SphinxBase  5prealpha
blas_lite.c
1 /*
2 NOTE: 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
10 extern 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 
18 extern doublereal slapy2_(real *, real *);
19 
20 
21 
22 /* Table of constant values */
23 
24 static integer c__1 = 1;
25 
26 logical 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 
135 doublereal 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 
198 L20:
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  }
211 L40:
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  }
220 L60:
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 
922 L20:
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  }
935 L40:
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 
cilist
Definition: f2c.h:45