|
2329
|
1 SUBROUTINE ZTRSM ( 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*16 ALPHA |
|
|
7 * .. Array Arguments .. |
|
|
8 COMPLEX*16 A( LDA, * ), B( LDB, * ) |
|
|
9 * .. |
|
|
10 * |
|
|
11 * Purpose |
|
|
12 * ======= |
|
|
13 * |
|
|
14 * ZTRSM 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*16 . |
|
|
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*16 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*16 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 DCONJG, MAX |
|
|
136 * .. Local Scalars .. |
|
|
137 LOGICAL LSIDE, NOCONJ, NOUNIT, UPPER |
|
|
138 INTEGER I, INFO, J, K, NROWA |
|
|
139 COMPLEX*16 TEMP |
|
|
140 * .. Parameters .. |
|
|
141 COMPLEX*16 ONE |
|
|
142 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) |
|
|
143 COMPLEX*16 ZERO |
|
|
144 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+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( 'ZTRSM ', 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 - DCONJG( A( K, I ) )*B( K, J ) |
|
|
264 120 CONTINUE |
|
|
265 IF( NOUNIT ) |
|
|
266 $ TEMP = TEMP/DCONJG( 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 - DCONJG( A( K, I ) )*B( K, J ) |
|
|
284 160 CONTINUE |
|
|
285 IF( NOUNIT ) |
|
|
286 $ TEMP = TEMP/DCONJG( 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/DCONJG( 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 = DCONJG( 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/DCONJG( 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 = DCONJG( 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 ZTRSM . |
|
|
413 * |
|
|
414 END |
