OB.DAAC Logo
NASA Logo
Ocean Color Science Software

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


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster