src/libsphinxbase/util/blas_lite.c

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