Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/blas/ctrmm.f
5182 views
1
SUBROUTINE CTRMM ( 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
* CTRMM performs one of the matrix-matrix operations
15
*
16
* B := alpha*op( A )*B, or B := alpha*B*op( A )
17
*
18
* where alpha is a scalar, B is an m by n matrix, 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
* Parameters
24
* ==========
25
*
26
* SIDE - CHARACTER*1.
27
* On entry, SIDE specifies whether op( A ) multiplies B from
28
* the left or right as follows:
29
*
30
* SIDE = 'L' or 'l' B := alpha*op( A )*B.
31
*
32
* SIDE = 'R' or 'r' B := alpha*B*op( A ).
33
*
34
* Unchanged on exit.
35
*
36
* UPLO - CHARACTER*1.
37
* On entry, UPLO specifies whether the matrix A is an upper or
38
* lower triangular matrix as follows:
39
*
40
* UPLO = 'U' or 'u' A is an upper triangular matrix.
41
*
42
* UPLO = 'L' or 'l' A is a lower triangular matrix.
43
*
44
* Unchanged on exit.
45
*
46
* TRANSA - CHARACTER*1.
47
* On entry, TRANSA specifies the form of op( A ) to be used in
48
* the matrix multiplication as follows:
49
*
50
* TRANSA = 'N' or 'n' op( A ) = A.
51
*
52
* TRANSA = 'T' or 't' op( A ) = A'.
53
*
54
* TRANSA = 'C' or 'c' op( A ) = conjg( A' ).
55
*
56
* Unchanged on exit.
57
*
58
* DIAG - CHARACTER*1.
59
* On entry, DIAG specifies whether or not A is unit triangular
60
* as follows:
61
*
62
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
63
*
64
* DIAG = 'N' or 'n' A is not assumed to be unit
65
* triangular.
66
*
67
* Unchanged on exit.
68
*
69
* M - INTEGER.
70
* On entry, M specifies the number of rows of B. M must be at
71
* least zero.
72
* Unchanged on exit.
73
*
74
* N - INTEGER.
75
* On entry, N specifies the number of columns of B. N must be
76
* at least zero.
77
* Unchanged on exit.
78
*
79
* ALPHA - COMPLEX .
80
* On entry, ALPHA specifies the scalar alpha. When alpha is
81
* zero then A is not referenced and B need not be set before
82
* entry.
83
* Unchanged on exit.
84
*
85
* A - COMPLEX array of DIMENSION ( LDA, k ), where k is m
86
* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
87
* Before entry with UPLO = 'U' or 'u', the leading k by k
88
* upper triangular part of the array A must contain the upper
89
* triangular matrix and the strictly lower triangular part of
90
* A is not referenced.
91
* Before entry with UPLO = 'L' or 'l', the leading k by k
92
* lower triangular part of the array A must contain the lower
93
* triangular matrix and the strictly upper triangular part of
94
* A is not referenced.
95
* Note that when DIAG = 'U' or 'u', the diagonal elements of
96
* A are not referenced either, but are assumed to be unity.
97
* Unchanged on exit.
98
*
99
* LDA - INTEGER.
100
* On entry, LDA specifies the first dimension of A as declared
101
* in the calling (sub) program. When SIDE = 'L' or 'l' then
102
* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
103
* then LDA must be at least max( 1, n ).
104
* Unchanged on exit.
105
*
106
* B - COMPLEX array of DIMENSION ( LDB, n ).
107
* Before entry, the leading m by n part of the array B must
108
* contain the matrix B, and on exit is overwritten by the
109
* transformed matrix.
110
*
111
* LDB - INTEGER.
112
* On entry, LDB specifies the first dimension of B as declared
113
* in the calling (sub) program. LDB must be at least
114
* max( 1, m ).
115
* Unchanged on exit.
116
*
117
*
118
* Level 3 Blas routine.
119
*
120
* -- Written on 8-February-1989.
121
* Jack Dongarra, Argonne National Laboratory.
122
* Iain Duff, AERE Harwell.
123
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
124
* Sven Hammarling, Numerical Algorithms Group Ltd.
125
*
126
*
127
* .. External Functions ..
128
LOGICAL LSAME
129
EXTERNAL LSAME
130
* .. External Subroutines ..
131
EXTERNAL XERBLA
132
* .. Intrinsic Functions ..
133
INTRINSIC CONJG, MAX
134
* .. Local Scalars ..
135
LOGICAL LSIDE, NOCONJ, NOUNIT, UPPER
136
INTEGER I, INFO, J, K, NROWA
137
COMPLEX TEMP
138
* .. Parameters ..
139
COMPLEX ONE
140
PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
141
COMPLEX ZERO
142
PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ) )
143
* ..
144
* .. Executable Statements ..
145
*
146
* Test the input parameters.
147
*
148
LSIDE = LSAME( SIDE , 'L' )
149
IF( LSIDE )THEN
150
NROWA = M
151
ELSE
152
NROWA = N
153
END IF
154
NOCONJ = LSAME( TRANSA, 'T' )
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( 'CTRMM ', 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*A*B.
208
*
209
IF( UPPER )THEN
210
DO 50, J = 1, N
211
DO 40, K = 1, M
212
IF( B( K, J ).NE.ZERO )THEN
213
TEMP = ALPHA*B( K, J )
214
DO 30, I = 1, K - 1
215
B( I, J ) = B( I, J ) + TEMP*A( I, K )
216
30 CONTINUE
217
IF( NOUNIT )
218
$ TEMP = TEMP*A( K, K )
219
B( K, J ) = TEMP
220
END IF
221
40 CONTINUE
222
50 CONTINUE
223
ELSE
224
DO 80, J = 1, N
225
DO 70 K = M, 1, -1
226
IF( B( K, J ).NE.ZERO )THEN
227
TEMP = ALPHA*B( K, J )
228
B( K, J ) = TEMP
229
IF( NOUNIT )
230
$ B( K, J ) = B( K, J )*A( K, K )
231
DO 60, I = K + 1, M
232
B( I, J ) = B( I, J ) + TEMP*A( I, K )
233
60 CONTINUE
234
END IF
235
70 CONTINUE
236
80 CONTINUE
237
END IF
238
ELSE
239
*
240
* Form B := alpha*A'*B or B := alpha*conjg( A' )*B.
241
*
242
IF( UPPER )THEN
243
DO 120, J = 1, N
244
DO 110, I = M, 1, -1
245
TEMP = B( I, J )
246
IF( NOCONJ )THEN
247
IF( NOUNIT )
248
$ TEMP = TEMP*A( I, I )
249
DO 90, K = 1, I - 1
250
TEMP = TEMP + A( K, I )*B( K, J )
251
90 CONTINUE
252
ELSE
253
IF( NOUNIT )
254
$ TEMP = TEMP*CONJG( A( I, I ) )
255
DO 100, K = 1, I - 1
256
TEMP = TEMP + CONJG( A( K, I ) )*B( K, J )
257
100 CONTINUE
258
END IF
259
B( I, J ) = ALPHA*TEMP
260
110 CONTINUE
261
120 CONTINUE
262
ELSE
263
DO 160, J = 1, N
264
DO 150, I = 1, M
265
TEMP = B( I, J )
266
IF( NOCONJ )THEN
267
IF( NOUNIT )
268
$ TEMP = TEMP*A( I, I )
269
DO 130, K = I + 1, M
270
TEMP = TEMP + A( K, I )*B( K, J )
271
130 CONTINUE
272
ELSE
273
IF( NOUNIT )
274
$ TEMP = TEMP*CONJG( A( I, I ) )
275
DO 140, K = I + 1, M
276
TEMP = TEMP + CONJG( A( K, I ) )*B( K, J )
277
140 CONTINUE
278
END IF
279
B( I, J ) = ALPHA*TEMP
280
150 CONTINUE
281
160 CONTINUE
282
END IF
283
END IF
284
ELSE
285
IF( LSAME( TRANSA, 'N' ) )THEN
286
*
287
* Form B := alpha*B*A.
288
*
289
IF( UPPER )THEN
290
DO 200, J = N, 1, -1
291
TEMP = ALPHA
292
IF( NOUNIT )
293
$ TEMP = TEMP*A( J, J )
294
DO 170, I = 1, M
295
B( I, J ) = TEMP*B( I, J )
296
170 CONTINUE
297
DO 190, K = 1, J - 1
298
IF( A( K, J ).NE.ZERO )THEN
299
TEMP = ALPHA*A( K, J )
300
DO 180, I = 1, M
301
B( I, J ) = B( I, J ) + TEMP*B( I, K )
302
180 CONTINUE
303
END IF
304
190 CONTINUE
305
200 CONTINUE
306
ELSE
307
DO 240, J = 1, N
308
TEMP = ALPHA
309
IF( NOUNIT )
310
$ TEMP = TEMP*A( J, J )
311
DO 210, I = 1, M
312
B( I, J ) = TEMP*B( I, J )
313
210 CONTINUE
314
DO 230, K = J + 1, N
315
IF( A( K, J ).NE.ZERO )THEN
316
TEMP = ALPHA*A( K, J )
317
DO 220, I = 1, M
318
B( I, J ) = B( I, J ) + TEMP*B( I, K )
319
220 CONTINUE
320
END IF
321
230 CONTINUE
322
240 CONTINUE
323
END IF
324
ELSE
325
*
326
* Form B := alpha*B*A' or B := alpha*B*conjg( A' ).
327
*
328
IF( UPPER )THEN
329
DO 280, K = 1, N
330
DO 260, J = 1, K - 1
331
IF( A( J, K ).NE.ZERO )THEN
332
IF( NOCONJ )THEN
333
TEMP = ALPHA*A( J, K )
334
ELSE
335
TEMP = ALPHA*CONJG( A( J, K ) )
336
END IF
337
DO 250, I = 1, M
338
B( I, J ) = B( I, J ) + TEMP*B( I, K )
339
250 CONTINUE
340
END IF
341
260 CONTINUE
342
TEMP = ALPHA
343
IF( NOUNIT )THEN
344
IF( NOCONJ )THEN
345
TEMP = TEMP*A( K, K )
346
ELSE
347
TEMP = TEMP*CONJG( A( K, K ) )
348
END IF
349
END IF
350
IF( TEMP.NE.ONE )THEN
351
DO 270, I = 1, M
352
B( I, K ) = TEMP*B( I, K )
353
270 CONTINUE
354
END IF
355
280 CONTINUE
356
ELSE
357
DO 320, K = N, 1, -1
358
DO 300, J = K + 1, N
359
IF( A( J, K ).NE.ZERO )THEN
360
IF( NOCONJ )THEN
361
TEMP = ALPHA*A( J, K )
362
ELSE
363
TEMP = ALPHA*CONJG( A( J, K ) )
364
END IF
365
DO 290, I = 1, M
366
B( I, J ) = B( I, J ) + TEMP*B( I, K )
367
290 CONTINUE
368
END IF
369
300 CONTINUE
370
TEMP = ALPHA
371
IF( NOUNIT )THEN
372
IF( NOCONJ )THEN
373
TEMP = TEMP*A( K, K )
374
ELSE
375
TEMP = TEMP*CONJG( A( K, K ) )
376
END IF
377
END IF
378
IF( TEMP.NE.ONE )THEN
379
DO 310, I = 1, M
380
B( I, K ) = TEMP*B( I, K )
381
310 CONTINUE
382
END IF
383
320 CONTINUE
384
END IF
385
END IF
386
END IF
387
*
388
RETURN
389
*
390
* End of CTRMM .
391
*
392
END
393
394