SphinxBase  0.6
src/libsphinxbase/util/blas_lite.c
00001 /*
00002 NOTE: This is generated code. Look in README.python for information on
00003       remaking this file.
00004 */
00005 #include "sphinxbase/f2c.h"
00006 
00007 #ifdef HAVE_CONFIG
00008 #include "config.h"
00009 #else
00010 extern doublereal slamch_(char *);
00011 #define EPSILON slamch_("Epsilon")
00012 #define SAFEMINIMUM slamch_("Safe minimum")
00013 #define PRECISION slamch_("Precision")
00014 #define BASE slamch_("Base")
00015 #endif
00016 
00017 
00018 extern doublereal slapy2_(real *, real *);
00019 
00020 
00021 
00022 /* Table of constant values */
00023 
00024 static integer c__1 = 1;
00025 
00026 logical lsame_(char *ca, char *cb)
00027 {
00028     /* System generated locals */
00029     logical ret_val;
00030 
00031     /* Local variables */
00032     static integer inta, intb, zcode;
00033 
00034 
00035 /*
00036     -- LAPACK auxiliary routine (version 3.0) --
00037        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
00038        Courant Institute, Argonne National Lab, and Rice University
00039        September 30, 1994
00040 
00041 
00042     Purpose
00043     =======
00044 
00045     LSAME returns .TRUE. if CA is the same letter as CB regardless of
00046     case.
00047 
00048     Arguments
00049     =========
00050 
00051     CA      (input) CHARACTER*1
00052     CB      (input) CHARACTER*1
00053             CA and CB specify the single characters to be compared.
00054 
00055    =====================================================================
00056 
00057 
00058        Test if the characters are equal
00059 */
00060 
00061     ret_val = *(unsigned char *)ca == *(unsigned char *)cb;
00062     if (ret_val) {
00063         return ret_val;
00064     }
00065 
00066 /*     Now test for equivalence if both characters are alphabetic. */
00067 
00068     zcode = 'Z';
00069 
00070 /*
00071        Use 'Z' rather than 'A' so that ASCII can be detected on Prime
00072        machines, on which ICHAR returns a value with bit 8 set.
00073        ICHAR('A') on Prime machines returns 193 which is the same as
00074        ICHAR('A') on an EBCDIC machine.
00075 */
00076 
00077     inta = *(unsigned char *)ca;
00078     intb = *(unsigned char *)cb;
00079 
00080     if (zcode == 90 || zcode == 122) {
00081 
00082 /*
00083           ASCII is assumed - ZCODE is the ASCII code of either lower or
00084           upper case 'Z'.
00085 */
00086 
00087         if (inta >= 97 && inta <= 122) {
00088             inta += -32;
00089         }
00090         if (intb >= 97 && intb <= 122) {
00091             intb += -32;
00092         }
00093 
00094     } else if (zcode == 233 || zcode == 169) {
00095 
00096 /*
00097           EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
00098           upper case 'Z'.
00099 */
00100 
00101         if (inta >= 129 && inta <= 137 || inta >= 145 && inta <= 153 || inta
00102                 >= 162 && inta <= 169) {
00103             inta += 64;
00104         }
00105         if (intb >= 129 && intb <= 137 || intb >= 145 && intb <= 153 || intb
00106                 >= 162 && intb <= 169) {
00107             intb += 64;
00108         }
00109 
00110     } else if (zcode == 218 || zcode == 250) {
00111 
00112 /*
00113           ASCII is assumed, on Prime machines - ZCODE is the ASCII code
00114           plus 128 of either lower or upper case 'Z'.
00115 */
00116 
00117         if (inta >= 225 && inta <= 250) {
00118             inta += -32;
00119         }
00120         if (intb >= 225 && intb <= 250) {
00121             intb += -32;
00122         }
00123     }
00124     ret_val = inta == intb;
00125 
00126 /*
00127        RETURN
00128 
00129        End of LSAME
00130 */
00131 
00132     return ret_val;
00133 } /* lsame_ */
00134 
00135 doublereal sdot_(integer *n, real *sx, integer *incx, real *sy, integer *incy)
00136 {
00137     /* System generated locals */
00138     integer i__1;
00139     real ret_val;
00140 
00141     /* Local variables */
00142     static integer i__, m, ix, iy, mp1;
00143     static real stemp;
00144 
00145 
00146 /*
00147        forms the dot product of two vectors.
00148        uses unrolled loops for increments equal to one.
00149        jack dongarra, linpack, 3/11/78.
00150        modified 12/3/93, array(1) declarations changed to array(*)
00151 */
00152 
00153 
00154     /* Parameter adjustments */
00155     --sy;
00156     --sx;
00157 
00158     /* Function Body */
00159     stemp = 0.f;
00160     ret_val = 0.f;
00161     if (*n <= 0) {
00162         return ret_val;
00163     }
00164     if (*incx == 1 && *incy == 1) {
00165         goto L20;
00166     }
00167 
00168 /*
00169           code for unequal increments or equal increments
00170             not equal to 1
00171 */
00172 
00173     ix = 1;
00174     iy = 1;
00175     if (*incx < 0) {
00176         ix = (-(*n) + 1) * *incx + 1;
00177     }
00178     if (*incy < 0) {
00179         iy = (-(*n) + 1) * *incy + 1;
00180     }
00181     i__1 = *n;
00182     for (i__ = 1; i__ <= i__1; ++i__) {
00183         stemp += sx[ix] * sy[iy];
00184         ix += *incx;
00185         iy += *incy;
00186 /* L10: */
00187     }
00188     ret_val = stemp;
00189     return ret_val;
00190 
00191 /*
00192           code for both increments equal to 1
00193 
00194 
00195           clean-up loop
00196 */
00197 
00198 L20:
00199     m = *n % 5;
00200     if (m == 0) {
00201         goto L40;
00202     }
00203     i__1 = m;
00204     for (i__ = 1; i__ <= i__1; ++i__) {
00205         stemp += sx[i__] * sy[i__];
00206 /* L30: */
00207     }
00208     if (*n < 5) {
00209         goto L60;
00210     }
00211 L40:
00212     mp1 = m + 1;
00213     i__1 = *n;
00214     for (i__ = mp1; i__ <= i__1; i__ += 5) {
00215         stemp = stemp + sx[i__] * sy[i__] + sx[i__ + 1] * sy[i__ + 1] + sx[
00216                 i__ + 2] * sy[i__ + 2] + sx[i__ + 3] * sy[i__ + 3] + sx[i__ +
00217                 4] * sy[i__ + 4];
00218 /* L50: */
00219     }
00220 L60:
00221     ret_val = stemp;
00222     return ret_val;
00223 } /* sdot_ */
00224 
00225 /* Subroutine */ int sgemm_(char *transa, char *transb, integer *m, integer *
00226         n, integer *k, real *alpha, real *a, integer *lda, real *b, integer *
00227         ldb, real *beta, real *c__, integer *ldc)
00228 {
00229     /* System generated locals */
00230     integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
00231             i__3;
00232 
00233     /* Local variables */
00234     static integer i__, j, l, info;
00235     static logical nota, notb;
00236     static real temp;
00237     static integer ncola;
00238     extern logical lsame_(char *, char *);
00239     static integer nrowa, nrowb;
00240     extern /* Subroutine */ int xerbla_(char *, integer *);
00241 
00242 
00243 /*
00244     Purpose
00245     =======
00246 
00247     SGEMM  performs one of the matrix-matrix operations
00248 
00249        C := alpha*op( A )*op( B ) + beta*C,
00250 
00251     where  op( X ) is one of
00252 
00253        op( X ) = X   or   op( X ) = X',
00254 
00255     alpha and beta are scalars, and A, B and C are matrices, with op( A )
00256     an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
00257 
00258     Parameters
00259     ==========
00260 
00261     TRANSA - CHARACTER*1.
00262              On entry, TRANSA specifies the form of op( A ) to be used in
00263              the matrix multiplication as follows:
00264 
00265                 TRANSA = 'N' or 'n',  op( A ) = A.
00266 
00267                 TRANSA = 'T' or 't',  op( A ) = A'.
00268 
00269                 TRANSA = 'C' or 'c',  op( A ) = A'.
00270 
00271              Unchanged on exit.
00272 
00273     TRANSB - CHARACTER*1.
00274              On entry, TRANSB specifies the form of op( B ) to be used in
00275              the matrix multiplication as follows:
00276 
00277                 TRANSB = 'N' or 'n',  op( B ) = B.
00278 
00279                 TRANSB = 'T' or 't',  op( B ) = B'.
00280 
00281                 TRANSB = 'C' or 'c',  op( B ) = B'.
00282 
00283              Unchanged on exit.
00284 
00285     M      - INTEGER.
00286              On entry,  M  specifies  the number  of rows  of the  matrix
00287              op( A )  and of the  matrix  C.  M  must  be at least  zero.
00288              Unchanged on exit.
00289 
00290     N      - INTEGER.
00291              On entry,  N  specifies the number  of columns of the matrix
00292              op( B ) and the number of columns of the matrix C. N must be
00293              at least zero.
00294              Unchanged on exit.
00295 
00296     K      - INTEGER.
00297              On entry,  K  specifies  the number of columns of the matrix
00298              op( A ) and the number of rows of the matrix op( B ). K must
00299              be at least  zero.
00300              Unchanged on exit.
00301 
00302     ALPHA  - REAL            .
00303              On entry, ALPHA specifies the scalar alpha.
00304              Unchanged on exit.
00305 
00306     A      - REAL             array of DIMENSION ( LDA, ka ), where ka is
00307              k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
00308              Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
00309              part of the array  A  must contain the matrix  A,  otherwise
00310              the leading  k by m  part of the array  A  must contain  the
00311              matrix A.
00312              Unchanged on exit.
00313 
00314     LDA    - INTEGER.
00315              On entry, LDA specifies the first dimension of A as declared
00316              in the calling (sub) program. When  TRANSA = 'N' or 'n' then
00317              LDA must be at least  max( 1, m ), otherwise  LDA must be at
00318              least  max( 1, k ).
00319              Unchanged on exit.
00320 
00321     B      - REAL             array of DIMENSION ( LDB, kb ), where kb is
00322              n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
00323              Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
00324              part of the array  B  must contain the matrix  B,  otherwise
00325              the leading  n by k  part of the array  B  must contain  the
00326              matrix B.
00327              Unchanged on exit.
00328 
00329     LDB    - INTEGER.
00330              On entry, LDB specifies the first dimension of B as declared
00331              in the calling (sub) program. When  TRANSB = 'N' or 'n' then
00332              LDB must be at least  max( 1, k ), otherwise  LDB must be at
00333              least  max( 1, n ).
00334              Unchanged on exit.
00335 
00336     BETA   - REAL            .
00337              On entry,  BETA  specifies the scalar  beta.  When  BETA  is
00338              supplied as zero then C need not be set on input.
00339              Unchanged on exit.
00340 
00341     C      - REAL             array of DIMENSION ( LDC, n ).
00342              Before entry, the leading  m by n  part of the array  C must
00343              contain the matrix  C,  except when  beta  is zero, in which
00344              case C need not be set on entry.
00345              On exit, the array  C  is overwritten by the  m by n  matrix
00346              ( alpha*op( A )*op( B ) + beta*C ).
00347 
00348     LDC    - INTEGER.
00349              On entry, LDC specifies the first dimension of C as declared
00350              in  the  calling  (sub)  program.   LDC  must  be  at  least
00351              max( 1, m ).
00352              Unchanged on exit.
00353 
00354 
00355     Level 3 Blas routine.
00356 
00357     -- Written on 8-February-1989.
00358        Jack Dongarra, Argonne National Laboratory.
00359        Iain Duff, AERE Harwell.
00360        Jeremy Du Croz, Numerical Algorithms Group Ltd.
00361        Sven Hammarling, Numerical Algorithms Group Ltd.
00362 
00363 
00364        Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
00365        transposed and set  NROWA, NCOLA and  NROWB  as the number of rows
00366        and  columns of  A  and the  number of  rows  of  B  respectively.
00367 */
00368 
00369     /* Parameter adjustments */
00370     a_dim1 = *lda;
00371     a_offset = 1 + a_dim1;
00372     a -= a_offset;
00373     b_dim1 = *ldb;
00374     b_offset = 1 + b_dim1;
00375     b -= b_offset;
00376     c_dim1 = *ldc;
00377     c_offset = 1 + c_dim1;
00378     c__ -= c_offset;
00379 
00380     /* Function Body */
00381     nota = lsame_(transa, "N");
00382     notb = lsame_(transb, "N");
00383     if (nota) {
00384         nrowa = *m;
00385         ncola = *k;
00386     } else {
00387         nrowa = *k;
00388         ncola = *m;
00389     }
00390     if (notb) {
00391         nrowb = *k;
00392     } else {
00393         nrowb = *n;
00394     }
00395 
00396 /*     Test the input parameters. */
00397 
00398     info = 0;
00399     if (! nota && ! lsame_(transa, "C") && ! lsame_(
00400             transa, "T")) {
00401         info = 1;
00402     } else if (! notb && ! lsame_(transb, "C") && !
00403             lsame_(transb, "T")) {
00404         info = 2;
00405     } else if (*m < 0) {
00406         info = 3;
00407     } else if (*n < 0) {
00408         info = 4;
00409     } else if (*k < 0) {
00410         info = 5;
00411     } else if (*lda < max(1,nrowa)) {
00412         info = 8;
00413     } else if (*ldb < max(1,nrowb)) {
00414         info = 10;
00415     } else if (*ldc < max(1,*m)) {
00416         info = 13;
00417     }
00418     if (info != 0) {
00419         xerbla_("SGEMM ", &info);
00420         return 0;
00421     }
00422 
00423 /*     Quick return if possible. */
00424 
00425     if (*m == 0 || *n == 0 || (*alpha == 0.f || *k == 0) && *beta == 1.f) {
00426         return 0;
00427     }
00428 
00429 /*     And if  alpha.eq.zero. */
00430 
00431     if (*alpha == 0.f) {
00432         if (*beta == 0.f) {
00433             i__1 = *n;
00434             for (j = 1; j <= i__1; ++j) {
00435                 i__2 = *m;
00436                 for (i__ = 1; i__ <= i__2; ++i__) {
00437                     c__[i__ + j * c_dim1] = 0.f;
00438 /* L10: */
00439                 }
00440 /* L20: */
00441             }
00442         } else {
00443             i__1 = *n;
00444             for (j = 1; j <= i__1; ++j) {
00445                 i__2 = *m;
00446                 for (i__ = 1; i__ <= i__2; ++i__) {
00447                     c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
00448 /* L30: */
00449                 }
00450 /* L40: */
00451             }
00452         }
00453         return 0;
00454     }
00455 
00456 /*     Start the operations. */
00457 
00458     if (notb) {
00459         if (nota) {
00460 
00461 /*           Form  C := alpha*A*B + beta*C. */
00462 
00463             i__1 = *n;
00464             for (j = 1; j <= i__1; ++j) {
00465                 if (*beta == 0.f) {
00466                     i__2 = *m;
00467                     for (i__ = 1; i__ <= i__2; ++i__) {
00468                         c__[i__ + j * c_dim1] = 0.f;
00469 /* L50: */
00470                     }
00471                 } else if (*beta != 1.f) {
00472                     i__2 = *m;
00473                     for (i__ = 1; i__ <= i__2; ++i__) {
00474                         c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
00475 /* L60: */
00476                     }
00477                 }
00478                 i__2 = *k;
00479                 for (l = 1; l <= i__2; ++l) {
00480                     if (b[l + j * b_dim1] != 0.f) {
00481                         temp = *alpha * b[l + j * b_dim1];
00482                         i__3 = *m;
00483                         for (i__ = 1; i__ <= i__3; ++i__) {
00484                             c__[i__ + j * c_dim1] += temp * a[i__ + l *
00485                                     a_dim1];
00486 /* L70: */
00487                         }
00488                     }
00489 /* L80: */
00490                 }
00491 /* L90: */
00492             }
00493         } else {
00494 
00495 /*           Form  C := alpha*A'*B + beta*C */
00496 
00497             i__1 = *n;
00498             for (j = 1; j <= i__1; ++j) {
00499                 i__2 = *m;
00500                 for (i__ = 1; i__ <= i__2; ++i__) {
00501                     temp = 0.f;
00502                     i__3 = *k;
00503                     for (l = 1; l <= i__3; ++l) {
00504                         temp += a[l + i__ * a_dim1] * b[l + j * b_dim1];
00505 /* L100: */
00506                     }
00507                     if (*beta == 0.f) {
00508                         c__[i__ + j * c_dim1] = *alpha * temp;
00509                     } else {
00510                         c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
00511                                 i__ + j * c_dim1];
00512                     }
00513 /* L110: */
00514                 }
00515 /* L120: */
00516             }
00517         }
00518     } else {
00519         if (nota) {
00520 
00521 /*           Form  C := alpha*A*B' + beta*C */
00522 
00523             i__1 = *n;
00524             for (j = 1; j <= i__1; ++j) {
00525                 if (*beta == 0.f) {
00526                     i__2 = *m;
00527                     for (i__ = 1; i__ <= i__2; ++i__) {
00528                         c__[i__ + j * c_dim1] = 0.f;
00529 /* L130: */
00530                     }
00531                 } else if (*beta != 1.f) {
00532                     i__2 = *m;
00533                     for (i__ = 1; i__ <= i__2; ++i__) {
00534                         c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
00535 /* L140: */
00536                     }
00537                 }
00538                 i__2 = *k;
00539                 for (l = 1; l <= i__2; ++l) {
00540                     if (b[j + l * b_dim1] != 0.f) {
00541                         temp = *alpha * b[j + l * b_dim1];
00542                         i__3 = *m;
00543                         for (i__ = 1; i__ <= i__3; ++i__) {
00544                             c__[i__ + j * c_dim1] += temp * a[i__ + l *
00545                                     a_dim1];
00546 /* L150: */
00547                         }
00548                     }
00549 /* L160: */
00550                 }
00551 /* L170: */
00552             }
00553         } else {
00554 
00555 /*           Form  C := alpha*A'*B' + beta*C */
00556 
00557             i__1 = *n;
00558             for (j = 1; j <= i__1; ++j) {
00559                 i__2 = *m;
00560                 for (i__ = 1; i__ <= i__2; ++i__) {
00561                     temp = 0.f;
00562                     i__3 = *k;
00563                     for (l = 1; l <= i__3; ++l) {
00564                         temp += a[l + i__ * a_dim1] * b[j + l * b_dim1];
00565 /* L180: */
00566                     }
00567                     if (*beta == 0.f) {
00568                         c__[i__ + j * c_dim1] = *alpha * temp;
00569                     } else {
00570                         c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
00571                                 i__ + j * c_dim1];
00572                     }
00573 /* L190: */
00574                 }
00575 /* L200: */
00576             }
00577         }
00578     }
00579 
00580     return 0;
00581 
00582 /*     End of SGEMM . */
00583 
00584 } /* sgemm_ */
00585 
00586 /* Subroutine */ int sgemv_(char *trans, integer *m, integer *n, real *alpha,
00587         real *a, integer *lda, real *x, integer *incx, real *beta, real *y,
00588         integer *incy)
00589 {
00590     /* System generated locals */
00591     integer a_dim1, a_offset, i__1, i__2;
00592 
00593     /* Local variables */
00594     static integer i__, j, ix, iy, jx, jy, kx, ky, info;
00595     static real temp;
00596     static integer lenx, leny;
00597     extern logical lsame_(char *, char *);
00598     extern /* Subroutine */ int xerbla_(char *, integer *);
00599 
00600 
00601 /*
00602     Purpose
00603     =======
00604 
00605     SGEMV  performs one of the matrix-vector operations
00606 
00607        y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
00608 
00609     where alpha and beta are scalars, x and y are vectors and A is an
00610     m by n matrix.
00611 
00612     Parameters
00613     ==========
00614 
00615     TRANS  - CHARACTER*1.
00616              On entry, TRANS specifies the operation to be performed as
00617              follows:
00618 
00619                 TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
00620 
00621                 TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
00622 
00623                 TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.
00624 
00625              Unchanged on exit.
00626 
00627     M      - INTEGER.
00628              On entry, M specifies the number of rows of the matrix A.
00629              M must be at least zero.
00630              Unchanged on exit.
00631 
00632     N      - INTEGER.
00633              On entry, N specifies the number of columns of the matrix A.
00634              N must be at least zero.
00635              Unchanged on exit.
00636 
00637     ALPHA  - REAL            .
00638              On entry, ALPHA specifies the scalar alpha.
00639              Unchanged on exit.
00640 
00641     A      - REAL             array of DIMENSION ( LDA, n ).
00642              Before entry, the leading m by n part of the array A must
00643              contain the matrix of coefficients.
00644              Unchanged on exit.
00645 
00646     LDA    - INTEGER.
00647              On entry, LDA specifies the first dimension of A as declared
00648              in the calling (sub) program. LDA must be at least
00649              max( 1, m ).
00650              Unchanged on exit.
00651 
00652     X      - REAL             array of DIMENSION at least
00653              ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
00654              and at least
00655              ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
00656              Before entry, the incremented array X must contain the
00657              vector x.
00658              Unchanged on exit.
00659 
00660     INCX   - INTEGER.
00661              On entry, INCX specifies the increment for the elements of
00662              X. INCX must not be zero.
00663              Unchanged on exit.
00664 
00665     BETA   - REAL            .
00666              On entry, BETA specifies the scalar beta. When BETA is
00667              supplied as zero then Y need not be set on input.
00668              Unchanged on exit.
00669 
00670     Y      - REAL             array of DIMENSION at least
00671              ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
00672              and at least
00673              ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
00674              Before entry with BETA non-zero, the incremented array Y
00675              must contain the vector y. On exit, Y is overwritten by the
00676              updated vector y.
00677 
00678     INCY   - INTEGER.
00679              On entry, INCY specifies the increment for the elements of
00680              Y. INCY must not be zero.
00681              Unchanged on exit.
00682 
00683 
00684     Level 2 Blas routine.
00685 
00686     -- Written on 22-October-1986.
00687        Jack Dongarra, Argonne National Lab.
00688        Jeremy Du Croz, Nag Central Office.
00689        Sven Hammarling, Nag Central Office.
00690        Richard Hanson, Sandia National Labs.
00691 
00692 
00693        Test the input parameters.
00694 */
00695 
00696     /* Parameter adjustments */
00697     a_dim1 = *lda;
00698     a_offset = 1 + a_dim1;
00699     a -= a_offset;
00700     --x;
00701     --y;
00702 
00703     /* Function Body */
00704     info = 0;
00705     if (! lsame_(trans, "N") && ! lsame_(trans, "T") && ! lsame_(trans, "C")
00706             ) {
00707         info = 1;
00708     } else if (*m < 0) {
00709         info = 2;
00710     } else if (*n < 0) {
00711         info = 3;
00712     } else if (*lda < max(1,*m)) {
00713         info = 6;
00714     } else if (*incx == 0) {
00715         info = 8;
00716     } else if (*incy == 0) {
00717         info = 11;
00718     }
00719     if (info != 0) {
00720         xerbla_("SGEMV ", &info);
00721         return 0;
00722     }
00723 
00724 /*     Quick return if possible. */
00725 
00726     if (*m == 0 || *n == 0 || *alpha == 0.f && *beta == 1.f) {
00727         return 0;
00728     }
00729 
00730 /*
00731        Set  LENX  and  LENY, the lengths of the vectors x and y, and set
00732        up the start points in  X  and  Y.
00733 */
00734 
00735     if (lsame_(trans, "N")) {
00736         lenx = *n;
00737         leny = *m;
00738     } else {
00739         lenx = *m;
00740         leny = *n;
00741     }
00742     if (*incx > 0) {
00743         kx = 1;
00744     } else {
00745         kx = 1 - (lenx - 1) * *incx;
00746     }
00747     if (*incy > 0) {
00748         ky = 1;
00749     } else {
00750         ky = 1 - (leny - 1) * *incy;
00751     }
00752 
00753 /*
00754        Start the operations. In this version the elements of A are
00755        accessed sequentially with one pass through A.
00756 
00757        First form  y := beta*y.
00758 */
00759 
00760     if (*beta != 1.f) {
00761         if (*incy == 1) {
00762             if (*beta == 0.f) {
00763                 i__1 = leny;
00764                 for (i__ = 1; i__ <= i__1; ++i__) {
00765                     y[i__] = 0.f;
00766 /* L10: */
00767                 }
00768             } else {
00769                 i__1 = leny;
00770                 for (i__ = 1; i__ <= i__1; ++i__) {
00771                     y[i__] = *beta * y[i__];
00772 /* L20: */
00773                 }
00774             }
00775         } else {
00776             iy = ky;
00777             if (*beta == 0.f) {
00778                 i__1 = leny;
00779                 for (i__ = 1; i__ <= i__1; ++i__) {
00780                     y[iy] = 0.f;
00781                     iy += *incy;
00782 /* L30: */
00783                 }
00784             } else {
00785                 i__1 = leny;
00786                 for (i__ = 1; i__ <= i__1; ++i__) {
00787                     y[iy] = *beta * y[iy];
00788                     iy += *incy;
00789 /* L40: */
00790                 }
00791             }
00792         }
00793     }
00794     if (*alpha == 0.f) {
00795         return 0;
00796     }
00797     if (lsame_(trans, "N")) {
00798 
00799 /*        Form  y := alpha*A*x + y. */
00800 
00801         jx = kx;
00802         if (*incy == 1) {
00803             i__1 = *n;
00804             for (j = 1; j <= i__1; ++j) {
00805                 if (x[jx] != 0.f) {
00806                     temp = *alpha * x[jx];
00807                     i__2 = *m;
00808                     for (i__ = 1; i__ <= i__2; ++i__) {
00809                         y[i__] += temp * a[i__ + j * a_dim1];
00810 /* L50: */
00811                     }
00812                 }
00813                 jx += *incx;
00814 /* L60: */
00815             }
00816         } else {
00817             i__1 = *n;
00818             for (j = 1; j <= i__1; ++j) {
00819                 if (x[jx] != 0.f) {
00820                     temp = *alpha * x[jx];
00821                     iy = ky;
00822                     i__2 = *m;
00823                     for (i__ = 1; i__ <= i__2; ++i__) {
00824                         y[iy] += temp * a[i__ + j * a_dim1];
00825                         iy += *incy;
00826 /* L70: */
00827                     }
00828                 }
00829                 jx += *incx;
00830 /* L80: */
00831             }
00832         }
00833     } else {
00834 
00835 /*        Form  y := alpha*A'*x + y. */
00836 
00837         jy = ky;
00838         if (*incx == 1) {
00839             i__1 = *n;
00840             for (j = 1; j <= i__1; ++j) {
00841                 temp = 0.f;
00842                 i__2 = *m;
00843                 for (i__ = 1; i__ <= i__2; ++i__) {
00844                     temp += a[i__ + j * a_dim1] * x[i__];
00845 /* L90: */
00846                 }
00847                 y[jy] += *alpha * temp;
00848                 jy += *incy;
00849 /* L100: */
00850             }
00851         } else {
00852             i__1 = *n;
00853             for (j = 1; j <= i__1; ++j) {
00854                 temp = 0.f;
00855                 ix = kx;
00856                 i__2 = *m;
00857                 for (i__ = 1; i__ <= i__2; ++i__) {
00858                     temp += a[i__ + j * a_dim1] * x[ix];
00859                     ix += *incx;
00860 /* L110: */
00861                 }
00862                 y[jy] += *alpha * temp;
00863                 jy += *incy;
00864 /* L120: */
00865             }
00866         }
00867     }
00868 
00869     return 0;
00870 
00871 /*     End of SGEMV . */
00872 
00873 } /* sgemv_ */
00874 
00875 /* Subroutine */ int sscal_(integer *n, real *sa, real *sx, integer *incx)
00876 {
00877     /* System generated locals */
00878     integer i__1, i__2;
00879 
00880     /* Local variables */
00881     static integer i__, m, mp1, nincx;
00882 
00883 
00884 /*
00885        scales a vector by a constant.
00886        uses unrolled loops for increment equal to 1.
00887        jack dongarra, linpack, 3/11/78.
00888        modified 3/93 to return if incx .le. 0.
00889        modified 12/3/93, array(1) declarations changed to array(*)
00890 */
00891 
00892 
00893     /* Parameter adjustments */
00894     --sx;
00895 
00896     /* Function Body */
00897     if (*n <= 0 || *incx <= 0) {
00898         return 0;
00899     }
00900     if (*incx == 1) {
00901         goto L20;
00902     }
00903 
00904 /*        code for increment not equal to 1 */
00905 
00906     nincx = *n * *incx;
00907     i__1 = nincx;
00908     i__2 = *incx;
00909     for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
00910         sx[i__] = *sa * sx[i__];
00911 /* L10: */
00912     }
00913     return 0;
00914 
00915 /*
00916           code for increment equal to 1
00917 
00918 
00919           clean-up loop
00920 */
00921 
00922 L20:
00923     m = *n % 5;
00924     if (m == 0) {
00925         goto L40;
00926     }
00927     i__2 = m;
00928     for (i__ = 1; i__ <= i__2; ++i__) {
00929         sx[i__] = *sa * sx[i__];
00930 /* L30: */
00931     }
00932     if (*n < 5) {
00933         return 0;
00934     }
00935 L40:
00936     mp1 = m + 1;
00937     i__2 = *n;
00938     for (i__ = mp1; i__ <= i__2; i__ += 5) {
00939         sx[i__] = *sa * sx[i__];
00940         sx[i__ + 1] = *sa * sx[i__ + 1];
00941         sx[i__ + 2] = *sa * sx[i__ + 2];
00942         sx[i__ + 3] = *sa * sx[i__ + 3];
00943         sx[i__ + 4] = *sa * sx[i__ + 4];
00944 /* L50: */
00945     }
00946     return 0;
00947 } /* sscal_ */
00948 
00949 /* Subroutine */ int ssymm_(char *side, char *uplo, integer *m, integer *n,
00950         real *alpha, real *a, integer *lda, real *b, integer *ldb, real *beta,
00951          real *c__, integer *ldc)
00952 {
00953     /* System generated locals */
00954     integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
00955             i__3;
00956 
00957     /* Local variables */
00958     static integer i__, j, k, info;
00959     static real temp1, temp2;
00960     extern logical lsame_(char *, char *);
00961     static integer nrowa;
00962     static logical upper;
00963     extern /* Subroutine */ int xerbla_(char *, integer *);
00964 
00965 
00966 /*
00967     Purpose
00968     =======
00969 
00970     SSYMM  performs one of the matrix-matrix operations
00971 
00972        C := alpha*A*B + beta*C,
00973 
00974     or
00975 
00976        C := alpha*B*A + beta*C,
00977 
00978     where alpha and beta are scalars,  A is a symmetric matrix and  B and
00979     C are  m by n matrices.
00980 
00981     Parameters
00982     ==========
00983 
00984     SIDE   - CHARACTER*1.
00985              On entry,  SIDE  specifies whether  the  symmetric matrix  A
00986              appears on the  left or right  in the  operation as follows:
00987 
00988                 SIDE = 'L' or 'l'   C := alpha*A*B + beta*C,
00989 
00990                 SIDE = 'R' or 'r'   C := alpha*B*A + beta*C,
00991 
00992              Unchanged on exit.
00993 
00994     UPLO   - CHARACTER*1.
00995              On  entry,   UPLO  specifies  whether  the  upper  or  lower
00996              triangular  part  of  the  symmetric  matrix   A  is  to  be
00997              referenced as follows:
00998 
00999                 UPLO = 'U' or 'u'   Only the upper triangular part of the
01000                                     symmetric matrix is to be referenced.
01001 
01002                 UPLO = 'L' or 'l'   Only the lower triangular part of the
01003                                     symmetric matrix is to be referenced.
01004 
01005              Unchanged on exit.
01006 
01007     M      - INTEGER.
01008              On entry,  M  specifies the number of rows of the matrix  C.
01009              M  must be at least zero.
01010              Unchanged on exit.
01011 
01012     N      - INTEGER.
01013              On entry, N specifies the number of columns of the matrix C.
01014              N  must be at least zero.
01015              Unchanged on exit.
01016 
01017     ALPHA  - REAL            .
01018              On entry, ALPHA specifies the scalar alpha.
01019              Unchanged on exit.
01020 
01021     A      - REAL             array of DIMENSION ( LDA, ka ), where ka is
01022              m  when  SIDE = 'L' or 'l'  and is  n otherwise.
01023              Before entry  with  SIDE = 'L' or 'l',  the  m by m  part of
01024              the array  A  must contain the  symmetric matrix,  such that
01025              when  UPLO = 'U' or 'u', the leading m by m upper triangular
01026              part of the array  A  must contain the upper triangular part
01027              of the  symmetric matrix and the  strictly  lower triangular
01028              part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
01029              the leading  m by m  lower triangular part  of the  array  A
01030              must  contain  the  lower triangular part  of the  symmetric
01031              matrix and the  strictly upper triangular part of  A  is not
01032              referenced.
01033              Before entry  with  SIDE = 'R' or 'r',  the  n by n  part of
01034              the array  A  must contain the  symmetric matrix,  such that
01035              when  UPLO = 'U' or 'u', the leading n by n upper triangular
01036              part of the array  A  must contain the upper triangular part
01037              of the  symmetric matrix and the  strictly  lower triangular
01038              part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
01039              the leading  n by n  lower triangular part  of the  array  A
01040              must  contain  the  lower triangular part  of the  symmetric
01041              matrix and the  strictly upper triangular part of  A  is not
01042              referenced.
01043              Unchanged on exit.
01044 
01045     LDA    - INTEGER.
01046              On entry, LDA specifies the first dimension of A as declared
01047              in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
01048              LDA must be at least  max( 1, m ), otherwise  LDA must be at
01049              least  max( 1, n ).
01050              Unchanged on exit.
01051 
01052     B      - REAL             array of DIMENSION ( LDB, n ).
01053              Before entry, the leading  m by n part of the array  B  must
01054              contain the matrix B.
01055              Unchanged on exit.
01056 
01057     LDB    - INTEGER.
01058              On entry, LDB specifies the first dimension of B as declared
01059              in  the  calling  (sub)  program.   LDB  must  be  at  least
01060              max( 1, m ).
01061              Unchanged on exit.
01062 
01063     BETA   - REAL            .
01064              On entry,  BETA  specifies the scalar  beta.  When  BETA  is
01065              supplied as zero then C need not be set on input.
01066              Unchanged on exit.
01067 
01068     C      - REAL             array of DIMENSION ( LDC, n ).
01069              Before entry, the leading  m by n  part of the array  C must
01070              contain the matrix  C,  except when  beta  is zero, in which
01071              case C need not be set on entry.
01072              On exit, the array  C  is overwritten by the  m by n updated
01073              matrix.
01074 
01075     LDC    - INTEGER.
01076              On entry, LDC specifies the first dimension of C as declared
01077              in  the  calling  (sub)  program.   LDC  must  be  at  least
01078              max( 1, m ).
01079              Unchanged on exit.
01080 
01081 
01082     Level 3 Blas routine.
01083 
01084     -- Written on 8-February-1989.
01085        Jack Dongarra, Argonne National Laboratory.
01086        Iain Duff, AERE Harwell.
01087        Jeremy Du Croz, Numerical Algorithms Group Ltd.
01088        Sven Hammarling, Numerical Algorithms Group Ltd.
01089 
01090 
01091        Set NROWA as the number of rows of A.
01092 */
01093 
01094     /* Parameter adjustments */
01095     a_dim1 = *lda;
01096     a_offset = 1 + a_dim1;
01097     a -= a_offset;
01098     b_dim1 = *ldb;
01099     b_offset = 1 + b_dim1;
01100     b -= b_offset;
01101     c_dim1 = *ldc;
01102     c_offset = 1 + c_dim1;
01103     c__ -= c_offset;
01104 
01105     /* Function Body */
01106     if (lsame_(side, "L")) {
01107         nrowa = *m;
01108     } else {
01109         nrowa = *n;
01110     }
01111     upper = lsame_(uplo, "U");
01112 
01113 /*     Test the input parameters. */
01114 
01115     info = 0;
01116     if (! lsame_(side, "L") && ! lsame_(side, "R")) {
01117         info = 1;
01118     } else if (! upper && ! lsame_(uplo, "L")) {
01119         info = 2;
01120     } else if (*m < 0) {
01121         info = 3;
01122     } else if (*n < 0) {
01123         info = 4;
01124     } else if (*lda < max(1,nrowa)) {
01125         info = 7;
01126     } else if (*ldb < max(1,*m)) {
01127         info = 9;
01128     } else if (*ldc < max(1,*m)) {
01129         info = 12;
01130     }
01131     if (info != 0) {
01132         xerbla_("SSYMM ", &info);
01133         return 0;
01134     }
01135 
01136 /*     Quick return if possible. */
01137 
01138     if (*m == 0 || *n == 0 || *alpha == 0.f && *beta == 1.f) {
01139         return 0;
01140     }
01141 
01142 /*     And when  alpha.eq.zero. */
01143 
01144     if (*alpha == 0.f) {
01145         if (*beta == 0.f) {
01146             i__1 = *n;
01147             for (j = 1; j <= i__1; ++j) {
01148                 i__2 = *m;
01149                 for (i__ = 1; i__ <= i__2; ++i__) {
01150                     c__[i__ + j * c_dim1] = 0.f;
01151 /* L10: */
01152                 }
01153 /* L20: */
01154             }
01155         } else {
01156             i__1 = *n;
01157             for (j = 1; j <= i__1; ++j) {
01158                 i__2 = *m;
01159                 for (i__ = 1; i__ <= i__2; ++i__) {
01160                     c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
01161 /* L30: */
01162                 }
01163 /* L40: */
01164             }
01165         }
01166         return 0;
01167     }
01168 
01169 /*     Start the operations. */
01170 
01171     if (lsame_(side, "L")) {
01172 
01173 /*        Form  C := alpha*A*B + beta*C. */
01174 
01175         if (upper) {
01176             i__1 = *n;
01177             for (j = 1; j <= i__1; ++j) {
01178                 i__2 = *m;
01179                 for (i__ = 1; i__ <= i__2; ++i__) {
01180                     temp1 = *alpha * b[i__ + j * b_dim1];
01181                     temp2 = 0.f;
01182                     i__3 = i__ - 1;
01183                     for (k = 1; k <= i__3; ++k) {
01184                         c__[k + j * c_dim1] += temp1 * a[k + i__ * a_dim1];
01185                         temp2 += b[k + j * b_dim1] * a[k + i__ * a_dim1];
01186 /* L50: */
01187                     }
01188                     if (*beta == 0.f) {
01189                         c__[i__ + j * c_dim1] = temp1 * a[i__ + i__ * a_dim1]
01190                                 + *alpha * temp2;
01191                     } else {
01192                         c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1]
01193                                 + temp1 * a[i__ + i__ * a_dim1] + *alpha *
01194                                 temp2;
01195                     }
01196 /* L60: */
01197                 }
01198 /* L70: */
01199             }
01200         } else {
01201             i__1 = *n;
01202             for (j = 1; j <= i__1; ++j) {
01203                 for (i__ = *m; i__ >= 1; --i__) {
01204                     temp1 = *alpha * b[i__ + j * b_dim1];
01205                     temp2 = 0.f;
01206                     i__2 = *m;
01207                     for (k = i__ + 1; k <= i__2; ++k) {
01208                         c__[k + j * c_dim1] += temp1 * a[k + i__ * a_dim1];
01209                         temp2 += b[k + j * b_dim1] * a[k + i__ * a_dim1];
01210 /* L80: */
01211                     }
01212                     if (*beta == 0.f) {
01213                         c__[i__ + j * c_dim1] = temp1 * a[i__ + i__ * a_dim1]
01214                                 + *alpha * temp2;
01215                     } else {
01216                         c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1]
01217                                 + temp1 * a[i__ + i__ * a_dim1] + *alpha *
01218                                 temp2;
01219                     }
01220 /* L90: */
01221                 }
01222 /* L100: */
01223             }
01224         }
01225     } else {
01226 
01227 /*        Form  C := alpha*B*A + beta*C. */
01228 
01229         i__1 = *n;
01230         for (j = 1; j <= i__1; ++j) {
01231             temp1 = *alpha * a[j + j * a_dim1];
01232             if (*beta == 0.f) {
01233                 i__2 = *m;
01234                 for (i__ = 1; i__ <= i__2; ++i__) {
01235                     c__[i__ + j * c_dim1] = temp1 * b[i__ + j * b_dim1];
01236 /* L110: */
01237                 }
01238             } else {
01239                 i__2 = *m;
01240                 for (i__ = 1; i__ <= i__2; ++i__) {
01241                     c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1] +
01242                             temp1 * b[i__ + j * b_dim1];
01243 /* L120: */
01244                 }
01245             }
01246             i__2 = j - 1;
01247             for (k = 1; k <= i__2; ++k) {
01248                 if (upper) {
01249                     temp1 = *alpha * a[k + j * a_dim1];
01250                 } else {
01251                     temp1 = *alpha * a[j + k * a_dim1];
01252                 }
01253                 i__3 = *m;
01254                 for (i__ = 1; i__ <= i__3; ++i__) {
01255                     c__[i__ + j * c_dim1] += temp1 * b[i__ + k * b_dim1];
01256 /* L130: */
01257                 }
01258 /* L140: */
01259             }
01260             i__2 = *n;
01261             for (k = j + 1; k <= i__2; ++k) {
01262                 if (upper) {
01263                     temp1 = *alpha * a[j + k * a_dim1];
01264                 } else {
01265                     temp1 = *alpha * a[k + j * a_dim1];
01266                 }
01267                 i__3 = *m;
01268                 for (i__ = 1; i__ <= i__3; ++i__) {
01269                     c__[i__ + j * c_dim1] += temp1 * b[i__ + k * b_dim1];
01270 /* L150: */
01271                 }
01272 /* L160: */
01273             }
01274 /* L170: */
01275         }
01276     }
01277 
01278     return 0;
01279 
01280 /*     End of SSYMM . */
01281 
01282 } /* ssymm_ */
01283 
01284 /* Subroutine */ int ssyrk_(char *uplo, char *trans, integer *n, integer *k,
01285         real *alpha, real *a, integer *lda, real *beta, real *c__, integer *
01286         ldc)
01287 {
01288     /* System generated locals */
01289     integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2, i__3;
01290 
01291     /* Local variables */
01292     static integer i__, j, l, info;
01293     static real temp;
01294     extern logical lsame_(char *, char *);
01295     static integer nrowa;
01296     static logical upper;
01297     extern /* Subroutine */ int xerbla_(char *, integer *);
01298 
01299 
01300 /*
01301     Purpose
01302     =======
01303 
01304     SSYRK  performs one of the symmetric rank k operations
01305 
01306        C := alpha*A*A' + beta*C,
01307 
01308     or
01309 
01310        C := alpha*A'*A + beta*C,
01311 
01312     where  alpha and beta  are scalars, C is an  n by n  symmetric matrix
01313     and  A  is an  n by k  matrix in the first case and a  k by n  matrix
01314     in the second case.
01315 
01316     Parameters
01317     ==========
01318 
01319     UPLO   - CHARACTER*1.
01320              On  entry,   UPLO  specifies  whether  the  upper  or  lower
01321              triangular  part  of the  array  C  is to be  referenced  as
01322              follows:
01323 
01324                 UPLO = 'U' or 'u'   Only the  upper triangular part of  C
01325                                     is to be referenced.
01326 
01327                 UPLO = 'L' or 'l'   Only the  lower triangular part of  C
01328                                     is to be referenced.
01329 
01330              Unchanged on exit.
01331 
01332     TRANS  - CHARACTER*1.
01333              On entry,  TRANS  specifies the operation to be performed as
01334              follows:
01335 
01336                 TRANS = 'N' or 'n'   C := alpha*A*A' + beta*C.
01337 
01338                 TRANS = 'T' or 't'   C := alpha*A'*A + beta*C.
01339 
01340                 TRANS = 'C' or 'c'   C := alpha*A'*A + beta*C.
01341 
01342              Unchanged on exit.
01343 
01344     N      - INTEGER.
01345              On entry,  N specifies the order of the matrix C.  N must be
01346              at least zero.
01347              Unchanged on exit.
01348 
01349     K      - INTEGER.
01350              On entry with  TRANS = 'N' or 'n',  K  specifies  the number
01351              of  columns   of  the   matrix   A,   and  on   entry   with
01352              TRANS = 'T' or 't' or 'C' or 'c',  K  specifies  the  number
01353              of rows of the matrix  A.  K must be at least zero.
01354              Unchanged on exit.
01355 
01356     ALPHA  - REAL            .
01357              On entry, ALPHA specifies the scalar alpha.
01358              Unchanged on exit.
01359 
01360     A      - REAL             array of DIMENSION ( LDA, ka ), where ka is
01361              k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
01362              Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
01363              part of the array  A  must contain the matrix  A,  otherwise
01364              the leading  k by n  part of the array  A  must contain  the
01365              matrix A.
01366              Unchanged on exit.
01367 
01368     LDA    - INTEGER.
01369              On entry, LDA specifies the first dimension of A as declared
01370              in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
01371              then  LDA must be at least  max( 1, n ), otherwise  LDA must
01372              be at least  max( 1, k ).
01373              Unchanged on exit.
01374 
01375     BETA   - REAL            .
01376              On entry, BETA specifies the scalar beta.
01377              Unchanged on exit.
01378 
01379     C      - REAL             array of DIMENSION ( LDC, n ).
01380              Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
01381              upper triangular part of the array C must contain the upper
01382              triangular part  of the  symmetric matrix  and the strictly
01383              lower triangular part of C is not referenced.  On exit, the
01384              upper triangular part of the array  C is overwritten by the
01385              upper triangular part of the updated matrix.
01386              Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
01387              lower triangular part of the array C must contain the lower
01388              triangular part  of the  symmetric matrix  and the strictly
01389              upper triangular part of C is not referenced.  On exit, the
01390              lower triangular part of the array  C is overwritten by the
01391              lower triangular part of the updated matrix.
01392 
01393     LDC    - INTEGER.
01394              On entry, LDC specifies the first dimension of C as declared
01395              in  the  calling  (sub)  program.   LDC  must  be  at  least
01396              max( 1, n ).
01397              Unchanged on exit.
01398 
01399 
01400     Level 3 Blas routine.
01401 
01402     -- Written on 8-February-1989.
01403        Jack Dongarra, Argonne National Laboratory.
01404        Iain Duff, AERE Harwell.
01405        Jeremy Du Croz, Numerical Algorithms Group Ltd.
01406        Sven Hammarling, Numerical Algorithms Group Ltd.
01407 
01408 
01409        Test the input parameters.
01410 */
01411 
01412     /* Parameter adjustments */
01413     a_dim1 = *lda;
01414     a_offset = 1 + a_dim1;
01415     a -= a_offset;
01416     c_dim1 = *ldc;
01417     c_offset = 1 + c_dim1;
01418     c__ -= c_offset;
01419 
01420     /* Function Body */
01421     if (lsame_(trans, "N")) {
01422         nrowa = *n;
01423     } else {
01424         nrowa = *k;
01425     }
01426     upper = lsame_(uplo, "U");
01427 
01428     info = 0;
01429     if (! upper && ! lsame_(uplo, "L")) {
01430         info = 1;
01431     } else if (! lsame_(trans, "N") && ! lsame_(trans,
01432             "T") && ! lsame_(trans, "C")) {
01433         info = 2;
01434     } else if (*n < 0) {
01435         info = 3;
01436     } else if (*k < 0) {
01437         info = 4;
01438     } else if (*lda < max(1,nrowa)) {
01439         info = 7;
01440     } else if (*ldc < max(1,*n)) {
01441         info = 10;
01442     }
01443     if (info != 0) {
01444         xerbla_("SSYRK ", &info);
01445         return 0;
01446     }
01447 
01448 /*     Quick return if possible. */
01449 
01450     if (*n == 0 || (*alpha == 0.f || *k == 0) && *beta == 1.f) {
01451         return 0;
01452     }
01453 
01454 /*     And when  alpha.eq.zero. */
01455 
01456     if (*alpha == 0.f) {
01457         if (upper) {
01458             if (*beta == 0.f) {
01459                 i__1 = *n;
01460                 for (j = 1; j <= i__1; ++j) {
01461                     i__2 = j;
01462                     for (i__ = 1; i__ <= i__2; ++i__) {
01463                         c__[i__ + j * c_dim1] = 0.f;
01464 /* L10: */
01465                     }
01466 /* L20: */
01467                 }
01468             } else {
01469                 i__1 = *n;
01470                 for (j = 1; j <= i__1; ++j) {
01471                     i__2 = j;
01472                     for (i__ = 1; i__ <= i__2; ++i__) {
01473                         c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
01474 /* L30: */
01475                     }
01476 /* L40: */
01477                 }
01478             }
01479         } else {
01480             if (*beta == 0.f) {
01481                 i__1 = *n;
01482                 for (j = 1; j <= i__1; ++j) {
01483                     i__2 = *n;
01484                     for (i__ = j; i__ <= i__2; ++i__) {
01485                         c__[i__ + j * c_dim1] = 0.f;
01486 /* L50: */
01487                     }
01488 /* L60: */
01489                 }
01490             } else {
01491                 i__1 = *n;
01492                 for (j = 1; j <= i__1; ++j) {
01493                     i__2 = *n;
01494                     for (i__ = j; i__ <= i__2; ++i__) {
01495                         c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
01496 /* L70: */
01497                     }
01498 /* L80: */
01499                 }
01500             }
01501         }
01502         return 0;
01503     }
01504 
01505 /*     Start the operations. */
01506 
01507     if (lsame_(trans, "N")) {
01508 
01509 /*        Form  C := alpha*A*A' + beta*C. */
01510 
01511         if (upper) {
01512             i__1 = *n;
01513             for (j = 1; j <= i__1; ++j) {
01514                 if (*beta == 0.f) {
01515                     i__2 = j;
01516                     for (i__ = 1; i__ <= i__2; ++i__) {
01517                         c__[i__ + j * c_dim1] = 0.f;
01518 /* L90: */
01519                     }
01520                 } else if (*beta != 1.f) {
01521                     i__2 = j;
01522                     for (i__ = 1; i__ <= i__2; ++i__) {
01523                         c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
01524 /* L100: */
01525                     }
01526                 }
01527                 i__2 = *k;
01528                 for (l = 1; l <= i__2; ++l) {
01529                     if (a[j + l * a_dim1] != 0.f) {
01530                         temp = *alpha * a[j + l * a_dim1];
01531                         i__3 = j;
01532                         for (i__ = 1; i__ <= i__3; ++i__) {
01533                             c__[i__ + j * c_dim1] += temp * a[i__ + l *
01534                                     a_dim1];
01535 /* L110: */
01536                         }
01537                     }
01538 /* L120: */
01539                 }
01540 /* L130: */
01541             }
01542         } else {
01543             i__1 = *n;
01544             for (j = 1; j <= i__1; ++j) {
01545                 if (*beta == 0.f) {
01546                     i__2 = *n;
01547                     for (i__ = j; i__ <= i__2; ++i__) {
01548                         c__[i__ + j * c_dim1] = 0.f;
01549 /* L140: */
01550                     }
01551                 } else if (*beta != 1.f) {
01552                     i__2 = *n;
01553                     for (i__ = j; i__ <= i__2; ++i__) {
01554                         c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
01555 /* L150: */
01556                     }
01557                 }
01558                 i__2 = *k;
01559                 for (l = 1; l <= i__2; ++l) {
01560                     if (a[j + l * a_dim1] != 0.f) {
01561                         temp = *alpha * a[j + l * a_dim1];
01562                         i__3 = *n;
01563                         for (i__ = j; i__ <= i__3; ++i__) {
01564                             c__[i__ + j * c_dim1] += temp * a[i__ + l *
01565                                     a_dim1];
01566 /* L160: */
01567                         }
01568                     }
01569 /* L170: */
01570                 }
01571 /* L180: */
01572             }
01573         }
01574     } else {
01575 
01576 /*        Form  C := alpha*A'*A + beta*C. */
01577 
01578         if (upper) {
01579             i__1 = *n;
01580             for (j = 1; j <= i__1; ++j) {
01581                 i__2 = j;
01582                 for (i__ = 1; i__ <= i__2; ++i__) {
01583                     temp = 0.f;
01584                     i__3 = *k;
01585                     for (l = 1; l <= i__3; ++l) {
01586                         temp += a[l + i__ * a_dim1] * a[l + j * a_dim1];
01587 /* L190: */
01588                     }
01589                     if (*beta == 0.f) {
01590                         c__[i__ + j * c_dim1] = *alpha * temp;
01591                     } else {
01592                         c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
01593                                 i__ + j * c_dim1];
01594                     }
01595 /* L200: */
01596                 }
01597 /* L210: */
01598             }
01599         } else {
01600             i__1 = *n;
01601             for (j = 1; j <= i__1; ++j) {
01602                 i__2 = *n;
01603                 for (i__ = j; i__ <= i__2; ++i__) {
01604                     temp = 0.f;
01605                     i__3 = *k;
01606                     for (l = 1; l <= i__3; ++l) {
01607                         temp += a[l + i__ * a_dim1] * a[l + j * a_dim1];
01608 /* L220: */
01609                     }
01610                     if (*beta == 0.f) {
01611                         c__[i__ + j * c_dim1] = *alpha * temp;
01612                     } else {
01613                         c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
01614                                 i__ + j * c_dim1];
01615                     }
01616 /* L230: */
01617                 }
01618 /* L240: */
01619             }
01620         }
01621     }
01622 
01623     return 0;
01624 
01625 /*     End of SSYRK . */
01626 
01627 } /* ssyrk_ */
01628 
01629 /* Subroutine */ int strsm_(char *side, char *uplo, char *transa, char *diag,
01630         integer *m, integer *n, real *alpha, real *a, integer *lda, real *b,
01631         integer *ldb)
01632 {
01633     /* System generated locals */
01634     integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
01635 
01636     /* Local variables */
01637     static integer i__, j, k, info;
01638     static real temp;
01639     static logical lside;
01640     extern logical lsame_(char *, char *);
01641     static integer nrowa;
01642     static logical upper;
01643     extern /* Subroutine */ int xerbla_(char *, integer *);
01644     static logical nounit;
01645 
01646 
01647 /*
01648     Purpose
01649     =======
01650 
01651     STRSM  solves one of the matrix equations
01652 
01653        op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
01654 
01655     where alpha is a scalar, X and B are m by n matrices, A is a unit, or
01656     non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
01657 
01658        op( A ) = A   or   op( A ) = A'.
01659 
01660     The matrix X is overwritten on B.
01661 
01662     Parameters
01663     ==========
01664 
01665     SIDE   - CHARACTER*1.
01666              On entry, SIDE specifies whether op( A ) appears on the left
01667              or right of X as follows:
01668 
01669                 SIDE = 'L' or 'l'   op( A )*X = alpha*B.
01670 
01671                 SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
01672 
01673              Unchanged on exit.
01674 
01675     UPLO   - CHARACTER*1.
01676              On entry, UPLO specifies whether the matrix A is an upper or
01677              lower triangular matrix as follows:
01678 
01679                 UPLO = 'U' or 'u'   A is an upper triangular matrix.
01680 
01681                 UPLO = 'L' or 'l'   A is a lower triangular matrix.
01682 
01683              Unchanged on exit.
01684 
01685     TRANSA - CHARACTER*1.
01686              On entry, TRANSA specifies the form of op( A ) to be used in
01687              the matrix multiplication as follows:
01688 
01689                 TRANSA = 'N' or 'n'   op( A ) = A.
01690 
01691                 TRANSA = 'T' or 't'   op( A ) = A'.
01692 
01693                 TRANSA = 'C' or 'c'   op( A ) = A'.
01694 
01695              Unchanged on exit.
01696 
01697     DIAG   - CHARACTER*1.
01698              On entry, DIAG specifies whether or not A is unit triangular
01699              as follows:
01700 
01701                 DIAG = 'U' or 'u'   A is assumed to be unit triangular.
01702 
01703                 DIAG = 'N' or 'n'   A is not assumed to be unit
01704                                     triangular.
01705 
01706              Unchanged on exit.
01707 
01708     M      - INTEGER.
01709              On entry, M specifies the number of rows of B. M must be at
01710              least zero.
01711              Unchanged on exit.
01712 
01713     N      - INTEGER.
01714              On entry, N specifies the number of columns of B.  N must be
01715              at least zero.
01716              Unchanged on exit.
01717 
01718     ALPHA  - REAL            .
01719              On entry,  ALPHA specifies the scalar  alpha. When  alpha is
01720              zero then  A is not referenced and  B need not be set before
01721              entry.
01722              Unchanged on exit.
01723 
01724     A      - REAL             array of DIMENSION ( LDA, k ), where k is m
01725              when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
01726              Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
01727              upper triangular part of the array  A must contain the upper
01728              triangular matrix  and the strictly lower triangular part of
01729              A is not referenced.
01730              Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
01731              lower triangular part of the array  A must contain the lower
01732              triangular matrix  and the strictly upper triangular part of
01733              A is not referenced.
01734              Note that when  DIAG = 'U' or 'u',  the diagonal elements of
01735              A  are not referenced either,  but are assumed to be  unity.
01736              Unchanged on exit.
01737 
01738     LDA    - INTEGER.
01739              On entry, LDA specifies the first dimension of A as declared
01740              in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
01741              LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
01742              then LDA must be at least max( 1, n ).
01743              Unchanged on exit.
01744 
01745     B      - REAL             array of DIMENSION ( LDB, n ).
01746              Before entry,  the leading  m by n part of the array  B must
01747              contain  the  right-hand  side  matrix  B,  and  on exit  is
01748              overwritten by the solution matrix  X.
01749 
01750     LDB    - INTEGER.
01751              On entry, LDB specifies the first dimension of B as declared
01752              in  the  calling  (sub)  program.   LDB  must  be  at  least
01753              max( 1, m ).
01754              Unchanged on exit.
01755 
01756 
01757     Level 3 Blas routine.
01758 
01759 
01760     -- Written on 8-February-1989.
01761        Jack Dongarra, Argonne National Laboratory.
01762        Iain Duff, AERE Harwell.
01763        Jeremy Du Croz, Numerical Algorithms Group Ltd.
01764        Sven Hammarling, Numerical Algorithms Group Ltd.
01765 
01766 
01767        Test the input parameters.
01768 */
01769 
01770     /* Parameter adjustments */
01771     a_dim1 = *lda;
01772     a_offset = 1 + a_dim1;
01773     a -= a_offset;
01774     b_dim1 = *ldb;
01775     b_offset = 1 + b_dim1;
01776     b -= b_offset;
01777 
01778     /* Function Body */
01779     lside = lsame_(side, "L");
01780     if (lside) {
01781         nrowa = *m;
01782     } else {
01783         nrowa = *n;
01784     }
01785     nounit = lsame_(diag, "N");
01786     upper = lsame_(uplo, "U");
01787 
01788     info = 0;
01789     if (! lside && ! lsame_(side, "R")) {
01790         info = 1;
01791     } else if (! upper && ! lsame_(uplo, "L")) {
01792         info = 2;
01793     } else if (! lsame_(transa, "N") && ! lsame_(transa,
01794              "T") && ! lsame_(transa, "C")) {
01795         info = 3;
01796     } else if (! lsame_(diag, "U") && ! lsame_(diag,
01797             "N")) {
01798         info = 4;
01799     } else if (*m < 0) {
01800         info = 5;
01801     } else if (*n < 0) {
01802         info = 6;
01803     } else if (*lda < max(1,nrowa)) {
01804         info = 9;
01805     } else if (*ldb < max(1,*m)) {
01806         info = 11;
01807     }
01808     if (info != 0) {
01809         xerbla_("STRSM ", &info);
01810         return 0;
01811     }
01812 
01813 /*     Quick return if possible. */
01814 
01815     if (*n == 0) {
01816         return 0;
01817     }
01818 
01819 /*     And when  alpha.eq.zero. */
01820 
01821     if (*alpha == 0.f) {
01822         i__1 = *n;
01823         for (j = 1; j <= i__1; ++j) {
01824             i__2 = *m;
01825             for (i__ = 1; i__ <= i__2; ++i__) {
01826                 b[i__ + j * b_dim1] = 0.f;
01827 /* L10: */
01828             }
01829 /* L20: */
01830         }
01831         return 0;
01832     }
01833 
01834 /*     Start the operations. */
01835 
01836     if (lside) {
01837         if (lsame_(transa, "N")) {
01838 
01839 /*           Form  B := alpha*inv( A )*B. */
01840 
01841             if (upper) {
01842                 i__1 = *n;
01843                 for (j = 1; j <= i__1; ++j) {
01844                     if (*alpha != 1.f) {
01845                         i__2 = *m;
01846                         for (i__ = 1; i__ <= i__2; ++i__) {
01847                             b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
01848                                     ;
01849 /* L30: */
01850                         }
01851                     }
01852                     for (k = *m; k >= 1; --k) {
01853                         if (b[k + j * b_dim1] != 0.f) {
01854                             if (nounit) {
01855                                 b[k + j * b_dim1] /= a[k + k * a_dim1];
01856                             }
01857                             i__2 = k - 1;
01858                             for (i__ = 1; i__ <= i__2; ++i__) {
01859                                 b[i__ + j * b_dim1] -= b[k + j * b_dim1] * a[
01860                                         i__ + k * a_dim1];
01861 /* L40: */
01862                             }
01863                         }
01864 /* L50: */
01865                     }
01866 /* L60: */
01867                 }
01868             } else {
01869                 i__1 = *n;
01870                 for (j = 1; j <= i__1; ++j) {
01871                     if (*alpha != 1.f) {
01872                         i__2 = *m;
01873                         for (i__ = 1; i__ <= i__2; ++i__) {
01874                             b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
01875                                     ;
01876 /* L70: */
01877                         }
01878                     }
01879                     i__2 = *m;
01880                     for (k = 1; k <= i__2; ++k) {
01881                         if (b[k + j * b_dim1] != 0.f) {
01882                             if (nounit) {
01883                                 b[k + j * b_dim1] /= a[k + k * a_dim1];
01884                             }
01885                             i__3 = *m;
01886                             for (i__ = k + 1; i__ <= i__3; ++i__) {
01887                                 b[i__ + j * b_dim1] -= b[k + j * b_dim1] * a[
01888                                         i__ + k * a_dim1];
01889 /* L80: */
01890                             }
01891                         }
01892 /* L90: */
01893                     }
01894 /* L100: */
01895                 }
01896             }
01897         } else {
01898 
01899 /*           Form  B := alpha*inv( A' )*B. */
01900 
01901             if (upper) {
01902                 i__1 = *n;
01903                 for (j = 1; j <= i__1; ++j) {
01904                     i__2 = *m;
01905                     for (i__ = 1; i__ <= i__2; ++i__) {
01906                         temp = *alpha * b[i__ + j * b_dim1];
01907                         i__3 = i__ - 1;
01908                         for (k = 1; k <= i__3; ++k) {
01909                             temp -= a[k + i__ * a_dim1] * b[k + j * b_dim1];
01910 /* L110: */
01911                         }
01912                         if (nounit) {
01913                             temp /= a[i__ + i__ * a_dim1];
01914                         }
01915                         b[i__ + j * b_dim1] = temp;
01916 /* L120: */
01917                     }
01918 /* L130: */
01919                 }
01920             } else {
01921                 i__1 = *n;
01922                 for (j = 1; j <= i__1; ++j) {
01923                     for (i__ = *m; i__ >= 1; --i__) {
01924                         temp = *alpha * b[i__ + j * b_dim1];
01925                         i__2 = *m;
01926                         for (k = i__ + 1; k <= i__2; ++k) {
01927                             temp -= a[k + i__ * a_dim1] * b[k + j * b_dim1];
01928 /* L140: */
01929                         }
01930                         if (nounit) {
01931                             temp /= a[i__ + i__ * a_dim1];
01932                         }
01933                         b[i__ + j * b_dim1] = temp;
01934 /* L150: */
01935                     }
01936 /* L160: */
01937                 }
01938             }
01939         }
01940     } else {
01941         if (lsame_(transa, "N")) {
01942 
01943 /*           Form  B := alpha*B*inv( A ). */
01944 
01945             if (upper) {
01946                 i__1 = *n;
01947                 for (j = 1; j <= i__1; ++j) {
01948                     if (*alpha != 1.f) {
01949                         i__2 = *m;
01950                         for (i__ = 1; i__ <= i__2; ++i__) {
01951                             b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
01952                                     ;
01953 /* L170: */
01954                         }
01955                     }
01956                     i__2 = j - 1;
01957                     for (k = 1; k <= i__2; ++k) {
01958                         if (a[k + j * a_dim1] != 0.f) {
01959                             i__3 = *m;
01960                             for (i__ = 1; i__ <= i__3; ++i__) {
01961                                 b[i__ + j * b_dim1] -= a[k + j * a_dim1] * b[
01962                                         i__ + k * b_dim1];
01963 /* L180: */
01964                             }
01965                         }
01966 /* L190: */
01967                     }
01968                     if (nounit) {
01969                         temp = 1.f / a[j + j * a_dim1];
01970                         i__2 = *m;
01971                         for (i__ = 1; i__ <= i__2; ++i__) {
01972                             b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
01973 /* L200: */
01974                         }
01975                     }
01976 /* L210: */
01977                 }
01978             } else {
01979                 for (j = *n; j >= 1; --j) {
01980                     if (*alpha != 1.f) {
01981                         i__1 = *m;
01982                         for (i__ = 1; i__ <= i__1; ++i__) {
01983                             b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
01984                                     ;
01985 /* L220: */
01986                         }
01987                     }
01988                     i__1 = *n;
01989                     for (k = j + 1; k <= i__1; ++k) {
01990                         if (a[k + j * a_dim1] != 0.f) {
01991                             i__2 = *m;
01992                             for (i__ = 1; i__ <= i__2; ++i__) {
01993                                 b[i__ + j * b_dim1] -= a[k + j * a_dim1] * b[
01994                                         i__ + k * b_dim1];
01995 /* L230: */
01996                             }
01997                         }
01998 /* L240: */
01999                     }
02000                     if (nounit) {
02001                         temp = 1.f / a[j + j * a_dim1];
02002                         i__1 = *m;
02003                         for (i__ = 1; i__ <= i__1; ++i__) {
02004                             b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
02005 /* L250: */
02006                         }
02007                     }
02008 /* L260: */
02009                 }
02010             }
02011         } else {
02012 
02013 /*           Form  B := alpha*B*inv( A' ). */
02014 
02015             if (upper) {
02016                 for (k = *n; k >= 1; --k) {
02017                     if (nounit) {
02018                         temp = 1.f / a[k + k * a_dim1];
02019                         i__1 = *m;
02020                         for (i__ = 1; i__ <= i__1; ++i__) {
02021                             b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
02022 /* L270: */
02023                         }
02024                     }
02025                     i__1 = k - 1;
02026                     for (j = 1; j <= i__1; ++j) {
02027                         if (a[j + k * a_dim1] != 0.f) {
02028                             temp = a[j + k * a_dim1];
02029                             i__2 = *m;
02030                             for (i__ = 1; i__ <= i__2; ++i__) {
02031                                 b[i__ + j * b_dim1] -= temp * b[i__ + k *
02032                                         b_dim1];
02033 /* L280: */
02034                             }
02035                         }
02036 /* L290: */
02037                     }
02038                     if (*alpha != 1.f) {
02039                         i__1 = *m;
02040                         for (i__ = 1; i__ <= i__1; ++i__) {
02041                             b[i__ + k * b_dim1] = *alpha * b[i__ + k * b_dim1]
02042                                     ;
02043 /* L300: */
02044                         }
02045                     }
02046 /* L310: */
02047                 }
02048             } else {
02049                 i__1 = *n;
02050                 for (k = 1; k <= i__1; ++k) {
02051                     if (nounit) {
02052                         temp = 1.f / a[k + k * a_dim1];
02053                         i__2 = *m;
02054                         for (i__ = 1; i__ <= i__2; ++i__) {
02055                             b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
02056 /* L320: */
02057                         }
02058                     }
02059                     i__2 = *n;
02060                     for (j = k + 1; j <= i__2; ++j) {
02061                         if (a[j + k * a_dim1] != 0.f) {
02062                             temp = a[j + k * a_dim1];
02063                             i__3 = *m;
02064                             for (i__ = 1; i__ <= i__3; ++i__) {
02065                                 b[i__ + j * b_dim1] -= temp * b[i__ + k *
02066                                         b_dim1];
02067 /* L330: */
02068                             }
02069                         }
02070 /* L340: */
02071                     }
02072                     if (*alpha != 1.f) {
02073                         i__2 = *m;
02074                         for (i__ = 1; i__ <= i__2; ++i__) {
02075                             b[i__ + k * b_dim1] = *alpha * b[i__ + k * b_dim1]
02076                                     ;
02077 /* L350: */
02078                         }
02079                     }
02080 /* L360: */
02081                 }
02082             }
02083         }
02084     }
02085 
02086     return 0;
02087 
02088 /*     End of STRSM . */
02089 
02090 } /* strsm_ */
02091 
02092 /* Subroutine */ int xerbla_(char *srname, integer *info)
02093 {
02094     /* Format strings */
02095     static char fmt_9999[] = "(\002 ** On entry to \002,a6,\002 parameter nu"
02096             "mber \002,i2,\002 had \002,\002an illegal value\002)";
02097 
02098     /* Builtin functions */
02099     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
02100     /* Subroutine */ int s_stop(char *, ftnlen);
02101 
02102     /* Fortran I/O blocks */
02103     static cilist io___60 = { 0, 6, 0, fmt_9999, 0 };
02104 
02105 
02106 /*
02107     -- LAPACK auxiliary routine (preliminary version) --
02108        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
02109        Courant Institute, Argonne National Lab, and Rice University
02110        February 29, 1992
02111 
02112 
02113     Purpose
02114     =======
02115 
02116     XERBLA  is an error handler for the LAPACK routines.
02117     It is called by an LAPACK routine if an input parameter has an
02118     invalid value.  A message is printed and execution stops.
02119 
02120     Installers may consider modifying the STOP statement in order to
02121     call system-specific exception-handling facilities.
02122 
02123     Arguments
02124     =========
02125 
02126     SRNAME  (input) CHARACTER*6
02127             The name of the routine which called XERBLA.
02128 
02129     INFO    (input) INTEGER
02130             The position of the invalid parameter in the parameter list
02131             of the calling routine.
02132 */
02133 
02134 
02135     s_wsfe(&io___60);
02136     do_fio(&c__1, srname, (ftnlen)6);
02137     do_fio(&c__1, (char *)&(*info), (ftnlen)sizeof(integer));
02138     e_wsfe();
02139 
02140     s_stop("", (ftnlen)0);
02141 
02142 
02143 /*     End of XERBLA */
02144 
02145     return 0;
02146 } /* xerbla_ */
02147