SphinxBase 5prealpha
blas_lite.c
1/*
2NOTE: This is generated code. Look in README.python for information on
3 remaking this file.
4*/
5#include "sphinxbase/f2c.h"
6
7#ifdef HAVE_CONFIG
8#include "config.h"
9#else
10extern doublereal slamch_(char *);
11#define EPSILON slamch_("Epsilon")
12#define SAFEMINIMUM slamch_("Safe minimum")
13#define PRECISION slamch_("Precision")
14#define BASE slamch_("Base")
15#endif
16
17
18extern doublereal slapy2_(real *, real *);
19
20
21
22/* Table of constant values */
23
24static integer c__1 = 1;
25
26logical lsame_(char *ca, char *cb)
27{
28 /* System generated locals */
29 logical ret_val;
30
31 /* Local variables */
32 static integer inta, intb, zcode;
33
34
35/*
36 -- LAPACK auxiliary routine (version 3.0) --
37 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
38 Courant Institute, Argonne National Lab, and Rice University
39 September 30, 1994
40
41
42 Purpose
43 =======
44
45 LSAME returns .TRUE. if CA is the same letter as CB regardless of
46 case.
47
48 Arguments
49 =========
50
51 CA (input) CHARACTER*1
52 CB (input) CHARACTER*1
53 CA and CB specify the single characters to be compared.
54
55 =====================================================================
56
57
58 Test if the characters are equal
59*/
60
61 ret_val = *(unsigned char *)ca == *(unsigned char *)cb;
62 if (ret_val) {
63 return ret_val;
64 }
65
66/* Now test for equivalence if both characters are alphabetic. */
67
68 zcode = 'Z';
69
70/*
71 Use 'Z' rather than 'A' so that ASCII can be detected on Prime
72 machines, on which ICHAR returns a value with bit 8 set.
73 ICHAR('A') on Prime machines returns 193 which is the same as
74 ICHAR('A') on an EBCDIC machine.
75*/
76
77 inta = *(unsigned char *)ca;
78 intb = *(unsigned char *)cb;
79
80 if (zcode == 90 || zcode == 122) {
81
82/*
83 ASCII is assumed - ZCODE is the ASCII code of either lower or
84 upper case 'Z'.
85*/
86
87 if (inta >= 97 && inta <= 122) {
88 inta += -32;
89 }
90 if (intb >= 97 && intb <= 122) {
91 intb += -32;
92 }
93
94 } else if (zcode == 233 || zcode == 169) {
95
96/*
97 EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
98 upper case 'Z'.
99*/
100
101 if (inta >= 129 && inta <= 137 || inta >= 145 && inta <= 153 || inta
102 >= 162 && inta <= 169) {
103 inta += 64;
104 }
105 if (intb >= 129 && intb <= 137 || intb >= 145 && intb <= 153 || intb
106 >= 162 && intb <= 169) {
107 intb += 64;
108 }
109
110 } else if (zcode == 218 || zcode == 250) {
111
112/*
113 ASCII is assumed, on Prime machines - ZCODE is the ASCII code
114 plus 128 of either lower or upper case 'Z'.
115*/
116
117 if (inta >= 225 && inta <= 250) {
118 inta += -32;
119 }
120 if (intb >= 225 && intb <= 250) {
121 intb += -32;
122 }
123 }
124 ret_val = inta == intb;
125
126/*
127 RETURN
128
129 End of LSAME
130*/
131
132 return ret_val;
133} /* lsame_ */
134
135doublereal sdot_(integer *n, real *sx, integer *incx, real *sy, integer *incy)
136{
137 /* System generated locals */
138 integer i__1;
139 real ret_val;
140
141 /* Local variables */
142 static integer i__, m, ix, iy, mp1;
143 static real stemp;
144
145
146/*
147 forms the dot product of two vectors.
148 uses unrolled loops for increments equal to one.
149 jack dongarra, linpack, 3/11/78.
150 modified 12/3/93, array(1) declarations changed to array(*)
151*/
152
153
154 /* Parameter adjustments */
155 --sy;
156 --sx;
157
158 /* Function Body */
159 stemp = 0.f;
160 ret_val = 0.f;
161 if (*n <= 0) {
162 return ret_val;
163 }
164 if (*incx == 1 && *incy == 1) {
165 goto L20;
166 }
167
168/*
169 code for unequal increments or equal increments
170 not equal to 1
171*/
172
173 ix = 1;
174 iy = 1;
175 if (*incx < 0) {
176 ix = (-(*n) + 1) * *incx + 1;
177 }
178 if (*incy < 0) {
179 iy = (-(*n) + 1) * *incy + 1;
180 }
181 i__1 = *n;
182 for (i__ = 1; i__ <= i__1; ++i__) {
183 stemp += sx[ix] * sy[iy];
184 ix += *incx;
185 iy += *incy;
186/* L10: */
187 }
188 ret_val = stemp;
189 return ret_val;
190
191/*
192 code for both increments equal to 1
193
194
195 clean-up loop
196*/
197
198L20:
199 m = *n % 5;
200 if (m == 0) {
201 goto L40;
202 }
203 i__1 = m;
204 for (i__ = 1; i__ <= i__1; ++i__) {
205 stemp += sx[i__] * sy[i__];
206/* L30: */
207 }
208 if (*n < 5) {
209 goto L60;
210 }
211L40:
212 mp1 = m + 1;
213 i__1 = *n;
214 for (i__ = mp1; i__ <= i__1; i__ += 5) {
215 stemp = stemp + sx[i__] * sy[i__] + sx[i__ + 1] * sy[i__ + 1] + sx[
216 i__ + 2] * sy[i__ + 2] + sx[i__ + 3] * sy[i__ + 3] + sx[i__ +
217 4] * sy[i__ + 4];
218/* L50: */
219 }
220L60:
221 ret_val = stemp;
222 return ret_val;
223} /* sdot_ */
224
225/* Subroutine */ int sgemm_(char *transa, char *transb, integer *m, integer *
226 n, integer *k, real *alpha, real *a, integer *lda, real *b, integer *
227 ldb, real *beta, real *c__, integer *ldc)
228{
229 /* System generated locals */
230 integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
231 i__3;
232
233 /* Local variables */
234 static integer i__, j, l, info;
235 static logical nota, notb;
236 static real temp;
237 static integer ncola;
238 extern logical lsame_(char *, char *);
239 static integer nrowa, nrowb;
240 extern /* Subroutine */ int xerbla_(char *, integer *);
241
242
243/*
244 Purpose
245 =======
246
247 SGEMM performs one of the matrix-matrix operations
248
249 C := alpha*op( A )*op( B ) + beta*C,
250
251 where op( X ) is one of
252
253 op( X ) = X or op( X ) = X',
254
255 alpha and beta are scalars, and A, B and C are matrices, with op( A )
256 an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
257
258 Parameters
259 ==========
260
261 TRANSA - CHARACTER*1.
262 On entry, TRANSA specifies the form of op( A ) to be used in
263 the matrix multiplication as follows:
264
265 TRANSA = 'N' or 'n', op( A ) = A.
266
267 TRANSA = 'T' or 't', op( A ) = A'.
268
269 TRANSA = 'C' or 'c', op( A ) = A'.
270
271 Unchanged on exit.
272
273 TRANSB - CHARACTER*1.
274 On entry, TRANSB specifies the form of op( B ) to be used in
275 the matrix multiplication as follows:
276
277 TRANSB = 'N' or 'n', op( B ) = B.
278
279 TRANSB = 'T' or 't', op( B ) = B'.
280
281 TRANSB = 'C' or 'c', op( B ) = B'.
282
283 Unchanged on exit.
284
285 M - INTEGER.
286 On entry, M specifies the number of rows of the matrix
287 op( A ) and of the matrix C. M must be at least zero.
288 Unchanged on exit.
289
290 N - INTEGER.
291 On entry, N specifies the number of columns of the matrix
292 op( B ) and the number of columns of the matrix C. N must be
293 at least zero.
294 Unchanged on exit.
295
296 K - INTEGER.
297 On entry, K specifies the number of columns of the matrix
298 op( A ) and the number of rows of the matrix op( B ). K must
299 be at least zero.
300 Unchanged on exit.
301
302 ALPHA - REAL .
303 On entry, ALPHA specifies the scalar alpha.
304 Unchanged on exit.
305
306 A - REAL array of DIMENSION ( LDA, ka ), where ka is
307 k when TRANSA = 'N' or 'n', and is m otherwise.
308 Before entry with TRANSA = 'N' or 'n', the leading m by k
309 part of the array A must contain the matrix A, otherwise
310 the leading k by m part of the array A must contain the
311 matrix A.
312 Unchanged on exit.
313
314 LDA - INTEGER.
315 On entry, LDA specifies the first dimension of A as declared
316 in the calling (sub) program. When TRANSA = 'N' or 'n' then
317 LDA must be at least max( 1, m ), otherwise LDA must be at
318 least max( 1, k ).
319 Unchanged on exit.
320
321 B - REAL array of DIMENSION ( LDB, kb ), where kb is
322 n when TRANSB = 'N' or 'n', and is k otherwise.
323 Before entry with TRANSB = 'N' or 'n', the leading k by n
324 part of the array B must contain the matrix B, otherwise
325 the leading n by k part of the array B must contain the
326 matrix B.
327 Unchanged on exit.
328
329 LDB - INTEGER.
330 On entry, LDB specifies the first dimension of B as declared
331 in the calling (sub) program. When TRANSB = 'N' or 'n' then
332 LDB must be at least max( 1, k ), otherwise LDB must be at
333 least max( 1, n ).
334 Unchanged on exit.
335
336 BETA - REAL .
337 On entry, BETA specifies the scalar beta. When BETA is
338 supplied as zero then C need not be set on input.
339 Unchanged on exit.
340
341 C - REAL array of DIMENSION ( LDC, n ).
342 Before entry, the leading m by n part of the array C must
343 contain the matrix C, except when beta is zero, in which
344 case C need not be set on entry.
345 On exit, the array C is overwritten by the m by n matrix
346 ( alpha*op( A )*op( B ) + beta*C ).
347
348 LDC - INTEGER.
349 On entry, LDC specifies the first dimension of C as declared
350 in the calling (sub) program. LDC must be at least
351 max( 1, m ).
352 Unchanged on exit.
353
354
355 Level 3 Blas routine.
356
357 -- Written on 8-February-1989.
358 Jack Dongarra, Argonne National Laboratory.
359 Iain Duff, AERE Harwell.
360 Jeremy Du Croz, Numerical Algorithms Group Ltd.
361 Sven Hammarling, Numerical Algorithms Group Ltd.
362
363
364 Set NOTA and NOTB as true if A and B respectively are not
365 transposed and set NROWA, NCOLA and NROWB as the number of rows
366 and columns of A and the number of rows of B respectively.
367*/
368
369 /* Parameter adjustments */
370 a_dim1 = *lda;
371 a_offset = 1 + a_dim1;
372 a -= a_offset;
373 b_dim1 = *ldb;
374 b_offset = 1 + b_dim1;
375 b -= b_offset;
376 c_dim1 = *ldc;
377 c_offset = 1 + c_dim1;
378 c__ -= c_offset;
379
380 /* Function Body */
381 nota = lsame_(transa, "N");
382 notb = lsame_(transb, "N");
383 if (nota) {
384 nrowa = *m;
385 ncola = *k;
386 } else {
387 nrowa = *k;
388 ncola = *m;
389 }
390 if (notb) {
391 nrowb = *k;
392 } else {
393 nrowb = *n;
394 }
395
396/* Test the input parameters. */
397
398 info = 0;
399 if (! nota && ! lsame_(transa, "C") && ! lsame_(
400 transa, "T")) {
401 info = 1;
402 } else if (! notb && ! lsame_(transb, "C") && !
403 lsame_(transb, "T")) {
404 info = 2;
405 } else if (*m < 0) {
406 info = 3;
407 } else if (*n < 0) {
408 info = 4;
409 } else if (*k < 0) {
410 info = 5;
411 } else if (*lda < max(1,nrowa)) {
412 info = 8;
413 } else if (*ldb < max(1,nrowb)) {
414 info = 10;
415 } else if (*ldc < max(1,*m)) {
416 info = 13;
417 }
418 if (info != 0) {
419 xerbla_("SGEMM ", &info);
420 return 0;
421 }
422
423/* Quick return if possible. */
424
425 if (*m == 0 || *n == 0 || (*alpha == 0.f || *k == 0) && *beta == 1.f) {
426 return 0;
427 }
428
429/* And if alpha.eq.zero. */
430
431 if (*alpha == 0.f) {
432 if (*beta == 0.f) {
433 i__1 = *n;
434 for (j = 1; j <= i__1; ++j) {
435 i__2 = *m;
436 for (i__ = 1; i__ <= i__2; ++i__) {
437 c__[i__ + j * c_dim1] = 0.f;
438/* L10: */
439 }
440/* L20: */
441 }
442 } else {
443 i__1 = *n;
444 for (j = 1; j <= i__1; ++j) {
445 i__2 = *m;
446 for (i__ = 1; i__ <= i__2; ++i__) {
447 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
448/* L30: */
449 }
450/* L40: */
451 }
452 }
453 return 0;
454 }
455
456/* Start the operations. */
457
458 if (notb) {
459 if (nota) {
460
461/* Form C := alpha*A*B + beta*C. */
462
463 i__1 = *n;
464 for (j = 1; j <= i__1; ++j) {
465 if (*beta == 0.f) {
466 i__2 = *m;
467 for (i__ = 1; i__ <= i__2; ++i__) {
468 c__[i__ + j * c_dim1] = 0.f;
469/* L50: */
470 }
471 } else if (*beta != 1.f) {
472 i__2 = *m;
473 for (i__ = 1; i__ <= i__2; ++i__) {
474 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
475/* L60: */
476 }
477 }
478 i__2 = *k;
479 for (l = 1; l <= i__2; ++l) {
480 if (b[l + j * b_dim1] != 0.f) {
481 temp = *alpha * b[l + j * b_dim1];
482 i__3 = *m;
483 for (i__ = 1; i__ <= i__3; ++i__) {
484 c__[i__ + j * c_dim1] += temp * a[i__ + l *
485 a_dim1];
486/* L70: */
487 }
488 }
489/* L80: */
490 }
491/* L90: */
492 }
493 } else {
494
495/* Form C := alpha*A'*B + beta*C */
496
497 i__1 = *n;
498 for (j = 1; j <= i__1; ++j) {
499 i__2 = *m;
500 for (i__ = 1; i__ <= i__2; ++i__) {
501 temp = 0.f;
502 i__3 = *k;
503 for (l = 1; l <= i__3; ++l) {
504 temp += a[l + i__ * a_dim1] * b[l + j * b_dim1];
505/* L100: */
506 }
507 if (*beta == 0.f) {
508 c__[i__ + j * c_dim1] = *alpha * temp;
509 } else {
510 c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
511 i__ + j * c_dim1];
512 }
513/* L110: */
514 }
515/* L120: */
516 }
517 }
518 } else {
519 if (nota) {
520
521/* Form C := alpha*A*B' + beta*C */
522
523 i__1 = *n;
524 for (j = 1; j <= i__1; ++j) {
525 if (*beta == 0.f) {
526 i__2 = *m;
527 for (i__ = 1; i__ <= i__2; ++i__) {
528 c__[i__ + j * c_dim1] = 0.f;
529/* L130: */
530 }
531 } else if (*beta != 1.f) {
532 i__2 = *m;
533 for (i__ = 1; i__ <= i__2; ++i__) {
534 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
535/* L140: */
536 }
537 }
538 i__2 = *k;
539 for (l = 1; l <= i__2; ++l) {
540 if (b[j + l * b_dim1] != 0.f) {
541 temp = *alpha * b[j + l * b_dim1];
542 i__3 = *m;
543 for (i__ = 1; i__ <= i__3; ++i__) {
544 c__[i__ + j * c_dim1] += temp * a[i__ + l *
545 a_dim1];
546/* L150: */
547 }
548 }
549/* L160: */
550 }
551/* L170: */
552 }
553 } else {
554
555/* Form C := alpha*A'*B' + beta*C */
556
557 i__1 = *n;
558 for (j = 1; j <= i__1; ++j) {
559 i__2 = *m;
560 for (i__ = 1; i__ <= i__2; ++i__) {
561 temp = 0.f;
562 i__3 = *k;
563 for (l = 1; l <= i__3; ++l) {
564 temp += a[l + i__ * a_dim1] * b[j + l * b_dim1];
565/* L180: */
566 }
567 if (*beta == 0.f) {
568 c__[i__ + j * c_dim1] = *alpha * temp;
569 } else {
570 c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
571 i__ + j * c_dim1];
572 }
573/* L190: */
574 }
575/* L200: */
576 }
577 }
578 }
579
580 return 0;
581
582/* End of SGEMM . */
583
584} /* sgemm_ */
585
586/* Subroutine */ int sgemv_(char *trans, integer *m, integer *n, real *alpha,
587 real *a, integer *lda, real *x, integer *incx, real *beta, real *y,
588 integer *incy)
589{
590 /* System generated locals */
591 integer a_dim1, a_offset, i__1, i__2;
592
593 /* Local variables */
594 static integer i__, j, ix, iy, jx, jy, kx, ky, info;
595 static real temp;
596 static integer lenx, leny;
597 extern logical lsame_(char *, char *);
598 extern /* Subroutine */ int xerbla_(char *, integer *);
599
600
601/*
602 Purpose
603 =======
604
605 SGEMV performs one of the matrix-vector operations
606
607 y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y,
608
609 where alpha and beta are scalars, x and y are vectors and A is an
610 m by n matrix.
611
612 Parameters
613 ==========
614
615 TRANS - CHARACTER*1.
616 On entry, TRANS specifies the operation to be performed as
617 follows:
618
619 TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
620
621 TRANS = 'T' or 't' y := alpha*A'*x + beta*y.
622
623 TRANS = 'C' or 'c' y := alpha*A'*x + beta*y.
624
625 Unchanged on exit.
626
627 M - INTEGER.
628 On entry, M specifies the number of rows of the matrix A.
629 M must be at least zero.
630 Unchanged on exit.
631
632 N - INTEGER.
633 On entry, N specifies the number of columns of the matrix A.
634 N must be at least zero.
635 Unchanged on exit.
636
637 ALPHA - REAL .
638 On entry, ALPHA specifies the scalar alpha.
639 Unchanged on exit.
640
641 A - REAL array of DIMENSION ( LDA, n ).
642 Before entry, the leading m by n part of the array A must
643 contain the matrix of coefficients.
644 Unchanged on exit.
645
646 LDA - INTEGER.
647 On entry, LDA specifies the first dimension of A as declared
648 in the calling (sub) program. LDA must be at least
649 max( 1, m ).
650 Unchanged on exit.
651
652 X - REAL array of DIMENSION at least
653 ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
654 and at least
655 ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
656 Before entry, the incremented array X must contain the
657 vector x.
658 Unchanged on exit.
659
660 INCX - INTEGER.
661 On entry, INCX specifies the increment for the elements of
662 X. INCX must not be zero.
663 Unchanged on exit.
664
665 BETA - REAL .
666 On entry, BETA specifies the scalar beta. When BETA is
667 supplied as zero then Y need not be set on input.
668 Unchanged on exit.
669
670 Y - REAL array of DIMENSION at least
671 ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
672 and at least
673 ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
674 Before entry with BETA non-zero, the incremented array Y
675 must contain the vector y. On exit, Y is overwritten by the
676 updated vector y.
677
678 INCY - INTEGER.
679 On entry, INCY specifies the increment for the elements of
680 Y. INCY must not be zero.
681 Unchanged on exit.
682
683
684 Level 2 Blas routine.
685
686 -- Written on 22-October-1986.
687 Jack Dongarra, Argonne National Lab.
688 Jeremy Du Croz, Nag Central Office.
689 Sven Hammarling, Nag Central Office.
690 Richard Hanson, Sandia National Labs.
691
692
693 Test the input parameters.
694*/
695
696 /* Parameter adjustments */
697 a_dim1 = *lda;
698 a_offset = 1 + a_dim1;
699 a -= a_offset;
700 --x;
701 --y;
702
703 /* Function Body */
704 info = 0;
705 if (! lsame_(trans, "N") && ! lsame_(trans, "T") && ! lsame_(trans, "C")
706 ) {
707 info = 1;
708 } else if (*m < 0) {
709 info = 2;
710 } else if (*n < 0) {
711 info = 3;
712 } else if (*lda < max(1,*m)) {
713 info = 6;
714 } else if (*incx == 0) {
715 info = 8;
716 } else if (*incy == 0) {
717 info = 11;
718 }
719 if (info != 0) {
720 xerbla_("SGEMV ", &info);
721 return 0;
722 }
723
724/* Quick return if possible. */
725
726 if (*m == 0 || *n == 0 || *alpha == 0.f && *beta == 1.f) {
727 return 0;
728 }
729
730/*
731 Set LENX and LENY, the lengths of the vectors x and y, and set
732 up the start points in X and Y.
733*/
734
735 if (lsame_(trans, "N")) {
736 lenx = *n;
737 leny = *m;
738 } else {
739 lenx = *m;
740 leny = *n;
741 }
742 if (*incx > 0) {
743 kx = 1;
744 } else {
745 kx = 1 - (lenx - 1) * *incx;
746 }
747 if (*incy > 0) {
748 ky = 1;
749 } else {
750 ky = 1 - (leny - 1) * *incy;
751 }
752
753/*
754 Start the operations. In this version the elements of A are
755 accessed sequentially with one pass through A.
756
757 First form y := beta*y.
758*/
759
760 if (*beta != 1.f) {
761 if (*incy == 1) {
762 if (*beta == 0.f) {
763 i__1 = leny;
764 for (i__ = 1; i__ <= i__1; ++i__) {
765 y[i__] = 0.f;
766/* L10: */
767 }
768 } else {
769 i__1 = leny;
770 for (i__ = 1; i__ <= i__1; ++i__) {
771 y[i__] = *beta * y[i__];
772/* L20: */
773 }
774 }
775 } else {
776 iy = ky;
777 if (*beta == 0.f) {
778 i__1 = leny;
779 for (i__ = 1; i__ <= i__1; ++i__) {
780 y[iy] = 0.f;
781 iy += *incy;
782/* L30: */
783 }
784 } else {
785 i__1 = leny;
786 for (i__ = 1; i__ <= i__1; ++i__) {
787 y[iy] = *beta * y[iy];
788 iy += *incy;
789/* L40: */
790 }
791 }
792 }
793 }
794 if (*alpha == 0.f) {
795 return 0;
796 }
797 if (lsame_(trans, "N")) {
798
799/* Form y := alpha*A*x + y. */
800
801 jx = kx;
802 if (*incy == 1) {
803 i__1 = *n;
804 for (j = 1; j <= i__1; ++j) {
805 if (x[jx] != 0.f) {
806 temp = *alpha * x[jx];
807 i__2 = *m;
808 for (i__ = 1; i__ <= i__2; ++i__) {
809 y[i__] += temp * a[i__ + j * a_dim1];
810/* L50: */
811 }
812 }
813 jx += *incx;
814/* L60: */
815 }
816 } else {
817 i__1 = *n;
818 for (j = 1; j <= i__1; ++j) {
819 if (x[jx] != 0.f) {
820 temp = *alpha * x[jx];
821 iy = ky;
822 i__2 = *m;
823 for (i__ = 1; i__ <= i__2; ++i__) {
824 y[iy] += temp * a[i__ + j * a_dim1];
825 iy += *incy;
826/* L70: */
827 }
828 }
829 jx += *incx;
830/* L80: */
831 }
832 }
833 } else {
834
835/* Form y := alpha*A'*x + y. */
836
837 jy = ky;
838 if (*incx == 1) {
839 i__1 = *n;
840 for (j = 1; j <= i__1; ++j) {
841 temp = 0.f;
842 i__2 = *m;
843 for (i__ = 1; i__ <= i__2; ++i__) {
844 temp += a[i__ + j * a_dim1] * x[i__];
845/* L90: */
846 }
847 y[jy] += *alpha * temp;
848 jy += *incy;
849/* L100: */
850 }
851 } else {
852 i__1 = *n;
853 for (j = 1; j <= i__1; ++j) {
854 temp = 0.f;
855 ix = kx;
856 i__2 = *m;
857 for (i__ = 1; i__ <= i__2; ++i__) {
858 temp += a[i__ + j * a_dim1] * x[ix];
859 ix += *incx;
860/* L110: */
861 }
862 y[jy] += *alpha * temp;
863 jy += *incy;
864/* L120: */
865 }
866 }
867 }
868
869 return 0;
870
871/* End of SGEMV . */
872
873} /* sgemv_ */
874
875/* Subroutine */ int sscal_(integer *n, real *sa, real *sx, integer *incx)
876{
877 /* System generated locals */
878 integer i__1, i__2;
879
880 /* Local variables */
881 static integer i__, m, mp1, nincx;
882
883
884/*
885 scales a vector by a constant.
886 uses unrolled loops for increment equal to 1.
887 jack dongarra, linpack, 3/11/78.
888 modified 3/93 to return if incx .le. 0.
889 modified 12/3/93, array(1) declarations changed to array(*)
890*/
891
892
893 /* Parameter adjustments */
894 --sx;
895
896 /* Function Body */
897 if (*n <= 0 || *incx <= 0) {
898 return 0;
899 }
900 if (*incx == 1) {
901 goto L20;
902 }
903
904/* code for increment not equal to 1 */
905
906 nincx = *n * *incx;
907 i__1 = nincx;
908 i__2 = *incx;
909 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
910 sx[i__] = *sa * sx[i__];
911/* L10: */
912 }
913 return 0;
914
915/*
916 code for increment equal to 1
917
918
919 clean-up loop
920*/
921
922L20:
923 m = *n % 5;
924 if (m == 0) {
925 goto L40;
926 }
927 i__2 = m;
928 for (i__ = 1; i__ <= i__2; ++i__) {
929 sx[i__] = *sa * sx[i__];
930/* L30: */
931 }
932 if (*n < 5) {
933 return 0;
934 }
935L40:
936 mp1 = m + 1;
937 i__2 = *n;
938 for (i__ = mp1; i__ <= i__2; i__ += 5) {
939 sx[i__] = *sa * sx[i__];
940 sx[i__ + 1] = *sa * sx[i__ + 1];
941 sx[i__ + 2] = *sa * sx[i__ + 2];
942 sx[i__ + 3] = *sa * sx[i__ + 3];
943 sx[i__ + 4] = *sa * sx[i__ + 4];
944/* L50: */
945 }
946 return 0;
947} /* sscal_ */
948
949/* Subroutine */ int ssymm_(char *side, char *uplo, integer *m, integer *n,
950 real *alpha, real *a, integer *lda, real *b, integer *ldb, real *beta,
951 real *c__, integer *ldc)
952{
953 /* System generated locals */
954 integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
955 i__3;
956
957 /* Local variables */
958 static integer i__, j, k, info;
959 static real temp1, temp2;
960 extern logical lsame_(char *, char *);
961 static integer nrowa;
962 static logical upper;
963 extern /* Subroutine */ int xerbla_(char *, integer *);
964
965
966/*
967 Purpose
968 =======
969
970 SSYMM performs one of the matrix-matrix operations
971
972 C := alpha*A*B + beta*C,
973
974 or
975
976 C := alpha*B*A + beta*C,
977
978 where alpha and beta are scalars, A is a symmetric matrix and B and
979 C are m by n matrices.
980
981 Parameters
982 ==========
983
984 SIDE - CHARACTER*1.
985 On entry, SIDE specifies whether the symmetric matrix A
986 appears on the left or right in the operation as follows:
987
988 SIDE = 'L' or 'l' C := alpha*A*B + beta*C,
989
990 SIDE = 'R' or 'r' C := alpha*B*A + beta*C,
991
992 Unchanged on exit.
993
994 UPLO - CHARACTER*1.
995 On entry, UPLO specifies whether the upper or lower
996 triangular part of the symmetric matrix A is to be
997 referenced as follows:
998
999 UPLO = 'U' or 'u' Only the upper triangular part of the
1000 symmetric matrix is to be referenced.
1001
1002 UPLO = 'L' or 'l' Only the lower triangular part of the
1003 symmetric matrix is to be referenced.
1004
1005 Unchanged on exit.
1006
1007 M - INTEGER.
1008 On entry, M specifies the number of rows of the matrix C.
1009 M must be at least zero.
1010 Unchanged on exit.
1011
1012 N - INTEGER.
1013 On entry, N specifies the number of columns of the matrix C.
1014 N must be at least zero.
1015 Unchanged on exit.
1016
1017 ALPHA - REAL .
1018 On entry, ALPHA specifies the scalar alpha.
1019 Unchanged on exit.
1020
1021 A - REAL array of DIMENSION ( LDA, ka ), where ka is
1022 m when SIDE = 'L' or 'l' and is n otherwise.
1023 Before entry with SIDE = 'L' or 'l', the m by m part of
1024 the array A must contain the symmetric matrix, such that
1025 when UPLO = 'U' or 'u', the leading m by m upper triangular
1026 part of the array A must contain the upper triangular part
1027 of the symmetric matrix and the strictly lower triangular
1028 part of A is not referenced, and when UPLO = 'L' or 'l',
1029 the leading m by m lower triangular part of the array A
1030 must contain the lower triangular part of the symmetric
1031 matrix and the strictly upper triangular part of A is not
1032 referenced.
1033 Before entry with SIDE = 'R' or 'r', the n by n part of
1034 the array A must contain the symmetric matrix, such that
1035 when UPLO = 'U' or 'u', the leading n by n upper triangular
1036 part of the array A must contain the upper triangular part
1037 of the symmetric matrix and the strictly lower triangular
1038 part of A is not referenced, and when UPLO = 'L' or 'l',
1039 the leading n by n lower triangular part of the array A
1040 must contain the lower triangular part of the symmetric
1041 matrix and the strictly upper triangular part of A is not
1042 referenced.
1043 Unchanged on exit.
1044
1045 LDA - INTEGER.
1046 On entry, LDA specifies the first dimension of A as declared
1047 in the calling (sub) program. When SIDE = 'L' or 'l' then
1048 LDA must be at least max( 1, m ), otherwise LDA must be at
1049 least max( 1, n ).
1050 Unchanged on exit.
1051
1052 B - REAL array of DIMENSION ( LDB, n ).
1053 Before entry, the leading m by n part of the array B must
1054 contain the matrix B.
1055 Unchanged on exit.
1056
1057 LDB - INTEGER.
1058 On entry, LDB specifies the first dimension of B as declared
1059 in the calling (sub) program. LDB must be at least
1060 max( 1, m ).
1061 Unchanged on exit.
1062
1063 BETA - REAL .
1064 On entry, BETA specifies the scalar beta. When BETA is
1065 supplied as zero then C need not be set on input.
1066 Unchanged on exit.
1067
1068 C - REAL array of DIMENSION ( LDC, n ).
1069 Before entry, the leading m by n part of the array C must
1070 contain the matrix C, except when beta is zero, in which
1071 case C need not be set on entry.
1072 On exit, the array C is overwritten by the m by n updated
1073 matrix.
1074
1075 LDC - INTEGER.
1076 On entry, LDC specifies the first dimension of C as declared
1077 in the calling (sub) program. LDC must be at least
1078 max( 1, m ).
1079 Unchanged on exit.
1080
1081
1082 Level 3 Blas routine.
1083
1084 -- Written on 8-February-1989.
1085 Jack Dongarra, Argonne National Laboratory.
1086 Iain Duff, AERE Harwell.
1087 Jeremy Du Croz, Numerical Algorithms Group Ltd.
1088 Sven Hammarling, Numerical Algorithms Group Ltd.
1089
1090
1091 Set NROWA as the number of rows of A.
1092*/
1093
1094 /* Parameter adjustments */
1095 a_dim1 = *lda;
1096 a_offset = 1 + a_dim1;
1097 a -= a_offset;
1098 b_dim1 = *ldb;
1099 b_offset = 1 + b_dim1;
1100 b -= b_offset;
1101 c_dim1 = *ldc;
1102 c_offset = 1 + c_dim1;
1103 c__ -= c_offset;
1104
1105 /* Function Body */
1106 if (lsame_(side, "L")) {
1107 nrowa = *m;
1108 } else {
1109 nrowa = *n;
1110 }
1111 upper = lsame_(uplo, "U");
1112
1113/* Test the input parameters. */
1114
1115 info = 0;
1116 if (! lsame_(side, "L") && ! lsame_(side, "R")) {
1117 info = 1;
1118 } else if (! upper && ! lsame_(uplo, "L")) {
1119 info = 2;
1120 } else if (*m < 0) {
1121 info = 3;
1122 } else if (*n < 0) {
1123 info = 4;
1124 } else if (*lda < max(1,nrowa)) {
1125 info = 7;
1126 } else if (*ldb < max(1,*m)) {
1127 info = 9;
1128 } else if (*ldc < max(1,*m)) {
1129 info = 12;
1130 }
1131 if (info != 0) {
1132 xerbla_("SSYMM ", &info);
1133 return 0;
1134 }
1135
1136/* Quick return if possible. */
1137
1138 if (*m == 0 || *n == 0 || *alpha == 0.f && *beta == 1.f) {
1139 return 0;
1140 }
1141
1142/* And when alpha.eq.zero. */
1143
1144 if (*alpha == 0.f) {
1145 if (*beta == 0.f) {
1146 i__1 = *n;
1147 for (j = 1; j <= i__1; ++j) {
1148 i__2 = *m;
1149 for (i__ = 1; i__ <= i__2; ++i__) {
1150 c__[i__ + j * c_dim1] = 0.f;
1151/* L10: */
1152 }
1153/* L20: */
1154 }
1155 } else {
1156 i__1 = *n;
1157 for (j = 1; j <= i__1; ++j) {
1158 i__2 = *m;
1159 for (i__ = 1; i__ <= i__2; ++i__) {
1160 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
1161/* L30: */
1162 }
1163/* L40: */
1164 }
1165 }
1166 return 0;
1167 }
1168
1169/* Start the operations. */
1170
1171 if (lsame_(side, "L")) {
1172
1173/* Form C := alpha*A*B + beta*C. */
1174
1175 if (upper) {
1176 i__1 = *n;
1177 for (j = 1; j <= i__1; ++j) {
1178 i__2 = *m;
1179 for (i__ = 1; i__ <= i__2; ++i__) {
1180 temp1 = *alpha * b[i__ + j * b_dim1];
1181 temp2 = 0.f;
1182 i__3 = i__ - 1;
1183 for (k = 1; k <= i__3; ++k) {
1184 c__[k + j * c_dim1] += temp1 * a[k + i__ * a_dim1];
1185 temp2 += b[k + j * b_dim1] * a[k + i__ * a_dim1];
1186/* L50: */
1187 }
1188 if (*beta == 0.f) {
1189 c__[i__ + j * c_dim1] = temp1 * a[i__ + i__ * a_dim1]
1190 + *alpha * temp2;
1191 } else {
1192 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1]
1193 + temp1 * a[i__ + i__ * a_dim1] + *alpha *
1194 temp2;
1195 }
1196/* L60: */
1197 }
1198/* L70: */
1199 }
1200 } else {
1201 i__1 = *n;
1202 for (j = 1; j <= i__1; ++j) {
1203 for (i__ = *m; i__ >= 1; --i__) {
1204 temp1 = *alpha * b[i__ + j * b_dim1];
1205 temp2 = 0.f;
1206 i__2 = *m;
1207 for (k = i__ + 1; k <= i__2; ++k) {
1208 c__[k + j * c_dim1] += temp1 * a[k + i__ * a_dim1];
1209 temp2 += b[k + j * b_dim1] * a[k + i__ * a_dim1];
1210/* L80: */
1211 }
1212 if (*beta == 0.f) {
1213 c__[i__ + j * c_dim1] = temp1 * a[i__ + i__ * a_dim1]
1214 + *alpha * temp2;
1215 } else {
1216 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1]
1217 + temp1 * a[i__ + i__ * a_dim1] + *alpha *
1218 temp2;
1219 }
1220/* L90: */
1221 }
1222/* L100: */
1223 }
1224 }
1225 } else {
1226
1227/* Form C := alpha*B*A + beta*C. */
1228
1229 i__1 = *n;
1230 for (j = 1; j <= i__1; ++j) {
1231 temp1 = *alpha * a[j + j * a_dim1];
1232 if (*beta == 0.f) {
1233 i__2 = *m;
1234 for (i__ = 1; i__ <= i__2; ++i__) {
1235 c__[i__ + j * c_dim1] = temp1 * b[i__ + j * b_dim1];
1236/* L110: */
1237 }
1238 } else {
1239 i__2 = *m;
1240 for (i__ = 1; i__ <= i__2; ++i__) {
1241 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1] +
1242 temp1 * b[i__ + j * b_dim1];
1243/* L120: */
1244 }
1245 }
1246 i__2 = j - 1;
1247 for (k = 1; k <= i__2; ++k) {
1248 if (upper) {
1249 temp1 = *alpha * a[k + j * a_dim1];
1250 } else {
1251 temp1 = *alpha * a[j + k * a_dim1];
1252 }
1253 i__3 = *m;
1254 for (i__ = 1; i__ <= i__3; ++i__) {
1255 c__[i__ + j * c_dim1] += temp1 * b[i__ + k * b_dim1];
1256/* L130: */
1257 }
1258/* L140: */
1259 }
1260 i__2 = *n;
1261 for (k = j + 1; k <= i__2; ++k) {
1262 if (upper) {
1263 temp1 = *alpha * a[j + k * a_dim1];
1264 } else {
1265 temp1 = *alpha * a[k + j * a_dim1];
1266 }
1267 i__3 = *m;
1268 for (i__ = 1; i__ <= i__3; ++i__) {
1269 c__[i__ + j * c_dim1] += temp1 * b[i__ + k * b_dim1];
1270/* L150: */
1271 }
1272/* L160: */
1273 }
1274/* L170: */
1275 }
1276 }
1277
1278 return 0;
1279
1280/* End of SSYMM . */
1281
1282} /* ssymm_ */
1283
1284/* Subroutine */ int ssyrk_(char *uplo, char *trans, integer *n, integer *k,
1285 real *alpha, real *a, integer *lda, real *beta, real *c__, integer *
1286 ldc)
1287{
1288 /* System generated locals */
1289 integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2, i__3;
1290
1291 /* Local variables */
1292 static integer i__, j, l, info;
1293 static real temp;
1294 extern logical lsame_(char *, char *);
1295 static integer nrowa;
1296 static logical upper;
1297 extern /* Subroutine */ int xerbla_(char *, integer *);
1298
1299
1300/*
1301 Purpose
1302 =======
1303
1304 SSYRK performs one of the symmetric rank k operations
1305
1306 C := alpha*A*A' + beta*C,
1307
1308 or
1309
1310 C := alpha*A'*A + beta*C,
1311
1312 where alpha and beta are scalars, C is an n by n symmetric matrix
1313 and A is an n by k matrix in the first case and a k by n matrix
1314 in the second case.
1315
1316 Parameters
1317 ==========
1318
1319 UPLO - CHARACTER*1.
1320 On entry, UPLO specifies whether the upper or lower
1321 triangular part of the array C is to be referenced as
1322 follows:
1323
1324 UPLO = 'U' or 'u' Only the upper triangular part of C
1325 is to be referenced.
1326
1327 UPLO = 'L' or 'l' Only the lower triangular part of C
1328 is to be referenced.
1329
1330 Unchanged on exit.
1331
1332 TRANS - CHARACTER*1.
1333 On entry, TRANS specifies the operation to be performed as
1334 follows:
1335
1336 TRANS = 'N' or 'n' C := alpha*A*A' + beta*C.
1337
1338 TRANS = 'T' or 't' C := alpha*A'*A + beta*C.
1339
1340 TRANS = 'C' or 'c' C := alpha*A'*A + beta*C.
1341
1342 Unchanged on exit.
1343
1344 N - INTEGER.
1345 On entry, N specifies the order of the matrix C. N must be
1346 at least zero.
1347 Unchanged on exit.
1348
1349 K - INTEGER.
1350 On entry with TRANS = 'N' or 'n', K specifies the number
1351 of columns of the matrix A, and on entry with
1352 TRANS = 'T' or 't' or 'C' or 'c', K specifies the number
1353 of rows of the matrix A. K must be at least zero.
1354 Unchanged on exit.
1355
1356 ALPHA - REAL .
1357 On entry, ALPHA specifies the scalar alpha.
1358 Unchanged on exit.
1359
1360 A - REAL array of DIMENSION ( LDA, ka ), where ka is
1361 k when TRANS = 'N' or 'n', and is n otherwise.
1362 Before entry with TRANS = 'N' or 'n', the leading n by k
1363 part of the array A must contain the matrix A, otherwise
1364 the leading k by n part of the array A must contain the
1365 matrix A.
1366 Unchanged on exit.
1367
1368 LDA - INTEGER.
1369 On entry, LDA specifies the first dimension of A as declared
1370 in the calling (sub) program. When TRANS = 'N' or 'n'
1371 then LDA must be at least max( 1, n ), otherwise LDA must
1372 be at least max( 1, k ).
1373 Unchanged on exit.
1374
1375 BETA - REAL .
1376 On entry, BETA specifies the scalar beta.
1377 Unchanged on exit.
1378
1379 C - REAL array of DIMENSION ( LDC, n ).
1380 Before entry with UPLO = 'U' or 'u', the leading n by n
1381 upper triangular part of the array C must contain the upper
1382 triangular part of the symmetric matrix and the strictly
1383 lower triangular part of C is not referenced. On exit, the
1384 upper triangular part of the array C is overwritten by the
1385 upper triangular part of the updated matrix.
1386 Before entry with UPLO = 'L' or 'l', the leading n by n
1387 lower triangular part of the array C must contain the lower
1388 triangular part of the symmetric matrix and the strictly
1389 upper triangular part of C is not referenced. On exit, the
1390 lower triangular part of the array C is overwritten by the
1391 lower triangular part of the updated matrix.
1392
1393 LDC - INTEGER.
1394 On entry, LDC specifies the first dimension of C as declared
1395 in the calling (sub) program. LDC must be at least
1396 max( 1, n ).
1397 Unchanged on exit.
1398
1399
1400 Level 3 Blas routine.
1401
1402 -- Written on 8-February-1989.
1403 Jack Dongarra, Argonne National Laboratory.
1404 Iain Duff, AERE Harwell.
1405 Jeremy Du Croz, Numerical Algorithms Group Ltd.
1406 Sven Hammarling, Numerical Algorithms Group Ltd.
1407
1408
1409 Test the input parameters.
1410*/
1411
1412 /* Parameter adjustments */
1413 a_dim1 = *lda;
1414 a_offset = 1 + a_dim1;
1415 a -= a_offset;
1416 c_dim1 = *ldc;
1417 c_offset = 1 + c_dim1;
1418 c__ -= c_offset;
1419
1420 /* Function Body */
1421 if (lsame_(trans, "N")) {
1422 nrowa = *n;
1423 } else {
1424 nrowa = *k;
1425 }
1426 upper = lsame_(uplo, "U");
1427
1428 info = 0;
1429 if (! upper && ! lsame_(uplo, "L")) {
1430 info = 1;
1431 } else if (! lsame_(trans, "N") && ! lsame_(trans,
1432 "T") && ! lsame_(trans, "C")) {
1433 info = 2;
1434 } else if (*n < 0) {
1435 info = 3;
1436 } else if (*k < 0) {
1437 info = 4;
1438 } else if (*lda < max(1,nrowa)) {
1439 info = 7;
1440 } else if (*ldc < max(1,*n)) {
1441 info = 10;
1442 }
1443 if (info != 0) {
1444 xerbla_("SSYRK ", &info);
1445 return 0;
1446 }
1447
1448/* Quick return if possible. */
1449
1450 if (*n == 0 || (*alpha == 0.f || *k == 0) && *beta == 1.f) {
1451 return 0;
1452 }
1453
1454/* And when alpha.eq.zero. */
1455
1456 if (*alpha == 0.f) {
1457 if (upper) {
1458 if (*beta == 0.f) {
1459 i__1 = *n;
1460 for (j = 1; j <= i__1; ++j) {
1461 i__2 = j;
1462 for (i__ = 1; i__ <= i__2; ++i__) {
1463 c__[i__ + j * c_dim1] = 0.f;
1464/* L10: */
1465 }
1466/* L20: */
1467 }
1468 } else {
1469 i__1 = *n;
1470 for (j = 1; j <= i__1; ++j) {
1471 i__2 = j;
1472 for (i__ = 1; i__ <= i__2; ++i__) {
1473 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
1474/* L30: */
1475 }
1476/* L40: */
1477 }
1478 }
1479 } else {
1480 if (*beta == 0.f) {
1481 i__1 = *n;
1482 for (j = 1; j <= i__1; ++j) {
1483 i__2 = *n;
1484 for (i__ = j; i__ <= i__2; ++i__) {
1485 c__[i__ + j * c_dim1] = 0.f;
1486/* L50: */
1487 }
1488/* L60: */
1489 }
1490 } else {
1491 i__1 = *n;
1492 for (j = 1; j <= i__1; ++j) {
1493 i__2 = *n;
1494 for (i__ = j; i__ <= i__2; ++i__) {
1495 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
1496/* L70: */
1497 }
1498/* L80: */
1499 }
1500 }
1501 }
1502 return 0;
1503 }
1504
1505/* Start the operations. */
1506
1507 if (lsame_(trans, "N")) {
1508
1509/* Form C := alpha*A*A' + beta*C. */
1510
1511 if (upper) {
1512 i__1 = *n;
1513 for (j = 1; j <= i__1; ++j) {
1514 if (*beta == 0.f) {
1515 i__2 = j;
1516 for (i__ = 1; i__ <= i__2; ++i__) {
1517 c__[i__ + j * c_dim1] = 0.f;
1518/* L90: */
1519 }
1520 } else if (*beta != 1.f) {
1521 i__2 = j;
1522 for (i__ = 1; i__ <= i__2; ++i__) {
1523 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
1524/* L100: */
1525 }
1526 }
1527 i__2 = *k;
1528 for (l = 1; l <= i__2; ++l) {
1529 if (a[j + l * a_dim1] != 0.f) {
1530 temp = *alpha * a[j + l * a_dim1];
1531 i__3 = j;
1532 for (i__ = 1; i__ <= i__3; ++i__) {
1533 c__[i__ + j * c_dim1] += temp * a[i__ + l *
1534 a_dim1];
1535/* L110: */
1536 }
1537 }
1538/* L120: */
1539 }
1540/* L130: */
1541 }
1542 } else {
1543 i__1 = *n;
1544 for (j = 1; j <= i__1; ++j) {
1545 if (*beta == 0.f) {
1546 i__2 = *n;
1547 for (i__ = j; i__ <= i__2; ++i__) {
1548 c__[i__ + j * c_dim1] = 0.f;
1549/* L140: */
1550 }
1551 } else if (*beta != 1.f) {
1552 i__2 = *n;
1553 for (i__ = j; i__ <= i__2; ++i__) {
1554 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
1555/* L150: */
1556 }
1557 }
1558 i__2 = *k;
1559 for (l = 1; l <= i__2; ++l) {
1560 if (a[j + l * a_dim1] != 0.f) {
1561 temp = *alpha * a[j + l * a_dim1];
1562 i__3 = *n;
1563 for (i__ = j; i__ <= i__3; ++i__) {
1564 c__[i__ + j * c_dim1] += temp * a[i__ + l *
1565 a_dim1];
1566/* L160: */
1567 }
1568 }
1569/* L170: */
1570 }
1571/* L180: */
1572 }
1573 }
1574 } else {
1575
1576/* Form C := alpha*A'*A + beta*C. */
1577
1578 if (upper) {
1579 i__1 = *n;
1580 for (j = 1; j <= i__1; ++j) {
1581 i__2 = j;
1582 for (i__ = 1; i__ <= i__2; ++i__) {
1583 temp = 0.f;
1584 i__3 = *k;
1585 for (l = 1; l <= i__3; ++l) {
1586 temp += a[l + i__ * a_dim1] * a[l + j * a_dim1];
1587/* L190: */
1588 }
1589 if (*beta == 0.f) {
1590 c__[i__ + j * c_dim1] = *alpha * temp;
1591 } else {
1592 c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
1593 i__ + j * c_dim1];
1594 }
1595/* L200: */
1596 }
1597/* L210: */
1598 }
1599 } else {
1600 i__1 = *n;
1601 for (j = 1; j <= i__1; ++j) {
1602 i__2 = *n;
1603 for (i__ = j; i__ <= i__2; ++i__) {
1604 temp = 0.f;
1605 i__3 = *k;
1606 for (l = 1; l <= i__3; ++l) {
1607 temp += a[l + i__ * a_dim1] * a[l + j * a_dim1];
1608/* L220: */
1609 }
1610 if (*beta == 0.f) {
1611 c__[i__ + j * c_dim1] = *alpha * temp;
1612 } else {
1613 c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
1614 i__ + j * c_dim1];
1615 }
1616/* L230: */
1617 }
1618/* L240: */
1619 }
1620 }
1621 }
1622
1623 return 0;
1624
1625/* End of SSYRK . */
1626
1627} /* ssyrk_ */
1628
1629/* Subroutine */ int strsm_(char *side, char *uplo, char *transa, char *diag,
1630 integer *m, integer *n, real *alpha, real *a, integer *lda, real *b,
1631 integer *ldb)
1632{
1633 /* System generated locals */
1634 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
1635
1636 /* Local variables */
1637 static integer i__, j, k, info;
1638 static real temp;
1639 static logical lside;
1640 extern logical lsame_(char *, char *);
1641 static integer nrowa;
1642 static logical upper;
1643 extern /* Subroutine */ int xerbla_(char *, integer *);
1644 static logical nounit;
1645
1646
1647/*
1648 Purpose
1649 =======
1650
1651 STRSM solves one of the matrix equations
1652
1653 op( A )*X = alpha*B, or X*op( A ) = alpha*B,
1654
1655 where alpha is a scalar, X and B are m by n matrices, A is a unit, or
1656 non-unit, upper or lower triangular matrix and op( A ) is one of
1657
1658 op( A ) = A or op( A ) = A'.
1659
1660 The matrix X is overwritten on B.
1661
1662 Parameters
1663 ==========
1664
1665 SIDE - CHARACTER*1.
1666 On entry, SIDE specifies whether op( A ) appears on the left
1667 or right of X as follows:
1668
1669 SIDE = 'L' or 'l' op( A )*X = alpha*B.
1670
1671 SIDE = 'R' or 'r' X*op( A ) = alpha*B.
1672
1673 Unchanged on exit.
1674
1675 UPLO - CHARACTER*1.
1676 On entry, UPLO specifies whether the matrix A is an upper or
1677 lower triangular matrix as follows:
1678
1679 UPLO = 'U' or 'u' A is an upper triangular matrix.
1680
1681 UPLO = 'L' or 'l' A is a lower triangular matrix.
1682
1683 Unchanged on exit.
1684
1685 TRANSA - CHARACTER*1.
1686 On entry, TRANSA specifies the form of op( A ) to be used in
1687 the matrix multiplication as follows:
1688
1689 TRANSA = 'N' or 'n' op( A ) = A.
1690
1691 TRANSA = 'T' or 't' op( A ) = A'.
1692
1693 TRANSA = 'C' or 'c' op( A ) = A'.
1694
1695 Unchanged on exit.
1696
1697 DIAG - CHARACTER*1.
1698 On entry, DIAG specifies whether or not A is unit triangular
1699 as follows:
1700
1701 DIAG = 'U' or 'u' A is assumed to be unit triangular.
1702
1703 DIAG = 'N' or 'n' A is not assumed to be unit
1704 triangular.
1705
1706 Unchanged on exit.
1707
1708 M - INTEGER.
1709 On entry, M specifies the number of rows of B. M must be at
1710 least zero.
1711 Unchanged on exit.
1712
1713 N - INTEGER.
1714 On entry, N specifies the number of columns of B. N must be
1715 at least zero.
1716 Unchanged on exit.
1717
1718 ALPHA - REAL .
1719 On entry, ALPHA specifies the scalar alpha. When alpha is
1720 zero then A is not referenced and B need not be set before
1721 entry.
1722 Unchanged on exit.
1723
1724 A - REAL array of DIMENSION ( LDA, k ), where k is m
1725 when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
1726 Before entry with UPLO = 'U' or 'u', the leading k by k
1727 upper triangular part of the array A must contain the upper
1728 triangular matrix and the strictly lower triangular part of
1729 A is not referenced.
1730 Before entry with UPLO = 'L' or 'l', the leading k by k
1731 lower triangular part of the array A must contain the lower
1732 triangular matrix and the strictly upper triangular part of
1733 A is not referenced.
1734 Note that when DIAG = 'U' or 'u', the diagonal elements of
1735 A are not referenced either, but are assumed to be unity.
1736 Unchanged on exit.
1737
1738 LDA - INTEGER.
1739 On entry, LDA specifies the first dimension of A as declared
1740 in the calling (sub) program. When SIDE = 'L' or 'l' then
1741 LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
1742 then LDA must be at least max( 1, n ).
1743 Unchanged on exit.
1744
1745 B - REAL array of DIMENSION ( LDB, n ).
1746 Before entry, the leading m by n part of the array B must
1747 contain the right-hand side matrix B, and on exit is
1748 overwritten by the solution matrix X.
1749
1750 LDB - INTEGER.
1751 On entry, LDB specifies the first dimension of B as declared
1752 in the calling (sub) program. LDB must be at least
1753 max( 1, m ).
1754 Unchanged on exit.
1755
1756
1757 Level 3 Blas routine.
1758
1759
1760 -- Written on 8-February-1989.
1761 Jack Dongarra, Argonne National Laboratory.
1762 Iain Duff, AERE Harwell.
1763 Jeremy Du Croz, Numerical Algorithms Group Ltd.
1764 Sven Hammarling, Numerical Algorithms Group Ltd.
1765
1766
1767 Test the input parameters.
1768*/
1769
1770 /* Parameter adjustments */
1771 a_dim1 = *lda;
1772 a_offset = 1 + a_dim1;
1773 a -= a_offset;
1774 b_dim1 = *ldb;
1775 b_offset = 1 + b_dim1;
1776 b -= b_offset;
1777
1778 /* Function Body */
1779 lside = lsame_(side, "L");
1780 if (lside) {
1781 nrowa = *m;
1782 } else {
1783 nrowa = *n;
1784 }
1785 nounit = lsame_(diag, "N");
1786 upper = lsame_(uplo, "U");
1787
1788 info = 0;
1789 if (! lside && ! lsame_(side, "R")) {
1790 info = 1;
1791 } else if (! upper && ! lsame_(uplo, "L")) {
1792 info = 2;
1793 } else if (! lsame_(transa, "N") && ! lsame_(transa,
1794 "T") && ! lsame_(transa, "C")) {
1795 info = 3;
1796 } else if (! lsame_(diag, "U") && ! lsame_(diag,
1797 "N")) {
1798 info = 4;
1799 } else if (*m < 0) {
1800 info = 5;
1801 } else if (*n < 0) {
1802 info = 6;
1803 } else if (*lda < max(1,nrowa)) {
1804 info = 9;
1805 } else if (*ldb < max(1,*m)) {
1806 info = 11;
1807 }
1808 if (info != 0) {
1809 xerbla_("STRSM ", &info);
1810 return 0;
1811 }
1812
1813/* Quick return if possible. */
1814
1815 if (*n == 0) {
1816 return 0;
1817 }
1818
1819/* And when alpha.eq.zero. */
1820
1821 if (*alpha == 0.f) {
1822 i__1 = *n;
1823 for (j = 1; j <= i__1; ++j) {
1824 i__2 = *m;
1825 for (i__ = 1; i__ <= i__2; ++i__) {
1826 b[i__ + j * b_dim1] = 0.f;
1827/* L10: */
1828 }
1829/* L20: */
1830 }
1831 return 0;
1832 }
1833
1834/* Start the operations. */
1835
1836 if (lside) {
1837 if (lsame_(transa, "N")) {
1838
1839/* Form B := alpha*inv( A )*B. */
1840
1841 if (upper) {
1842 i__1 = *n;
1843 for (j = 1; j <= i__1; ++j) {
1844 if (*alpha != 1.f) {
1845 i__2 = *m;
1846 for (i__ = 1; i__ <= i__2; ++i__) {
1847 b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
1848 ;
1849/* L30: */
1850 }
1851 }
1852 for (k = *m; k >= 1; --k) {
1853 if (b[k + j * b_dim1] != 0.f) {
1854 if (nounit) {
1855 b[k + j * b_dim1] /= a[k + k * a_dim1];
1856 }
1857 i__2 = k - 1;
1858 for (i__ = 1; i__ <= i__2; ++i__) {
1859 b[i__ + j * b_dim1] -= b[k + j * b_dim1] * a[
1860 i__ + k * a_dim1];
1861/* L40: */
1862 }
1863 }
1864/* L50: */
1865 }
1866/* L60: */
1867 }
1868 } else {
1869 i__1 = *n;
1870 for (j = 1; j <= i__1; ++j) {
1871 if (*alpha != 1.f) {
1872 i__2 = *m;
1873 for (i__ = 1; i__ <= i__2; ++i__) {
1874 b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
1875 ;
1876/* L70: */
1877 }
1878 }
1879 i__2 = *m;
1880 for (k = 1; k <= i__2; ++k) {
1881 if (b[k + j * b_dim1] != 0.f) {
1882 if (nounit) {
1883 b[k + j * b_dim1] /= a[k + k * a_dim1];
1884 }
1885 i__3 = *m;
1886 for (i__ = k + 1; i__ <= i__3; ++i__) {
1887 b[i__ + j * b_dim1] -= b[k + j * b_dim1] * a[
1888 i__ + k * a_dim1];
1889/* L80: */
1890 }
1891 }
1892/* L90: */
1893 }
1894/* L100: */
1895 }
1896 }
1897 } else {
1898
1899/* Form B := alpha*inv( A' )*B. */
1900
1901 if (upper) {
1902 i__1 = *n;
1903 for (j = 1; j <= i__1; ++j) {
1904 i__2 = *m;
1905 for (i__ = 1; i__ <= i__2; ++i__) {
1906 temp = *alpha * b[i__ + j * b_dim1];
1907 i__3 = i__ - 1;
1908 for (k = 1; k <= i__3; ++k) {
1909 temp -= a[k + i__ * a_dim1] * b[k + j * b_dim1];
1910/* L110: */
1911 }
1912 if (nounit) {
1913 temp /= a[i__ + i__ * a_dim1];
1914 }
1915 b[i__ + j * b_dim1] = temp;
1916/* L120: */
1917 }
1918/* L130: */
1919 }
1920 } else {
1921 i__1 = *n;
1922 for (j = 1; j <= i__1; ++j) {
1923 for (i__ = *m; i__ >= 1; --i__) {
1924 temp = *alpha * b[i__ + j * b_dim1];
1925 i__2 = *m;
1926 for (k = i__ + 1; k <= i__2; ++k) {
1927 temp -= a[k + i__ * a_dim1] * b[k + j * b_dim1];
1928/* L140: */
1929 }
1930 if (nounit) {
1931 temp /= a[i__ + i__ * a_dim1];
1932 }
1933 b[i__ + j * b_dim1] = temp;
1934/* L150: */
1935 }
1936/* L160: */
1937 }
1938 }
1939 }
1940 } else {
1941 if (lsame_(transa, "N")) {
1942
1943/* Form B := alpha*B*inv( A ). */
1944
1945 if (upper) {
1946 i__1 = *n;
1947 for (j = 1; j <= i__1; ++j) {
1948 if (*alpha != 1.f) {
1949 i__2 = *m;
1950 for (i__ = 1; i__ <= i__2; ++i__) {
1951 b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
1952 ;
1953/* L170: */
1954 }
1955 }
1956 i__2 = j - 1;
1957 for (k = 1; k <= i__2; ++k) {
1958 if (a[k + j * a_dim1] != 0.f) {
1959 i__3 = *m;
1960 for (i__ = 1; i__ <= i__3; ++i__) {
1961 b[i__ + j * b_dim1] -= a[k + j * a_dim1] * b[
1962 i__ + k * b_dim1];
1963/* L180: */
1964 }
1965 }
1966/* L190: */
1967 }
1968 if (nounit) {
1969 temp = 1.f / a[j + j * a_dim1];
1970 i__2 = *m;
1971 for (i__ = 1; i__ <= i__2; ++i__) {
1972 b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
1973/* L200: */
1974 }
1975 }
1976/* L210: */
1977 }
1978 } else {
1979 for (j = *n; j >= 1; --j) {
1980 if (*alpha != 1.f) {
1981 i__1 = *m;
1982 for (i__ = 1; i__ <= i__1; ++i__) {
1983 b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
1984 ;
1985/* L220: */
1986 }
1987 }
1988 i__1 = *n;
1989 for (k = j + 1; k <= i__1; ++k) {
1990 if (a[k + j * a_dim1] != 0.f) {
1991 i__2 = *m;
1992 for (i__ = 1; i__ <= i__2; ++i__) {
1993 b[i__ + j * b_dim1] -= a[k + j * a_dim1] * b[
1994 i__ + k * b_dim1];
1995/* L230: */
1996 }
1997 }
1998/* L240: */
1999 }
2000 if (nounit) {
2001 temp = 1.f / a[j + j * a_dim1];
2002 i__1 = *m;
2003 for (i__ = 1; i__ <= i__1; ++i__) {
2004 b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
2005/* L250: */
2006 }
2007 }
2008/* L260: */
2009 }
2010 }
2011 } else {
2012
2013/* Form B := alpha*B*inv( A' ). */
2014
2015 if (upper) {
2016 for (k = *n; k >= 1; --k) {
2017 if (nounit) {
2018 temp = 1.f / a[k + k * a_dim1];
2019 i__1 = *m;
2020 for (i__ = 1; i__ <= i__1; ++i__) {
2021 b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
2022/* L270: */
2023 }
2024 }
2025 i__1 = k - 1;
2026 for (j = 1; j <= i__1; ++j) {
2027 if (a[j + k * a_dim1] != 0.f) {
2028 temp = a[j + k * a_dim1];
2029 i__2 = *m;
2030 for (i__ = 1; i__ <= i__2; ++i__) {
2031 b[i__ + j * b_dim1] -= temp * b[i__ + k *
2032 b_dim1];
2033/* L280: */
2034 }
2035 }
2036/* L290: */
2037 }
2038 if (*alpha != 1.f) {
2039 i__1 = *m;
2040 for (i__ = 1; i__ <= i__1; ++i__) {
2041 b[i__ + k * b_dim1] = *alpha * b[i__ + k * b_dim1]
2042 ;
2043/* L300: */
2044 }
2045 }
2046/* L310: */
2047 }
2048 } else {
2049 i__1 = *n;
2050 for (k = 1; k <= i__1; ++k) {
2051 if (nounit) {
2052 temp = 1.f / a[k + k * a_dim1];
2053 i__2 = *m;
2054 for (i__ = 1; i__ <= i__2; ++i__) {
2055 b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
2056/* L320: */
2057 }
2058 }
2059 i__2 = *n;
2060 for (j = k + 1; j <= i__2; ++j) {
2061 if (a[j + k * a_dim1] != 0.f) {
2062 temp = a[j + k * a_dim1];
2063 i__3 = *m;
2064 for (i__ = 1; i__ <= i__3; ++i__) {
2065 b[i__ + j * b_dim1] -= temp * b[i__ + k *
2066 b_dim1];
2067/* L330: */
2068 }
2069 }
2070/* L340: */
2071 }
2072 if (*alpha != 1.f) {
2073 i__2 = *m;
2074 for (i__ = 1; i__ <= i__2; ++i__) {
2075 b[i__ + k * b_dim1] = *alpha * b[i__ + k * b_dim1]
2076 ;
2077/* L350: */
2078 }
2079 }
2080/* L360: */
2081 }
2082 }
2083 }
2084 }
2085
2086 return 0;
2087
2088/* End of STRSM . */
2089
2090} /* strsm_ */
2091
2092/* Subroutine */ int xerbla_(char *srname, integer *info)
2093{
2094 /* Format strings */
2095 static char fmt_9999[] = "(\002 ** On entry to \002,a6,\002 parameter nu"
2096 "mber \002,i2,\002 had \002,\002an illegal value\002)";
2097
2098 /* Builtin functions */
2099 integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
2100 /* Subroutine */ int s_stop(char *, ftnlen);
2101
2102 /* Fortran I/O blocks */
2103 static cilist io___60 = { 0, 6, 0, fmt_9999, 0 };
2104
2105
2106/*
2107 -- LAPACK auxiliary routine (preliminary version) --
2108 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
2109 Courant Institute, Argonne National Lab, and Rice University
2110 February 29, 1992
2111
2112
2113 Purpose
2114 =======
2115
2116 XERBLA is an error handler for the LAPACK routines.
2117 It is called by an LAPACK routine if an input parameter has an
2118 invalid value. A message is printed and execution stops.
2119
2120 Installers may consider modifying the STOP statement in order to
2121 call system-specific exception-handling facilities.
2122
2123 Arguments
2124 =========
2125
2126 SRNAME (input) CHARACTER*6
2127 The name of the routine which called XERBLA.
2128
2129 INFO (input) INTEGER
2130 The position of the invalid parameter in the parameter list
2131 of the calling routine.
2132*/
2133
2134
2135 s_wsfe(&io___60);
2136 do_fio(&c__1, srname, (ftnlen)6);
2137 do_fio(&c__1, (char *)&(*info), (ftnlen)sizeof(integer));
2138 e_wsfe();
2139
2140 s_stop("", (ftnlen)0);
2141
2142
2143/* End of XERBLA */
2144
2145 return 0;
2146} /* xerbla_ */
2147
Definition f2c.h:46