 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CTFSM( 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 COMPLEX ALPHA 00017 * .. 00018 * .. Array Arguments .. 00019 COMPLEX 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 * CTFSM 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**H. 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 * = 'C': The Conjugate-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 = 'C' or 'c' op( A ) = conjg( 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) COMPLEX 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) COMPLEX array, dimension (N*(N+1)/2) 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 = 'C' then RFP is the Conjugate-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 * conjugate-transpose Format. If UPLO = 'L' the RFP A contains 00112 * the NT elements of lower packed A either in normal or 00113 * conjugate-transpose Format. The LDA of RFP A is (N+1)/2 when 00114 * TRANSR = 'C'. 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) COMPLEX 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 Standard Packed Format when N is even. 00133 * 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 * conjugate-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 * conjugate-transpose of the last three columns of AP lower. 00152 * To denote conjugate we place -- above the element. This covers the 00153 * case N even and TRANSR = 'N'. 00154 * 00155 * RFP A RFP A 00156 * 00157 * -- -- -- 00158 * 03 04 05 33 43 53 00159 * -- -- 00160 * 13 14 15 00 44 54 00161 * -- 00162 * 23 24 25 10 11 55 00163 * 00164 * 33 34 35 20 21 22 00165 * -- 00166 * 00 44 45 30 31 32 00167 * -- -- 00168 * 01 11 55 40 41 42 00169 * -- -- -- 00170 * 02 12 22 50 51 52 00171 * 00172 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- 00173 * transpose of RFP A above. One therefore gets: 00174 * 00175 * 00176 * RFP A RFP A 00177 * 00178 * -- -- -- -- -- -- -- -- -- -- 00179 * 03 13 23 33 00 01 02 33 00 10 20 30 40 50 00180 * -- -- -- -- -- -- -- -- -- -- 00181 * 04 14 24 34 44 11 12 43 44 11 21 31 41 51 00182 * -- -- -- -- -- -- -- -- -- -- 00183 * 05 15 25 35 45 55 22 53 54 55 22 32 42 52 00184 * 00185 * 00186 * We next consider Standard Packed Format when N is odd. 00187 * We give an example where N = 5. 00188 * 00189 * AP is Upper AP is Lower 00190 * 00191 * 00 01 02 03 04 00 00192 * 11 12 13 14 10 11 00193 * 22 23 24 20 21 22 00194 * 33 34 30 31 32 33 00195 * 44 40 41 42 43 44 00196 * 00197 * 00198 * Let TRANSR = 'N'. RFP holds AP as follows: 00199 * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last 00200 * three columns of AP upper. The lower triangle A(3:4,0:1) consists of 00201 * conjugate-transpose of the first two columns of AP upper. 00202 * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first 00203 * three columns of AP lower. The upper triangle A(0:1,1:2) consists of 00204 * conjugate-transpose of the last two columns of AP lower. 00205 * To denote conjugate we place -- above the element. This covers the 00206 * case N odd and TRANSR = 'N'. 00207 * 00208 * RFP A RFP A 00209 * 00210 * -- -- 00211 * 02 03 04 00 33 43 00212 * -- 00213 * 12 13 14 10 11 44 00214 * 00215 * 22 23 24 20 21 22 00216 * -- 00217 * 00 33 34 30 31 32 00218 * -- -- 00219 * 01 11 44 40 41 42 00220 * 00221 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- 00222 * transpose of RFP A above. One therefore gets: 00223 * 00224 * 00225 * RFP A RFP A 00226 * 00227 * -- -- -- -- -- -- -- -- -- 00228 * 02 12 22 00 01 00 10 20 30 40 50 00229 * -- -- -- -- -- -- -- -- -- 00230 * 03 13 23 33 11 33 11 21 31 41 51 00231 * -- -- -- -- -- -- -- -- -- 00232 * 04 14 24 34 44 43 44 22 32 42 52 00233 * 00234 * ===================================================================== 00235 * .. 00236 * .. Parameters .. 00237 COMPLEX CONE, CZERO 00238 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ), 00239 $ CZERO = ( 0.0E+0, 0.0E+0 ) ) 00240 * .. 00241 * .. Local Scalars .. 00242 LOGICAL LOWER, LSIDE, MISODD, NISODD, NORMALTRANSR, 00243 $ NOTRANS 00244 INTEGER M1, M2, N1, N2, K, INFO, I, J 00245 * .. 00246 * .. External Functions .. 00247 LOGICAL LSAME 00248 EXTERNAL LSAME 00249 * .. 00250 * .. External Subroutines .. 00251 EXTERNAL XERBLA, CGEMM, CTRSM 00252 * .. 00253 * .. Intrinsic Functions .. 00254 INTRINSIC MAX, MOD 00255 * .. 00256 * .. Executable Statements .. 00257 * 00258 * Test the input parameters. 00259 * 00260 INFO = 0 00261 NORMALTRANSR = LSAME( TRANSR, 'N' ) 00262 LSIDE = LSAME( SIDE, 'L' ) 00263 LOWER = LSAME( UPLO, 'L' ) 00264 NOTRANS = LSAME( TRANS, 'N' ) 00265 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN 00266 INFO = -1 00267 ELSE IF( .NOT.LSIDE .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN 00268 INFO = -2 00269 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 00270 INFO = -3 00271 ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN 00272 INFO = -4 00273 ELSE IF( .NOT.LSAME( DIAG, 'N' ) .AND. .NOT.LSAME( DIAG, 'U' ) ) 00274 $ THEN 00275 INFO = -5 00276 ELSE IF( M.LT.0 ) THEN 00277 INFO = -6 00278 ELSE IF( N.LT.0 ) THEN 00279 INFO = -7 00280 ELSE IF( LDB.LT.MAX( 1, M ) ) THEN 00281 INFO = -11 00282 END IF 00283 IF( INFO.NE.0 ) THEN 00284 CALL XERBLA( 'CTFSM ', -INFO ) 00285 RETURN 00286 END IF 00287 * 00288 * Quick return when ( (N.EQ.0).OR.(M.EQ.0) ) 00289 * 00290 IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) ) 00291 $ RETURN 00292 * 00293 * Quick return when ALPHA.EQ.(0E+0,0E+0) 00294 * 00295 IF( ALPHA.EQ.CZERO ) THEN 00296 DO 20 J = 0, N - 1 00297 DO 10 I = 0, M - 1 00298 B( I, J ) = CZERO 00299 10 CONTINUE 00300 20 CONTINUE 00301 RETURN 00302 END IF 00303 * 00304 IF( LSIDE ) THEN 00305 * 00306 * SIDE = 'L' 00307 * 00308 * A is M-by-M. 00309 * If M is odd, set NISODD = .TRUE., and M1 and M2. 00310 * If M is even, NISODD = .FALSE., and M. 00311 * 00312 IF( MOD( M, 2 ).EQ.0 ) THEN 00313 MISODD = .FALSE. 00314 K = M / 2 00315 ELSE 00316 MISODD = .TRUE. 00317 IF( LOWER ) THEN 00318 M2 = M / 2 00319 M1 = M - M2 00320 ELSE 00321 M1 = M / 2 00322 M2 = M - M1 00323 END IF 00324 END IF 00325 * 00326 IF( MISODD ) THEN 00327 * 00328 * SIDE = 'L' and N is odd 00329 * 00330 IF( NORMALTRANSR ) THEN 00331 * 00332 * SIDE = 'L', N is odd, and TRANSR = 'N' 00333 * 00334 IF( LOWER ) THEN 00335 * 00336 * SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'L' 00337 * 00338 IF( NOTRANS ) THEN 00339 * 00340 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and 00341 * TRANS = 'N' 00342 * 00343 IF( M.EQ.1 ) THEN 00344 CALL CTRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA, 00345 $ A, M, B, LDB ) 00346 ELSE 00347 CALL CTRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA, 00348 $ A( 0 ), M, B, LDB ) 00349 CALL CGEMM( 'N', 'N', M2, N, M1, -CONE, A( M1 ), 00350 $ M, B, LDB, ALPHA, B( M1, 0 ), LDB ) 00351 CALL CTRSM( 'L', 'U', 'C', DIAG, M2, N, CONE, 00352 $ A( M ), M, B( M1, 0 ), LDB ) 00353 END IF 00354 * 00355 ELSE 00356 * 00357 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and 00358 * TRANS = 'C' 00359 * 00360 IF( M.EQ.1 ) THEN 00361 CALL CTRSM( 'L', 'L', 'C', DIAG, M1, N, ALPHA, 00362 $ A( 0 ), M, B, LDB ) 00363 ELSE 00364 CALL CTRSM( 'L', 'U', 'N', DIAG, M2, N, ALPHA, 00365 $ A( M ), M, B( M1, 0 ), LDB ) 00366 CALL CGEMM( 'C', 'N', M1, N, M2, -CONE, A( M1 ), 00367 $ M, B( M1, 0 ), LDB, ALPHA, B, LDB ) 00368 CALL CTRSM( 'L', 'L', 'C', DIAG, M1, N, CONE, 00369 $ A( 0 ), M, B, LDB ) 00370 END IF 00371 * 00372 END IF 00373 * 00374 ELSE 00375 * 00376 * SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'U' 00377 * 00378 IF( .NOT.NOTRANS ) THEN 00379 * 00380 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and 00381 * TRANS = 'N' 00382 * 00383 CALL CTRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA, 00384 $ A( M2 ), M, B, LDB ) 00385 CALL CGEMM( 'C', 'N', M2, N, M1, -CONE, A( 0 ), M, 00386 $ B, LDB, ALPHA, B( M1, 0 ), LDB ) 00387 CALL CTRSM( 'L', 'U', 'C', DIAG, M2, N, CONE, 00388 $ A( M1 ), M, B( M1, 0 ), LDB ) 00389 * 00390 ELSE 00391 * 00392 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and 00393 * TRANS = 'C' 00394 * 00395 CALL CTRSM( 'L', 'U', 'N', DIAG, M2, N, ALPHA, 00396 $ A( M1 ), M, B( M1, 0 ), LDB ) 00397 CALL CGEMM( 'N', 'N', M1, N, M2, -CONE, A( 0 ), M, 00398 $ B( M1, 0 ), LDB, ALPHA, B, LDB ) 00399 CALL CTRSM( 'L', 'L', 'C', DIAG, M1, N, CONE, 00400 $ A( M2 ), M, B, LDB ) 00401 * 00402 END IF 00403 * 00404 END IF 00405 * 00406 ELSE 00407 * 00408 * SIDE = 'L', N is odd, and TRANSR = 'C' 00409 * 00410 IF( LOWER ) THEN 00411 * 00412 * SIDE ='L', N is odd, TRANSR = 'C', and UPLO = 'L' 00413 * 00414 IF( NOTRANS ) THEN 00415 * 00416 * SIDE ='L', N is odd, TRANSR = 'C', UPLO = 'L', and 00417 * TRANS = 'N' 00418 * 00419 IF( M.EQ.1 ) THEN 00420 CALL CTRSM( 'L', 'U', 'C', DIAG, M1, N, ALPHA, 00421 $ A( 0 ), M1, B, LDB ) 00422 ELSE 00423 CALL CTRSM( 'L', 'U', 'C', DIAG, M1, N, ALPHA, 00424 $ A( 0 ), M1, B, LDB ) 00425 CALL CGEMM( 'C', 'N', M2, N, M1, -CONE, 00426 $ A( M1*M1 ), M1, B, LDB, ALPHA, 00427 $ B( M1, 0 ), LDB ) 00428 CALL CTRSM( 'L', 'L', 'N', DIAG, M2, N, CONE, 00429 $ A( 1 ), M1, B( M1, 0 ), LDB ) 00430 END IF 00431 * 00432 ELSE 00433 * 00434 * SIDE ='L', N is odd, TRANSR = 'C', UPLO = 'L', and 00435 * TRANS = 'C' 00436 * 00437 IF( M.EQ.1 ) THEN 00438 CALL CTRSM( 'L', 'U', 'N', DIAG, M1, N, ALPHA, 00439 $ A( 0 ), M1, B, LDB ) 00440 ELSE 00441 CALL CTRSM( 'L', 'L', 'C', DIAG, M2, N, ALPHA, 00442 $ A( 1 ), M1, B( M1, 0 ), LDB ) 00443 CALL CGEMM( 'N', 'N', M1, N, M2, -CONE, 00444 $ A( M1*M1 ), M1, B( M1, 0 ), LDB, 00445 $ ALPHA, B, LDB ) 00446 CALL CTRSM( 'L', 'U', 'N', DIAG, M1, N, CONE, 00447 $ A( 0 ), M1, B, LDB ) 00448 END IF 00449 * 00450 END IF 00451 * 00452 ELSE 00453 * 00454 * SIDE ='L', N is odd, TRANSR = 'C', and UPLO = 'U' 00455 * 00456 IF( .NOT.NOTRANS ) THEN 00457 * 00458 * SIDE ='L', N is odd, TRANSR = 'C', UPLO = 'U', and 00459 * TRANS = 'N' 00460 * 00461 CALL CTRSM( 'L', 'U', 'C', DIAG, M1, N, ALPHA, 00462 $ A( M2*M2 ), M2, B, LDB ) 00463 CALL CGEMM( 'N', 'N', M2, N, M1, -CONE, A( 0 ), M2, 00464 $ B, LDB, ALPHA, B( M1, 0 ), LDB ) 00465 CALL CTRSM( 'L', 'L', 'N', DIAG, M2, N, CONE, 00466 $ A( M1*M2 ), M2, B( M1, 0 ), LDB ) 00467 * 00468 ELSE 00469 * 00470 * SIDE ='L', N is odd, TRANSR = 'C', UPLO = 'U', and 00471 * TRANS = 'C' 00472 * 00473 CALL CTRSM( 'L', 'L', 'C', DIAG, M2, N, ALPHA, 00474 $ A( M1*M2 ), M2, B( M1, 0 ), LDB ) 00475 CALL CGEMM( 'C', 'N', M1, N, M2, -CONE, A( 0 ), M2, 00476 $ B( M1, 0 ), LDB, ALPHA, B, LDB ) 00477 CALL CTRSM( 'L', 'U', 'N', DIAG, M1, N, CONE, 00478 $ A( M2*M2 ), M2, B, LDB ) 00479 * 00480 END IF 00481 * 00482 END IF 00483 * 00484 END IF 00485 * 00486 ELSE 00487 * 00488 * SIDE = 'L' and N is even 00489 * 00490 IF( NORMALTRANSR ) THEN 00491 * 00492 * SIDE = 'L', N is even, and TRANSR = 'N' 00493 * 00494 IF( LOWER ) THEN 00495 * 00496 * SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'L' 00497 * 00498 IF( NOTRANS ) THEN 00499 * 00500 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L', 00501 * and TRANS = 'N' 00502 * 00503 CALL CTRSM( 'L', 'L', 'N', DIAG, K, N, ALPHA, 00504 $ A( 1 ), M+1, B, LDB ) 00505 CALL CGEMM( 'N', 'N', K, N, K, -CONE, A( K+1 ), 00506 $ M+1, B, LDB, ALPHA, B( K, 0 ), LDB ) 00507 CALL CTRSM( 'L', 'U', 'C', DIAG, K, N, CONE, 00508 $ A( 0 ), M+1, B( K, 0 ), LDB ) 00509 * 00510 ELSE 00511 * 00512 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L', 00513 * and TRANS = 'C' 00514 * 00515 CALL CTRSM( 'L', 'U', 'N', DIAG, K, N, ALPHA, 00516 $ A( 0 ), M+1, B( K, 0 ), LDB ) 00517 CALL CGEMM( 'C', 'N', K, N, K, -CONE, A( K+1 ), 00518 $ M+1, B( K, 0 ), LDB, ALPHA, B, LDB ) 00519 CALL CTRSM( 'L', 'L', 'C', DIAG, K, N, CONE, 00520 $ A( 1 ), M+1, B, LDB ) 00521 * 00522 END IF 00523 * 00524 ELSE 00525 * 00526 * SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'U' 00527 * 00528 IF( .NOT.NOTRANS ) THEN 00529 * 00530 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U', 00531 * and TRANS = 'N' 00532 * 00533 CALL CTRSM( 'L', 'L', 'N', DIAG, K, N, ALPHA, 00534 $ A( K+1 ), M+1, B, LDB ) 00535 CALL CGEMM( 'C', 'N', K, N, K, -CONE, A( 0 ), M+1, 00536 $ B, LDB, ALPHA, B( K, 0 ), LDB ) 00537 CALL CTRSM( 'L', 'U', 'C', DIAG, K, N, CONE, 00538 $ A( K ), M+1, B( K, 0 ), LDB ) 00539 * 00540 ELSE 00541 * 00542 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U', 00543 * and TRANS = 'C' 00544 CALL CTRSM( 'L', 'U', 'N', DIAG, K, N, ALPHA, 00545 $ A( K ), M+1, B( K, 0 ), LDB ) 00546 CALL CGEMM( 'N', 'N', K, N, K, -CONE, A( 0 ), M+1, 00547 $ B( K, 0 ), LDB, ALPHA, B, LDB ) 00548 CALL CTRSM( 'L', 'L', 'C', DIAG, K, N, CONE, 00549 $ A( K+1 ), M+1, B, LDB ) 00550 * 00551 END IF 00552 * 00553 END IF 00554 * 00555 ELSE 00556 * 00557 * SIDE = 'L', N is even, and TRANSR = 'C' 00558 * 00559 IF( LOWER ) THEN 00560 * 00561 * SIDE ='L', N is even, TRANSR = 'C', and UPLO = 'L' 00562 * 00563 IF( NOTRANS ) THEN 00564 * 00565 * SIDE ='L', N is even, TRANSR = 'C', UPLO = 'L', 00566 * and TRANS = 'N' 00567 * 00568 CALL CTRSM( 'L', 'U', 'C', DIAG, K, N, ALPHA, 00569 $ A( K ), K, B, LDB ) 00570 CALL CGEMM( 'C', 'N', K, N, K, -CONE, 00571 $ A( K*( K+1 ) ), K, B, LDB, ALPHA, 00572 $ B( K, 0 ), LDB ) 00573 CALL CTRSM( 'L', 'L', 'N', DIAG, K, N, CONE, 00574 $ A( 0 ), K, B( K, 0 ), LDB ) 00575 * 00576 ELSE 00577 * 00578 * SIDE ='L', N is even, TRANSR = 'C', UPLO = 'L', 00579 * and TRANS = 'C' 00580 * 00581 CALL CTRSM( 'L', 'L', 'C', DIAG, K, N, ALPHA, 00582 $ A( 0 ), K, B( K, 0 ), LDB ) 00583 CALL CGEMM( 'N', 'N', K, N, K, -CONE, 00584 $ A( K*( K+1 ) ), K, B( K, 0 ), LDB, 00585 $ ALPHA, B, LDB ) 00586 CALL CTRSM( 'L', 'U', 'N', DIAG, K, N, CONE, 00587 $ A( K ), K, B, LDB ) 00588 * 00589 END IF 00590 * 00591 ELSE 00592 * 00593 * SIDE ='L', N is even, TRANSR = 'C', and UPLO = 'U' 00594 * 00595 IF( .NOT.NOTRANS ) THEN 00596 * 00597 * SIDE ='L', N is even, TRANSR = 'C', UPLO = 'U', 00598 * and TRANS = 'N' 00599 * 00600 CALL CTRSM( 'L', 'U', 'C', DIAG, K, N, ALPHA, 00601 $ A( K*( K+1 ) ), K, B, LDB ) 00602 CALL CGEMM( 'N', 'N', K, N, K, -CONE, A( 0 ), K, B, 00603 $ LDB, ALPHA, B( K, 0 ), LDB ) 00604 CALL CTRSM( 'L', 'L', 'N', DIAG, K, N, CONE, 00605 $ A( K*K ), K, B( K, 0 ), LDB ) 00606 * 00607 ELSE 00608 * 00609 * SIDE ='L', N is even, TRANSR = 'C', UPLO = 'U', 00610 * and TRANS = 'C' 00611 * 00612 CALL CTRSM( 'L', 'L', 'C', DIAG, K, N, ALPHA, 00613 $ A( K*K ), K, B( K, 0 ), LDB ) 00614 CALL CGEMM( 'C', 'N', K, N, K, -CONE, A( 0 ), K, 00615 $ B( K, 0 ), LDB, ALPHA, B, LDB ) 00616 CALL CTRSM( 'L', 'U', 'N', DIAG, K, N, CONE, 00617 $ A( K*( K+1 ) ), K, B, LDB ) 00618 * 00619 END IF 00620 * 00621 END IF 00622 * 00623 END IF 00624 * 00625 END IF 00626 * 00627 ELSE 00628 * 00629 * SIDE = 'R' 00630 * 00631 * A is N-by-N. 00632 * If N is odd, set NISODD = .TRUE., and N1 and N2. 00633 * If N is even, NISODD = .FALSE., and K. 00634 * 00635 IF( MOD( N, 2 ).EQ.0 ) THEN 00636 NISODD = .FALSE. 00637 K = N / 2 00638 ELSE 00639 NISODD = .TRUE. 00640 IF( LOWER ) THEN 00641 N2 = N / 2 00642 N1 = N - N2 00643 ELSE 00644 N1 = N / 2 00645 N2 = N - N1 00646 END IF 00647 END IF 00648 * 00649 IF( NISODD ) THEN 00650 * 00651 * SIDE = 'R' and N is odd 00652 * 00653 IF( NORMALTRANSR ) THEN 00654 * 00655 * SIDE = 'R', N is odd, and TRANSR = 'N' 00656 * 00657 IF( LOWER ) THEN 00658 * 00659 * SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'L' 00660 * 00661 IF( NOTRANS ) THEN 00662 * 00663 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and 00664 * TRANS = 'N' 00665 * 00666 CALL CTRSM( 'R', 'U', 'C', DIAG, M, N2, ALPHA, 00667 $ A( N ), N, B( 0, N1 ), LDB ) 00668 CALL CGEMM( 'N', 'N', M, N1, N2, -CONE, B( 0, N1 ), 00669 $ LDB, A( N1 ), N, ALPHA, B( 0, 0 ), 00670 $ LDB ) 00671 CALL CTRSM( 'R', 'L', 'N', DIAG, M, N1, CONE, 00672 $ A( 0 ), N, B( 0, 0 ), LDB ) 00673 * 00674 ELSE 00675 * 00676 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and 00677 * TRANS = 'C' 00678 * 00679 CALL CTRSM( 'R', 'L', 'C', DIAG, M, N1, ALPHA, 00680 $ A( 0 ), N, B( 0, 0 ), LDB ) 00681 CALL CGEMM( 'N', 'C', M, N2, N1, -CONE, B( 0, 0 ), 00682 $ LDB, A( N1 ), N, ALPHA, B( 0, N1 ), 00683 $ LDB ) 00684 CALL CTRSM( 'R', 'U', 'N', DIAG, M, N2, CONE, 00685 $ A( N ), N, B( 0, N1 ), LDB ) 00686 * 00687 END IF 00688 * 00689 ELSE 00690 * 00691 * SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'U' 00692 * 00693 IF( NOTRANS ) THEN 00694 * 00695 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and 00696 * TRANS = 'N' 00697 * 00698 CALL CTRSM( 'R', 'L', 'C', DIAG, M, N1, ALPHA, 00699 $ A( N2 ), N, B( 0, 0 ), LDB ) 00700 CALL CGEMM( 'N', 'N', M, N2, N1, -CONE, B( 0, 0 ), 00701 $ LDB, A( 0 ), N, ALPHA, B( 0, N1 ), 00702 $ LDB ) 00703 CALL CTRSM( 'R', 'U', 'N', DIAG, M, N2, CONE, 00704 $ A( N1 ), N, B( 0, N1 ), LDB ) 00705 * 00706 ELSE 00707 * 00708 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and 00709 * TRANS = 'C' 00710 * 00711 CALL CTRSM( 'R', 'U', 'C', DIAG, M, N2, ALPHA, 00712 $ A( N1 ), N, B( 0, N1 ), LDB ) 00713 CALL CGEMM( 'N', 'C', M, N1, N2, -CONE, B( 0, N1 ), 00714 $ LDB, A( 0 ), N, ALPHA, B( 0, 0 ), LDB ) 00715 CALL CTRSM( 'R', 'L', 'N', DIAG, M, N1, CONE, 00716 $ A( N2 ), N, B( 0, 0 ), LDB ) 00717 * 00718 END IF 00719 * 00720 END IF 00721 * 00722 ELSE 00723 * 00724 * SIDE = 'R', N is odd, and TRANSR = 'C' 00725 * 00726 IF( LOWER ) THEN 00727 * 00728 * SIDE ='R', N is odd, TRANSR = 'C', and UPLO = 'L' 00729 * 00730 IF( NOTRANS ) THEN 00731 * 00732 * SIDE ='R', N is odd, TRANSR = 'C', UPLO = 'L', and 00733 * TRANS = 'N' 00734 * 00735 CALL CTRSM( 'R', 'L', 'N', DIAG, M, N2, ALPHA, 00736 $ A( 1 ), N1, B( 0, N1 ), LDB ) 00737 CALL CGEMM( 'N', 'C', M, N1, N2, -CONE, B( 0, N1 ), 00738 $ LDB, A( N1*N1 ), N1, ALPHA, B( 0, 0 ), 00739 $ LDB ) 00740 CALL CTRSM( 'R', 'U', 'C', DIAG, M, N1, CONE, 00741 $ A( 0 ), N1, B( 0, 0 ), LDB ) 00742 * 00743 ELSE 00744 * 00745 * SIDE ='R', N is odd, TRANSR = 'C', UPLO = 'L', and 00746 * TRANS = 'C' 00747 * 00748 CALL CTRSM( 'R', 'U', 'N', DIAG, M, N1, ALPHA, 00749 $ A( 0 ), N1, B( 0, 0 ), LDB ) 00750 CALL CGEMM( 'N', 'N', M, N2, N1, -CONE, B( 0, 0 ), 00751 $ LDB, A( N1*N1 ), N1, ALPHA, B( 0, N1 ), 00752 $ LDB ) 00753 CALL CTRSM( 'R', 'L', 'C', DIAG, M, N2, CONE, 00754 $ A( 1 ), N1, B( 0, N1 ), LDB ) 00755 * 00756 END IF 00757 * 00758 ELSE 00759 * 00760 * SIDE ='R', N is odd, TRANSR = 'C', and UPLO = 'U' 00761 * 00762 IF( NOTRANS ) THEN 00763 * 00764 * SIDE ='R', N is odd, TRANSR = 'C', UPLO = 'U', and 00765 * TRANS = 'N' 00766 * 00767 CALL CTRSM( 'R', 'U', 'N', DIAG, M, N1, ALPHA, 00768 $ A( N2*N2 ), N2, B( 0, 0 ), LDB ) 00769 CALL CGEMM( 'N', 'C', M, N2, N1, -CONE, B( 0, 0 ), 00770 $ LDB, A( 0 ), N2, ALPHA, B( 0, N1 ), 00771 $ LDB ) 00772 CALL CTRSM( 'R', 'L', 'C', DIAG, M, N2, CONE, 00773 $ A( N1*N2 ), N2, B( 0, N1 ), LDB ) 00774 * 00775 ELSE 00776 * 00777 * SIDE ='R', N is odd, TRANSR = 'C', UPLO = 'U', and 00778 * TRANS = 'C' 00779 * 00780 CALL CTRSM( 'R', 'L', 'N', DIAG, M, N2, ALPHA, 00781 $ A( N1*N2 ), N2, B( 0, N1 ), LDB ) 00782 CALL CGEMM( 'N', 'N', M, N1, N2, -CONE, B( 0, N1 ), 00783 $ LDB, A( 0 ), N2, ALPHA, B( 0, 0 ), 00784 $ LDB ) 00785 CALL CTRSM( 'R', 'U', 'C', DIAG, M, N1, CONE, 00786 $ A( N2*N2 ), N2, B( 0, 0 ), LDB ) 00787 * 00788 END IF 00789 * 00790 END IF 00791 * 00792 END IF 00793 * 00794 ELSE 00795 * 00796 * SIDE = 'R' and N is even 00797 * 00798 IF( NORMALTRANSR ) THEN 00799 * 00800 * SIDE = 'R', N is even, and TRANSR = 'N' 00801 * 00802 IF( LOWER ) THEN 00803 * 00804 * SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'L' 00805 * 00806 IF( NOTRANS ) THEN 00807 * 00808 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L', 00809 * and TRANS = 'N' 00810 * 00811 CALL CTRSM( 'R', 'U', 'C', DIAG, M, K, ALPHA, 00812 $ A( 0 ), N+1, B( 0, K ), LDB ) 00813 CALL CGEMM( 'N', 'N', M, K, K, -CONE, B( 0, K ), 00814 $ LDB, A( K+1 ), N+1, ALPHA, B( 0, 0 ), 00815 $ LDB ) 00816 CALL CTRSM( 'R', 'L', 'N', DIAG, M, K, CONE, 00817 $ A( 1 ), N+1, B( 0, 0 ), LDB ) 00818 * 00819 ELSE 00820 * 00821 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L', 00822 * and TRANS = 'C' 00823 * 00824 CALL CTRSM( 'R', 'L', 'C', DIAG, M, K, ALPHA, 00825 $ A( 1 ), N+1, B( 0, 0 ), LDB ) 00826 CALL CGEMM( 'N', 'C', M, K, K, -CONE, B( 0, 0 ), 00827 $ LDB, A( K+1 ), N+1, ALPHA, B( 0, K ), 00828 $ LDB ) 00829 CALL CTRSM( 'R', 'U', 'N', DIAG, M, K, CONE, 00830 $ A( 0 ), N+1, B( 0, K ), LDB ) 00831 * 00832 END IF 00833 * 00834 ELSE 00835 * 00836 * SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'U' 00837 * 00838 IF( NOTRANS ) THEN 00839 * 00840 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U', 00841 * and TRANS = 'N' 00842 * 00843 CALL CTRSM( 'R', 'L', 'C', DIAG, M, K, ALPHA, 00844 $ A( K+1 ), N+1, B( 0, 0 ), LDB ) 00845 CALL CGEMM( 'N', 'N', M, K, K, -CONE, B( 0, 0 ), 00846 $ LDB, A( 0 ), N+1, ALPHA, B( 0, K ), 00847 $ LDB ) 00848 CALL CTRSM( 'R', 'U', 'N', DIAG, M, K, CONE, 00849 $ A( K ), N+1, B( 0, K ), LDB ) 00850 * 00851 ELSE 00852 * 00853 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U', 00854 * and TRANS = 'C' 00855 * 00856 CALL CTRSM( 'R', 'U', 'C', DIAG, M, K, ALPHA, 00857 $ A( K ), N+1, B( 0, K ), LDB ) 00858 CALL CGEMM( 'N', 'C', M, K, K, -CONE, B( 0, K ), 00859 $ LDB, A( 0 ), N+1, ALPHA, B( 0, 0 ), 00860 $ LDB ) 00861 CALL CTRSM( 'R', 'L', 'N', DIAG, M, K, CONE, 00862 $ A( K+1 ), N+1, B( 0, 0 ), LDB ) 00863 * 00864 END IF 00865 * 00866 END IF 00867 * 00868 ELSE 00869 * 00870 * SIDE = 'R', N is even, and TRANSR = 'C' 00871 * 00872 IF( LOWER ) THEN 00873 * 00874 * SIDE ='R', N is even, TRANSR = 'C', and UPLO = 'L' 00875 * 00876 IF( NOTRANS ) THEN 00877 * 00878 * SIDE ='R', N is even, TRANSR = 'C', UPLO = 'L', 00879 * and TRANS = 'N' 00880 * 00881 CALL CTRSM( 'R', 'L', 'N', DIAG, M, K, ALPHA, 00882 $ A( 0 ), K, B( 0, K ), LDB ) 00883 CALL CGEMM( 'N', 'C', M, K, K, -CONE, B( 0, K ), 00884 $ LDB, A( ( K+1 )*K ), K, ALPHA, 00885 $ B( 0, 0 ), LDB ) 00886 CALL CTRSM( 'R', 'U', 'C', DIAG, M, K, CONE, 00887 $ A( K ), K, B( 0, 0 ), LDB ) 00888 * 00889 ELSE 00890 * 00891 * SIDE ='R', N is even, TRANSR = 'C', UPLO = 'L', 00892 * and TRANS = 'C' 00893 * 00894 CALL CTRSM( 'R', 'U', 'N', DIAG, M, K, ALPHA, 00895 $ A( K ), K, B( 0, 0 ), LDB ) 00896 CALL CGEMM( 'N', 'N', M, K, K, -CONE, B( 0, 0 ), 00897 $ LDB, A( ( K+1 )*K ), K, ALPHA, 00898 $ B( 0, K ), LDB ) 00899 CALL CTRSM( 'R', 'L', 'C', DIAG, M, K, CONE, 00900 $ A( 0 ), K, B( 0, K ), LDB ) 00901 * 00902 END IF 00903 * 00904 ELSE 00905 * 00906 * SIDE ='R', N is even, TRANSR = 'C', and UPLO = 'U' 00907 * 00908 IF( NOTRANS ) THEN 00909 * 00910 * SIDE ='R', N is even, TRANSR = 'C', UPLO = 'U', 00911 * and TRANS = 'N' 00912 * 00913 CALL CTRSM( 'R', 'U', 'N', DIAG, M, K, ALPHA, 00914 $ A( ( K+1 )*K ), K, B( 0, 0 ), LDB ) 00915 CALL CGEMM( 'N', 'C', M, K, K, -CONE, B( 0, 0 ), 00916 $ LDB, A( 0 ), K, ALPHA, B( 0, K ), LDB ) 00917 CALL CTRSM( 'R', 'L', 'C', DIAG, M, K, CONE, 00918 $ A( K*K ), K, B( 0, K ), LDB ) 00919 * 00920 ELSE 00921 * 00922 * SIDE ='R', N is even, TRANSR = 'C', UPLO = 'U', 00923 * and TRANS = 'C' 00924 * 00925 CALL CTRSM( 'R', 'L', 'N', DIAG, M, K, ALPHA, 00926 $ A( K*K ), K, B( 0, K ), LDB ) 00927 CALL CGEMM( 'N', 'N', M, K, K, -CONE, B( 0, K ), 00928 $ LDB, A( 0 ), K, ALPHA, B( 0, 0 ), LDB ) 00929 CALL CTRSM( 'R', 'U', 'C', DIAG, M, K, CONE, 00930 $ A( ( K+1 )*K ), K, B( 0, 0 ), LDB ) 00931 * 00932 END IF 00933 * 00934 END IF 00935 * 00936 END IF 00937 * 00938 END IF 00939 END IF 00940 * 00941 RETURN 00942 * 00943 * End of CTFSM 00944 * 00945 END
1.7.3 
