 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DTFSM( TRANSR, SIDE, UPLO, TRANS, DIAG, M, N, ALPHA, A, 00002 $ B, LDB ) 00003 * 00004 * -- LAPACK routine (version 3.3.1) -- 00005 * 00006 * -- Contributed by Fred Gustavson of the IBM Watson Research Center -- 00007 * -- April 2011 -- 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 DOUBLE PRECISION ALPHA 00017 * .. 00018 * .. Array Arguments .. 00019 DOUBLE PRECISION 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 * DTFSM 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**T. 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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 DOUBLE PRECISION ONE, ZERO 00221 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+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 XERBLA, DGEMM, DTRSM 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( 'DTFSM ', -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 * 00309 IF( MISODD ) THEN 00310 * 00311 * SIDE = 'L' and N is odd 00312 * 00313 IF( NORMALTRANSR ) THEN 00314 * 00315 * SIDE = 'L', N is odd, and TRANSR = 'N' 00316 * 00317 IF( LOWER ) THEN 00318 * 00319 * SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'L' 00320 * 00321 IF( NOTRANS ) THEN 00322 * 00323 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and 00324 * TRANS = 'N' 00325 * 00326 IF( M.EQ.1 ) THEN 00327 CALL DTRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA, 00328 $ A, M, B, LDB ) 00329 ELSE 00330 CALL DTRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA, 00331 $ A( 0 ), M, B, LDB ) 00332 CALL DGEMM( 'N', 'N', M2, N, M1, -ONE, A( M1 ), 00333 $ M, B, LDB, ALPHA, B( M1, 0 ), LDB ) 00334 CALL DTRSM( 'L', 'U', 'T', DIAG, M2, N, ONE, 00335 $ A( M ), M, B( M1, 0 ), LDB ) 00336 END IF 00337 * 00338 ELSE 00339 * 00340 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and 00341 * TRANS = 'T' 00342 * 00343 IF( M.EQ.1 ) THEN 00344 CALL DTRSM( 'L', 'L', 'T', DIAG, M1, N, ALPHA, 00345 $ A( 0 ), M, B, LDB ) 00346 ELSE 00347 CALL DTRSM( 'L', 'U', 'N', DIAG, M2, N, ALPHA, 00348 $ A( M ), M, B( M1, 0 ), LDB ) 00349 CALL DGEMM( 'T', 'N', M1, N, M2, -ONE, A( M1 ), 00350 $ M, B( M1, 0 ), LDB, ALPHA, B, LDB ) 00351 CALL DTRSM( 'L', 'L', 'T', DIAG, M1, N, ONE, 00352 $ A( 0 ), M, B, LDB ) 00353 END IF 00354 * 00355 END IF 00356 * 00357 ELSE 00358 * 00359 * SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'U' 00360 * 00361 IF( .NOT.NOTRANS ) THEN 00362 * 00363 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and 00364 * TRANS = 'N' 00365 * 00366 CALL DTRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA, 00367 $ A( M2 ), M, B, LDB ) 00368 CALL DGEMM( 'T', 'N', M2, N, M1, -ONE, A( 0 ), M, 00369 $ B, LDB, ALPHA, B( M1, 0 ), LDB ) 00370 CALL DTRSM( 'L', 'U', 'T', DIAG, M2, N, ONE, 00371 $ A( M1 ), M, B( M1, 0 ), LDB ) 00372 * 00373 ELSE 00374 * 00375 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and 00376 * TRANS = 'T' 00377 * 00378 CALL DTRSM( 'L', 'U', 'N', DIAG, M2, N, ALPHA, 00379 $ A( M1 ), M, B( M1, 0 ), LDB ) 00380 CALL DGEMM( 'N', 'N', M1, N, M2, -ONE, A( 0 ), M, 00381 $ B( M1, 0 ), LDB, ALPHA, B, LDB ) 00382 CALL DTRSM( 'L', 'L', 'T', DIAG, M1, N, ONE, 00383 $ A( M2 ), M, B, LDB ) 00384 * 00385 END IF 00386 * 00387 END IF 00388 * 00389 ELSE 00390 * 00391 * SIDE = 'L', N is odd, and TRANSR = 'T' 00392 * 00393 IF( LOWER ) THEN 00394 * 00395 * SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'L' 00396 * 00397 IF( NOTRANS ) THEN 00398 * 00399 * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and 00400 * TRANS = 'N' 00401 * 00402 IF( M.EQ.1 ) THEN 00403 CALL DTRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA, 00404 $ A( 0 ), M1, B, LDB ) 00405 ELSE 00406 CALL DTRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA, 00407 $ A( 0 ), M1, B, LDB ) 00408 CALL DGEMM( 'T', 'N', M2, N, M1, -ONE, 00409 $ A( M1*M1 ), M1, B, LDB, ALPHA, 00410 $ B( M1, 0 ), LDB ) 00411 CALL DTRSM( 'L', 'L', 'N', DIAG, M2, N, ONE, 00412 $ A( 1 ), M1, B( M1, 0 ), LDB ) 00413 END IF 00414 * 00415 ELSE 00416 * 00417 * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and 00418 * TRANS = 'T' 00419 * 00420 IF( M.EQ.1 ) THEN 00421 CALL DTRSM( 'L', 'U', 'N', DIAG, M1, N, ALPHA, 00422 $ A( 0 ), M1, B, LDB ) 00423 ELSE 00424 CALL DTRSM( 'L', 'L', 'T', DIAG, M2, N, ALPHA, 00425 $ A( 1 ), M1, B( M1, 0 ), LDB ) 00426 CALL DGEMM( 'N', 'N', M1, N, M2, -ONE, 00427 $ A( M1*M1 ), M1, B( M1, 0 ), LDB, 00428 $ ALPHA, B, LDB ) 00429 CALL DTRSM( 'L', 'U', 'N', DIAG, M1, N, ONE, 00430 $ A( 0 ), M1, B, LDB ) 00431 END IF 00432 * 00433 END IF 00434 * 00435 ELSE 00436 * 00437 * SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'U' 00438 * 00439 IF( .NOT.NOTRANS ) THEN 00440 * 00441 * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and 00442 * TRANS = 'N' 00443 * 00444 CALL DTRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA, 00445 $ A( M2*M2 ), M2, B, LDB ) 00446 CALL DGEMM( 'N', 'N', M2, N, M1, -ONE, A( 0 ), M2, 00447 $ B, LDB, ALPHA, B( M1, 0 ), LDB ) 00448 CALL DTRSM( 'L', 'L', 'N', DIAG, M2, N, ONE, 00449 $ A( M1*M2 ), M2, B( M1, 0 ), LDB ) 00450 * 00451 ELSE 00452 * 00453 * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and 00454 * TRANS = 'T' 00455 * 00456 CALL DTRSM( 'L', 'L', 'T', DIAG, M2, N, ALPHA, 00457 $ A( M1*M2 ), M2, B( M1, 0 ), LDB ) 00458 CALL DGEMM( 'T', 'N', M1, N, M2, -ONE, A( 0 ), M2, 00459 $ B( M1, 0 ), LDB, ALPHA, B, LDB ) 00460 CALL DTRSM( 'L', 'U', 'N', DIAG, M1, N, ONE, 00461 $ A( M2*M2 ), M2, B, LDB ) 00462 * 00463 END IF 00464 * 00465 END IF 00466 * 00467 END IF 00468 * 00469 ELSE 00470 * 00471 * SIDE = 'L' and N is even 00472 * 00473 IF( NORMALTRANSR ) THEN 00474 * 00475 * SIDE = 'L', N is even, and TRANSR = 'N' 00476 * 00477 IF( LOWER ) THEN 00478 * 00479 * SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'L' 00480 * 00481 IF( NOTRANS ) THEN 00482 * 00483 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L', 00484 * and TRANS = 'N' 00485 * 00486 CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ALPHA, 00487 $ A( 1 ), M+1, B, LDB ) 00488 CALL DGEMM( 'N', 'N', K, N, K, -ONE, A( K+1 ), 00489 $ M+1, B, LDB, ALPHA, B( K, 0 ), LDB ) 00490 CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ONE, 00491 $ A( 0 ), M+1, B( K, 0 ), LDB ) 00492 * 00493 ELSE 00494 * 00495 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L', 00496 * and TRANS = 'T' 00497 * 00498 CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ALPHA, 00499 $ A( 0 ), M+1, B( K, 0 ), LDB ) 00500 CALL DGEMM( 'T', 'N', K, N, K, -ONE, A( K+1 ), 00501 $ M+1, B( K, 0 ), LDB, ALPHA, B, LDB ) 00502 CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ONE, 00503 $ A( 1 ), M+1, B, LDB ) 00504 * 00505 END IF 00506 * 00507 ELSE 00508 * 00509 * SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'U' 00510 * 00511 IF( .NOT.NOTRANS ) THEN 00512 * 00513 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U', 00514 * and TRANS = 'N' 00515 * 00516 CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ALPHA, 00517 $ A( K+1 ), M+1, B, LDB ) 00518 CALL DGEMM( 'T', 'N', K, N, K, -ONE, A( 0 ), M+1, 00519 $ B, LDB, ALPHA, B( K, 0 ), LDB ) 00520 CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ONE, 00521 $ A( K ), M+1, B( K, 0 ), LDB ) 00522 * 00523 ELSE 00524 * 00525 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U', 00526 * and TRANS = 'T' 00527 CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ALPHA, 00528 $ A( K ), M+1, B( K, 0 ), LDB ) 00529 CALL DGEMM( 'N', 'N', K, N, K, -ONE, A( 0 ), M+1, 00530 $ B( K, 0 ), LDB, ALPHA, B, LDB ) 00531 CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ONE, 00532 $ A( K+1 ), M+1, B, LDB ) 00533 * 00534 END IF 00535 * 00536 END IF 00537 * 00538 ELSE 00539 * 00540 * SIDE = 'L', N is even, and TRANSR = 'T' 00541 * 00542 IF( LOWER ) THEN 00543 * 00544 * SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'L' 00545 * 00546 IF( NOTRANS ) THEN 00547 * 00548 * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L', 00549 * and TRANS = 'N' 00550 * 00551 CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ALPHA, 00552 $ A( K ), K, B, LDB ) 00553 CALL DGEMM( 'T', 'N', K, N, K, -ONE, 00554 $ A( K*( K+1 ) ), K, B, LDB, ALPHA, 00555 $ B( K, 0 ), LDB ) 00556 CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ONE, 00557 $ A( 0 ), K, B( K, 0 ), LDB ) 00558 * 00559 ELSE 00560 * 00561 * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L', 00562 * and TRANS = 'T' 00563 * 00564 CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ALPHA, 00565 $ A( 0 ), K, B( K, 0 ), LDB ) 00566 CALL DGEMM( 'N', 'N', K, N, K, -ONE, 00567 $ A( K*( K+1 ) ), K, B( K, 0 ), LDB, 00568 $ ALPHA, B, LDB ) 00569 CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ONE, 00570 $ A( K ), K, B, LDB ) 00571 * 00572 END IF 00573 * 00574 ELSE 00575 * 00576 * SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'U' 00577 * 00578 IF( .NOT.NOTRANS ) THEN 00579 * 00580 * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U', 00581 * and TRANS = 'N' 00582 * 00583 CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ALPHA, 00584 $ A( K*( K+1 ) ), K, B, LDB ) 00585 CALL DGEMM( 'N', 'N', K, N, K, -ONE, A( 0 ), K, B, 00586 $ LDB, ALPHA, B( K, 0 ), LDB ) 00587 CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ONE, 00588 $ A( K*K ), K, B( K, 0 ), LDB ) 00589 * 00590 ELSE 00591 * 00592 * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U', 00593 * and TRANS = 'T' 00594 * 00595 CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ALPHA, 00596 $ A( K*K ), K, B( K, 0 ), LDB ) 00597 CALL DGEMM( 'T', 'N', K, N, K, -ONE, A( 0 ), K, 00598 $ B( K, 0 ), LDB, ALPHA, B, LDB ) 00599 CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ONE, 00600 $ A( K*( K+1 ) ), K, B, LDB ) 00601 * 00602 END IF 00603 * 00604 END IF 00605 * 00606 END IF 00607 * 00608 END IF 00609 * 00610 ELSE 00611 * 00612 * SIDE = 'R' 00613 * 00614 * A is N-by-N. 00615 * If N is odd, set NISODD = .TRUE., and N1 and N2. 00616 * If N is even, NISODD = .FALSE., and K. 00617 * 00618 IF( MOD( N, 2 ).EQ.0 ) THEN 00619 NISODD = .FALSE. 00620 K = N / 2 00621 ELSE 00622 NISODD = .TRUE. 00623 IF( LOWER ) THEN 00624 N2 = N / 2 00625 N1 = N - N2 00626 ELSE 00627 N1 = N / 2 00628 N2 = N - N1 00629 END IF 00630 END IF 00631 * 00632 IF( NISODD ) THEN 00633 * 00634 * SIDE = 'R' and N is odd 00635 * 00636 IF( NORMALTRANSR ) THEN 00637 * 00638 * SIDE = 'R', N is odd, and TRANSR = 'N' 00639 * 00640 IF( LOWER ) THEN 00641 * 00642 * SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'L' 00643 * 00644 IF( NOTRANS ) THEN 00645 * 00646 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and 00647 * TRANS = 'N' 00648 * 00649 CALL DTRSM( 'R', 'U', 'T', DIAG, M, N2, ALPHA, 00650 $ A( N ), N, B( 0, N1 ), LDB ) 00651 CALL DGEMM( 'N', 'N', M, N1, N2, -ONE, B( 0, N1 ), 00652 $ LDB, A( N1 ), N, ALPHA, B( 0, 0 ), 00653 $ LDB ) 00654 CALL DTRSM( 'R', 'L', 'N', DIAG, M, N1, ONE, 00655 $ A( 0 ), N, B( 0, 0 ), LDB ) 00656 * 00657 ELSE 00658 * 00659 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and 00660 * TRANS = 'T' 00661 * 00662 CALL DTRSM( 'R', 'L', 'T', DIAG, M, N1, ALPHA, 00663 $ A( 0 ), N, B( 0, 0 ), LDB ) 00664 CALL DGEMM( 'N', 'T', M, N2, N1, -ONE, B( 0, 0 ), 00665 $ LDB, A( N1 ), N, ALPHA, B( 0, N1 ), 00666 $ LDB ) 00667 CALL DTRSM( 'R', 'U', 'N', DIAG, M, N2, ONE, 00668 $ A( N ), N, B( 0, N1 ), LDB ) 00669 * 00670 END IF 00671 * 00672 ELSE 00673 * 00674 * SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'U' 00675 * 00676 IF( NOTRANS ) THEN 00677 * 00678 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and 00679 * TRANS = 'N' 00680 * 00681 CALL DTRSM( 'R', 'L', 'T', DIAG, M, N1, ALPHA, 00682 $ A( N2 ), N, B( 0, 0 ), LDB ) 00683 CALL DGEMM( 'N', 'N', M, N2, N1, -ONE, B( 0, 0 ), 00684 $ LDB, A( 0 ), N, ALPHA, B( 0, N1 ), 00685 $ LDB ) 00686 CALL DTRSM( 'R', 'U', 'N', DIAG, M, N2, ONE, 00687 $ A( N1 ), N, B( 0, N1 ), LDB ) 00688 * 00689 ELSE 00690 * 00691 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and 00692 * TRANS = 'T' 00693 * 00694 CALL DTRSM( 'R', 'U', 'T', DIAG, M, N2, ALPHA, 00695 $ A( N1 ), N, B( 0, N1 ), LDB ) 00696 CALL DGEMM( 'N', 'T', M, N1, N2, -ONE, B( 0, N1 ), 00697 $ LDB, A( 0 ), N, ALPHA, B( 0, 0 ), LDB ) 00698 CALL DTRSM( 'R', 'L', 'N', DIAG, M, N1, ONE, 00699 $ A( N2 ), N, B( 0, 0 ), LDB ) 00700 * 00701 END IF 00702 * 00703 END IF 00704 * 00705 ELSE 00706 * 00707 * SIDE = 'R', N is odd, and TRANSR = 'T' 00708 * 00709 IF( LOWER ) THEN 00710 * 00711 * SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'L' 00712 * 00713 IF( NOTRANS ) THEN 00714 * 00715 * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and 00716 * TRANS = 'N' 00717 * 00718 CALL DTRSM( 'R', 'L', 'N', DIAG, M, N2, ALPHA, 00719 $ A( 1 ), N1, B( 0, N1 ), LDB ) 00720 CALL DGEMM( 'N', 'T', M, N1, N2, -ONE, B( 0, N1 ), 00721 $ LDB, A( N1*N1 ), N1, ALPHA, B( 0, 0 ), 00722 $ LDB ) 00723 CALL DTRSM( 'R', 'U', 'T', DIAG, M, N1, ONE, 00724 $ A( 0 ), N1, B( 0, 0 ), LDB ) 00725 * 00726 ELSE 00727 * 00728 * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and 00729 * TRANS = 'T' 00730 * 00731 CALL DTRSM( 'R', 'U', 'N', DIAG, M, N1, ALPHA, 00732 $ A( 0 ), N1, B( 0, 0 ), LDB ) 00733 CALL DGEMM( 'N', 'N', M, N2, N1, -ONE, B( 0, 0 ), 00734 $ LDB, A( N1*N1 ), N1, ALPHA, B( 0, N1 ), 00735 $ LDB ) 00736 CALL DTRSM( 'R', 'L', 'T', DIAG, M, N2, ONE, 00737 $ A( 1 ), N1, B( 0, N1 ), LDB ) 00738 * 00739 END IF 00740 * 00741 ELSE 00742 * 00743 * SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'U' 00744 * 00745 IF( NOTRANS ) THEN 00746 * 00747 * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and 00748 * TRANS = 'N' 00749 * 00750 CALL DTRSM( 'R', 'U', 'N', DIAG, M, N1, ALPHA, 00751 $ A( N2*N2 ), N2, B( 0, 0 ), LDB ) 00752 CALL DGEMM( 'N', 'T', M, N2, N1, -ONE, B( 0, 0 ), 00753 $ LDB, A( 0 ), N2, ALPHA, B( 0, N1 ), 00754 $ LDB ) 00755 CALL DTRSM( 'R', 'L', 'T', DIAG, M, N2, ONE, 00756 $ A( N1*N2 ), N2, B( 0, N1 ), LDB ) 00757 * 00758 ELSE 00759 * 00760 * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and 00761 * TRANS = 'T' 00762 * 00763 CALL DTRSM( 'R', 'L', 'N', DIAG, M, N2, ALPHA, 00764 $ A( N1*N2 ), N2, B( 0, N1 ), LDB ) 00765 CALL DGEMM( 'N', 'N', M, N1, N2, -ONE, B( 0, N1 ), 00766 $ LDB, A( 0 ), N2, ALPHA, B( 0, 0 ), 00767 $ LDB ) 00768 CALL DTRSM( 'R', 'U', 'T', DIAG, M, N1, ONE, 00769 $ A( N2*N2 ), N2, B( 0, 0 ), LDB ) 00770 * 00771 END IF 00772 * 00773 END IF 00774 * 00775 END IF 00776 * 00777 ELSE 00778 * 00779 * SIDE = 'R' and N is even 00780 * 00781 IF( NORMALTRANSR ) THEN 00782 * 00783 * SIDE = 'R', N is even, and TRANSR = 'N' 00784 * 00785 IF( LOWER ) THEN 00786 * 00787 * SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'L' 00788 * 00789 IF( NOTRANS ) THEN 00790 * 00791 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L', 00792 * and TRANS = 'N' 00793 * 00794 CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ALPHA, 00795 $ A( 0 ), N+1, B( 0, K ), LDB ) 00796 CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, K ), 00797 $ LDB, A( K+1 ), N+1, ALPHA, B( 0, 0 ), 00798 $ LDB ) 00799 CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ONE, 00800 $ A( 1 ), N+1, B( 0, 0 ), LDB ) 00801 * 00802 ELSE 00803 * 00804 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L', 00805 * and TRANS = 'T' 00806 * 00807 CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ALPHA, 00808 $ A( 1 ), N+1, B( 0, 0 ), LDB ) 00809 CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, 0 ), 00810 $ LDB, A( K+1 ), N+1, ALPHA, B( 0, K ), 00811 $ LDB ) 00812 CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ONE, 00813 $ A( 0 ), N+1, B( 0, K ), LDB ) 00814 * 00815 END IF 00816 * 00817 ELSE 00818 * 00819 * SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'U' 00820 * 00821 IF( NOTRANS ) THEN 00822 * 00823 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U', 00824 * and TRANS = 'N' 00825 * 00826 CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ALPHA, 00827 $ A( K+1 ), N+1, B( 0, 0 ), LDB ) 00828 CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, 0 ), 00829 $ LDB, A( 0 ), N+1, ALPHA, B( 0, K ), 00830 $ LDB ) 00831 CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ONE, 00832 $ A( K ), N+1, B( 0, K ), LDB ) 00833 * 00834 ELSE 00835 * 00836 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U', 00837 * and TRANS = 'T' 00838 * 00839 CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ALPHA, 00840 $ A( K ), N+1, B( 0, K ), LDB ) 00841 CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, K ), 00842 $ LDB, A( 0 ), N+1, ALPHA, B( 0, 0 ), 00843 $ LDB ) 00844 CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ONE, 00845 $ A( K+1 ), N+1, B( 0, 0 ), LDB ) 00846 * 00847 END IF 00848 * 00849 END IF 00850 * 00851 ELSE 00852 * 00853 * SIDE = 'R', N is even, and TRANSR = 'T' 00854 * 00855 IF( LOWER ) THEN 00856 * 00857 * SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'L' 00858 * 00859 IF( NOTRANS ) THEN 00860 * 00861 * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L', 00862 * and TRANS = 'N' 00863 * 00864 CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ALPHA, 00865 $ A( 0 ), K, B( 0, K ), LDB ) 00866 CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, K ), 00867 $ LDB, A( ( K+1 )*K ), K, ALPHA, 00868 $ B( 0, 0 ), LDB ) 00869 CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ONE, 00870 $ A( K ), K, B( 0, 0 ), LDB ) 00871 * 00872 ELSE 00873 * 00874 * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L', 00875 * and TRANS = 'T' 00876 * 00877 CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ALPHA, 00878 $ A( K ), K, B( 0, 0 ), LDB ) 00879 CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, 0 ), 00880 $ LDB, A( ( K+1 )*K ), K, ALPHA, 00881 $ B( 0, K ), LDB ) 00882 CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ONE, 00883 $ A( 0 ), K, B( 0, K ), LDB ) 00884 * 00885 END IF 00886 * 00887 ELSE 00888 * 00889 * SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'U' 00890 * 00891 IF( NOTRANS ) THEN 00892 * 00893 * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U', 00894 * and TRANS = 'N' 00895 * 00896 CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ALPHA, 00897 $ A( ( K+1 )*K ), K, B( 0, 0 ), LDB ) 00898 CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, 0 ), 00899 $ LDB, A( 0 ), K, ALPHA, B( 0, K ), LDB ) 00900 CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ONE, 00901 $ A( K*K ), K, B( 0, K ), LDB ) 00902 * 00903 ELSE 00904 * 00905 * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U', 00906 * and TRANS = 'T' 00907 * 00908 CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ALPHA, 00909 $ A( K*K ), K, B( 0, K ), LDB ) 00910 CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, K ), 00911 $ LDB, A( 0 ), K, ALPHA, B( 0, 0 ), LDB ) 00912 CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ONE, 00913 $ A( ( K+1 )*K ), K, B( 0, 0 ), LDB ) 00914 * 00915 END IF 00916 * 00917 END IF 00918 * 00919 END IF 00920 * 00921 END IF 00922 END IF 00923 * 00924 RETURN 00925 * 00926 * End of DTFSM 00927 * 00928 END
1.7.3 
