[Bio] / Numeric / Src / blas_lite.c Repository:
ViewVC logotype

Annotation of /Numeric/Src/blas_lite.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (view) (download) (as text)

1 : efrank 1.1 #include "Numeric/f2c.h"
2 :    
3 :     /* Subroutine */ int daxpy_(integer *n, doublereal *da, doublereal *dx,
4 :     integer *incx, doublereal *dy, integer *incy)
5 :     {
6 :     /* System generated locals */
7 :     integer i__1;
8 :     /* Local variables */
9 :     static integer i__, m, ix, iy, mp1;
10 :     /* constant times a vector plus a vector.
11 :     uses unrolled loops for increments equal to one.
12 :     jack dongarra, linpack, 3/11/78.
13 :     modified 12/3/93, array(1) declarations changed to array(*)
14 :     Parameter adjustments */
15 :     --dy;
16 :     --dx;
17 :     /* Function Body */
18 :     if (*n <= 0) {
19 :     return 0;
20 :     }
21 :     if (*da == 0.) {
22 :     return 0;
23 :     }
24 :     if (*incx == 1 && *incy == 1) {
25 :     goto L20;
26 :     }
27 :     /* code for unequal increments or equal increments
28 :     not equal to 1 */
29 :     ix = 1;
30 :     iy = 1;
31 :     if (*incx < 0) {
32 :     ix = (-(*n) + 1) * *incx + 1;
33 :     }
34 :     if (*incy < 0) {
35 :     iy = (-(*n) + 1) * *incy + 1;
36 :     }
37 :     i__1 = *n;
38 :     for (i__ = 1; i__ <= i__1; ++i__) {
39 :     dy[iy] += *da * dx[ix];
40 :     ix += *incx;
41 :     iy += *incy;
42 :     /* L10: */
43 :     }
44 :     return 0;
45 :     /* code for both increments equal to 1
46 :     clean-up loop */
47 :     L20:
48 :     m = *n % 4;
49 :     if (m == 0) {
50 :     goto L40;
51 :     }
52 :     i__1 = m;
53 :     for (i__ = 1; i__ <= i__1; ++i__) {
54 :     dy[i__] += *da * dx[i__];
55 :     /* L30: */
56 :     }
57 :     if (*n < 4) {
58 :     return 0;
59 :     }
60 :     L40:
61 :     mp1 = m + 1;
62 :     i__1 = *n;
63 :     for (i__ = mp1; i__ <= i__1; i__ += 4) {
64 :     dy[i__] += *da * dx[i__];
65 :     dy[i__ + 1] += *da * dx[i__ + 1];
66 :     dy[i__ + 2] += *da * dx[i__ + 2];
67 :     dy[i__ + 3] += *da * dx[i__ + 3];
68 :     /* L50: */
69 :     }
70 :     return 0;
71 :     } /* daxpy_ */
72 :    
73 :    
74 :     doublereal dcabs1_(doublecomplex *z__)
75 :     {
76 :     /* System generated locals */
77 :     doublereal ret_val;
78 :     static doublecomplex equiv_0[1];
79 :     /* Local variables */
80 :     #define t ((doublereal *)equiv_0)
81 :     #define zz (equiv_0)
82 :     zz->r = z__->r, zz->i = z__->i;
83 :     ret_val = abs(t[0]) + abs(t[1]);
84 :     return ret_val;
85 :     } /* dcabs1_ */
86 :     #undef zz
87 :     #undef t
88 :    
89 :    
90 :     /* Subroutine */ int dcopy_(integer *n, doublereal *dx, integer *incx,
91 :     doublereal *dy, integer *incy)
92 :     {
93 :     /* System generated locals */
94 :     integer i__1;
95 :     /* Local variables */
96 :     static integer i__, m, ix, iy, mp1;
97 :     /* copies a vector, x, to a vector, y.
98 :     uses unrolled loops for increments equal to one.
99 :     jack dongarra, linpack, 3/11/78.
100 :     modified 12/3/93, array(1) declarations changed to array(*)
101 :     Parameter adjustments */
102 :     --dy;
103 :     --dx;
104 :     /* Function Body */
105 :     if (*n <= 0) {
106 :     return 0;
107 :     }
108 :     if (*incx == 1 && *incy == 1) {
109 :     goto L20;
110 :     }
111 :     /* code for unequal increments or equal increments
112 :     not equal to 1 */
113 :     ix = 1;
114 :     iy = 1;
115 :     if (*incx < 0) {
116 :     ix = (-(*n) + 1) * *incx + 1;
117 :     }
118 :     if (*incy < 0) {
119 :     iy = (-(*n) + 1) * *incy + 1;
120 :     }
121 :     i__1 = *n;
122 :     for (i__ = 1; i__ <= i__1; ++i__) {
123 :     dy[iy] = dx[ix];
124 :     ix += *incx;
125 :     iy += *incy;
126 :     /* L10: */
127 :     }
128 :     return 0;
129 :     /* code for both increments equal to 1
130 :     clean-up loop */
131 :     L20:
132 :     m = *n % 7;
133 :     if (m == 0) {
134 :     goto L40;
135 :     }
136 :     i__1 = m;
137 :     for (i__ = 1; i__ <= i__1; ++i__) {
138 :     dy[i__] = dx[i__];
139 :     /* L30: */
140 :     }
141 :     if (*n < 7) {
142 :     return 0;
143 :     }
144 :     L40:
145 :     mp1 = m + 1;
146 :     i__1 = *n;
147 :     for (i__ = mp1; i__ <= i__1; i__ += 7) {
148 :     dy[i__] = dx[i__];
149 :     dy[i__ + 1] = dx[i__ + 1];
150 :     dy[i__ + 2] = dx[i__ + 2];
151 :     dy[i__ + 3] = dx[i__ + 3];
152 :     dy[i__ + 4] = dx[i__ + 4];
153 :     dy[i__ + 5] = dx[i__ + 5];
154 :     dy[i__ + 6] = dx[i__ + 6];
155 :     /* L50: */
156 :     }
157 :     return 0;
158 :     } /* dcopy_ */
159 :    
160 :    
161 :     doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy,
162 :     integer *incy)
163 :     {
164 :     /* System generated locals */
165 :     integer i__1;
166 :     doublereal ret_val;
167 :     /* Local variables */
168 :     static integer i__, m;
169 :     static doublereal dtemp;
170 :     static integer ix, iy, mp1;
171 :     /* forms the dot product of two vectors.
172 :     uses unrolled loops for increments equal to one.
173 :     jack dongarra, linpack, 3/11/78.
174 :     modified 12/3/93, array(1) declarations changed to array(*)
175 :     Parameter adjustments */
176 :     --dy;
177 :     --dx;
178 :     /* Function Body */
179 :     ret_val = 0.;
180 :     dtemp = 0.;
181 :     if (*n <= 0) {
182 :     return ret_val;
183 :     }
184 :     if (*incx == 1 && *incy == 1) {
185 :     goto L20;
186 :     }
187 :     /* code for unequal increments or equal increments
188 :     not equal to 1 */
189 :     ix = 1;
190 :     iy = 1;
191 :     if (*incx < 0) {
192 :     ix = (-(*n) + 1) * *incx + 1;
193 :     }
194 :     if (*incy < 0) {
195 :     iy = (-(*n) + 1) * *incy + 1;
196 :     }
197 :     i__1 = *n;
198 :     for (i__ = 1; i__ <= i__1; ++i__) {
199 :     dtemp += dx[ix] * dy[iy];
200 :     ix += *incx;
201 :     iy += *incy;
202 :     /* L10: */
203 :     }
204 :     ret_val = dtemp;
205 :     return ret_val;
206 :     /* code for both increments equal to 1
207 :     clean-up loop */
208 :     L20:
209 :     m = *n % 5;
210 :     if (m == 0) {
211 :     goto L40;
212 :     }
213 :     i__1 = m;
214 :     for (i__ = 1; i__ <= i__1; ++i__) {
215 :     dtemp += dx[i__] * dy[i__];
216 :     /* L30: */
217 :     }
218 :     if (*n < 5) {
219 :     goto L60;
220 :     }
221 :     L40:
222 :     mp1 = m + 1;
223 :     i__1 = *n;
224 :     for (i__ = mp1; i__ <= i__1; i__ += 5) {
225 :     dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[
226 :     i__ + 2] * dy[i__ + 2] + dx[i__ + 3] * dy[i__ + 3] + dx[i__ +
227 :     4] * dy[i__ + 4];
228 :     /* L50: */
229 :     }
230 :     L60:
231 :     ret_val = dtemp;
232 :     return ret_val;
233 :     } /* ddot_ */
234 :    
235 :    
236 :     /* Subroutine */ int dgemm_(char *transa, char *transb, integer *m, integer *
237 :     n, integer *k, doublereal *alpha, doublereal *a, integer *lda,
238 :     doublereal *b, integer *ldb, doublereal *beta, doublereal *c__,
239 :     integer *ldc)
240 :     {
241 :     /* System generated locals */
242 :     integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
243 :     i__3;
244 :     /* Local variables */
245 :     static integer info;
246 :     static logical nota, notb;
247 :     static doublereal temp;
248 :     static integer i__, j, l, ncola;
249 :     extern logical lsame_(char *, char *);
250 :     static integer nrowa, nrowb;
251 :     extern /* Subroutine */ int xerbla_(char *, integer *);
252 :     #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
253 :     #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
254 :     #define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1]
255 :     /* Purpose
256 :     =======
257 :     DGEMM performs one of the matrix-matrix operations
258 :     C := alpha*op( A )*op( B ) + beta*C,
259 :     where op( X ) is one of
260 :     op( X ) = X or op( X ) = X',
261 :     alpha and beta are scalars, and A, B and C are matrices, with op( A )
262 :     an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
263 :     Parameters
264 :     ==========
265 :     TRANSA - CHARACTER*1.
266 :     On entry, TRANSA specifies the form of op( A ) to be used in
267 :     the matrix multiplication as follows:
268 :     TRANSA = 'N' or 'n', op( A ) = A.
269 :     TRANSA = 'T' or 't', op( A ) = A'.
270 :     TRANSA = 'C' or 'c', op( A ) = A'.
271 :     Unchanged on exit.
272 :     TRANSB - CHARACTER*1.
273 :     On entry, TRANSB specifies the form of op( B ) to be used in
274 :     the matrix multiplication as follows:
275 :     TRANSB = 'N' or 'n', op( B ) = B.
276 :     TRANSB = 'T' or 't', op( B ) = B'.
277 :     TRANSB = 'C' or 'c', op( B ) = B'.
278 :     Unchanged on exit.
279 :     M - INTEGER.
280 :     On entry, M specifies the number of rows of the matrix
281 :     op( A ) and of the matrix C. M must be at least zero.
282 :     Unchanged on exit.
283 :     N - INTEGER.
284 :     On entry, N specifies the number of columns of the matrix
285 :     op( B ) and the number of columns of the matrix C. N must be
286 :     at least zero.
287 :     Unchanged on exit.
288 :     K - INTEGER.
289 :     On entry, K specifies the number of columns of the matrix
290 :     op( A ) and the number of rows of the matrix op( B ). K must
291 :     be at least zero.
292 :     Unchanged on exit.
293 :     ALPHA - DOUBLE PRECISION.
294 :     On entry, ALPHA specifies the scalar alpha.
295 :     Unchanged on exit.
296 :     A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
297 :     k when TRANSA = 'N' or 'n', and is m otherwise.
298 :     Before entry with TRANSA = 'N' or 'n', the leading m by k
299 :     part of the array A must contain the matrix A, otherwise
300 :     the leading k by m part of the array A must contain the
301 :     matrix A.
302 :     Unchanged on exit.
303 :     LDA - INTEGER.
304 :     On entry, LDA specifies the first dimension of A as declared
305 :     in the calling (sub) program. When TRANSA = 'N' or 'n' then
306 :     LDA must be at least max( 1, m ), otherwise LDA must be at
307 :     least max( 1, k ).
308 :     Unchanged on exit.
309 :     B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
310 :     n when TRANSB = 'N' or 'n', and is k otherwise.
311 :     Before entry with TRANSB = 'N' or 'n', the leading k by n
312 :     part of the array B must contain the matrix B, otherwise
313 :     the leading n by k part of the array B must contain the
314 :     matrix B.
315 :     Unchanged on exit.
316 :     LDB - INTEGER.
317 :     On entry, LDB specifies the first dimension of B as declared
318 :     in the calling (sub) program. When TRANSB = 'N' or 'n' then
319 :     LDB must be at least max( 1, k ), otherwise LDB must be at
320 :     least max( 1, n ).
321 :     Unchanged on exit.
322 :     BETA - DOUBLE PRECISION.
323 :     On entry, BETA specifies the scalar beta. When BETA is
324 :     supplied as zero then C need not be set on input.
325 :     Unchanged on exit.
326 :     C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
327 :     Before entry, the leading m by n part of the array C must
328 :     contain the matrix C, except when beta is zero, in which
329 :     case C need not be set on entry.
330 :     On exit, the array C is overwritten by the m by n matrix
331 :     ( alpha*op( A )*op( B ) + beta*C ).
332 :     LDC - INTEGER.
333 :     On entry, LDC specifies the first dimension of C as declared
334 :     in the calling (sub) program. LDC must be at least
335 :     max( 1, m ).
336 :     Unchanged on exit.
337 :     Level 3 Blas routine.
338 :     -- Written on 8-February-1989.
339 :     Jack Dongarra, Argonne National Laboratory.
340 :     Iain Duff, AERE Harwell.
341 :     Jeremy Du Croz, Numerical Algorithms Group Ltd.
342 :     Sven Hammarling, Numerical Algorithms Group Ltd.
343 :     Set NOTA and NOTB as true if A and B respectively are not
344 :     transposed and set NROWA, NCOLA and NROWB as the number of rows
345 :     and columns of A and the number of rows of B respectively.
346 :     Parameter adjustments */
347 :     a_dim1 = *lda;
348 :     a_offset = 1 + a_dim1 * 1;
349 :     a -= a_offset;
350 :     b_dim1 = *ldb;
351 :     b_offset = 1 + b_dim1 * 1;
352 :     b -= b_offset;
353 :     c_dim1 = *ldc;
354 :     c_offset = 1 + c_dim1 * 1;
355 :     c__ -= c_offset;
356 :     /* Function Body */
357 :     nota = lsame_(transa, "N");
358 :     notb = lsame_(transb, "N");
359 :     if (nota) {
360 :     nrowa = *m;
361 :     ncola = *k;
362 :     } else {
363 :     nrowa = *k;
364 :     ncola = *m;
365 :     }
366 :     if (notb) {
367 :     nrowb = *k;
368 :     } else {
369 :     nrowb = *n;
370 :     }
371 :     /* Test the input parameters. */
372 :     info = 0;
373 :     if (! nota && ! lsame_(transa, "C") && ! lsame_(
374 :     transa, "T")) {
375 :     info = 1;
376 :     } else if (! notb && ! lsame_(transb, "C") && !
377 :     lsame_(transb, "T")) {
378 :     info = 2;
379 :     } else if (*m < 0) {
380 :     info = 3;
381 :     } else if (*n < 0) {
382 :     info = 4;
383 :     } else if (*k < 0) {
384 :     info = 5;
385 :     } else if (*lda < max(1,nrowa)) {
386 :     info = 8;
387 :     } else if (*ldb < max(1,nrowb)) {
388 :     info = 10;
389 :     } else if (*ldc < max(1,*m)) {
390 :     info = 13;
391 :     }
392 :     if (info != 0) {
393 :     xerbla_("DGEMM ", &info);
394 :     return 0;
395 :     }
396 :     /* Quick return if possible. */
397 :     if (*m == 0 || *n == 0 || (*alpha == 0. || *k == 0) && *beta == 1.) {
398 :     return 0;
399 :     }
400 :     /* And if alpha.eq.zero. */
401 :     if (*alpha == 0.) {
402 :     if (*beta == 0.) {
403 :     i__1 = *n;
404 :     for (j = 1; j <= i__1; ++j) {
405 :     i__2 = *m;
406 :     for (i__ = 1; i__ <= i__2; ++i__) {
407 :     c___ref(i__, j) = 0.;
408 :     /* L10: */
409 :     }
410 :     /* L20: */
411 :     }
412 :     } else {
413 :     i__1 = *n;
414 :     for (j = 1; j <= i__1; ++j) {
415 :     i__2 = *m;
416 :     for (i__ = 1; i__ <= i__2; ++i__) {
417 :     c___ref(i__, j) = *beta * c___ref(i__, j);
418 :     /* L30: */
419 :     }
420 :     /* L40: */
421 :     }
422 :     }
423 :     return 0;
424 :     }
425 :     /* Start the operations. */
426 :     if (notb) {
427 :     if (nota) {
428 :     /* Form C := alpha*A*B + beta*C. */
429 :     i__1 = *n;
430 :     for (j = 1; j <= i__1; ++j) {
431 :     if (*beta == 0.) {
432 :     i__2 = *m;
433 :     for (i__ = 1; i__ <= i__2; ++i__) {
434 :     c___ref(i__, j) = 0.;
435 :     /* L50: */
436 :     }
437 :     } else if (*beta != 1.) {
438 :     i__2 = *m;
439 :     for (i__ = 1; i__ <= i__2; ++i__) {
440 :     c___ref(i__, j) = *beta * c___ref(i__, j);
441 :     /* L60: */
442 :     }
443 :     }
444 :     i__2 = *k;
445 :     for (l = 1; l <= i__2; ++l) {
446 :     if (b_ref(l, j) != 0.) {
447 :     temp = *alpha * b_ref(l, j);
448 :     i__3 = *m;
449 :     for (i__ = 1; i__ <= i__3; ++i__) {
450 :     c___ref(i__, j) = c___ref(i__, j) + temp * a_ref(
451 :     i__, l);
452 :     /* L70: */
453 :     }
454 :     }
455 :     /* L80: */
456 :     }
457 :     /* L90: */
458 :     }
459 :     } else {
460 :     /* Form C := alpha*A'*B + beta*C */
461 :     i__1 = *n;
462 :     for (j = 1; j <= i__1; ++j) {
463 :     i__2 = *m;
464 :     for (i__ = 1; i__ <= i__2; ++i__) {
465 :     temp = 0.;
466 :     i__3 = *k;
467 :     for (l = 1; l <= i__3; ++l) {
468 :     temp += a_ref(l, i__) * b_ref(l, j);
469 :     /* L100: */
470 :     }
471 :     if (*beta == 0.) {
472 :     c___ref(i__, j) = *alpha * temp;
473 :     } else {
474 :     c___ref(i__, j) = *alpha * temp + *beta * c___ref(i__,
475 :     j);
476 :     }
477 :     /* L110: */
478 :     }
479 :     /* L120: */
480 :     }
481 :     }
482 :     } else {
483 :     if (nota) {
484 :     /* Form C := alpha*A*B' + beta*C */
485 :     i__1 = *n;
486 :     for (j = 1; j <= i__1; ++j) {
487 :     if (*beta == 0.) {
488 :     i__2 = *m;
489 :     for (i__ = 1; i__ <= i__2; ++i__) {
490 :     c___ref(i__, j) = 0.;
491 :     /* L130: */
492 :     }
493 :     } else if (*beta != 1.) {
494 :     i__2 = *m;
495 :     for (i__ = 1; i__ <= i__2; ++i__) {
496 :     c___ref(i__, j) = *beta * c___ref(i__, j);
497 :     /* L140: */
498 :     }
499 :     }
500 :     i__2 = *k;
501 :     for (l = 1; l <= i__2; ++l) {
502 :     if (b_ref(j, l) != 0.) {
503 :     temp = *alpha * b_ref(j, l);
504 :     i__3 = *m;
505 :     for (i__ = 1; i__ <= i__3; ++i__) {
506 :     c___ref(i__, j) = c___ref(i__, j) + temp * a_ref(
507 :     i__, l);
508 :     /* L150: */
509 :     }
510 :     }
511 :     /* L160: */
512 :     }
513 :     /* L170: */
514 :     }
515 :     } else {
516 :     /* Form C := alpha*A'*B' + beta*C */
517 :     i__1 = *n;
518 :     for (j = 1; j <= i__1; ++j) {
519 :     i__2 = *m;
520 :     for (i__ = 1; i__ <= i__2; ++i__) {
521 :     temp = 0.;
522 :     i__3 = *k;
523 :     for (l = 1; l <= i__3; ++l) {
524 :     temp += a_ref(l, i__) * b_ref(j, l);
525 :     /* L180: */
526 :     }
527 :     if (*beta == 0.) {
528 :     c___ref(i__, j) = *alpha * temp;
529 :     } else {
530 :     c___ref(i__, j) = *alpha * temp + *beta * c___ref(i__,
531 :     j);
532 :     }
533 :     /* L190: */
534 :     }
535 :     /* L200: */
536 :     }
537 :     }
538 :     }
539 :     return 0;
540 :     /* End of DGEMM . */
541 :     } /* dgemm_ */
542 :     #undef c___ref
543 :     #undef b_ref
544 :     #undef a_ref
545 :    
546 :    
547 :     /* Subroutine */ int dgemv_(char *trans, integer *m, integer *n, doublereal *
548 :     alpha, doublereal *a, integer *lda, doublereal *x, integer *incx,
549 :     doublereal *beta, doublereal *y, integer *incy)
550 :     {
551 :     /* System generated locals */
552 :     integer a_dim1, a_offset, i__1, i__2;
553 :     /* Local variables */
554 :     static integer info;
555 :     static doublereal temp;
556 :     static integer lenx, leny, i__, j;
557 :     extern logical lsame_(char *, char *);
558 :     static integer ix, iy, jx, jy, kx, ky;
559 :     extern /* Subroutine */ int xerbla_(char *, integer *);
560 :     #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
561 :     /* Purpose
562 :     =======
563 :     DGEMV performs one of the matrix-vector operations
564 :     y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y,
565 :     where alpha and beta are scalars, x and y are vectors and A is an
566 :     m by n matrix.
567 :     Parameters
568 :     ==========
569 :     TRANS - CHARACTER*1.
570 :     On entry, TRANS specifies the operation to be performed as
571 :     follows:
572 :     TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
573 :     TRANS = 'T' or 't' y := alpha*A'*x + beta*y.
574 :     TRANS = 'C' or 'c' y := alpha*A'*x + beta*y.
575 :     Unchanged on exit.
576 :     M - INTEGER.
577 :     On entry, M specifies the number of rows of the matrix A.
578 :     M must be at least zero.
579 :     Unchanged on exit.
580 :     N - INTEGER.
581 :     On entry, N specifies the number of columns of the matrix A.
582 :     N must be at least zero.
583 :     Unchanged on exit.
584 :     ALPHA - DOUBLE PRECISION.
585 :     On entry, ALPHA specifies the scalar alpha.
586 :     Unchanged on exit.
587 :     A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
588 :     Before entry, the leading m by n part of the array A must
589 :     contain the matrix of coefficients.
590 :     Unchanged on exit.
591 :     LDA - INTEGER.
592 :     On entry, LDA specifies the first dimension of A as declared
593 :     in the calling (sub) program. LDA must be at least
594 :     max( 1, m ).
595 :     Unchanged on exit.
596 :     X - DOUBLE PRECISION array of DIMENSION at least
597 :     ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
598 :     and at least
599 :     ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
600 :     Before entry, the incremented array X must contain the
601 :     vector x.
602 :     Unchanged on exit.
603 :     INCX - INTEGER.
604 :     On entry, INCX specifies the increment for the elements of
605 :     X. INCX must not be zero.
606 :     Unchanged on exit.
607 :     BETA - DOUBLE PRECISION.
608 :     On entry, BETA specifies the scalar beta. When BETA is
609 :     supplied as zero then Y need not be set on input.
610 :     Unchanged on exit.
611 :     Y - DOUBLE PRECISION array of DIMENSION at least
612 :     ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
613 :     and at least
614 :     ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
615 :     Before entry with BETA non-zero, the incremented array Y
616 :     must contain the vector y. On exit, Y is overwritten by the
617 :     updated vector y.
618 :     INCY - INTEGER.
619 :     On entry, INCY specifies the increment for the elements of
620 :     Y. INCY must not be zero.
621 :     Unchanged on exit.
622 :     Level 2 Blas routine.
623 :     -- Written on 22-October-1986.
624 :     Jack Dongarra, Argonne National Lab.
625 :     Jeremy Du Croz, Nag Central Office.
626 :     Sven Hammarling, Nag Central Office.
627 :     Richard Hanson, Sandia National Labs.
628 :     Test the input parameters.
629 :     Parameter adjustments */
630 :     a_dim1 = *lda;
631 :     a_offset = 1 + a_dim1 * 1;
632 :     a -= a_offset;
633 :     --x;
634 :     --y;
635 :     /* Function Body */
636 :     info = 0;
637 :     if (! lsame_(trans, "N") && ! lsame_(trans, "T") && ! lsame_(trans, "C")
638 :     ) {
639 :     info = 1;
640 :     } else if (*m < 0) {
641 :     info = 2;
642 :     } else if (*n < 0) {
643 :     info = 3;
644 :     } else if (*lda < max(1,*m)) {
645 :     info = 6;
646 :     } else if (*incx == 0) {
647 :     info = 8;
648 :     } else if (*incy == 0) {
649 :     info = 11;
650 :     }
651 :     if (info != 0) {
652 :     xerbla_("DGEMV ", &info);
653 :     return 0;
654 :     }
655 :     /* Quick return if possible. */
656 :     if (*m == 0 || *n == 0 || *alpha == 0. && *beta == 1.) {
657 :     return 0;
658 :     }
659 :     /* Set LENX and LENY, the lengths of the vectors x and y, and set
660 :     up the start points in X and Y. */
661 :     if (lsame_(trans, "N")) {
662 :     lenx = *n;
663 :     leny = *m;
664 :     } else {
665 :     lenx = *m;
666 :     leny = *n;
667 :     }
668 :     if (*incx > 0) {
669 :     kx = 1;
670 :     } else {
671 :     kx = 1 - (lenx - 1) * *incx;
672 :     }
673 :     if (*incy > 0) {
674 :     ky = 1;
675 :     } else {
676 :     ky = 1 - (leny - 1) * *incy;
677 :     }
678 :     /* Start the operations. In this version the elements of A are
679 :     accessed sequentially with one pass through A.
680 :     First form y := beta*y. */
681 :     if (*beta != 1.) {
682 :     if (*incy == 1) {
683 :     if (*beta == 0.) {
684 :     i__1 = leny;
685 :     for (i__ = 1; i__ <= i__1; ++i__) {
686 :     y[i__] = 0.;
687 :     /* L10: */
688 :     }
689 :     } else {
690 :     i__1 = leny;
691 :     for (i__ = 1; i__ <= i__1; ++i__) {
692 :     y[i__] = *beta * y[i__];
693 :     /* L20: */
694 :     }
695 :     }
696 :     } else {
697 :     iy = ky;
698 :     if (*beta == 0.) {
699 :     i__1 = leny;
700 :     for (i__ = 1; i__ <= i__1; ++i__) {
701 :     y[iy] = 0.;
702 :     iy += *incy;
703 :     /* L30: */
704 :     }
705 :     } else {
706 :     i__1 = leny;
707 :     for (i__ = 1; i__ <= i__1; ++i__) {
708 :     y[iy] = *beta * y[iy];
709 :     iy += *incy;
710 :     /* L40: */
711 :     }
712 :     }
713 :     }
714 :     }
715 :     if (*alpha == 0.) {
716 :     return 0;
717 :     }
718 :     if (lsame_(trans, "N")) {
719 :     /* Form y := alpha*A*x + y. */
720 :     jx = kx;
721 :     if (*incy == 1) {
722 :     i__1 = *n;
723 :     for (j = 1; j <= i__1; ++j) {
724 :     if (x[jx] != 0.) {
725 :     temp = *alpha * x[jx];
726 :     i__2 = *m;
727 :     for (i__ = 1; i__ <= i__2; ++i__) {
728 :     y[i__] += temp * a_ref(i__, j);
729 :     /* L50: */
730 :     }
731 :     }
732 :     jx += *incx;
733 :     /* L60: */
734 :     }
735 :     } else {
736 :     i__1 = *n;
737 :     for (j = 1; j <= i__1; ++j) {
738 :     if (x[jx] != 0.) {
739 :     temp = *alpha * x[jx];
740 :     iy = ky;
741 :     i__2 = *m;
742 :     for (i__ = 1; i__ <= i__2; ++i__) {
743 :     y[iy] += temp * a_ref(i__, j);
744 :     iy += *incy;
745 :     /* L70: */
746 :     }
747 :     }
748 :     jx += *incx;
749 :     /* L80: */
750 :     }
751 :     }
752 :     } else {
753 :     /* Form y := alpha*A'*x + y. */
754 :     jy = ky;
755 :     if (*incx == 1) {
756 :     i__1 = *n;
757 :     for (j = 1; j <= i__1; ++j) {
758 :     temp = 0.;
759 :     i__2 = *m;
760 :     for (i__ = 1; i__ <= i__2; ++i__) {
761 :     temp += a_ref(i__, j) * x[i__];
762 :     /* L90: */
763 :     }
764 :     y[jy] += *alpha * temp;
765 :     jy += *incy;
766 :     /* L100: */
767 :     }
768 :     } else {
769 :     i__1 = *n;
770 :     for (j = 1; j <= i__1; ++j) {
771 :     temp = 0.;
772 :     ix = kx;
773 :     i__2 = *m;
774 :     for (i__ = 1; i__ <= i__2; ++i__) {
775 :     temp += a_ref(i__, j) * x[ix];
776 :     ix += *incx;
777 :     /* L110: */
778 :     }
779 :     y[jy] += *alpha * temp;
780 :     jy += *incy;
781 :     /* L120: */
782 :     }
783 :     }
784 :     }
785 :     return 0;
786 :     /* End of DGEMV . */
787 :     } /* dgemv_ */
788 :     #undef a_ref
789 :    
790 :    
791 :     /* Subroutine */ int dger_(integer *m, integer *n, doublereal *alpha,
792 :     doublereal *x, integer *incx, doublereal *y, integer *incy,
793 :     doublereal *a, integer *lda)
794 :     {
795 :     /* System generated locals */
796 :     integer a_dim1, a_offset, i__1, i__2;
797 :     /* Local variables */
798 :     static integer info;
799 :     static doublereal temp;
800 :     static integer i__, j, ix, jy, kx;
801 :     extern /* Subroutine */ int xerbla_(char *, integer *);
802 :     #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
803 :     /* Purpose
804 :     =======
805 :     DGER performs the rank 1 operation
806 :     A := alpha*x*y' + A,
807 :     where alpha is a scalar, x is an m element vector, y is an n element
808 :     vector and A is an m by n matrix.
809 :     Parameters
810 :     ==========
811 :     M - INTEGER.
812 :     On entry, M specifies the number of rows of the matrix A.
813 :     M must be at least zero.
814 :     Unchanged on exit.
815 :     N - INTEGER.
816 :     On entry, N specifies the number of columns of the matrix A.
817 :     N must be at least zero.
818 :     Unchanged on exit.
819 :     ALPHA - DOUBLE PRECISION.
820 :     On entry, ALPHA specifies the scalar alpha.
821 :     Unchanged on exit.
822 :     X - DOUBLE PRECISION array of dimension at least
823 :     ( 1 + ( m - 1 )*abs( INCX ) ).
824 :     Before entry, the incremented array X must contain the m
825 :     element vector x.
826 :     Unchanged on exit.
827 :     INCX - INTEGER.
828 :     On entry, INCX specifies the increment for the elements of
829 :     X. INCX must not be zero.
830 :     Unchanged on exit.
831 :     Y - DOUBLE PRECISION array of dimension at least
832 :     ( 1 + ( n - 1 )*abs( INCY ) ).
833 :     Before entry, the incremented array Y must contain the n
834 :     element vector y.
835 :     Unchanged on exit.
836 :     INCY - INTEGER.
837 :     On entry, INCY specifies the increment for the elements of
838 :     Y. INCY must not be zero.
839 :     Unchanged on exit.
840 :     A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
841 :     Before entry, the leading m by n part of the array A must
842 :     contain the matrix of coefficients. On exit, A is
843 :     overwritten by the updated matrix.
844 :     LDA - INTEGER.
845 :     On entry, LDA specifies the first dimension of A as declared
846 :     in the calling (sub) program. LDA must be at least
847 :     max( 1, m ).
848 :     Unchanged on exit.
849 :     Level 2 Blas routine.
850 :     -- Written on 22-October-1986.
851 :     Jack Dongarra, Argonne National Lab.
852 :     Jeremy Du Croz, Nag Central Office.
853 :     Sven Hammarling, Nag Central Office.
854 :     Richard Hanson, Sandia National Labs.
855 :     Test the input parameters.
856 :     Parameter adjustments */
857 :     --x;
858 :     --y;
859 :     a_dim1 = *lda;
860 :     a_offset = 1 + a_dim1 * 1;
861 :     a -= a_offset;
862 :     /* Function Body */
863 :     info = 0;
864 :     if (*m < 0) {
865 :     info = 1;
866 :     } else if (*n < 0) {
867 :     info = 2;
868 :     } else if (*incx == 0) {
869 :     info = 5;
870 :     } else if (*incy == 0) {
871 :     info = 7;
872 :     } else if (*lda < max(1,*m)) {
873 :     info = 9;
874 :     }
875 :     if (info != 0) {
876 :     xerbla_("DGER ", &info);
877 :     return 0;
878 :     }
879 :     /* Quick return if possible. */
880 :     if (*m == 0 || *n == 0 || *alpha == 0.) {
881 :     return 0;
882 :     }
883 :     /* Start the operations. In this version the elements of A are
884 :     accessed sequentially with one pass through A. */
885 :     if (*incy > 0) {
886 :     jy = 1;
887 :     } else {
888 :     jy = 1 - (*n - 1) * *incy;
889 :     }
890 :     if (*incx == 1) {
891 :     i__1 = *n;
892 :     for (j = 1; j <= i__1; ++j) {
893 :     if (y[jy] != 0.) {
894 :     temp = *alpha * y[jy];
895 :     i__2 = *m;
896 :     for (i__ = 1; i__ <= i__2; ++i__) {
897 :     a_ref(i__, j) = a_ref(i__, j) + x[i__] * temp;
898 :     /* L10: */
899 :     }
900 :     }
901 :     jy += *incy;
902 :     /* L20: */
903 :     }
904 :     } else {
905 :     if (*incx > 0) {
906 :     kx = 1;
907 :     } else {
908 :     kx = 1 - (*m - 1) * *incx;
909 :     }
910 :     i__1 = *n;
911 :     for (j = 1; j <= i__1; ++j) {
912 :     if (y[jy] != 0.) {
913 :     temp = *alpha * y[jy];
914 :     ix = kx;
915 :     i__2 = *m;
916 :     for (i__ = 1; i__ <= i__2; ++i__) {
917 :     a_ref(i__, j) = a_ref(i__, j) + x[ix] * temp;
918 :     ix += *incx;
919 :     /* L30: */
920 :     }
921 :     }
922 :     jy += *incy;
923 :     /* L40: */
924 :     }
925 :     }
926 :     return 0;
927 :     /* End of DGER . */
928 :     } /* dger_ */
929 :     #undef a_ref
930 :    
931 :    
932 :     doublereal dnrm2_(integer *n, doublereal *x, integer *incx)
933 :     {
934 :     /* The following loop is equivalent to this call to the LAPACK
935 :     auxiliary routine:
936 :     CALL DLASSQ( N, X, INCX, SCALE, SSQ ) */
937 :     /* System generated locals */
938 :     integer i__1, i__2;
939 :     doublereal ret_val, d__1;
940 :     /* Builtin functions */
941 :     double sqrt(doublereal);
942 :     /* Local variables */
943 :     static doublereal norm, scale, absxi;
944 :     static integer ix;
945 :     static doublereal ssq;
946 :     /* DNRM2 returns the euclidean norm of a vector via the function
947 :     name, so that
948 :     DNRM2 := sqrt( x'*x )
949 :     -- This version written on 25-October-1982.
950 :     Modified on 14-October-1993 to inline the call to DLASSQ.
951 :     Sven Hammarling, Nag Ltd.
952 :     Parameter adjustments */
953 :     --x;
954 :     /* Function Body */
955 :     if (*n < 1 || *incx < 1) {
956 :     norm = 0.;
957 :     } else if (*n == 1) {
958 :     norm = abs(x[1]);
959 :     } else {
960 :     scale = 0.;
961 :     ssq = 1.;
962 :    
963 :    
964 :     i__1 = (*n - 1) * *incx + 1;
965 :     i__2 = *incx;
966 :     for (ix = 1; i__2 < 0 ? ix >= i__1 : ix <= i__1; ix += i__2) {
967 :     if (x[ix] != 0.) {
968 :     absxi = (d__1 = x[ix], abs(d__1));
969 :     if (scale < absxi) {
970 :     /* Computing 2nd power */
971 :     d__1 = scale / absxi;
972 :     ssq = ssq * (d__1 * d__1) + 1.;
973 :     scale = absxi;
974 :     } else {
975 :     /* Computing 2nd power */
976 :     d__1 = absxi / scale;
977 :     ssq += d__1 * d__1;
978 :     }
979 :     }
980 :     /* L10: */
981 :     }
982 :     norm = scale * sqrt(ssq);
983 :     }
984 :    
985 :     ret_val = norm;
986 :     return ret_val;
987 :    
988 :     /* End of DNRM2. */
989 :    
990 :     } /* dnrm2_ */
991 :    
992 :    
993 :     /* Subroutine */ int drot_(integer *n, doublereal *dx, integer *incx,
994 :     doublereal *dy, integer *incy, doublereal *c__, doublereal *s)
995 :     {
996 :     /* System generated locals */
997 :     integer i__1;
998 :     /* Local variables */
999 :     static integer i__;
1000 :     static doublereal dtemp;
1001 :     static integer ix, iy;
1002 :     /* applies a plane rotation.
1003 :     jack dongarra, linpack, 3/11/78.
1004 :     modified 12/3/93, array(1) declarations changed to array(*)
1005 :     Parameter adjustments */
1006 :     --dy;
1007 :     --dx;
1008 :     /* Function Body */
1009 :     if (*n <= 0) {
1010 :     return 0;
1011 :     }
1012 :     if (*incx == 1 && *incy == 1) {
1013 :     goto L20;
1014 :     }
1015 :     /* code for unequal increments or equal increments not equal
1016 :     to 1 */
1017 :     ix = 1;
1018 :     iy = 1;
1019 :     if (*incx < 0) {
1020 :     ix = (-(*n) + 1) * *incx + 1;
1021 :     }
1022 :     if (*incy < 0) {
1023 :     iy = (-(*n) + 1) * *incy + 1;
1024 :     }
1025 :     i__1 = *n;
1026 :     for (i__ = 1; i__ <= i__1; ++i__) {
1027 :     dtemp = *c__ * dx[ix] + *s * dy[iy];
1028 :     dy[iy] = *c__ * dy[iy] - *s * dx[ix];
1029 :     dx[ix] = dtemp;
1030 :     ix += *incx;
1031 :     iy += *incy;
1032 :     /* L10: */
1033 :     }
1034 :     return 0;
1035 :     /* code for both increments equal to 1 */
1036 :     L20:
1037 :     i__1 = *n;
1038 :     for (i__ = 1; i__ <= i__1; ++i__) {
1039 :     dtemp = *c__ * dx[i__] + *s * dy[i__];
1040 :     dy[i__] = *c__ * dy[i__] - *s * dx[i__];
1041 :     dx[i__] = dtemp;
1042 :     /* L30: */
1043 :     }
1044 :     return 0;
1045 :     } /* drot_ */
1046 :    
1047 :    
1048 :     /* Subroutine */ int dscal_(integer *n, doublereal *da, doublereal *dx,
1049 :     integer *incx)
1050 :     {
1051 :     /* System generated locals */
1052 :     integer i__1, i__2;
1053 :     /* Local variables */
1054 :     static integer i__, m, nincx, mp1;
1055 :     /* scales a vector by a constant.
1056 :     uses unrolled loops for increment equal to one.
1057 :     jack dongarra, linpack, 3/11/78.
1058 :     modified 3/93 to return if incx .le. 0.
1059 :     modified 12/3/93, array(1) declarations changed to array(*)
1060 :     Parameter adjustments */
1061 :     --dx;
1062 :     /* Function Body */
1063 :     if (*n <= 0 || *incx <= 0) {
1064 :     return 0;
1065 :     }
1066 :     if (*incx == 1) {
1067 :     goto L20;
1068 :     }
1069 :     /* code for increment not equal to 1 */
1070 :     nincx = *n * *incx;
1071 :     i__1 = nincx;
1072 :     i__2 = *incx;
1073 :     for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
1074 :     dx[i__] = *da * dx[i__];
1075 :     /* L10: */
1076 :     }
1077 :     return 0;
1078 :     /* code for increment equal to 1
1079 :     clean-up loop */
1080 :     L20:
1081 :     m = *n % 5;
1082 :     if (m == 0) {
1083 :     goto L40;
1084 :     }
1085 :     i__2 = m;
1086 :     for (i__ = 1; i__ <= i__2; ++i__) {
1087 :     dx[i__] = *da * dx[i__];
1088 :     /* L30: */
1089 :     }
1090 :     if (*n < 5) {
1091 :     return 0;
1092 :     }
1093 :     L40:
1094 :     mp1 = m + 1;
1095 :     i__2 = *n;
1096 :     for (i__ = mp1; i__ <= i__2; i__ += 5) {
1097 :     dx[i__] = *da * dx[i__];
1098 :     dx[i__ + 1] = *da * dx[i__ + 1];
1099 :     dx[i__ + 2] = *da * dx[i__ + 2];
1100 :     dx[i__ + 3] = *da * dx[i__ + 3];
1101 :     dx[i__ + 4] = *da * dx[i__ + 4];
1102 :     /* L50: */
1103 :     }
1104 :     return 0;
1105 :     } /* dscal_ */
1106 :    
1107 :    
1108 :     /* Subroutine */ int dswap_(integer *n, doublereal *dx, integer *incx,
1109 :     doublereal *dy, integer *incy)
1110 :     {
1111 :     /* System generated locals */
1112 :     integer i__1;
1113 :     /* Local variables */
1114 :     static integer i__, m;
1115 :     static doublereal dtemp;
1116 :     static integer ix, iy, mp1;
1117 :     /* interchanges two vectors.
1118 :     uses unrolled loops for increments equal one.
1119 :     jack dongarra, linpack, 3/11/78.
1120 :     modified 12/3/93, array(1) declarations changed to array(*)
1121 :     Parameter adjustments */
1122 :     --dy;
1123 :     --dx;
1124 :     /* Function Body */
1125 :     if (*n <= 0) {
1126 :     return 0;
1127 :     }
1128 :     if (*incx == 1 && *incy == 1) {
1129 :     goto L20;
1130 :     }
1131 :     /* code for unequal increments or equal increments not equal
1132 :     to 1 */
1133 :     ix = 1;
1134 :     iy = 1;
1135 :     if (*incx < 0) {
1136 :     ix = (-(*n) + 1) * *incx + 1;
1137 :     }
1138 :     if (*incy < 0) {
1139 :     iy = (-(*n) + 1) * *incy + 1;
1140 :     }
1141 :     i__1 = *n;
1142 :     for (i__ = 1; i__ <= i__1; ++i__) {
1143 :     dtemp = dx[ix];
1144 :     dx[ix] = dy[iy];
1145 :     dy[iy] = dtemp;
1146 :     ix += *incx;
1147 :     iy += *incy;
1148 :     /* L10: */
1149 :     }
1150 :     return 0;
1151 :     /* code for both increments equal to 1
1152 :     clean-up loop */
1153 :     L20:
1154 :     m = *n % 3;
1155 :     if (m == 0) {
1156 :     goto L40;
1157 :     }
1158 :     i__1 = m;
1159 :     for (i__ = 1; i__ <= i__1; ++i__) {
1160 :     dtemp = dx[i__];
1161 :     dx[i__] = dy[i__];
1162 :     dy[i__] = dtemp;
1163 :     /* L30: */
1164 :     }
1165 :     if (*n < 3) {
1166 :     return 0;
1167 :     }
1168 :     L40:
1169 :     mp1 = m + 1;
1170 :     i__1 = *n;
1171 :     for (i__ = mp1; i__ <= i__1; i__ += 3) {
1172 :     dtemp = dx[i__];
1173 :     dx[i__] = dy[i__];
1174 :     dy[i__] = dtemp;
1175 :     dtemp = dx[i__ + 1];
1176 :     dx[i__ + 1] = dy[i__ + 1];
1177 :     dy[i__ + 1] = dtemp;
1178 :     dtemp = dx[i__ + 2];
1179 :     dx[i__ + 2] = dy[i__ + 2];
1180 :     dy[i__ + 2] = dtemp;
1181 :     /* L50: */
1182 :     }
1183 :     return 0;
1184 :     } /* dswap_ */
1185 :    
1186 :    
1187 :     /* Subroutine */ int dsymv_(char *uplo, integer *n, doublereal *alpha,
1188 :     doublereal *a, integer *lda, doublereal *x, integer *incx, doublereal
1189 :     *beta, doublereal *y, integer *incy)
1190 :     {
1191 :     /* System generated locals */
1192 :     integer a_dim1, a_offset, i__1, i__2;
1193 :     /* Local variables */
1194 :     static integer info;
1195 :     static doublereal temp1, temp2;
1196 :     static integer i__, j;
1197 :     extern logical lsame_(char *, char *);
1198 :     static integer ix, iy, jx, jy, kx, ky;
1199 :     extern /* Subroutine */ int xerbla_(char *, integer *);
1200 :     #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
1201 :     /* Purpose
1202 :     =======
1203 :     DSYMV performs the matrix-vector operation
1204 :     y := alpha*A*x + beta*y,
1205 :     where alpha and beta are scalars, x and y are n element vectors and
1206 :     A is an n by n symmetric matrix.
1207 :     Parameters
1208 :     ==========
1209 :     UPLO - CHARACTER*1.
1210 :     On entry, UPLO specifies whether the upper or lower
1211 :     triangular part of the array A is to be referenced as
1212 :     follows:
1213 :     UPLO = 'U' or 'u' Only the upper triangular part of A
1214 :     is to be referenced.
1215 :     UPLO = 'L' or 'l' Only the lower triangular part of A
1216 :     is to be referenced.
1217 :     Unchanged on exit.
1218 :     N - INTEGER.
1219 :     On entry, N specifies the order of the matrix A.
1220 :     N must be at least zero.
1221 :     Unchanged on exit.
1222 :     ALPHA - DOUBLE PRECISION.
1223 :     On entry, ALPHA specifies the scalar alpha.
1224 :     Unchanged on exit.
1225 :     A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
1226 :     Before entry with UPLO = 'U' or 'u', the leading n by n
1227 :     upper triangular part of the array A must contain the upper
1228 :     triangular part of the symmetric matrix and the strictly
1229 :     lower triangular part of A is not referenced.
1230 :     Before entry with UPLO = 'L' or 'l', the leading n by n
1231 :     lower triangular part of the array A must contain the lower
1232 :     triangular part of the symmetric matrix and the strictly
1233 :     upper triangular part of A is not referenced.
1234 :     Unchanged on exit.
1235 :     LDA - INTEGER.
1236 :     On entry, LDA specifies the first dimension of A as declared
1237 :     in the calling (sub) program. LDA must be at least
1238 :     max( 1, n ).
1239 :     Unchanged on exit.
1240 :     X - DOUBLE PRECISION array of dimension at least
1241 :     ( 1 + ( n - 1 )*abs( INCX ) ).
1242 :     Before entry, the incremented array X must contain the n
1243 :     element vector x.
1244 :     Unchanged on exit.
1245 :     INCX - INTEGER.
1246 :     On entry, INCX specifies the increment for the elements of
1247 :     X. INCX must not be zero.
1248 :     Unchanged on exit.
1249 :     BETA - DOUBLE PRECISION.
1250 :     On entry, BETA specifies the scalar beta. When BETA is
1251 :     supplied as zero then Y need not be set on input.
1252 :     Unchanged on exit.
1253 :     Y - DOUBLE PRECISION array of dimension at least
1254 :     ( 1 + ( n - 1 )*abs( INCY ) ).
1255 :     Before entry, the incremented array Y must contain the n
1256 :     element vector y. On exit, Y is overwritten by the updated
1257 :     vector y.
1258 :     INCY - INTEGER.
1259 :     On entry, INCY specifies the increment for the elements of
1260 :     Y. INCY must not be zero.
1261 :     Unchanged on exit.
1262 :     Level 2 Blas routine.
1263 :     -- Written on 22-October-1986.
1264 :     Jack Dongarra, Argonne National Lab.
1265 :     Jeremy Du Croz, Nag Central Office.
1266 :     Sven Hammarling, Nag Central Office.
1267 :     Richard Hanson, Sandia National Labs.
1268 :     Test the input parameters.
1269 :     Parameter adjustments */
1270 :     a_dim1 = *lda;
1271 :     a_offset = 1 + a_dim1 * 1;
1272 :     a -= a_offset;
1273 :     --x;
1274 :     --y;
1275 :     /* Function Body */
1276 :     info = 0;
1277 :     if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
1278 :     info = 1;
1279 :     } else if (*n < 0) {
1280 :     info = 2;
1281 :     } else if (*lda < max(1,*n)) {
1282 :     info = 5;
1283 :     } else if (*incx == 0) {
1284 :     info = 7;
1285 :     } else if (*incy == 0) {
1286 :     info = 10;
1287 :     }
1288 :     if (info != 0) {
1289 :     xerbla_("DSYMV ", &info);
1290 :     return 0;
1291 :     }
1292 :     /* Quick return if possible. */
1293 :     if (*n == 0 || *alpha == 0. && *beta == 1.) {
1294 :     return 0;
1295 :     }
1296 :     /* Set up the start points in X and Y. */
1297 :     if (*incx > 0) {
1298 :     kx = 1;
1299 :     } else {
1300 :     kx = 1 - (*n - 1) * *incx;
1301 :     }
1302 :     if (*incy > 0) {
1303 :     ky = 1;
1304 :     } else {
1305 :     ky = 1 - (*n - 1) * *incy;
1306 :     }
1307 :     /* Start the operations. In this version the elements of A are
1308 :     accessed sequentially with one pass through the triangular part
1309 :     of A.
1310 :     First form y := beta*y. */
1311 :     if (*beta != 1.) {
1312 :     if (*incy == 1) {
1313 :     if (*beta == 0.) {
1314 :     i__1 = *n;
1315 :     for (i__ = 1; i__ <= i__1; ++i__) {
1316 :     y[i__] = 0.;
1317 :     /* L10: */
1318 :     }
1319 :     } else {
1320 :     i__1 = *n;
1321 :     for (i__ = 1; i__ <= i__1; ++i__) {
1322 :     y[i__] = *beta * y[i__];
1323 :     /* L20: */
1324 :     }
1325 :     }
1326 :     } else {
1327 :     iy = ky;
1328 :     if (*beta == 0.) {
1329 :     i__1 = *n;
1330 :     for (i__ = 1; i__ <= i__1; ++i__) {
1331 :     y[iy] = 0.;
1332 :     iy += *incy;
1333 :     /* L30: */
1334 :     }
1335 :     } else {
1336 :     i__1 = *n;
1337 :     for (i__ = 1; i__ <= i__1; ++i__) {
1338 :     y[iy] = *beta * y[iy];
1339 :     iy += *incy;
1340 :     /* L40: */
1341 :     }
1342 :     }
1343 :     }
1344 :     }
1345 :     if (*alpha == 0.) {
1346 :     return 0;
1347 :     }
1348 :     if (lsame_(uplo, "U")) {
1349 :     /* Form y when A is stored in upper triangle. */
1350 :     if (*incx == 1 && *incy == 1) {
1351 :     i__1 = *n;
1352 :     for (j = 1; j <= i__1; ++j) {
1353 :     temp1 = *alpha * x[j];
1354 :     temp2 = 0.;
1355 :     i__2 = j - 1;
1356 :     for (i__ = 1; i__ <= i__2; ++i__) {
1357 :     y[i__] += temp1 * a_ref(i__, j);
1358 :     temp2 += a_ref(i__, j) * x[i__];
1359 :     /* L50: */
1360 :     }
1361 :     y[j] = y[j] + temp1 * a_ref(j, j) + *alpha * temp2;
1362 :     /* L60: */
1363 :     }
1364 :     } else {
1365 :     jx = kx;
1366 :     jy = ky;
1367 :     i__1 = *n;
1368 :     for (j = 1; j <= i__1; ++j) {
1369 :     temp1 = *alpha * x[jx];
1370 :     temp2 = 0.;
1371 :     ix = kx;
1372 :     iy = ky;
1373 :     i__2 = j - 1;
1374 :     for (i__ = 1; i__ <= i__2; ++i__) {
1375 :     y[iy] += temp1 * a_ref(i__, j);
1376 :     temp2 += a_ref(i__, j) * x[ix];
1377 :     ix += *incx;
1378 :     iy += *incy;
1379 :     /* L70: */
1380 :     }
1381 :     y[jy] = y[jy] + temp1 * a_ref(j, j) + *alpha * temp2;
1382 :     jx += *incx;
1383 :     jy += *incy;
1384 :     /* L80: */
1385 :     }
1386 :     }
1387 :     } else {
1388 :     /* Form y when A is stored in lower triangle. */
1389 :     if (*incx == 1 && *incy == 1) {
1390 :     i__1 = *n;
1391 :     for (j = 1; j <= i__1; ++j) {
1392 :     temp1 = *alpha * x[j];
1393 :     temp2 = 0.;
1394 :     y[j] += temp1 * a_ref(j, j);
1395 :     i__2 = *n;
1396 :     for (i__ = j + 1; i__ <= i__2; ++i__) {
1397 :     y[i__] += temp1 * a_ref(i__, j);
1398 :     temp2 += a_ref(i__, j) * x[i__];
1399 :     /* L90: */
1400 :     }
1401 :     y[j] += *alpha * temp2;
1402 :     /* L100: */
1403 :     }
1404 :     } else {
1405 :     jx = kx;
1406 :     jy = ky;
1407 :     i__1 = *n;
1408 :     for (j = 1; j <= i__1; ++j) {
1409 :     temp1 = *alpha * x[jx];
1410 :     temp2 = 0.;
1411 :     y[jy] += temp1 * a_ref(j, j);
1412 :     ix = jx;
1413 :     iy = jy;
1414 :     i__2 = *n;
1415 :     for (i__ = j + 1; i__ <= i__2; ++i__) {
1416 :     ix += *incx;
1417 :     iy += *incy;
1418 :     y[iy] += temp1 * a_ref(i__, j);
1419 :     temp2 += a_ref(i__, j) * x[ix];
1420 :     /* L110: */
1421 :     }
1422 :     y[jy] += *alpha * temp2;
1423 :     jx += *incx;
1424 :     jy += *incy;
1425 :     /* L120: */
1426 :     }
1427 :     }
1428 :     }
1429 :     return 0;
1430 :     /* End of DSYMV . */
1431 :     } /* dsymv_ */
1432 :     #undef a_ref
1433 :    
1434 :    
1435 :     /* Subroutine */ int dsyr2_(char *uplo, integer *n, doublereal *alpha,
1436 :     doublereal *x, integer *incx, doublereal *y, integer *incy,
1437 :     doublereal *a, integer *lda)
1438 :     {
1439 :     /* System generated locals */
1440 :     integer a_dim1, a_offset, i__1, i__2;
1441 :     /* Local variables */
1442 :     static integer info;
1443 :     static doublereal temp1, temp2;
1444 :     static integer i__, j;
1445 :     extern logical lsame_(char *, char *);
1446 :     static integer ix, iy, jx, jy, kx, ky;
1447 :     extern /* Subroutine */ int xerbla_(char *, integer *);
1448 :     #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
1449 :     /* Purpose
1450 :     =======
1451 :     DSYR2 performs the symmetric rank 2 operation
1452 :     A := alpha*x*y' + alpha*y*x' + A,
1453 :     where alpha is a scalar, x and y are n element vectors and A is an n
1454 :     by n symmetric matrix.
1455 :     Parameters
1456 :     ==========
1457 :     UPLO - CHARACTER*1.
1458 :     On entry, UPLO specifies whether the upper or lower
1459 :     triangular part of the array A is to be referenced as
1460 :     follows:
1461 :     UPLO = 'U' or 'u' Only the upper triangular part of A
1462 :     is to be referenced.
1463 :     UPLO = 'L' or 'l' Only the lower triangular part of A
1464 :     is to be referenced.
1465 :     Unchanged on exit.
1466 :     N - INTEGER.
1467 :     On entry, N specifies the order of the matrix A.
1468 :     N must be at least zero.
1469 :     Unchanged on exit.
1470 :     ALPHA - DOUBLE PRECISION.
1471 :     On entry, ALPHA specifies the scalar alpha.
1472 :     Unchanged on exit.
1473 :     X - DOUBLE PRECISION array of dimension at least
1474 :     ( 1 + ( n - 1 )*abs( INCX ) ).
1475 :     Before entry, the incremented array X must contain the n
1476 :     element vector x.
1477 :     Unchanged on exit.
1478 :     INCX - INTEGER.
1479 :     On entry, INCX specifies the increment for the elements of
1480 :     X. INCX must not be zero.
1481 :     Unchanged on exit.
1482 :     Y - DOUBLE PRECISION array of dimension at least
1483 :     ( 1 + ( n - 1 )*abs( INCY ) ).
1484 :     Before entry, the incremented array Y must contain the n
1485 :     element vector y.
1486 :     Unchanged on exit.
1487 :     INCY - INTEGER.
1488 :     On entry, INCY specifies the increment for the elements of
1489 :     Y. INCY must not be zero.
1490 :     Unchanged on exit.
1491 :     A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
1492 :     Before entry with UPLO = 'U' or 'u', the leading n by n
1493 :     upper triangular part of the array A must contain the upper
1494 :     triangular part of the symmetric matrix and the strictly
1495 :     lower triangular part of A is not referenced. On exit, the
1496 :     upper triangular part of the array A is overwritten by the
1497 :     upper triangular part of the updated matrix.
1498 :     Before entry with UPLO = 'L' or 'l', the leading n by n
1499 :     lower triangular part of the array A must contain the lower
1500 :     triangular part of the symmetric matrix and the strictly
1501 :     upper triangular part of A is not referenced. On exit, the
1502 :     lower triangular part of the array A is overwritten by the
1503 :     lower triangular part of the updated matrix.
1504 :     LDA - INTEGER.
1505 :     On entry, LDA specifies the first dimension of A as declared
1506 :     in the calling (sub) program. LDA must be at least
1507 :     max( 1, n ).
1508 :     Unchanged on exit.
1509 :     Level 2 Blas routine.
1510 :     -- Written on 22-October-1986.
1511 :     Jack Dongarra, Argonne National Lab.
1512 :     Jeremy Du Croz, Nag Central Office.
1513 :     Sven Hammarling, Nag Central Office.
1514 :     Richard Hanson, Sandia National Labs.
1515 :     Test the input parameters.
1516 :     Parameter adjustments */
1517 :     --x;
1518 :     --y;
1519 :     a_dim1 = *lda;
1520 :     a_offset = 1 + a_dim1 * 1;
1521 :     a -= a_offset;
1522 :     /* Function Body */
1523 :     info = 0;
1524 :     if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
1525 :     info = 1;
1526 :     } else if (*n < 0) {
1527 :     info = 2;
1528 :     } else if (*incx == 0) {
1529 :     info = 5;
1530 :     } else if (*incy == 0) {
1531 :     info = 7;
1532 :     } else if (*lda < max(1,*n)) {
1533 :     info = 9;
1534 :     }
1535 :     if (info != 0) {
1536 :     xerbla_("DSYR2 ", &info);
1537 :     return 0;
1538 :     }
1539 :     /* Quick return if possible. */
1540 :     if (*n == 0 || *alpha == 0.) {
1541 :     return 0;
1542 :     }
1543 :     /* Set up the start points in X and Y if the increments are not both
1544 :     unity. */
1545 :     if (*incx != 1 || *incy != 1) {
1546 :     if (*incx > 0) {
1547 :     kx = 1;
1548 :     } else {
1549 :     kx = 1 - (*n - 1) * *incx;
1550 :     }
1551 :     if (*incy > 0) {
1552 :     ky = 1;
1553 :     } else {
1554 :     ky = 1 - (*n - 1) * *incy;
1555 :     }
1556 :     jx = kx;
1557 :     jy = ky;
1558 :     }
1559 :     /* Start the operations. In this version the elements of A are
1560 :     accessed sequentially with one pass through the triangular part
1561 :     of A. */
1562 :     if (lsame_(uplo, "U")) {
1563 :     /* Form A when A is stored in the upper triangle. */
1564 :     if (*incx == 1 && *incy == 1) {
1565 :     i__1 = *n;
1566 :     for (j = 1; j <= i__1; ++j) {
1567 :     if (x[j] != 0. || y[j] != 0.) {
1568 :     temp1 = *alpha * y[j];
1569 :     temp2 = *alpha * x[j];
1570 :     i__2 = j;
1571 :     for (i__ = 1; i__ <= i__2; ++i__) {
1572 :     a_ref(i__, j) = a_ref(i__, j) + x[i__] * temp1 + y[
1573 :     i__] * temp2;
1574 :     /* L10: */
1575 :     }
1576 :     }
1577 :     /* L20: */
1578 :     }
1579 :     } else {
1580 :     i__1 = *n;
1581 :     for (j = 1; j <= i__1; ++j) {
1582 :     if (x[jx] != 0. || y[jy] != 0.) {
1583 :     temp1 = *alpha * y[jy];
1584 :     temp2 = *alpha * x[jx];
1585 :     ix = kx;
1586 :     iy = ky;
1587 :     i__2 = j;
1588 :     for (i__ = 1; i__ <= i__2; ++i__) {
1589 :     a_ref(i__, j) = a_ref(i__, j) + x[ix] * temp1 + y[iy]
1590 :     * temp2;
1591 :     ix += *incx;
1592 :     iy += *incy;
1593 :     /* L30: */
1594 :     }
1595 :     }
1596 :     jx += *incx;
1597 :     jy += *incy;
1598 :     /* L40: */
1599 :     }
1600 :     }
1601 :     } else {
1602 :     /* Form A when A is stored in the lower triangle. */
1603 :     if (*incx == 1 && *incy == 1) {
1604 :     i__1 = *n;
1605 :     for (j = 1; j <= i__1; ++j) {
1606 :     if (x[j] != 0. || y[j] != 0.) {
1607 :     temp1 = *alpha * y[j];
1608 :     temp2 = *alpha * x[j];
1609 :     i__2 = *n;
1610 :     for (i__ = j; i__ <= i__2; ++i__) {
1611 :     a_ref(i__, j) = a_ref(i__, j) + x[i__] * temp1 + y[
1612 :     i__] * temp2;
1613 :     /* L50: */
1614 :     }
1615 :     }
1616 :     /* L60: */
1617 :     }
1618 :     } else {
1619 :     i__1 = *n;
1620 :     for (j = 1; j <= i__1; ++j) {
1621 :     if (x[jx] != 0. || y[jy] != 0.) {
1622 :     temp1 = *alpha * y[jy];
1623 :     temp2 = *alpha * x[jx];
1624 :     ix = jx;
1625 :     iy = jy;
1626 :     i__2 = *n;
1627 :     for (i__ = j; i__ <= i__2; ++i__) {
1628 :     a_ref(i__, j) = a_ref(i__, j) + x[ix] * temp1 + y[iy]
1629 :     * temp2;
1630 :     ix += *incx;
1631 :     iy += *incy;
1632 :     /* L70: */
1633 :     }
1634 :     }
1635 :     jx += *incx;
1636 :     jy += *incy;
1637 :     /* L80: */
1638 :     }
1639 :     }
1640 :     }
1641 :     return 0;
1642 :     /* End of DSYR2 . */
1643 :     } /* dsyr2_ */
1644 :     #undef a_ref
1645 :    
1646 :    
1647 :     /* Subroutine */ int dsyr2k_(char *uplo, char *trans, integer *n, integer *k,
1648 :     doublereal *alpha, doublereal *a, integer *lda, doublereal *b,
1649 :     integer *ldb, doublereal *beta, doublereal *c__, integer *ldc)
1650 :     {
1651 :     /* System generated locals */
1652 :     integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
1653 :     i__3;
1654 :     /* Local variables */
1655 :     static integer info;
1656 :     static doublereal temp1, temp2;
1657 :     static integer i__, j, l;
1658 :     extern logical lsame_(char *, char *);
1659 :     static integer nrowa;
1660 :     static logical upper;
1661 :     extern /* Subroutine */ int xerbla_(char *, integer *);
1662 :     #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
1663 :     #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
1664 :     #define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1]
1665 :     /* Purpose
1666 :     =======
1667 :     DSYR2K performs one of the symmetric rank 2k operations
1668 :     C := alpha*A*B' + alpha*B*A' + beta*C,
1669 :     or
1670 :     C := alpha*A'*B + alpha*B'*A + beta*C,
1671 :     where alpha and beta are scalars, C is an n by n symmetric matrix
1672 :     and A and B are n by k matrices in the first case and k by n
1673 :     matrices in the second case.
1674 :     Parameters
1675 :     ==========
1676 :     UPLO - CHARACTER*1.
1677 :     On entry, UPLO specifies whether the upper or lower
1678 :     triangular part of the array C is to be referenced as
1679 :     follows:
1680 :     UPLO = 'U' or 'u' Only the upper triangular part of C
1681 :     is to be referenced.
1682 :     UPLO = 'L' or 'l' Only the lower triangular part of C
1683 :     is to be referenced.
1684 :     Unchanged on exit.
1685 :     TRANS - CHARACTER*1.
1686 :     On entry, TRANS specifies the operation to be performed as
1687 :     follows:
1688 :     TRANS = 'N' or 'n' C := alpha*A*B' + alpha*B*A' +
1689 :     beta*C.
1690 :     TRANS = 'T' or 't' C := alpha*A'*B + alpha*B'*A +
1691 :     beta*C.
1692 :     TRANS = 'C' or 'c' C := alpha*A'*B + alpha*B'*A +
1693 :     beta*C.
1694 :     Unchanged on exit.
1695 :     N - INTEGER.
1696 :     On entry, N specifies the order of the matrix C. N must be
1697 :     at least zero.
1698 :     Unchanged on exit.
1699 :     K - INTEGER.
1700 :     On entry with TRANS = 'N' or 'n', K specifies the number
1701 :     of columns of the matrices A and B, and on entry with
1702 :     TRANS = 'T' or 't' or 'C' or 'c', K specifies the number
1703 :     of rows of the matrices A and B. K must be at least zero.
1704 :     Unchanged on exit.
1705 :     ALPHA - DOUBLE PRECISION.
1706 :     On entry, ALPHA specifies the scalar alpha.
1707 :     Unchanged on exit.
1708 :     A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
1709 :     k when TRANS = 'N' or 'n', and is n otherwise.
1710 :     Before entry with TRANS = 'N' or 'n', the leading n by k
1711 :     part of the array A must contain the matrix A, otherwise
1712 :     the leading k by n part of the array A must contain the
1713 :     matrix A.
1714 :     Unchanged on exit.
1715 :     LDA - INTEGER.
1716 :     On entry, LDA specifies the first dimension of A as declared
1717 :     in the calling (sub) program. When TRANS = 'N' or 'n'
1718 :     then LDA must be at least max( 1, n ), otherwise LDA must
1719 :     be at least max( 1, k ).
1720 :     Unchanged on exit.
1721 :     B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
1722 :     k when TRANS = 'N' or 'n', and is n otherwise.
1723 :     Before entry with TRANS = 'N' or 'n', the leading n by k
1724 :     part of the array B must contain the matrix B, otherwise
1725 :     the leading k by n part of the array B must contain the
1726 :     matrix B.
1727 :     Unchanged on exit.
1728 :     LDB - INTEGER.
1729 :     On entry, LDB specifies the first dimension of B as declared
1730 :     in the calling (sub) program. When TRANS = 'N' or 'n'
1731 :     then LDB must be at least max( 1, n ), otherwise LDB must
1732 :     be at least max( 1, k ).
1733 :     Unchanged on exit.
1734 :     BETA - DOUBLE PRECISION.
1735 :     On entry, BETA specifies the scalar beta.
1736 :     Unchanged on exit.
1737 :     C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
1738 :     Before entry with UPLO = 'U' or 'u', the leading n by n
1739 :     upper triangular part of the array C must contain the upper
1740 :     triangular part of the symmetric matrix and the strictly
1741 :     lower triangular part of C is not referenced. On exit, the
1742 :     upper triangular part of the array C is overwritten by the
1743 :     upper triangular part of the updated matrix.
1744 :     Before entry with UPLO = 'L' or 'l', the leading n by n
1745 :     lower triangular part of the array C must contain the lower
1746 :     triangular part of the symmetric matrix and the strictly
1747 :     upper triangular part of C is not referenced. On exit, the
1748 :     lower triangular part of the array C is overwritten by the
1749 :     lower triangular part of the updated matrix.
1750 :     LDC - INTEGER.
1751 :     On entry, LDC specifies the first dimension of C as declared
1752 :     in the calling (sub) program. LDC must be at least
1753 :     max( 1, n ).
1754 :     Unchanged on exit.
1755 :     Level 3 Blas routine.
1756 :     -- Written on 8-February-1989.
1757 :     Jack Dongarra, Argonne National Laboratory.
1758 :     Iain Duff, AERE Harwell.
1759 :     Jeremy Du Croz, Numerical Algorithms Group Ltd.
1760 :     Sven Hammarling, Numerical Algorithms Group Ltd.
1761 :     Test the input parameters.
1762 :     Parameter adjustments */
1763 :     a_dim1 = *lda;
1764 :     a_offset = 1 + a_dim1 * 1;
1765 :     a -= a_offset;
1766 :     b_dim1 = *ldb;
1767 :     b_offset = 1 + b_dim1 * 1;
1768 :     b -= b_offset;
1769 :     c_dim1 = *ldc;
1770 :     c_offset = 1 + c_dim1 * 1;
1771 :     c__ -= c_offset;
1772 :     /* Function Body */
1773 :     if (lsame_(trans, "N")) {
1774 :     nrowa = *n;
1775 :     } else {
1776 :     nrowa = *k;
1777 :     }
1778 :     upper = lsame_(uplo, "U");
1779 :     info = 0;
1780 :     if (! upper && ! lsame_(uplo, "L")) {
1781 :     info = 1;
1782 :     } else if (! lsame_(trans, "N") && ! lsame_(trans,
1783 :     "T") && ! lsame_(trans, "C")) {
1784 :     info = 2;
1785 :     } else if (*n < 0) {
1786 :     info = 3;
1787 :     } else if (*k < 0) {
1788 :     info = 4;
1789 :     } else if (*lda < max(1,nrowa)) {
1790 :     info = 7;
1791 :     } else if (*ldb < max(1,nrowa)) {
1792 :     info = 9;
1793 :     } else if (*ldc < max(1,*n)) {
1794 :     info = 12;
1795 :     }
1796 :     if (info != 0) {
1797 :     xerbla_("DSYR2K", &info);
1798 :     return 0;
1799 :     }
1800 :     /* Quick return if possible. */
1801 :     if (*n == 0 || (*alpha == 0. || *k == 0) && *beta == 1.) {
1802 :     return 0;
1803 :     }
1804 :     /* And when alpha.eq.zero. */
1805 :     if (*alpha == 0.) {
1806 :     if (upper) {
1807 :     if (*beta == 0.) {
1808 :     i__1 = *n;
1809 :     for (j = 1; j <= i__1; ++j) {
1810 :     i__2 = j;
1811 :     for (i__ = 1; i__ <= i__2; ++i__) {
1812 :     c___ref(i__, j) = 0.;
1813 :     /* L10: */
1814 :     }
1815 :     /* L20: */
1816 :     }
1817 :     } else {
1818 :     i__1 = *n;
1819 :     for (j = 1; j <= i__1; ++j) {
1820 :     i__2 = j;
1821 :     for (i__ = 1; i__ <= i__2; ++i__) {
1822 :     c___ref(i__, j) = *beta * c___ref(i__, j);
1823 :     /* L30: */
1824 :     }
1825 :     /* L40: */
1826 :     }
1827 :     }
1828 :     } else {
1829 :     if (*beta == 0.) {
1830 :     i__1 = *n;
1831 :     for (j = 1; j <= i__1; ++j) {
1832 :     i__2 = *n;
1833 :     for (i__ = j; i__ <= i__2; ++i__) {
1834 :     c___ref(i__, j) = 0.;
1835 :     /* L50: */
1836 :     }
1837 :     /* L60: */
1838 :     }
1839 :     } else {
1840 :     i__1 = *n;
1841 :     for (j = 1; j <= i__1; ++j) {
1842 :     i__2 = *n;
1843 :     for (i__ = j; i__ <= i__2; ++i__) {
1844 :     c___ref(i__, j) = *beta * c___ref(i__, j);
1845 :     /* L70: */
1846 :     }
1847 :     /* L80: */
1848 :     }
1849 :     }
1850 :     }
1851 :     return 0;
1852 :     }
1853 :     /* Start the operations. */
1854 :     if (lsame_(trans, "N")) {
1855 :     /* Form C := alpha*A*B' + alpha*B*A' + C. */
1856 :     if (upper) {
1857 :     i__1 = *n;
1858 :     for (j = 1; j <= i__1; ++j) {
1859 :     if (*beta == 0.) {
1860 :     i__2 = j;
1861 :     for (i__ = 1; i__ <= i__2; ++i__) {
1862 :     c___ref(i__, j) = 0.;
1863 :     /* L90: */
1864 :     }
1865 :     } else if (*beta != 1.) {
1866 :     i__2 = j;
1867 :     for (i__ = 1; i__ <= i__2; ++i__) {
1868 :     c___ref(i__, j) = *beta * c___ref(i__, j);
1869 :     /* L100: */
1870 :     }
1871 :     }
1872 :     i__2 = *k;
1873 :     for (l = 1; l <= i__2; ++l) {
1874 :     if (a_ref(j, l) != 0. || b_ref(j, l) != 0.) {
1875 :     temp1 = *alpha * b_ref(j, l);
1876 :     temp2 = *alpha * a_ref(j, l);
1877 :     i__3 = j;
1878 :     for (i__ = 1; i__ <= i__3; ++i__) {
1879 :     c___ref(i__, j) = c___ref(i__, j) + a_ref(i__, l)
1880 :     * temp1 + b_ref(i__, l) * temp2;
1881 :     /* L110: */
1882 :     }
1883 :     }
1884 :     /* L120: */
1885 :     }
1886 :     /* L130: */
1887 :     }
1888 :     } else {
1889 :     i__1 = *n;
1890 :     for (j = 1; j <= i__1; ++j) {
1891 :     if (*beta == 0.) {
1892 :     i__2 = *n;
1893 :     for (i__ = j; i__ <= i__2; ++i__) {
1894 :     c___ref(i__, j) = 0.;
1895 :     /* L140: */
1896 :     }
1897 :     } else if (*beta != 1.) {
1898 :     i__2 = *n;
1899 :     for (i__ = j; i__ <= i__2; ++i__) {
1900 :     c___ref(i__, j) = *beta * c___ref(i__, j);
1901 :     /* L150: */
1902 :     }
1903 :     }
1904 :     i__2 = *k;
1905 :     for (l = 1; l <= i__2; ++l) {
1906 :     if (a_ref(j, l) != 0. || b_ref(j, l) != 0.) {
1907 :     temp1 = *alpha * b_ref(j, l);
1908 :     temp2 = *alpha * a_ref(j, l);
1909 :     i__3 = *n;
1910 :     for (i__ = j; i__ <= i__3; ++i__) {
1911 :     c___ref(i__, j) = c___ref(i__, j) + a_ref(i__, l)
1912 :     * temp1 + b_ref(i__, l) * temp2;
1913 :     /* L160: */
1914 :     }
1915 :     }
1916 :     /* L170: */
1917 :     }
1918 :     /* L180: */
1919 :     }
1920 :     }
1921 :     } else {
1922 :     /* Form C := alpha*A'*B + alpha*B'*A + C. */
1923 :     if (upper) {
1924 :     i__1 = *n;
1925 :     for (j = 1; j <= i__1; ++j) {
1926 :     i__2 = j;
1927 :     for (i__ = 1; i__ <= i__2; ++i__) {
1928 :     temp1 = 0.;
1929 :     temp2 = 0.;
1930 :     i__3 = *k;
1931 :     for (l = 1; l <= i__3; ++l) {
1932 :     temp1 += a_ref(l, i__) * b_ref(l, j);
1933 :     temp2 += b_ref(l, i__) * a_ref(l, j);
1934 :     /* L190: */
1935 :     }
1936 :     if (*beta == 0.) {
1937 :     c___ref(i__, j) = *alpha * temp1 + *alpha * temp2;
1938 :     } else {
1939 :     c___ref(i__, j) = *beta * c___ref(i__, j) + *alpha *
1940 :     temp1 + *alpha * temp2;
1941 :     }
1942 :     /* L200: */
1943 :     }
1944 :     /* L210: */
1945 :     }
1946 :     } else {
1947 :     i__1 = *n;
1948 :     for (j = 1; j <= i__1; ++j) {
1949 :     i__2 = *n;
1950 :     for (i__ = j; i__ <= i__2; ++i__) {
1951 :     temp1 = 0.;
1952 :     temp2 = 0.;
1953 :     i__3 = *k;
1954 :     for (l = 1; l <= i__3; ++l) {
1955 :     temp1 += a_ref(l, i__) * b_ref(l, j);
1956 :     temp2 += b_ref(l, i__) * a_ref(l, j);
1957 :     /* L220: */
1958 :     }
1959 :     if (*beta == 0.) {
1960 :     c___ref(i__, j) = *alpha * temp1 + *alpha * temp2;
1961 :     } else {
1962 :     c___ref(i__, j) = *beta * c___ref(i__, j) + *alpha *
1963 :     temp1 + *alpha * temp2;
1964 :     }
1965 :     /* L230: */
1966 :     }
1967 :     /* L240: */
1968 :     }
1969 :     }
1970 :     }
1971 :     return 0;
1972 :     /* End of DSYR2K. */
1973 :     } /* dsyr2k_ */
1974 :     #undef c___ref
1975 :     #undef b_ref
1976 :     #undef a_ref
1977 :    
1978 :    
1979 :     /* Subroutine */ int dtrmm_(char *side, char *uplo, char *transa, char *diag,
1980 :     integer *m, integer *n, doublereal *alpha, doublereal *a, integer *
1981 :     lda, doublereal *b, integer *ldb)
1982 :     {
1983 :     /* System generated locals */
1984 :     integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
1985 :     /* Local variables */
1986 :     static integer info;
1987 :     static doublereal temp;
1988 :     static integer i__, j, k;
1989 :     static logical lside;
1990 :     extern logical lsame_(char *, char *);
1991 :     static integer nrowa;
1992 :     static logical upper;
1993 :     extern /* Subroutine */ int xerbla_(char *, integer *);
1994 :     static logical nounit;
1995 :     #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
1996 :     #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
1997 :     /* Purpose
1998 :     =======
1999 :     DTRMM performs one of the matrix-matrix operations
2000 :     B := alpha*op( A )*B, or B := alpha*B*op( A ),
2001 :     where alpha is a scalar, B is an m by n matrix, A is a unit, or
2002 :     non-unit, upper or lower triangular matrix and op( A ) is one of
2003 :     op( A ) = A or op( A ) = A'.
2004 :     Parameters
2005 :     ==========
2006 :     SIDE - CHARACTER*1.
2007 :     On entry, SIDE specifies whether op( A ) multiplies B from
2008 :     the left or right as follows:
2009 :     SIDE = 'L' or 'l' B := alpha*op( A )*B.
2010 :     SIDE = 'R' or 'r' B := alpha*B*op( A ).
2011 :     Unchanged on exit.
2012 :     UPLO - CHARACTER*1.
2013 :     On entry, UPLO specifies whether the matrix A is an upper or
2014 :     lower triangular matrix as follows:
2015 :     UPLO = 'U' or 'u' A is an upper triangular matrix.
2016 :     UPLO = 'L' or 'l' A is a lower triangular matrix.
2017 :     Unchanged on exit.
2018 :     TRANSA - CHARACTER*1.
2019 :     On entry, TRANSA specifies the form of op( A ) to be used in
2020 :     the matrix multiplication as follows:
2021 :     TRANSA = 'N' or 'n' op( A ) = A.
2022 :     TRANSA = 'T' or 't' op( A ) = A'.
2023 :     TRANSA = 'C' or 'c' op( A ) = A'.
2024 :     Unchanged on exit.
2025 :     DIAG - CHARACTER*1.
2026 :     On entry, DIAG specifies whether or not A is unit triangular
2027 :     as follows:
2028 :     DIAG = 'U' or 'u' A is assumed to be unit triangular.
2029 :     DIAG = 'N' or 'n' A is not assumed to be unit
2030 :     triangular.
2031 :     Unchanged on exit.
2032 :     M - INTEGER.
2033 :     On entry, M specifies the number of rows of B. M must be at
2034 :     least zero.
2035 :     Unchanged on exit.
2036 :     N - INTEGER.
2037 :     On entry, N specifies the number of columns of B. N must be
2038 :     at least zero.
2039 :     Unchanged on exit.
2040 :     ALPHA - DOUBLE PRECISION.
2041 :     On entry, ALPHA specifies the scalar alpha. When alpha is
2042 :     zero then A is not referenced and B need not be set before
2043 :     entry.
2044 :     Unchanged on exit.
2045 :     A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
2046 :     when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
2047 :     Before entry with UPLO = 'U' or 'u', the leading k by k
2048 :     upper triangular part of the array A must contain the upper
2049 :     triangular matrix and the strictly lower triangular part of
2050 :     A is not referenced.
2051 :     Before entry with UPLO = 'L' or 'l', the leading k by k
2052 :     lower triangular part of the array A must contain the lower
2053 :     triangular matrix and the strictly upper triangular part of
2054 :     A is not referenced.
2055 :     Note that when DIAG = 'U' or 'u', the diagonal elements of
2056 :     A are not referenced either, but are assumed to be unity.
2057 :     Unchanged on exit.
2058 :     LDA - INTEGER.
2059 :     On entry, LDA specifies the first dimension of A as declared
2060 :     in the calling (sub) program. When SIDE = 'L' or 'l' then
2061 :     LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
2062 :     then LDA must be at least max( 1, n ).
2063 :     Unchanged on exit.
2064 :     B - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
2065 :     Before entry, the leading m by n part of the array B must
2066 :     contain the matrix B, and on exit is overwritten by the
2067 :     transformed matrix.
2068 :     LDB - INTEGER.
2069 :     On entry, LDB specifies the first dimension of B as declared
2070 :     in the calling (sub) program. LDB must be at least
2071 :     max( 1, m ).
2072 :     Unchanged on exit.
2073 :     Level 3 Blas routine.
2074 :     -- Written on 8-February-1989.
2075 :     Jack Dongarra, Argonne National Laboratory.
2076 :     Iain Duff, AERE Harwell.
2077 :     Jeremy Du Croz, Numerical Algorithms Group Ltd.
2078 :     Sven Hammarling, Numerical Algorithms Group Ltd.
2079 :     Test the input parameters.
2080 :     Parameter adjustments */
2081 :     a_dim1 = *lda;
2082 :     a_offset = 1 + a_dim1 * 1;
2083 :     a -= a_offset;
2084 :     b_dim1 = *ldb;
2085 :     b_offset = 1 + b_dim1 * 1;
2086 :     b -= b_offset;
2087 :     /* Function Body */
2088 :     lside = lsame_(side, "L");
2089 :     if (lside) {
2090 :     nrowa = *m;
2091 :     } else {
2092 :     nrowa = *n;
2093 :     }
2094 :     nounit = lsame_(diag, "N");
2095 :     upper = lsame_(uplo, "U");
2096 :     info = 0;
2097 :     if (! lside && ! lsame_(side, "R")) {
2098 :     info = 1;
2099 :     } else if (! upper && ! lsame_(uplo, "L")) {
2100 :     info = 2;
2101 :     } else if (! lsame_(transa, "N") && ! lsame_(transa,
2102 :     "T") && ! lsame_(transa, "C")) {
2103 :     info = 3;
2104 :     } else if (! lsame_(diag, "U") && ! lsame_(diag,
2105 :     "N")) {
2106 :     info = 4;
2107 :     } else if (*m < 0) {
2108 :     info = 5;
2109 :     } else if (*n < 0) {
2110 :     info = 6;
2111 :     } else if (*lda < max(1,nrowa)) {
2112 :     info = 9;
2113 :     } else if (*ldb < max(1,*m)) {
2114 :     info = 11;
2115 :     }
2116 :     if (info != 0) {
2117 :     xerbla_("DTRMM ", &info);
2118 :     return 0;
2119 :     }
2120 :     /* Quick return if possible. */
2121 :     if (*n == 0) {
2122 :     return 0;
2123 :     }
2124 :     /* And when alpha.eq.zero. */
2125 :     if (*alpha == 0.) {
2126 :     i__1 = *n;
2127 :     for (j = 1; j <= i__1; ++j) {
2128 :     i__2 = *m;
2129 :     for (i__ = 1; i__ <= i__2; ++i__) {
2130 :     b_ref(i__, j) = 0.;
2131 :     /* L10: */
2132 :     }
2133 :     /* L20: */
2134 :     }
2135 :     return 0;
2136 :     }
2137 :     /* Start the operations. */
2138 :     if (lside) {
2139 :     if (lsame_(transa, "N")) {
2140 :     /* Form B := alpha*A*B. */
2141 :     if (upper) {
2142 :     i__1 = *n;
2143 :     for (j = 1; j <= i__1; ++j) {
2144 :     i__2 = *m;
2145 :     for (k = 1; k <= i__2; ++k) {
2146 :     if (b_ref(k, j) != 0.) {
2147 :     temp = *alpha * b_ref(k, j);
2148 :     i__3 = k - 1;
2149 :     for (i__ = 1; i__ <= i__3; ++i__) {
2150 :     b_ref(i__, j) = b_ref(i__, j) + temp * a_ref(
2151 :     i__, k);
2152 :     /* L30: */
2153 :     }
2154 :     if (nounit) {
2155 :     temp *= a_ref(k, k);
2156 :     }
2157 :     b_ref(k, j) = temp;
2158 :     }
2159 :     /* L40: */
2160 :     }
2161 :     /* L50: */
2162 :     }
2163 :     } else {
2164 :     i__1 = *n;
2165 :     for (j = 1; j <= i__1; ++j) {
2166 :     for (k = *m; k >= 1; --k) {
2167 :     if (b_ref(k, j) != 0.) {
2168 :     temp = *alpha * b_ref(k, j);
2169 :     b_ref(k, j) = temp;
2170 :     if (nounit) {
2171 :     b_ref(k, j) = b_ref(k, j) * a_ref(k, k);
2172 :     }
2173 :     i__2 = *m;
2174 :     for (i__ = k + 1; i__ <= i__2; ++i__) {
2175 :     b_ref(i__, j) = b_ref(i__, j) + temp * a_ref(
2176 :     i__, k);
2177 :     /* L60: */
2178 :     }
2179 :     }
2180 :     /* L70: */
2181 :     }
2182 :     /* L80: */
2183 :     }
2184 :     }
2185 :     } else {
2186 :     /* Form B := alpha*A'*B. */
2187 :     if (upper) {
2188 :     i__1 = *n;
2189 :     for (j = 1; j <= i__1; ++j) {
2190 :     for (i__ = *m; i__ >= 1; --i__) {
2191 :     temp = b_ref(i__, j);
2192 :     if (nounit) {
2193 :     temp *= a_ref(i__, i__);
2194 :     }
2195 :     i__2 = i__ - 1;
2196 :     for (k = 1; k <= i__2; ++k) {
2197 :     temp += a_ref(k, i__) * b_ref(k, j);
2198 :     /* L90: */
2199 :     }
2200 :     b_ref(i__, j) = *alpha * temp;
2201 :     /* L100: */
2202 :     }
2203 :     /* L110: */
2204 :     }
2205 :     } else {
2206 :     i__1 = *n;
2207 :     for (j = 1; j <= i__1; ++j) {
2208 :     i__2 = *m;
2209 :     for (i__ = 1; i__ <= i__2; ++i__) {
2210 :     temp = b_ref(i__, j);
2211 :     if (nounit) {
2212 :     temp *= a_ref(i__, i__);
2213 :     }
2214 :     i__3 = *m;
2215 :     for (k = i__ + 1; k <= i__3; ++k) {
2216 :     temp += a_ref(k, i__) * b_ref(k, j);
2217 :     /* L120: */
2218 :     }
2219 :     b_ref(i__, j) = *alpha * temp;
2220 :     /* L130: */
2221 :     }
2222 :     /* L140: */
2223 :     }
2224 :     }
2225 :     }
2226 :     } else {
2227 :     if (lsame_(transa, "N")) {
2228 :     /* Form B := alpha*B*A. */
2229 :     if (upper) {
2230 :     for (j = *n; j >= 1; --j) {
2231 :     temp = *alpha;
2232 :     if (nounit) {
2233 :     temp *= a_ref(j, j);
2234 :     }
2235 :     i__1 = *m;
2236 :     for (i__ = 1; i__ <= i__1; ++i__) {
2237 :     b_ref(i__, j) = temp * b_ref(i__, j);
2238 :     /* L150: */
2239 :     }
2240 :     i__1 = j - 1;
2241 :     for (k = 1; k <= i__1; ++k) {
2242 :     if (a_ref(k, j) != 0.) {
2243 :     temp = *alpha * a_ref(k, j);
2244 :     i__2 = *m;
2245 :     for (i__ = 1; i__ <= i__2; ++i__) {
2246 :     b_ref(i__, j) = b_ref(i__, j) + temp * b_ref(
2247 :     i__, k);
2248 :     /* L160: */
2249 :     }
2250 :     }
2251 :     /* L170: */
2252 :     }
2253 :     /* L180: */
2254 :     }
2255 :     } else {
2256 :     i__1 = *n;
2257 :     for (j = 1; j <= i__1; ++j) {
2258 :     temp = *alpha;
2259 :     if (nounit) {
2260 :     temp *= a_ref(j, j);
2261 :     }
2262 :     i__2 = *m;
2263 :     for (i__ = 1; i__ <= i__2; ++i__) {
2264 :     b_ref(i__, j) = temp * b_ref(i__, j);
2265 :     /* L190: */
2266 :     }
2267 :     i__2 = *n;
2268 :     for (k = j + 1; k <= i__2; ++k) {
2269 :     if (a_ref(k, j) != 0.) {
2270 :     temp = *alpha * a_ref(k, j);
2271 :     i__3 = *m;
2272 :     for (i__ = 1; i__ <= i__3; ++i__) {
2273 :     b_ref(i__, j) = b_ref(i__, j) + temp * b_ref(
2274 :     i__, k);
2275 :     /* L200: */
2276 :     }
2277 :     }
2278 :     /* L210: */
2279 :     }
2280 :     /* L220: */
2281 :     }
2282 :     }
2283 :     } else {
2284 :     /* Form B := alpha*B*A'. */
2285 :     if (upper) {
2286 :     i__1 = *n;
2287 :     for (k = 1; k <= i__1; ++k) {
2288 :     i__2 = k - 1;
2289 :     for (j = 1; j <= i__2; ++j) {
2290 :     if (a_ref(j, k) != 0.) {
2291 :     temp = *alpha * a_ref(j, k);
2292 :     i__3 = *m;
2293 :     for (i__ = 1; i__ <= i__3; ++i__) {
2294 :     b_ref(i__, j) = b_ref(i__, j) + temp * b_ref(
2295 :     i__, k);
2296 :     /* L230: */
2297 :     }
2298 :     }
2299 :     /* L240: */
2300 :     }
2301 :     temp = *alpha;
2302 :     if (nounit) {
2303 :     temp *= a_ref(k, k);
2304 :     }
2305 :     if (temp != 1.) {
2306 :     i__2 = *m;
2307 :     for (i__ = 1; i__ <= i__2; ++i__) {
2308 :     b_ref(i__, k) = temp * b_ref(i__, k);
2309 :     /* L250: */
2310 :     }
2311 :     }
2312 :     /* L260: */
2313 :     }
2314 :     } else {
2315 :     for (k = *n; k >= 1; --k) {
2316 :     i__1 = *n;
2317 :     for (j = k + 1; j <= i__1; ++j) {
2318 :     if (a_ref(j, k) != 0.) {
2319 :     temp = *alpha * a_ref(j, k);
2320 :     i__2 = *m;
2321 :     for (i__ = 1; i__ <= i__2; ++i__) {
2322 :     b_ref(i__, j) = b_ref(i__, j) + temp * b_ref(
2323 :     i__, k);
2324 :     /* L270: */
2325 :     }
2326 :     }
2327 :     /* L280: */
2328 :     }
2329 :     temp = *alpha;
2330 :     if (nounit) {
2331 :     temp *= a_ref(k, k);
2332 :     }
2333 :     if (temp != 1.) {
2334 :     i__1 = *m;
2335 :     for (i__ = 1; i__ <= i__1; ++i__) {
2336 :     b_ref(i__, k) = temp * b_ref(i__, k);
2337 :     /* L290: */
2338 :     }
2339 :     }
2340 :     /* L300: */
2341 :     }
2342 :     }
2343 :     }
2344 :     }
2345 :     return 0;
2346 :     /* End of DTRMM . */
2347 :     } /* dtrmm_ */
2348 :     #undef b_ref
2349 :     #undef a_ref
2350 :    
2351 :    
2352 :     /* Subroutine */ int dtrmv_(char *uplo, char *trans, char *diag, integer *n,
2353 :     doublereal *a, integer *lda, doublereal *x, integer *incx)
2354 :     {
2355 :     /* System generated locals */
2356 :     integer a_dim1, a_offset, i__1, i__2;
2357 :     /* Local variables */
2358 :     static integer info;
2359 :     static doublereal temp;
2360 :     static integer i__, j;
2361 :     extern logical lsame_(char *, char *);
2362 :     static integer ix, jx, kx;
2363 :     extern /* Subroutine */ int xerbla_(char *, integer *);
2364 :     static logical nounit;
2365 :     #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
2366 :     /* Purpose
2367 :     =======
2368 :     DTRMV performs one of the matrix-vector operations
2369 :     x := A*x, or x := A'*x,
2370 :     where x is an n element vector and A is an n by n unit, or non-unit,
2371 :     upper or lower triangular matrix.
2372 :     Parameters
2373 :     ==========
2374 :     UPLO - CHARACTER*1.
2375 :     On entry, UPLO specifies whether the matrix is an upper or
2376 :     lower triangular matrix as follows:
2377 :     UPLO = 'U' or 'u' A is an upper triangular matrix.
2378 :     UPLO = 'L' or 'l' A is a lower triangular matrix.
2379 :     Unchanged on exit.
2380 :     TRANS - CHARACTER*1.
2381 :     On entry, TRANS specifies the operation to be performed as
2382 :     follows:
2383 :     TRANS = 'N' or 'n' x := A*x.
2384 :     TRANS = 'T' or 't' x := A'*x.
2385 :     TRANS = 'C' or 'c' x := A'*x.
2386 :     Unchanged on exit.
2387 :     DIAG - CHARACTER*1.
2388 :     On entry, DIAG specifies whether or not A is unit
2389 :     triangular as follows:
2390 :     DIAG = 'U' or 'u' A is assumed to be unit triangular.
2391 :     DIAG = 'N' or 'n' A is not assumed to be unit
2392 :     triangular.
2393 :     Unchanged on exit.
2394 :     N - INTEGER.
2395 :     On entry, N specifies the order of the matrix A.
2396 :     N must be at least zero.
2397 :     Unchanged on exit.
2398 :     A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
2399 :     Before entry with UPLO = 'U' or 'u', the leading n by n
2400 :     upper triangular part of the array A must contain the upper
2401 :     triangular matrix and the strictly lower triangular part of
2402 :     A is not referenced.
2403 :     Before entry with UPLO = 'L' or 'l', the leading n by n
2404 :     lower triangular part of the array A must contain the lower
2405 :     triangular matrix and the strictly upper triangular part of
2406 :     A is not referenced.
2407 :     Note that when DIAG = 'U' or 'u', the diagonal elements of
2408 :     A are not referenced either, but are assumed to be unity.
2409 :     Unchanged on exit.
2410 :     LDA - INTEGER.
2411 :     On entry, LDA specifies the first dimension of A as declared
2412 :     in the calling (sub) program. LDA must be at least
2413 :     max( 1, n ).
2414 :     Unchanged on exit.
2415 :     X - DOUBLE PRECISION array of dimension at least
2416 :     ( 1 + ( n - 1 )*abs( INCX ) ).
2417 :     Before entry, the incremented array X must contain the n
2418 :     element vector x. On exit, X is overwritten with the
2419 :     tranformed vector x.
2420 :     INCX - INTEGER.
2421 :     On entry, INCX specifies the increment for the elements of
2422 :     X. INCX must not be zero.
2423 :     Unchanged on exit.
2424 :     Level 2 Blas routine.
2425 :     -- Written on 22-October-1986.
2426 :     Jack Dongarra, Argonne National Lab.
2427 :     Jeremy Du Croz, Nag Central Office.
2428 :     Sven Hammarling, Nag Central Office.
2429 :     Richard Hanson, Sandia National Labs.
2430 :     Test the input parameters.
2431 :     Parameter adjustments */
2432 :     a_dim1 = *lda;
2433 :     a_offset = 1 + a_dim1 * 1;
2434 :     a -= a_offset;
2435 :     --x;
2436 :     /* Function Body */
2437 :     info = 0;
2438 :     if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
2439 :     info = 1;
2440 :     } else if (! lsame_(trans, "N") && ! lsame_(trans,
2441 :     "T") && ! lsame_(trans, "C")) {
2442 :     info = 2;
2443 :     } else if (! lsame_(diag, "U") && ! lsame_(diag,
2444 :     "N")) {
2445 :     info = 3;
2446 :     } else if (*n < 0) {
2447 :     info = 4;
2448 :     } else if (*lda < max(1,*n)) {
2449 :     info = 6;
2450 :     } else if (*incx == 0) {
2451 :     info = 8;
2452 :     }
2453 :     if (info != 0) {
2454 :     xerbla_("DTRMV ", &info);
2455 :     return 0;
2456 :     }
2457 :     /* Quick return if possible. */
2458 :     if (*n == 0) {
2459 :     return 0;
2460 :     }
2461 :     nounit = lsame_(diag, "N");
2462 :     /* Set up the start point in X if the increment is not unity. This
2463 :     will be ( N - 1 )*INCX too small for descending loops. */
2464 :     if (*incx <= 0) {
2465 :     kx = 1 - (*n - 1) * *incx;
2466 :     } else if (*incx != 1) {
2467 :     kx = 1;
2468 :     }
2469 :     /* Start the operations. In this version the elements of A are
2470 :     accessed sequentially with one pass through A. */
2471 :     if (lsame_(trans, "N")) {
2472 :     /* Form x := A*x. */
2473 :     if (lsame_(uplo, "U")) {
2474 :     if (*incx == 1) {
2475 :     i__1 = *n;
2476 :     for (j = 1; j <= i__1; ++j) {
2477 :     if (x[j] != 0.) {
2478 :     temp = x[j];
2479 :     i__2 = j - 1;
2480 :     for (i__ = 1; i__ <= i__2; ++i__) {
2481 :     x[i__] += temp * a_ref(i__, j);
2482 :     /* L10: */
2483 :     }
2484 :     if (nounit) {
2485 :     x[j] *= a_ref(j, j);
2486 :     }
2487 :     }
2488 :     /* L20: */
2489 :     }
2490 :     } else {
2491 :     jx = kx;
2492 :     i__1 = *n;
2493 :     for (j = 1; j <= i__1; ++j) {
2494 :     if (x[jx] != 0.) {
2495 :     temp = x[jx];
2496 :     ix = kx;
2497 :     i__2 = j - 1;
2498 :     for (i__ = 1; i__ <= i__2; ++i__) {
2499 :     x[ix] += temp * a_ref(i__, j);
2500 :     ix += *incx;
2501 :     /* L30: */
2502 :     }
2503 :     if (nounit) {
2504 :     x[jx] *= a_ref(j, j);
2505 :     }
2506 :     }
2507 :     jx += *incx;
2508 :     /* L40: */
2509 :     }
2510 :     }
2511 :     } else {
2512 :     if (*incx == 1) {
2513 :     for (j = *n; j >= 1; --j) {
2514 :     if (x[j] != 0.) {
2515 :     temp = x[j];
2516 :     i__1 = j + 1;
2517 :     for (i__ = *n; i__ >= i__1; --i__) {
2518 :     x[i__] += temp * a_ref(i__, j);
2519 :     /* L50: */
2520 :     }
2521 :     if (nounit) {
2522 :     x[j] *= a_ref(j, j);
2523 :     }
2524 :     }
2525 :     /* L60: */
2526 :     }
2527 :     } else {
2528 :     kx += (*n - 1) * *incx;
2529 :     jx = kx;
2530 :     for (j = *n; j >= 1; --j) {
2531 :     if (x[jx] != 0.) {
2532 :     temp = x[jx];
2533 :     ix = kx;
2534 :     i__1 = j + 1;
2535 :     for (i__ = *n; i__ >= i__1; --i__) {
2536 :     x[ix] += temp * a_ref(i__, j);
2537 :     ix -= *incx;
2538 :     /* L70: */
2539 :     }
2540 :     if (nounit) {
2541 :     x[jx] *= a_ref(j, j);
2542 :     }
2543 :     }
2544 :     jx -= *incx;
2545 :     /* L80: */
2546 :     }
2547 :     }
2548 :     }
2549 :     } else {
2550 :     /* Form x := A'*x. */
2551 :     if (lsame_(uplo, "U")) {
2552 :     if (*incx == 1) {
2553 :     for (j = *n; j >= 1; --j) {
2554 :     temp = x[j];
2555 :     if (nounit) {
2556 :     temp *= a_ref(j, j);
2557 :     }
2558 :     for (i__ = j - 1; i__ >= 1; --i__) {
2559 :     temp += a_ref(i__, j) * x[i__];
2560 :     /* L90: */
2561 :     }
2562 :     x[j] = temp;
2563 :     /* L100: */
2564 :     }
2565 :     } else {
2566 :     jx = kx + (*n - 1) * *incx;
2567 :     for (j = *n; j >= 1; --j) {
2568 :     temp = x[jx];
2569 :     ix = jx;
2570 :     if (nounit) {
2571 :     temp *= a_ref(j, j);
2572 :     }
2573 :     for (i__ = j - 1; i__ >= 1; --i__) {
2574 :     ix -= *incx;
2575 :     temp += a_ref(i__, j) * x[ix];
2576 :     /* L110: */
2577 :     }
2578 :     x[jx] = temp;
2579 :     jx -= *incx;
2580 :     /* L120: */
2581 :     }
2582 :     }
2583 :     } else {
2584 :     if (*incx == 1) {
2585 :     i__1 = *n;
2586 :     for (j = 1; j <= i__1; ++j) {
2587 :     temp = x[j];
2588 :     if (nounit) {
2589 :     temp *= a_ref(j, j);
2590 :     }
2591 :     i__2 = *n;
2592 :     for (i__ = j + 1; i__ <= i__2; ++i__) {
2593 :     temp += a_ref(i__, j) * x[i__];
2594 :     /* L130: */
2595 :     }
2596 :     x[j] = temp;
2597 :     /* L140: */
2598 :     }
2599 :     } else {
2600 :     jx = kx;
2601 :     i__1 = *n;
2602 :     for (j = 1; j <= i__1; ++j) {
2603 :     temp = x[jx];
2604 :     ix = jx;
2605 :     if (nounit) {
2606 :     temp *= a_ref(j, j);
2607 :     }
2608 :     i__2 = *n;
2609 :     for (i__ = j + 1; i__ <= i__2; ++i__) {
2610 :     ix += *incx;
2611 :     temp += a_ref(i__, j) * x[ix];
2612 :     /* L150: */
2613 :     }
2614 :     x[jx] = temp;
2615 :     jx += *incx;
2616 :     /* L160: */
2617 :     }
2618 :     }
2619 :     }
2620 :     }
2621 :     return 0;
2622 :     /* End of DTRMV . */
2623 :     } /* dtrmv_ */
2624 :     #undef a_ref
2625 :    
2626 :    
2627 :     /* Subroutine */ int dtrsm_(char *side, char *uplo, char *transa, char *diag,
2628 :     integer *m, integer *n, doublereal *alpha, doublereal *a, integer *
2629 :     lda, doublereal *b, integer *ldb)
2630 :     {
2631 :     /* System generated locals */
2632 :     integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
2633 :     /* Local variables */
2634 :     static integer info;
2635 :     static doublereal temp;
2636 :     static integer i__, j, k;
2637 :     static logical lside;
2638 :     extern logical lsame_(char *, char *);
2639 :     static integer nrowa;
2640 :     static logical upper;
2641 :     extern /* Subroutine */ int xerbla_(char *, integer *);
2642 :     static logical nounit;
2643 :     #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
2644 :     #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
2645 :     /* Purpose
2646 :     =======
2647 :     DTRSM solves one of the matrix equations
2648 :     op( A )*X = alpha*B, or X*op( A ) = alpha*B,
2649 :     where alpha is a scalar, X and B are m by n matrices, A is a unit, or
2650 :     non-unit, upper or lower triangular matrix and op( A ) is one of
2651 :     op( A ) = A or op( A ) = A'.
2652 :     The matrix X is overwritten on B.
2653 :     Parameters
2654 :     ==========
2655 :     SIDE - CHARACTER*1.
2656 :     On entry, SIDE specifies whether op( A ) appears on the left
2657 :     or right of X as follows:
2658 :     SIDE = 'L' or 'l' op( A )*X = alpha*B.
2659 :     SIDE = 'R' or 'r' X*op( A ) = alpha*B.
2660 :     Unchanged on exit.
2661 :     UPLO - CHARACTER*1.
2662 :     On entry, UPLO specifies whether the matrix A is an upper or
2663 :     lower triangular matrix as follows:
2664 :     UPLO = 'U' or 'u' A is an upper triangular matrix.
2665 :     UPLO = 'L' or 'l' A is a lower triangular matrix.
2666 :     Unchanged on exit.
2667 :     TRANSA - CHARACTER*1.
2668 :     On entry, TRANSA specifies the form of op( A ) to be used in
2669 :     the matrix multiplication as follows:
2670 :     TRANSA = 'N' or 'n' op( A ) = A.
2671 :     TRANSA = 'T' or 't' op( A ) = A'.
2672 :     TRANSA = 'C' or 'c' op( A ) = A'.
2673 :     Unchanged on exit.
2674 :     DIAG - CHARACTER*1.
2675 :     On entry, DIAG specifies whether or not A is unit triangular
2676 :     as follows:
2677 :     DIAG = 'U' or 'u' A is assumed to be unit triangular.
2678 :     DIAG = 'N' or 'n' A is not assumed to be unit
2679 :     triangular.
2680 :     Unchanged on exit.
2681 :     M - INTEGER.
2682 :     On entry, M specifies the number of rows of B. M must be at
2683 :     least zero.
2684 :     Unchanged on exit.
2685 :     N - INTEGER.
2686 :     On entry, N specifies the number of columns of B. N must be
2687 :     at least zero.
2688 :     Unchanged on exit.
2689 :     ALPHA - DOUBLE PRECISION.
2690 :     On entry, ALPHA specifies the scalar alpha. When alpha is
2691 :     zero then A is not referenced and B need not be set before
2692 :     entry.
2693 :     Unchanged on exit.
2694 :     A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
2695 :     when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
2696 :     Before entry with UPLO = 'U' or 'u', the leading k by k
2697 :     upper triangular part of the array A must contain the upper
2698 :     triangular matrix and the strictly lower triangular part of
2699 :     A is not referenced.
2700 :     Before entry with UPLO = 'L' or 'l', the leading k by k
2701 :     lower triangular part of the array A must contain the lower
2702 :     triangular matrix and the strictly upper triangular part of
2703 :     A is not referenced.
2704 :     Note that when DIAG = 'U' or 'u', the diagonal elements of
2705 :     A are not referenced either, but are assumed to be unity.
2706 :     Unchanged on exit.
2707 :     LDA - INTEGER.
2708 :     On entry, LDA specifies the first dimension of A as declared
2709 :     in the calling (sub) program. When SIDE = 'L' or 'l' then
2710 :     LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
2711 :     then LDA must be at least max( 1, n ).
2712 :     Unchanged on exit.
2713 :     B - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
2714 :     Before entry, the leading m by n part of the array B must
2715 :     contain the right-hand side matrix B, and on exit is
2716 :     overwritten by the solution matrix X.
2717 :     LDB - INTEGER.
2718 :     On entry, LDB specifies the first dimension of B as declared
2719 :     in the calling (sub) program. LDB must be at least
2720 :     max( 1, m ).
2721 :     Unchanged on exit.
2722 :     Level 3 Blas routine.
2723 :     -- Written on 8-February-1989.
2724 :     Jack Dongarra, Argonne National Laboratory.
2725 :     Iain Duff, AERE Harwell.
2726 :     Jeremy Du Croz, Numerical Algorithms Group Ltd.
2727 :     Sven Hammarling, Numerical Algorithms Group Ltd.
2728 :     Test the input parameters.
2729 :     Parameter adjustments */
2730 :     a_dim1 = *lda;
2731 :     a_offset = 1 + a_dim1 * 1;
2732 :     a -= a_offset;
2733 :     b_dim1 = *ldb;
2734 :     b_offset = 1 + b_dim1 * 1;
2735 :     b -= b_offset;
2736 :     /* Function Body */
2737 :     lside = lsame_(side, "L");
2738 :     if (lside) {
2739 :     nrowa = *m;
2740 :     } else {
2741 :     nrowa = *n;
2742 :     }
2743 :     nounit = lsame_(diag, "N");
2744 :     upper = lsame_(uplo, "U");
2745 :     info = 0;
2746 :     if (! lside && ! lsame_(side, "R")) {
2747 :     info = 1;
2748 :     } else if (! upper && ! lsame_(uplo, "L")) {
2749 :     info = 2;
2750 :     } else if (! lsame_(transa, "N") && ! lsame_(transa,
2751 :     "T") && ! lsame_(transa, "C")) {
2752 :     info = 3;
2753 :     } else if (! lsame_(diag, "U") && ! lsame_(diag,
2754 :     "N")) {
2755 :     info = 4;
2756 :     } else if (*m < 0) {
2757 :     info = 5;
2758 :     } else if (*n < 0) {
2759 :     info = 6;
2760 :     } else if (*lda < max(1,nrowa)) {
2761 :     info = 9;
2762 :     } else if (*ldb < max(1,*m)) {
2763 :     info = 11;
2764 :     }
2765 :     if (info != 0) {
2766 :     xerbla_("DTRSM ", &info);
2767 :     return 0;
2768 :     }
2769 :     /* Quick return if possible. */
2770 :     if (*n == 0) {
2771 :     return 0;
2772 :     }
2773 :     /* And when alpha.eq.zero. */
2774 :     if (*alpha == 0.) {
2775 :     i__1 = *n;
2776 :     for (j = 1; j <= i__1; ++j) {
2777 :     i__2 = *m;
2778 :     for (i__ = 1; i__ <= i__2; ++i__) {
2779 :     b_ref(i__, j) = 0.;
2780 :     /* L10: */
2781 :     }
2782 :     /* L20: */
2783 :     }
2784 :     return 0;
2785 :     }
2786 :     /* Start the operations. */
2787 :     if (lside) {
2788 :     if (lsame_(transa, "N")) {
2789 :     /* Form B := alpha*inv( A )*B. */
2790 :     if (upper) {
2791 :     i__1 = *n;
2792 :     for (j = 1; j <= i__1; ++j) {
2793 :     if (*alpha != 1.) {
2794 :     i__2 = *m;
2795 :     for (i__ = 1; i__ <= i__2; ++i__) {
2796 :     b_ref(i__, j) = *alpha * b_ref(i__, j);
2797 :     /* L30: */
2798 :     }
2799 :     }
2800 :     for (k = *m; k >= 1; --k) {
2801 :     if (b_ref(k, j) != 0.) {
2802 :     if (nounit) {
2803 :     b_ref(k, j) = b_ref(k, j) / a_ref(k, k);
2804 :     }
2805 :     i__2 = k - 1;
2806 :     for (i__ = 1; i__ <= i__2; ++i__) {
2807 :     b_ref(i__, j) = b_ref(i__, j) - b_ref(k, j) *
2808 :     a_ref(i__, k);
2809 :     /* L40: */
2810 :     }
2811 :     }
2812 :     /* L50: */
2813 :     }
2814 :     /* L60: */
2815 :     }
2816 :     } else {
2817 :     i__1 = *n;
2818 :     for (j = 1; j <= i__1; ++j) {
2819 :     if (*alpha != 1.) {
2820 :     i__2 = *m;
2821 :     for (i__ = 1; i__ <= i__2; ++i__) {
2822 :     b_ref(i__, j) = *alpha * b_ref(i__, j);
2823 :     /* L70: */
2824 :     }
2825 :     }
2826 :     i__2 = *m;
2827 :     for (k = 1; k <= i__2; ++k) {
2828 :     if (b_ref(k, j) != 0.) {
2829 :     if (nounit) {
2830 :     b_ref(k, j) = b_ref(k, j) / a_ref(k, k);
2831 :     }
2832 :     i__3 = *m;
2833 :     for (i__ = k + 1; i__ <= i__3; ++i__) {
2834 :     b_ref(i__, j) = b_ref(i__, j) - b_ref(k, j) *
2835 :     a_ref(i__, k);
2836 :     /* L80: */
2837 :     }
2838 :     }
2839 :     /* L90: */
2840 :     }
2841 :     /* L100: */
2842 :     }
2843 :     }
2844 :     } else {
2845 :     /* Form B := alpha*inv( A' )*B. */
2846 :     if (upper) {
2847 :     i__1 = *n;
2848 :     for (j = 1; j <= i__1; ++j) {
2849 :     i__2 = *m;
2850 :     for (i__ = 1; i__ <= i__2; ++i__) {
2851 :     temp = *alpha * b_ref(i__, j);
2852 :     i__3 = i__ - 1;
2853 :     for (k = 1; k <= i__3; ++k) {
2854 :     temp -= a_ref(k, i__) * b_ref(k, j);
2855 :     /* L110: */
2856 :     }
2857 :     if (nounit) {
2858 :     temp /= a_ref(i__, i__);
2859 :     }
2860 :     b_ref(i__, j) = temp;
2861 :     /* L120: */
2862 :     }
2863 :     /* L130: */
2864 :     }
2865 :     } else {
2866 :     i__1 = *n;
2867 :     for (j = 1; j <= i__1; ++j) {
2868 :     for (i__ = *m; i__ >= 1; --i__) {
2869 :     temp = *alpha * b_ref(i__, j);
2870 :     i__2 = *m;
2871 :     for (k = i__ + 1; k <= i__2; ++k) {
2872 :     temp -= a_ref(k, i__) * b_ref(k, j);
2873 :     /* L140: */
2874 :     }
2875 :     if (nounit) {
2876 :     temp /= a_ref(i__, i__);
2877 :     }
2878 :     b_ref(i__, j) = temp;
2879 :     /* L150: */
2880 :     }
2881 :     /* L160: */
2882 :     }
2883 :     }
2884 :     }
2885 :     } else {
2886 :     if (lsame_(transa, "N")) {
2887 :     /* Form B := alpha*B*inv( A ). */
2888 :     if (upper) {
2889 :     i__1 = *n;
2890 :     for (j = 1; j <= i__1; ++j) {
2891 :     if (*alpha != 1.) {
2892 :     i__2 = *m;
2893 :     for (i__ = 1; i__ <= i__2; ++i__) {
2894 :     b_ref(i__, j) = *alpha * b_ref(i__, j);
2895 :     /* L170: */
2896 :     }
2897 :     }
2898 :     i__2 = j - 1;
2899 :     for (k = 1; k <= i__2; ++k) {
2900 :     if (a_ref(k, j) != 0.) {
2901 :     i__3 = *m;
2902 :     for (i__ = 1; i__ <= i__3; ++i__) {
2903 :     b_ref(i__, j) = b_ref(i__, j) - a_ref(k, j) *
2904 :     b_ref(i__, k);
2905 :     /* L180: */
2906 :     }
2907 :     }
2908 :     /* L190: */
2909 :     }
2910 :     if (nounit) {
2911 :     temp = 1. / a_ref(j, j);
2912 :     i__2 = *m;
2913 :     for (i__ = 1; i__ <= i__2; ++i__) {
2914 :     b_ref(i__, j) = temp * b_ref(i__, j);
2915 :     /* L200: */
2916 :     }
2917 :     }
2918 :     /* L210: */
2919 :     }
2920 :     } else {
2921 :     for (j = *n; j >= 1; --j) {
2922 :     if (*alpha != 1.) {
2923 :     i__1 = *m;
2924 :     for (i__ = 1; i__ <= i__1; ++i__) {
2925 :     b_ref(i__, j) = *alpha * b_ref(i__, j);
2926 :     /* L220: */
2927 :     }
2928 :     }
2929 :     i__1 = *n;
2930 :     for (k = j + 1; k <= i__1; ++k) {
2931 :     if (a_ref(k, j) != 0.) {
2932 :     i__2 = *m;
2933 :     for (i__ = 1; i__ <= i__2; ++i__) {
2934 :     b_ref(i__, j) = b_ref(i__, j) - a_ref(k, j) *
2935 :     b_ref(i__, k);
2936 :     /* L230: */
2937 :     }
2938 :     }
2939 :     /* L240: */
2940 :     }
2941 :     if (nounit) {
2942 :     temp = 1. / a_ref(j, j);
2943 :     i__1 = *m;
2944 :     for (i__ = 1; i__ <= i__1; ++i__) {
2945 :     b_ref(i__, j) = temp * b_ref(i__, j);
2946 :     /* L250: */
2947 :     }
2948 :     }
2949 :     /* L260: */
2950 :     }
2951 :     }
2952 :     } else {
2953 :     /* Form B := alpha*B*inv( A' ). */
2954 :     if (upper) {
2955 :     for (k = *n; k >= 1; --k) {
2956 :     if (nounit) {
2957 :     temp = 1. / a_ref(k, k);
2958 :     i__1 = *m;
2959 :     for (i__ = 1; i__ <= i__1; ++i__) {
2960 :     b_ref(i__, k) = temp * b_ref(i__, k);
2961 :     /* L270: */
2962 :     }
2963 :     }
2964 :     i__1 = k - 1;
2965 :     for (j = 1; j <= i__1; ++j) {
2966 :     if (a_ref(j, k) != 0.) {
2967 :     temp = a_ref(j, k);
2968 :     i__2 = *m;
2969 :     for (i__ = 1; i__ <= i__2; ++i__) {
2970 :     b_ref(i__, j) = b_ref(i__, j) - temp * b_ref(
2971 :     i__, k);
2972 :     /* L280: */
2973 :     }
2974 :     }
2975 :     /* L290: */
2976 :     }
2977 :     if (*alpha != 1.) {
2978 :     i__1 = *m;
2979 :     for (i__ = 1; i__ <= i__1; ++i__) {
2980 :     b_ref(i__, k) = *alpha * b_ref(i__, k);
2981 :     /* L300: */
2982 :     }
2983 :     }
2984 :     /* L310: */
2985 :     }
2986 :     } else {
2987 :     i__1 = *n;
2988 :     for (k = 1; k <= i__1; ++k) {
2989 :     if (nounit) {
2990 :     temp = 1. / a_ref(k, k);
2991 :     i__2 = *m;
2992 :     for (i__ = 1; i__ <= i__2; ++i__) {
2993 :     b_ref(i__, k) = temp * b_ref(i__, k);
2994 :     /* L320: */
2995 :     }
2996 :     }
2997 :     i__2 = *n;
2998 :     for (j = k + 1; j <= i__2; ++j) {
2999 :     if (a_ref(j, k) != 0.) {
3000 :     temp = a_ref(j, k);
3001 :     i__3 = *m;
3002 :     for (i__ = 1; i__ <= i__3; ++i__) {
3003 :     b_ref(i__, j) = b_ref(i__, j) - temp * b_ref(
3004 :     i__, k);
3005 :     /* L330: */
3006 :     }
3007 :     }
3008 :     /* L340: */
3009 :     }
3010 :     if (*alpha != 1.) {
3011 :     i__2 = *m;
3012 :     for (i__ = 1; i__ <= i__2; ++i__) {
3013 :     b_ref(i__, k) = *alpha * b_ref(i__, k);
3014 :     /* L350: */
3015 :     }
3016 :     }
3017 :     /* L360: */
3018 :     }
3019 :     }
3020 :     }
3021 :     }
3022 :     return 0;
3023 :     /* End of DTRSM . */
3024 :     } /* dtrsm_ */
3025 :     #undef b_ref
3026 :     #undef a_ref
3027 :    
3028 :    
3029 :     doublereal dzasum_(integer *n, doublecomplex *zx, integer *incx)
3030 :     {
3031 :     /* System generated locals */
3032 :     integer i__1;
3033 :     doublereal ret_val;
3034 :     /* Local variables */
3035 :     static integer i__;
3036 :     static doublereal stemp;
3037 :     extern doublereal dcabs1_(doublecomplex *);
3038 :     static integer ix;
3039 :     /* takes the sum of the absolute values.
3040 :     jack dongarra, 3/11/78.
3041 :     modified 3/93 to return if incx .le. 0.
3042 :     modified 12/3/93, array(1) declarations changed to array(*)
3043 :     Parameter adjustments */
3044 :     --zx;
3045 :     /* Function Body */
3046 :     ret_val = 0.;
3047 :     stemp = 0.;
3048 :     if (*n <= 0 || *incx <= 0) {
3049 :     return ret_val;
3050 :     }
3051 :     if (*incx == 1) {
3052 :     goto L20;
3053 :     }
3054 :     /* code for increment not equal to 1 */
3055 :     ix = 1;
3056 :     i__1 = *n;
3057 :     for (i__ = 1; i__ <= i__1; ++i__) {
3058 :     stemp += dcabs1_(&zx[ix]);
3059 :     ix += *incx;
3060 :     /* L10: */
3061 :     }
3062 :     ret_val = stemp;
3063 :     return ret_val;
3064 :     /* code for increment equal to 1 */
3065 :     L20:
3066 :     i__1 = *n;
3067 :     for (i__ = 1; i__ <= i__1; ++i__) {
3068 :     stemp += dcabs1_(&zx[i__]);
3069 :     /* L30: */
3070 :     }
3071 :     ret_val = stemp;
3072 :     return ret_val;
3073 :     } /* dzasum_ */
3074 :    
3075 :    
3076 :     doublereal dznrm2_(integer *n, doublecomplex *x, integer *incx)
3077 :     {
3078 :     /* The following loop is equivalent to this call to the LAPACK
3079 :     auxiliary routine:
3080 :     CALL ZLASSQ( N, X, INCX, SCALE, SSQ ) */
3081 :     /* System generated locals */
3082 :     integer i__1, i__2, i__3;
3083 :     doublereal ret_val, d__1;
3084 :     /* Builtin functions */
3085 :     double d_imag(doublecomplex *), sqrt(doublereal);
3086 :     /* Local variables */
3087 :     static doublereal temp, norm, scale;
3088 :     static integer ix;
3089 :     static doublereal ssq;
3090 :     /* DZNRM2 returns the euclidean norm of a vector via the function
3091 :     name, so that
3092 :     DZNRM2 := sqrt( conjg( x' )*x )
3093 :     -- This version written on 25-October-1982.
3094 :     Modified on 14-October-1993 to inline the call to ZLASSQ.
3095 :     Sven Hammarling, Nag Ltd.
3096 :     Parameter adjustments */
3097 :     --x;
3098 :     /* Function Body */
3099 :     if (*n < 1 || *incx < 1) {
3100 :     norm = 0.;
3101 :     } else {
3102 :     scale = 0.;
3103 :     ssq = 1.;
3104 :    
3105 :    
3106 :     i__1 = (*n - 1) * *incx + 1;
3107 :     i__2 = *incx;
3108 :     for (ix = 1; i__2 < 0 ? ix >= i__1 : ix <= i__1; ix += i__2) {
3109 :     i__3 = ix;
3110 :     if (x[i__3].r != 0.) {
3111 :     i__3 = ix;
3112 :     temp = (d__1 = x[i__3].r, abs(d__1));
3113 :     if (scale < temp) {
3114 :     /* Computing 2nd power */
3115 :     d__1 = scale / temp;
3116 :     ssq = ssq * (d__1 * d__1) + 1.;
3117 :     scale = temp;
3118 :     } else {
3119 :     /* Computing 2nd power */
3120 :     d__1 = temp / scale;
3121 :     ssq += d__1 * d__1;
3122 :     }
3123 :     }
3124 :     if (d_imag(&x[ix]) != 0.) {
3125 :     temp = (d__1 = d_imag(&x[ix]), abs(d__1));
3126 :     if (scale < temp) {
3127 :     /* Computing 2nd power */
3128 :     d__1 = scale / temp;
3129 :     ssq = ssq * (d__1 * d__1) + 1.;
3130 :     scale = temp;
3131 :     } else {
3132 :     /* Computing 2nd power */
3133 :     d__1 = temp / scale;
3134 :     ssq += d__1 * d__1;
3135 :     }
3136 :     }
3137 :     /* L10: */
3138 :     }
3139 :     norm = scale * sqrt(ssq);
3140 :     }
3141 :    
3142 :     ret_val = norm;
3143 :     return ret_val;
3144 :    
3145 :     /* End of DZNRM2. */
3146 :    
3147 :     } /* dznrm2_ */
3148 :    
3149 :    
3150 :     integer idamax_(integer *n, doublereal *dx, integer *incx)
3151 :     {
3152 :     /* System generated locals */
3153 :     integer ret_val, i__1;
3154 :     doublereal d__1;
3155 :     /* Local variables */
3156 :     static doublereal dmax__;
3157 :     static integer i__, ix;
3158 :     /* finds the index of element having max. absolute value.
3159 :     jack dongarra, linpack, 3/11/78.
3160 :     modified 3/93 to return if incx .le. 0.
3161 :     modified 12/3/93, array(1) declarations changed to array(*)
3162 :     Parameter adjustments */
3163 :     --dx;
3164 :     /* Function Body */
3165 :     ret_val = 0;
3166 :     if (*n < 1 || *incx <= 0) {
3167 :     return ret_val;
3168 :     }
3169 :     ret_val = 1;
3170 :     if (*n == 1) {
3171 :     return ret_val;
3172 :     }
3173 :     if (*incx == 1) {
3174 :     goto L20;
3175 :     }
3176 :     /* code for increment not equal to 1 */
3177 :     ix = 1;
3178 :     dmax__ = abs(dx[1]);
3179 :     ix += *incx;
3180 :     i__1 = *n;
3181 :     for (i__ = 2; i__ <= i__1; ++i__) {
3182 :     if ((d__1 = dx[ix], abs(d__1)) <= dmax__) {
3183 :     goto L5;
3184 :     }
3185 :     ret_val = i__;
3186 :     dmax__ = (d__1 = dx[ix], abs(d__1));
3187 :     L5:
3188 :     ix += *incx;
3189 :     /* L10: */
3190 :     }
3191 :     return ret_val;
3192 :     /* code for increment equal to 1 */
3193 :     L20:
3194 :     dmax__ = abs(dx[1]);
3195 :     i__1 = *n;
3196 :     for (i__ = 2; i__ <= i__1; ++i__) {
3197 :     if ((d__1 = dx[i__], abs(d__1)) <= dmax__) {
3198 :     goto L30;
3199 :     }
3200 :     ret_val = i__;
3201 :     dmax__ = (d__1 = dx[i__], abs(d__1));
3202 :     L30:
3203 :     ;
3204 :     }
3205 :     return ret_val;
3206 :     } /* idamax_ */
3207 :    
3208 :    
3209 :     integer izamax_(integer *n, doublecomplex *zx, integer *incx)
3210 :     {
3211 :     /* System generated locals */
3212 :     integer ret_val, i__1;
3213 :     /* Local variables */
3214 :     static doublereal smax;
3215 :     static integer i__;
3216 :     extern doublereal dcabs1_(doublecomplex *);
3217 :     static integer ix;
3218 :     /* finds the index of element having max. absolute value.
3219 :     jack dongarra, 1/15/85.
3220 :     modified 3/93 to return if incx .le. 0.
3221 :     modified 12/3/93, array(1) declarations changed to array(*)
3222 :     Parameter adjustments */
3223 :     --zx;
3224 :     /* Function Body */
3225 :     ret_val = 0;
3226 :     if (*n < 1 || *incx <= 0) {
3227 :     return ret_val;
3228 :     }
3229 :     ret_val = 1;
3230 :     if (*n == 1) {
3231 :     return ret_val;
3232 :     }
3233 :     if (*incx == 1) {
3234 :     goto L20;
3235 :     }
3236 :     /* code for increment not equal to 1 */
3237 :     ix = 1;
3238 :     smax = dcabs1_(&zx[1]);
3239 :     ix += *incx;
3240 :     i__1 = *n;
3241 :     for (i__ = 2; i__ <= i__1; ++i__) {
3242 :     if (dcabs1_(&zx[ix]) <= smax) {
3243 :     goto L5;
3244 :     }
3245 :     ret_val = i__;
3246 :     smax = dcabs1_(&zx[ix]);
3247 :     L5:
3248 :     ix += *incx;
3249 :     /* L10: */
3250 :     }
3251 :     return ret_val;
3252 :     /* code for increment equal to 1 */
3253 :     L20:
3254 :     smax = dcabs1_(&zx[1]);
3255 :     i__1 = *n;
3256 :     for (i__ = 2; i__ <= i__1; ++i__) {
3257 :     if (dcabs1_(&zx[i__]) <= smax) {
3258 :     goto L30;
3259 :     }
3260 :     ret_val = i__;
3261 :     smax = dcabs1_(&zx[i__]);
3262 :     L30:
3263 :     ;
3264 :     }
3265 :     return ret_val;
3266 :     } /* izamax_ */
3267 :    
3268 :    
3269 :     /* Subroutine */ int xerbla_(char *srname, integer *info)
3270 :     {
3271 :     /* -- LAPACK auxiliary routine (preliminary version) --
3272 :     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
3273 :     Courant Institute, Argonne National Lab, and Rice University
3274 :     February 29, 1992
3275 :    
3276 :    
3277 :     Purpose
3278 :     =======
3279 :    
3280 :     XERBLA is an error handler for the LAPACK routines.
3281 :     It is called by an LAPACK routine if an input parameter has an
3282 :     invalid value. A message is printed and execution stops.
3283 :    
3284 :     Installers may consider modifying the STOP statement in order to
3285 :     call system-specific exception-handling facilities.
3286 :    
3287 :     Arguments
3288 :     =========
3289 :    
3290 :     SRNAME (input) CHARACTER*6
3291 :     The name of the routine which called XERBLA.
3292 :    
3293 :     INFO (input) INTEGER
3294 :     The position of the invalid parameter in the parameter list
3295 :     of the calling routine. */
3296 :     /* Table of constant values */
3297 :     static integer c__1 = 1;
3298 :    
3299 :     /* Format strings */
3300 :     static char fmt_9999[] = "(\002 ** On entry to \002,a6,\002 parameter nu"
3301 :     "mber \002,i2,\002 had \002,\002an illegal value\002)";
3302 :     /* Builtin functions */
3303 :     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
3304 :     /* Subroutine */ int s_stop(char *, ftnlen);
3305 :     /* Fortran I/O blocks */
3306 :     static cilist io___1 = { 0, 6, 0, fmt_9999, 0 };
3307 :    
3308 :    
3309 :    
3310 :    
3311 :     s_wsfe(&io___1);
3312 :     do_fio(&c__1, srname, (ftnlen)6);
3313 :     do_fio(&c__1, (char *)&(*info), (ftnlen)sizeof(integer));
3314 :     e_wsfe();
3315 :    
3316 :     s_stop("", (ftnlen)0);
3317 :    
3318 :    
3319 :     /* End of XERBLA */
3320 :    
3321 :     return 0;
3322 :     } /* xerbla_ */
3323 :    
3324 :    
3325 :     /* Subroutine */ int zaxpy_(integer *n, doublecomplex *za, doublecomplex *zx,
3326 :     integer *incx, doublecomplex *zy, integer *incy)
3327 :     {
3328 :     /* System generated locals */
3329 :     integer i__1, i__2, i__3, i__4;
3330 :     doublecomplex z__1, z__2;
3331 :     /* Local variables */
3332 :     static integer i__;
3333 :     extern doublereal dcabs1_(doublecomplex *);
3334 :     static integer ix, iy;
3335 :     /* constant times a vector plus a vector.
3336 :     jack dongarra, 3/11/78.
3337 :     modified 12/3/93, array(1) declarations changed to array(*)
3338 :     Parameter adjustments */
3339 :     --zy;
3340 :     --zx;
3341 :     /* Function Body */
3342 :     if (*n <= 0) {
3343 :     return 0;
3344 :     }
3345 :     if (dcabs1_(za) == 0.) {
3346 :     return 0;
3347 :     }
3348 :     if (*incx == 1 && *incy == 1) {
3349 :     goto L20;
3350 :     }
3351 :     /* code for unequal increments or equal increments
3352 :     not equal to 1 */
3353 :     ix = 1;
3354 :     iy = 1;
3355 :     if (*incx < 0) {
3356 :     ix = (-(*n) + 1) * *incx + 1;
3357 :     }
3358 :     if (*incy < 0) {
3359 :     iy = (-(*n) + 1) * *incy + 1;
3360 :     }
3361 :     i__1 = *n;
3362 :     for (i__ = 1; i__ <= i__1; ++i__) {
3363 :     i__2 = iy;
3364 :     i__3 = iy;
3365 :     i__4 = ix;
3366 :     z__2.r = za->r * zx[i__4].r - za->i * zx[i__4].i, z__2.i = za->r * zx[
3367 :     i__4].i + za->i * zx[i__4].r;
3368 :     z__1.r = zy[i__3].r + z__2.r, z__1.i = zy[i__3].i + z__2.i;
3369 :     zy[i__2].r = z__1.r, zy[i__2].i = z__1.i;
3370 :     ix += *incx;
3371 :     iy += *incy;
3372 :     /* L10: */
3373 :     }
3374 :     return 0;
3375 :     /* code for both increments equal to 1 */
3376 :     L20:
3377 :     i__1 = *n;
3378 :     for (i__ = 1; i__ <= i__1; ++i__) {
3379 :     i__2 = i__;
3380 :     i__3 = i__;
3381 :     i__4 = i__;
3382 :     z__2.r = za->r * zx[i__4].r - za->i * zx[i__4].i, z__2.i = za->r * zx[
3383 :     i__4].i + za->i * zx[i__4].r;
3384 :     z__1.r = zy[i__3].r + z__2.r, z__1.i = zy[i__3].i + z__2.i;
3385 :     zy[i__2].r = z__1.r, zy[i__2].i = z__1.i;
3386 :     /* L30: */
3387 :     }
3388 :     return 0;
3389 :     } /* zaxpy_ */
3390 :    
3391 :    
3392 :     /* Subroutine */ int zcopy_(integer *n, doublecomplex *zx, integer *incx,
3393 :     doublecomplex *zy, integer *incy)
3394 :     {
3395 :     /* System generated locals */
3396 :     integer i__1, i__2, i__3;
3397 :     /* Local variables */
3398 :     static integer i__, ix, iy;
3399 :     /* copies a vector, x, to a vector, y.
3400 :     jack dongarra, linpack, 4/11/78.
3401 :     modified 12/3/93, array(1) declarations changed to array(*)
3402 :     Parameter adjustments */
3403 :     --zy;
3404 :     --zx;
3405 :     /* Function Body */
3406 :     if (*n <= 0) {
3407 :     return 0;
3408 :     }
3409 :     if (*incx == 1 && *incy == 1) {
3410 :     goto L20;
3411 :     }
3412 :     /* code for unequal increments or equal increments
3413 :     not equal to 1 */
3414 :     ix = 1;
3415 :     iy = 1;
3416 :     if (*incx < 0) {
3417 :     ix = (-(*n) + 1) * *incx + 1;
3418 :     }
3419 :     if (*incy < 0) {
3420 :     iy = (-(*n) + 1) * *incy + 1;
3421 :     }
3422 :     i__1 = *n;
3423 :     for (i__ = 1; i__ <= i__1; ++i__) {
3424 :     i__2 = iy;
3425 :     i__3 = ix;
3426 :     zy[i__2].r = zx[i__3].r, zy[i__2].i = zx[i__3].i;
3427 :     ix += *incx;
3428 :     iy += *incy;
3429 :     /* L10: */
3430 :     }
3431 :     return 0;
3432 :     /* code for both increments equal to 1 */
3433 :     L20:
3434 :     i__1 = *n;
3435 :     for (i__ = 1; i__ <= i__1; ++i__) {
3436 :     i__2 = i__;
3437 :     i__3 = i__;
3438 :     zy[i__2].r = zx[i__3].r, zy[i__2].i = zx[i__3].i;
3439 :     /* L30: */
3440 :     }
3441 :     return 0;
3442 :     } /* zcopy_ */
3443 :    
3444 :    
3445 :     /* Double Complex */ VOID zdotc_(doublecomplex * ret_val, integer *n,
3446 :     doublecomplex *zx, integer *incx, doublecomplex *zy, integer *incy)
3447 :     {
3448 :     /* System generated locals */
3449 :     integer i__1, i__2;
3450 :     doublecomplex z__1, z__2, z__3;
3451 :     /* Builtin functions */
3452 :     void d_cnjg(doublecomplex *, doublecomplex *);
3453 :     /* Local variables */
3454 :     static integer i__;
3455 :     static doublecomplex ztemp;
3456 :     static integer ix, iy;
3457 :     /* forms the dot product of a vector.
3458 :     jack dongarra, 3/11/78.
3459 :     modified 12/3/93, array(1) declarations changed to array(*)
3460 :     Parameter adjustments */
3461 :     --zy;
3462 :     --zx;
3463 :     /* Function Body */
3464 :     ztemp.r = 0., ztemp.i = 0.;
3465 :     ret_val->r = 0., ret_val->i = 0.;
3466 :     if (*n <= 0) {
3467 :     return ;
3468 :     }
3469 :     if (*incx == 1 && *incy == 1) {
3470 :     goto L20;
3471 :     }
3472 :     /* code for unequal increments or equal increments
3473 :     not equal to 1 */
3474 :     ix = 1;
3475 :     iy = 1;
3476 :     if (*incx < 0) {
3477 :     ix = (-(*n) + 1) * *incx + 1;
3478 :     }
3479 :     if (*incy < 0) {
3480 :     iy = (-(*n) + 1) * *incy + 1;
3481 :     }
3482 :     i__1 = *n;
3483 :     for (i__ = 1; i__ <= i__1; ++i__) {
3484 :     d_cnjg(&z__3, &zx[ix]);
3485 :     i__2 = iy;
3486 :     z__2.r = z__3.r * zy[i__2].r - z__3.i * zy[i__2].i, z__2.i = z__3.r *
3487 :     zy[i__2].i + z__3.i * zy[i__2].r;
3488 :     z__1.r = ztemp.r + z__2.r, z__1.i = ztemp.i + z__2.i;
3489 :     ztemp.r = z__1.r, ztemp.i = z__1.i;
3490 :     ix += *incx;
3491 :     iy += *incy;
3492 :     /* L10: */
3493 :     }
3494 :     ret_val->r = ztemp.r, ret_val->i = ztemp.i;
3495 :     return ;
3496 :     /* code for both increments equal to 1 */
3497 :     L20:
3498 :     i__1 = *n;
3499 :     for (i__ = 1; i__ <= i__1; ++i__) {
3500 :     d_cnjg(&z__3, &zx[i__]);
3501 :     i__2 = i__;
3502 :     z__2.r = z__3.r * zy[i__2].r - z__3.i * zy[i__2].i, z__2.i = z__3.r *
3503 :     zy[i__2].i + z__3.i * zy[i__2].r;
3504 :     z__1.r = ztemp.r + z__2.r, z__1.i = ztemp.i + z__2.i;
3505 :     ztemp.r = z__1.r, ztemp.i = z__1.i;
3506 :     /* L30: */
3507 :     }
3508 :     ret_val->r = ztemp.r, ret_val->i = ztemp.i;
3509 :     return ;
3510 :     } /* zdotc_ */
3511 :    
3512 :    
3513 :     /* Double Complex */ VOID zdotu_(doublecomplex * ret_val, integer *n,
3514 :     doublecomplex *zx, integer *incx, doublecomplex *zy, integer *incy)
3515 :     {
3516 :     /* System generated locals */
3517 :     integer i__1, i__2, i__3;
3518 :     doublecomplex z__1, z__2;
3519 :     /* Local variables */
3520 :     static integer i__;
3521 :     static doublecomplex ztemp;
3522 :     static integer ix, iy;
3523 :     /* forms the dot product of two vectors.
3524 :     jack dongarra, 3/11/78.
3525 :     modified 12/3/93, array(1) declarations changed to array(*)
3526 :     Parameter adjustments */
3527 :     --zy;
3528 :     --zx;
3529 :     /* Function Body */
3530 :     ztemp.r = 0., ztemp.i = 0.;
3531 :     ret_val->r = 0., ret_val->i = 0.;
3532 :     if (*n <= 0) {
3533 :     return ;
3534 :     }
3535 :     if (*incx == 1 && *incy == 1) {
3536 :     goto L20;
3537 :     }
3538 :     /* code for unequal increments or equal increments
3539 :     not equal to 1 */
3540 :     ix = 1;
3541 :     iy = 1;
3542 :     if (*incx < 0) {
3543 :     ix = (-(*n) + 1) * *incx + 1;
3544 :     }
3545 :     if (*incy < 0) {
3546 :     iy = (-(*n) + 1) * *incy + 1;
3547 :     }
3548 :     i__1 = *n;
3549 :     for (i__ = 1; i__ <= i__1; ++i__) {
3550 :     i__2 = ix;
3551 :     i__3 = iy;
3552 :     z__2.r = zx[i__2].r * zy[i__3].r - zx[i__2].i * zy[i__3].i, z__2.i =
3553 :     zx[i__2].r * zy[i__3].i + zx[i__2].i * zy[i__3].r;
3554 :     z__1.r = ztemp.r + z__2.r, z__1.i = ztemp.i + z__2.i;
3555 :     ztemp.r = z__1.r, ztemp.i = z__1.i;
3556 :     ix += *incx;
3557 :     iy += *incy;
3558 :     /* L10: */
3559 :     }
3560 :     ret_val->r = ztemp.r, ret_val->i = ztemp.i;
3561 :     return ;
3562 :     /* code for both increments equal to 1 */
3563 :     L20:
3564 :     i__1 = *n;
3565 :     for (i__ = 1; i__ <= i__1; ++i__) {
3566 :     i__2 = i__;
3567 :     i__3 = i__;
3568 :     z__2.r = zx[i__2].r * zy[i__3].r - zx[i__2].i * zy[i__3].i, z__2.i =
3569 :     zx[i__2].r * zy[i__3].i + zx[i__2].i * zy[i__3].r;
3570 :     z__1.r = ztemp.r + z__2.r, z__1.i = ztemp.i + z__2.i;
3571 :     ztemp.r = z__1.r, ztemp.i = z__1.i;
3572 :     /* L30: */
3573 :     }
3574 :     ret_val->r = ztemp.r, ret_val->i = ztemp.i;
3575 :     return ;
3576 :     } /* zdotu_ */
3577 :    
3578 :    
3579 :     /* Subroutine */ int zdscal_(integer *n, doublereal *da, doublecomplex *zx,
3580 :     integer *incx)
3581 :     {
3582 :     /* System generated locals */
3583 :     integer i__1, i__2, i__3;
3584 :     doublecomplex z__1, z__2;
3585 :     /* Local variables */
3586 :     static integer i__, ix;
3587 :     /* scales a vector by a constant.
3588 :     jack dongarra, 3/11/78.
3589 :     modified 3/93 to return if incx .le. 0.
3590 :     modified 12/3/93, array(1) declarations changed to array(*)
3591 :     Parameter adjustments */
3592 :     --zx;
3593 :     /* Function Body */
3594 :     if (*n <= 0 || *incx <= 0) {
3595 :     return 0;
3596 :     }
3597 :     if (*incx == 1) {
3598 :     goto L20;
3599 :     }
3600 :     /* code for increment not equal to 1 */
3601 :     ix = 1;
3602 :     i__1 = *n;
3603 :     for (i__ = 1; i__ <= i__1; ++i__) {
3604 :     i__2 = ix;
3605 :     z__2.r = *da, z__2.i = 0.;
3606 :     i__3 = ix;
3607 :     z__1.r = z__2.r * zx[i__3].r - z__2.i * zx[i__3].i, z__1.i = z__2.r *
3608 :     zx[i__3].i + z__2.i * zx[i__3].r;
3609 :     zx[i__2].r = z__1.r, zx[i__2].i = z__1.i;
3610 :     ix += *incx;
3611 :     /* L10: */
3612 :     }
3613 :     return 0;
3614 :     /* code for increment equal to 1 */
3615 :     L20:
3616 :     i__1 = *n;
3617 :     for (i__ = 1; i__ <= i__1; ++i__) {
3618 :     i__2 = i__;
3619 :     z__2.r = *da, z__2.i = 0.;
3620 :     i__3 = i__;
3621 :     z__1.r = z__2.r * zx[i__3].r - z__2.i * zx[i__3].i, z__1.i = z__2.r *
3622 :     zx[i__3].i + z__2.i * zx[i__3].r;
3623 :     zx[i__2].r = z__1.r, zx[i__2].i = z__1.i;
3624 :     /* L30: */
3625 :     }
3626 :     return 0;
3627 :     } /* zdscal_ */
3628 :    
3629 :    
3630 :     /* Subroutine */ int zgemm_(char *transa, char *transb, integer *m, integer *
3631 :     n, integer *k, doublecomplex *alpha, doublecomplex *a, integer *lda,
3632 :     doublecomplex *b, integer *ldb, doublecomplex *beta, doublecomplex *
3633 :     c__, integer *ldc)
3634 :     {
3635 :     /* System generated locals */
3636 :     integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
3637 :     i__3, i__4, i__5, i__6;
3638 :     doublecomplex z__1, z__2, z__3, z__4;
3639 :     /* Builtin functions */
3640 :     void d_cnjg(doublecomplex *, doublecomplex *);
3641 :     /* Local variables */
3642 :     static integer info;
3643 :     static logical nota, notb;
3644 :     static doublecomplex temp;
3645 :     static integer i__, j, l;
3646 :     static logical conja, conjb;
3647 :     static integer ncola;
3648 :     extern logical lsame_(char *, char *);
3649 :     static integer nrowa, nrowb;
3650 :     extern /* Subroutine */ int xerbla_(char *, integer *);
3651 :     #define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1
3652 :     #define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]
3653 :     #define b_subscr(a_1,a_2) (a_2)*b_dim1 + a_1
3654 :     #define b_ref(a_1,a_2) b[b_subscr(a_1,a_2)]
3655 :     #define c___subscr(a_1,a_2) (a_2)*c_dim1 + a_1
3656 :     #define c___ref(a_1,a_2) c__[c___subscr(a_1,a_2)]
3657 :     /* Purpose
3658 :     =======
3659 :     ZGEMM performs one of the matrix-matrix operations
3660 :     C := alpha*op( A )*op( B ) + beta*C,
3661 :     where op( X ) is one of
3662 :     op( X ) = X or op( X ) = X' or op( X ) = conjg( X' ),
3663 :     alpha and beta are scalars, and A, B and C are matrices, with op( A )
3664 :     an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
3665 :     Parameters
3666 :     ==========
3667 :     TRANSA - CHARACTER*1.
3668 :     On entry, TRANSA specifies the form of op( A ) to be used in
3669 :     the matrix multiplication as follows:
3670 :     TRANSA = 'N' or 'n', op( A ) = A.
3671 :     TRANSA = 'T' or 't', op( A ) = A'.
3672 :     TRANSA = 'C' or 'c', op( A ) = conjg( A' ).
3673 :     Unchanged on exit.
3674 :     TRANSB - CHARACTER*1.
3675 :     On entry, TRANSB specifies the form of op( B ) to be used in
3676 :     the matrix multiplication as follows:
3677 :     TRANSB = 'N' or 'n', op( B ) = B.
3678 :     TRANSB = 'T' or 't', op( B ) = B'.
3679 :     TRANSB = 'C' or 'c', op( B ) = conjg( B' ).
3680 :     Unchanged on exit.
3681 :     M - INTEGER.
3682 :     On entry, M specifies the number of rows of the matrix
3683 :     op( A ) and of the matrix C. M must be at least zero.
3684 :     Unchanged on exit.
3685 :     N - INTEGER.
3686 :     On entry, N specifies the number of columns of the matrix
3687 :     op( B ) and the number of columns of the matrix C. N must be
3688 :     at least zero.
3689 :     Unchanged on exit.
3690 :     K - INTEGER.
3691 :     On entry, K specifies the number of columns of the matrix
3692 :     op( A ) and the number of rows of the matrix op( B ). K must
3693 :     be at least zero.
3694 :     Unchanged on exit.
3695 :     ALPHA - COMPLEX*16 .
3696 :     On entry, ALPHA specifies the scalar alpha.
3697 :     Unchanged on exit.
3698 :     A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is
3699 :     k when TRANSA = 'N' or 'n', and is m otherwise.
3700 :     Before entry with TRANSA = 'N' or 'n', the leading m by k
3701 :     part of the array A must contain the matrix A, otherwise
3702 :     the leading k by m part of the array A must contain the
3703 :     matrix A.
3704 :     Unchanged on exit.
3705 :     LDA - INTEGER.
3706 :     On entry, LDA specifies the first dimension of A as declared
3707 :     in the calling (sub) program. When TRANSA = 'N' or 'n' then
3708 :     LDA must be at least max( 1, m ), otherwise LDA must be at
3709 :     least max( 1, k ).
3710 :     Unchanged on exit.
3711 :     B - COMPLEX*16 array of DIMENSION ( LDB, kb ), where kb is
3712 :     n when TRANSB = 'N' or 'n', and is k otherwise.
3713 :     Before entry with TRANSB = 'N' or 'n', the leading k by n
3714 :     part of the array B must contain the matrix B, otherwise
3715 :     the leading n by k part of the array B must contain the
3716 :     matrix B.
3717 :     Unchanged on exit.
3718 :     LDB - INTEGER.
3719 :     On entry, LDB specifies the first dimension of B as declared
3720 :     in the calling (sub) program. When TRANSB = 'N' or 'n' then
3721 :     LDB must be at least max( 1, k ), otherwise LDB must be at
3722 :     least max( 1, n ).
3723 :     Unchanged on exit.
3724 :     BETA - COMPLEX*16 .
3725 :     On entry, BETA specifies the scalar beta. When BETA is
3726 :     supplied as zero then C need not be set on input.
3727 :     Unchanged on exit.
3728 :     C - COMPLEX*16 array of DIMENSION ( LDC, n ).
3729 :     Before entry, the leading m by n part of the array C must
3730 :     contain the matrix C, except when beta is zero, in which
3731 :     case C need not be set on entry.
3732 :     On exit, the array C is overwritten by the m by n matrix
3733 :     ( alpha*op( A )*op( B ) + beta*C ).
3734 :     LDC - INTEGER.
3735 :     On entry, LDC specifies the first dimension of C as declared
3736 :     in the calling (sub) program. LDC must be at least
3737 :     max( 1, m ).
3738 :     Unchanged on exit.
3739 :     Level 3 Blas routine.
3740 :     -- Written on 8-February-1989.
3741 :     Jack Dongarra, Argonne National Laboratory.
3742 :     Iain Duff, AERE Harwell.
3743 :     Jeremy Du Croz, Numerical Algorithms Group Ltd.
3744 :     Sven Hammarling, Numerical Algorithms Group Ltd.
3745 :     Set NOTA and NOTB as true if A and B respectively are not
3746 :     conjugated or transposed, set CONJA and CONJB as true if A and
3747 :     B respectively are to be transposed but not conjugated and set
3748 :     NROWA, NCOLA and NROWB as the number of rows and columns of A
3749 :     and the number of rows of B respectively.
3750 :     Parameter adjustments */
3751 :     a_dim1 = *lda;
3752 :     a_offset = 1 + a_dim1 * 1;
3753 :     a -= a_offset;
3754 :     b_dim1 = *ldb;
3755 :     b_offset = 1 + b_dim1 * 1;
3756 :     b -= b_offset;
3757 :     c_dim1 = *ldc;
3758 :     c_offset = 1 + c_dim1 * 1;
3759 :     c__ -= c_offset;
3760 :     /* Function Body */
3761 :     nota = lsame_(transa, "N");
3762 :     notb = lsame_(transb, "N");
3763 :     conja = lsame_(transa, "C");
3764 :     conjb = lsame_(transb, "C");
3765 :     if (nota) {
3766 :     nrowa = *m;
3767 :     ncola = *k;
3768 :     } else {
3769 :     nrowa = *k;
3770 :     ncola = *m;
3771 :     }
3772 :     if (notb) {
3773 :     nrowb = *k;
3774 :     } else {
3775 :     nrowb = *n;
3776 :     }
3777 :     /* Test the input parameters. */
3778 :     info = 0;
3779 :     if (! nota && ! conja && ! lsame_(transa, "T")) {
3780 :     info = 1;
3781 :     } else if (! notb && ! conjb && ! lsame_(transb, "T")) {
3782 :     info = 2;
3783 :     } else if (*m < 0) {
3784 :     info = 3;
3785 :     } else if (*n < 0) {
3786 :     info = 4;
3787 :     } else if (*k < 0) {
3788 :     info = 5;
3789 :     } else if (*lda < max(1,nrowa)) {
3790 :     info = 8;
3791 :     } else if (*ldb < max(1,nrowb)) {
3792 :     info = 10;
3793 :     } else if (*ldc < max(1,*m)) {
3794 :     info = 13;
3795 :     }
3796 :     if (info != 0) {
3797 :     xerbla_("ZGEMM ", &info);
3798 :     return 0;
3799 :     }
3800 :     /* Quick return if possible. */
3801 :     if (*m == 0 || *n == 0 || (alpha->r == 0. && alpha->i == 0. || *k == 0) &&
3802 :     (beta->r == 1. && beta->i == 0.)) {
3803 :     return 0;
3804 :     }
3805 :     /* And when alpha.eq.zero. */
3806 :     if (alpha->r == 0. && alpha->i == 0.) {
3807 :     if (beta->r == 0. && beta->i == 0.) {
3808 :     i__1 = *n;
3809 :     for (j = 1; j <= i__1; ++j) {
3810 :     i__2 = *m;
3811 :     for (i__ = 1; i__ <= i__2; ++i__) {
3812 :     i__3 = c___subscr(i__, j);
3813 :     c__[i__3].r = 0., c__[i__3].i = 0.;
3814 :     /* L10: */
3815 :     }
3816 :     /* L20: */
3817 :     }
3818 :     } else {
3819 :     i__1 = *n;
3820 :     for (j = 1; j <= i__1; ++j) {
3821 :     i__2 = *m;
3822 :     for (i__ = 1; i__ <= i__2; ++i__) {
3823 :     i__3 = c___subscr(i__, j);
3824 :     i__4 = c___subscr(i__, j);
3825 :     z__1.r = beta->r * c__[i__4].r - beta->i * c__[i__4].i,
3826 :     z__1.i = beta->r * c__[i__4].i + beta->i * c__[
3827 :     i__4].r;
3828 :     c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
3829 :     /* L30: */
3830 :     }
3831 :     /* L40: */
3832 :     }
3833 :     }
3834 :     return 0;
3835 :     }
3836 :     /* Start the operations. */
3837 :     if (notb) {
3838 :     if (nota) {
3839 :     /* Form C := alpha*A*B + beta*C. */
3840 :     i__1 = *n;
3841 :     for (j = 1; j <= i__1; ++j) {
3842 :     if (beta->r == 0. && beta->i == 0.) {
3843 :     i__2 = *m;
3844 :     for (i__ = 1; i__ <= i__2; ++i__) {
3845 :     i__3 = c___subscr(i__, j);
3846 :     c__[i__3].r = 0., c__[i__3].i = 0.;
3847 :     /* L50: */
3848 :     }
3849 :     } else if (beta->r != 1. || beta->i != 0.) {
3850 :     i__2 = *m;
3851 :     for (i__ = 1; i__ <= i__2; ++i__) {
3852 :     i__3 = c___subscr(i__, j);
3853 :     i__4 = c___subscr(i__, j);
3854 :     z__1.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
3855 :     .i, z__1.i = beta->r * c__[i__4].i + beta->i *
3856 :     c__[i__4].r;
3857 :     c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
3858 :     /* L60: */
3859 :     }
3860 :     }
3861 :     i__2 = *k;
3862 :     for (l = 1; l <= i__2; ++l) {
3863 :     i__3 = b_subscr(l, j);
3864 :     if (b[i__3].r != 0. || b[i__3].i != 0.) {
3865 :     i__3 = b_subscr(l, j);
3866 :     z__1.r = alpha->r * b[i__3].r - alpha->i * b[i__3].i,
3867 :     z__1.i = alpha->r * b[i__3].i + alpha->i * b[
3868 :     i__3].r;
3869 :     temp.r = z__1.r, temp.i = z__1.i;
3870 :     i__3 = *m;
3871 :     for (i__ = 1; i__ <= i__3; ++i__) {
3872 :     i__4 = c___subscr(i__, j);
3873 :     i__5 = c___subscr(i__, j);
3874 :     i__6 = a_subscr(i__, l);
3875 :     z__2.r = temp.r * a[i__6].r - temp.i * a[i__6].i,
3876 :     z__2.i = temp.r * a[i__6].i + temp.i * a[
3877 :     i__6].r;
3878 :     z__1.r = c__[i__5].r + z__2.r, z__1.i = c__[i__5]
3879 :     .i + z__2.i;
3880 :     c__[i__4].r = z__1.r, c__[i__4].i = z__1.i;
3881 :     /* L70: */
3882 :     }
3883 :     }
3884 :     /* L80: */
3885 :     }
3886 :     /* L90: */
3887 :     }
3888 :     } else if (conja) {
3889 :     /* Form C := alpha*conjg( A' )*B + beta*C. */
3890 :     i__1 = *n;
3891 :     for (j = 1; j <= i__1; ++j) {
3892 :     i__2 = *m;
3893 :     for (i__ = 1; i__ <= i__2; ++i__) {
3894 :     temp.r = 0., temp.i = 0.;
3895 :     i__3 = *k;
3896 :     for (l = 1; l <= i__3; ++l) {
3897 :     d_cnjg(&z__3, &a_ref(l, i__));
3898 :     i__4 = b_subscr(l, j);
3899 :     z__2.r = z__3.r * b[i__4].r - z__3.i * b[i__4].i,
3900 :     z__2.i = z__3.r * b[i__4].i + z__3.i * b[i__4]
3901 :     .r;
3902 :     z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
3903 :     temp.r = z__1.r, temp.i = z__1.i;
3904 :     /* L100: */
3905 :     }
3906 :     if (beta->r == 0. && beta->i == 0.) {
3907 :     i__3 = c___subscr(i__, j);
3908 :     z__1.r = alpha->r * temp.r - alpha->i * temp.i,
3909 :     z__1.i = alpha->r * temp.i + alpha->i *
3910 :     temp.r;
3911 :     c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
3912 :     } else {
3913 :     i__3 = c___subscr(i__, j);
3914 :     z__2.r = alpha->r * temp.r - alpha->i * temp.i,
3915 :     z__2.i = alpha->r * temp.i + alpha->i *
3916 :     temp.r;
3917 :     i__4 = c___subscr(i__, j);
3918 :     z__3.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
3919 :     .i, z__3.i = beta->r * c__[i__4].i + beta->i *
3920 :     c__[i__4].r;
3921 :     z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
3922 :     c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
3923 :     }
3924 :     /* L110: */
3925 :     }
3926 :     /* L120: */
3927 :     }
3928 :     } else {
3929 :     /* Form C := alpha*A'*B + beta*C */
3930 :     i__1 = *n;
3931 :     for (j = 1; j <= i__1; ++j) {
3932 :     i__2 = *m;
3933 :     for (i__ = 1; i__ <= i__2; ++i__) {
3934 :     temp.r = 0., temp.i = 0.;
3935 :     i__3 = *k;
3936 :     for (l = 1; l <= i__3; ++l) {
3937 :     i__4 = a_subscr(l, i__);
3938 :     i__5 = b_subscr(l, j);
3939 :     z__2.r = a[i__4].r * b[i__5].r - a[i__4].i * b[i__5]
3940 :     .i, z__2.i = a[i__4].r * b[i__5].i + a[i__4]
3941 :     .i * b[i__5].r;
3942 :     z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
3943 :     temp.r = z__1.r, temp.i = z__1.i;
3944 :     /* L130: */
3945 :     }
3946 :     if (beta->r == 0. && beta->i == 0.) {
3947 :     i__3 = c___subscr(i__, j);
3948 :     z__1.r = alpha->r * temp.r - alpha->i * temp.i,
3949 :     z__1.i = alpha->r * temp.i + alpha->i *
3950 :     temp.r;
3951 :     c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
3952 :     } else {
3953 :     i__3 = c___subscr(i__, j);
3954 :     z__2.r = alpha->r * temp.r - alpha->i * temp.i,
3955 :     z__2.i = alpha->r * temp.i + alpha->i *
3956 :     temp.r;
3957 :     i__4 = c___subscr(i__, j);
3958 :     z__3.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
3959 :     .i, z__3.i = beta->r * c__[i__4].i + beta->i *
3960 :     c__[i__4].r;
3961 :     z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
3962 :     c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
3963 :     }
3964 :     /* L140: */
3965 :     }
3966 :     /* L150: */
3967 :     }
3968 :     }
3969 :     } else if (nota) {
3970 :     if (conjb) {
3971 :     /* Form C := alpha*A*conjg( B' ) + beta*C. */
3972 :     i__1 = *n;
3973 :     for (j = 1; j <= i__1; ++j) {
3974 :     if (beta->r == 0. && beta->i == 0.) {
3975 :     i__2 = *m;
3976 :     for (i__ = 1; i__ <= i__2; ++i__) {
3977 :     i__3 = c___subscr(i__, j);
3978 :     c__[i__3].r = 0., c__[i__3].i = 0.;
3979 :     /* L160: */
3980 :     }
3981 :     } else if (beta->r != 1. || beta->i != 0.) {
3982 :     i__2 = *m;
3983 :     for (i__ = 1; i__ <= i__2; ++i__) {
3984 :     i__3 = c___subscr(i__, j);
3985 :     i__4 = c___subscr(i__, j);
3986 :     z__1.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
3987 :     .i, z__1.i = beta->r * c__[i__4].i + beta->i *
3988 :     c__[i__4].r;
3989 :     c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
3990 :     /* L170: */
3991 :     }
3992 :     }
3993 :     i__2 = *k;
3994 :     for (l = 1; l <= i__2; ++l) {
3995 :     i__3 = b_subscr(j, l);
3996 :     if (b[i__3].r != 0. || b[i__3].i != 0.) {
3997 :     d_cnjg(&z__2, &b_ref(j, l));
3998 :     z__1.r = alpha->r * z__2.r - alpha->i * z__2.i,
3999 :     z__1.i = alpha->r * z__2.i + alpha->i *
4000 :     z__2.r;
4001 :     temp.r = z__1.r, temp.i = z__1.i;
4002 :     i__3 = *m;
4003 :     for (i__ = 1; i__ <= i__3; ++i__) {
4004 :     i__4 = c___subscr(i__, j);
4005 :     i__5 = c___subscr(i__, j);
4006 :     i__6 = a_subscr(i__, l);
4007 :     z__2.r = temp.r * a[i__6].r - temp.i * a[i__6].i,
4008 :     z__2.i = temp.r * a[i__6].i + temp.i * a[
4009 :     i__6].r;
4010 :     z__1.r = c__[i__5].r + z__2.r, z__1.i = c__[i__5]
4011 :     .i + z__2.i;
4012 :     c__[i__4].r = z__1.r, c__[i__4].i = z__1.i;
4013 :     /* L180: */
4014 :     }
4015 :     }
4016 :     /* L190: */
4017 :     }
4018 :     /* L200: */
4019 :     }
4020 :     } else {
4021 :     /* Form C := alpha*A*B' + beta*C */
4022 :     i__1 = *n;
4023 :     for (j = 1; j <= i__1; ++j) {
4024 :     if (beta->r == 0. && beta->i == 0.) {
4025 :     i__2 = *m;
4026 :     for (i__ = 1; i__ <= i__2; ++i__) {
4027 :     i__3 = c___subscr(i__, j);
4028 :     c__[i__3].r = 0., c__[i__3].i = 0.;
4029 :     /* L210: */
4030 :     }
4031 :     } else if (beta->r != 1. || beta->i != 0.) {
4032 :     i__2 = *m;
4033 :     for (i__ = 1; i__ <= i__2; ++i__) {
4034 :     i__3 = c___subscr(i__, j);
4035 :     i__4 = c___subscr(i__, j);
4036 :     z__1.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
4037 :     .i, z__1.i = beta->r * c__[i__4].i + beta->i *
4038 :     c__[i__4].r;
4039 :     c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
4040 :     /* L220: */
4041 :     }
4042 :     }
4043 :     i__2 = *k;
4044 :     for (l = 1; l <= i__2; ++l) {
4045 :     i__3 = b_subscr(j, l);
4046 :     if (b[i__3].r != 0. || b[i__3].i != 0.) {
4047 :     i__3 = b_subscr(j, l);
4048 :     z__1.r = alpha->r * b[i__3].r - alpha->i * b[i__3].i,
4049 :     z__1.i = alpha->r * b[i__3].i + alpha->i * b[
4050 :     i__3].r;
4051 :     temp.r = z__1.r, temp.i = z__1.i;
4052 :     i__3 = *m;
4053 :     for (i__ = 1; i__ <= i__3; ++i__) {
4054 :     i__4 = c___subscr(i__, j);
4055 :     i__5 = c___subscr(i__, j);
4056 :     i__6 = a_subscr(i__, l);
4057 :     z__2.r = temp.r * a[i__6].r - temp.i * a[i__6].i,
4058 :     z__2.i = temp.r * a[i__6].i + temp.i * a[
4059 :     i__6].r;
4060 :     z__1.r = c__[i__5].r + z__2.r, z__1.i = c__[i__5]
4061 :     .i + z__2.i;
4062 :     c__[i__4].r = z__1.r, c__[i__4].i = z__1.i;
4063 :     /* L230: */
4064 :     }
4065 :     }
4066 :     /* L240: */
4067 :     }
4068 :     /* L250: */
4069 :     }
4070 :     }
4071 :     } else if (conja) {
4072 :     if (conjb) {
4073 :     /* Form C := alpha*conjg( A' )*conjg( B' ) + beta*C. */
4074 :     i__1 = *n;
4075 :     for (j = 1; j <= i__1; ++j) {
4076 :     i__2 = *m;
4077 :     for (i__ = 1; i__ <= i__2; ++i__) {
4078 :     temp.r = 0., temp.i = 0.;
4079 :     i__3 = *k;
4080 :     for (l = 1; l <= i__3; ++l) {
4081 :     d_cnjg(&z__3, &a_ref(l, i__));
4082 :     d_cnjg(&z__4, &b_ref(j, l));
4083 :     z__2.r = z__3.r * z__4.r - z__3.i * z__4.i, z__2.i =
4084 :     z__3.r * z__4.i + z__3.i * z__4.r;
4085 :     z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
4086 :     temp.r = z__1.r, temp.i = z__1.i;
4087 :     /* L260: */
4088 :     }
4089 :     if (beta->r == 0. && beta->i == 0.) {
4090 :     i__3 = c___subscr(i__, j);
4091 :     z__1.r = alpha->r * temp.r - alpha->i * temp.i,
4092 :     z__1.i = alpha->r * temp.i + alpha->i *
4093 :     temp.r;
4094 :     c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
4095 :     } else {
4096 :     i__3 = c___subscr(i__, j);
4097 :     z__2.r = alpha->r * temp.r - alpha->i * temp.i,
4098 :     z__2.i = alpha->r * temp.i + alpha->i *
4099 :     temp.r;
4100 :     i__4 = c___subscr(i__, j);
4101 :     z__3.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
4102 :     .i, z__3.i = beta->r * c__[i__4].i + beta->i *
4103 :     c__[i__4].r;
4104 :     z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
4105 :     c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
4106 :     }
4107 :     /* L270: */
4108 :     }
4109 :     /* L280: */
4110 :     }
4111 :     } else {
4112 :     /* Form C := alpha*conjg( A' )*B' + beta*C */
4113 :     i__1 = *n;
4114 :     for (j = 1; j <= i__1; ++j) {
4115 :     i__2 = *m;
4116 :     for (i__ = 1; i__ <= i__2; ++i__) {
4117 :     temp.r = 0., temp.i = 0.;
4118 :     i__3 = *k;
4119 :     for (l = 1; l <= i__3; ++l) {
4120 :     d_cnjg(&z__3, &a_ref(l, i__));
4121 :     i__4 = b_subscr(j, l);
4122 :     z__2.r = z__3.r * b[i__4].r - z__3.i * b[i__4].i,
4123 :     z__2.i = z__3.r * b[i__4].i + z__3.i * b[i__4]
4124 :     .r;
4125 :     z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
4126 :     temp.r = z__1.r, temp.i = z__1.i;
4127 :     /* L290: */
4128 :     }
4129 :     if (beta->r == 0. && beta->i == 0.) {
4130 :     i__3 = c___subscr(i__, j);
4131 :     z__1.r = alpha->r * temp.r - alpha->i * temp.i,
4132 :     z__1.i = alpha->r * temp.i + alpha->i *
4133 :     temp.r;
4134 :     c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
4135 :     } else {
4136 :     i__3 = c___subscr(i__, j);
4137 :     z__2.r = alpha->r * temp.r - alpha->i * temp.i,
4138 :     z__2.i = alpha->r * temp.i + alpha->i *
4139 :     temp.r;
4140 :     i__4 = c___subscr(i__, j);
4141 :     z__3.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
4142 :     .i, z__3.i = beta->r * c__[i__4].i + beta->i *
4143 :     c__[i__4].r;
4144 :     z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
4145 :     c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
4146 :     }
4147 :     /* L300: */
4148 :     }
4149 :     /* L310: */
4150 :     }
4151 :     }
4152 :     } else {
4153 :     if (conjb) {
4154 :     /* Form C := alpha*A'*conjg( B' ) + beta*C */
4155 :     i__1 = *n;
4156 :     for (j = 1; j <= i__1; ++j) {
4157 :     i__2 = *m;
4158 :     for (i__ = 1; i__ <= i__2; ++i__) {
4159 :     temp.r = 0., temp.i = 0.;
4160 :     i__3 = *k;
4161 :     for (l = 1; l <= i__3; ++l) {
4162 :     i__4 = a_subscr(l, i__);
4163 :     d_cnjg(&z__3, &b_ref(j, l));
4164 :     z__2.r = a[i__4].r * z__3.r - a[i__4].i * z__3.i,
4165 :     z__2.i = a[i__4].r * z__3.i + a[i__4].i *
4166 :     z__3.r;
4167 :     z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
4168 :     temp.r = z__1.r, temp.i = z__1.i;
4169 :     /* L320: */
4170 :     }
4171 :     if (beta->r == 0. && beta->i == 0.) {
4172 :     i__3 = c___subscr(i__, j);
4173 :     z__1.r = alpha->r * temp.r - alpha->i * temp.i,
4174 :     z__1.i = alpha->r * temp.i + alpha->i *
4175 :     temp.r;
4176 :     c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
4177 :     } else {
4178 :     i__3 = c___subscr(i__, j);
4179 :     z__2.r = alpha->r * temp.r - alpha->i * temp.i,
4180 :     z__2.i = alpha->r * temp.i + alpha->i *
4181 :     temp.r;
4182 :     i__4 = c___subscr(i__, j);
4183 :     z__3.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
4184 :     .i, z__3.i = beta->r * c__[i__4].i + beta->i *
4185 :     c__[i__4].r;
4186 :     z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
4187 :     c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
4188 :     }
4189 :     /* L330: */
4190 :     }
4191 :     /* L340: */
4192 :     }
4193 :     } else {
4194 :     /* Form C := alpha*A'*B' + beta*C */
4195 :     i__1 = *n;
4196 :     for (j = 1; j <= i__1; ++j) {
4197 :     i__2 = *m;
4198 :     for (i__ = 1; i__ <= i__2; ++i__) {
4199 :     temp.r = 0., temp.i = 0.;
4200 :     i__3 = *k;
4201 :     for (l = 1; l <= i__3; ++l) {
4202 :     i__4 = a_subscr(l, i__);
4203 :     i__5 = b_subscr(j, l);
4204 :     z__2.r = a[i__4].r * b[i__5].r - a[i__4].i * b[i__5]
4205 :     .i, z__2.i = a[i__4].r * b[i__5].i + a[i__4]
4206 :     .i * b[i__5].r;
4207 :     z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
4208 :     temp.r = z__1.r, temp.i = z__1.i;
4209 :     /* L350: */
4210 :     }
4211 :     if (beta->r == 0. && beta->i == 0.) {
4212 :     i__3 = c___subscr(i__, j);
4213 :     z__1.r = alpha->r * temp.r - alpha->i * temp.i,
4214 :     z__1.i = alpha->r * temp.i + alpha->i *
4215 :     temp.r;
4216 :     c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
4217 :     } else {
4218 :     i__3 = c___subscr(i__, j);
4219 :     z__2.r = alpha->r * temp.r - alpha->i * temp.i,
4220 :     z__2.i = alpha->r * temp.i + alpha->i *
4221 :     temp.r;
4222 :     i__4 = c___subscr(i__, j);
4223 :     z__3.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
4224 :     .i, z__3.i = beta->r * c__[i__4].i + beta->i *
4225 :     c__[i__4].r;
4226 :     z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
4227 :     c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
4228 :     }
4229 :     /* L360: */
4230 :     }
4231 :     /* L370: */
4232 :     }
4233 :     }
4234 :     }
4235 :     return 0;
4236 :     /* End of ZGEMM . */
4237 :     } /* zgemm_ */
4238 :     #undef c___ref
4239 :     #undef c___subscr
4240 :     #undef b_ref
4241 :     #undef b_subscr
4242 :     #undef a_ref
4243 :     #undef a_subscr
4244 :    
4245 :    
4246 :     /* Subroutine */ int zgemv_(char *trans, integer *m, integer *n,
4247 :     doublecomplex *alpha, doublecomplex *a, integer *lda, doublecomplex *
4248 :     x, integer *incx, doublecomplex *beta, doublecomplex *y, integer *
4249 :     incy)
4250 :     {
4251 :     /* System generated locals */
4252 :     integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
4253 :     doublecomplex z__1, z__2, z__3;
4254 :     /* Builtin functions */
4255 :     void d_cnjg(doublecomplex *, doublecomplex *);
4256 :     /* Local variables */
4257 :     static integer info;
4258 :     static doublecomplex temp;
4259 :     static integer lenx, leny, i__, j;
4260 :     extern logical lsame_(char *, char *);
4261 :     static integer ix, iy, jx, jy, kx, ky;
4262 :     extern /* Subroutine */ int xerbla_(char *, integer *);
4263 :     static logical noconj;
4264 :     #define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1
4265 :     #define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]
4266 :     /* Purpose
4267 :     =======
4268 :     ZGEMV performs one of the matrix-vector operations
4269 :     y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, or
4270 :     y := alpha*conjg( A' )*x + beta*y,
4271 :     where alpha and beta are scalars, x and y are vectors and A is an
4272 :     m by n matrix.
4273 :     Parameters
4274 :     ==========
4275 :     TRANS - CHARACTER*1.
4276 :     On entry, TRANS specifies the operation to be performed as
4277 :     follows:
4278 :     TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
4279 :     TRANS = 'T' or 't' y := alpha*A'*x + beta*y.
4280 :     TRANS = 'C' or 'c' y := alpha*conjg( A' )*x + beta*y.
4281 :     Unchanged on exit.
4282 :     M - INTEGER.
4283 :     On entry, M specifies the number of rows of the matrix A.
4284 :     M must be at least zero.
4285 :     Unchanged on exit.
4286 :     N - INTEGER.
4287 :     On entry, N specifies the number of columns of the matrix A.
4288 :     N must be at least zero.
4289 :     Unchanged on exit.
4290 :     ALPHA - COMPLEX*16 .
4291 :     On entry, ALPHA specifies the scalar alpha.
4292 :     Unchanged on exit.
4293 :     A - COMPLEX*16 array of DIMENSION ( LDA, n ).
4294 :     Before entry, the leading m by n part of the array A must
4295 :     contain the matrix of coefficients.
4296 :     Unchanged on exit.
4297 :     LDA - INTEGER.
4298 :     On entry, LDA specifies the first dimension of A as declared
4299 :     in the calling (sub) program. LDA must be at least
4300 :     max( 1, m ).
4301 :     Unchanged on exit.
4302 :     X - COMPLEX*16 array of DIMENSION at least
4303 :     ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
4304 :     and at least
4305 :     ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
4306 :     Before entry, the incremented array X must contain the
4307 :     vector x.
4308 :     Unchanged on exit.
4309 :     INCX - INTEGER.
4310 :     On entry, INCX specifies the increment for the elements of
4311 :     X. INCX must not be zero.
4312 :     Unchanged on exit.
4313 :     BETA - COMPLEX*16 .
4314 :     On entry, BETA specifies the scalar beta. When BETA is
4315 :     supplied as zero then Y need not be set on input.
4316 :     Unchanged on exit.
4317 :     Y - COMPLEX*16 array of DIMENSION at least
4318 :     ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
4319 :     and at least
4320 :     ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
4321 :     Before entry with BETA non-zero, the incremented array Y
4322 :     must contain the vector y. On exit, Y is overwritten by the
4323 :     updated vector y.
4324 :     INCY - INTEGER.
4325 :     On entry, INCY specifies the increment for the elements of
4326 :     Y. INCY must not be zero.
4327 :     Unchanged on exit.
4328 :     Level 2 Blas routine.
4329 :     -- Written on 22-October-1986.
4330 :     Jack Dongarra, Argonne National Lab.
4331 :     Jeremy Du Croz, Nag Central Office.
4332 :     Sven Hammarling, Nag Central Office.
4333 :     Richard Hanson, Sandia National Labs.
4334 :     Test the input parameters.
4335 :     Parameter adjustments */
4336 :     a_dim1 = *lda;
4337 :     a_offset = 1 + a_dim1 * 1;
4338 :     a -= a_offset;
4339 :     --x;
4340 :     --y;
4341 :     /* Function Body */
4342 :     info = 0;
4343 :     if (! lsame_(trans, "N") && ! lsame_(trans, "T") && ! lsame_(trans, "C")
4344 :     ) {
4345 :     info = 1;
4346 :     } else if (*m < 0) {
4347 :     info = 2;
4348 :     } else if (*n < 0) {
4349 :     info = 3;
4350 :     } else if (*lda < max(1,*m)) {
4351 :     info = 6;
4352 :     } else if (*incx == 0) {
4353 :     info = 8;
4354 :     } else if (*incy == 0) {
4355 :     info = 11;
4356 :     }
4357 :     if (info != 0) {
4358 :     xerbla_("ZGEMV ", &info);
4359 :     return 0;
4360 :     }
4361 :     /* Quick return if possible. */
4362 :     if (*m == 0 || *n == 0 || alpha->r == 0. && alpha->i == 0. && (beta->r ==
4363 :     1. && beta->i == 0.)) {
4364 :     return 0;
4365 :     }
4366 :     noconj = lsame_(trans, "T");
4367 :     /* Set LENX and LENY, the lengths of the vectors x and y, and set
4368 :     up the start points in X and Y. */
4369 :     if (lsame_(trans, "N")) {
4370 :     lenx = *n;
4371 :     leny = *m;
4372 :     } else {
4373 :     lenx = *m;
4374 :     leny = *n;
4375 :     }
4376 :     if (*incx > 0) {
4377 :     kx = 1;
4378 :     } else {
4379 :     kx = 1 - (lenx - 1) * *incx;
4380 :     }
4381 :     if (*incy > 0) {
4382 :     ky = 1;
4383 :     } else {
4384 :     ky = 1 - (leny - 1) * *incy;
4385 :     }
4386 :     /* Start the operations. In this version the elements of A are
4387 :     accessed sequentially with one pass through A.
4388 :     First form y := beta*y. */
4389 :     if (beta->r != 1. || beta->i != 0.) {
4390 :     if (*incy == 1) {
4391 :     if (beta->r == 0. && beta->i == 0.) {
4392 :     i__1 = leny;
4393 :     for (i__ = 1; i__ <= i__1; ++i__) {
4394 :     i__2 = i__;
4395 :     y[i__2].r = 0., y[i__2].i = 0.;
4396 :     /* L10: */
4397 :     }
4398 :     } else {
4399 :     i__1 = leny;
4400 :     for (i__ = 1; i__ <= i__1; ++i__) {
4401 :     i__2 = i__;
4402 :     i__3 = i__;
4403 :     z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i,
4404 :     z__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
4405 :     .r;
4406 :     y[i__2].r = z__1.r, y[i__2].i = z__1.i;
4407 :     /* L20: */
4408 :     }
4409 :     }
4410 :     } else {
4411 :     iy = ky;
4412 :     if (beta->r == 0. && beta->i == 0.) {
4413 :     i__1 = leny;
4414 :     for (i__ = 1; i__ <= i__1; ++i__) {
4415 :     i__2 = iy;
4416 :     y[i__2].r = 0., y[i__2].i = 0.;
4417 :     iy += *incy;
4418 :     /* L30: */
4419 :     }
4420 :     } else {
4421 :     i__1 = leny;
4422 :     for (i__ = 1; i__ <= i__1; ++i__) {
4423 :     i__2 = iy;
4424 :     i__3 = iy;
4425 :     z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i,
4426 :     z__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
4427 :     .r;
4428 :     y[i__2].r = z__1.r, y[i__2].i = z__1.i;
4429 :     iy += *incy;
4430 :     /* L40: */
4431 :     }
4432 :     }
4433 :     }
4434 :     }
4435 :     if (alpha->r == 0. && alpha->i == 0.) {
4436 :     return 0;
4437 :     }
4438 :     if (lsame_(trans, "N")) {
4439 :     /* Form y := alpha*A*x + y. */
4440 :     jx = kx;
4441 :     if (*incy == 1) {
4442 :     i__1 = *n;
4443 :     for (j = 1; j <= i__1; ++j) {
4444 :     i__2 = jx;
4445 :     if (x[i__2].r != 0. || x[i__2].i != 0.) {
4446 :     i__2 = jx;
4447 :     z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i,
4448 :     z__1.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
4449 :     .r;
4450 :     temp.r = z__1.r, temp.i = z__1.i;
4451 :     i__2 = *m;
4452 :     for (i__ = 1; i__ <= i__2; ++i__) {
4453 :     i__3 = i__;
4454 :     i__4 = i__;
4455 :     i__5 = a_subscr(i__, j);
4456 :     z__2.r = temp.r * a[i__5].r - temp.i * a[i__5].i,
4457 :     z__2.i = temp.r * a[i__5].i + temp.i * a[i__5]
4458 :     .r;
4459 :     z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i +
4460 :     z__2.i;
4461 :     y[i__3].r = z__1.r, y[i__3].i = z__1.i;
4462 :     /* L50: */
4463 :     }
4464 :     }
4465 :     jx += *incx;
4466 :     /* L60: */
4467 :     }
4468 :     } else {
4469 :     i__1 = *n;
4470 :     for (j = 1; j <= i__1; ++j) {
4471 :     i__2 = jx;
4472 :     if (x[i__2].r != 0. || x[i__2].i != 0.) {
4473 :     i__2 = jx;
4474 :     z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i,
4475 :     z__1.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
4476 :     .r;
4477 :     temp.r = z__1.r, temp.i = z__1.i;
4478 :     iy = ky;
4479 :     i__2 = *m;
4480 :     for (i__ = 1; i__ <= i__2; ++i__) {
4481 :     i__3 = iy;
4482 :     i__4 = iy;
4483 :     i__5 = a_subscr(i__, j);
4484 :     z__2.r = temp.r * a[i__5].r - temp.i * a[i__5].i,
4485 :     z__2.i = temp.r * a[i__5].i + temp.i * a[i__5]
4486 :     .r;
4487 :     z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i +
4488 :     z__2.i;
4489 :     y[i__3].r = z__1.r, y[i__3].i = z__1.i;
4490 :     iy += *incy;
4491 :     /* L70: */
4492 :     }
4493 :     }
4494 :     jx += *incx;
4495 :     /* L80: */
4496 :     }
4497 :     }
4498 :     } else {
4499 :     /* Form y := alpha*A'*x + y or y := alpha*conjg( A' )*x + y. */
4500 :     jy = ky;
4501 :     if (*incx == 1) {
4502 :     i__1 = *n;
4503 :     for (j = 1; j <= i__1; ++j) {
4504 :     temp.r = 0., temp.i = 0.;
4505 :     if (noconj) {
4506 :     i__2 = *m;
4507 :     for (i__ = 1; i__ <= i__2; ++i__) {
4508 :     i__3 = a_subscr(i__, j);
4509 :     i__4 = i__;
4510 :     z__2.r = a[i__3].r * x[i__4].r - a[i__3].i * x[i__4]
4511 :     .i, z__2.i = a[i__3].r * x[i__4].i + a[i__3]
4512 :     .i * x[i__4].r;
4513 :     z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
4514 :     temp.r = z__1.r, temp.i = z__1.i;
4515 :     /* L90: */
4516 :     }
4517 :     } else {
4518 :     i__2 = *m;
4519 :     for (i__ = 1; i__ <= i__2; ++i__) {
4520 :     d_cnjg(&z__3, &a_ref(i__, j));
4521 :     i__3 = i__;
4522 :     z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i,
4523 :     z__2.i = z__3.r * x[i__3].i + z__3.i * x[i__3]
4524 :     .r;
4525 :     z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
4526 :     temp.r = z__1.r, temp.i = z__1.i;
4527 :     /* L100: */
4528 :     }
4529 :     }
4530 :     i__2 = jy;
4531 :     i__3 = jy;
4532 :     z__2.r = alpha->r * temp.r - alpha->i * temp.i, z__2.i =
4533 :     alpha->r * temp.i + alpha->i * temp.r;
4534 :     z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i;
4535 :     y[i__2].r = z__1.r, y[i__2].i = z__1.i;
4536 :     jy += *incy;
4537 :     /* L110: */
4538 :     }
4539 :     } else {
4540 :     i__1 = *n;
4541 :     for (j = 1; j <= i__1; ++j) {
4542 :     temp.r = 0., temp.i = 0.;
4543 :     ix = kx;
4544 :     if (noconj) {
4545 :     i__2 = *m;
4546 :     for (i__ = 1; i__ <= i__2; ++i__) {
4547 :     i__3 = a_subscr(i__, j);
4548 :     i__4 = ix;
4549 :     z__2.r = a[i__3].r * x[i__4].r - a[i__3].i * x[i__4]
4550 :     .i, z__2.i = a[i__3].r * x[i__4].i + a[i__3]
4551 :     .i * x[i__4].r;
4552 :     z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
4553 :     temp.r = z__1.r, temp.i = z__1.i;
4554 :     ix += *incx;
4555 :     /* L120: */
4556 :     }
4557 :     } else {
4558 :     i__2 = *m;
4559 :     for (i__ = 1; i__ <= i__2; ++i__) {
4560 :     d_cnjg(&z__3, &a_ref(i__, j));
4561 :     i__3 = ix;
4562 :     z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i,
4563 :     z__2.i = z__3.r * x[i__3].i + z__3.i * x[i__3]
4564 :     .r;
4565 :     z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
4566 :     temp.r = z__1.r, temp.i = z__1.i;
4567 :     ix += *incx;
4568 :     /* L130: */
4569 :     }
4570 :     }
4571 :     i__2 = jy;
4572 :     i__3 = jy;
4573 :     z__2.r = alpha->r * temp.r - alpha->i * temp.i, z__2.i =
4574 :     alpha->r * temp.i + alpha->i * temp.r;
4575 :     z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i;
4576 :     y[i__2].r = z__1.r, y[i__2].i = z__1.i;
4577 :     jy += *incy;
4578 :     /* L140: */
4579 :     }
4580 :     }
4581 :     }
4582 :     return 0;
4583 :     /* End of ZGEMV . */
4584 :     } /* zgemv_ */
4585 :     #undef a_ref
4586 :     #undef a_subscr
4587 :    
4588 :    
4589 :     /* Subroutine */ int zgerc_(integer *m, integer *n, doublecomplex *alpha,
4590 :     doublecomplex *x, integer *incx, doublecomplex *y, integer *incy,
4591 :     doublecomplex *a, integer *lda)
4592 :     {
4593 :     /* System generated locals */
4594 :     integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
4595 :     doublecomplex z__1, z__2;
4596 :     /* Builtin functions */
4597 :     void d_cnjg(doublecomplex *, doublecomplex *);
4598 :     /* Local variables */
4599 :     static integer info;
4600 :     static doublecomplex temp;
4601 :     static integer i__, j, ix, jy, kx;
4602 :     extern /* Subroutine */ int xerbla_(char *, integer *);
4603 :     #define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1
4604 :     #define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]
4605 :     /* Purpose
4606 :     =======
4607 :     ZGERC performs the rank 1 operation
4608 :     A := alpha*x*conjg( y' ) + A,
4609 :     where alpha is a scalar, x is an m element vector, y is an n element
4610 :     vector and A is an m by n matrix.
4611 :     Parameters
4612 :     ==========
4613 :     M - INTEGER.
4614 :     On entry, M specifies the number of rows of the matrix A.
4615 :     M must be at least zero.
4616 :     Unchanged on exit.
4617 :     N - INTEGER.
4618 :     On entry, N specifies the number of columns of the matrix A.
4619 :     N must be at least zero.
4620 :     Unchanged on exit.
4621 :     ALPHA - COMPLEX*16 .
4622 :     On entry, ALPHA specifies the scalar alpha.
4623 :     Unchanged on exit.
4624 :     X - COMPLEX*16 array of dimension at least
4625 :     ( 1 + ( m - 1 )*abs( INCX ) ).
4626 :     Before entry, the incremented array X must contain the m
4627 :     element vector x.
4628 :     Unchanged on exit.
4629 :     INCX - INTEGER.
4630 :     On entry, INCX specifies the increment for the elements of
4631 :     X. INCX must not be zero.
4632 :     Unchanged on exit.
4633 :     Y - COMPLEX*16 array of dimension at least
4634 :     ( 1 + ( n - 1 )*abs( INCY ) ).
4635 :     Before entry, the incremented array Y must contain the n
4636 :     element vector y.
4637 :     Unchanged on exit.
4638 :     INCY - INTEGER.
4639 :     On entry, INCY specifies the increment for the elements of
4640 :     Y. INCY must not be zero.
4641 :     Unchanged on exit.
4642 :     A - COMPLEX*16 array of DIMENSION ( LDA, n ).
4643 :     Before entry, the leading m by n part of the array A must
4644 :     contain the matrix of coefficients. On exit, A is
4645 :     overwritten by the updated matrix.
4646 :     LDA - INTEGER.
4647 :     On entry, LDA specifies the first dimension of A as declared
4648 :     in the calling (sub) program. LDA must be at least
4649 :     max( 1, m ).
4650 :     Unchanged on exit.
4651 :     Level 2 Blas routine.
4652 :     -- Written on 22-October-1986.
4653 :     Jack Dongarra, Argonne National Lab.
4654 :     Jeremy Du Croz, Nag Central Office.
4655 :     Sven Hammarling, Nag Central Office.
4656 :     Richard Hanson, Sandia National Labs.
4657 :     Test the input parameters.
4658 :     Parameter adjustments */
4659 :     --x;
4660 :     --y;
4661 :     a_dim1 = *lda;
4662 :     a_offset = 1 + a_dim1 * 1;
4663 :     a -= a_offset;
4664 :     /* Function Body */
4665 :     info = 0;
4666 :     if (*m < 0) {
4667 :     info = 1;
4668 :     } else if (*n < 0) {
4669 :     info = 2;
4670 :     } else if (*incx == 0) {
4671 :     info = 5;
4672 :     } else if (*incy == 0) {
4673 :     info = 7;
4674 :     } else if (*lda < max(1,*m)) {
4675 :     info = 9;
4676 :     }
4677 :     if (info != 0) {
4678 :     xerbla_("ZGERC ", &info);
4679 :     return 0;
4680 :     }
4681 :     /* Quick return if possible. */
4682 :     if (*m == 0 || *n == 0 || alpha->r == 0. && alpha->i == 0.) {
4683 :     return 0;
4684 :     }
4685 :     /* Start the operations. In this version the elements of A are
4686 :     accessed sequentially with one pass through A. */
4687 :     if (*incy > 0) {
4688 :     jy = 1;
4689 :     } else {
4690 :     jy = 1 - (*n - 1) * *incy;
4691 :     }
4692 :     if (*incx == 1) {
4693 :     i__1 = *n;
4694 :     for (j = 1; j <= i__1; ++j) {
4695 :     i__2 = jy;
4696 :     if (y[i__2].r != 0. || y[i__2].i != 0.) {
4697 :     d_cnjg(&z__2, &y[jy]);
4698 :     z__1.r = alpha->r * z__2.r - alpha->i * z__2.i, z__1.i =
4699 :     alpha->r * z__2.i + alpha->i * z__2.r;
4700 :     temp.r = z__1.r, temp.i = z__1.i;
4701 :     i__2 = *m;
4702 :     for (i__ = 1; i__ <= i__2; ++i__) {
4703 :     i__3 = a_subscr(i__, j);
4704 :     i__4 = a_subscr(i__, j);
4705 :     i__5 = i__;
4706 :     z__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, z__2.i =
4707 :     x[i__5].r * temp.i + x[i__5].i * temp.r;
4708 :     z__1.r = a[i__4].r + z__2.r, z__1.i = a[i__4].i + z__2.i;
4709 :     a[i__3].r = z__1.r, a[i__3].i = z__1.i;
4710 :     /* L10: */
4711 :     }
4712 :     }
4713 :     jy += *incy;
4714 :     /* L20: */
4715 :     }
4716 :     } else {
4717 :     if (*incx > 0) {
4718 :     kx = 1;
4719 :     } else {
4720 :     kx = 1 - (*m - 1) * *incx;
4721 :     }
4722 :     i__1 = *n;
4723 :     for (j = 1; j <= i__1; ++j) {
4724 :     i__2 = jy;
4725 :     if (y[i__2].r != 0. || y[i__2].i != 0.) {
4726 :     d_cnjg(&z__2, &y[jy]);
4727 :     z__1.r = alpha->r * z__2.r - alpha->i * z__2.i, z__1.i =
4728 :     alpha->r * z__2.i + alpha->i * z__2.r;
4729 :     temp.r = z__1.r, temp.i = z__1.i;
4730 :     ix = kx;
4731 :     i__2 = *m;
4732 :     for (i__ = 1; i__ <= i__2; ++i__) {
4733 :     i__3 = a_subscr(i__, j);
4734 :     i__4 = a_subscr(i__, j);
4735 :     i__5 = ix;
4736 :     z__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, z__2.i =
4737 :     x[i__5].r * temp.i + x[i__5].i * temp.r;
4738 :     z__1.r = a[i__4].r + z__2.r, z__1.i = a[i__4].i + z__2.i;
4739 :     a[i__3].r = z__1.r, a[i__3].i = z__1.i;
4740 :     ix += *incx;
4741 :     /* L30: */
4742 :     }
4743 :     }
4744 :     jy += *incy;
4745 :     /* L40: */
4746 :     }
4747 :     }
4748 :     return 0;
4749 :     /* End of ZGERC . */
4750 :     } /* zgerc_ */
4751 :     #undef a_ref
4752 :     #undef a_subscr
4753 :    
4754 :    
4755 :     /* Subroutine */ int zgeru_(integer *m, integer *n, doublecomplex *alpha,
4756 :     doublecomplex *x, integer *incx, doublecomplex *y, integer *incy,
4757 :     doublecomplex *a, integer *lda)
4758 :     {
4759 :     /* System generated locals */