summaryrefslogtreecommitdiff
path: root/media/sphinxbase/src/libsphinxbase/util/blas_lite.c
diff options
context:
space:
mode:
Diffstat (limited to 'media/sphinxbase/src/libsphinxbase/util/blas_lite.c')
-rw-r--r--media/sphinxbase/src/libsphinxbase/util/blas_lite.c2147
1 files changed, 0 insertions, 2147 deletions
diff --git a/media/sphinxbase/src/libsphinxbase/util/blas_lite.c b/media/sphinxbase/src/libsphinxbase/util/blas_lite.c
deleted file mode 100644
index c175eaa528..0000000000
--- a/media/sphinxbase/src/libsphinxbase/util/blas_lite.c
+++ /dev/null
@@ -1,2147 +0,0 @@
-/*
-NOTE: This is generated code. Look in README.python for information on
- remaking this file.
-*/
-#include "sphinxbase/f2c.h"
-
-#ifdef HAVE_CONFIG
-#include "config.h"
-#else
-extern doublereal slamch_(char *);
-#define EPSILON slamch_("Epsilon")
-#define SAFEMINIMUM slamch_("Safe minimum")
-#define PRECISION slamch_("Precision")
-#define BASE slamch_("Base")
-#endif
-
-
-extern doublereal slapy2_(real *, real *);
-
-
-
-/* Table of constant values */
-
-static integer c__1 = 1;
-
-logical lsame_(char *ca, char *cb)
-{
- /* System generated locals */
- logical ret_val;
-
- /* Local variables */
- static integer inta, intb, zcode;
-
-
-/*
- -- LAPACK auxiliary routine (version 3.0) --
- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
- Courant Institute, Argonne National Lab, and Rice University
- September 30, 1994
-
-
- Purpose
- =======
-
- LSAME returns .TRUE. if CA is the same letter as CB regardless of
- case.
-
- Arguments
- =========
-
- CA (input) CHARACTER*1
- CB (input) CHARACTER*1
- CA and CB specify the single characters to be compared.
-
- =====================================================================
-
-
- Test if the characters are equal
-*/
-
- ret_val = *(unsigned char *)ca == *(unsigned char *)cb;
- if (ret_val) {
- return ret_val;
- }
-
-/* Now test for equivalence if both characters are alphabetic. */
-
- zcode = 'Z';
-
-/*
- Use 'Z' rather than 'A' so that ASCII can be detected on Prime
- machines, on which ICHAR returns a value with bit 8 set.
- ICHAR('A') on Prime machines returns 193 which is the same as
- ICHAR('A') on an EBCDIC machine.
-*/
-
- inta = *(unsigned char *)ca;
- intb = *(unsigned char *)cb;
-
- if (zcode == 90 || zcode == 122) {
-
-/*
- ASCII is assumed - ZCODE is the ASCII code of either lower or
- upper case 'Z'.
-*/
-
- if (inta >= 97 && inta <= 122) {
- inta += -32;
- }
- if (intb >= 97 && intb <= 122) {
- intb += -32;
- }
-
- } else if (zcode == 233 || zcode == 169) {
-
-/*
- EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
- upper case 'Z'.
-*/
-
- if (inta >= 129 && inta <= 137 || inta >= 145 && inta <= 153 || inta
- >= 162 && inta <= 169) {
- inta += 64;
- }
- if (intb >= 129 && intb <= 137 || intb >= 145 && intb <= 153 || intb
- >= 162 && intb <= 169) {
- intb += 64;
- }
-
- } else if (zcode == 218 || zcode == 250) {
-
-/*
- ASCII is assumed, on Prime machines - ZCODE is the ASCII code
- plus 128 of either lower or upper case 'Z'.
-*/
-
- if (inta >= 225 && inta <= 250) {
- inta += -32;
- }
- if (intb >= 225 && intb <= 250) {
- intb += -32;
- }
- }
- ret_val = inta == intb;
-
-/*
- RETURN
-
- End of LSAME
-*/
-
- return ret_val;
-} /* lsame_ */
-
-doublereal sdot_(integer *n, real *sx, integer *incx, real *sy, integer *incy)
-{
- /* System generated locals */
- integer i__1;
- real ret_val;
-
- /* Local variables */
- static integer i__, m, ix, iy, mp1;
- static real stemp;
-
-
-/*
- forms the dot product of two vectors.
- uses unrolled loops for increments equal to one.
- jack dongarra, linpack, 3/11/78.
- modified 12/3/93, array(1) declarations changed to array(*)
-*/
-
-
- /* Parameter adjustments */
- --sy;
- --sx;
-
- /* Function Body */
- stemp = 0.f;
- ret_val = 0.f;
- if (*n <= 0) {
- return ret_val;
- }
- if (*incx == 1 && *incy == 1) {
- goto L20;
- }
-
-/*
- code for unequal increments or equal increments
- not equal to 1
-*/
-
- ix = 1;
- iy = 1;
- if (*incx < 0) {
- ix = (-(*n) + 1) * *incx + 1;
- }
- if (*incy < 0) {
- iy = (-(*n) + 1) * *incy + 1;
- }
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__) {
- stemp += sx[ix] * sy[iy];
- ix += *incx;
- iy += *incy;
-/* L10: */
- }
- ret_val = stemp;
- return ret_val;
-
-/*
- code for both increments equal to 1
-
-
- clean-up loop
-*/
-
-L20:
- m = *n % 5;
- if (m == 0) {
- goto L40;
- }
- i__1 = m;
- for (i__ = 1; i__ <= i__1; ++i__) {
- stemp += sx[i__] * sy[i__];
-/* L30: */
- }
- if (*n < 5) {
- goto L60;
- }
-L40:
- mp1 = m + 1;
- i__1 = *n;
- for (i__ = mp1; i__ <= i__1; i__ += 5) {
- stemp = stemp + sx[i__] * sy[i__] + sx[i__ + 1] * sy[i__ + 1] + sx[
- i__ + 2] * sy[i__ + 2] + sx[i__ + 3] * sy[i__ + 3] + sx[i__ +
- 4] * sy[i__ + 4];
-/* L50: */
- }
-L60:
- ret_val = stemp;
- return ret_val;
-} /* sdot_ */
-
-/* Subroutine */ int sgemm_(char *transa, char *transb, integer *m, integer *
- n, integer *k, real *alpha, real *a, integer *lda, real *b, integer *
- ldb, real *beta, real *c__, integer *ldc)
-{
- /* System generated locals */
- integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
- i__3;
-
- /* Local variables */
- static integer i__, j, l, info;
- static logical nota, notb;
- static real temp;
- static integer ncola;
- extern logical lsame_(char *, char *);
- static integer nrowa, nrowb;
- extern /* Subroutine */ int xerbla_(char *, integer *);
-
-
-/*
- Purpose
- =======
-
- SGEMM performs one of the matrix-matrix operations
-
- C := alpha*op( A )*op( B ) + beta*C,
-
- where op( X ) is one of
-
- op( X ) = X or op( X ) = X',
-
- alpha and beta are scalars, and A, B and C are matrices, with op( A )
- an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
-
- Parameters
- ==========
-
- TRANSA - CHARACTER*1.
- On entry, TRANSA specifies the form of op( A ) to be used in
- the matrix multiplication as follows:
-
- TRANSA = 'N' or 'n', op( A ) = A.
-
- TRANSA = 'T' or 't', op( A ) = A'.
-
- TRANSA = 'C' or 'c', op( A ) = A'.
-
- Unchanged on exit.
-
- TRANSB - CHARACTER*1.
- On entry, TRANSB specifies the form of op( B ) to be used in
- the matrix multiplication as follows:
-
- TRANSB = 'N' or 'n', op( B ) = B.
-
- TRANSB = 'T' or 't', op( B ) = B'.
-
- TRANSB = 'C' or 'c', op( B ) = B'.
-
- Unchanged on exit.
-
- M - INTEGER.
- On entry, M specifies the number of rows of the matrix
- op( A ) and of the matrix C. M must be at least zero.
- Unchanged on exit.
-
- N - INTEGER.
- On entry, N specifies the number of columns of the matrix
- op( B ) and the number of columns of the matrix C. N must be
- at least zero.
- Unchanged on exit.
-
- K - INTEGER.
- On entry, K specifies the number of columns of the matrix
- op( A ) and the number of rows of the matrix op( B ). K must
- be at least zero.
- Unchanged on exit.
-
- ALPHA - REAL .
- On entry, ALPHA specifies the scalar alpha.
- Unchanged on exit.
-
- A - REAL array of DIMENSION ( LDA, ka ), where ka is
- k when TRANSA = 'N' or 'n', and is m otherwise.
- Before entry with TRANSA = 'N' or 'n', the leading m by k
- part of the array A must contain the matrix A, otherwise
- the leading k by m part of the array A must contain the
- matrix A.
- Unchanged on exit.
-
- LDA - INTEGER.
- On entry, LDA specifies the first dimension of A as declared
- in the calling (sub) program. When TRANSA = 'N' or 'n' then
- LDA must be at least max( 1, m ), otherwise LDA must be at
- least max( 1, k ).
- Unchanged on exit.
-
- B - REAL array of DIMENSION ( LDB, kb ), where kb is
- n when TRANSB = 'N' or 'n', and is k otherwise.
- Before entry with TRANSB = 'N' or 'n', the leading k by n
- part of the array B must contain the matrix B, otherwise
- the leading n by k part of the array B must contain the
- matrix B.
- Unchanged on exit.
-
- LDB - INTEGER.
- On entry, LDB specifies the first dimension of B as declared
- in the calling (sub) program. When TRANSB = 'N' or 'n' then
- LDB must be at least max( 1, k ), otherwise LDB must be at
- least max( 1, n ).
- Unchanged on exit.
-
- BETA - REAL .
- On entry, BETA specifies the scalar beta. When BETA is
- supplied as zero then C need not be set on input.
- Unchanged on exit.
-
- C - REAL array of DIMENSION ( LDC, n ).
- Before entry, the leading m by n part of the array C must
- contain the matrix C, except when beta is zero, in which
- case C need not be set on entry.
- On exit, the array C is overwritten by the m by n matrix
- ( alpha*op( A )*op( B ) + beta*C ).
-
- LDC - INTEGER.
- On entry, LDC specifies the first dimension of C as declared
- in the calling (sub) program. LDC must be at least
- max( 1, m ).
- Unchanged on exit.
-
-
- Level 3 Blas routine.
-
- -- Written on 8-February-1989.
- Jack Dongarra, Argonne National Laboratory.
- Iain Duff, AERE Harwell.
- Jeremy Du Croz, Numerical Algorithms Group Ltd.
- Sven Hammarling, Numerical Algorithms Group Ltd.
-
-
- Set NOTA and NOTB as true if A and B respectively are not
- transposed and set NROWA, NCOLA and NROWB as the number of rows
- and columns of A and the number of rows of B respectively.
-*/
-
- /* Parameter adjustments */
- a_dim1 = *lda;
- a_offset = 1 + a_dim1;
- a -= a_offset;
- b_dim1 = *ldb;
- b_offset = 1 + b_dim1;
- b -= b_offset;
- c_dim1 = *ldc;
- c_offset = 1 + c_dim1;
- c__ -= c_offset;
-
- /* Function Body */
- nota = lsame_(transa, "N");
- notb = lsame_(transb, "N");
- if (nota) {
- nrowa = *m;
- ncola = *k;
- } else {
- nrowa = *k;
- ncola = *m;
- }
- if (notb) {
- nrowb = *k;
- } else {
- nrowb = *n;
- }
-
-/* Test the input parameters. */
-
- info = 0;
- if (! nota && ! lsame_(transa, "C") && ! lsame_(
- transa, "T")) {
- info = 1;
- } else if (! notb && ! lsame_(transb, "C") && !
- lsame_(transb, "T")) {
- info = 2;
- } else if (*m < 0) {
- info = 3;
- } else if (*n < 0) {
- info = 4;
- } else if (*k < 0) {
- info = 5;
- } else if (*lda < max(1,nrowa)) {
- info = 8;
- } else if (*ldb < max(1,nrowb)) {
- info = 10;
- } else if (*ldc < max(1,*m)) {
- info = 13;
- }
- if (info != 0) {
- xerbla_("SGEMM ", &info);
- return 0;
- }
-
-/* Quick return if possible. */
-
- if (*m == 0 || *n == 0 || (*alpha == 0.f || *k == 0) && *beta == 1.f) {
- return 0;
- }
-
-/* And if alpha.eq.zero. */
-
- if (*alpha == 0.f) {
- if (*beta == 0.f) {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = 0.f;
-/* L10: */
- }
-/* L20: */
- }
- } else {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
-/* L30: */
- }
-/* L40: */
- }
- }
- return 0;
- }
-
-/* Start the operations. */
-
- if (notb) {
- if (nota) {
-
-/* Form C := alpha*A*B + beta*C. */
-
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- if (*beta == 0.f) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = 0.f;
-/* L50: */
- }
- } else if (*beta != 1.f) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
-/* L60: */
- }
- }
- i__2 = *k;
- for (l = 1; l <= i__2; ++l) {
- if (b[l + j * b_dim1] != 0.f) {
- temp = *alpha * b[l + j * b_dim1];
- i__3 = *m;
- for (i__ = 1; i__ <= i__3; ++i__) {
- c__[i__ + j * c_dim1] += temp * a[i__ + l *
- a_dim1];
-/* L70: */
- }
- }
-/* L80: */
- }
-/* L90: */
- }
- } else {
-
-/* Form C := alpha*A'*B + beta*C */
-
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- temp = 0.f;
- i__3 = *k;
- for (l = 1; l <= i__3; ++l) {
- temp += a[l + i__ * a_dim1] * b[l + j * b_dim1];
-/* L100: */
- }
- if (*beta == 0.f) {
- c__[i__ + j * c_dim1] = *alpha * temp;
- } else {
- c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
- i__ + j * c_dim1];
- }
-/* L110: */
- }
-/* L120: */
- }
- }
- } else {
- if (nota) {
-
-/* Form C := alpha*A*B' + beta*C */
-
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- if (*beta == 0.f) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = 0.f;
-/* L130: */
- }
- } else if (*beta != 1.f) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
-/* L140: */
- }
- }
- i__2 = *k;
- for (l = 1; l <= i__2; ++l) {
- if (b[j + l * b_dim1] != 0.f) {
- temp = *alpha * b[j + l * b_dim1];
- i__3 = *m;
- for (i__ = 1; i__ <= i__3; ++i__) {
- c__[i__ + j * c_dim1] += temp * a[i__ + l *
- a_dim1];
-/* L150: */
- }
- }
-/* L160: */
- }
-/* L170: */
- }
- } else {
-
-/* Form C := alpha*A'*B' + beta*C */
-
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- temp = 0.f;
- i__3 = *k;
- for (l = 1; l <= i__3; ++l) {
- temp += a[l + i__ * a_dim1] * b[j + l * b_dim1];
-/* L180: */
- }
- if (*beta == 0.f) {
- c__[i__ + j * c_dim1] = *alpha * temp;
- } else {
- c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
- i__ + j * c_dim1];
- }
-/* L190: */
- }
-/* L200: */
- }
- }
- }
-
- return 0;
-
-/* End of SGEMM . */
-
-} /* sgemm_ */
-
-/* Subroutine */ int sgemv_(char *trans, integer *m, integer *n, real *alpha,
- real *a, integer *lda, real *x, integer *incx, real *beta, real *y,
- integer *incy)
-{
- /* System generated locals */
- integer a_dim1, a_offset, i__1, i__2;
-
- /* Local variables */
- static integer i__, j, ix, iy, jx, jy, kx, ky, info;
- static real temp;
- static integer lenx, leny;
- extern logical lsame_(char *, char *);
- extern /* Subroutine */ int xerbla_(char *, integer *);
-
-
-/*
- Purpose
- =======
-
- SGEMV performs one of the matrix-vector operations
-
- y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y,
-
- where alpha and beta are scalars, x and y are vectors and A is an
- m by n matrix.
-
- Parameters
- ==========
-
- TRANS - CHARACTER*1.
- On entry, TRANS specifies the operation to be performed as
- follows:
-
- TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
-
- TRANS = 'T' or 't' y := alpha*A'*x + beta*y.
-
- TRANS = 'C' or 'c' y := alpha*A'*x + beta*y.
-
- Unchanged on exit.
-
- M - INTEGER.
- On entry, M specifies the number of rows of the matrix A.
- M must be at least zero.
- Unchanged on exit.
-
- N - INTEGER.
- On entry, N specifies the number of columns of the matrix A.
- N must be at least zero.
- Unchanged on exit.
-
- ALPHA - REAL .
- On entry, ALPHA specifies the scalar alpha.
- Unchanged on exit.
-
- A - REAL array of DIMENSION ( LDA, n ).
- Before entry, the leading m by n part of the array A must
- contain the matrix of coefficients.
- Unchanged on exit.
-
- LDA - INTEGER.
- On entry, LDA specifies the first dimension of A as declared
- in the calling (sub) program. LDA must be at least
- max( 1, m ).
- Unchanged on exit.
-
- X - REAL array of DIMENSION at least
- ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
- and at least
- ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
- Before entry, the incremented array X must contain the
- vector x.
- Unchanged on exit.
-
- INCX - INTEGER.
- On entry, INCX specifies the increment for the elements of
- X. INCX must not be zero.
- Unchanged on exit.
-
- BETA - REAL .
- On entry, BETA specifies the scalar beta. When BETA is
- supplied as zero then Y need not be set on input.
- Unchanged on exit.
-
- Y - REAL array of DIMENSION at least
- ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
- and at least
- ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
- Before entry with BETA non-zero, the incremented array Y
- must contain the vector y. On exit, Y is overwritten by the
- updated vector y.
-
- INCY - INTEGER.
- On entry, INCY specifies the increment for the elements of
- Y. INCY must not be zero.
- Unchanged on exit.
-
-
- Level 2 Blas routine.
-
- -- Written on 22-October-1986.
- Jack Dongarra, Argonne National Lab.
- Jeremy Du Croz, Nag Central Office.
- Sven Hammarling, Nag Central Office.
- Richard Hanson, Sandia National Labs.
-
-
- Test the input parameters.
-*/
-
- /* Parameter adjustments */
- a_dim1 = *lda;
- a_offset = 1 + a_dim1;
- a -= a_offset;
- --x;
- --y;
-
- /* Function Body */
- info = 0;
- if (! lsame_(trans, "N") && ! lsame_(trans, "T") && ! lsame_(trans, "C")
- ) {
- info = 1;
- } else if (*m < 0) {
- info = 2;
- } else if (*n < 0) {
- info = 3;
- } else if (*lda < max(1,*m)) {
- info = 6;
- } else if (*incx == 0) {
- info = 8;
- } else if (*incy == 0) {
- info = 11;
- }
- if (info != 0) {
- xerbla_("SGEMV ", &info);
- return 0;
- }
-
-/* Quick return if possible. */
-
- if (*m == 0 || *n == 0 || *alpha == 0.f && *beta == 1.f) {
- return 0;
- }
-
-/*
- Set LENX and LENY, the lengths of the vectors x and y, and set
- up the start points in X and Y.
-*/
-
- if (lsame_(trans, "N")) {
- lenx = *n;
- leny = *m;
- } else {
- lenx = *m;
- leny = *n;
- }
- if (*incx > 0) {
- kx = 1;
- } else {
- kx = 1 - (lenx - 1) * *incx;
- }
- if (*incy > 0) {
- ky = 1;
- } else {
- ky = 1 - (leny - 1) * *incy;
- }
-
-/*
- Start the operations. In this version the elements of A are
- accessed sequentially with one pass through A.
-
- First form y := beta*y.
-*/
-
- if (*beta != 1.f) {
- if (*incy == 1) {
- if (*beta == 0.f) {
- i__1 = leny;
- for (i__ = 1; i__ <= i__1; ++i__) {
- y[i__] = 0.f;
-/* L10: */
- }
- } else {
- i__1 = leny;
- for (i__ = 1; i__ <= i__1; ++i__) {
- y[i__] = *beta * y[i__];
-/* L20: */
- }
- }
- } else {
- iy = ky;
- if (*beta == 0.f) {
- i__1 = leny;
- for (i__ = 1; i__ <= i__1; ++i__) {
- y[iy] = 0.f;
- iy += *incy;
-/* L30: */
- }
- } else {
- i__1 = leny;
- for (i__ = 1; i__ <= i__1; ++i__) {
- y[iy] = *beta * y[iy];
- iy += *incy;
-/* L40: */
- }
- }
- }
- }
- if (*alpha == 0.f) {
- return 0;
- }
- if (lsame_(trans, "N")) {
-
-/* Form y := alpha*A*x + y. */
-
- jx = kx;
- if (*incy == 1) {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- if (x[jx] != 0.f) {
- temp = *alpha * x[jx];
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- y[i__] += temp * a[i__ + j * a_dim1];
-/* L50: */
- }
- }
- jx += *incx;
-/* L60: */
- }
- } else {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- if (x[jx] != 0.f) {
- temp = *alpha * x[jx];
- iy = ky;
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- y[iy] += temp * a[i__ + j * a_dim1];
- iy += *incy;
-/* L70: */
- }
- }
- jx += *incx;
-/* L80: */
- }
- }
- } else {
-
-/* Form y := alpha*A'*x + y. */
-
- jy = ky;
- if (*incx == 1) {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- temp = 0.f;
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- temp += a[i__ + j * a_dim1] * x[i__];
-/* L90: */
- }
- y[jy] += *alpha * temp;
- jy += *incy;
-/* L100: */
- }
- } else {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- temp = 0.f;
- ix = kx;
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- temp += a[i__ + j * a_dim1] * x[ix];
- ix += *incx;
-/* L110: */
- }
- y[jy] += *alpha * temp;
- jy += *incy;
-/* L120: */
- }
- }
- }
-
- return 0;
-
-/* End of SGEMV . */
-
-} /* sgemv_ */
-
-/* Subroutine */ int sscal_(integer *n, real *sa, real *sx, integer *incx)
-{
- /* System generated locals */
- integer i__1, i__2;
-
- /* Local variables */
- static integer i__, m, mp1, nincx;
-
-
-/*
- scales a vector by a constant.
- uses unrolled loops for increment equal to 1.
- jack dongarra, linpack, 3/11/78.
- modified 3/93 to return if incx .le. 0.
- modified 12/3/93, array(1) declarations changed to array(*)
-*/
-
-
- /* Parameter adjustments */
- --sx;
-
- /* Function Body */
- if (*n <= 0 || *incx <= 0) {
- return 0;
- }
- if (*incx == 1) {
- goto L20;
- }
-
-/* code for increment not equal to 1 */
-
- nincx = *n * *incx;
- i__1 = nincx;
- i__2 = *incx;
- for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
- sx[i__] = *sa * sx[i__];
-/* L10: */
- }
- return 0;
-
-/*
- code for increment equal to 1
-
-
- clean-up loop
-*/
-
-L20:
- m = *n % 5;
- if (m == 0) {
- goto L40;
- }
- i__2 = m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- sx[i__] = *sa * sx[i__];
-/* L30: */
- }
- if (*n < 5) {
- return 0;
- }
-L40:
- mp1 = m + 1;
- i__2 = *n;
- for (i__ = mp1; i__ <= i__2; i__ += 5) {
- sx[i__] = *sa * sx[i__];
- sx[i__ + 1] = *sa * sx[i__ + 1];
- sx[i__ + 2] = *sa * sx[i__ + 2];
- sx[i__ + 3] = *sa * sx[i__ + 3];
- sx[i__ + 4] = *sa * sx[i__ + 4];
-/* L50: */
- }
- return 0;
-} /* sscal_ */
-
-/* Subroutine */ int ssymm_(char *side, char *uplo, integer *m, integer *n,
- real *alpha, real *a, integer *lda, real *b, integer *ldb, real *beta,
- real *c__, integer *ldc)
-{
- /* System generated locals */
- integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
- i__3;
-
- /* Local variables */
- static integer i__, j, k, info;
- static real temp1, temp2;
- extern logical lsame_(char *, char *);
- static integer nrowa;
- static logical upper;
- extern /* Subroutine */ int xerbla_(char *, integer *);
-
-
-/*
- Purpose
- =======
-
- SSYMM performs one of the matrix-matrix operations
-
- C := alpha*A*B + beta*C,
-
- or
-
- C := alpha*B*A + beta*C,
-
- where alpha and beta are scalars, A is a symmetric matrix and B and
- C are m by n matrices.
-
- Parameters
- ==========
-
- SIDE - CHARACTER*1.
- On entry, SIDE specifies whether the symmetric matrix A
- appears on the left or right in the operation as follows:
-
- SIDE = 'L' or 'l' C := alpha*A*B + beta*C,
-
- SIDE = 'R' or 'r' C := alpha*B*A + beta*C,
-
- Unchanged on exit.
-
- UPLO - CHARACTER*1.
- On entry, UPLO specifies whether the upper or lower
- triangular part of the symmetric matrix A is to be
- referenced as follows:
-
- UPLO = 'U' or 'u' Only the upper triangular part of the
- symmetric matrix is to be referenced.
-
- UPLO = 'L' or 'l' Only the lower triangular part of the
- symmetric matrix is to be referenced.
-
- Unchanged on exit.
-
- M - INTEGER.
- On entry, M specifies the number of rows of the matrix C.
- M must be at least zero.
- Unchanged on exit.
-
- N - INTEGER.
- On entry, N specifies the number of columns of the matrix C.
- N must be at least zero.
- Unchanged on exit.
-
- ALPHA - REAL .
- On entry, ALPHA specifies the scalar alpha.
- Unchanged on exit.
-
- A - REAL array of DIMENSION ( LDA, ka ), where ka is
- m when SIDE = 'L' or 'l' and is n otherwise.
- Before entry with SIDE = 'L' or 'l', the m by m part of
- the array A must contain the symmetric matrix, such that
- when UPLO = 'U' or 'u', the leading m by m upper triangular
- part of the array A must contain the upper triangular part
- of the symmetric matrix and the strictly lower triangular
- part of A is not referenced, and when UPLO = 'L' or 'l',
- the leading m by m lower triangular part of the array A
- must contain the lower triangular part of the symmetric
- matrix and the strictly upper triangular part of A is not
- referenced.
- Before entry with SIDE = 'R' or 'r', the n by n part of
- the array A must contain the symmetric matrix, such that
- when UPLO = 'U' or 'u', the leading n by n upper triangular
- part of the array A must contain the upper triangular part
- of the symmetric matrix and the strictly lower triangular
- part of A is not referenced, and when UPLO = 'L' or 'l',
- the leading n by n lower triangular part of the array A
- must contain the lower triangular part of the symmetric
- matrix and the strictly upper triangular part of A is not
- referenced.
- Unchanged on exit.
-
- LDA - INTEGER.
- On entry, LDA specifies the first dimension of A as declared
- in the calling (sub) program. When SIDE = 'L' or 'l' then
- LDA must be at least max( 1, m ), otherwise LDA must be at
- least max( 1, n ).
- Unchanged on exit.
-
- B - REAL array of DIMENSION ( LDB, n ).
- Before entry, the leading m by n part of the array B must
- contain the matrix B.
- Unchanged on exit.
-
- LDB - INTEGER.
- On entry, LDB specifies the first dimension of B as declared
- in the calling (sub) program. LDB must be at least
- max( 1, m ).
- Unchanged on exit.
-
- BETA - REAL .
- On entry, BETA specifies the scalar beta. When BETA is
- supplied as zero then C need not be set on input.
- Unchanged on exit.
-
- C - REAL array of DIMENSION ( LDC, n ).
- Before entry, the leading m by n part of the array C must
- contain the matrix C, except when beta is zero, in which
- case C need not be set on entry.
- On exit, the array C is overwritten by the m by n updated
- matrix.
-
- LDC - INTEGER.
- On entry, LDC specifies the first dimension of C as declared
- in the calling (sub) program. LDC must be at least
- max( 1, m ).
- Unchanged on exit.
-
-
- Level 3 Blas routine.
-
- -- Written on 8-February-1989.
- Jack Dongarra, Argonne National Laboratory.
- Iain Duff, AERE Harwell.
- Jeremy Du Croz, Numerical Algorithms Group Ltd.
- Sven Hammarling, Numerical Algorithms Group Ltd.
-
-
- Set NROWA as the number of rows of A.
-*/
-
- /* Parameter adjustments */
- a_dim1 = *lda;
- a_offset = 1 + a_dim1;
- a -= a_offset;
- b_dim1 = *ldb;
- b_offset = 1 + b_dim1;
- b -= b_offset;
- c_dim1 = *ldc;
- c_offset = 1 + c_dim1;
- c__ -= c_offset;
-
- /* Function Body */
- if (lsame_(side, "L")) {
- nrowa = *m;
- } else {
- nrowa = *n;
- }
- upper = lsame_(uplo, "U");
-
-/* Test the input parameters. */
-
- info = 0;
- if (! lsame_(side, "L") && ! lsame_(side, "R")) {
- info = 1;
- } else if (! upper && ! lsame_(uplo, "L")) {
- info = 2;
- } else if (*m < 0) {
- info = 3;
- } else if (*n < 0) {
- info = 4;
- } else if (*lda < max(1,nrowa)) {
- info = 7;
- } else if (*ldb < max(1,*m)) {
- info = 9;
- } else if (*ldc < max(1,*m)) {
- info = 12;
- }
- if (info != 0) {
- xerbla_("SSYMM ", &info);
- return 0;
- }
-
-/* Quick return if possible. */
-
- if (*m == 0 || *n == 0 || *alpha == 0.f && *beta == 1.f) {
- return 0;
- }
-
-/* And when alpha.eq.zero. */
-
- if (*alpha == 0.f) {
- if (*beta == 0.f) {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = 0.f;
-/* L10: */
- }
-/* L20: */
- }
- } else {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
-/* L30: */
- }
-/* L40: */
- }
- }
- return 0;
- }
-
-/* Start the operations. */
-
- if (lsame_(side, "L")) {
-
-/* Form C := alpha*A*B + beta*C. */
-
- if (upper) {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- temp1 = *alpha * b[i__ + j * b_dim1];
- temp2 = 0.f;
- i__3 = i__ - 1;
- for (k = 1; k <= i__3; ++k) {
- c__[k + j * c_dim1] += temp1 * a[k + i__ * a_dim1];
- temp2 += b[k + j * b_dim1] * a[k + i__ * a_dim1];
-/* L50: */
- }
- if (*beta == 0.f) {
- c__[i__ + j * c_dim1] = temp1 * a[i__ + i__ * a_dim1]
- + *alpha * temp2;
- } else {
- c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1]
- + temp1 * a[i__ + i__ * a_dim1] + *alpha *
- temp2;
- }
-/* L60: */
- }
-/* L70: */
- }
- } else {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- for (i__ = *m; i__ >= 1; --i__) {
- temp1 = *alpha * b[i__ + j * b_dim1];
- temp2 = 0.f;
- i__2 = *m;
- for (k = i__ + 1; k <= i__2; ++k) {
- c__[k + j * c_dim1] += temp1 * a[k + i__ * a_dim1];
- temp2 += b[k + j * b_dim1] * a[k + i__ * a_dim1];
-/* L80: */
- }
- if (*beta == 0.f) {
- c__[i__ + j * c_dim1] = temp1 * a[i__ + i__ * a_dim1]
- + *alpha * temp2;
- } else {
- c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1]
- + temp1 * a[i__ + i__ * a_dim1] + *alpha *
- temp2;
- }
-/* L90: */
- }
-/* L100: */
- }
- }
- } else {
-
-/* Form C := alpha*B*A + beta*C. */
-
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- temp1 = *alpha * a[j + j * a_dim1];
- if (*beta == 0.f) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = temp1 * b[i__ + j * b_dim1];
-/* L110: */
- }
- } else {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1] +
- temp1 * b[i__ + j * b_dim1];
-/* L120: */
- }
- }
- i__2 = j - 1;
- for (k = 1; k <= i__2; ++k) {
- if (upper) {
- temp1 = *alpha * a[k + j * a_dim1];
- } else {
- temp1 = *alpha * a[j + k * a_dim1];
- }
- i__3 = *m;
- for (i__ = 1; i__ <= i__3; ++i__) {
- c__[i__ + j * c_dim1] += temp1 * b[i__ + k * b_dim1];
-/* L130: */
- }
-/* L140: */
- }
- i__2 = *n;
- for (k = j + 1; k <= i__2; ++k) {
- if (upper) {
- temp1 = *alpha * a[j + k * a_dim1];
- } else {
- temp1 = *alpha * a[k + j * a_dim1];
- }
- i__3 = *m;
- for (i__ = 1; i__ <= i__3; ++i__) {
- c__[i__ + j * c_dim1] += temp1 * b[i__ + k * b_dim1];
-/* L150: */
- }
-/* L160: */
- }
-/* L170: */
- }
- }
-
- return 0;
-
-/* End of SSYMM . */
-
-} /* ssymm_ */
-
-/* Subroutine */ int ssyrk_(char *uplo, char *trans, integer *n, integer *k,
- real *alpha, real *a, integer *lda, real *beta, real *c__, integer *
- ldc)
-{
- /* System generated locals */
- integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2, i__3;
-
- /* Local variables */
- static integer i__, j, l, info;
- static real temp;
- extern logical lsame_(char *, char *);
- static integer nrowa;
- static logical upper;
- extern /* Subroutine */ int xerbla_(char *, integer *);
-
-
-/*
- Purpose
- =======
-
- SSYRK performs one of the symmetric rank k operations
-
- C := alpha*A*A' + beta*C,
-
- or
-
- C := alpha*A'*A + beta*C,
-
- where alpha and beta are scalars, C is an n by n symmetric matrix
- and A is an n by k matrix in the first case and a k by n matrix
- in the second case.
-
- Parameters
- ==========
-
- UPLO - CHARACTER*1.
- On entry, UPLO specifies whether the upper or lower
- triangular part of the array C is to be referenced as
- follows:
-
- UPLO = 'U' or 'u' Only the upper triangular part of C
- is to be referenced.
-
- UPLO = 'L' or 'l' Only the lower triangular part of C
- is to be referenced.
-
- Unchanged on exit.
-
- TRANS - CHARACTER*1.
- On entry, TRANS specifies the operation to be performed as
- follows:
-
- TRANS = 'N' or 'n' C := alpha*A*A' + beta*C.
-
- TRANS = 'T' or 't' C := alpha*A'*A + beta*C.
-
- TRANS = 'C' or 'c' C := alpha*A'*A + beta*C.
-
- Unchanged on exit.
-
- N - INTEGER.
- On entry, N specifies the order of the matrix C. N must be
- at least zero.
- Unchanged on exit.
-
- K - INTEGER.
- On entry with TRANS = 'N' or 'n', K specifies the number
- of columns of the matrix A, and on entry with
- TRANS = 'T' or 't' or 'C' or 'c', K specifies the number
- of rows of the matrix A. K must be at least zero.
- Unchanged on exit.
-
- ALPHA - REAL .
- On entry, ALPHA specifies the scalar alpha.
- Unchanged on exit.
-
- A - REAL array of DIMENSION ( LDA, ka ), where ka is
- k when TRANS = 'N' or 'n', and is n otherwise.
- Before entry with TRANS = 'N' or 'n', the leading n by k
- part of the array A must contain the matrix A, otherwise
- the leading k by n part of the array A must contain the
- matrix A.
- Unchanged on exit.
-
- LDA - INTEGER.
- On entry, LDA specifies the first dimension of A as declared
- in the calling (sub) program. When TRANS = 'N' or 'n'
- then LDA must be at least max( 1, n ), otherwise LDA must
- be at least max( 1, k ).
- Unchanged on exit.
-
- BETA - REAL .
- On entry, BETA specifies the scalar beta.
- Unchanged on exit.
-
- C - REAL array of DIMENSION ( LDC, n ).
- Before entry with UPLO = 'U' or 'u', the leading n by n
- upper triangular part of the array C must contain the upper
- triangular part of the symmetric matrix and the strictly
- lower triangular part of C is not referenced. On exit, the
- upper triangular part of the array C is overwritten by the
- upper triangular part of the updated matrix.
- Before entry with UPLO = 'L' or 'l', the leading n by n
- lower triangular part of the array C must contain the lower
- triangular part of the symmetric matrix and the strictly
- upper triangular part of C is not referenced. On exit, the
- lower triangular part of the array C is overwritten by the
- lower triangular part of the updated matrix.
-
- LDC - INTEGER.
- On entry, LDC specifies the first dimension of C as declared
- in the calling (sub) program. LDC must be at least
- max( 1, n ).
- Unchanged on exit.
-
-
- Level 3 Blas routine.
-
- -- Written on 8-February-1989.
- Jack Dongarra, Argonne National Laboratory.
- Iain Duff, AERE Harwell.
- Jeremy Du Croz, Numerical Algorithms Group Ltd.
- Sven Hammarling, Numerical Algorithms Group Ltd.
-
-
- Test the input parameters.
-*/
-
- /* Parameter adjustments */
- a_dim1 = *lda;
- a_offset = 1 + a_dim1;
- a -= a_offset;
- c_dim1 = *ldc;
- c_offset = 1 + c_dim1;
- c__ -= c_offset;
-
- /* Function Body */
- if (lsame_(trans, "N")) {
- nrowa = *n;
- } else {
- nrowa = *k;
- }
- upper = lsame_(uplo, "U");
-
- info = 0;
- if (! upper && ! lsame_(uplo, "L")) {
- info = 1;
- } else if (! lsame_(trans, "N") && ! lsame_(trans,
- "T") && ! lsame_(trans, "C")) {
- info = 2;
- } else if (*n < 0) {
- info = 3;
- } else if (*k < 0) {
- info = 4;
- } else if (*lda < max(1,nrowa)) {
- info = 7;
- } else if (*ldc < max(1,*n)) {
- info = 10;
- }
- if (info != 0) {
- xerbla_("SSYRK ", &info);
- return 0;
- }
-
-/* Quick return if possible. */
-
- if (*n == 0 || (*alpha == 0.f || *k == 0) && *beta == 1.f) {
- return 0;
- }
-
-/* And when alpha.eq.zero. */
-
- if (*alpha == 0.f) {
- if (upper) {
- if (*beta == 0.f) {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__2 = j;
- for (i__ = 1; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = 0.f;
-/* L10: */
- }
-/* L20: */
- }
- } else {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__2 = j;
- for (i__ = 1; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
-/* L30: */
- }
-/* L40: */
- }
- }
- } else {
- if (*beta == 0.f) {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__2 = *n;
- for (i__ = j; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = 0.f;
-/* L50: */
- }
-/* L60: */
- }
- } else {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__2 = *n;
- for (i__ = j; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
-/* L70: */
- }
-/* L80: */
- }
- }
- }
- return 0;
- }
-
-/* Start the operations. */
-
- if (lsame_(trans, "N")) {
-
-/* Form C := alpha*A*A' + beta*C. */
-
- if (upper) {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- if (*beta == 0.f) {
- i__2 = j;
- for (i__ = 1; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = 0.f;
-/* L90: */
- }
- } else if (*beta != 1.f) {
- i__2 = j;
- for (i__ = 1; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
-/* L100: */
- }
- }
- i__2 = *k;
- for (l = 1; l <= i__2; ++l) {
- if (a[j + l * a_dim1] != 0.f) {
- temp = *alpha * a[j + l * a_dim1];
- i__3 = j;
- for (i__ = 1; i__ <= i__3; ++i__) {
- c__[i__ + j * c_dim1] += temp * a[i__ + l *
- a_dim1];
-/* L110: */
- }
- }
-/* L120: */
- }
-/* L130: */
- }
- } else {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- if (*beta == 0.f) {
- i__2 = *n;
- for (i__ = j; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = 0.f;
-/* L140: */
- }
- } else if (*beta != 1.f) {
- i__2 = *n;
- for (i__ = j; i__ <= i__2; ++i__) {
- c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
-/* L150: */
- }
- }
- i__2 = *k;
- for (l = 1; l <= i__2; ++l) {
- if (a[j + l * a_dim1] != 0.f) {
- temp = *alpha * a[j + l * a_dim1];
- i__3 = *n;
- for (i__ = j; i__ <= i__3; ++i__) {
- c__[i__ + j * c_dim1] += temp * a[i__ + l *
- a_dim1];
-/* L160: */
- }
- }
-/* L170: */
- }
-/* L180: */
- }
- }
- } else {
-
-/* Form C := alpha*A'*A + beta*C. */
-
- if (upper) {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__2 = j;
- for (i__ = 1; i__ <= i__2; ++i__) {
- temp = 0.f;
- i__3 = *k;
- for (l = 1; l <= i__3; ++l) {
- temp += a[l + i__ * a_dim1] * a[l + j * a_dim1];
-/* L190: */
- }
- if (*beta == 0.f) {
- c__[i__ + j * c_dim1] = *alpha * temp;
- } else {
- c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
- i__ + j * c_dim1];
- }
-/* L200: */
- }
-/* L210: */
- }
- } else {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__2 = *n;
- for (i__ = j; i__ <= i__2; ++i__) {
- temp = 0.f;
- i__3 = *k;
- for (l = 1; l <= i__3; ++l) {
- temp += a[l + i__ * a_dim1] * a[l + j * a_dim1];
-/* L220: */
- }
- if (*beta == 0.f) {
- c__[i__ + j * c_dim1] = *alpha * temp;
- } else {
- c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
- i__ + j * c_dim1];
- }
-/* L230: */
- }
-/* L240: */
- }
- }
- }
-
- return 0;
-
-/* End of SSYRK . */
-
-} /* ssyrk_ */
-
-/* Subroutine */ int strsm_(char *side, char *uplo, char *transa, char *diag,
- integer *m, integer *n, real *alpha, real *a, integer *lda, real *b,
- integer *ldb)
-{
- /* System generated locals */
- integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
-
- /* Local variables */
- static integer i__, j, k, info;
- static real temp;
- static logical lside;
- extern logical lsame_(char *, char *);
- static integer nrowa;
- static logical upper;
- extern /* Subroutine */ int xerbla_(char *, integer *);
- static logical nounit;
-
-
-/*
- Purpose
- =======
-
- STRSM solves one of the matrix equations
-
- op( A )*X = alpha*B, or X*op( A ) = alpha*B,
-
- where alpha is a scalar, X and B are m by n matrices, A is a unit, or
- non-unit, upper or lower triangular matrix and op( A ) is one of
-
- op( A ) = A or op( A ) = A'.
-
- The matrix X is overwritten on B.
-
- Parameters
- ==========
-
- SIDE - CHARACTER*1.
- On entry, SIDE specifies whether op( A ) appears on the left
- or right of X as follows:
-
- SIDE = 'L' or 'l' op( A )*X = alpha*B.
-
- SIDE = 'R' or 'r' X*op( A ) = alpha*B.
-
- Unchanged on exit.
-
- UPLO - CHARACTER*1.
- On entry, UPLO specifies whether the matrix A is an upper or
- lower triangular matrix as follows:
-
- UPLO = 'U' or 'u' A is an upper triangular matrix.
-
- UPLO = 'L' or 'l' A is a lower triangular matrix.
-
- Unchanged on exit.
-
- TRANSA - CHARACTER*1.
- On entry, TRANSA specifies the form of op( A ) to be used in
- the matrix multiplication as follows:
-
- TRANSA = 'N' or 'n' op( A ) = A.
-
- TRANSA = 'T' or 't' op( A ) = A'.
-
- TRANSA = 'C' or 'c' op( A ) = A'.
-
- Unchanged on exit.
-
- DIAG - CHARACTER*1.
- On entry, DIAG specifies whether or not A is unit triangular
- as follows:
-
- DIAG = 'U' or 'u' A is assumed to be unit triangular.
-
- DIAG = 'N' or 'n' A is not assumed to be unit
- triangular.
-
- Unchanged on exit.
-
- M - INTEGER.
- On entry, M specifies the number of rows of B. M must be at
- least zero.
- Unchanged on exit.
-
- N - INTEGER.
- On entry, N specifies the number of columns of B. N must be
- at least zero.
- Unchanged on exit.
-
- ALPHA - REAL .
- On entry, ALPHA specifies the scalar alpha. When alpha is
- zero then A is not referenced and B need not be set before
- entry.
- Unchanged on exit.
-
- A - REAL array of DIMENSION ( LDA, k ), where k is m
- when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
- Before entry with UPLO = 'U' or 'u', the leading k by k
- upper triangular part of the array A must contain the upper
- triangular matrix and the strictly lower triangular part of
- A is not referenced.
- Before entry with UPLO = 'L' or 'l', the leading k by k
- lower triangular part of the array A must contain the lower
- triangular matrix and the strictly upper triangular part of
- A is not referenced.
- Note that when DIAG = 'U' or 'u', the diagonal elements of
- A are not referenced either, but are assumed to be unity.
- Unchanged on exit.
-
- LDA - INTEGER.
- On entry, LDA specifies the first dimension of A as declared
- in the calling (sub) program. When SIDE = 'L' or 'l' then
- LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
- then LDA must be at least max( 1, n ).
- Unchanged on exit.
-
- B - REAL array of DIMENSION ( LDB, n ).
- Before entry, the leading m by n part of the array B must
- contain the right-hand side matrix B, and on exit is
- overwritten by the solution matrix X.
-
- LDB - INTEGER.
- On entry, LDB specifies the first dimension of B as declared
- in the calling (sub) program. LDB must be at least
- max( 1, m ).
- Unchanged on exit.
-
-
- Level 3 Blas routine.
-
-
- -- Written on 8-February-1989.
- Jack Dongarra, Argonne National Laboratory.
- Iain Duff, AERE Harwell.
- Jeremy Du Croz, Numerical Algorithms Group Ltd.
- Sven Hammarling, Numerical Algorithms Group Ltd.
-
-
- Test the input parameters.
-*/
-
- /* Parameter adjustments */
- a_dim1 = *lda;
- a_offset = 1 + a_dim1;
- a -= a_offset;
- b_dim1 = *ldb;
- b_offset = 1 + b_dim1;
- b -= b_offset;
-
- /* Function Body */
- lside = lsame_(side, "L");
- if (lside) {
- nrowa = *m;
- } else {
- nrowa = *n;
- }
- nounit = lsame_(diag, "N");
- upper = lsame_(uplo, "U");
-
- info = 0;
- if (! lside && ! lsame_(side, "R")) {
- info = 1;
- } else if (! upper && ! lsame_(uplo, "L")) {
- info = 2;
- } else if (! lsame_(transa, "N") && ! lsame_(transa,
- "T") && ! lsame_(transa, "C")) {
- info = 3;
- } else if (! lsame_(diag, "U") && ! lsame_(diag,
- "N")) {
- info = 4;
- } else if (*m < 0) {
- info = 5;
- } else if (*n < 0) {
- info = 6;
- } else if (*lda < max(1,nrowa)) {
- info = 9;
- } else if (*ldb < max(1,*m)) {
- info = 11;
- }
- if (info != 0) {
- xerbla_("STRSM ", &info);
- return 0;
- }
-
-/* Quick return if possible. */
-
- if (*n == 0) {
- return 0;
- }
-
-/* And when alpha.eq.zero. */
-
- if (*alpha == 0.f) {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- b[i__ + j * b_dim1] = 0.f;
-/* L10: */
- }
-/* L20: */
- }
- return 0;
- }
-
-/* Start the operations. */
-
- if (lside) {
- if (lsame_(transa, "N")) {
-
-/* Form B := alpha*inv( A )*B. */
-
- if (upper) {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- if (*alpha != 1.f) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
- ;
-/* L30: */
- }
- }
- for (k = *m; k >= 1; --k) {
- if (b[k + j * b_dim1] != 0.f) {
- if (nounit) {
- b[k + j * b_dim1] /= a[k + k * a_dim1];
- }
- i__2 = k - 1;
- for (i__ = 1; i__ <= i__2; ++i__) {
- b[i__ + j * b_dim1] -= b[k + j * b_dim1] * a[
- i__ + k * a_dim1];
-/* L40: */
- }
- }
-/* L50: */
- }
-/* L60: */
- }
- } else {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- if (*alpha != 1.f) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
- ;
-/* L70: */
- }
- }
- i__2 = *m;
- for (k = 1; k <= i__2; ++k) {
- if (b[k + j * b_dim1] != 0.f) {
- if (nounit) {
- b[k + j * b_dim1] /= a[k + k * a_dim1];
- }
- i__3 = *m;
- for (i__ = k + 1; i__ <= i__3; ++i__) {
- b[i__ + j * b_dim1] -= b[k + j * b_dim1] * a[
- i__ + k * a_dim1];
-/* L80: */
- }
- }
-/* L90: */
- }
-/* L100: */
- }
- }
- } else {
-
-/* Form B := alpha*inv( A' )*B. */
-
- if (upper) {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- temp = *alpha * b[i__ + j * b_dim1];
- i__3 = i__ - 1;
- for (k = 1; k <= i__3; ++k) {
- temp -= a[k + i__ * a_dim1] * b[k + j * b_dim1];
-/* L110: */
- }
- if (nounit) {
- temp /= a[i__ + i__ * a_dim1];
- }
- b[i__ + j * b_dim1] = temp;
-/* L120: */
- }
-/* L130: */
- }
- } else {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- for (i__ = *m; i__ >= 1; --i__) {
- temp = *alpha * b[i__ + j * b_dim1];
- i__2 = *m;
- for (k = i__ + 1; k <= i__2; ++k) {
- temp -= a[k + i__ * a_dim1] * b[k + j * b_dim1];
-/* L140: */
- }
- if (nounit) {
- temp /= a[i__ + i__ * a_dim1];
- }
- b[i__ + j * b_dim1] = temp;
-/* L150: */
- }
-/* L160: */
- }
- }
- }
- } else {
- if (lsame_(transa, "N")) {
-
-/* Form B := alpha*B*inv( A ). */
-
- if (upper) {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- if (*alpha != 1.f) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
- ;
-/* L170: */
- }
- }
- i__2 = j - 1;
- for (k = 1; k <= i__2; ++k) {
- if (a[k + j * a_dim1] != 0.f) {
- i__3 = *m;
- for (i__ = 1; i__ <= i__3; ++i__) {
- b[i__ + j * b_dim1] -= a[k + j * a_dim1] * b[
- i__ + k * b_dim1];
-/* L180: */
- }
- }
-/* L190: */
- }
- if (nounit) {
- temp = 1.f / a[j + j * a_dim1];
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
-/* L200: */
- }
- }
-/* L210: */
- }
- } else {
- for (j = *n; j >= 1; --j) {
- if (*alpha != 1.f) {
- i__1 = *m;
- for (i__ = 1; i__ <= i__1; ++i__) {
- b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
- ;
-/* L220: */
- }
- }
- i__1 = *n;
- for (k = j + 1; k <= i__1; ++k) {
- if (a[k + j * a_dim1] != 0.f) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- b[i__ + j * b_dim1] -= a[k + j * a_dim1] * b[
- i__ + k * b_dim1];
-/* L230: */
- }
- }
-/* L240: */
- }
- if (nounit) {
- temp = 1.f / a[j + j * a_dim1];
- i__1 = *m;
- for (i__ = 1; i__ <= i__1; ++i__) {
- b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
-/* L250: */
- }
- }
-/* L260: */
- }
- }
- } else {
-
-/* Form B := alpha*B*inv( A' ). */
-
- if (upper) {
- for (k = *n; k >= 1; --k) {
- if (nounit) {
- temp = 1.f / a[k + k * a_dim1];
- i__1 = *m;
- for (i__ = 1; i__ <= i__1; ++i__) {
- b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
-/* L270: */
- }
- }
- i__1 = k - 1;
- for (j = 1; j <= i__1; ++j) {
- if (a[j + k * a_dim1] != 0.f) {
- temp = a[j + k * a_dim1];
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- b[i__ + j * b_dim1] -= temp * b[i__ + k *
- b_dim1];
-/* L280: */
- }
- }
-/* L290: */
- }
- if (*alpha != 1.f) {
- i__1 = *m;
- for (i__ = 1; i__ <= i__1; ++i__) {
- b[i__ + k * b_dim1] = *alpha * b[i__ + k * b_dim1]
- ;
-/* L300: */
- }
- }
-/* L310: */
- }
- } else {
- i__1 = *n;
- for (k = 1; k <= i__1; ++k) {
- if (nounit) {
- temp = 1.f / a[k + k * a_dim1];
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
-/* L320: */
- }
- }
- i__2 = *n;
- for (j = k + 1; j <= i__2; ++j) {
- if (a[j + k * a_dim1] != 0.f) {
- temp = a[j + k * a_dim1];
- i__3 = *m;
- for (i__ = 1; i__ <= i__3; ++i__) {
- b[i__ + j * b_dim1] -= temp * b[i__ + k *
- b_dim1];
-/* L330: */
- }
- }
-/* L340: */
- }
- if (*alpha != 1.f) {
- i__2 = *m;
- for (i__ = 1; i__ <= i__2; ++i__) {
- b[i__ + k * b_dim1] = *alpha * b[i__ + k * b_dim1]
- ;
-/* L350: */
- }
- }
-/* L360: */
- }
- }
- }
- }
-
- return 0;
-
-/* End of STRSM . */
-
-} /* strsm_ */
-
-/* Subroutine */ int xerbla_(char *srname, integer *info)
-{
- /* Format strings */
- static char fmt_9999[] = "(\002 ** On entry to \002,a6,\002 parameter nu"
- "mber \002,i2,\002 had \002,\002an illegal value\002)";
-
- /* Builtin functions */
- integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
- /* Subroutine */ int s_stop(char *, ftnlen);
-
- /* Fortran I/O blocks */
- static cilist io___60 = { 0, 6, 0, fmt_9999, 0 };
-
-
-/*
- -- LAPACK auxiliary routine (preliminary version) --
- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
- Courant Institute, Argonne National Lab, and Rice University
- February 29, 1992
-
-
- Purpose
- =======
-
- XERBLA is an error handler for the LAPACK routines.
- It is called by an LAPACK routine if an input parameter has an
- invalid value. A message is printed and execution stops.
-
- Installers may consider modifying the STOP statement in order to
- call system-specific exception-handling facilities.
-
- Arguments
- =========
-
- SRNAME (input) CHARACTER*6
- The name of the routine which called XERBLA.
-
- INFO (input) INTEGER
- The position of the invalid parameter in the parameter list
- of the calling routine.
-*/
-
-
- s_wsfe(&io___60);
- do_fio(&c__1, srname, (ftnlen)6);
- do_fio(&c__1, (char *)&(*info), (ftnlen)sizeof(integer));
- e_wsfe();
-
- s_stop("", (ftnlen)0);
-
-
-/* End of XERBLA */
-
- return 0;
-} /* xerbla_ */
-
generated by cgit v1.2.3 (git 2.26.2) at 2024-09-22 11:28:45 +0000