1 SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, 2 $ B, LDB ) 3* .. Scalar Arguments .. 4 CHARACTER*1 SIDE, UPLO, TRANSA, DIAG 5 INTEGER M, N, LDA, LDB 6 DOUBLE PRECISION ALPHA 7* .. Array Arguments .. 8 DOUBLE PRECISION A( LDA, * ), B( LDB, * ) 9* .. 10* 11* Purpose 12* ======= 13* 14* DTRSM solves one of the matrix equations 15* 16* op( A )*X = alpha*B, or X*op( A ) = alpha*B, 17* 18* where alpha is a scalar, X and B are m by n matrices, A is a unit, or 19* non-unit, upper or lower triangular matrix and op( A ) is one of 20* 21* op( A ) = A or op( A ) = A'. 22* 23* The matrix X is overwritten on B. 24* 25* Parameters 26* ========== 27* 28* SIDE - CHARACTER*1. 29* On entry, SIDE specifies whether op( A ) appears on the left 30* or right of X as follows: 31* 32* SIDE = 'L' or 'l' op( A )*X = alpha*B. 33* 34* SIDE = 'R' or 'r' X*op( A ) = alpha*B. 35* 36* Unchanged on exit. 37* 38* UPLO - CHARACTER*1. 39* On entry, UPLO specifies whether the matrix A is an upper or 40* lower triangular matrix as follows: 41* 42* UPLO = 'U' or 'u' A is an upper triangular matrix. 43* 44* UPLO = 'L' or 'l' A is a lower triangular matrix. 45* 46* Unchanged on exit. 47* 48* TRANSA - CHARACTER*1. 49* On entry, TRANSA specifies the form of op( A ) to be used in 50* the matrix multiplication as follows: 51* 52* TRANSA = 'N' or 'n' op( A ) = A. 53* 54* TRANSA = 'T' or 't' op( A ) = A'. 55* 56* TRANSA = 'C' or 'c' op( A ) = A'. 57* 58* Unchanged on exit. 59* 60* DIAG - CHARACTER*1. 61* On entry, DIAG specifies whether or not A is unit triangular 62* as follows: 63* 64* DIAG = 'U' or 'u' A is assumed to be unit triangular. 65* 66* DIAG = 'N' or 'n' A is not assumed to be unit 67* triangular. 68* 69* Unchanged on exit. 70* 71* M - INTEGER. 72* On entry, M specifies the number of rows of B. M must be at 73* least zero. 74* Unchanged on exit. 75* 76* N - INTEGER. 77* On entry, N specifies the number of columns of B. N must be 78* at least zero. 79* Unchanged on exit. 80* 81* ALPHA - DOUBLE PRECISION. 82* On entry, ALPHA specifies the scalar alpha. When alpha is 83* zero then A is not referenced and B need not be set before 84* entry. 85* Unchanged on exit. 86* 87* A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m 88* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. 89* Before entry with UPLO = 'U' or 'u', the leading k by k 90* upper triangular part of the array A must contain the upper 91* triangular matrix and the strictly lower triangular part of 92* A is not referenced. 93* Before entry with UPLO = 'L' or 'l', the leading k by k 94* lower triangular part of the array A must contain the lower 95* triangular matrix and the strictly upper triangular part of 96* A is not referenced. 97* Note that when DIAG = 'U' or 'u', the diagonal elements of 98* A are not referenced either, but are assumed to be unity. 99* Unchanged on exit. 100* 101* LDA - INTEGER. 102* On entry, LDA specifies the first dimension of A as declared 103* in the calling (sub) program. When SIDE = 'L' or 'l' then 104* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' 105* then LDA must be at least max( 1, n ). 106* Unchanged on exit. 107* 108* B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). 109* Before entry, the leading m by n part of the array B must 110* contain the right-hand side matrix B, and on exit is 111* overwritten by the solution matrix X. 112* 113* LDB - INTEGER. 114* On entry, LDB specifies the first dimension of B as declared 115* in the calling (sub) program. LDB must be at least 116* max( 1, m ). 117* Unchanged on exit. 118* 119* 120* Level 3 Blas routine. 121* 122* 123* -- Written on 8-February-1989. 124* Jack Dongarra, Argonne National Laboratory. 125* Iain Duff, AERE Harwell. 126* Jeremy Du Croz, Numerical Algorithms Group Ltd. 127* Sven Hammarling, Numerical Algorithms Group Ltd. 128* 129* 130* .. External Functions .. 131 LOGICAL LSAME 132 EXTERNAL LSAME 133* .. External Subroutines .. 134 EXTERNAL XERBLA 135* .. Intrinsic Functions .. 136 INTRINSIC MAX 137* .. Local Scalars .. 138 LOGICAL LSIDE, NOUNIT, UPPER 139 INTEGER I, INFO, J, K, NROWA 140 DOUBLE PRECISION TEMP 141* .. Parameters .. 142 DOUBLE PRECISION ONE , ZERO 143 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 144* .. 145* .. Executable Statements .. 146* 147* Test the input parameters. 148* 149 LSIDE = LSAME( SIDE , 'L' ) 150 IF( LSIDE )THEN 151 NROWA = M 152 ELSE 153 NROWA = N 154 END IF 155 NOUNIT = LSAME( DIAG , 'N' ) 156 UPPER = LSAME( UPLO , 'U' ) 157* 158 INFO = 0 159 IF( ( .NOT.LSIDE ).AND. 160 $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN 161 INFO = 1 162 ELSE IF( ( .NOT.UPPER ).AND. 163 $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN 164 INFO = 2 165 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND. 166 $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND. 167 $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN 168 INFO = 3 169 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND. 170 $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN 171 INFO = 4 172 ELSE IF( M .LT.0 )THEN 173 INFO = 5 174 ELSE IF( N .LT.0 )THEN 175 INFO = 6 176 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN 177 INFO = 9 178 ELSE IF( LDB.LT.MAX( 1, M ) )THEN 179 INFO = 11 180 END IF 181 IF( INFO.NE.0 )THEN 182 CALL XERBLA( 'DTRSM ', INFO ) 183 RETURN 184 END IF 185* 186* Quick return if possible. 187* 188 IF( N.EQ.0 ) 189 $ RETURN 190* 191* And when alpha.eq.zero. 192* 193 IF( ALPHA.EQ.ZERO )THEN 194 DO 20, J = 1, N 195 DO 10, I = 1, M 196 B( I, J ) = ZERO 197 10 CONTINUE 198 20 CONTINUE 199 RETURN 200 END IF 201* 202* Start the operations. 203* 204 IF( LSIDE )THEN 205 IF( LSAME( TRANSA, 'N' ) )THEN 206* 207* Form B := alpha*inv( A )*B. 208* 209 IF( UPPER )THEN 210 DO 60, J = 1, N 211 IF( ALPHA.NE.ONE )THEN 212 DO 30, I = 1, M 213 B( I, J ) = ALPHA*B( I, J ) 214 30 CONTINUE 215 END IF 216 DO 50, K = M, 1, -1 217 IF( B( K, J ).NE.ZERO )THEN 218 IF( NOUNIT ) 219 $ B( K, J ) = B( K, J )/A( K, K ) 220 DO 40, I = 1, K - 1 221 B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 222 40 CONTINUE 223 END IF 224 50 CONTINUE 225 60 CONTINUE 226 ELSE 227 DO 100, J = 1, N 228 IF( ALPHA.NE.ONE )THEN 229 DO 70, I = 1, M 230 B( I, J ) = ALPHA*B( I, J ) 231 70 CONTINUE 232 END IF 233 DO 90 K = 1, M 234 IF( B( K, J ).NE.ZERO )THEN 235 IF( NOUNIT ) 236 $ B( K, J ) = B( K, J )/A( K, K ) 237 DO 80, I = K + 1, M 238 B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 239 80 CONTINUE 240 END IF 241 90 CONTINUE 242 100 CONTINUE 243 END IF 244 ELSE 245* 246* Form B := alpha*inv( A' )*B. 247* 248 IF( UPPER )THEN 249 DO 130, J = 1, N 250 DO 120, I = 1, M 251 TEMP = ALPHA*B( I, J ) 252 DO 110, K = 1, I - 1 253 TEMP = TEMP - A( K, I )*B( K, J ) 254 110 CONTINUE 255 IF( NOUNIT ) 256 $ TEMP = TEMP/A( I, I ) 257 B( I, J ) = TEMP 258 120 CONTINUE 259 130 CONTINUE 260 ELSE 261 DO 160, J = 1, N 262 DO 150, I = M, 1, -1 263 TEMP = ALPHA*B( I, J ) 264 DO 140, K = I + 1, M 265 TEMP = TEMP - A( K, I )*B( K, J ) 266 140 CONTINUE 267 IF( NOUNIT ) 268 $ TEMP = TEMP/A( I, I ) 269 B( I, J ) = TEMP 270 150 CONTINUE 271 160 CONTINUE 272 END IF 273 END IF 274 ELSE 275 IF( LSAME( TRANSA, 'N' ) )THEN 276* 277* Form B := alpha*B*inv( A ). 278* 279 IF( UPPER )THEN 280 DO 210, J = 1, N 281 IF( ALPHA.NE.ONE )THEN 282 DO 170, I = 1, M 283 B( I, J ) = ALPHA*B( I, J ) 284 170 CONTINUE 285 END IF 286 DO 190, K = 1, J - 1 287 IF( A( K, J ).NE.ZERO )THEN 288 DO 180, I = 1, M 289 B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 290 180 CONTINUE 291 END IF 292 190 CONTINUE 293 IF( NOUNIT )THEN 294 TEMP = ONE/A( J, J ) 295 DO 200, I = 1, M 296 B( I, J ) = TEMP*B( I, J ) 297 200 CONTINUE 298 END IF 299 210 CONTINUE 300 ELSE 301 DO 260, J = N, 1, -1 302 IF( ALPHA.NE.ONE )THEN 303 DO 220, I = 1, M 304 B( I, J ) = ALPHA*B( I, J ) 305 220 CONTINUE 306 END IF 307 DO 240, K = J + 1, N 308 IF( A( K, J ).NE.ZERO )THEN 309 DO 230, I = 1, M 310 B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 311 230 CONTINUE 312 END IF 313 240 CONTINUE 314 IF( NOUNIT )THEN 315 TEMP = ONE/A( J, J ) 316 DO 250, I = 1, M 317 B( I, J ) = TEMP*B( I, J ) 318 250 CONTINUE 319 END IF 320 260 CONTINUE 321 END IF 322 ELSE 323* 324* Form B := alpha*B*inv( A' ). 325* 326 IF( UPPER )THEN 327 DO 310, K = N, 1, -1 328 IF( NOUNIT )THEN 329 TEMP = ONE/A( K, K ) 330 DO 270, I = 1, M 331 B( I, K ) = TEMP*B( I, K ) 332 270 CONTINUE 333 END IF 334 DO 290, J = 1, K - 1 335 IF( A( J, K ).NE.ZERO )THEN 336 TEMP = A( J, K ) 337 DO 280, I = 1, M 338 B( I, J ) = B( I, J ) - TEMP*B( I, K ) 339 280 CONTINUE 340 END IF 341 290 CONTINUE 342 IF( ALPHA.NE.ONE )THEN 343 DO 300, I = 1, M 344 B( I, K ) = ALPHA*B( I, K ) 345 300 CONTINUE 346 END IF 347 310 CONTINUE 348 ELSE 349 DO 360, K = 1, N 350 IF( NOUNIT )THEN 351 TEMP = ONE/A( K, K ) 352 DO 320, I = 1, M 353 B( I, K ) = TEMP*B( I, K ) 354 320 CONTINUE 355 END IF 356 DO 340, J = K + 1, N 357 IF( A( J, K ).NE.ZERO )THEN 358 TEMP = A( J, K ) 359 DO 330, I = 1, M 360 B( I, J ) = B( I, J ) - TEMP*B( I, K ) 361 330 CONTINUE 362 END IF 363 340 CONTINUE 364 IF( ALPHA.NE.ONE )THEN 365 DO 350, I = 1, M 366 B( I, K ) = ALPHA*B( I, K ) 367 350 CONTINUE 368 END IF 369 360 CONTINUE 370 END IF 371 END IF 372 END IF 373* 374 RETURN 375* 376* End of DTRSM . 377* 378 END 379