Annotation of rpl/lapack/lapack/dtfsm.f, revision 1.6
1.1 bertrand 1: SUBROUTINE DTFSM( TRANSR, SIDE, UPLO, TRANS, DIAG, M, N, ALPHA, A,
1.6 ! bertrand 2: $ B, LDB )
1.1 bertrand 3: *
1.6 ! bertrand 4: * -- LAPACK routine (version 3.3.1) --
1.1 bertrand 5: *
6: * -- Contributed by Fred Gustavson of the IBM Watson Research Center --
1.6 ! bertrand 7: * -- April 2011 --
1.1 bertrand 8: *
9: * -- LAPACK is a software package provided by Univ. of Tennessee, --
10: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
11: *
12: * ..
13: * .. Scalar Arguments ..
14: CHARACTER TRANSR, DIAG, SIDE, TRANS, UPLO
15: INTEGER LDB, M, N
16: DOUBLE PRECISION ALPHA
17: * ..
18: * .. Array Arguments ..
19: DOUBLE PRECISION A( 0: * ), B( 0: LDB-1, 0: * )
20: * ..
21: *
22: * Purpose
23: * =======
24: *
25: * Level 3 BLAS like routine for A in RFP Format.
26: *
27: * DTFSM solves the matrix equation
28: *
29: * op( A )*X = alpha*B or X*op( A ) = alpha*B
30: *
31: * where alpha is a scalar, X and B are m by n matrices, A is a unit, or
32: * non-unit, upper or lower triangular matrix and op( A ) is one of
33: *
1.6 ! bertrand 34: * op( A ) = A or op( A ) = A**T.
1.1 bertrand 35: *
36: * A is in Rectangular Full Packed (RFP) Format.
37: *
38: * The matrix X is overwritten on B.
39: *
40: * Arguments
41: * ==========
42: *
1.4 bertrand 43: * TRANSR (input) CHARACTER*1
1.1 bertrand 44: * = 'N': The Normal Form of RFP A is stored;
45: * = 'T': The Transpose Form of RFP A is stored.
46: *
1.4 bertrand 47: * SIDE (input) CHARACTER*1
1.1 bertrand 48: * On entry, SIDE specifies whether op( A ) appears on the left
49: * or right of X as follows:
50: *
51: * SIDE = 'L' or 'l' op( A )*X = alpha*B.
52: *
53: * SIDE = 'R' or 'r' X*op( A ) = alpha*B.
54: *
55: * Unchanged on exit.
56: *
1.4 bertrand 57: * UPLO (input) CHARACTER*1
1.1 bertrand 58: * On entry, UPLO specifies whether the RFP matrix A came from
59: * an upper or lower triangular matrix as follows:
60: * UPLO = 'U' or 'u' RFP A came from an upper triangular matrix
61: * UPLO = 'L' or 'l' RFP A came from a lower triangular matrix
62: *
63: * Unchanged on exit.
64: *
1.4 bertrand 65: * TRANS (input) CHARACTER*1
1.1 bertrand 66: * On entry, TRANS specifies the form of op( A ) to be used
67: * in the matrix multiplication as follows:
68: *
69: * TRANS = 'N' or 'n' op( A ) = A.
70: *
71: * TRANS = 'T' or 't' op( A ) = A'.
72: *
73: * Unchanged on exit.
74: *
1.4 bertrand 75: * DIAG (input) CHARACTER*1
1.1 bertrand 76: * On entry, DIAG specifies whether or not RFP A is unit
77: * triangular as follows:
78: *
79: * DIAG = 'U' or 'u' A is assumed to be unit triangular.
80: *
81: * DIAG = 'N' or 'n' A is not assumed to be unit
82: * triangular.
83: *
84: * Unchanged on exit.
85: *
86: * M (input) INTEGER
87: * On entry, M specifies the number of rows of B. M must be at
88: * least zero.
89: * Unchanged on exit.
90: *
91: * N (input) INTEGER
92: * On entry, N specifies the number of columns of B. N must be
93: * at least zero.
94: * Unchanged on exit.
95: *
96: * ALPHA (input) DOUBLE PRECISION
97: * On entry, ALPHA specifies the scalar alpha. When alpha is
98: * zero then A is not referenced and B need not be set before
99: * entry.
100: * Unchanged on exit.
101: *
102: * A (input) DOUBLE PRECISION array, dimension (NT)
103: * NT = N*(N+1)/2. On entry, the matrix A in RFP Format.
104: * RFP Format is described by TRANSR, UPLO and N as follows:
105: * If TRANSR='N' then RFP A is (0:N,0:K-1) when N is even;
106: * K=N/2. RFP A is (0:N-1,0:K) when N is odd; K=N/2. If
107: * TRANSR = 'T' then RFP is the transpose of RFP A as
108: * defined when TRANSR = 'N'. The contents of RFP A are defined
109: * by UPLO as follows: If UPLO = 'U' the RFP A contains the NT
110: * elements of upper packed A either in normal or
111: * transpose Format. If UPLO = 'L' the RFP A contains
112: * the NT elements of lower packed A either in normal or
113: * transpose Format. The LDA of RFP A is (N+1)/2 when
114: * TRANSR = 'T'. When TRANSR is 'N' the LDA is N+1 when N is
115: * even and is N when is odd.
116: * See the Note below for more details. Unchanged on exit.
117: *
118: * B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
119: * Before entry, the leading m by n part of the array B must
120: * contain the right-hand side matrix B, and on exit is
121: * overwritten by the solution matrix X.
122: *
123: * LDB (input) INTEGER
124: * On entry, LDB specifies the first dimension of B as declared
125: * in the calling (sub) program. LDB must be at least
126: * max( 1, m ).
127: * Unchanged on exit.
128: *
129: * Further Details
130: * ===============
131: *
132: * We first consider Rectangular Full Packed (RFP) Format when N is
133: * even. We give an example where N = 6.
134: *
135: * AP is Upper AP is Lower
136: *
137: * 00 01 02 03 04 05 00
138: * 11 12 13 14 15 10 11
139: * 22 23 24 25 20 21 22
140: * 33 34 35 30 31 32 33
141: * 44 45 40 41 42 43 44
142: * 55 50 51 52 53 54 55
143: *
144: *
145: * Let TRANSR = 'N'. RFP holds AP as follows:
146: * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
147: * three columns of AP upper. The lower triangle A(4:6,0:2) consists of
148: * the transpose of the first three columns of AP upper.
149: * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
150: * three columns of AP lower. The upper triangle A(0:2,0:2) consists of
151: * the transpose of the last three columns of AP lower.
152: * This covers the case N even and TRANSR = 'N'.
153: *
154: * RFP A RFP A
155: *
156: * 03 04 05 33 43 53
157: * 13 14 15 00 44 54
158: * 23 24 25 10 11 55
159: * 33 34 35 20 21 22
160: * 00 44 45 30 31 32
161: * 01 11 55 40 41 42
162: * 02 12 22 50 51 52
163: *
164: * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
165: * transpose of RFP A above. One therefore gets:
166: *
167: *
168: * RFP A RFP A
169: *
170: * 03 13 23 33 00 01 02 33 00 10 20 30 40 50
171: * 04 14 24 34 44 11 12 43 44 11 21 31 41 51
172: * 05 15 25 35 45 55 22 53 54 55 22 32 42 52
173: *
174: *
175: * We then consider Rectangular Full Packed (RFP) Format when N is
176: * odd. We give an example where N = 5.
177: *
178: * AP is Upper AP is Lower
179: *
180: * 00 01 02 03 04 00
181: * 11 12 13 14 10 11
182: * 22 23 24 20 21 22
183: * 33 34 30 31 32 33
184: * 44 40 41 42 43 44
185: *
186: *
187: * Let TRANSR = 'N'. RFP holds AP as follows:
188: * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
189: * three columns of AP upper. The lower triangle A(3:4,0:1) consists of
190: * the transpose of the first two columns of AP upper.
191: * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
192: * three columns of AP lower. The upper triangle A(0:1,1:2) consists of
193: * the transpose of the last two columns of AP lower.
194: * This covers the case N odd and TRANSR = 'N'.
195: *
196: * RFP A RFP A
197: *
198: * 02 03 04 00 33 43
199: * 12 13 14 10 11 44
200: * 22 23 24 20 21 22
201: * 00 33 34 30 31 32
202: * 01 11 44 40 41 42
203: *
204: * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
205: * transpose of RFP A above. One therefore gets:
206: *
207: * RFP A RFP A
208: *
209: * 02 12 22 00 01 00 10 20 30 40 50
210: * 03 13 23 33 11 33 11 21 31 41 51
211: * 04 14 24 34 44 43 44 22 32 42 52
212: *
213: * Reference
214: * =========
215: *
216: * =====================================================================
217: *
218: * ..
219: * .. Parameters ..
220: DOUBLE PRECISION ONE, ZERO
221: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
222: * ..
223: * .. Local Scalars ..
224: LOGICAL LOWER, LSIDE, MISODD, NISODD, NORMALTRANSR,
1.6 ! bertrand 225: $ NOTRANS
1.1 bertrand 226: INTEGER M1, M2, N1, N2, K, INFO, I, J
227: * ..
228: * .. External Functions ..
229: LOGICAL LSAME
230: EXTERNAL LSAME
231: * ..
232: * .. External Subroutines ..
233: EXTERNAL XERBLA, DGEMM, DTRSM
234: * ..
235: * .. Intrinsic Functions ..
236: INTRINSIC MAX, MOD
237: * ..
238: * .. Executable Statements ..
239: *
240: * Test the input parameters.
241: *
242: INFO = 0
243: NORMALTRANSR = LSAME( TRANSR, 'N' )
244: LSIDE = LSAME( SIDE, 'L' )
245: LOWER = LSAME( UPLO, 'L' )
246: NOTRANS = LSAME( TRANS, 'N' )
247: IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
248: INFO = -1
249: ELSE IF( .NOT.LSIDE .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
250: INFO = -2
251: ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
252: INFO = -3
253: ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
254: INFO = -4
255: ELSE IF( .NOT.LSAME( DIAG, 'N' ) .AND. .NOT.LSAME( DIAG, 'U' ) )
1.6 ! bertrand 256: $ THEN
1.1 bertrand 257: INFO = -5
258: ELSE IF( M.LT.0 ) THEN
259: INFO = -6
260: ELSE IF( N.LT.0 ) THEN
261: INFO = -7
262: ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
263: INFO = -11
264: END IF
265: IF( INFO.NE.0 ) THEN
266: CALL XERBLA( 'DTFSM ', -INFO )
267: RETURN
268: END IF
269: *
270: * Quick return when ( (N.EQ.0).OR.(M.EQ.0) )
271: *
272: IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )
1.6 ! bertrand 273: $ RETURN
1.1 bertrand 274: *
275: * Quick return when ALPHA.EQ.(0D+0)
276: *
277: IF( ALPHA.EQ.ZERO ) THEN
278: DO 20 J = 0, N - 1
279: DO 10 I = 0, M - 1
280: B( I, J ) = ZERO
281: 10 CONTINUE
282: 20 CONTINUE
283: RETURN
284: END IF
285: *
286: IF( LSIDE ) THEN
287: *
288: * SIDE = 'L'
289: *
290: * A is M-by-M.
291: * If M is odd, set NISODD = .TRUE., and M1 and M2.
292: * If M is even, NISODD = .FALSE., and M.
293: *
294: IF( MOD( M, 2 ).EQ.0 ) THEN
295: MISODD = .FALSE.
296: K = M / 2
297: ELSE
298: MISODD = .TRUE.
299: IF( LOWER ) THEN
300: M2 = M / 2
301: M1 = M - M2
302: ELSE
303: M1 = M / 2
304: M2 = M - M1
305: END IF
306: END IF
307: *
308: *
309: IF( MISODD ) THEN
310: *
311: * SIDE = 'L' and N is odd
312: *
313: IF( NORMALTRANSR ) THEN
314: *
315: * SIDE = 'L', N is odd, and TRANSR = 'N'
316: *
317: IF( LOWER ) THEN
318: *
319: * SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'L'
320: *
321: IF( NOTRANS ) THEN
322: *
323: * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
324: * TRANS = 'N'
325: *
326: IF( M.EQ.1 ) THEN
327: CALL DTRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA,
1.6 ! bertrand 328: $ A, M, B, LDB )
1.1 bertrand 329: ELSE
330: CALL DTRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA,
1.6 ! bertrand 331: $ A( 0 ), M, B, LDB )
1.1 bertrand 332: CALL DGEMM( 'N', 'N', M2, N, M1, -ONE, A( M1 ),
1.6 ! bertrand 333: $ M, B, LDB, ALPHA, B( M1, 0 ), LDB )
1.1 bertrand 334: CALL DTRSM( 'L', 'U', 'T', DIAG, M2, N, ONE,
1.6 ! bertrand 335: $ A( M ), M, B( M1, 0 ), LDB )
1.1 bertrand 336: END IF
337: *
338: ELSE
339: *
340: * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
341: * TRANS = 'T'
342: *
343: IF( M.EQ.1 ) THEN
344: CALL DTRSM( 'L', 'L', 'T', DIAG, M1, N, ALPHA,
1.6 ! bertrand 345: $ A( 0 ), M, B, LDB )
1.1 bertrand 346: ELSE
347: CALL DTRSM( 'L', 'U', 'N', DIAG, M2, N, ALPHA,
1.6 ! bertrand 348: $ A( M ), M, B( M1, 0 ), LDB )
1.1 bertrand 349: CALL DGEMM( 'T', 'N', M1, N, M2, -ONE, A( M1 ),
1.6 ! bertrand 350: $ M, B( M1, 0 ), LDB, ALPHA, B, LDB )
1.1 bertrand 351: CALL DTRSM( 'L', 'L', 'T', DIAG, M1, N, ONE,
1.6 ! bertrand 352: $ A( 0 ), M, B, LDB )
1.1 bertrand 353: END IF
354: *
355: END IF
356: *
357: ELSE
358: *
359: * SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'U'
360: *
361: IF( .NOT.NOTRANS ) THEN
362: *
363: * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
364: * TRANS = 'N'
365: *
366: CALL DTRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA,
1.6 ! bertrand 367: $ A( M2 ), M, B, LDB )
1.1 bertrand 368: CALL DGEMM( 'T', 'N', M2, N, M1, -ONE, A( 0 ), M,
1.6 ! bertrand 369: $ B, LDB, ALPHA, B( M1, 0 ), LDB )
1.1 bertrand 370: CALL DTRSM( 'L', 'U', 'T', DIAG, M2, N, ONE,
1.6 ! bertrand 371: $ A( M1 ), M, B( M1, 0 ), LDB )
1.1 bertrand 372: *
373: ELSE
374: *
375: * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
376: * TRANS = 'T'
377: *
378: CALL DTRSM( 'L', 'U', 'N', DIAG, M2, N, ALPHA,
1.6 ! bertrand 379: $ A( M1 ), M, B( M1, 0 ), LDB )
1.1 bertrand 380: CALL DGEMM( 'N', 'N', M1, N, M2, -ONE, A( 0 ), M,
1.6 ! bertrand 381: $ B( M1, 0 ), LDB, ALPHA, B, LDB )
1.1 bertrand 382: CALL DTRSM( 'L', 'L', 'T', DIAG, M1, N, ONE,
1.6 ! bertrand 383: $ A( M2 ), M, B, LDB )
1.1 bertrand 384: *
385: END IF
386: *
387: END IF
388: *
389: ELSE
390: *
391: * SIDE = 'L', N is odd, and TRANSR = 'T'
392: *
393: IF( LOWER ) THEN
394: *
395: * SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'L'
396: *
397: IF( NOTRANS ) THEN
398: *
399: * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and
400: * TRANS = 'N'
401: *
402: IF( M.EQ.1 ) THEN
403: CALL DTRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA,
1.6 ! bertrand 404: $ A( 0 ), M1, B, LDB )
1.1 bertrand 405: ELSE
406: CALL DTRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA,
1.6 ! bertrand 407: $ A( 0 ), M1, B, LDB )
1.1 bertrand 408: CALL DGEMM( 'T', 'N', M2, N, M1, -ONE,
1.6 ! bertrand 409: $ A( M1*M1 ), M1, B, LDB, ALPHA,
! 410: $ B( M1, 0 ), LDB )
1.1 bertrand 411: CALL DTRSM( 'L', 'L', 'N', DIAG, M2, N, ONE,
1.6 ! bertrand 412: $ A( 1 ), M1, B( M1, 0 ), LDB )
1.1 bertrand 413: END IF
414: *
415: ELSE
416: *
417: * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and
418: * TRANS = 'T'
419: *
420: IF( M.EQ.1 ) THEN
421: CALL DTRSM( 'L', 'U', 'N', DIAG, M1, N, ALPHA,
1.6 ! bertrand 422: $ A( 0 ), M1, B, LDB )
1.1 bertrand 423: ELSE
424: CALL DTRSM( 'L', 'L', 'T', DIAG, M2, N, ALPHA,
1.6 ! bertrand 425: $ A( 1 ), M1, B( M1, 0 ), LDB )
1.1 bertrand 426: CALL DGEMM( 'N', 'N', M1, N, M2, -ONE,
1.6 ! bertrand 427: $ A( M1*M1 ), M1, B( M1, 0 ), LDB,
! 428: $ ALPHA, B, LDB )
1.1 bertrand 429: CALL DTRSM( 'L', 'U', 'N', DIAG, M1, N, ONE,
1.6 ! bertrand 430: $ A( 0 ), M1, B, LDB )
1.1 bertrand 431: END IF
432: *
433: END IF
434: *
435: ELSE
436: *
437: * SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'U'
438: *
439: IF( .NOT.NOTRANS ) THEN
440: *
441: * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and
442: * TRANS = 'N'
443: *
444: CALL DTRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA,
1.6 ! bertrand 445: $ A( M2*M2 ), M2, B, LDB )
1.1 bertrand 446: CALL DGEMM( 'N', 'N', M2, N, M1, -ONE, A( 0 ), M2,
1.6 ! bertrand 447: $ B, LDB, ALPHA, B( M1, 0 ), LDB )
1.1 bertrand 448: CALL DTRSM( 'L', 'L', 'N', DIAG, M2, N, ONE,
1.6 ! bertrand 449: $ A( M1*M2 ), M2, B( M1, 0 ), LDB )
1.1 bertrand 450: *
451: ELSE
452: *
453: * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and
454: * TRANS = 'T'
455: *
456: CALL DTRSM( 'L', 'L', 'T', DIAG, M2, N, ALPHA,
1.6 ! bertrand 457: $ A( M1*M2 ), M2, B( M1, 0 ), LDB )
1.1 bertrand 458: CALL DGEMM( 'T', 'N', M1, N, M2, -ONE, A( 0 ), M2,
1.6 ! bertrand 459: $ B( M1, 0 ), LDB, ALPHA, B, LDB )
1.1 bertrand 460: CALL DTRSM( 'L', 'U', 'N', DIAG, M1, N, ONE,
1.6 ! bertrand 461: $ A( M2*M2 ), M2, B, LDB )
1.1 bertrand 462: *
463: END IF
464: *
465: END IF
466: *
467: END IF
468: *
469: ELSE
470: *
471: * SIDE = 'L' and N is even
472: *
473: IF( NORMALTRANSR ) THEN
474: *
475: * SIDE = 'L', N is even, and TRANSR = 'N'
476: *
477: IF( LOWER ) THEN
478: *
479: * SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'L'
480: *
481: IF( NOTRANS ) THEN
482: *
483: * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L',
484: * and TRANS = 'N'
485: *
486: CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ALPHA,
1.6 ! bertrand 487: $ A( 1 ), M+1, B, LDB )
1.1 bertrand 488: CALL DGEMM( 'N', 'N', K, N, K, -ONE, A( K+1 ),
1.6 ! bertrand 489: $ M+1, B, LDB, ALPHA, B( K, 0 ), LDB )
1.1 bertrand 490: CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ONE,
1.6 ! bertrand 491: $ A( 0 ), M+1, B( K, 0 ), LDB )
1.1 bertrand 492: *
493: ELSE
494: *
495: * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L',
496: * and TRANS = 'T'
497: *
498: CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ALPHA,
1.6 ! bertrand 499: $ A( 0 ), M+1, B( K, 0 ), LDB )
1.1 bertrand 500: CALL DGEMM( 'T', 'N', K, N, K, -ONE, A( K+1 ),
1.6 ! bertrand 501: $ M+1, B( K, 0 ), LDB, ALPHA, B, LDB )
1.1 bertrand 502: CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ONE,
1.6 ! bertrand 503: $ A( 1 ), M+1, B, LDB )
1.1 bertrand 504: *
505: END IF
506: *
507: ELSE
508: *
509: * SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'U'
510: *
511: IF( .NOT.NOTRANS ) THEN
512: *
513: * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U',
514: * and TRANS = 'N'
515: *
516: CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ALPHA,
1.6 ! bertrand 517: $ A( K+1 ), M+1, B, LDB )
1.1 bertrand 518: CALL DGEMM( 'T', 'N', K, N, K, -ONE, A( 0 ), M+1,
1.6 ! bertrand 519: $ B, LDB, ALPHA, B( K, 0 ), LDB )
1.1 bertrand 520: CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ONE,
1.6 ! bertrand 521: $ A( K ), M+1, B( K, 0 ), LDB )
1.1 bertrand 522: *
523: ELSE
524: *
525: * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U',
526: * and TRANS = 'T'
527: CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ALPHA,
1.6 ! bertrand 528: $ A( K ), M+1, B( K, 0 ), LDB )
1.1 bertrand 529: CALL DGEMM( 'N', 'N', K, N, K, -ONE, A( 0 ), M+1,
1.6 ! bertrand 530: $ B( K, 0 ), LDB, ALPHA, B, LDB )
1.1 bertrand 531: CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ONE,
1.6 ! bertrand 532: $ A( K+1 ), M+1, B, LDB )
1.1 bertrand 533: *
534: END IF
535: *
536: END IF
537: *
538: ELSE
539: *
540: * SIDE = 'L', N is even, and TRANSR = 'T'
541: *
542: IF( LOWER ) THEN
543: *
544: * SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'L'
545: *
546: IF( NOTRANS ) THEN
547: *
548: * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L',
549: * and TRANS = 'N'
550: *
551: CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ALPHA,
1.6 ! bertrand 552: $ A( K ), K, B, LDB )
1.1 bertrand 553: CALL DGEMM( 'T', 'N', K, N, K, -ONE,
1.6 ! bertrand 554: $ A( K*( K+1 ) ), K, B, LDB, ALPHA,
! 555: $ B( K, 0 ), LDB )
1.1 bertrand 556: CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ONE,
1.6 ! bertrand 557: $ A( 0 ), K, B( K, 0 ), LDB )
1.1 bertrand 558: *
559: ELSE
560: *
561: * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L',
562: * and TRANS = 'T'
563: *
564: CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ALPHA,
1.6 ! bertrand 565: $ A( 0 ), K, B( K, 0 ), LDB )
1.1 bertrand 566: CALL DGEMM( 'N', 'N', K, N, K, -ONE,
1.6 ! bertrand 567: $ A( K*( K+1 ) ), K, B( K, 0 ), LDB,
! 568: $ ALPHA, B, LDB )
1.1 bertrand 569: CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ONE,
1.6 ! bertrand 570: $ A( K ), K, B, LDB )
1.1 bertrand 571: *
572: END IF
573: *
574: ELSE
575: *
576: * SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'U'
577: *
578: IF( .NOT.NOTRANS ) THEN
579: *
580: * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U',
581: * and TRANS = 'N'
582: *
583: CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ALPHA,
1.6 ! bertrand 584: $ A( K*( K+1 ) ), K, B, LDB )
1.1 bertrand 585: CALL DGEMM( 'N', 'N', K, N, K, -ONE, A( 0 ), K, B,
1.6 ! bertrand 586: $ LDB, ALPHA, B( K, 0 ), LDB )
1.1 bertrand 587: CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ONE,
1.6 ! bertrand 588: $ A( K*K ), K, B( K, 0 ), LDB )
1.1 bertrand 589: *
590: ELSE
591: *
592: * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U',
593: * and TRANS = 'T'
594: *
595: CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ALPHA,
1.6 ! bertrand 596: $ A( K*K ), K, B( K, 0 ), LDB )
1.1 bertrand 597: CALL DGEMM( 'T', 'N', K, N, K, -ONE, A( 0 ), K,
1.6 ! bertrand 598: $ B( K, 0 ), LDB, ALPHA, B, LDB )
1.1 bertrand 599: CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ONE,
1.6 ! bertrand 600: $ A( K*( K+1 ) ), K, B, LDB )
1.1 bertrand 601: *
602: END IF
603: *
604: END IF
605: *
606: END IF
607: *
608: END IF
609: *
610: ELSE
611: *
612: * SIDE = 'R'
613: *
614: * A is N-by-N.
615: * If N is odd, set NISODD = .TRUE., and N1 and N2.
616: * If N is even, NISODD = .FALSE., and K.
617: *
618: IF( MOD( N, 2 ).EQ.0 ) THEN
619: NISODD = .FALSE.
620: K = N / 2
621: ELSE
622: NISODD = .TRUE.
623: IF( LOWER ) THEN
624: N2 = N / 2
625: N1 = N - N2
626: ELSE
627: N1 = N / 2
628: N2 = N - N1
629: END IF
630: END IF
631: *
632: IF( NISODD ) THEN
633: *
634: * SIDE = 'R' and N is odd
635: *
636: IF( NORMALTRANSR ) THEN
637: *
638: * SIDE = 'R', N is odd, and TRANSR = 'N'
639: *
640: IF( LOWER ) THEN
641: *
642: * SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'L'
643: *
644: IF( NOTRANS ) THEN
645: *
646: * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
647: * TRANS = 'N'
648: *
649: CALL DTRSM( 'R', 'U', 'T', DIAG, M, N2, ALPHA,
1.6 ! bertrand 650: $ A( N ), N, B( 0, N1 ), LDB )
1.1 bertrand 651: CALL DGEMM( 'N', 'N', M, N1, N2, -ONE, B( 0, N1 ),
1.6 ! bertrand 652: $ LDB, A( N1 ), N, ALPHA, B( 0, 0 ),
! 653: $ LDB )
1.1 bertrand 654: CALL DTRSM( 'R', 'L', 'N', DIAG, M, N1, ONE,
1.6 ! bertrand 655: $ A( 0 ), N, B( 0, 0 ), LDB )
1.1 bertrand 656: *
657: ELSE
658: *
659: * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
660: * TRANS = 'T'
661: *
662: CALL DTRSM( 'R', 'L', 'T', DIAG, M, N1, ALPHA,
1.6 ! bertrand 663: $ A( 0 ), N, B( 0, 0 ), LDB )
1.1 bertrand 664: CALL DGEMM( 'N', 'T', M, N2, N1, -ONE, B( 0, 0 ),
1.6 ! bertrand 665: $ LDB, A( N1 ), N, ALPHA, B( 0, N1 ),
! 666: $ LDB )
1.1 bertrand 667: CALL DTRSM( 'R', 'U', 'N', DIAG, M, N2, ONE,
1.6 ! bertrand 668: $ A( N ), N, B( 0, N1 ), LDB )
1.1 bertrand 669: *
670: END IF
671: *
672: ELSE
673: *
674: * SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'U'
675: *
676: IF( NOTRANS ) THEN
677: *
678: * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
679: * TRANS = 'N'
680: *
681: CALL DTRSM( 'R', 'L', 'T', DIAG, M, N1, ALPHA,
1.6 ! bertrand 682: $ A( N2 ), N, B( 0, 0 ), LDB )
1.1 bertrand 683: CALL DGEMM( 'N', 'N', M, N2, N1, -ONE, B( 0, 0 ),
1.6 ! bertrand 684: $ LDB, A( 0 ), N, ALPHA, B( 0, N1 ),
! 685: $ LDB )
1.1 bertrand 686: CALL DTRSM( 'R', 'U', 'N', DIAG, M, N2, ONE,
1.6 ! bertrand 687: $ A( N1 ), N, B( 0, N1 ), LDB )
1.1 bertrand 688: *
689: ELSE
690: *
691: * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
692: * TRANS = 'T'
693: *
694: CALL DTRSM( 'R', 'U', 'T', DIAG, M, N2, ALPHA,
1.6 ! bertrand 695: $ A( N1 ), N, B( 0, N1 ), LDB )
1.1 bertrand 696: CALL DGEMM( 'N', 'T', M, N1, N2, -ONE, B( 0, N1 ),
1.6 ! bertrand 697: $ LDB, A( 0 ), N, ALPHA, B( 0, 0 ), LDB )
1.1 bertrand 698: CALL DTRSM( 'R', 'L', 'N', DIAG, M, N1, ONE,
1.6 ! bertrand 699: $ A( N2 ), N, B( 0, 0 ), LDB )
1.1 bertrand 700: *
701: END IF
702: *
703: END IF
704: *
705: ELSE
706: *
707: * SIDE = 'R', N is odd, and TRANSR = 'T'
708: *
709: IF( LOWER ) THEN
710: *
711: * SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'L'
712: *
713: IF( NOTRANS ) THEN
714: *
715: * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and
716: * TRANS = 'N'
717: *
718: CALL DTRSM( 'R', 'L', 'N', DIAG, M, N2, ALPHA,
1.6 ! bertrand 719: $ A( 1 ), N1, B( 0, N1 ), LDB )
1.1 bertrand 720: CALL DGEMM( 'N', 'T', M, N1, N2, -ONE, B( 0, N1 ),
1.6 ! bertrand 721: $ LDB, A( N1*N1 ), N1, ALPHA, B( 0, 0 ),
! 722: $ LDB )
1.1 bertrand 723: CALL DTRSM( 'R', 'U', 'T', DIAG, M, N1, ONE,
1.6 ! bertrand 724: $ A( 0 ), N1, B( 0, 0 ), LDB )
1.1 bertrand 725: *
726: ELSE
727: *
728: * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and
729: * TRANS = 'T'
730: *
731: CALL DTRSM( 'R', 'U', 'N', DIAG, M, N1, ALPHA,
1.6 ! bertrand 732: $ A( 0 ), N1, B( 0, 0 ), LDB )
1.1 bertrand 733: CALL DGEMM( 'N', 'N', M, N2, N1, -ONE, B( 0, 0 ),
1.6 ! bertrand 734: $ LDB, A( N1*N1 ), N1, ALPHA, B( 0, N1 ),
! 735: $ LDB )
1.1 bertrand 736: CALL DTRSM( 'R', 'L', 'T', DIAG, M, N2, ONE,
1.6 ! bertrand 737: $ A( 1 ), N1, B( 0, N1 ), LDB )
1.1 bertrand 738: *
739: END IF
740: *
741: ELSE
742: *
743: * SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'U'
744: *
745: IF( NOTRANS ) THEN
746: *
747: * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and
748: * TRANS = 'N'
749: *
750: CALL DTRSM( 'R', 'U', 'N', DIAG, M, N1, ALPHA,
1.6 ! bertrand 751: $ A( N2*N2 ), N2, B( 0, 0 ), LDB )
1.1 bertrand 752: CALL DGEMM( 'N', 'T', M, N2, N1, -ONE, B( 0, 0 ),
1.6 ! bertrand 753: $ LDB, A( 0 ), N2, ALPHA, B( 0, N1 ),
! 754: $ LDB )
1.1 bertrand 755: CALL DTRSM( 'R', 'L', 'T', DIAG, M, N2, ONE,
1.6 ! bertrand 756: $ A( N1*N2 ), N2, B( 0, N1 ), LDB )
1.1 bertrand 757: *
758: ELSE
759: *
760: * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and
761: * TRANS = 'T'
762: *
763: CALL DTRSM( 'R', 'L', 'N', DIAG, M, N2, ALPHA,
1.6 ! bertrand 764: $ A( N1*N2 ), N2, B( 0, N1 ), LDB )
1.1 bertrand 765: CALL DGEMM( 'N', 'N', M, N1, N2, -ONE, B( 0, N1 ),
1.6 ! bertrand 766: $ LDB, A( 0 ), N2, ALPHA, B( 0, 0 ),
! 767: $ LDB )
1.1 bertrand 768: CALL DTRSM( 'R', 'U', 'T', DIAG, M, N1, ONE,
1.6 ! bertrand 769: $ A( N2*N2 ), N2, B( 0, 0 ), LDB )
1.1 bertrand 770: *
771: END IF
772: *
773: END IF
774: *
775: END IF
776: *
777: ELSE
778: *
779: * SIDE = 'R' and N is even
780: *
781: IF( NORMALTRANSR ) THEN
782: *
783: * SIDE = 'R', N is even, and TRANSR = 'N'
784: *
785: IF( LOWER ) THEN
786: *
787: * SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'L'
788: *
789: IF( NOTRANS ) THEN
790: *
791: * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L',
792: * and TRANS = 'N'
793: *
794: CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ALPHA,
1.6 ! bertrand 795: $ A( 0 ), N+1, B( 0, K ), LDB )
1.1 bertrand 796: CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, K ),
1.6 ! bertrand 797: $ LDB, A( K+1 ), N+1, ALPHA, B( 0, 0 ),
! 798: $ LDB )
1.1 bertrand 799: CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ONE,
1.6 ! bertrand 800: $ A( 1 ), N+1, B( 0, 0 ), LDB )
1.1 bertrand 801: *
802: ELSE
803: *
804: * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L',
805: * and TRANS = 'T'
806: *
807: CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ALPHA,
1.6 ! bertrand 808: $ A( 1 ), N+1, B( 0, 0 ), LDB )
1.1 bertrand 809: CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, 0 ),
1.6 ! bertrand 810: $ LDB, A( K+1 ), N+1, ALPHA, B( 0, K ),
! 811: $ LDB )
1.1 bertrand 812: CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ONE,
1.6 ! bertrand 813: $ A( 0 ), N+1, B( 0, K ), LDB )
1.1 bertrand 814: *
815: END IF
816: *
817: ELSE
818: *
819: * SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'U'
820: *
821: IF( NOTRANS ) THEN
822: *
823: * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U',
824: * and TRANS = 'N'
825: *
826: CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ALPHA,
1.6 ! bertrand 827: $ A( K+1 ), N+1, B( 0, 0 ), LDB )
1.1 bertrand 828: CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, 0 ),
1.6 ! bertrand 829: $ LDB, A( 0 ), N+1, ALPHA, B( 0, K ),
! 830: $ LDB )
1.1 bertrand 831: CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ONE,
1.6 ! bertrand 832: $ A( K ), N+1, B( 0, K ), LDB )
1.1 bertrand 833: *
834: ELSE
835: *
836: * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U',
837: * and TRANS = 'T'
838: *
839: CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ALPHA,
1.6 ! bertrand 840: $ A( K ), N+1, B( 0, K ), LDB )
1.1 bertrand 841: CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, K ),
1.6 ! bertrand 842: $ LDB, A( 0 ), N+1, ALPHA, B( 0, 0 ),
! 843: $ LDB )
1.1 bertrand 844: CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ONE,
1.6 ! bertrand 845: $ A( K+1 ), N+1, B( 0, 0 ), LDB )
1.1 bertrand 846: *
847: END IF
848: *
849: END IF
850: *
851: ELSE
852: *
853: * SIDE = 'R', N is even, and TRANSR = 'T'
854: *
855: IF( LOWER ) THEN
856: *
857: * SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'L'
858: *
859: IF( NOTRANS ) THEN
860: *
861: * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L',
862: * and TRANS = 'N'
863: *
864: CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ALPHA,
1.6 ! bertrand 865: $ A( 0 ), K, B( 0, K ), LDB )
1.1 bertrand 866: CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, K ),
1.6 ! bertrand 867: $ LDB, A( ( K+1 )*K ), K, ALPHA,
! 868: $ B( 0, 0 ), LDB )
1.1 bertrand 869: CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ONE,
1.6 ! bertrand 870: $ A( K ), K, B( 0, 0 ), LDB )
1.1 bertrand 871: *
872: ELSE
873: *
874: * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L',
875: * and TRANS = 'T'
876: *
877: CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ALPHA,
1.6 ! bertrand 878: $ A( K ), K, B( 0, 0 ), LDB )
1.1 bertrand 879: CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, 0 ),
1.6 ! bertrand 880: $ LDB, A( ( K+1 )*K ), K, ALPHA,
! 881: $ B( 0, K ), LDB )
1.1 bertrand 882: CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ONE,
1.6 ! bertrand 883: $ A( 0 ), K, B( 0, K ), LDB )
1.1 bertrand 884: *
885: END IF
886: *
887: ELSE
888: *
889: * SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'U'
890: *
891: IF( NOTRANS ) THEN
892: *
893: * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U',
894: * and TRANS = 'N'
895: *
896: CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ALPHA,
1.6 ! bertrand 897: $ A( ( K+1 )*K ), K, B( 0, 0 ), LDB )
1.1 bertrand 898: CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, 0 ),
1.6 ! bertrand 899: $ LDB, A( 0 ), K, ALPHA, B( 0, K ), LDB )
1.1 bertrand 900: CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ONE,
1.6 ! bertrand 901: $ A( K*K ), K, B( 0, K ), LDB )
1.1 bertrand 902: *
903: ELSE
904: *
905: * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U',
906: * and TRANS = 'T'
907: *
908: CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ALPHA,
1.6 ! bertrand 909: $ A( K*K ), K, B( 0, K ), LDB )
1.1 bertrand 910: CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, K ),
1.6 ! bertrand 911: $ LDB, A( 0 ), K, ALPHA, B( 0, 0 ), LDB )
1.1 bertrand 912: CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ONE,
1.6 ! bertrand 913: $ A( ( K+1 )*K ), K, B( 0, 0 ), LDB )
1.1 bertrand 914: *
915: END IF
916: *
917: END IF
918: *
919: END IF
920: *
921: END IF
922: END IF
923: *
924: RETURN
925: *
926: * End of DTFSM
927: *
928: END
CVSweb interface <joel.bertrand@systella.fr>