|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE STFSM( TRANSR, SIDE, UPLO, TRANS, DIAG, M, N, ALPHA, A, 00002 + B, LDB ) 00003 * 00004 * -- LAPACK routine (version 3.3.0) -- 00005 * 00006 * -- Contributed by Fred Gustavson of the IBM Watson Research Center -- 00007 * November 2010 00008 * 00009 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00010 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00011 * 00012 * .. 00013 * .. Scalar Arguments .. 00014 CHARACTER TRANSR, DIAG, SIDE, TRANS, UPLO 00015 INTEGER LDB, M, N 00016 REAL ALPHA 00017 * .. 00018 * .. Array Arguments .. 00019 REAL A( 0: * ), B( 0: LDB-1, 0: * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * Level 3 BLAS like routine for A in RFP Format. 00026 * 00027 * STFSM solves the matrix equation 00028 * 00029 * op( A )*X = alpha*B or X*op( A ) = alpha*B 00030 * 00031 * where alpha is a scalar, X and B are m by n matrices, A is a unit, or 00032 * non-unit, upper or lower triangular matrix and op( A ) is one of 00033 * 00034 * op( A ) = A or op( A ) = A'. 00035 * 00036 * A is in Rectangular Full Packed (RFP) Format. 00037 * 00038 * The matrix X is overwritten on B. 00039 * 00040 * Arguments 00041 * ========== 00042 * 00043 * TRANSR (input) CHARACTER*1 00044 * = 'N': The Normal Form of RFP A is stored; 00045 * = 'T': The Transpose Form of RFP A is stored. 00046 * 00047 * SIDE (input) CHARACTER*1 00048 * On entry, SIDE specifies whether op( A ) appears on the left 00049 * or right of X as follows: 00050 * 00051 * SIDE = 'L' or 'l' op( A )*X = alpha*B. 00052 * 00053 * SIDE = 'R' or 'r' X*op( A ) = alpha*B. 00054 * 00055 * Unchanged on exit. 00056 * 00057 * UPLO (input) CHARACTER*1 00058 * On entry, UPLO specifies whether the RFP matrix A came from 00059 * an upper or lower triangular matrix as follows: 00060 * UPLO = 'U' or 'u' RFP A came from an upper triangular matrix 00061 * UPLO = 'L' or 'l' RFP A came from a lower triangular matrix 00062 * 00063 * Unchanged on exit. 00064 * 00065 * TRANS (input) CHARACTER*1 00066 * On entry, TRANS specifies the form of op( A ) to be used 00067 * in the matrix multiplication as follows: 00068 * 00069 * TRANS = 'N' or 'n' op( A ) = A. 00070 * 00071 * TRANS = 'T' or 't' op( A ) = A'. 00072 * 00073 * Unchanged on exit. 00074 * 00075 * DIAG (input) CHARACTER*1 00076 * On entry, DIAG specifies whether or not RFP A is unit 00077 * triangular as follows: 00078 * 00079 * DIAG = 'U' or 'u' A is assumed to be unit triangular. 00080 * 00081 * DIAG = 'N' or 'n' A is not assumed to be unit 00082 * triangular. 00083 * 00084 * Unchanged on exit. 00085 * 00086 * M (input) INTEGER 00087 * On entry, M specifies the number of rows of B. M must be at 00088 * least zero. 00089 * Unchanged on exit. 00090 * 00091 * N (input) INTEGER 00092 * On entry, N specifies the number of columns of B. N must be 00093 * at least zero. 00094 * Unchanged on exit. 00095 * 00096 * ALPHA (input) REAL 00097 * On entry, ALPHA specifies the scalar alpha. When alpha is 00098 * zero then A is not referenced and B need not be set before 00099 * entry. 00100 * Unchanged on exit. 00101 * 00102 * A (input) REAL array, dimension (NT) 00103 * NT = N*(N+1)/2. On entry, the matrix A in RFP Format. 00104 * RFP Format is described by TRANSR, UPLO and N as follows: 00105 * If TRANSR='N' then RFP A is (0:N,0:K-1) when N is even; 00106 * K=N/2. RFP A is (0:N-1,0:K) when N is odd; K=N/2. If 00107 * TRANSR = 'T' then RFP is the transpose of RFP A as 00108 * defined when TRANSR = 'N'. The contents of RFP A are defined 00109 * by UPLO as follows: If UPLO = 'U' the RFP A contains the NT 00110 * elements of upper packed A either in normal or 00111 * transpose Format. If UPLO = 'L' the RFP A contains 00112 * the NT elements of lower packed A either in normal or 00113 * transpose Format. The LDA of RFP A is (N+1)/2 when 00114 * TRANSR = 'T'. When TRANSR is 'N' the LDA is N+1 when N is 00115 * even and is N when is odd. 00116 * See the Note below for more details. Unchanged on exit. 00117 * 00118 * B (input/output) REAL array, DIMENSION (LDB,N) 00119 * Before entry, the leading m by n part of the array B must 00120 * contain the right-hand side matrix B, and on exit is 00121 * overwritten by the solution matrix X. 00122 * 00123 * LDB (input) INTEGER 00124 * On entry, LDB specifies the first dimension of B as declared 00125 * in the calling (sub) program. LDB must be at least 00126 * max( 1, m ). 00127 * Unchanged on exit. 00128 * 00129 * Further Details 00130 * =============== 00131 * 00132 * We first consider Rectangular Full Packed (RFP) Format when N is 00133 * even. We give an example where N = 6. 00134 * 00135 * AP is Upper AP is Lower 00136 * 00137 * 00 01 02 03 04 05 00 00138 * 11 12 13 14 15 10 11 00139 * 22 23 24 25 20 21 22 00140 * 33 34 35 30 31 32 33 00141 * 44 45 40 41 42 43 44 00142 * 55 50 51 52 53 54 55 00143 * 00144 * 00145 * Let TRANSR = 'N'. RFP holds AP as follows: 00146 * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last 00147 * three columns of AP upper. The lower triangle A(4:6,0:2) consists of 00148 * the transpose of the first three columns of AP upper. 00149 * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first 00150 * three columns of AP lower. The upper triangle A(0:2,0:2) consists of 00151 * the transpose of the last three columns of AP lower. 00152 * This covers the case N even and TRANSR = 'N'. 00153 * 00154 * RFP A RFP A 00155 * 00156 * 03 04 05 33 43 53 00157 * 13 14 15 00 44 54 00158 * 23 24 25 10 11 55 00159 * 33 34 35 20 21 22 00160 * 00 44 45 30 31 32 00161 * 01 11 55 40 41 42 00162 * 02 12 22 50 51 52 00163 * 00164 * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the 00165 * transpose of RFP A above. One therefore gets: 00166 * 00167 * 00168 * RFP A RFP A 00169 * 00170 * 03 13 23 33 00 01 02 33 00 10 20 30 40 50 00171 * 04 14 24 34 44 11 12 43 44 11 21 31 41 51 00172 * 05 15 25 35 45 55 22 53 54 55 22 32 42 52 00173 * 00174 * 00175 * We then consider Rectangular Full Packed (RFP) Format when N is 00176 * odd. We give an example where N = 5. 00177 * 00178 * AP is Upper AP is Lower 00179 * 00180 * 00 01 02 03 04 00 00181 * 11 12 13 14 10 11 00182 * 22 23 24 20 21 22 00183 * 33 34 30 31 32 33 00184 * 44 40 41 42 43 44 00185 * 00186 * 00187 * Let TRANSR = 'N'. RFP holds AP as follows: 00188 * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last 00189 * three columns of AP upper. The lower triangle A(3:4,0:1) consists of 00190 * the transpose of the first two columns of AP upper. 00191 * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first 00192 * three columns of AP lower. The upper triangle A(0:1,1:2) consists of 00193 * the transpose of the last two columns of AP lower. 00194 * This covers the case N odd and TRANSR = 'N'. 00195 * 00196 * RFP A RFP A 00197 * 00198 * 02 03 04 00 33 43 00199 * 12 13 14 10 11 44 00200 * 22 23 24 20 21 22 00201 * 00 33 34 30 31 32 00202 * 01 11 44 40 41 42 00203 * 00204 * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the 00205 * transpose of RFP A above. One therefore gets: 00206 * 00207 * RFP A RFP A 00208 * 00209 * 02 12 22 00 01 00 10 20 30 40 50 00210 * 03 13 23 33 11 33 11 21 31 41 51 00211 * 04 14 24 34 44 43 44 22 32 42 52 00212 * 00213 * Reference 00214 * ========= 00215 * 00216 * ===================================================================== 00217 * 00218 * .. 00219 * .. Parameters .. 00220 REAL ONE, ZERO 00221 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00222 * .. 00223 * .. Local Scalars .. 00224 LOGICAL LOWER, LSIDE, MISODD, NISODD, NORMALTRANSR, 00225 + NOTRANS 00226 INTEGER M1, M2, N1, N2, K, INFO, I, J 00227 * .. 00228 * .. External Functions .. 00229 LOGICAL LSAME 00230 EXTERNAL LSAME 00231 * .. 00232 * .. External Subroutines .. 00233 EXTERNAL SGEMM, STRSM, XERBLA 00234 * .. 00235 * .. Intrinsic Functions .. 00236 INTRINSIC MAX, MOD 00237 * .. 00238 * .. Executable Statements .. 00239 * 00240 * Test the input parameters. 00241 * 00242 INFO = 0 00243 NORMALTRANSR = LSAME( TRANSR, 'N' ) 00244 LSIDE = LSAME( SIDE, 'L' ) 00245 LOWER = LSAME( UPLO, 'L' ) 00246 NOTRANS = LSAME( TRANS, 'N' ) 00247 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN 00248 INFO = -1 00249 ELSE IF( .NOT.LSIDE .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN 00250 INFO = -2 00251 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 00252 INFO = -3 00253 ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN 00254 INFO = -4 00255 ELSE IF( .NOT.LSAME( DIAG, 'N' ) .AND. .NOT.LSAME( DIAG, 'U' ) ) 00256 + THEN 00257 INFO = -5 00258 ELSE IF( M.LT.0 ) THEN 00259 INFO = -6 00260 ELSE IF( N.LT.0 ) THEN 00261 INFO = -7 00262 ELSE IF( LDB.LT.MAX( 1, M ) ) THEN 00263 INFO = -11 00264 END IF 00265 IF( INFO.NE.0 ) THEN 00266 CALL XERBLA( 'STFSM ', -INFO ) 00267 RETURN 00268 END IF 00269 * 00270 * Quick return when ( (N.EQ.0).OR.(M.EQ.0) ) 00271 * 00272 IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) ) 00273 + RETURN 00274 * 00275 * Quick return when ALPHA.EQ.(0D+0) 00276 * 00277 IF( ALPHA.EQ.ZERO ) THEN 00278 DO 20 J = 0, N - 1 00279 DO 10 I = 0, M - 1 00280 B( I, J ) = ZERO 00281 10 CONTINUE 00282 20 CONTINUE 00283 RETURN 00284 END IF 00285 * 00286 IF( LSIDE ) THEN 00287 * 00288 * SIDE = 'L' 00289 * 00290 * A is M-by-M. 00291 * If M is odd, set NISODD = .TRUE., and M1 and M2. 00292 * If M is even, NISODD = .FALSE., and M. 00293 * 00294 IF( MOD( M, 2 ).EQ.0 ) THEN 00295 MISODD = .FALSE. 00296 K = M / 2 00297 ELSE 00298 MISODD = .TRUE. 00299 IF( LOWER ) THEN 00300 M2 = M / 2 00301 M1 = M - M2 00302 ELSE 00303 M1 = M / 2 00304 M2 = M - M1 00305 END IF 00306 END IF 00307 * 00308 IF( MISODD ) THEN 00309 * 00310 * SIDE = 'L' and N is odd 00311 * 00312 IF( NORMALTRANSR ) THEN 00313 * 00314 * SIDE = 'L', N is odd, and TRANSR = 'N' 00315 * 00316 IF( LOWER ) THEN 00317 * 00318 * SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'L' 00319 * 00320 IF( NOTRANS ) THEN 00321 * 00322 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and 00323 * TRANS = 'N' 00324 * 00325 IF( M.EQ.1 ) THEN 00326 CALL STRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA, 00327 + A, M, B, LDB ) 00328 ELSE 00329 CALL STRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA, 00330 + A( 0 ), M, B, LDB ) 00331 CALL SGEMM( 'N', 'N', M2, N, M1, -ONE, A( M1 ), 00332 + M, B, LDB, ALPHA, B( M1, 0 ), LDB ) 00333 CALL STRSM( 'L', 'U', 'T', DIAG, M2, N, ONE, 00334 + A( M ), M, B( M1, 0 ), LDB ) 00335 END IF 00336 * 00337 ELSE 00338 * 00339 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and 00340 * TRANS = 'T' 00341 * 00342 IF( M.EQ.1 ) THEN 00343 CALL STRSM( 'L', 'L', 'T', DIAG, M1, N, ALPHA, 00344 + A( 0 ), M, B, LDB ) 00345 ELSE 00346 CALL STRSM( 'L', 'U', 'N', DIAG, M2, N, ALPHA, 00347 + A( M ), M, B( M1, 0 ), LDB ) 00348 CALL SGEMM( 'T', 'N', M1, N, M2, -ONE, A( M1 ), 00349 + M, B( M1, 0 ), LDB, ALPHA, B, LDB ) 00350 CALL STRSM( 'L', 'L', 'T', DIAG, M1, N, ONE, 00351 + A( 0 ), M, B, LDB ) 00352 END IF 00353 * 00354 END IF 00355 * 00356 ELSE 00357 * 00358 * SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'U' 00359 * 00360 IF( .NOT.NOTRANS ) THEN 00361 * 00362 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and 00363 * TRANS = 'N' 00364 * 00365 CALL STRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA, 00366 + A( M2 ), M, B, LDB ) 00367 CALL SGEMM( 'T', 'N', M2, N, M1, -ONE, A( 0 ), M, 00368 + B, LDB, ALPHA, B( M1, 0 ), LDB ) 00369 CALL STRSM( 'L', 'U', 'T', DIAG, M2, N, ONE, 00370 + A( M1 ), M, B( M1, 0 ), LDB ) 00371 * 00372 ELSE 00373 * 00374 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and 00375 * TRANS = 'T' 00376 * 00377 CALL STRSM( 'L', 'U', 'N', DIAG, M2, N, ALPHA, 00378 + A( M1 ), M, B( M1, 0 ), LDB ) 00379 CALL SGEMM( 'N', 'N', M1, N, M2, -ONE, A( 0 ), M, 00380 + B( M1, 0 ), LDB, ALPHA, B, LDB ) 00381 CALL STRSM( 'L', 'L', 'T', DIAG, M1, N, ONE, 00382 + A( M2 ), M, B, LDB ) 00383 * 00384 END IF 00385 * 00386 END IF 00387 * 00388 ELSE 00389 * 00390 * SIDE = 'L', N is odd, and TRANSR = 'T' 00391 * 00392 IF( LOWER ) THEN 00393 * 00394 * SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'L' 00395 * 00396 IF( NOTRANS ) THEN 00397 * 00398 * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and 00399 * TRANS = 'N' 00400 * 00401 IF( M.EQ.1 ) THEN 00402 CALL STRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA, 00403 + A( 0 ), M1, B, LDB ) 00404 ELSE 00405 CALL STRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA, 00406 + A( 0 ), M1, B, LDB ) 00407 CALL SGEMM( 'T', 'N', M2, N, M1, -ONE, 00408 + A( M1*M1 ), M1, B, LDB, ALPHA, 00409 + B( M1, 0 ), LDB ) 00410 CALL STRSM( 'L', 'L', 'N', DIAG, M2, N, ONE, 00411 + A( 1 ), M1, B( M1, 0 ), LDB ) 00412 END IF 00413 * 00414 ELSE 00415 * 00416 * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and 00417 * TRANS = 'T' 00418 * 00419 IF( M.EQ.1 ) THEN 00420 CALL STRSM( 'L', 'U', 'N', DIAG, M1, N, ALPHA, 00421 + A( 0 ), M1, B, LDB ) 00422 ELSE 00423 CALL STRSM( 'L', 'L', 'T', DIAG, M2, N, ALPHA, 00424 + A( 1 ), M1, B( M1, 0 ), LDB ) 00425 CALL SGEMM( 'N', 'N', M1, N, M2, -ONE, 00426 + A( M1*M1 ), M1, B( M1, 0 ), LDB, 00427 + ALPHA, B, LDB ) 00428 CALL STRSM( 'L', 'U', 'N', DIAG, M1, N, ONE, 00429 + A( 0 ), M1, B, LDB ) 00430 END IF 00431 * 00432 END IF 00433 * 00434 ELSE 00435 * 00436 * SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'U' 00437 * 00438 IF( .NOT.NOTRANS ) THEN 00439 * 00440 * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and 00441 * TRANS = 'N' 00442 * 00443 CALL STRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA, 00444 + A( M2*M2 ), M2, B, LDB ) 00445 CALL SGEMM( 'N', 'N', M2, N, M1, -ONE, A( 0 ), M2, 00446 + B, LDB, ALPHA, B( M1, 0 ), LDB ) 00447 CALL STRSM( 'L', 'L', 'N', DIAG, M2, N, ONE, 00448 + A( M1*M2 ), M2, B( M1, 0 ), LDB ) 00449 * 00450 ELSE 00451 * 00452 * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and 00453 * TRANS = 'T' 00454 * 00455 CALL STRSM( 'L', 'L', 'T', DIAG, M2, N, ALPHA, 00456 + A( M1*M2 ), M2, B( M1, 0 ), LDB ) 00457 CALL SGEMM( 'T', 'N', M1, N, M2, -ONE, A( 0 ), M2, 00458 + B( M1, 0 ), LDB, ALPHA, B, LDB ) 00459 CALL STRSM( 'L', 'U', 'N', DIAG, M1, N, ONE, 00460 + A( M2*M2 ), M2, B, LDB ) 00461 * 00462 END IF 00463 * 00464 END IF 00465 * 00466 END IF 00467 * 00468 ELSE 00469 * 00470 * SIDE = 'L' and N is even 00471 * 00472 IF( NORMALTRANSR ) THEN 00473 * 00474 * SIDE = 'L', N is even, and TRANSR = 'N' 00475 * 00476 IF( LOWER ) THEN 00477 * 00478 * SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'L' 00479 * 00480 IF( NOTRANS ) THEN 00481 * 00482 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L', 00483 * and TRANS = 'N' 00484 * 00485 CALL STRSM( 'L', 'L', 'N', DIAG, K, N, ALPHA, 00486 + A( 1 ), M+1, B, LDB ) 00487 CALL SGEMM( 'N', 'N', K, N, K, -ONE, A( K+1 ), 00488 + M+1, B, LDB, ALPHA, B( K, 0 ), LDB ) 00489 CALL STRSM( 'L', 'U', 'T', DIAG, K, N, ONE, 00490 + A( 0 ), M+1, B( K, 0 ), LDB ) 00491 * 00492 ELSE 00493 * 00494 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L', 00495 * and TRANS = 'T' 00496 * 00497 CALL STRSM( 'L', 'U', 'N', DIAG, K, N, ALPHA, 00498 + A( 0 ), M+1, B( K, 0 ), LDB ) 00499 CALL SGEMM( 'T', 'N', K, N, K, -ONE, A( K+1 ), 00500 + M+1, B( K, 0 ), LDB, ALPHA, B, LDB ) 00501 CALL STRSM( 'L', 'L', 'T', DIAG, K, N, ONE, 00502 + A( 1 ), M+1, B, LDB ) 00503 * 00504 END IF 00505 * 00506 ELSE 00507 * 00508 * SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'U' 00509 * 00510 IF( .NOT.NOTRANS ) THEN 00511 * 00512 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U', 00513 * and TRANS = 'N' 00514 * 00515 CALL STRSM( 'L', 'L', 'N', DIAG, K, N, ALPHA, 00516 + A( K+1 ), M+1, B, LDB ) 00517 CALL SGEMM( 'T', 'N', K, N, K, -ONE, A( 0 ), M+1, 00518 + B, LDB, ALPHA, B( K, 0 ), LDB ) 00519 CALL STRSM( 'L', 'U', 'T', DIAG, K, N, ONE, 00520 + A( K ), M+1, B( K, 0 ), LDB ) 00521 * 00522 ELSE 00523 * 00524 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U', 00525 * and TRANS = 'T' 00526 CALL STRSM( 'L', 'U', 'N', DIAG, K, N, ALPHA, 00527 + A( K ), M+1, B( K, 0 ), LDB ) 00528 CALL SGEMM( 'N', 'N', K, N, K, -ONE, A( 0 ), M+1, 00529 + B( K, 0 ), LDB, ALPHA, B, LDB ) 00530 CALL STRSM( 'L', 'L', 'T', DIAG, K, N, ONE, 00531 + A( K+1 ), M+1, B, LDB ) 00532 * 00533 END IF 00534 * 00535 END IF 00536 * 00537 ELSE 00538 * 00539 * SIDE = 'L', N is even, and TRANSR = 'T' 00540 * 00541 IF( LOWER ) THEN 00542 * 00543 * SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'L' 00544 * 00545 IF( NOTRANS ) THEN 00546 * 00547 * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L', 00548 * and TRANS = 'N' 00549 * 00550 CALL STRSM( 'L', 'U', 'T', DIAG, K, N, ALPHA, 00551 + A( K ), K, B, LDB ) 00552 CALL SGEMM( 'T', 'N', K, N, K, -ONE, 00553 + A( K*( K+1 ) ), K, B, LDB, ALPHA, 00554 + B( K, 0 ), LDB ) 00555 CALL STRSM( 'L', 'L', 'N', DIAG, K, N, ONE, 00556 + A( 0 ), K, B( K, 0 ), LDB ) 00557 * 00558 ELSE 00559 * 00560 * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L', 00561 * and TRANS = 'T' 00562 * 00563 CALL STRSM( 'L', 'L', 'T', DIAG, K, N, ALPHA, 00564 + A( 0 ), K, B( K, 0 ), LDB ) 00565 CALL SGEMM( 'N', 'N', K, N, K, -ONE, 00566 + A( K*( K+1 ) ), K, B( K, 0 ), LDB, 00567 + ALPHA, B, LDB ) 00568 CALL STRSM( 'L', 'U', 'N', DIAG, K, N, ONE, 00569 + A( K ), K, B, LDB ) 00570 * 00571 END IF 00572 * 00573 ELSE 00574 * 00575 * SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'U' 00576 * 00577 IF( .NOT.NOTRANS ) THEN 00578 * 00579 * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U', 00580 * and TRANS = 'N' 00581 * 00582 CALL STRSM( 'L', 'U', 'T', DIAG, K, N, ALPHA, 00583 + A( K*( K+1 ) ), K, B, LDB ) 00584 CALL SGEMM( 'N', 'N', K, N, K, -ONE, A( 0 ), K, B, 00585 + LDB, ALPHA, B( K, 0 ), LDB ) 00586 CALL STRSM( 'L', 'L', 'N', DIAG, K, N, ONE, 00587 + A( K*K ), K, B( K, 0 ), LDB ) 00588 * 00589 ELSE 00590 * 00591 * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U', 00592 * and TRANS = 'T' 00593 * 00594 CALL STRSM( 'L', 'L', 'T', DIAG, K, N, ALPHA, 00595 + A( K*K ), K, B( K, 0 ), LDB ) 00596 CALL SGEMM( 'T', 'N', K, N, K, -ONE, A( 0 ), K, 00597 + B( K, 0 ), LDB, ALPHA, B, LDB ) 00598 CALL STRSM( 'L', 'U', 'N', DIAG, K, N, ONE, 00599 + A( K*( K+1 ) ), K, B, LDB ) 00600 * 00601 END IF 00602 * 00603 END IF 00604 * 00605 END IF 00606 * 00607 END IF 00608 * 00609 ELSE 00610 * 00611 * SIDE = 'R' 00612 * 00613 * A is N-by-N. 00614 * If N is odd, set NISODD = .TRUE., and N1 and N2. 00615 * If N is even, NISODD = .FALSE., and K. 00616 * 00617 IF( MOD( N, 2 ).EQ.0 ) THEN 00618 NISODD = .FALSE. 00619 K = N / 2 00620 ELSE 00621 NISODD = .TRUE. 00622 IF( LOWER ) THEN 00623 N2 = N / 2 00624 N1 = N - N2 00625 ELSE 00626 N1 = N / 2 00627 N2 = N - N1 00628 END IF 00629 END IF 00630 * 00631 IF( NISODD ) THEN 00632 * 00633 * SIDE = 'R' and N is odd 00634 * 00635 IF( NORMALTRANSR ) THEN 00636 * 00637 * SIDE = 'R', N is odd, and TRANSR = 'N' 00638 * 00639 IF( LOWER ) THEN 00640 * 00641 * SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'L' 00642 * 00643 IF( NOTRANS ) THEN 00644 * 00645 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and 00646 * TRANS = 'N' 00647 * 00648 CALL STRSM( 'R', 'U', 'T', DIAG, M, N2, ALPHA, 00649 + A( N ), N, B( 0, N1 ), LDB ) 00650 CALL SGEMM( 'N', 'N', M, N1, N2, -ONE, B( 0, N1 ), 00651 + LDB, A( N1 ), N, ALPHA, B( 0, 0 ), 00652 + LDB ) 00653 CALL STRSM( 'R', 'L', 'N', DIAG, M, N1, ONE, 00654 + A( 0 ), N, B( 0, 0 ), LDB ) 00655 * 00656 ELSE 00657 * 00658 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and 00659 * TRANS = 'T' 00660 * 00661 CALL STRSM( 'R', 'L', 'T', DIAG, M, N1, ALPHA, 00662 + A( 0 ), N, B( 0, 0 ), LDB ) 00663 CALL SGEMM( 'N', 'T', M, N2, N1, -ONE, B( 0, 0 ), 00664 + LDB, A( N1 ), N, ALPHA, B( 0, N1 ), 00665 + LDB ) 00666 CALL STRSM( 'R', 'U', 'N', DIAG, M, N2, ONE, 00667 + A( N ), N, B( 0, N1 ), LDB ) 00668 * 00669 END IF 00670 * 00671 ELSE 00672 * 00673 * SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'U' 00674 * 00675 IF( NOTRANS ) THEN 00676 * 00677 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and 00678 * TRANS = 'N' 00679 * 00680 CALL STRSM( 'R', 'L', 'T', DIAG, M, N1, ALPHA, 00681 + A( N2 ), N, B( 0, 0 ), LDB ) 00682 CALL SGEMM( 'N', 'N', M, N2, N1, -ONE, B( 0, 0 ), 00683 + LDB, A( 0 ), N, ALPHA, B( 0, N1 ), 00684 + LDB ) 00685 CALL STRSM( 'R', 'U', 'N', DIAG, M, N2, ONE, 00686 + A( N1 ), N, B( 0, N1 ), LDB ) 00687 * 00688 ELSE 00689 * 00690 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and 00691 * TRANS = 'T' 00692 * 00693 CALL STRSM( 'R', 'U', 'T', DIAG, M, N2, ALPHA, 00694 + A( N1 ), N, B( 0, N1 ), LDB ) 00695 CALL SGEMM( 'N', 'T', M, N1, N2, -ONE, B( 0, N1 ), 00696 + LDB, A( 0 ), N, ALPHA, B( 0, 0 ), LDB ) 00697 CALL STRSM( 'R', 'L', 'N', DIAG, M, N1, ONE, 00698 + A( N2 ), N, B( 0, 0 ), LDB ) 00699 * 00700 END IF 00701 * 00702 END IF 00703 * 00704 ELSE 00705 * 00706 * SIDE = 'R', N is odd, and TRANSR = 'T' 00707 * 00708 IF( LOWER ) THEN 00709 * 00710 * SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'L' 00711 * 00712 IF( NOTRANS ) THEN 00713 * 00714 * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and 00715 * TRANS = 'N' 00716 * 00717 CALL STRSM( 'R', 'L', 'N', DIAG, M, N2, ALPHA, 00718 + A( 1 ), N1, B( 0, N1 ), LDB ) 00719 CALL SGEMM( 'N', 'T', M, N1, N2, -ONE, B( 0, N1 ), 00720 + LDB, A( N1*N1 ), N1, ALPHA, B( 0, 0 ), 00721 + LDB ) 00722 CALL STRSM( 'R', 'U', 'T', DIAG, M, N1, ONE, 00723 + A( 0 ), N1, B( 0, 0 ), LDB ) 00724 * 00725 ELSE 00726 * 00727 * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and 00728 * TRANS = 'T' 00729 * 00730 CALL STRSM( 'R', 'U', 'N', DIAG, M, N1, ALPHA, 00731 + A( 0 ), N1, B( 0, 0 ), LDB ) 00732 CALL SGEMM( 'N', 'N', M, N2, N1, -ONE, B( 0, 0 ), 00733 + LDB, A( N1*N1 ), N1, ALPHA, B( 0, N1 ), 00734 + LDB ) 00735 CALL STRSM( 'R', 'L', 'T', DIAG, M, N2, ONE, 00736 + A( 1 ), N1, B( 0, N1 ), LDB ) 00737 * 00738 END IF 00739 * 00740 ELSE 00741 * 00742 * SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'U' 00743 * 00744 IF( NOTRANS ) THEN 00745 * 00746 * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and 00747 * TRANS = 'N' 00748 * 00749 CALL STRSM( 'R', 'U', 'N', DIAG, M, N1, ALPHA, 00750 + A( N2*N2 ), N2, B( 0, 0 ), LDB ) 00751 CALL SGEMM( 'N', 'T', M, N2, N1, -ONE, B( 0, 0 ), 00752 + LDB, A( 0 ), N2, ALPHA, B( 0, N1 ), 00753 + LDB ) 00754 CALL STRSM( 'R', 'L', 'T', DIAG, M, N2, ONE, 00755 + A( N1*N2 ), N2, B( 0, N1 ), LDB ) 00756 * 00757 ELSE 00758 * 00759 * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and 00760 * TRANS = 'T' 00761 * 00762 CALL STRSM( 'R', 'L', 'N', DIAG, M, N2, ALPHA, 00763 + A( N1*N2 ), N2, B( 0, N1 ), LDB ) 00764 CALL SGEMM( 'N', 'N', M, N1, N2, -ONE, B( 0, N1 ), 00765 + LDB, A( 0 ), N2, ALPHA, B( 0, 0 ), 00766 + LDB ) 00767 CALL STRSM( 'R', 'U', 'T', DIAG, M, N1, ONE, 00768 + A( N2*N2 ), N2, B( 0, 0 ), LDB ) 00769 * 00770 END IF 00771 * 00772 END IF 00773 * 00774 END IF 00775 * 00776 ELSE 00777 * 00778 * SIDE = 'R' and N is even 00779 * 00780 IF( NORMALTRANSR ) THEN 00781 * 00782 * SIDE = 'R', N is even, and TRANSR = 'N' 00783 * 00784 IF( LOWER ) THEN 00785 * 00786 * SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'L' 00787 * 00788 IF( NOTRANS ) THEN 00789 * 00790 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L', 00791 * and TRANS = 'N' 00792 * 00793 CALL STRSM( 'R', 'U', 'T', DIAG, M, K, ALPHA, 00794 + A( 0 ), N+1, B( 0, K ), LDB ) 00795 CALL SGEMM( 'N', 'N', M, K, K, -ONE, B( 0, K ), 00796 + LDB, A( K+1 ), N+1, ALPHA, B( 0, 0 ), 00797 + LDB ) 00798 CALL STRSM( 'R', 'L', 'N', DIAG, M, K, ONE, 00799 + A( 1 ), N+1, B( 0, 0 ), LDB ) 00800 * 00801 ELSE 00802 * 00803 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L', 00804 * and TRANS = 'T' 00805 * 00806 CALL STRSM( 'R', 'L', 'T', DIAG, M, K, ALPHA, 00807 + A( 1 ), N+1, B( 0, 0 ), LDB ) 00808 CALL SGEMM( 'N', 'T', M, K, K, -ONE, B( 0, 0 ), 00809 + LDB, A( K+1 ), N+1, ALPHA, B( 0, K ), 00810 + LDB ) 00811 CALL STRSM( 'R', 'U', 'N', DIAG, M, K, ONE, 00812 + A( 0 ), N+1, B( 0, K ), LDB ) 00813 * 00814 END IF 00815 * 00816 ELSE 00817 * 00818 * SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'U' 00819 * 00820 IF( NOTRANS ) THEN 00821 * 00822 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U', 00823 * and TRANS = 'N' 00824 * 00825 CALL STRSM( 'R', 'L', 'T', DIAG, M, K, ALPHA, 00826 + A( K+1 ), N+1, B( 0, 0 ), LDB ) 00827 CALL SGEMM( 'N', 'N', M, K, K, -ONE, B( 0, 0 ), 00828 + LDB, A( 0 ), N+1, ALPHA, B( 0, K ), 00829 + LDB ) 00830 CALL STRSM( 'R', 'U', 'N', DIAG, M, K, ONE, 00831 + A( K ), N+1, B( 0, K ), LDB ) 00832 * 00833 ELSE 00834 * 00835 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U', 00836 * and TRANS = 'T' 00837 * 00838 CALL STRSM( 'R', 'U', 'T', DIAG, M, K, ALPHA, 00839 + A( K ), N+1, B( 0, K ), LDB ) 00840 CALL SGEMM( 'N', 'T', M, K, K, -ONE, B( 0, K ), 00841 + LDB, A( 0 ), N+1, ALPHA, B( 0, 0 ), 00842 + LDB ) 00843 CALL STRSM( 'R', 'L', 'N', DIAG, M, K, ONE, 00844 + A( K+1 ), N+1, B( 0, 0 ), LDB ) 00845 * 00846 END IF 00847 * 00848 END IF 00849 * 00850 ELSE 00851 * 00852 * SIDE = 'R', N is even, and TRANSR = 'T' 00853 * 00854 IF( LOWER ) THEN 00855 * 00856 * SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'L' 00857 * 00858 IF( NOTRANS ) THEN 00859 * 00860 * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L', 00861 * and TRANS = 'N' 00862 * 00863 CALL STRSM( 'R', 'L', 'N', DIAG, M, K, ALPHA, 00864 + A( 0 ), K, B( 0, K ), LDB ) 00865 CALL SGEMM( 'N', 'T', M, K, K, -ONE, B( 0, K ), 00866 + LDB, A( ( K+1 )*K ), K, ALPHA, 00867 + B( 0, 0 ), LDB ) 00868 CALL STRSM( 'R', 'U', 'T', DIAG, M, K, ONE, 00869 + A( K ), K, B( 0, 0 ), LDB ) 00870 * 00871 ELSE 00872 * 00873 * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L', 00874 * and TRANS = 'T' 00875 * 00876 CALL STRSM( 'R', 'U', 'N', DIAG, M, K, ALPHA, 00877 + A( K ), K, B( 0, 0 ), LDB ) 00878 CALL SGEMM( 'N', 'N', M, K, K, -ONE, B( 0, 0 ), 00879 + LDB, A( ( K+1 )*K ), K, ALPHA, 00880 + B( 0, K ), LDB ) 00881 CALL STRSM( 'R', 'L', 'T', DIAG, M, K, ONE, 00882 + A( 0 ), K, B( 0, K ), LDB ) 00883 * 00884 END IF 00885 * 00886 ELSE 00887 * 00888 * SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'U' 00889 * 00890 IF( NOTRANS ) THEN 00891 * 00892 * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U', 00893 * and TRANS = 'N' 00894 * 00895 CALL STRSM( 'R', 'U', 'N', DIAG, M, K, ALPHA, 00896 + A( ( K+1 )*K ), K, B( 0, 0 ), LDB ) 00897 CALL SGEMM( 'N', 'T', M, K, K, -ONE, B( 0, 0 ), 00898 + LDB, A( 0 ), K, ALPHA, B( 0, K ), LDB ) 00899 CALL STRSM( 'R', 'L', 'T', DIAG, M, K, ONE, 00900 + A( K*K ), K, B( 0, K ), LDB ) 00901 * 00902 ELSE 00903 * 00904 * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U', 00905 * and TRANS = 'T' 00906 * 00907 CALL STRSM( 'R', 'L', 'N', DIAG, M, K, ALPHA, 00908 + A( K*K ), K, B( 0, K ), LDB ) 00909 CALL SGEMM( 'N', 'N', M, K, K, -ONE, B( 0, K ), 00910 + LDB, A( 0 ), K, ALPHA, B( 0, 0 ), LDB ) 00911 CALL STRSM( 'R', 'U', 'T', DIAG, M, K, ONE, 00912 + A( ( K+1 )*K ), K, B( 0, 0 ), LDB ) 00913 * 00914 END IF 00915 * 00916 END IF 00917 * 00918 END IF 00919 * 00920 END IF 00921 END IF 00922 * 00923 RETURN 00924 * 00925 * End of STFSM 00926 * 00927 END
1.7.3 
