1
SUBROUTINE ZTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
4
* $Id: ztrsm.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 ZTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
13
* .. Scalar Arguments ..
16
* CHARACTER DIAG,SIDE,TRANSA,UPLO
18
* .. Array Arguments ..
19
* COMPLEX*16 A(LDA,*),B(LDB,*)
28
*> ZTRSM 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 or op( A ) = A**H.
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**H.
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 COMPLEX*16
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 COMPLEX*16 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 COMPLEX*16 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 complex16_blas_level3
166
*> \par Further Details:
167
* =====================
171
*> Level 3 Blas routine.
173
*> -- Written on 8-February-1989.
174
*> Jack Dongarra, Argonne National Laboratory.
175
*> Iain Duff, AERE Harwell.
176
*> Jeremy Du Croz, Numerical Algorithms Group Ltd.
177
*> Sven Hammarling, Numerical Algorithms Group Ltd.
180
* =====================================================================
181
SUBROUTINE ZTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
183
* -- Reference BLAS level3 routine (version 3.4.0) --
184
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
185
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6
188
* .. Scalar Arguments ..
7
CHARACTER*1 SIDE, UPLO, TRANSA, DIAG
191
CHARACTER DIAG,SIDE,TRANSA,UPLO
10
193
* .. Array Arguments ..
11
COMPLEX*16 A( LDA, * ), B( LDB, * )
194
COMPLEX*16 A(LDA,*),B(LDB,*)
17
* ZTRSM 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' or op( A ) = conjg( 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 ) = conjg( 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 - COMPLEX*16 .
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 - COMPLEX*16 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 - COMPLEX*16 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.
125
* -- Written on 8-February-1989.
126
* Jack Dongarra, Argonne National Laboratory.
127
* Iain Duff, AERE Harwell.
128
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
129
* Sven Hammarling, Numerical Algorithms Group Ltd.
197
* =====================================================================
132
199
* .. External Functions ..
135
203
* .. External Subroutines ..
137
206
* .. Intrinsic Functions ..
138
INTRINSIC DCONJG, MAX
139
209
* .. Local Scalars ..
140
LOGICAL LSIDE, NOCONJ, NOUNIT, UPPER
141
INTEGER I, INFO, J, K, NROWA
211
INTEGER I,INFO,J,K,NROWA
212
LOGICAL LSIDE,NOCONJ,NOUNIT,UPPER
143
214
* .. Parameters ..
145
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
147
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
216
PARAMETER (ONE= (1.0D+0,0.0D+0))
218
PARAMETER (ZERO= (0.0D+0,0.0D+0))
149
* .. Executable Statements ..
151
221
* Test the input parameters.
153
LSIDE = LSAME( SIDE , 'L' )
223
LSIDE = LSAME(SIDE,'L')
159
NOCONJ = LSAME( TRANSA, 'T' )
160
NOUNIT = LSAME( DIAG , 'N' )
161
UPPER = LSAME( UPLO , 'U' )
229
NOCONJ = LSAME(TRANSA,'T')
230
NOUNIT = LSAME(DIAG,'N')
231
UPPER = LSAME(UPLO,'U')
164
IF( ( .NOT.LSIDE ).AND.
165
$ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN
167
ELSE IF( ( .NOT.UPPER ).AND.
168
$ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN
170
ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
171
$ ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
172
$ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN
174
ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND.
175
$ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN
177
ELSE IF( M .LT.0 )THEN
179
ELSE IF( N .LT.0 )THEN
181
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
183
ELSE IF( LDB.LT.MAX( 1, M ) )THEN
234
IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
236
ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
238
ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND.
239
+ (.NOT.LSAME(TRANSA,'T')) .AND.
240
+ (.NOT.LSAME(TRANSA,'C'))) THEN
242
ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN
244
ELSE IF (M.LT.0) THEN
246
ELSE IF (N.LT.0) THEN
248
ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
250
ELSE IF (LDB.LT.MAX(1,M)) THEN
187
CALL XERBLA( 'ZTRSM ', INFO )
254
CALL XERBLA('ZTRSM ',INFO)
191
258
* Quick return if possible.
260
IF (M.EQ.0 .OR. N.EQ.0) RETURN
196
262
* And when alpha.eq.zero.
198
IF( ALPHA.EQ.ZERO )THEN
264
IF (ALPHA.EQ.ZERO) THEN
207
273
* Start the operations.
210
IF( LSAME( TRANSA, 'N' ) )THEN
276
IF (LSAME(TRANSA,'N')) THEN
212
278
* Form B := alpha*inv( A )*B.
216
IF( ALPHA.NE.ONE )THEN
218
B( I, J ) = ALPHA*B( I, J )
222
IF( B( K, J ).NE.ZERO )THEN
224
$ B( K, J ) = B( K, J )/A( K, K )
226
B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
233
IF( ALPHA.NE.ONE )THEN
235
B( I, J ) = ALPHA*B( I, J )
239
IF( B( K, J ).NE.ZERO )THEN
241
$ B( K, J ) = B( K, J )/A( K, K )
243
B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
251
* Form B := alpha*inv( A' )*B
252
* or B := alpha*inv( conjg( A' ) )*B.
257
TEMP = ALPHA*B( I, J )
260
TEMP = TEMP - A( K, I )*B( K, J )
263
$ TEMP = TEMP/A( I, I )
266
TEMP = TEMP - DCONJG( A( K, I ) )*B( K, J )
269
$ TEMP = TEMP/DCONJG( A( I, I ) )
277
TEMP = ALPHA*B( I, J )
280
TEMP = TEMP - A( K, I )*B( K, J )
283
$ TEMP = TEMP/A( I, I )
286
TEMP = TEMP - DCONJG( A( K, I ) )*B( K, J )
289
$ TEMP = TEMP/DCONJG( A( I, I ) )
282
IF (ALPHA.NE.ONE) THEN
284
B(I,J) = ALPHA*B(I,J)
288
IF (B(K,J).NE.ZERO) THEN
289
IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
291
B(I,J) = B(I,J) - B(K,J)*A(I,K)
298
IF (ALPHA.NE.ONE) THEN
300
B(I,J) = ALPHA*B(I,J)
304
IF (B(K,J).NE.ZERO) THEN
305
IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
307
B(I,J) = B(I,J) - B(K,J)*A(I,K)
315
* Form B := alpha*inv( A**T )*B
316
* or B := alpha*inv( A**H )*B.
324
TEMP = TEMP - A(K,I)*B(K,J)
326
IF (NOUNIT) TEMP = TEMP/A(I,I)
329
TEMP = TEMP - DCONJG(A(K,I))*B(K,J)
331
IF (NOUNIT) TEMP = TEMP/DCONJG(A(I,I))
342
TEMP = TEMP - A(K,I)*B(K,J)
344
IF (NOUNIT) TEMP = TEMP/A(I,I)
347
TEMP = TEMP - DCONJG(A(K,I))*B(K,J)
349
IF (NOUNIT) TEMP = TEMP/DCONJG(A(I,I))
297
IF( LSAME( TRANSA, 'N' ) )THEN
357
IF (LSAME(TRANSA,'N')) THEN
299
359
* Form B := alpha*B*inv( A ).
303
IF( ALPHA.NE.ONE )THEN
305
B( I, J ) = ALPHA*B( I, J )
309
IF( A( K, J ).NE.ZERO )THEN
311
B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
318
B( I, J ) = TEMP*B( I, J )
324
IF( ALPHA.NE.ONE )THEN
326
B( I, J ) = ALPHA*B( I, J )
330
IF( A( K, J ).NE.ZERO )THEN
332
B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
339
B( I, J ) = TEMP*B( I, J )
346
* Form B := alpha*B*inv( A' )
347
* or B := alpha*B*inv( conjg( A' ) ).
355
TEMP = ONE/DCONJG( A( K, K ) )
358
B( I, K ) = TEMP*B( I, K )
362
IF( A( J, K ).NE.ZERO )THEN
366
TEMP = DCONJG( A( J, K ) )
369
B( I, J ) = B( I, J ) - TEMP*B( I, K )
373
IF( ALPHA.NE.ONE )THEN
375
B( I, K ) = ALPHA*B( I, K )
385
TEMP = ONE/DCONJG( A( K, K ) )
388
B( I, K ) = TEMP*B( I, K )
392
IF( A( J, K ).NE.ZERO )THEN
396
TEMP = DCONJG( A( J, K ) )
399
B( I, J ) = B( I, J ) - TEMP*B( I, K )
403
IF( ALPHA.NE.ONE )THEN
405
B( I, K ) = ALPHA*B( I, K )
363
IF (ALPHA.NE.ONE) THEN
365
B(I,J) = ALPHA*B(I,J)
369
IF (A(K,J).NE.ZERO) THEN
371
B(I,J) = B(I,J) - A(K,J)*B(I,K)
384
IF (ALPHA.NE.ONE) THEN
386
B(I,J) = ALPHA*B(I,J)
390
IF (A(K,J).NE.ZERO) THEN
392
B(I,J) = B(I,J) - A(K,J)*B(I,K)
406
* Form B := alpha*B*inv( A**T )
407
* or B := alpha*B*inv( A**H ).
415
TEMP = ONE/DCONJG(A(K,K))
422
IF (A(J,K).NE.ZERO) THEN
426
TEMP = DCONJG(A(J,K))
429
B(I,J) = B(I,J) - TEMP*B(I,K)
433
IF (ALPHA.NE.ONE) THEN
435
B(I,K) = ALPHA*B(I,K)
445
TEMP = ONE/DCONJG(A(K,K))
452
IF (A(J,K).NE.ZERO) THEN
456
TEMP = DCONJG(A(J,K))
459
B(I,J) = B(I,J) - TEMP*B(I,K)
463
IF (ALPHA.NE.ONE) THEN
465
B(I,K) = ALPHA*B(I,K)