1
SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
4
* $Id: dtrsm.f 19695 2010-10-29 16:51:02Z d3y133 $
3
* =========== DOCUMENTATION ===========
5
* Online html documentation available at
6
* http://www.netlib.org/lapack/explore-html/
11
* SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
13
* .. Scalar Arguments ..
14
* DOUBLE PRECISION ALPHA
16
* CHARACTER DIAG,SIDE,TRANSA,UPLO
18
* .. Array Arguments ..
19
* DOUBLE PRECISION A(LDA,*),B(LDB,*)
28
*> DTRSM solves one of the matrix equations
30
*> op( A )*X = alpha*B, or X*op( A ) = alpha*B,
32
*> where alpha is a scalar, X and B are m by n matrices, A is a unit, or
33
*> non-unit, upper or lower triangular matrix and op( A ) is one of
35
*> op( A ) = A or op( A ) = A**T.
37
*> The matrix X is overwritten on B.
45
*> SIDE is CHARACTER*1
46
*> On entry, SIDE specifies whether op( A ) appears on the left
47
*> or right of X as follows:
49
*> SIDE = 'L' or 'l' op( A )*X = alpha*B.
51
*> SIDE = 'R' or 'r' X*op( A ) = alpha*B.
56
*> UPLO is CHARACTER*1
57
*> On entry, UPLO specifies whether the matrix A is an upper or
58
*> lower triangular matrix as follows:
60
*> UPLO = 'U' or 'u' A is an upper triangular matrix.
62
*> UPLO = 'L' or 'l' A is a lower triangular matrix.
67
*> TRANSA is CHARACTER*1
68
*> On entry, TRANSA specifies the form of op( A ) to be used in
69
*> the matrix multiplication as follows:
71
*> TRANSA = 'N' or 'n' op( A ) = A.
73
*> TRANSA = 'T' or 't' op( A ) = A**T.
75
*> TRANSA = 'C' or 'c' op( A ) = A**T.
80
*> DIAG is CHARACTER*1
81
*> On entry, DIAG specifies whether or not A is unit triangular
84
*> DIAG = 'U' or 'u' A is assumed to be unit triangular.
86
*> DIAG = 'N' or 'n' A is not assumed to be unit
93
*> On entry, M specifies the number of rows of B. M must be at
100
*> On entry, N specifies the number of columns of B. N must be
106
*> ALPHA is DOUBLE PRECISION.
107
*> On entry, ALPHA specifies the scalar alpha. When alpha is
108
*> zero then A is not referenced and B need not be set before
114
*> A is DOUBLE PRECISION array of DIMENSION ( LDA, k ),
115
*> where k is m when SIDE = 'L' or 'l'
116
*> and k is n when SIDE = 'R' or 'r'.
117
*> Before entry with UPLO = 'U' or 'u', the leading k by k
118
*> upper triangular part of the array A must contain the upper
119
*> triangular matrix and the strictly lower triangular part of
120
*> A is not referenced.
121
*> Before entry with UPLO = 'L' or 'l', the leading k by k
122
*> lower triangular part of the array A must contain the lower
123
*> triangular matrix and the strictly upper triangular part of
124
*> A is not referenced.
125
*> Note that when DIAG = 'U' or 'u', the diagonal elements of
126
*> A are not referenced either, but are assumed to be unity.
132
*> On entry, LDA specifies the first dimension of A as declared
133
*> in the calling (sub) program. When SIDE = 'L' or 'l' then
134
*> LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
135
*> then LDA must be at least max( 1, n ).
140
*> B is DOUBLE PRECISION array of DIMENSION ( LDB, n ).
141
*> Before entry, the leading m by n part of the array B must
142
*> contain the right-hand side matrix B, and on exit is
143
*> overwritten by the solution matrix X.
149
*> On entry, LDB specifies the first dimension of B as declared
150
*> in the calling (sub) program. LDB must be at least
157
*> \author Univ. of Tennessee
158
*> \author Univ. of California Berkeley
159
*> \author Univ. of Colorado Denver
162
*> \date November 2011
164
*> \ingroup double_blas_level3
166
*> \par Further Details:
167
* =====================
171
*> Level 3 Blas routine.
174
*> -- Written on 8-February-1989.
175
*> Jack Dongarra, Argonne National Laboratory.
176
*> Iain Duff, AERE Harwell.
177
*> Jeremy Du Croz, Numerical Algorithms Group Ltd.
178
*> Sven Hammarling, Numerical Algorithms Group Ltd.
181
* =====================================================================
182
SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
184
* -- Reference BLAS level3 routine (version 3.4.0) --
185
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
186
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6
189
* .. Scalar Arguments ..
7
CHARACTER*1 SIDE, UPLO, TRANSA, DIAG
190
DOUBLE PRECISION ALPHA
192
CHARACTER DIAG,SIDE,TRANSA,UPLO
10
194
* .. Array Arguments ..
11
DOUBLE PRECISION A( LDA, * ), B( LDB, * )
195
DOUBLE PRECISION A(LDA,*),B(LDB,*)
17
* DTRSM solves one of the matrix equations
19
* op( A )*X = alpha*B, or X*op( A ) = alpha*B,
21
* where alpha is a scalar, X and B are m by n matrices, A is a unit, or
22
* non-unit, upper or lower triangular matrix and op( A ) is one of
24
* op( A ) = A or op( A ) = A'.
26
* The matrix X is overwritten on B.
32
* On entry, SIDE specifies whether op( A ) appears on the left
33
* or right of X as follows:
35
* SIDE = 'L' or 'l' op( A )*X = alpha*B.
37
* SIDE = 'R' or 'r' X*op( A ) = alpha*B.
42
* On entry, UPLO specifies whether the matrix A is an upper or
43
* lower triangular matrix as follows:
45
* UPLO = 'U' or 'u' A is an upper triangular matrix.
47
* UPLO = 'L' or 'l' A is a lower triangular matrix.
51
* TRANSA - CHARACTER*1.
52
* On entry, TRANSA specifies the form of op( A ) to be used in
53
* the matrix multiplication as follows:
55
* TRANSA = 'N' or 'n' op( A ) = A.
57
* TRANSA = 'T' or 't' op( A ) = A'.
59
* TRANSA = 'C' or 'c' op( A ) = A'.
64
* On entry, DIAG specifies whether or not A is unit triangular
67
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
69
* DIAG = 'N' or 'n' A is not assumed to be unit
75
* On entry, M specifies the number of rows of B. M must be at
80
* On entry, N specifies the number of columns of B. N must be
84
* ALPHA - DOUBLE PRECISION.
85
* On entry, ALPHA specifies the scalar alpha. When alpha is
86
* zero then A is not referenced and B need not be set before
90
* A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
91
* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
92
* Before entry with UPLO = 'U' or 'u', the leading k by k
93
* upper triangular part of the array A must contain the upper
94
* triangular matrix and the strictly lower triangular part of
95
* A is not referenced.
96
* Before entry with UPLO = 'L' or 'l', the leading k by k
97
* lower triangular part of the array A must contain the lower
98
* triangular matrix and the strictly upper triangular part of
99
* A is not referenced.
100
* Note that when DIAG = 'U' or 'u', the diagonal elements of
101
* A are not referenced either, but are assumed to be unity.
105
* On entry, LDA specifies the first dimension of A as declared
106
* in the calling (sub) program. When SIDE = 'L' or 'l' then
107
* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
108
* then LDA must be at least max( 1, n ).
111
* B - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
112
* Before entry, the leading m by n part of the array B must
113
* contain the right-hand side matrix B, and on exit is
114
* overwritten by the solution matrix X.
117
* On entry, LDB specifies the first dimension of B as declared
118
* in the calling (sub) program. LDB must be at least
123
* Level 3 Blas routine.
126
* -- Written on 8-February-1989.
127
* Jack Dongarra, Argonne National Laboratory.
128
* Iain Duff, AERE Harwell.
129
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
130
* Sven Hammarling, Numerical Algorithms Group Ltd.
198
* =====================================================================
133
200
* .. External Functions ..
136
204
* .. External Subroutines ..
138
207
* .. Intrinsic Functions ..
140
210
* .. Local Scalars ..
141
LOGICAL LSIDE, NOUNIT, UPPER
142
INTEGER I, INFO, J, K, NROWA
143
DOUBLE PRECISION TEMP
211
DOUBLE PRECISION TEMP
212
INTEGER I,INFO,J,K,NROWA
213
LOGICAL LSIDE,NOUNIT,UPPER
144
215
* .. Parameters ..
145
DOUBLE PRECISION ONE , ZERO
146
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
216
DOUBLE PRECISION ONE,ZERO
217
PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
148
* .. Executable Statements ..
150
220
* Test the input parameters.
152
LSIDE = LSAME( SIDE , 'L' )
222
LSIDE = LSAME(SIDE,'L')
158
NOUNIT = LSAME( DIAG , 'N' )
159
UPPER = LSAME( UPLO , 'U' )
228
NOUNIT = LSAME(DIAG,'N')
229
UPPER = LSAME(UPLO,'U')
162
IF( ( .NOT.LSIDE ).AND.
163
$ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN
165
ELSE IF( ( .NOT.UPPER ).AND.
166
$ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN
168
ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
169
$ ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
170
$ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN
172
ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND.
173
$ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN
175
ELSE IF( M .LT.0 )THEN
177
ELSE IF( N .LT.0 )THEN
179
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
181
ELSE IF( LDB.LT.MAX( 1, M ) )THEN
232
IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
234
ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
236
ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND.
237
+ (.NOT.LSAME(TRANSA,'T')) .AND.
238
+ (.NOT.LSAME(TRANSA,'C'))) THEN
240
ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN
242
ELSE IF (M.LT.0) THEN
244
ELSE IF (N.LT.0) THEN
246
ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
248
ELSE IF (LDB.LT.MAX(1,M)) THEN
185
CALL XERBLA( 'DTRSM ', INFO )
252
CALL XERBLA('DTRSM ',INFO)
189
256
* Quick return if possible.
258
IF (M.EQ.0 .OR. N.EQ.0) RETURN
194
260
* And when alpha.eq.zero.
196
IF( ALPHA.EQ.ZERO )THEN
262
IF (ALPHA.EQ.ZERO) THEN
205
271
* Start the operations.
208
IF( LSAME( TRANSA, 'N' ) )THEN
274
IF (LSAME(TRANSA,'N')) THEN
210
276
* Form B := alpha*inv( A )*B.
214
IF( ALPHA.NE.ONE )THEN
216
B( I, J ) = ALPHA*B( I, J )
220
IF( B( K, J ).NE.ZERO )THEN
222
$ B( K, J ) = B( K, J )/A( K, K )
224
B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
231
IF( ALPHA.NE.ONE )THEN
233
B( I, J ) = ALPHA*B( I, J )
237
IF( B( K, J ).NE.ZERO )THEN
239
$ B( K, J ) = B( K, J )/A( K, K )
241
B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
249
* Form B := alpha*inv( A' )*B.
254
TEMP = ALPHA*B( I, J )
256
TEMP = TEMP - A( K, I )*B( K, J )
259
$ TEMP = TEMP/A( I, I )
266
TEMP = ALPHA*B( I, J )
268
TEMP = TEMP - A( K, I )*B( K, J )
271
$ TEMP = TEMP/A( I, I )
280
IF (ALPHA.NE.ONE) THEN
282
B(I,J) = ALPHA*B(I,J)
286
IF (B(K,J).NE.ZERO) THEN
287
IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
289
B(I,J) = B(I,J) - B(K,J)*A(I,K)
296
IF (ALPHA.NE.ONE) THEN
298
B(I,J) = ALPHA*B(I,J)
302
IF (B(K,J).NE.ZERO) THEN
303
IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
305
B(I,J) = B(I,J) - B(K,J)*A(I,K)
313
* Form B := alpha*inv( A**T )*B.
320
TEMP = TEMP - A(K,I)*B(K,J)
322
IF (NOUNIT) TEMP = TEMP/A(I,I)
331
TEMP = TEMP - A(K,I)*B(K,J)
333
IF (NOUNIT) TEMP = TEMP/A(I,I)
278
IF( LSAME( TRANSA, 'N' ) )THEN
340
IF (LSAME(TRANSA,'N')) THEN
280
342
* Form B := alpha*B*inv( A ).
284
IF( ALPHA.NE.ONE )THEN
286
B( I, J ) = ALPHA*B( I, J )
290
IF( A( K, J ).NE.ZERO )THEN
292
B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
299
B( I, J ) = TEMP*B( I, J )
305
IF( ALPHA.NE.ONE )THEN
307
B( I, J ) = ALPHA*B( I, J )
311
IF( A( K, J ).NE.ZERO )THEN
313
B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
320
B( I, J ) = TEMP*B( I, J )
327
* Form B := alpha*B*inv( A' ).
334
B( I, K ) = TEMP*B( I, K )
338
IF( A( J, K ).NE.ZERO )THEN
341
B( I, J ) = B( I, J ) - TEMP*B( I, K )
345
IF( ALPHA.NE.ONE )THEN
347
B( I, K ) = ALPHA*B( I, K )
356
B( I, K ) = TEMP*B( I, K )
360
IF( A( J, K ).NE.ZERO )THEN
363
B( I, J ) = B( I, J ) - TEMP*B( I, K )
367
IF( ALPHA.NE.ONE )THEN
369
B( I, K ) = ALPHA*B( I, K )
346
IF (ALPHA.NE.ONE) THEN
348
B(I,J) = ALPHA*B(I,J)
352
IF (A(K,J).NE.ZERO) THEN
354
B(I,J) = B(I,J) - A(K,J)*B(I,K)
367
IF (ALPHA.NE.ONE) THEN
369
B(I,J) = ALPHA*B(I,J)
373
IF (A(K,J).NE.ZERO) THEN
375
B(I,J) = B(I,J) - A(K,J)*B(I,K)
389
* Form B := alpha*B*inv( A**T ).
400
IF (A(J,K).NE.ZERO) THEN
403
B(I,J) = B(I,J) - TEMP*B(I,K)
407
IF (ALPHA.NE.ONE) THEN
409
B(I,K) = ALPHA*B(I,K)
422
IF (A(J,K).NE.ZERO) THEN
425
B(I,J) = B(I,J) - TEMP*B(I,K)
429
IF (ALPHA.NE.ONE) THEN
431
B(I,K) = ALPHA*B(I,K)