SphinxBase 0.6
|
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