File: C:\NOAA\NEMS_11731\src\atmos\gfs\dyn\noblas.f

```1           SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
2          \$                   BETA, C, LDC )
3     *     .. Scalar Arguments ..
4           CHARACTER*1        TRANSA, TRANSB
5           INTEGER            M, N, K, LDA, LDB, LDC
6           REAL               ALPHA, BETA
7     *     .. Array Arguments ..
8           REAL               A( LDA, * ), B( LDB, * ), C( LDC, * )
9     *     ..
10     *
11     *  Purpose
12     *  =======
13     *
14     *  DGEMM  performs one of the matrix-matrix operations
15     *
16     *     C := alpha*op( A )*op( B ) + beta*C,
17     *
18     *  where  op( X ) is one of
19     *
20     *     op( X ) = X   or   op( X ) = X',
21     *
22     *  alpha and beta are scalars, and A, B and C are matrices, with op( A )
23     *  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
24     *
25     *  Parameters
26     *  ==========
27     *
28     *  TRANSA - CHARACTER*1.
29     *           On entry, TRANSA specifies the form of op( A ) to be used in
30     *           the matrix multiplication as follows:
31     *
32     *              TRANSA = 'N' or 'n',  op( A ) = A.
33     *
34     *              TRANSA = 'T' or 't',  op( A ) = A'.
35     *
36     *              TRANSA = 'C' or 'c',  op( A ) = A'.
37     *
38     *           Unchanged on exit.
39     *
40     *  TRANSB - CHARACTER*1.
41     *           On entry, TRANSB specifies the form of op( B ) to be used in
42     *           the matrix multiplication as follows:
43     *
44     *              TRANSB = 'N' or 'n',  op( B ) = B.
45     *
46     *              TRANSB = 'T' or 't',  op( B ) = B'.
47     *
48     *              TRANSB = 'C' or 'c',  op( B ) = B'.
49     *
50     *           Unchanged on exit.
51     *
52     *  M      - INTEGER.
53     *           On entry,  M  specifies  the number  of rows  of the  matrix
54     *           op( A )  and of the  matrix  C.  M  must  be at least  zero.
55     *           Unchanged on exit.
56     *
57     *  N      - INTEGER.
58     *           On entry,  N  specifies the number  of columns of the matrix
59     *           op( B ) and the number of columns of the matrix C. N must be
60     *           at least zero.
61     *           Unchanged on exit.
62     *
63     *  K      - INTEGER.
64     *           On entry,  K  specifies  the number of columns of the matrix
65     *           op( A ) and the number of rows of the matrix op( B ). K must
66     *           be at least  zero.
67     *           Unchanged on exit.
68     *
69     *  ALPHA  - REAL            .
70     *           On entry, ALPHA specifies the scalar alpha.
71     *           Unchanged on exit.
72     *
73     *  A      - REAL             array of DIMENSION ( LDA, ka ), where ka is
74     *           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
75     *           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
76     *           part of the array  A  must contain the matrix  A,  otherwise
77     *           the leading  k by m  part of the array  A  must contain  the
78     *           matrix A.
79     *           Unchanged on exit.
80     *
81     *  LDA    - INTEGER.
82     *           On entry, LDA specifies the first dimension of A as declared
83     *           in the calling (sub) program. When  TRANSA = 'N' or 'n' then
84     *           LDA must be at least  max( 1, m ), otherwise  LDA must be at
85     *           least  max( 1, k ).
86     *           Unchanged on exit.
87     *
88     *  B      - REAL             array of DIMENSION ( LDB, kb ), where kb is
89     *           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
90     *           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
91     *           part of the array  B  must contain the matrix  B,  otherwise
92     *           the leading  n by k  part of the array  B  must contain  the
93     *           matrix B.
94     *           Unchanged on exit.
95     *
96     *  LDB    - INTEGER.
97     *           On entry, LDB specifies the first dimension of B as declared
98     *           in the calling (sub) program. When  TRANSB = 'N' or 'n' then
99     *           LDB must be at least  max( 1, k ), otherwise  LDB must be at
100     *           least  max( 1, n ).
101     *           Unchanged on exit.
102     *
103     *  BETA   - REAL            .
104     *           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
105     *           supplied as zero then C need not be set on input.
106     *           Unchanged on exit.
107     *
108     *  C      - REAL             array of DIMENSION ( LDC, n ).
109     *           Before entry, the leading  m by n  part of the array  C must
110     *           contain the matrix  C,  except when  beta  is zero, in which
111     *           case C need not be set on entry.
112     *           On exit, the array  C  is overwritten by the  m by n  matrix
113     *           ( alpha*op( A )*op( B ) + beta*C ).
114     *
115     *  LDC    - INTEGER.
116     *           On entry, LDC specifies the first dimension of C as declared
117     *           in  the  calling  (sub)  program.   LDC  must  be  at  least
118     *           max( 1, m ).
119     *           Unchanged on exit.
120     *
121     *
122     *  Level 3 Blas routine.
123     *
124     *  -- Written on 8-February-1989.
125     *     Jack Dongarra, Argonne National Laboratory.
126     *     Iain Duff, AERE Harwell.
127     *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
128     *     Sven Hammarling, Numerical Algorithms Group Ltd.
129     *
130     *
131     *     .. External Functions ..
132           LOGICAL            LSAME
133           EXTERNAL           LSAME
134     *     .. External Subroutines ..
135           EXTERNAL           XERBLA
136     *     .. Intrinsic Functions ..
137           INTRINSIC          MAX
138     *     .. Local Scalars ..
139           LOGICAL            NOTA, NOTB
140           INTEGER            I, INFO, J, L, NCOLA, NROWA, NROWB
141           REAL               TEMP
142     *     .. Parameters ..
143           REAL               ONE         , ZERO
144           PARAMETER        ( ONE = 1.0E+0, ZERO = 0.0E+0 )
145     *     ..
146     *     .. Executable Statements ..
147     *
148     *     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
149     *     transposed and set  NROWA, NCOLA and  NROWB  as the number of rows
150     *     and  columns of  A  and the  number of  rows  of  B  respectively.
151     *
152           NOTA  = LSAME( TRANSA, 'N' )
153           NOTB  = LSAME( TRANSB, 'N' )
154           IF( NOTA )THEN
155              NROWA = M
156              NCOLA = K
157           ELSE
158              NROWA = K
159              NCOLA = M
160           END IF
161           IF( NOTB )THEN
162              NROWB = K
163           ELSE
164              NROWB = N
165           END IF
166     *
167     *     Test the input parameters.
168     *
169           INFO = 0
170           IF(      ( .NOT.NOTA                 ).AND.
171          \$         ( .NOT.LSAME( TRANSA, 'C' ) ).AND.
172          \$         ( .NOT.LSAME( TRANSA, 'T' ) )      )THEN
173              INFO = 1
174           ELSE IF( ( .NOT.NOTB                 ).AND.
175          \$         ( .NOT.LSAME( TRANSB, 'C' ) ).AND.
176          \$         ( .NOT.LSAME( TRANSB, 'T' ) )      )THEN
177              INFO = 2
178           ELSE IF( M  .LT.0               )THEN
179              INFO = 3
180           ELSE IF( N  .LT.0               )THEN
181              INFO = 4
182           ELSE IF( K  .LT.0               )THEN
183              INFO = 5
184           ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
185              INFO = 8
186           ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN
187              INFO = 10
188           ELSE IF( LDC.LT.MAX( 1, M     ) )THEN
189              INFO = 13
190           END IF
191           IF( INFO.NE.0 )THEN
192              CALL XERBLA( 'DGEMM ', INFO )
193              RETURN
194           END IF
195     *
196     *     Quick return if possible.
197     *
198           IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
199          \$    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
200          \$   RETURN
201     *
202     *     And if  alpha.eq.zero.
203     *
204           IF( ALPHA.EQ.ZERO )THEN
205              IF( BETA.EQ.ZERO )THEN
206                 DO 20, J = 1, N
207                    DO 10, I = 1, M
208                       C( I, J ) = ZERO
209        10          CONTINUE
210        20       CONTINUE
211              ELSE
212                 DO 40, J = 1, N
213                    DO 30, I = 1, M
214                       C( I, J ) = BETA*C( I, J )
215        30          CONTINUE
216        40       CONTINUE
217              END IF
218              RETURN
219           END IF
220     *
221     *     Start the operations.
222     *
223           IF( NOTB )THEN
224              IF( NOTA )THEN
225     *
226     *           Form  C := alpha*A*B + beta*C.
227     *
228                 DO 90, J = 1, N
229                    IF( BETA.EQ.ZERO )THEN
230                       DO 50, I = 1, M
231                          C( I, J ) = ZERO
232        50             CONTINUE
233                    ELSE IF( BETA.NE.ONE )THEN
234                       DO 60, I = 1, M
235                          C( I, J ) = BETA*C( I, J )
236        60             CONTINUE
237                    END IF
238                    DO 80, L = 1, K
239                       IF( B( L, J ).NE.ZERO )THEN
240                          TEMP = ALPHA*B( L, J )
241                          DO 70, I = 1, M
242                             C( I, J ) = C( I, J ) + TEMP*A( I, L )
243        70                CONTINUE
244                       END IF
245        80          CONTINUE
246        90       CONTINUE
247              ELSE
248     *
249     *           Form  C := alpha*A'*B + beta*C
250     *
251                 DO 120, J = 1, N
252                    DO 110, I = 1, M
253                       TEMP = ZERO
254                       DO 100, L = 1, K
255                          TEMP = TEMP + A( L, I )*B( L, J )
256       100             CONTINUE
257                       IF( BETA.EQ.ZERO )THEN
258                          C( I, J ) = ALPHA*TEMP
259                       ELSE
260                          C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
261                       END IF
262       110          CONTINUE
263       120       CONTINUE
264              END IF
265           ELSE
266              IF( NOTA )THEN
267     *
268     *           Form  C := alpha*A*B' + beta*C
269     *
270                 DO 170, J = 1, N
271                    IF( BETA.EQ.ZERO )THEN
272                       DO 130, I = 1, M
273                          C( I, J ) = ZERO
274       130             CONTINUE
275                    ELSE IF( BETA.NE.ONE )THEN
276                       DO 140, I = 1, M
277                          C( I, J ) = BETA*C( I, J )
278       140             CONTINUE
279                    END IF
280                    DO 160, L = 1, K
281                       IF( B( J, L ).NE.ZERO )THEN
282                          TEMP = ALPHA*B( J, L )
283                          DO 150, I = 1, M
284                             C( I, J ) = C( I, J ) + TEMP*A( I, L )
285       150                CONTINUE
286                       END IF
287       160          CONTINUE
288       170       CONTINUE
289              ELSE
290     *
291     *           Form  C := alpha*A'*B' + beta*C
292     *
293                 DO 200, J = 1, N
294                    DO 190, I = 1, M
295                       TEMP = ZERO
296                       DO 180, L = 1, K
297                          TEMP = TEMP + A( L, I )*B( J, L )
298       180             CONTINUE
299                       IF( BETA.EQ.ZERO )THEN
300                          C( I, J ) = ALPHA*TEMP
301                       ELSE
302                          C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
303                       END IF
304       190          CONTINUE
305       200       CONTINUE
306              END IF
307           END IF
308     *
309           RETURN
310     *
311     *     End of DGEMM .
312     *
313           END
314           SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
315     *
316     *  -- LAPACK routine (version 2.0) --
317     *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
318     *     Courant Institute, Argonne National Lab, and Rice University
319     *     September 30, 1994
320     *
321     *     .. Scalar Arguments ..
322           INTEGER            INFO, LDA, LWORK, N
323     *     ..
324     *     .. Array Arguments ..
325           INTEGER            IPIV( * )
326           REAL               A( LDA, * ), WORK( LWORK )
327     *     ..
328     *
329     *  Purpose
330     *  =======
331     *
332     *  DGETRI computes the inverse of a matrix using the LU factorization
333     *  computed by DGETRF.
334     *
335     *  This method inverts U and then computes inv(A) by solving the system
336     *  inv(A)*L = inv(U) for inv(A).
337     *
338     *  Arguments
339     *  =========
340     *
341     *  N       (input) INTEGER
342     *          The order of the matrix A.  N >= 0.
343     *
344     *  A       (input/output) REAL array, dimension (LDA,N)
345     *          On entry, the factors L and U from the factorization
346     *          A = P*L*U as computed by DGETRF.
347     *          On exit, if INFO = 0, the inverse of the original matrix A.
348     *
349     *  LDA     (input) INTEGER
350     *          The leading dimension of the array A.  LDA >= max(1,N).
351     *
352     *  IPIV    (input) INTEGER array, dimension (N)
353     *          The pivot indices from DGETRF; for 1<=i<=N, row i of the
354     *          matrix was interchanged with row IPIV(i).
355     *
356     *  WORK    (workspace/output) REAL array, dimension (LWORK)
357     *          On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
358     *
359     *  LWORK   (input) INTEGER
360     *          The dimension of the array WORK.  LWORK >= max(1,N).
361     *          For optimal performance LWORK >= N*NB, where NB is
362     *          the optimal blocksize returned by ILAENV.
363     *
364     *  INFO    (output) INTEGER
365     *          = 0:  successful exit
366     *          < 0:  if INFO = -i, the i-th argument had an illegal value
367     *          > 0:  if INFO = i, U(i,i) is exactly zero; the matrix is
368     *                singular and its inverse could not be computed.
369     *
370     *  =====================================================================
371     *
372     *     .. Parameters ..
373           REAL               ZERO, ONE
374           PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
375     *     ..
376     *     .. Local Scalars ..
377           INTEGER            I, IWS, J, JB, JJ, JP, LDWORK, NB, NBMIN, NN
378     *     ..
379     *     .. External Functions ..
380           INTEGER            ILAENV
381           EXTERNAL           ILAENV
382     *     ..
383     *     .. External Subroutines ..
384           EXTERNAL           DGEMM, DGEMV, DSWAP, DTRSM, DTRTRI, XERBLA
385     *     ..
386     *     .. Intrinsic Functions ..
387           INTRINSIC          MAX, MIN
388     *     ..
389     *     .. Executable Statements ..
390     *
391     *     Test the input parameters.
392     *
393           INFO = 0
394           WORK( 1 ) = MAX( N, 1 )
395           IF( N.LT.0 ) THEN
396              INFO = -1
397           ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
398              INFO = -3
399           ELSE IF( LWORK.LT.MAX( 1, N ) ) THEN
400              INFO = -6
401           END IF
402           IF( INFO.NE.0 ) THEN
403              CALL XERBLA( 'DGETRI', -INFO )
404              RETURN
405           END IF
406     *
407     *     Quick return if possible
408     *
409           IF( N.EQ.0 )
410          \$   RETURN
411     *
412     *     Form inv(U).  If INFO > 0 from DTRTRI, then U is singular,
413     *     and the inverse is not computed.
414     *
415           CALL DTRTRI( 'Upper', 'Non-unit', N, A, LDA, INFO )
416           IF( INFO.GT.0 )
417          \$   RETURN
418     *
419     *     Determine the block size for this environment.
420     *
421           NB = ILAENV( 1, 'DGETRI', ' ', N, -1, -1, -1 )
422           NBMIN = 2
423           LDWORK = N
424           IF( NB.GT.1 .AND. NB.LT.N ) THEN
425              IWS = MAX( LDWORK*NB, 1 )
426              IF( LWORK.LT.IWS ) THEN
427                 NB = LWORK / LDWORK
428                 NBMIN = MAX( 2, ILAENV( 2, 'DGETRI', ' ', N, -1, -1, -1 ) )
429              END IF
430           ELSE
431              IWS = N
432           END IF
433     *
434     *     Solve the equation inv(A)*L = inv(U) for inv(A).
435     *
436           IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN
437     *
438     *        Use unblocked code.
439     *
440              DO 20 J = N, 1, -1
441     *
442     *           Copy current column of L to WORK and replace with zeros.
443     *
444                 DO 10 I = J + 1, N
445                    WORK( I ) = A( I, J )
446                    A( I, J ) = ZERO
447        10       CONTINUE
448     *
449     *           Compute current column of inv(A).
450     *
451                 IF( J.LT.N )
452          \$         CALL DGEMV( 'No transpose', N, N-J, -ONE, A( 1, J+1 ),
453          \$                     LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 )
454        20    CONTINUE
455           ELSE
456     *
457     *        Use blocked code.
458     *
459              NN = ( ( N-1 ) / NB )*NB + 1
460              DO 50 J = NN, 1, -NB
461                 JB = MIN( NB, N-J+1 )
462     *
463     *           Copy current block column of L to WORK and replace with
464     *           zeros.
465     *
466                 DO 40 JJ = J, J + JB - 1
467                    DO 30 I = JJ + 1, N
468                       WORK( I+( JJ-J )*LDWORK ) = A( I, JJ )
469                       A( I, JJ ) = ZERO
470        30          CONTINUE
471        40       CONTINUE
472     *
473     *           Compute current block column of inv(A).
474     *
475                 IF( J+JB.LE.N )
476          \$         CALL DGEMM( 'No transpose', 'No transpose', N, JB,
477          \$                     N-J-JB+1, -ONE, A( 1, J+JB ), LDA,
478          \$                     WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA )
479                 CALL DTRSM( 'Right', 'Lower', 'No transpose', 'Unit', N, JB,
480          \$                  ONE, WORK( J ), LDWORK, A( 1, J ), LDA )
481        50    CONTINUE
482           END IF
483     *
484     *     Apply column interchanges.
485     *
486           DO 60 J = N - 1, 1, -1
487              JP = IPIV( J )
488              IF( JP.NE.J )
489          \$      CALL DSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 )
490        60 CONTINUE
491     *
492           WORK( 1 ) = IWS
493           RETURN
494     *
495     *     End of DGETRI
496     *
497           END
498           SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO )
499     *
500     *  -- LAPACK routine (version 2.0) --
501     *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
502     *     Courant Institute, Argonne National Lab, and Rice University
503     *     March 31, 1993
504     *
505     *     .. Scalar Arguments ..
506           INTEGER            INFO, LDA, M, N
507     *     ..
508     *     .. Array Arguments ..
509           INTEGER            IPIV( * )
510           REAL               A( LDA, * )
511     *     ..
512     *
513     *  Purpose
514     *  =======
515     *
516     *  DGETRF computes an LU factorization of a general M-by-N matrix A
517     *  using partial pivoting with row interchanges.
518     *
519     *  The factorization has the form
520     *     A = P * L * U
521     *  where P is a permutation matrix, L is lower triangular with unit
522     *  diagonal elements (lower trapezoidal if m > n), and U is upper
523     *  triangular (upper trapezoidal if m < n).
524     *
525     *  This is the right-looking Level 3 BLAS version of the algorithm.
526     *
527     *  Arguments
528     *  =========
529     *
530     *  M       (input) INTEGER
531     *          The number of rows of the matrix A.  M >= 0.
532     *
533     *  N       (input) INTEGER
534     *          The number of columns of the matrix A.  N >= 0.
535     *
536     *  A       (input/output) REAL array, dimension (LDA,N)
537     *          On entry, the M-by-N matrix to be factored.
538     *          On exit, the factors L and U from the factorization
539     *          A = P*L*U; the unit diagonal elements of L are not stored.
540     *
541     *  LDA     (input) INTEGER
542     *          The leading dimension of the array A.  LDA >= max(1,M).
543     *
544     *  IPIV    (output) INTEGER array, dimension (min(M,N))
545     *          The pivot indices; for 1 <= i <= min(M,N), row i of the
546     *          matrix was interchanged with row IPIV(i).
547     *
548     *  INFO    (output) INTEGER
549     *          = 0:  successful exit
550     *          < 0:  if INFO = -i, the i-th argument had an illegal value
551     *          > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
552     *                has been completed, but the factor U is exactly
553     *                singular, and division by zero will occur if it is used
554     *                to solve a system of equations.
555     *
556     *  =====================================================================
557     *
558     *     .. Parameters ..
559           REAL               ONE
560           PARAMETER          ( ONE = 1.0E+0 )
561     *     ..
562     *     .. Local Scalars ..
563           INTEGER            I, IINFO, J, JB, NB
564     *     ..
565     *     .. External Subroutines ..
566           EXTERNAL           DGEMM, DGETF2, DLASWP, DTRSM, XERBLA
567     *     ..
568     *     .. External Functions ..
569           INTEGER            ILAENV
570           EXTERNAL           ILAENV
571     *     ..
572     *     .. Intrinsic Functions ..
573           INTRINSIC          MAX, MIN
574     *     ..
575     *     .. Executable Statements ..
576     *
577     *     Test the input parameters.
578     *
579           INFO = 0
580           IF( M.LT.0 ) THEN
581              INFO = -1
582           ELSE IF( N.LT.0 ) THEN
583              INFO = -2
584           ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
585              INFO = -4
586           END IF
587           IF( INFO.NE.0 ) THEN
588              CALL XERBLA( 'DGETRF', -INFO )
589              RETURN
590           END IF
591     *
592     *     Quick return if possible
593     *
594           IF( M.EQ.0 .OR. N.EQ.0 )
595          \$   RETURN
596     *
597     *     Determine the block size for this environment.
598     *
599           NB = ILAENV( 1, 'DGETRF', ' ', M, N, -1, -1 )
600           IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN
601     *
602     *        Use unblocked code.
603     *
604              CALL DGETF2( M, N, A, LDA, IPIV, INFO )
605           ELSE
606     *
607     *        Use blocked code.
608     *
609              DO 20 J = 1, MIN( M, N ), NB
610                 JB = MIN( MIN( M, N )-J+1, NB )
611     *
612     *           Factor diagonal and subdiagonal blocks and test for exact
613     *           singularity.
614     *
615                 CALL DGETF2( M-J+1, JB, A( J, J ), LDA, IPIV( J ), IINFO )
616     *
617     *           Adjust INFO and the pivot indices.
618     *
619                 IF( INFO.EQ.0 .AND. IINFO.GT.0 )
620          \$         INFO = IINFO + J - 1
621                 DO 10 I = J, MIN( M, J+JB-1 )
622                    IPIV( I ) = J - 1 + IPIV( I )
623        10       CONTINUE
624     *
625     *           Apply interchanges to columns 1:J-1.
626     *
627                 CALL DLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 )
628     *
629                 IF( J+JB.LE.N ) THEN
630     *
631     *              Apply interchanges to columns J+JB:N.
632     *
633                    CALL DLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1,
634          \$                      IPIV, 1 )
635     *
636     *              Compute block row of U.
637     *
638                    CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB,
639          \$                     N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ),
640          \$                     LDA )
641                    IF( J+JB.LE.M ) THEN
642     *
643     *                 Update trailing submatrix.
644     *
645                       CALL DGEMM( 'No transpose', 'No transpose', M-J-JB+1,
646          \$                        N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA,
647          \$                        A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ),
648          \$                        LDA )
649                    END IF
650                 END IF
651        20    CONTINUE
652           END IF
653           RETURN
654     *
655     *     End of DGETRF
656     *
657           END
658           SUBROUTINE XERBLA( SRNAME, INFO )
659     *
660     *  -- LAPACK auxiliary routine (preliminary version) --
661     *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
662     *     Courant Institute, Argonne National Lab, and Rice University
663     *     February 29, 1992
664     *
665     *     .. Scalar Arguments ..
666           CHARACTER*6        SRNAME
667           INTEGER            INFO
668     *     ..
669     *
670     *  Purpose
671     *  =======
672     *
673     *  XERBLA  is an error handler for the LAPACK routines.
674     *  It is called by an LAPACK routine if an input parameter has an
675     *  invalid value.  A message is printed and execution stops.
676     *
677     *  Installers may consider modifying the STOP statement in order to
678     *  call system-specific exception-handling facilities.
679     *
680     *  Arguments
681     *  =========
682     *
683     *  SRNAME  (input) CHARACTER*6
684     *          The name of the routine which called XERBLA.
685     *
686     *  INFO    (input) INTEGER
687     *          The position of the invalid parameter in the parameter list
688     *          of the calling routine.
689     *
690     *
691           WRITE( *, FMT = 9999 )SRNAME, INFO
692     *
693           STOP
694     *
695      9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ',
696          \$      'an illegal value' )
697     *
698     *     End of XERBLA
699     *
700           END
701           SUBROUTINE DTRTRI( UPLO, DIAG, N, A, LDA, INFO )
702     *
703     *  -- LAPACK routine (version 2.0) --
704     *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
705     *     Courant Institute, Argonne National Lab, and Rice University
706     *     March 31, 1993
707     *
708     *     .. Scalar Arguments ..
709           CHARACTER          DIAG, UPLO
710           INTEGER            INFO, LDA, N
711     *     ..
712     *     .. Array Arguments ..
713           REAL               A( LDA, * )
714     *     ..
715     *
716     *  Purpose
717     *  =======
718     *
719     *  DTRTRI computes the inverse of a real upper or lower triangular
720     *  matrix A.
721     *
722     *  This is the Level 3 BLAS version of the algorithm.
723     *
724     *  Arguments
725     *  =========
726     *
727     *  UPLO    (input) CHARACTER*1
728     *          = 'U':  A is upper triangular;
729     *          = 'L':  A is lower triangular.
730     *
731     *  DIAG    (input) CHARACTER*1
732     *          = 'N':  A is non-unit triangular;
733     *          = 'U':  A is unit triangular.
734     *
735     *  N       (input) INTEGER
736     *          The order of the matrix A.  N >= 0.
737     *
738     *  A       (input/output) REAL array, dimension (LDA,N)
739     *          On entry, the triangular matrix A.  If UPLO = 'U', the
740     *          leading N-by-N upper triangular part of the array A contains
741     *          the upper triangular matrix, and the strictly lower
742     *          triangular part of A is not referenced.  If UPLO = 'L', the
743     *          leading N-by-N lower triangular part of the array A contains
744     *          the lower triangular matrix, and the strictly upper
745     *          triangular part of A is not referenced.  If DIAG = 'U', the
746     *          diagonal elements of A are also not referenced and are
747     *          assumed to be 1.
748     *          On exit, the (triangular) inverse of the original matrix, in
749     *          the same storage format.
750     *
751     *  LDA     (input) INTEGER
752     *          The leading dimension of the array A.  LDA >= max(1,N).
753     *
754     *  INFO    (output) INTEGER
755     *          = 0: successful exit
756     *          < 0: if INFO = -i, the i-th argument had an illegal value
757     *          > 0: if INFO = i, A(i,i) is exactly zero.  The triangular
758     *               matrix is singular and its inverse can not be computed.
759     *
760     *  =====================================================================
761     *
762     *     .. Parameters ..
763           REAL               ONE, ZERO
764           PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
765     *     ..
766     *     .. Local Scalars ..
767           LOGICAL            NOUNIT, UPPER
768           INTEGER            J, JB, NB, NN
769     *     ..
770     *     .. External Functions ..
771           LOGICAL            LSAME
772           INTEGER            ILAENV
773           EXTERNAL           LSAME, ILAENV
774     *     ..
775     *     .. External Subroutines ..
776           EXTERNAL           DTRMM, DTRSM, DTRTI2, XERBLA
777     *     ..
778     *     .. Intrinsic Functions ..
779           INTRINSIC          MAX, MIN
780     *     ..
781     *     .. Executable Statements ..
782     *
783     *     Test the input parameters.
784     *
785           INFO = 0
786           UPPER = LSAME( UPLO, 'U' )
787           NOUNIT = LSAME( DIAG, 'N' )
788           IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
789              INFO = -1
790           ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
791              INFO = -2
792           ELSE IF( N.LT.0 ) THEN
793              INFO = -3
794           ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
795              INFO = -5
796           END IF
797           IF( INFO.NE.0 ) THEN
798              CALL XERBLA( 'DTRTRI', -INFO )
799              RETURN
800           END IF
801     *
802     *     Quick return if possible
803     *
804           IF( N.EQ.0 )
805          \$   RETURN
806     *
807     *     Check for singularity if non-unit.
808     *
809           IF( NOUNIT ) THEN
810              DO 10 INFO = 1, N
811                 IF( A( INFO, INFO ).EQ.ZERO )
812          \$         RETURN
813        10    CONTINUE
814              INFO = 0
815           END IF
816     *
817     *     Determine the block size for this environment.
818     *
819           NB = ILAENV( 1, 'DTRTRI', UPLO // DIAG, N, -1, -1, -1 )
820           IF( NB.LE.1 .OR. NB.GE.N ) THEN
821     *
822     *        Use unblocked code
823     *
824              CALL DTRTI2( UPLO, DIAG, N, A, LDA, INFO )
825           ELSE
826     *
827     *        Use blocked code
828     *
829              IF( UPPER ) THEN
830     *
831     *           Compute inverse of upper triangular matrix
832     *
833                 DO 20 J = 1, N, NB
834                    JB = MIN( NB, N-J+1 )
835     *
836     *              Compute rows 1:j-1 of current block column
837     *
838                    CALL DTRMM( 'Left', 'Upper', 'No transpose', DIAG, J-1,
839          \$                     JB, ONE, A, LDA, A( 1, J ), LDA )
840                    CALL DTRSM( 'Right', 'Upper', 'No transpose', DIAG, J-1,
841          \$                     JB, -ONE, A( J, J ), LDA, A( 1, J ), LDA )
842     *
843     *              Compute inverse of current diagonal block
844     *
845                    CALL DTRTI2( 'Upper', DIAG, JB, A( J, J ), LDA, INFO )
846        20       CONTINUE
847              ELSE
848     *
849     *           Compute inverse of lower triangular matrix
850     *
851                 NN = ( ( N-1 ) / NB )*NB + 1
852                 DO 30 J = NN, 1, -NB
853                    JB = MIN( NB, N-J+1 )
854                    IF( J+JB.LE.N ) THEN
855     *
856     *                 Compute rows j+jb:n of current block column
857     *
858                       CALL DTRMM( 'Left', 'Lower', 'No transpose', DIAG,
859          \$                        N-J-JB+1, JB, ONE, A( J+JB, J+JB ), LDA,
860          \$                        A( J+JB, J ), LDA )
861                       CALL DTRSM( 'Right', 'Lower', 'No transpose', DIAG,
862          \$                        N-J-JB+1, JB, -ONE, A( J, J ), LDA,
863          \$                        A( J+JB, J ), LDA )
864                    END IF
865     *
866     *              Compute inverse of current diagonal block
867     *
868                    CALL DTRTI2( 'Lower', DIAG, JB, A( J, J ), LDA, INFO )
869        30       CONTINUE
870              END IF
871           END IF
872     *
873           RETURN
874     *
875     *     End of DTRTRI
876     *
877           END
878           SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
879          \$                   BETA, Y, INCY )
880     *     .. Scalar Arguments ..
881           REAL               ALPHA, BETA
882           INTEGER            INCX, INCY, LDA, M, N
883           CHARACTER*1        TRANS
884     *     .. Array Arguments ..
885           REAL               A( LDA, * ), X( * ), Y( * )
886     *     ..
887     *
888     *  Purpose
889     *  =======
890     *
891     *  DGEMV  performs one of the matrix-vector operations
892     *
893     *     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
894     *
895     *  where alpha and beta are scalars, x and y are vectors and A is an
896     *  m by n matrix.
897     *
898     *  Parameters
899     *  ==========
900     *
901     *  TRANS  - CHARACTER*1.
902     *           On entry, TRANS specifies the operation to be performed as
903     *           follows:
904     *
905     *              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
906     *
907     *              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
908     *
909     *              TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.
910     *
911     *           Unchanged on exit.
912     *
913     *  M      - INTEGER.
914     *           On entry, M specifies the number of rows of the matrix A.
915     *           M must be at least zero.
916     *           Unchanged on exit.
917     *
918     *  N      - INTEGER.
919     *           On entry, N specifies the number of columns of the matrix A.
920     *           N must be at least zero.
921     *           Unchanged on exit.
922     *
923     *  ALPHA  - REAL            .
924     *           On entry, ALPHA specifies the scalar alpha.
925     *           Unchanged on exit.
926     *
927     *  A      - REAL             array of DIMENSION ( LDA, n ).
928     *           Before entry, the leading m by n part of the array A must
929     *           contain the matrix of coefficients.
930     *           Unchanged on exit.
931     *
932     *  LDA    - INTEGER.
933     *           On entry, LDA specifies the first dimension of A as declared
934     *           in the calling (sub) program. LDA must be at least
935     *           max( 1, m ).
936     *           Unchanged on exit.
937     *
938     *  X      - REAL             array of DIMENSION at least
939     *           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
940     *           and at least
941     *           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
942     *           Before entry, the incremented array X must contain the
943     *           vector x.
944     *           Unchanged on exit.
945     *
946     *  INCX   - INTEGER.
947     *           On entry, INCX specifies the increment for the elements of
948     *           X. INCX must not be zero.
949     *           Unchanged on exit.
950     *
951     *  BETA   - REAL            .
952     *           On entry, BETA specifies the scalar beta. When BETA is
953     *           supplied as zero then Y need not be set on input.
954     *           Unchanged on exit.
955     *
956     *  Y      - REAL             array of DIMENSION at least
957     *           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
958     *           and at least
959     *           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
960     *           Before entry with BETA non-zero, the incremented array Y
961     *           must contain the vector y. On exit, Y is overwritten by the
962     *           updated vector y.
963     *
964     *  INCY   - INTEGER.
965     *           On entry, INCY specifies the increment for the elements of
966     *           Y. INCY must not be zero.
967     *           Unchanged on exit.
968     *
969     *
970     *  Level 2 Blas routine.
971     *
972     *  -- Written on 22-October-1986.
973     *     Jack Dongarra, Argonne National Lab.
974     *     Jeremy Du Croz, Nag Central Office.
975     *     Sven Hammarling, Nag Central Office.
976     *     Richard Hanson, Sandia National Labs.
977     *
978     *
979     *     .. Parameters ..
980           REAL               ONE         , ZERO
981           PARAMETER        ( ONE = 1.0E+0, ZERO = 0.0E+0 )
982     *     .. Local Scalars ..
983           REAL               TEMP
984           INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY
985     *     .. External Functions ..
986           LOGICAL            LSAME
987           EXTERNAL           LSAME
988     *     .. External Subroutines ..
989           EXTERNAL           XERBLA
990     *     .. Intrinsic Functions ..
991           INTRINSIC          MAX
992     *     ..
993     *     .. Executable Statements ..
994     *
995     *     Test the input parameters.
996     *
997           INFO = 0
998           IF     ( .NOT.LSAME( TRANS, 'N' ).AND.
999          \$         .NOT.LSAME( TRANS, 'T' ).AND.
1000          \$         .NOT.LSAME( TRANS, 'C' )      )THEN
1001              INFO = 1
1002           ELSE IF( M.LT.0 )THEN
1003              INFO = 2
1004           ELSE IF( N.LT.0 )THEN
1005              INFO = 3
1006           ELSE IF( LDA.LT.MAX( 1, M ) )THEN
1007              INFO = 6
1008           ELSE IF( INCX.EQ.0 )THEN
1009              INFO = 8
1010           ELSE IF( INCY.EQ.0 )THEN
1011              INFO = 11
1012           END IF
1013           IF( INFO.NE.0 )THEN
1014              CALL XERBLA( 'DGEMV ', INFO )
1015              RETURN
1016           END IF
1017     *
1018     *     Quick return if possible.
1019     *
1020           IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
1021          \$    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
1022          \$   RETURN
1023     *
1024     *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
1025     *     up the start points in  X  and  Y.
1026     *
1027           IF( LSAME( TRANS, 'N' ) )THEN
1028              LENX = N
1029              LENY = M
1030           ELSE
1031              LENX = M
1032              LENY = N
1033           END IF
1034           IF( INCX.GT.0 )THEN
1035              KX = 1
1036           ELSE
1037              KX = 1 - ( LENX - 1 )*INCX
1038           END IF
1039           IF( INCY.GT.0 )THEN
1040              KY = 1
1041           ELSE
1042              KY = 1 - ( LENY - 1 )*INCY
1043           END IF
1044     *
1045     *     Start the operations. In this version the elements of A are
1046     *     accessed sequentially with one pass through A.
1047     *
1048     *     First form  y := beta*y.
1049     *
1050           IF( BETA.NE.ONE )THEN
1051              IF( INCY.EQ.1 )THEN
1052                 IF( BETA.EQ.ZERO )THEN
1053                    DO 10, I = 1, LENY
1054                       Y( I ) = ZERO
1055        10          CONTINUE
1056                 ELSE
1057                    DO 20, I = 1, LENY
1058                       Y( I ) = BETA*Y( I )
1059        20          CONTINUE
1060                 END IF
1061              ELSE
1062                 IY = KY
1063                 IF( BETA.EQ.ZERO )THEN
1064                    DO 30, I = 1, LENY
1065                       Y( IY ) = ZERO
1066                       IY      = IY   + INCY
1067        30          CONTINUE
1068                 ELSE
1069                    DO 40, I = 1, LENY
1070                       Y( IY ) = BETA*Y( IY )
1071                       IY      = IY           + INCY
1072        40          CONTINUE
1073                 END IF
1074              END IF
1075           END IF
1076           IF( ALPHA.EQ.ZERO )
1077          \$   RETURN
1078           IF( LSAME( TRANS, 'N' ) )THEN
1079     *
1080     *        Form  y := alpha*A*x + y.
1081     *
1082              JX = KX
1083              IF( INCY.EQ.1 )THEN
1084                 DO 60, J = 1, N
1085                    IF( X( JX ).NE.ZERO )THEN
1086                       TEMP = ALPHA*X( JX )
1087                       DO 50, I = 1, M
1088                          Y( I ) = Y( I ) + TEMP*A( I, J )
1089        50             CONTINUE
1090                    END IF
1091                    JX = JX + INCX
1092        60       CONTINUE
1093              ELSE
1094                 DO 80, J = 1, N
1095                    IF( X( JX ).NE.ZERO )THEN
1096                       TEMP = ALPHA*X( JX )
1097                       IY   = KY
1098                       DO 70, I = 1, M
1099                          Y( IY ) = Y( IY ) + TEMP*A( I, J )
1100                          IY      = IY      + INCY
1101        70             CONTINUE
1102                    END IF
1103                    JX = JX + INCX
1104        80       CONTINUE
1105              END IF
1106           ELSE
1107     *
1108     *        Form  y := alpha*A'*x + y.
1109     *
1110              JY = KY
1111              IF( INCX.EQ.1 )THEN
1112                 DO 100, J = 1, N
1113                    TEMP = ZERO
1114                    DO 90, I = 1, M
1115                       TEMP = TEMP + A( I, J )*X( I )
1116        90          CONTINUE
1117                    Y( JY ) = Y( JY ) + ALPHA*TEMP
1118                    JY      = JY      + INCY
1119       100       CONTINUE
1120              ELSE
1121                 DO 120, J = 1, N
1122                    TEMP = ZERO
1123                    IX   = KX
1124                    DO 110, I = 1, M
1125                       TEMP = TEMP + A( I, J )*X( IX )
1126                       IX   = IX   + INCX
1127       110          CONTINUE
1128                    Y( JY ) = Y( JY ) + ALPHA*TEMP
1129                    JY      = JY      + INCY
1130       120       CONTINUE
1131              END IF
1132           END IF
1133     *
1134           RETURN
1135     *
1136     *     End of DGEMV .
1137     *
1138           END
1139           SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
1140          \$                   B, LDB )
1141     *     .. Scalar Arguments ..
1142           CHARACTER*1        SIDE, UPLO, TRANSA, DIAG
1143           INTEGER            M, N, LDA, LDB
1144           REAL               ALPHA
1145     *     .. Array Arguments ..
1146           REAL               A( LDA, * ), B( LDB, * )
1147     *     ..
1148     *
1149     *  Purpose
1150     *  =======
1151     *
1152     *  DTRSM  solves one of the matrix equations
1153     *
1154     *     op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
1155     *
1156     *  where alpha is a scalar, X and B are m by n matrices, A is a unit, or
1157     *  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
1158     *
1159     *     op( A ) = A   or   op( A ) = A'.
1160     *
1161     *  The matrix X is overwritten on B.
1162     *
1163     *  Parameters
1164     *  ==========
1165     *
1166     *  SIDE   - CHARACTER*1.
1167     *           On entry, SIDE specifies whether op( A ) appears on the left
1168     *           or right of X as follows:
1169     *
1170     *              SIDE = 'L' or 'l'   op( A )*X = alpha*B.
1171     *
1172     *              SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
1173     *
1174     *           Unchanged on exit.
1175     *
1176     *  UPLO   - CHARACTER*1.
1177     *           On entry, UPLO specifies whether the matrix A is an upper or
1178     *           lower triangular matrix as follows:
1179     *
1180     *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
1181     *
1182     *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
1183     *
1184     *           Unchanged on exit.
1185     *
1186     *  TRANSA - CHARACTER*1.
1187     *           On entry, TRANSA specifies the form of op( A ) to be used in
1188     *           the matrix multiplication as follows:
1189     *
1190     *              TRANSA = 'N' or 'n'   op( A ) = A.
1191     *
1192     *              TRANSA = 'T' or 't'   op( A ) = A'.
1193     *
1194     *              TRANSA = 'C' or 'c'   op( A ) = A'.
1195     *
1196     *           Unchanged on exit.
1197     *
1198     *  DIAG   - CHARACTER*1.
1199     *           On entry, DIAG specifies whether or not A is unit triangular
1200     *           as follows:
1201     *
1202     *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
1203     *
1204     *              DIAG = 'N' or 'n'   A is not assumed to be unit
1205     *                                  triangular.
1206     *
1207     *           Unchanged on exit.
1208     *
1209     *  M      - INTEGER.
1210     *           On entry, M specifies the number of rows of B. M must be at
1211     *           least zero.
1212     *           Unchanged on exit.
1213     *
1214     *  N      - INTEGER.
1215     *           On entry, N specifies the number of columns of B.  N must be
1216     *           at least zero.
1217     *           Unchanged on exit.
1218     *
1219     *  ALPHA  - REAL            .
1220     *           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
1221     *           zero then  A is not referenced and  B need not be set before
1222     *           entry.
1223     *           Unchanged on exit.
1224     *
1225     *  A      - REAL             array of DIMENSION ( LDA, k ), where k is m
1226     *           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
1227     *           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
1228     *           upper triangular part of the array  A must contain the upper
1229     *           triangular matrix  and the strictly lower triangular part of
1230     *           A is not referenced.
1231     *           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
1232     *           lower triangular part of the array  A must contain the lower
1233     *           triangular matrix  and the strictly upper triangular part of
1234     *           A is not referenced.
1235     *           Note that when  DIAG = 'U' or 'u',  the diagonal elements of
1236     *           A  are not referenced either,  but are assumed to be  unity.
1237     *           Unchanged on exit.
1238     *
1239     *  LDA    - INTEGER.
1240     *           On entry, LDA specifies the first dimension of A as declared
1241     *           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
1242     *           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
1243     *           then LDA must be at least max( 1, n ).
1244     *           Unchanged on exit.
1245     *
1246     *  B      - REAL             array of DIMENSION ( LDB, n ).
1247     *           Before entry,  the leading  m by n part of the array  B must
1248     *           contain  the  right-hand  side  matrix  B,  and  on exit  is
1249     *           overwritten by the solution matrix  X.
1250     *
1251     *  LDB    - INTEGER.
1252     *           On entry, LDB specifies the first dimension of B as declared
1253     *           in  the  calling  (sub)  program.   LDB  must  be  at  least
1254     *           max( 1, m ).
1255     *           Unchanged on exit.
1256     *
1257     *
1258     *  Level 3 Blas routine.
1259     *
1260     *
1261     *  -- Written on 8-February-1989.
1262     *     Jack Dongarra, Argonne National Laboratory.
1263     *     Iain Duff, AERE Harwell.
1264     *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
1265     *     Sven Hammarling, Numerical Algorithms Group Ltd.
1266     *
1267     *
1268     *     .. External Functions ..
1269           LOGICAL            LSAME
1270           EXTERNAL           LSAME
1271     *     .. External Subroutines ..
1272           EXTERNAL           XERBLA
1273     *     .. Intrinsic Functions ..
1274           INTRINSIC          MAX
1275     *     .. Local Scalars ..
1276           LOGICAL            LSIDE, NOUNIT, UPPER
1277           INTEGER            I, INFO, J, K, NROWA
1278           REAL               TEMP
1279     *     .. Parameters ..
1280           REAL               ONE         , ZERO
1281           PARAMETER        ( ONE = 1.0E+0, ZERO = 0.0E+0 )
1282     *     ..
1283     *     .. Executable Statements ..
1284     *
1285     *     Test the input parameters.
1286     *
1287           LSIDE  = LSAME( SIDE  , 'L' )
1288           IF( LSIDE )THEN
1289              NROWA = M
1290           ELSE
1291              NROWA = N
1292           END IF
1293           NOUNIT = LSAME( DIAG  , 'N' )
1294           UPPER  = LSAME( UPLO  , 'U' )
1295     *
1296           INFO   = 0
1297           IF(      ( .NOT.LSIDE                ).AND.
1298          \$         ( .NOT.LSAME( SIDE  , 'R' ) )      )THEN
1299              INFO = 1
1300           ELSE IF( ( .NOT.UPPER                ).AND.
1301          \$         ( .NOT.LSAME( UPLO  , 'L' ) )      )THEN
1302              INFO = 2
1303           ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
1304          \$         ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
1305          \$         ( .NOT.LSAME( TRANSA, 'C' ) )      )THEN
1306              INFO = 3
1307           ELSE IF( ( .NOT.LSAME( DIAG  , 'U' ) ).AND.
1308          \$         ( .NOT.LSAME( DIAG  , 'N' ) )      )THEN
1309              INFO = 4
1310           ELSE IF( M  .LT.0               )THEN
1311              INFO = 5
1312           ELSE IF( N  .LT.0               )THEN
1313              INFO = 6
1314           ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
1315              INFO = 9
1316           ELSE IF( LDB.LT.MAX( 1, M     ) )THEN
1317              INFO = 11
1318           END IF
1319           IF( INFO.NE.0 )THEN
1320              CALL XERBLA( 'DTRSM ', INFO )
1321              RETURN
1322           END IF
1323     *
1324     *     Quick return if possible.
1325     *
1326           IF( N.EQ.0 )
1327          \$   RETURN
1328     *
1329     *     And when  alpha.eq.zero.
1330     *
1331           IF( ALPHA.EQ.ZERO )THEN
1332              DO 20, J = 1, N
1333                 DO 10, I = 1, M
1334                    B( I, J ) = ZERO
1335        10       CONTINUE
1336        20    CONTINUE
1337              RETURN
1338           END IF
1339     *
1340     *     Start the operations.
1341     *
1342           IF( LSIDE )THEN
1343              IF( LSAME( TRANSA, 'N' ) )THEN
1344     *
1345     *           Form  B := alpha*inv( A )*B.
1346     *
1347                 IF( UPPER )THEN
1348                    DO 60, J = 1, N
1349                       IF( ALPHA.NE.ONE )THEN
1350                          DO 30, I = 1, M
1351                             B( I, J ) = ALPHA*B( I, J )
1352        30                CONTINUE
1353                       END IF
1354                       DO 50, K = M, 1, -1
1355                          IF( B( K, J ).NE.ZERO )THEN
1356                             IF( NOUNIT )
1357          \$                     B( K, J ) = B( K, J )/A( K, K )
1358                             DO 40, I = 1, K - 1
1359                                B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
1360        40                   CONTINUE
1361                          END IF
1362        50             CONTINUE
1363        60          CONTINUE
1364                 ELSE
1365                    DO 100, J = 1, N
1366                       IF( ALPHA.NE.ONE )THEN
1367                          DO 70, I = 1, M
1368                             B( I, J ) = ALPHA*B( I, J )
1369        70                CONTINUE
1370                       END IF
1371                       DO 90 K = 1, M
1372                          IF( B( K, J ).NE.ZERO )THEN
1373                             IF( NOUNIT )
1374          \$                     B( K, J ) = B( K, J )/A( K, K )
1375                             DO 80, I = K + 1, M
1376                                B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
1377        80                   CONTINUE
1378                          END IF
1379        90             CONTINUE
1380       100          CONTINUE
1381                 END IF
1382              ELSE
1383     *
1384     *           Form  B := alpha*inv( A' )*B.
1385     *
1386                 IF( UPPER )THEN
1387                    DO 130, J = 1, N
1388                       DO 120, I = 1, M
1389                          TEMP = ALPHA*B( I, J )
1390                          DO 110, K = 1, I - 1
1391                             TEMP = TEMP - A( K, I )*B( K, J )
1392       110                CONTINUE
1393                          IF( NOUNIT )
1394          \$                  TEMP = TEMP/A( I, I )
1395                          B( I, J ) = TEMP
1396       120             CONTINUE
1397       130          CONTINUE
1398                 ELSE
1399                    DO 160, J = 1, N
1400                       DO 150, I = M, 1, -1
1401                          TEMP = ALPHA*B( I, J )
1402                          DO 140, K = I + 1, M
1403                             TEMP = TEMP - A( K, I )*B( K, J )
1404       140                CONTINUE
1405                          IF( NOUNIT )
1406          \$                  TEMP = TEMP/A( I, I )
1407                          B( I, J ) = TEMP
1408       150             CONTINUE
1409       160          CONTINUE
1410                 END IF
1411              END IF
1412           ELSE
1413              IF( LSAME( TRANSA, 'N' ) )THEN
1414     *
1415     *           Form  B := alpha*B*inv( A ).
1416     *
1417                 IF( UPPER )THEN
1418                    DO 210, J = 1, N
1419                       IF( ALPHA.NE.ONE )THEN
1420                          DO 170, I = 1, M
1421                             B( I, J ) = ALPHA*B( I, J )
1422       170                CONTINUE
1423                       END IF
1424                       DO 190, K = 1, J - 1
1425                          IF( A( K, J ).NE.ZERO )THEN
1426                             DO 180, I = 1, M
1427                                B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
1428       180                   CONTINUE
1429                          END IF
1430       190             CONTINUE
1431                       IF( NOUNIT )THEN
1432                          TEMP = ONE/A( J, J )
1433                          DO 200, I = 1, M
1434                             B( I, J ) = TEMP*B( I, J )
1435       200                CONTINUE
1436                       END IF
1437       210          CONTINUE
1438                 ELSE
1439                    DO 260, J = N, 1, -1
1440                       IF( ALPHA.NE.ONE )THEN
1441                          DO 220, I = 1, M
1442                             B( I, J ) = ALPHA*B( I, J )
1443       220                CONTINUE
1444                       END IF
1445                       DO 240, K = J + 1, N
1446                          IF( A( K, J ).NE.ZERO )THEN
1447                             DO 230, I = 1, M
1448                                B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
1449       230                   CONTINUE
1450                          END IF
1451       240             CONTINUE
1452                       IF( NOUNIT )THEN
1453                          TEMP = ONE/A( J, J )
1454                          DO 250, I = 1, M
1455                            B( I, J ) = TEMP*B( I, J )
1456       250                CONTINUE
1457                       END IF
1458       260          CONTINUE
1459                 END IF
1460              ELSE
1461     *
1462     *           Form  B := alpha*B*inv( A' ).
1463     *
1464                 IF( UPPER )THEN
1465                    DO 310, K = N, 1, -1
1466                       IF( NOUNIT )THEN
1467                          TEMP = ONE/A( K, K )
1468                          DO 270, I = 1, M
1469                             B( I, K ) = TEMP*B( I, K )
1470       270                CONTINUE
1471                       END IF
1472                       DO 290, J = 1, K - 1
1473                          IF( A( J, K ).NE.ZERO )THEN
1474                             TEMP = A( J, K )
1475                             DO 280, I = 1, M
1476                                B( I, J ) = B( I, J ) - TEMP*B( I, K )
1477       280                   CONTINUE
1478                          END IF
1479       290             CONTINUE
1480                       IF( ALPHA.NE.ONE )THEN
1481                          DO 300, I = 1, M
1482                             B( I, K ) = ALPHA*B( I, K )
1483       300                CONTINUE
1484                       END IF
1485       310          CONTINUE
1486                 ELSE
1487                    DO 360, K = 1, N
1488                       IF( NOUNIT )THEN
1489                          TEMP = ONE/A( K, K )
1490                          DO 320, I = 1, M
1491                             B( I, K ) = TEMP*B( I, K )
1492       320                CONTINUE
1493                       END IF
1494                       DO 340, J = K + 1, N
1495                          IF( A( J, K ).NE.ZERO )THEN
1496                             TEMP = A( J, K )
1497                             DO 330, I = 1, M
1498                                B( I, J ) = B( I, J ) - TEMP*B( I, K )
1499       330                   CONTINUE
1500                          END IF
1501       340             CONTINUE
1502                       IF( ALPHA.NE.ONE )THEN
1503                          DO 350, I = 1, M
1504                             B( I, K ) = ALPHA*B( I, K )
1505       350                CONTINUE
1506                       END IF
1507       360          CONTINUE
1508                 END IF
1509              END IF
1510           END IF
1511     *
1512           RETURN
1513     *
1514     *     End of DTRSM .
1515     *
1516           END
1517
1518           SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO )
1519     *
1520     *  -- LAPACK routine (version 2.0) --
1521     *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1522     *     Courant Institute, Argonne National Lab, and Rice University
1523     *     June 30, 1992
1524     *
1525     *     .. Scalar Arguments ..
1526           INTEGER            INFO, LDA, M, N
1527     *     ..
1528     *     .. Array Arguments ..
1529           INTEGER            IPIV( * )
1530           REAL               A( LDA, * )
1531     *     ..
1532     *
1533     *  Purpose
1534     *  =======
1535     *
1536     *  DGETF2 computes an LU factorization of a general m-by-n matrix A
1537     *  using partial pivoting with row interchanges.
1538     *
1539     *  The factorization has the form
1540     *     A = P * L * U
1541     *  where P is a permutation matrix, L is lower triangular with unit
1542     *  diagonal elements (lower trapezoidal if m > n), and U is upper
1543     *  triangular (upper trapezoidal if m < n).
1544     *
1545     *  This is the right-looking Level 2 BLAS version of the algorithm.
1546     *
1547     *  Arguments
1548     *  =========
1549     *
1550     *  M       (input) INTEGER
1551     *          The number of rows of the matrix A.  M >= 0.
1552     *
1553     *  N       (input) INTEGER
1554     *          The number of columns of the matrix A.  N >= 0.
1555     *
1556     *  A       (input/output) REAL array, dimension (LDA,N)
1557     *          On entry, the m by n matrix to be factored.
1558     *          On exit, the factors L and U from the factorization
1559     *          A = P*L*U; the unit diagonal elements of L are not stored.
1560     *
1561     *  LDA     (input) INTEGER
1562     *          The leading dimension of the array A.  LDA >= max(1,M).
1563     *
1564     *  IPIV    (output) INTEGER array, dimension (min(M,N))
1565     *          The pivot indices; for 1 <= i <= min(M,N), row i of the
1566     *          matrix was interchanged with row IPIV(i).
1567     *
1568     *  INFO    (output) INTEGER
1569     *          = 0: successful exit
1570     *          < 0: if INFO = -k, the k-th argument had an illegal value
1571     *          > 0: if INFO = k, U(k,k) is exactly zero. The factorization
1572     *               has been completed, but the factor U is exactly
1573     *               singular, and division by zero will occur if it is used
1574     *               to solve a system of equations.
1575     *
1576     *  =====================================================================
1577     *
1578     *     .. Parameters ..
1579           REAL               ONE, ZERO
1580           PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
1581     *     ..
1582     *     .. Local Scalars ..
1583           INTEGER            J, JP
1584     *     ..
1585     *     .. External Functions ..
1586           INTEGER            ISAMAX
1587           EXTERNAL           ISAMAX
1588     *     ..
1589     *     .. External Subroutines ..
1590           EXTERNAL           DGER, DSCAL, DSWAP, XERBLA
1591     *     ..
1592     *     .. Intrinsic Functions ..
1593           INTRINSIC          MAX, MIN
1594     *     ..
1595     *     .. Executable Statements ..
1596     *
1597     *     Test the input parameters.
1598     *
1599           INFO = 0
1600           IF( M.LT.0 ) THEN
1601              INFO = -1
1602           ELSE IF( N.LT.0 ) THEN
1603              INFO = -2
1604           ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
1605              INFO = -4
1606           END IF
1607           IF( INFO.NE.0 ) THEN
1608              CALL XERBLA( 'DGETF2', -INFO )
1609              RETURN
1610           END IF
1611     *
1612     *     Quick return if possible
1613     *
1614           IF( M.EQ.0 .OR. N.EQ.0 )
1615          \$   RETURN
1616     *
1617           DO 10 J = 1, MIN( M, N )
1618     *
1619     *        Find pivot and test for singularity.
1620     *
1621              JP = J - 1 + ISAMAX( M-J+1, A( J, J ), 1 )
1622              IPIV( J ) = JP
1623              IF( A( JP, J ).NE.ZERO ) THEN
1624     *
1625     *           Apply the interchange to columns 1:N.
1626     *
1627                 IF( JP.NE.J )
1628          \$         CALL DSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA )
1629     *
1630     *           Compute elements J+1:M of J-th column.
1631     *
1632                 IF( J.LT.M )
1633          \$         CALL DSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 )
1634     *
1635              ELSE IF( INFO.EQ.0 ) THEN
1636     *
1637                 INFO = J
1638              END IF
1639     *
1640              IF( J.LT.MIN( M, N ) ) THEN
1641     *
1642     *           Update trailing submatrix.
1643     *
1644                 CALL DGER( M-J, N-J, -ONE, A( J+1, J ), 1, A( J, J+1 ), LDA,
1645          \$                 A( J+1, J+1 ), LDA )
1646              END IF
1647        10 CONTINUE
1648           RETURN
1649     *
1650     *     End of DGETF2
1651     *
1652           END
1653           SUBROUTINE DLASWP( N, A, LDA, K1, K2, IPIV, INCX )
1654     *
1655     *  -- LAPACK auxiliary routine (version 2.0) --
1656     *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1657     *     Courant Institute, Argonne National Lab, and Rice University
1658     *     October 31, 1992
1659     *
1660     *     .. Scalar Arguments ..
1661           INTEGER            INCX, K1, K2, LDA, N
1662     *     ..
1663     *     .. Array Arguments ..
1664           INTEGER            IPIV( * )
1665           REAL               A( LDA, * )
1666     *     ..
1667     *
1668     *  Purpose
1669     *  =======
1670     *
1671     *  DLASWP performs a series of row interchanges on the matrix A.
1672     *  One row interchange is initiated for each of rows K1 through K2 of A.
1673     *
1674     *  Arguments
1675     *  =========
1676     *
1677     *  N       (input) INTEGER
1678     *          The number of columns of the matrix A.
1679     *
1680     *  A       (input/output) REAL array, dimension (LDA,N)
1681     *          On entry, the matrix of column dimension N to which the row
1682     *          interchanges will be applied.
1683     *          On exit, the permuted matrix.
1684     *
1685     *  LDA     (input) INTEGER
1686     *          The leading dimension of the array A.
1687     *
1688     *  K1      (input) INTEGER
1689     *          The first element of IPIV for which a row interchange will
1690     *          be done.
1691     *
1692     *  K2      (input) INTEGER
1693     *          The last element of IPIV for which a row interchange will
1694     *          be done.
1695     *
1696     *  IPIV    (input) INTEGER array, dimension (M*abs(INCX))
1697     *          The vector of pivot indices.  Only the elements in positions
1698     *          K1 through K2 of IPIV are accessed.
1699     *          IPIV(K) = L implies rows K and L are to be interchanged.
1700     *
1701     *  INCX    (input) INTEGER
1702     *          The increment between successive values of IPIV.  If IPIV
1703     *          is negative, the pivots are applied in reverse order.
1704     *
1705     * =====================================================================
1706     *
1707     *     .. Local Scalars ..
1708           INTEGER            I, IP, IX
1709     *     ..
1710     *     .. External Subroutines ..
1711           EXTERNAL           DSWAP
1712     *     ..
1713     *     .. Executable Statements ..
1714     *
1715     *     Interchange row I with row IPIV(I) for each of rows K1 through K2.
1716     *
1717           IF( INCX.EQ.0 )
1718          \$   RETURN
1719           IF( INCX.GT.0 ) THEN
1720              IX = K1
1721           ELSE
1722              IX = 1 + ( 1-K2 )*INCX
1723           END IF
1724           IF( INCX.EQ.1 ) THEN
1725              DO 10 I = K1, K2
1726                 IP = IPIV( I )
1727                 IF( IP.NE.I )
1728          \$         CALL DSWAP( N, A( I, 1 ), LDA, A( IP, 1 ), LDA )
1729        10    CONTINUE
1730           ELSE IF( INCX.GT.1 ) THEN
1731              DO 20 I = K1, K2
1732                 IP = IPIV( IX )
1733                 IF( IP.NE.I )
1734          \$         CALL DSWAP( N, A( I, 1 ), LDA, A( IP, 1 ), LDA )
1735                 IX = IX + INCX
1736        20    CONTINUE
1737           ELSE IF( INCX.LT.0 ) THEN
1738              DO 30 I = K2, K1, -1
1739                 IP = IPIV( IX )
1740                 IF( IP.NE.I )
1741          \$         CALL DSWAP( N, A( I, 1 ), LDA, A( IP, 1 ), LDA )
1742                 IX = IX + INCX
1743        30    CONTINUE
1744           END IF
1745     *
1746           RETURN
1747     *
1748     *     End of DLASWP
1749     *
1750           END
1751           subroutine dswap (n,sx,incx,sy,incy)
1752     c
1753     c     interchanges two vectors.
1754     c     uses unrolled loops for increments equal to 1.
1755     c     jack dongarra, linpack, 3/11/78.
1756     c     modified 12/3/93, array(1) declarations changed to array(*)
1757     c
1758           real sx(*),sy(*),stemp
1759           integer i,incx,incy,ix,iy,m,mp1,n
1760     c
1761           if(n.le.0)return
1762           if(incx.eq.1.and.incy.eq.1)go to 20
1763     c
1764     c       code for unequal increments or equal increments not equal
1765     c         to 1
1766     c
1767           ix = 1
1768           iy = 1
1769           if(incx.lt.0)ix = (-n+1)*incx + 1
1770           if(incy.lt.0)iy = (-n+1)*incy + 1
1771           do 10 i = 1,n
1772             stemp = sx(ix)
1773             sx(ix) = sy(iy)
1774             sy(iy) = stemp
1775             ix = ix + incx
1776             iy = iy + incy
1777        10 continue
1778           return
1779     c
1780     c       code for both increments equal to 1
1781     c
1782     c
1783     c       clean-up loop
1784     c
1785        20 m = mod(n,3)
1786           if( m .eq. 0 ) go to 40
1787           do 30 i = 1,m
1788             stemp = sx(i)
1789             sx(i) = sy(i)
1790             sy(i) = stemp
1791        30 continue
1792           if( n .lt. 3 ) return
1793        40 mp1 = m + 1
1794           do 50 i = mp1,n,3
1795             stemp = sx(i)
1796             sx(i) = sy(i)
1797             sy(i) = stemp
1798             stemp = sx(i + 1)
1799             sx(i + 1) = sy(i + 1)
1800             sy(i + 1) = stemp
1801             stemp = sx(i + 2)
1802             sx(i + 2) = sy(i + 2)
1803             sy(i + 2) = stemp
1804        50 continue
1805           return
1806           end
1807           subroutine dscal(n,sa,sx,incx)
1808     c
1809     c     scales a vector by a constant.
1810     c     uses unrolled loops for increment equal to 1.
1811     c     jack dongarra, linpack, 3/11/78.
1812     c     modified 3/93 to return if incx .le. 0.
1813     c     modified 12/3/93, array(1) declarations changed to array(*)
1814     c
1815           real sa,sx(*)
1816           integer i,incx,m,mp1,n,nincx
1817     c
1818           if( n.le.0 .or. incx.le.0 )return
1819           if(incx.eq.1)go to 20
1820     c
1821     c        code for increment not equal to 1
1822     c
1823           nincx = n*incx
1824           do 10 i = 1,nincx,incx
1825             sx(i) = sa*sx(i)
1826        10 continue
1827           return
1828     c
1829     c        code for increment equal to 1
1830     c
1831     c
1832     c        clean-up loop
1833     c
1834        20 m = mod(n,5)
1835           if( m .eq. 0 ) go to 40
1836           do 30 i = 1,m
1837             sx(i) = sa*sx(i)
1838        30 continue
1839           if( n .lt. 5 ) return
1840        40 mp1 = m + 1
1841           do 50 i = mp1,n,5
1842             sx(i) = sa*sx(i)
1843             sx(i + 1) = sa*sx(i + 1)
1844             sx(i + 2) = sa*sx(i + 2)
1845             sx(i + 3) = sa*sx(i + 3)
1846             sx(i + 4) = sa*sx(i + 4)
1847        50 continue
1848           return
1849           end
1850           SUBROUTINE DGER  ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
1851     *     .. Scalar Arguments ..
1852           REAL               ALPHA
1853           INTEGER            INCX, INCY, LDA, M, N
1854     *     .. Array Arguments ..
1855           REAL               A( LDA, * ), X( * ), Y( * )
1856     *     ..
1857     *
1858     *  Purpose
1859     *  =======
1860     *
1861     *  DGER   performs the rank 1 operation
1862     *
1863     *     A := alpha*x*y' + A,
1864     *
1865     *  where alpha is a scalar, x is an m element vector, y is an n element
1866     *  vector and A is an m by n matrix.
1867     *
1868     *  Parameters
1869     *  ==========
1870     *
1871     *  M      - INTEGER.
1872     *           On entry, M specifies the number of rows of the matrix A.
1873     *           M must be at least zero.
1874     *           Unchanged on exit.
1875     *
1876     *  N      - INTEGER.
1877     *           On entry, N specifies the number of columns of the matrix A.
1878     *           N must be at least zero.
1879     *           Unchanged on exit.
1880     *
1881     *  ALPHA  - REAL            .
1882     *           On entry, ALPHA specifies the scalar alpha.
1883     *           Unchanged on exit.
1884     *
1885     *  X      - REAL             array of dimension at least
1886     *           ( 1 + ( m - 1 )*abs( INCX ) ).
1887     *           Before entry, the incremented array X must contain the m
1888     *           element vector x.
1889     *           Unchanged on exit.
1890     *
1891     *  INCX   - INTEGER.
1892     *           On entry, INCX specifies the increment for the elements of
1893     *           X. INCX must not be zero.
1894     *           Unchanged on exit.
1895     *
1896     *  Y      - REAL             array of dimension at least
1897     *           ( 1 + ( n - 1 )*abs( INCY ) ).
1898     *           Before entry, the incremented array Y must contain the n
1899     *           element vector y.
1900     *           Unchanged on exit.
1901     *
1902     *  INCY   - INTEGER.
1903     *           On entry, INCY specifies the increment for the elements of
1904     *           Y. INCY must not be zero.
1905     *           Unchanged on exit.
1906     *
1907     *  A      - REAL             array of DIMENSION ( LDA, n ).
1908     *           Before entry, the leading m by n part of the array A must
1909     *           contain the matrix of coefficients. On exit, A is
1910     *           overwritten by the updated matrix.
1911     *
1912     *  LDA    - INTEGER.
1913     *           On entry, LDA specifies the first dimension of A as declared
1914     *           in the calling (sub) program. LDA must be at least
1915     *           max( 1, m ).
1916     *           Unchanged on exit.
1917     *
1918     *
1919     *  Level 2 Blas routine.
1920     *
1921     *  -- Written on 22-October-1986.
1922     *     Jack Dongarra, Argonne National Lab.
1923     *     Jeremy Du Croz, Nag Central Office.
1924     *     Sven Hammarling, Nag Central Office.
1925     *     Richard Hanson, Sandia National Labs.
1926     *
1927     *
1928     *     .. Parameters ..
1929           REAL               ZERO
1930           PARAMETER        ( ZERO = 0.0E+0 )
1931     *     .. Local Scalars ..
1932           REAL               TEMP
1933           INTEGER            I, INFO, IX, J, JY, KX
1934     *     .. External Subroutines ..
1935           EXTERNAL           XERBLA
1936     *     .. Intrinsic Functions ..
1937           INTRINSIC          MAX
1938     *     ..
1939     *     .. Executable Statements ..
1940     *
1941     *     Test the input parameters.
1942     *
1943           INFO = 0
1944           IF     ( M.LT.0 )THEN
1945              INFO = 1
1946           ELSE IF( N.LT.0 )THEN
1947              INFO = 2
1948           ELSE IF( INCX.EQ.0 )THEN
1949              INFO = 5
1950           ELSE IF( INCY.EQ.0 )THEN
1951              INFO = 7
1952           ELSE IF( LDA.LT.MAX( 1, M ) )THEN
1953              INFO = 9
1954           END IF
1955           IF( INFO.NE.0 )THEN
1956              CALL XERBLA( 'DGER  ', INFO )
1957              RETURN
1958           END IF
1959     *
1960     *     Quick return if possible.
1961     *
1962           IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
1963          \$   RETURN
1964     *
1965     *     Start the operations. In this version the elements of A are
1966     *     accessed sequentially with one pass through A.
1967     *
1968           IF( INCY.GT.0 )THEN
1969              JY = 1
1970           ELSE
1971              JY = 1 - ( N - 1 )*INCY
1972           END IF
1973           IF( INCX.EQ.1 )THEN
1974              DO 20, J = 1, N
1975                 IF( Y( JY ).NE.ZERO )THEN
1976                    TEMP = ALPHA*Y( JY )
1977                    DO 10, I = 1, M
1978                       A( I, J ) = A( I, J ) + X( I )*TEMP
1979        10          CONTINUE
1980                 END IF
1981                 JY = JY + INCY
1982        20    CONTINUE
1983           ELSE
1984              IF( INCX.GT.0 )THEN
1985                 KX = 1
1986              ELSE
1987                 KX = 1 - ( M - 1 )*INCX
1988              END IF
1989              DO 40, J = 1, N
1990                 IF( Y( JY ).NE.ZERO )THEN
1991                    TEMP = ALPHA*Y( JY )
1992                    IX   = KX
1993                    DO 30, I = 1, M
1994                       A( I, J ) = A( I, J ) + X( IX )*TEMP
1995                       IX        = IX        + INCX
1996        30          CONTINUE
1997                 END IF
1998                 JY = JY + INCY
1999        40    CONTINUE
2000           END IF
2001     *
2002           RETURN
2003     *
2004     *     End of DGER  .
2005     *
2006           END
2007           SUBROUTINE DTRTI2( UPLO, DIAG, N, A, LDA, INFO )
2008     *
2009     *  -- LAPACK routine (version 2.0) --
2010     *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
2011     *     Courant Institute, Argonne National Lab, and Rice University
2012     *     February 29, 1992
2013     *
2014     *     .. Scalar Arguments ..
2015           CHARACTER          DIAG, UPLO
2016           INTEGER            INFO, LDA, N
2017     *     ..
2018     *     .. Array Arguments ..
2019           REAL               A( LDA, * )
2020     *     ..
2021     *
2022     *  Purpose
2023     *  =======
2024     *
2025     *  DTRTI2 computes the inverse of a real upper or lower triangular
2026     *  matrix.
2027     *
2028     *  This is the Level 2 BLAS version of the algorithm.
2029     *
2030     *  Arguments
2031     *  =========
2032     *
2033     *  UPLO    (input) CHARACTER*1
2034     *          Specifies whether the matrix A is upper or lower triangular.
2035     *          = 'U':  Upper triangular
2036     *          = 'L':  Lower triangular
2037     *
2038     *  DIAG    (input) CHARACTER*1
2039     *          Specifies whether or not the matrix A is unit triangular.
2040     *          = 'N':  Non-unit triangular
2041     *          = 'U':  Unit triangular
2042     *
2043     *  N       (input) INTEGER
2044     *          The order of the matrix A.  N >= 0.
2045     *
2046     *  A       (input/output) REAL array, dimension (LDA,N)
2047     *          On entry, the triangular matrix A.  If UPLO = 'U', the
2048     *          leading n by n upper triangular part of the array A contains
2049     *          the upper triangular matrix, and the strictly lower
2050     *          triangular part of A is not referenced.  If UPLO = 'L', the
2051     *          leading n by n lower triangular part of the array A contains
2052     *          the lower triangular matrix, and the strictly upper
2053     *          triangular part of A is not referenced.  If DIAG = 'U', the
2054     *          diagonal elements of A are also not referenced and are
2055     *          assumed to be 1.
2056     *
2057     *          On exit, the (triangular) inverse of the original matrix, in
2058     *          the same storage format.
2059     *
2060     *  LDA     (input) INTEGER
2061     *          The leading dimension of the array A.  LDA >= max(1,N).
2062     *
2063     *  INFO    (output) INTEGER
2064     *          = 0: successful exit
2065     *          < 0: if INFO = -k, the k-th argument had an illegal value
2066     *
2067     *  =====================================================================
2068     *
2069     *     .. Parameters ..
2070           REAL               ONE
2071           PARAMETER          ( ONE = 1.0E+0 )
2072     *     ..
2073     *     .. Local Scalars ..
2074           LOGICAL            NOUNIT, UPPER
2075           INTEGER            J
2076           REAL               AJJ
2077     *     ..
2078     *     .. External Functions ..
2079           LOGICAL            LSAME
2080           EXTERNAL           LSAME
2081     *     ..
2082     *     .. External Subroutines ..
2083           EXTERNAL           DSCAL, DTRMV, XERBLA
2084     *     ..
2085     *     .. Intrinsic Functions ..
2086           INTRINSIC          MAX
2087     *     ..
2088     *     .. Executable Statements ..
2089     *
2090     *     Test the input parameters.
2091     *
2092           INFO = 0
2093           UPPER = LSAME( UPLO, 'U' )
2094           NOUNIT = LSAME( DIAG, 'N' )
2095           IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
2096              INFO = -1
2097           ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
2098              INFO = -2
2099           ELSE IF( N.LT.0 ) THEN
2100              INFO = -3
2101           ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
2102              INFO = -5
2103           END IF
2104           IF( INFO.NE.0 ) THEN
2105              CALL XERBLA( 'DTRTI2', -INFO )
2106              RETURN
2107           END IF
2108     *
2109           IF( UPPER ) THEN
2110     *
2111     *        Compute inverse of upper triangular matrix.
2112     *
2113              DO 10 J = 1, N
2114                 IF( NOUNIT ) THEN
2115                    A( J, J ) = ONE / A( J, J )
2116                    AJJ = -A( J, J )
2117                 ELSE
2118                    AJJ = -ONE
2119                 END IF
2120     *
2121     *           Compute elements 1:j-1 of j-th column.
2122     *
2123                 CALL DTRMV( 'Upper', 'No transpose', DIAG, J-1, A, LDA,
2124          \$                  A( 1, J ), 1 )
2125                 CALL DSCAL( J-1, AJJ, A( 1, J ), 1 )
2126        10    CONTINUE
2127           ELSE
2128     *
2129     *        Compute inverse of lower triangular matrix.
2130     *
2131              DO 20 J = N, 1, -1
2132                 IF( NOUNIT ) THEN
2133                    A( J, J ) = ONE / A( J, J )
2134                    AJJ = -A( J, J )
2135                 ELSE
2136                    AJJ = -ONE
2137                 END IF
2138                 IF( J.LT.N ) THEN
2139     *
2140     *              Compute elements j+1:n of j-th column.
2141     *
2142                    CALL DTRMV( 'Lower', 'No transpose', DIAG, N-J,
2143          \$                     A( J+1, J+1 ), LDA, A( J+1, J ), 1 )
2144                    CALL DSCAL( N-J, AJJ, A( J+1, J ), 1 )
2145                 END IF
2146        20    CONTINUE
2147           END IF
2148     *
2149           RETURN
2150     *
2151     *     End of DTRTI2
2152     *
2153           END
2154           SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
2155          \$                   B, LDB )
2156     *     .. Scalar Arguments ..
2157           CHARACTER*1        SIDE, UPLO, TRANSA, DIAG
2158           INTEGER            M, N, LDA, LDB
2159           REAL               ALPHA
2160     *     .. Array Arguments ..
2161           REAL               A( LDA, * ), B( LDB, * )
2162     *     ..
2163     *
2164     *  Purpose
2165     *  =======
2166     *
2167     *  DTRMM  performs one of the matrix-matrix operations
2168     *
2169     *     B := alpha*op( A )*B,   or   B := alpha*B*op( A ),
2170     *
2171     *  where  alpha  is a scalar,  B  is an m by n matrix,  A  is a unit, or
2172     *  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
2173     *
2174     *     op( A ) = A   or   op( A ) = A'.
2175     *
2176     *  Parameters
2177     *  ==========
2178     *
2179     *  SIDE   - CHARACTER*1.
2180     *           On entry,  SIDE specifies whether  op( A ) multiplies B from
2181     *           the left or right as follows:
2182     *
2183     *              SIDE = 'L' or 'l'   B := alpha*op( A )*B.
2184     *
2185     *              SIDE = 'R' or 'r'   B := alpha*B*op( A ).
2186     *
2187     *           Unchanged on exit.
2188     *
2189     *  UPLO   - CHARACTER*1.
2190     *           On entry, UPLO specifies whether the matrix A is an upper or
2191     *           lower triangular matrix as follows:
2192     *
2193     *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
2194     *
2195     *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
2196     *
2197     *           Unchanged on exit.
2198     *
2199     *  TRANSA - CHARACTER*1.
2200     *           On entry, TRANSA specifies the form of op( A ) to be used in
2201     *           the matrix multiplication as follows:
2202     *
2203     *              TRANSA = 'N' or 'n'   op( A ) = A.
2204     *
2205     *              TRANSA = 'T' or 't'   op( A ) = A'.
2206     *
2207     *              TRANSA = 'C' or 'c'   op( A ) = A'.
2208     *
2209     *           Unchanged on exit.
2210     *
2211     *  DIAG   - CHARACTER*1.
2212     *           On entry, DIAG specifies whether or not A is unit triangular
2213     *           as follows:
2214     *
2215     *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
2216     *
2217     *              DIAG = 'N' or 'n'   A is not assumed to be unit
2218     *                                  triangular.
2219     *
2220     *           Unchanged on exit.
2221     *
2222     *  M      - INTEGER.
2223     *           On entry, M specifies the number of rows of B. M must be at
2224     *           least zero.
2225     *           Unchanged on exit.
2226     *
2227     *  N      - INTEGER.
2228     *           On entry, N specifies the number of columns of B.  N must be
2229     *           at least zero.
2230     *           Unchanged on exit.
2231     *
2232     *  ALPHA  - REAL            .
2233     *           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
2234     *           zero then  A is not referenced and  B need not be set before
2235     *           entry.
2236     *           Unchanged on exit.
2237     *
2238     *  A      - REAL             array of DIMENSION ( LDA, k ), where k is m
2239     *           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
2240     *           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
2241     *           upper triangular part of the array  A must contain the upper
2242     *           triangular matrix  and the strictly lower triangular part of
2243     *           A is not referenced.
2244     *           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
2245     *           lower triangular part of the array  A must contain the lower
2246     *           triangular matrix  and the strictly upper triangular part of
2247     *           A is not referenced.
2248     *           Note that when  DIAG = 'U' or 'u',  the diagonal elements of
2249     *           A  are not referenced either,  but are assumed to be  unity.
2250     *           Unchanged on exit.
2251     *
2252     *  LDA    - INTEGER.
2253     *           On entry, LDA specifies the first dimension of A as declared
2254     *           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
2255     *           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
2256     *           then LDA must be at least max( 1, n ).
2257     *           Unchanged on exit.
2258     *
2259     *  B      - REAL             array of DIMENSION ( LDB, n ).
2260     *           Before entry,  the leading  m by n part of the array  B must
2261     *           contain the matrix  B,  and  on exit  is overwritten  by the
2262     *           transformed matrix.
2263     *
2264     *  LDB    - INTEGER.
2265     *           On entry, LDB specifies the first dimension of B as declared
2266     *           in  the  calling  (sub)  program.   LDB  must  be  at  least
2267     *           max( 1, m ).
2268     *           Unchanged on exit.
2269     *
2270     *
2271     *  Level 3 Blas routine.
2272     *
2273     *  -- Written on 8-February-1989.
2274     *     Jack Dongarra, Argonne National Laboratory.
2275     *     Iain Duff, AERE Harwell.
2276     *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
2277     *     Sven Hammarling, Numerical Algorithms Group Ltd.
2278     *
2279     *
2280     *     .. External Functions ..
2281           LOGICAL            LSAME
2282           EXTERNAL           LSAME
2283     *     .. External Subroutines ..
2284           EXTERNAL           XERBLA
2285     *     .. Intrinsic Functions ..
2286           INTRINSIC          MAX
2287     *     .. Local Scalars ..
2288           LOGICAL            LSIDE, NOUNIT, UPPER
2289           INTEGER            I, INFO, J, K, NROWA
2290           REAL               TEMP
2291     *     .. Parameters ..
2292           REAL               ONE         , ZERO
2293           PARAMETER        ( ONE = 1.0E+0, ZERO = 0.0E+0 )
2294     *     ..
2295     *     .. Executable Statements ..
2296     *
2297     *     Test the input parameters.
2298     *
2299           LSIDE  = LSAME( SIDE  , 'L' )
2300           IF( LSIDE )THEN
2301              NROWA = M
2302           ELSE
2303              NROWA = N
2304           END IF
2305           NOUNIT = LSAME( DIAG  , 'N' )
2306           UPPER  = LSAME( UPLO  , 'U' )
2307     *
2308           INFO   = 0
2309           IF(      ( .NOT.LSIDE                ).AND.
2310          \$         ( .NOT.LSAME( SIDE  , 'R' ) )      )THEN
2311              INFO = 1
2312           ELSE IF( ( .NOT.UPPER                ).AND.
2313          \$         ( .NOT.LSAME( UPLO  , 'L' ) )      )THEN
2314              INFO = 2
2315           ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
2316          \$         ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
2317          \$         ( .NOT.LSAME( TRANSA, 'C' ) )      )THEN
2318              INFO = 3
2319           ELSE IF( ( .NOT.LSAME( DIAG  , 'U' ) ).AND.
2320          \$         ( .NOT.LSAME( DIAG  , 'N' ) )      )THEN
2321              INFO = 4
2322           ELSE IF( M  .LT.0               )THEN
2323              INFO = 5
2324           ELSE IF( N  .LT.0               )THEN
2325              INFO = 6
2326           ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
2327              INFO = 9
2328           ELSE IF( LDB.LT.MAX( 1, M     ) )THEN
2329              INFO = 11
2330           END IF
2331           IF( INFO.NE.0 )THEN
2332              CALL XERBLA( 'DTRMM ', INFO )
2333              RETURN
2334           END IF
2335     *
2336     *     Quick return if possible.
2337     *
2338           IF( N.EQ.0 )
2339          \$   RETURN
2340     *
2341     *     And when  alpha.eq.zero.
2342     *
2343           IF( ALPHA.EQ.ZERO )THEN
2344              DO 20, J = 1, N
2345                 DO 10, I = 1, M
2346                    B( I, J ) = ZERO
2347        10       CONTINUE
2348        20    CONTINUE
2349              RETURN
2350           END IF
2351     *
2352     *     Start the operations.
2353     *
2354           IF( LSIDE )THEN
2355              IF( LSAME( TRANSA, 'N' ) )THEN
2356     *
2357     *           Form  B := alpha*A*B.
2358     *
2359                 IF( UPPER )THEN
2360                    DO 50, J = 1, N
2361                       DO 40, K = 1, M
2362                          IF( B( K, J ).NE.ZERO )THEN
2363                             TEMP = ALPHA*B( K, J )
2364                             DO 30, I = 1, K - 1
2365                                B( I, J ) = B( I, J ) + TEMP*A( I, K )
2366        30                   CONTINUE
2367                             IF( NOUNIT )
2368          \$                     TEMP = TEMP*A( K, K )
2369                             B( K, J ) = TEMP
2370                          END IF
2371        40             CONTINUE
2372        50          CONTINUE
2373                 ELSE
2374                    DO 80, J = 1, N
2375                       DO 70 K = M, 1, -1
2376                          IF( B( K, J ).NE.ZERO )THEN
2377                             TEMP      = ALPHA*B( K, J )
2378                             B( K, J ) = TEMP
2379                             IF( NOUNIT )
2380          \$                     B( K, J ) = B( K, J )*A( K, K )
2381                             DO 60, I = K + 1, M
2382                                B( I, J ) = B( I, J ) + TEMP*A( I, K )
2383        60                   CONTINUE
2384                          END IF
2385        70             CONTINUE
2386        80          CONTINUE
2387                 END IF
2388              ELSE
2389     *
2390     *           Form  B := alpha*A'*B.
2391     *
2392                 IF( UPPER )THEN
2393                    DO 110, J = 1, N
2394                       DO 100, I = M, 1, -1
2395                          TEMP = B( I, J )
2396                          IF( NOUNIT )
2397          \$                  TEMP = TEMP*A( I, I )
2398                          DO 90, K = 1, I - 1
2399                             TEMP = TEMP + A( K, I )*B( K, J )
2400        90                CONTINUE
2401                          B( I, J ) = ALPHA*TEMP
2402       100             CONTINUE
2403       110          CONTINUE
2404                 ELSE
2405                    DO 140, J = 1, N
2406                       DO 130, I = 1, M
2407                          TEMP = B( I, J )
2408                          IF( NOUNIT )
2409          \$                  TEMP = TEMP*A( I, I )
2410                          DO 120, K = I + 1, M
2411                             TEMP = TEMP + A( K, I )*B( K, J )
2412       120                CONTINUE
2413                          B( I, J ) = ALPHA*TEMP
2414       130             CONTINUE
2415       140          CONTINUE
2416                 END IF
2417              END IF
2418           ELSE
2419              IF( LSAME( TRANSA, 'N' ) )THEN
2420     *
2421     *           Form  B := alpha*B*A.
2422     *
2423                 IF( UPPER )THEN
2424                    DO 180, J = N, 1, -1
2425                       TEMP = ALPHA
2426                       IF( NOUNIT )
2427          \$               TEMP = TEMP*A( J, J )
2428                       DO 150, I = 1, M
2429                          B( I, J ) = TEMP*B( I, J )
2430       150             CONTINUE
2431                       DO 170, K = 1, J - 1
2432                          IF( A( K, J ).NE.ZERO )THEN
2433                             TEMP = ALPHA*A( K, J )
2434                             DO 160, I = 1, M
2435                                B( I, J ) = B( I, J ) + TEMP*B( I, K )
2436       160                   CONTINUE
2437                          END IF
2438       170             CONTINUE
2439       180          CONTINUE
2440                 ELSE
2441                    DO 220, J = 1, N
2442                       TEMP = ALPHA
2443                       IF( NOUNIT )
2444          \$               TEMP = TEMP*A( J, J )
2445                       DO 190, I = 1, M
2446                          B( I, J ) = TEMP*B( I, J )
2447       190             CONTINUE
2448                       DO 210, K = J + 1, N
2449                          IF( A( K, J ).NE.ZERO )THEN
2450                             TEMP = ALPHA*A( K, J )
2451                             DO 200, I = 1, M
2452                                B( I, J ) = B( I, J ) + TEMP*B( I, K )
2453       200                   CONTINUE
2454                          END IF
2455       210             CONTINUE
2456       220          CONTINUE
2457                 END IF
2458              ELSE
2459     *
2460     *           Form  B := alpha*B*A'.
2461     *
2462                 IF( UPPER )THEN
2463                    DO 260, K = 1, N
2464                       DO 240, J = 1, K - 1
2465                          IF( A( J, K ).NE.ZERO )THEN
2466                             TEMP = ALPHA*A( J, K )
2467                             DO 230, I = 1, M
2468                                B( I, J ) = B( I, J ) + TEMP*B( I, K )
2469       230                   CONTINUE
2470                          END IF
2471       240             CONTINUE
2472                       TEMP = ALPHA
2473                       IF( NOUNIT )
2474          \$               TEMP = TEMP*A( K, K )
2475                       IF( TEMP.NE.ONE )THEN
2476                          DO 250, I = 1, M
2477                             B( I, K ) = TEMP*B( I, K )
2478       250                CONTINUE
2479                       END IF
2480       260          CONTINUE
2481                 ELSE
2482                    DO 300, K = N, 1, -1
2483                       DO 280, J = K + 1, N
2484                          IF( A( J, K ).NE.ZERO )THEN
2485                             TEMP = ALPHA*A( J, K )
2486                             DO 270, I = 1, M
2487                                B( I, J ) = B( I, J ) + TEMP*B( I, K )
2488       270                   CONTINUE
2489                          END IF
2490       280             CONTINUE
2491                       TEMP = ALPHA
2492                       IF( NOUNIT )
2493          \$               TEMP = TEMP*A( K, K )
2494                       IF( TEMP.NE.ONE )THEN
2495                          DO 290, I = 1, M
2496                             B( I, K ) = TEMP*B( I, K )
2497       290                CONTINUE
2498                       END IF
2499       300          CONTINUE
2500                 END IF
2501              END IF
2502           END IF
2503     *
2504           RETURN
2505     *
2506     *     End of DTRMM .
2507     *
2508           END
2509           SUBROUTINE DTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
2510     *     .. Scalar Arguments ..
2511           INTEGER            INCX, LDA, N
2512           CHARACTER*1        DIAG, TRANS, UPLO
2513     *     .. Array Arguments ..
2514           REAL               A( LDA, * ), X( * )
2515     *     ..
2516     *
2517     *  Purpose
2518     *  =======
2519     *
2520     *  DTRMV  performs one of the matrix-vector operations
2521     *
2522     *     x := A*x,   or   x := A'*x,
2523     *
2524     *  where x is an n element vector and  A is an n by n unit, or non-unit,
2525     *  upper or lower triangular matrix.
2526     *
2527     *  Parameters
2528     *  ==========
2529     *
2530     *  UPLO   - CHARACTER*1.
2531     *           On entry, UPLO specifies whether the matrix is an upper or
2532     *           lower triangular matrix as follows:
2533     *
2534     *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
2535     *
2536     *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
2537     *
2538     *           Unchanged on exit.
2539     *
2540     *  TRANS  - CHARACTER*1.
2541     *           On entry, TRANS specifies the operation to be performed as
2542     *           follows:
2543     *
2544     *              TRANS = 'N' or 'n'   x := A*x.
2545     *
2546     *              TRANS = 'T' or 't'   x := A'*x.
2547     *
2548     *              TRANS = 'C' or 'c'   x := A'*x.
2549     *
2550     *           Unchanged on exit.
2551     *
2552     *  DIAG   - CHARACTER*1.
2553     *           On entry, DIAG specifies whether or not A is unit
2554     *           triangular as follows:
2555     *
2556     *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
2557     *
2558     *              DIAG = 'N' or 'n'   A is not assumed to be unit
2559     *                                  triangular.
2560     *
2561     *           Unchanged on exit.
2562     *
2563     *  N      - INTEGER.
2564     *           On entry, N specifies the order of the matrix A.
2565     *           N must be at least zero.
2566     *           Unchanged on exit.
2567     *
2568     *  A      - REAL             array of DIMENSION ( LDA, n ).
2569     *           Before entry with  UPLO = 'U' or 'u', the leading n by n
2570     *           upper triangular part of the array A must contain the upper
2571     *           triangular matrix and the strictly lower triangular part of
2572     *           A is not referenced.
2573     *           Before entry with UPLO = 'L' or 'l', the leading n by n
2574     *           lower triangular part of the array A must contain the lower
2575     *           triangular matrix and the strictly upper triangular part of
2576     *           A is not referenced.
2577     *           Note that when  DIAG = 'U' or 'u', the diagonal elements of
2578     *           A are not referenced either, but are assumed to be unity.
2579     *           Unchanged on exit.
2580     *
2581     *  LDA    - INTEGER.
2582     *           On entry, LDA specifies the first dimension of A as declared
2583     *           in the calling (sub) program. LDA must be at least
2584     *           max( 1, n ).
2585     *           Unchanged on exit.
2586     *
2587     *  X      - REAL             array of dimension at least
2588     *           ( 1 + ( n - 1 )*abs( INCX ) ).
2589     *           Before entry, the incremented array X must contain the n
2590     *           element vector x. On exit, X is overwritten with the
2591     *           tranformed vector x.
2592     *
2593     *  INCX   - INTEGER.
2594     *           On entry, INCX specifies the increment for the elements of
2595     *           X. INCX must not be zero.
2596     *           Unchanged on exit.
2597     *
2598     *
2599     *  Level 2 Blas routine.
2600     *
2601     *  -- Written on 22-October-1986.
2602     *     Jack Dongarra, Argonne National Lab.
2603     *     Jeremy Du Croz, Nag Central Office.
2604     *     Sven Hammarling, Nag Central Office.
2605     *     Richard Hanson, Sandia National Labs.
2606     *
2607     *
2608     *     .. Parameters ..
2609           REAL               ZERO
2610           PARAMETER        ( ZERO = 0.0E+0 )
2611     *     .. Local Scalars ..
2612           REAL               TEMP
2613           INTEGER            I, INFO, IX, J, JX, KX
2614           LOGICAL            NOUNIT
2615     *     .. External Functions ..
2616           LOGICAL            LSAME
2617           EXTERNAL           LSAME
2618     *     .. External Subroutines ..
2619           EXTERNAL           XERBLA
2620     *     .. Intrinsic Functions ..
2621           INTRINSIC          MAX
2622     *     ..
2623     *     .. Executable Statements ..
2624     *
2625     *     Test the input parameters.
2626     *
2627           INFO = 0
2628           IF     ( .NOT.LSAME( UPLO , 'U' ).AND.
2629          \$         .NOT.LSAME( UPLO , 'L' )      )THEN
2630              INFO = 1
2631           ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
2632          \$         .NOT.LSAME( TRANS, 'T' ).AND.
2633          \$         .NOT.LSAME( TRANS, 'C' )      )THEN
2634              INFO = 2
2635           ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
2636          \$         .NOT.LSAME( DIAG , 'N' )      )THEN
2637              INFO = 3
2638           ELSE IF( N.LT.0 )THEN
2639              INFO = 4
2640           ELSE IF( LDA.LT.MAX( 1, N ) )THEN
2641              INFO = 6
2642           ELSE IF( INCX.EQ.0 )THEN
2643              INFO = 8
2644           END IF
2645           IF( INFO.NE.0 )THEN
2646              CALL XERBLA( 'DTRMV ', INFO )
2647              RETURN
2648           END IF
2649     *
2650     *     Quick return if possible.
2651     *
2652           IF( N.EQ.0 )
2653          \$   RETURN
2654     *
2655           NOUNIT = LSAME( DIAG, 'N' )
2656     *
2657     *     Set up the start point in X if the increment is not unity. This
2658     *     will be  ( N - 1 )*INCX  too small for descending loops.
2659     *
2660           IF( INCX.LE.0 )THEN
2661              KX = 1 - ( N - 1 )*INCX
2662           ELSE IF( INCX.NE.1 )THEN
2663              KX = 1
2664           END IF
2665     *
2666     *     Start the operations. In this version the elements of A are
2667     *     accessed sequentially with one pass through A.
2668     *
2669           IF( LSAME( TRANS, 'N' ) )THEN
2670     *
2671     *        Form  x := A*x.
2672     *
2673              IF( LSAME( UPLO, 'U' ) )THEN
2674                 IF( INCX.EQ.1 )THEN
2675                    DO 20, J = 1, N
2676                       IF( X( J ).NE.ZERO )THEN
2677                          TEMP = X( J )
2678                          DO 10, I = 1, J - 1
2679                             X( I ) = X( I ) + TEMP*A( I, J )
2680        10                CONTINUE
2681                          IF( NOUNIT )
2682          \$                  X( J ) = X( J )*A( J, J )
2683                       END IF
2684        20          CONTINUE
2685                 ELSE
2686                    JX = KX
2687                    DO 40, J = 1, N
2688                       IF( X( JX ).NE.ZERO )THEN
2689                          TEMP = X( JX )
2690                          IX   = KX
2691                          DO 30, I = 1, J - 1
2692                             X( IX ) = X( IX ) + TEMP*A( I, J )
2693                             IX      = IX      + INCX
2694        30                CONTINUE
2695                          IF( NOUNIT )
2696          \$                  X( JX ) = X( JX )*A( J, J )
2697                       END IF
2698                       JX = JX + INCX
2699        40          CONTINUE
2700                 END IF
2701              ELSE
2702                 IF( INCX.EQ.1 )THEN
2703                    DO 60, J = N, 1, -1
2704                       IF( X( J ).NE.ZERO )THEN
2705                          TEMP = X( J )
2706                          DO 50, I = N, J + 1, -1
2707                             X( I ) = X( I ) + TEMP*A( I, J )
2708        50                CONTINUE
2709                          IF( NOUNIT )
2710          \$                  X( J ) = X( J )*A( J, J )
2711                       END IF
2712        60          CONTINUE
2713                 ELSE
2714                    KX = KX + ( N - 1 )*INCX
2715                    JX = KX
2716                    DO 80, J = N, 1, -1
2717                       IF( X( JX ).NE.ZERO )THEN
2718                          TEMP = X( JX )
2719                          IX   = KX
2720                          DO 70, I = N, J + 1, -1
2721                             X( IX ) = X( IX ) + TEMP*A( I, J )
2722                             IX      = IX      - INCX
2723        70                CONTINUE
2724                          IF( NOUNIT )
2725          \$                  X( JX ) = X( JX )*A( J, J )
2726                       END IF
2727                       JX = JX - INCX
2728        80          CONTINUE
2729                 END IF
2730              END IF
2731           ELSE
2732     *
2733     *        Form  x := A'*x.
2734     *
2735              IF( LSAME( UPLO, 'U' ) )THEN
2736                 IF( INCX.EQ.1 )THEN
2737                    DO 100, J = N, 1, -1
2738                       TEMP = X( J )
2739                       IF( NOUNIT )
2740          \$               TEMP = TEMP*A( J, J )
2741                       DO 90, I = J - 1, 1, -1
2742                          TEMP = TEMP + A( I, J )*X( I )
2743        90             CONTINUE
2744                       X( J ) = TEMP
2745       100          CONTINUE
2746                 ELSE
2747                    JX = KX + ( N - 1 )*INCX
2748                    DO 120, J = N, 1, -1
2749                       TEMP = X( JX )
2750                       IX   = JX
2751                       IF( NOUNIT )
2752          \$               TEMP = TEMP*A( J, J )
2753                       DO 110, I = J - 1, 1, -1
2754                          IX   = IX   - INCX
2755                          TEMP = TEMP + A( I, J )*X( IX )
2756       110             CONTINUE
2757                       X( JX ) = TEMP
2758                       JX      = JX   - INCX
2759       120          CONTINUE
2760                 END IF
2761              ELSE
2762                 IF( INCX.EQ.1 )THEN
2763                    DO 140, J = 1, N
2764                       TEMP = X( J )
2765                       IF( NOUNIT )
2766          \$               TEMP = TEMP*A( J, J )
2767                       DO 130, I = J + 1, N
2768                          TEMP = TEMP + A( I, J )*X( I )
2769       130             CONTINUE
2770                       X( J ) = TEMP
2771       140          CONTINUE
2772                 ELSE
2773                    JX = KX
2774                    DO 160, J = 1, N
2775                       TEMP = X( JX )
2776                       IX   = JX
2777                       IF( NOUNIT )
2778          \$               TEMP = TEMP*A( J, J )
2779                       DO 150, I = J + 1, N
2780                          IX   = IX   + INCX
2781                          TEMP = TEMP + A( I, J )*X( IX )
2782       150             CONTINUE
2783                       X( JX ) = TEMP
2784                       JX      = JX   + INCX
2785       160          CONTINUE
2786                 END IF
2787              END IF
2788           END IF
2789     *
2790           RETURN
2791     *
2792     *     End of DTRMV .
2793     *
2794           END
2795
2796
2797
2798           INTEGER          FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3,
2799          \$                 N4 )
2800     *
2801     *  -- LAPACK auxiliary routine (version 2.0) --
2802     *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
2803     *     Courant Institute, Argonne National Lab, and Rice University
2804     *     September 30, 1994
2805     *
2806     *     .. Scalar Arguments ..
2807           CHARACTER*( * )    NAME, OPTS
2808           INTEGER            ISPEC, N1, N2, N3, N4
2809     *     ..
2810     *
2811     *  Purpose
2812     *  =======
2813     *
2814     *  ILAENV is called from the LAPACK routines to choose problem-dependent
2815     *  parameters for the local environment.  See ISPEC for a description of
2816     *  the parameters.
2817     *
2818     *  This version provides a set of parameters which should give good,
2819     *  but not optimal, performance on many of the currently available
2820     *  computers.  Users are encouraged to modify this subroutine to set
2821     *  the tuning parameters for their particular machine using the option
2822     *  and problem size information in the arguments.
2823     *
2824     *  This routine will not function correctly if it is converted to all
2825     *  lower case.  Converting it to all upper case is allowed.
2826     *
2827     *  Arguments
2828     *  =========
2829     *
2830     *  ISPEC   (input) INTEGER
2831     *          Specifies the parameter to be returned as the value of
2832     *          ILAENV.
2833     *          = 1: the optimal blocksize; if this value is 1, an unblocked
2834     *               algorithm will give the best performance.
2835     *          = 2: the minimum block size for which the block routine
2836     *               should be used; if the usable block size is less than
2837     *               this value, an unblocked routine should be used.
2838     *          = 3: the crossover point (in a block routine, for N less
2839     *               than this value, an unblocked routine should be used)
2840     *          = 4: the number of shifts, used in the nonsymmetric
2841     *               eigenvalue routines
2842     *          = 5: the minimum column dimension for blocking to be used;
2843     *               rectangular blocks must have dimension at least k by m,
2844     *               where k is given by ILAENV(2,...) and m by ILAENV(5,...)
2845     *          = 6: the crossover point for the SVD (when reducing an m by n
2846     *               matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds
2847     *               this value, a QR factorization is used first to reduce
2848     *               the matrix to a triangular form.)
2849     *          = 7: the number of processors
2850     *          = 8: the crossover point for the multishift QR and QZ methods
2851     *               for nonsymmetric eigenvalue problems.
2852     *
2853     *  NAME    (input) CHARACTER*(*)
2854     *          The name of the calling subroutine, in either upper case or
2855     *          lower case.
2856     *
2857     *  OPTS    (input) CHARACTER*(*)
2858     *          The character options to the subroutine NAME, concatenated
2859     *          into a single character string.  For example, UPLO = 'U',
2860     *          TRANS = 'T', and DIAG = 'N' for a triangular routine would
2861     *          be specified as OPTS = 'UTN'.
2862     *
2863     *  N1      (input) INTEGER
2864     *  N2      (input) INTEGER
2865     *  N3      (input) INTEGER
2866     *  N4      (input) INTEGER
2867     *          Problem dimensions for the subroutine NAME; these may not all
2868     *          be required.
2869     *
2870     * (ILAENV) (output) INTEGER
2871     *          >= 0: the value of the parameter specified by ISPEC
2872     *          < 0:  if ILAENV = -k, the k-th argument had an illegal value.
2873     *
2874     *  Further Details
2875     *  ===============
2876     *
2877     *  The following conventions have been used when calling ILAENV from the
2878     *  LAPACK routines:
2879     *  1)  OPTS is a concatenation of all of the character options to
2880     *      subroutine NAME, in the same order that they appear in the
2881     *      argument list for NAME, even if they are not used in determining
2882     *      the value of the parameter specified by ISPEC.
2883     *  2)  The problem dimensions N1, N2, N3, N4 are specified in the order
2884     *      that they appear in the argument list for NAME.  N1 is used
2885     *      first, N2 second, and so on, and unused problem dimensions are
2886     *      passed a value of -1.
2887     *  3)  The parameter value returned by ILAENV is checked for validity in
2888     *      the calling subroutine.  For example, ILAENV is used to retrieve
2889     *      the optimal blocksize for DTRTRI as follows:
2890     *
2891     *      NB = ILAENV( 1, 'DTRTRI', UPLO // DIAG, N, -1, -1, -1 )
2892     *      IF( NB.LE.1 ) NB = MAX( 1, N )
2893     *
2894     *  =====================================================================
2895     *
2896     *     .. Local Scalars ..
2897           LOGICAL            CNAME, SNAME
2898           CHARACTER*1        C1
2899           CHARACTER*2        C2, C4
2900           CHARACTER*3        C3
2901           CHARACTER*6        SUBNAM
2902           INTEGER            I, IC, IZ, NB, NBMIN, NX
2903     *     ..
2904     *     .. Intrinsic Functions ..
2905           INTRINSIC          CHAR, ICHAR, INT, MIN, REAL
2906     *     ..
2907     *     .. Executable Statements ..
2908     *
2909           GO TO ( 100, 100, 100, 400, 500, 600, 700, 800 ) ISPEC
2910     *
2911     *     Invalid value for ISPEC
2912     *
2913           ILAENV = -1
2914           RETURN
2915     *
2916       100 CONTINUE
2917     *
2918     *     Convert NAME to upper case if the first character is lower case.
2919     *
2920           ILAENV = 1
2921           SUBNAM = NAME
2922           IC = MOVA2I( SUBNAM( 1:1 ) )
2923           IZ = MOVA2I( 'Z' )
2924           IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN
2925     *
2926     *        ASCII character set
2927     *
2928              IF( IC.GE.97 .AND. IC.LE.122 ) THEN
2929                 SUBNAM( 1:1 ) = CHAR( IC-32 )
2930                 DO 10 I = 2, 6
2931                    IC = MOVA2I( SUBNAM( I:I ) )
2932                    IF( IC.GE.97 .AND. IC.LE.122 )
2933          \$            SUBNAM( I:I ) = CHAR( IC-32 )
2934        10       CONTINUE
2935              END IF
2936     *
2937           ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN
2938     *
2939     *        EBCDIC character set
2940     *
2941              IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR.
2942          \$       ( IC.GE.145 .AND. IC.LE.153 ) .OR.
2943          \$       ( IC.GE.162 .AND. IC.LE.169 ) ) THEN
2944                 SUBNAM( 1:1 ) = CHAR( IC+64 )
2945                 DO 20 I = 2, 6
2946                    IC = MOVA2I( SUBNAM( I:I ) )
2947                    IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR.
2948          \$             ( IC.GE.145 .AND. IC.LE.153 ) .OR.
2949          \$             ( IC.GE.162 .AND. IC.LE.169 ) )
2950          \$            SUBNAM( I:I ) = CHAR( IC+64 )
2951        20       CONTINUE
2952              END IF
2953     *
2954           ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN
2955     *
2956     *        Prime machines:  ASCII+128
2957     *
2958              IF( IC.GE.225 .AND. IC.LE.250 ) THEN
2959                 SUBNAM( 1:1 ) = CHAR( IC-32 )
2960                 DO 30 I = 2, 6
2961                    IC = MOVA2I( SUBNAM( I:I ) )
2962                    IF( IC.GE.225 .AND. IC.LE.250 )
2963          \$            SUBNAM( I:I ) = CHAR( IC-32 )
2964        30       CONTINUE
2965              END IF
2966           END IF
2967     *
2968           C1 = SUBNAM( 1:1 )
2969           SNAME = C1.EQ.'S' .OR. C1.EQ.'D'
2970           CNAME = C1.EQ.'C' .OR. C1.EQ.'Z'
2971           IF( .NOT.( CNAME .OR. SNAME ) )
2972          \$   RETURN
2973           C2 = SUBNAM( 2:3 )
2974           C3 = SUBNAM( 4:6 )
2975           C4 = C3( 2:3 )
2976     *
2977           GO TO ( 110, 200, 300 ) ISPEC
2978     *
2979       110 CONTINUE
2980     *
2981     *     ISPEC = 1:  block size
2982     *
2983     *     In these examples, separate code is provided for setting NB for
2984     *     real and complex.  We assume that NB will take the same value in
2985     *     single or double precision.
2986     *
2987           NB = 1
2988     *
2989           IF( C2.EQ.'GE' ) THEN
2990              IF( C3.EQ.'TRF' ) THEN
2991                 IF( SNAME ) THEN
2992                    NB = 64
2993                 ELSE
2994                    NB = 64
2995                 END IF
2996              ELSE IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
2997          \$            C3.EQ.'QLF' ) THEN
2998                 IF( SNAME ) THEN
2999                    NB = 32
3000                 ELSE
3001                    NB = 32
3002                 END IF
3003              ELSE IF( C3.EQ.'HRD' ) THEN
3004                 IF( SNAME ) THEN
3005                    NB = 32
3006                 ELSE
3007                    NB = 32
3008                 END IF
3009              ELSE IF( C3.EQ.'BRD' ) THEN
3010                 IF( SNAME ) THEN
3011                    NB = 32
3012                 ELSE
3013                    NB = 32
3014                 END IF
3015              ELSE IF( C3.EQ.'TRI' ) THEN
3016                 IF( SNAME ) THEN
3017                    NB = 64
3018                 ELSE
3019                    NB = 64
3020                 END IF
3021              END IF
3022           ELSE IF( C2.EQ.'PO' ) THEN
3023              IF( C3.EQ.'TRF' ) THEN
3024                 IF( SNAME ) THEN
3025                    NB = 64
3026                 ELSE
3027                    NB = 64
3028                 END IF
3029              END IF
3030           ELSE IF( C2.EQ.'SY' ) THEN
3031              IF( C3.EQ.'TRF' ) THEN
3032                 IF( SNAME ) THEN
3033                    NB = 64
3034                 ELSE
3035                    NB = 64
3036                 END IF
3037              ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN
3038                 NB = 1
3039              ELSE IF( SNAME .AND. C3.EQ.'GST' ) THEN
3040                 NB = 64
3041              END IF
3042           ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
3043              IF( C3.EQ.'TRF' ) THEN
3044                 NB = 64
3045              ELSE IF( C3.EQ.'TRD' ) THEN
3046                 NB = 1
3047              ELSE IF( C3.EQ.'GST' ) THEN
3048                 NB = 64
3049              END IF
3050           ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
3051              IF( C3( 1:1 ).EQ.'G' ) THEN
3052                 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
3053          \$          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
3054          \$          C4.EQ.'BR' ) THEN
3055                    NB = 32
3056                 END IF
3057              ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
3058                 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
3059          \$          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
3060          \$          C4.EQ.'BR' ) THEN
3061                    NB = 32
3062                 END IF
3063              END IF
3064           ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
3065              IF( C3( 1:1 ).EQ.'G' ) THEN
3066                 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
3067          \$          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
3068          \$          C4.EQ.'BR' ) THEN
3069                    NB = 32
3070                 END IF
3071              ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
3072                 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
3073          \$          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
3074          \$          C4.EQ.'BR' ) THEN
3075                    NB = 32
3076                 END IF
3077              END IF
3078           ELSE IF( C2.EQ.'GB' ) THEN
3079              IF( C3.EQ.'TRF' ) THEN
3080                 IF( SNAME ) THEN
3081                    IF( N4.LE.64 ) THEN
3082                       NB = 1
3083                    ELSE
3084                       NB = 32
3085                    END IF
3086                 ELSE
3087                    IF( N4.LE.64 ) THEN
3088                       NB = 1
3089                    ELSE
3090                       NB = 32
3091                    END IF
3092                 END IF
3093              END IF
3094           ELSE IF( C2.EQ.'PB' ) THEN
3095              IF( C3.EQ.'TRF' ) THEN
3096                 IF( SNAME ) THEN
3097                    IF( N2.LE.64 ) THEN
3098                       NB = 1
3099                    ELSE
3100                       NB = 32
3101                    END IF
3102                 ELSE
3103                    IF( N2.LE.64 ) THEN
3104                       NB = 1
3105                    ELSE
3106                       NB = 32
3107                    END IF
3108                 END IF
3109              END IF
3110           ELSE IF( C2.EQ.'TR' ) THEN
3111              IF( C3.EQ.'TRI' ) THEN
3112                 IF( SNAME ) THEN
3113                    NB = 64
3114                 ELSE
3115                    NB = 64
3116                 END IF
3117              END IF
3118           ELSE IF( C2.EQ.'LA' ) THEN
3119              IF( C3.EQ.'UUM' ) THEN
3120                 IF( SNAME ) THEN
3121                    NB = 64
3122                 ELSE
3123                    NB = 64
3124                 END IF
3125              END IF
3126           ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN
3127              IF( C3.EQ.'EBZ' ) THEN
3128                 NB = 1
3129              END IF
3130           END IF
3131           ILAENV = NB
3132           RETURN
3133     *
3134       200 CONTINUE
3135     *
3136     *     ISPEC = 2:  minimum block size
3137     *
3138           NBMIN = 2
3139           IF( C2.EQ.'GE' ) THEN
3140              IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
3141          \$       C3.EQ.'QLF' ) THEN
3142                 IF( SNAME ) THEN
3143                    NBMIN = 2
3144                 ELSE
3145                    NBMIN = 2
3146                 END IF
3147              ELSE IF( C3.EQ.'HRD' ) THEN
3148                 IF( SNAME ) THEN
3149                    NBMIN = 2
3150                 ELSE
3151                    NBMIN = 2
3152                 END IF
3153              ELSE IF( C3.EQ.'BRD' ) THEN
3154                 IF( SNAME ) THEN
3155                    NBMIN = 2
3156                 ELSE
3157                    NBMIN = 2
3158                 END IF
3159              ELSE IF( C3.EQ.'TRI' ) THEN
3160                 IF( SNAME ) THEN
3161                    NBMIN = 2
3162                 ELSE
3163                    NBMIN = 2
3164                 END IF
3165              END IF
3166           ELSE IF( C2.EQ.'SY' ) THEN
3167              IF( C3.EQ.'TRF' ) THEN
3168                 IF( SNAME ) THEN
3169                    NBMIN = 8
3170                 ELSE
3171                    NBMIN = 8
3172                 END IF
3173              ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN
3174                 NBMIN = 2
3175              END IF
3176           ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
3177              IF( C3.EQ.'TRD' ) THEN
3178                 NBMIN = 2
3179              END IF
3180           ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
3181              IF( C3( 1:1 ).EQ.'G' ) THEN
3182                 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
3183          \$          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
3184          \$          C4.EQ.'BR' ) THEN
3185                    NBMIN = 2
3186                 END IF
3187              ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
3188                 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
3189          \$          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
3190          \$          C4.EQ.'BR' ) THEN
3191                    NBMIN = 2
3192                 END IF
3193              END IF
3194           ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
3195              IF( C3( 1:1 ).EQ.'G' ) THEN
3196                 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
3197          \$          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
3198          \$          C4.EQ.'BR' ) THEN
3199                    NBMIN = 2
3200                 END IF
3201              ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
3202                 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
3203          \$          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
3204          \$          C4.EQ.'BR' ) THEN
3205                    NBMIN = 2
3206                 END IF
3207              END IF
3208           END IF
3209           ILAENV = NBMIN
3210           RETURN
3211     *
3212       300 CONTINUE
3213     *
3214     *     ISPEC = 3:  crossover point
3215     *
3216           NX = 0
3217           IF( C2.EQ.'GE' ) THEN
3218              IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
3219          \$       C3.EQ.'QLF' ) THEN
3220                 IF( SNAME ) THEN
3221                    NX = 128
3222                 ELSE
3223                    NX = 128
3224                 END IF
3225              ELSE IF( C3.EQ.'HRD' ) THEN
3226                 IF( SNAME ) THEN
3227                    NX = 128
3228                 ELSE
3229                    NX = 128
3230                 END IF
3231              ELSE IF( C3.EQ.'BRD' ) THEN
3232                 IF( SNAME ) THEN
3233                    NX = 128
3234                 ELSE
3235                    NX = 128
3236                 END IF
3237              END IF
3238           ELSE IF( C2.EQ.'SY' ) THEN
3239              IF( SNAME .AND. C3.EQ.'TRD' ) THEN
3240                 NX = 1
3241              END IF
3242           ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
3243              IF( C3.EQ.'TRD' ) THEN
3244                 NX = 1
3245              END IF
3246           ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
3247              IF( C3( 1:1 ).EQ.'G' ) THEN
3248                 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
3249          \$          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
3250          \$          C4.EQ.'BR' ) THEN
3251                    NX = 128
3252                 END IF
3253              END IF
3254           ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
3255              IF( C3( 1:1 ).EQ.'G' ) THEN
3256                 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
3257          \$          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
3258          \$          C4.EQ.'BR' ) THEN
3259                    NX = 128
3260                 END IF
3261              END IF
3262           END IF
3263           ILAENV = NX
3264           RETURN
3265     *
3266       400 CONTINUE
3267     *
3268     *     ISPEC = 4:  number of shifts (used by xHSEQR)
3269     *
3270           ILAENV = 6
3271           RETURN
3272     *
3273       500 CONTINUE
3274     *
3275     *     ISPEC = 5:  minimum column dimension (not used)
3276     *
3277           ILAENV = 2
3278           RETURN
3279     *
3280       600 CONTINUE
3281     *
3282     *     ISPEC = 6:  crossover point for SVD (used by xGELSS and xGESVD)
3283     *
3284           ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 )
3285           RETURN
3286     *
3287       700 CONTINUE
3288     *
3289     *     ISPEC = 7:  number of processors (not used)
3290     *
3291           ILAENV = 1
3292           RETURN
3293     *
3294       800 CONTINUE
3295     *
3296     *     ISPEC = 8:  crossover point for multishift (used by xHSEQR)
3297     *
3298           ILAENV = 50
3299           RETURN
3300     *
3301     *     End of ILAENV
3302     *
3303           END
3304
3305
3306
3307
3308
3309           LOGICAL          FUNCTION LSAME( CA, CB )
3310     *
3311     *  -- LAPACK auxiliary routine (version 2.0) --
3312     *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
3313     *     Courant Institute, Argonne National Lab, and Rice University
3314     *     September 30, 1994
3315     *
3316     *     .. Scalar Arguments ..
3317           CHARACTER          CA, CB
3318     *     ..
3319     *
3320     *  Purpose
3321     *  =======
3322     *
3323     *  LSAME returns .TRUE. if CA is the same letter as CB regardless of
3324     *  case.
3325     *
3326     *  Arguments
3327     *  =========
3328     *
3329     *  CA      (input) CHARACTER*1
3330     *  CB      (input) CHARACTER*1
3331     *          CA and CB specify the single characters to be compared.
3332     *
3333     * =====================================================================
3334     *
3335     *     .. Intrinsic Functions ..
3336           INTRINSIC          ICHAR
3337     *     ..
3338     *     .. Local Scalars ..
3339           INTEGER            INTA, INTB, ZCODE
3340     *     ..
3341     *     .. Executable Statements ..
3342     *
3343     *     Test if the characters are equal
3344     *
3345           LSAME = CA.EQ.CB
3346           IF( LSAME )
3347          \$   RETURN
3348     *
3349     *     Now test for equivalence if both characters are alphabetic.
3350     *
3351           ZCODE = MOVA2I( 'Z' )
3352     *
3353     *     Use 'Z' rather than 'A' so that ASCII can be detected on Prime
3354     *     machines, on which ICHAR returns a value with bit 8 set.
3355     *     MOVA2I('A') on Prime machines returns 193 which is the same as
3356     *     MOVA2I('A') on an EBCDIC machine.
3357     *
3358           INTA = MOVA2I( CA )
3359           INTB = MOVA2I( CB )
3360     *
3361           IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN
3362     *
3363     *        ASCII is assumed - ZCODE is the ASCII code of either lower or
3364     *        upper case 'Z'.
3365     *
3366              IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32
3367              IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32
3368     *
3369           ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN
3370     *
3371     *        EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
3372     *        upper case 'Z'.
3373     *
3374              IF( INTA.GE.129 .AND. INTA.LE.137 .OR.
3375          \$       INTA.GE.145 .AND. INTA.LE.153 .OR.
3376          \$       INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64
3377              IF( INTB.GE.129 .AND. INTB.LE.137 .OR.
3378          \$       INTB.GE.145 .AND. INTB.LE.153 .OR.
3379          \$       INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64
3380     *
3381           ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN
3382     *
3383     *        ASCII is assumed, on Prime machines - ZCODE is the ASCII code
3384     *        plus 128 of either lower or upper case 'Z'.
3385     *
3386              IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32
3387              IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32
3388           END IF
3389           LSAME = INTA.EQ.INTB
3390     *
3391     *     RETURN
3392     *
3393     *     End of LSAME
3394     *
3395           END
3396
3397
3398
3399           integer function isamax(n,sx,incx)
3400     c
3401     c     finds the index of element having max. absolute value.
3402     c     jack dongarra, linpack, 3/11/78.
3403     c     modified 3/93 to return if incx .le. 0.
3404     c     modified 12/3/93, array(1) declarations changed to array(*)
3405     c
3406           real sx(*),smax
3407           integer i,incx,ix,n
3408     c
3409           isamax = 0
3410           if( n.lt.1 .or. incx.le.0 ) return
3411           isamax = 1
3412           if(n.eq.1)return
3413           if(incx.eq.1)go to 20
3414     c
3415     c        code for increment not equal to 1
3416     c
3417           ix = 1
3418           smax = abs(sx(1))
3419           ix = ix + incx
3420           do 10 i = 2,n
3421              if(abs(sx(ix)).le.smax) go to 5
3422              isamax = i
3423              smax = abs(sx(ix))
3424         5    ix = ix + incx
3425        10 continue
3426           return
3427     c
3428     c        code for increment equal to 1
3429     c
3430        20 smax = abs(sx(1))
3431           do 30 i = 2,n
3432              if(abs(sx(i)).le.smax) go to 30
3433              isamax = i
3434              smax = abs(sx(i))
3435        30 continue
3436           return
3437           end
3438
3439
3440           subroutine dcopy(n,sx,incx,sy,incy)
3441     c
3442     c     copies a vector, x, to a vector, y.
3443     c     uses unrolled loops for increments equal to 1.
3444     c     jack dongarra, linpack, 3/11/78.
3445     c     modified 12/3/93, array(1) declarations changed to array(*)
3446     c
3447           real sx(*),sy(*)
3448           integer i,incx,incy,ix,iy,m,mp1,n
3449     c
3450           if(n.le.0)return
3451           if(incx.eq.1.and.incy.eq.1)go to 20
3452     c
3453     c        code for unequal increments or equal increments
3454     c          not equal to 1
3455     c
3456           ix = 1
3457           iy = 1
3458           if(incx.lt.0)ix = (-n+1)*incx + 1
3459           if(incy.lt.0)iy = (-n+1)*incy + 1
3460           do 10 i = 1,n
3461             sy(iy) = sx(ix)
3462             ix = ix + incx
3463             iy = iy + incy
3464        10 continue
3465           return
3466     c
3467     c        code for both increments equal to 1
3468     c
3469     c
3470     c        clean-up loop
3471     c
3472        20 m = mod(n,7)
3473           if( m .eq. 0 ) go to 40
3474           do 30 i = 1,m
3475             sy(i) = sx(i)
3476        30 continue
3477           if( n .lt. 7 ) return
3478        40 mp1 = m + 1
3479           do 50 i = mp1,n,7
3480             sy(i) = sx(i)
3481             sy(i + 1) = sx(i + 1)
3482             sy(i + 2) = sx(i + 2)
3483             sy(i + 3) = sx(i + 3)
3484             sy(i + 4) = sx(i + 4)
3485             sy(i + 5) = sx(i + 5)
3486             sy(i + 6) = sx(i + 6)
3487        50 continue
3488           return
3489           end
3490     ```