Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/blas/ctrsm.f
5218 views
1
SUBROUTINE CTRSM ( 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
COMPLEX ALPHA
7
* .. Array Arguments ..
8
COMPLEX A( LDA, * ), B( LDB, * )
9
* ..
10
*
11
* Purpose
12
* =======
13
*
14
* CTRSM 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' or op( A ) = conjg( 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 ) = conjg( 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 - COMPLEX .
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 - COMPLEX 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 - COMPLEX 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
* -- Written on 8-February-1989.
123
* Jack Dongarra, Argonne National Laboratory.
124
* Iain Duff, AERE Harwell.
125
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
126
* Sven Hammarling, Numerical Algorithms Group Ltd.
127
*
128
*
129
* .. External Functions ..
130
LOGICAL LSAME
131
EXTERNAL LSAME
132
* .. External Subroutines ..
133
EXTERNAL XERBLA
134
* .. Intrinsic Functions ..
135
INTRINSIC CONJG, MAX
136
* .. Local Scalars ..
137
LOGICAL LSIDE, NOCONJ, NOUNIT, UPPER
138
INTEGER I, INFO, J, K, NROWA
139
COMPLEX TEMP
140
* .. Parameters ..
141
COMPLEX ONE
142
PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
143
COMPLEX ZERO
144
PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ) )
145
* ..
146
* .. Executable Statements ..
147
*
148
* Test the input parameters.
149
*
150
LSIDE = LSAME( SIDE , 'L' )
151
IF( LSIDE )THEN
152
NROWA = M
153
ELSE
154
NROWA = N
155
END IF
156
NOCONJ = LSAME( TRANSA, 'T' )
157
NOUNIT = LSAME( DIAG , 'N' )
158
UPPER = LSAME( UPLO , 'U' )
159
*
160
INFO = 0
161
IF( ( .NOT.LSIDE ).AND.
162
$ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN
163
INFO = 1
164
ELSE IF( ( .NOT.UPPER ).AND.
165
$ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN
166
INFO = 2
167
ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
168
$ ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
169
$ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN
170
INFO = 3
171
ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND.
172
$ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN
173
INFO = 4
174
ELSE IF( M .LT.0 )THEN
175
INFO = 5
176
ELSE IF( N .LT.0 )THEN
177
INFO = 6
178
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
179
INFO = 9
180
ELSE IF( LDB.LT.MAX( 1, M ) )THEN
181
INFO = 11
182
END IF
183
IF( INFO.NE.0 )THEN
184
CALL XERBLA( 'CTRSM ', INFO )
185
RETURN
186
END IF
187
*
188
* Quick return if possible.
189
*
190
IF( N.EQ.0 )
191
$ RETURN
192
*
193
* And when alpha.eq.zero.
194
*
195
IF( ALPHA.EQ.ZERO )THEN
196
DO 20, J = 1, N
197
DO 10, I = 1, M
198
B( I, J ) = ZERO
199
10 CONTINUE
200
20 CONTINUE
201
RETURN
202
END IF
203
*
204
* Start the operations.
205
*
206
IF( LSIDE )THEN
207
IF( LSAME( TRANSA, 'N' ) )THEN
208
*
209
* Form B := alpha*inv( A )*B.
210
*
211
IF( UPPER )THEN
212
DO 60, J = 1, N
213
IF( ALPHA.NE.ONE )THEN
214
DO 30, I = 1, M
215
B( I, J ) = ALPHA*B( I, J )
216
30 CONTINUE
217
END IF
218
DO 50, K = M, 1, -1
219
IF( B( K, J ).NE.ZERO )THEN
220
IF( NOUNIT )
221
$ B( K, J ) = B( K, J )/A( K, K )
222
DO 40, I = 1, K - 1
223
B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
224
40 CONTINUE
225
END IF
226
50 CONTINUE
227
60 CONTINUE
228
ELSE
229
DO 100, J = 1, N
230
IF( ALPHA.NE.ONE )THEN
231
DO 70, I = 1, M
232
B( I, J ) = ALPHA*B( I, J )
233
70 CONTINUE
234
END IF
235
DO 90 K = 1, M
236
IF( B( K, J ).NE.ZERO )THEN
237
IF( NOUNIT )
238
$ B( K, J ) = B( K, J )/A( K, K )
239
DO 80, I = K + 1, M
240
B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
241
80 CONTINUE
242
END IF
243
90 CONTINUE
244
100 CONTINUE
245
END IF
246
ELSE
247
*
248
* Form B := alpha*inv( A' )*B
249
* or B := alpha*inv( conjg( A' ) )*B.
250
*
251
IF( UPPER )THEN
252
DO 140, J = 1, N
253
DO 130, I = 1, M
254
TEMP = ALPHA*B( I, J )
255
IF( NOCONJ )THEN
256
DO 110, K = 1, I - 1
257
TEMP = TEMP - A( K, I )*B( K, J )
258
110 CONTINUE
259
IF( NOUNIT )
260
$ TEMP = TEMP/A( I, I )
261
ELSE
262
DO 120, K = 1, I - 1
263
TEMP = TEMP - CONJG( A( K, I ) )*B( K, J )
264
120 CONTINUE
265
IF( NOUNIT )
266
$ TEMP = TEMP/CONJG( A( I, I ) )
267
END IF
268
B( I, J ) = TEMP
269
130 CONTINUE
270
140 CONTINUE
271
ELSE
272
DO 180, J = 1, N
273
DO 170, I = M, 1, -1
274
TEMP = ALPHA*B( I, J )
275
IF( NOCONJ )THEN
276
DO 150, K = I + 1, M
277
TEMP = TEMP - A( K, I )*B( K, J )
278
150 CONTINUE
279
IF( NOUNIT )
280
$ TEMP = TEMP/A( I, I )
281
ELSE
282
DO 160, K = I + 1, M
283
TEMP = TEMP - CONJG( A( K, I ) )*B( K, J )
284
160 CONTINUE
285
IF( NOUNIT )
286
$ TEMP = TEMP/CONJG( A( I, I ) )
287
END IF
288
B( I, J ) = TEMP
289
170 CONTINUE
290
180 CONTINUE
291
END IF
292
END IF
293
ELSE
294
IF( LSAME( TRANSA, 'N' ) )THEN
295
*
296
* Form B := alpha*B*inv( A ).
297
*
298
IF( UPPER )THEN
299
DO 230, J = 1, N
300
IF( ALPHA.NE.ONE )THEN
301
DO 190, I = 1, M
302
B( I, J ) = ALPHA*B( I, J )
303
190 CONTINUE
304
END IF
305
DO 210, K = 1, J - 1
306
IF( A( K, J ).NE.ZERO )THEN
307
DO 200, I = 1, M
308
B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
309
200 CONTINUE
310
END IF
311
210 CONTINUE
312
IF( NOUNIT )THEN
313
TEMP = ONE/A( J, J )
314
DO 220, I = 1, M
315
B( I, J ) = TEMP*B( I, J )
316
220 CONTINUE
317
END IF
318
230 CONTINUE
319
ELSE
320
DO 280, J = N, 1, -1
321
IF( ALPHA.NE.ONE )THEN
322
DO 240, I = 1, M
323
B( I, J ) = ALPHA*B( I, J )
324
240 CONTINUE
325
END IF
326
DO 260, K = J + 1, N
327
IF( A( K, J ).NE.ZERO )THEN
328
DO 250, I = 1, M
329
B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
330
250 CONTINUE
331
END IF
332
260 CONTINUE
333
IF( NOUNIT )THEN
334
TEMP = ONE/A( J, J )
335
DO 270, I = 1, M
336
B( I, J ) = TEMP*B( I, J )
337
270 CONTINUE
338
END IF
339
280 CONTINUE
340
END IF
341
ELSE
342
*
343
* Form B := alpha*B*inv( A' )
344
* or B := alpha*B*inv( conjg( A' ) ).
345
*
346
IF( UPPER )THEN
347
DO 330, K = N, 1, -1
348
IF( NOUNIT )THEN
349
IF( NOCONJ )THEN
350
TEMP = ONE/A( K, K )
351
ELSE
352
TEMP = ONE/CONJG( A( K, K ) )
353
END IF
354
DO 290, I = 1, M
355
B( I, K ) = TEMP*B( I, K )
356
290 CONTINUE
357
END IF
358
DO 310, J = 1, K - 1
359
IF( A( J, K ).NE.ZERO )THEN
360
IF( NOCONJ )THEN
361
TEMP = A( J, K )
362
ELSE
363
TEMP = CONJG( A( J, K ) )
364
END IF
365
DO 300, I = 1, M
366
B( I, J ) = B( I, J ) - TEMP*B( I, K )
367
300 CONTINUE
368
END IF
369
310 CONTINUE
370
IF( ALPHA.NE.ONE )THEN
371
DO 320, I = 1, M
372
B( I, K ) = ALPHA*B( I, K )
373
320 CONTINUE
374
END IF
375
330 CONTINUE
376
ELSE
377
DO 380, K = 1, N
378
IF( NOUNIT )THEN
379
IF( NOCONJ )THEN
380
TEMP = ONE/A( K, K )
381
ELSE
382
TEMP = ONE/CONJG( A( K, K ) )
383
END IF
384
DO 340, I = 1, M
385
B( I, K ) = TEMP*B( I, K )
386
340 CONTINUE
387
END IF
388
DO 360, J = K + 1, N
389
IF( A( J, K ).NE.ZERO )THEN
390
IF( NOCONJ )THEN
391
TEMP = A( J, K )
392
ELSE
393
TEMP = CONJG( A( J, K ) )
394
END IF
395
DO 350, I = 1, M
396
B( I, J ) = B( I, J ) - TEMP*B( I, K )
397
350 CONTINUE
398
END IF
399
360 CONTINUE
400
IF( ALPHA.NE.ONE )THEN
401
DO 370, I = 1, M
402
B( I, K ) = ALPHA*B( I, K )
403
370 CONTINUE
404
END IF
405
380 CONTINUE
406
END IF
407
END IF
408
END IF
409
*
410
RETURN
411
*
412
* End of CTRSM .
413
*
414
END
415
416