1: SUBROUTINE DTFSM( TRANSR, SIDE, UPLO, TRANS, DIAG, M, N, ALPHA, A,
2: + B, LDB )
3: *
4: * -- LAPACK routine (version 3.2.2) --
5: *
6: * -- Contributed by Fred Gustavson of the IBM Watson Research Center --
7: * -- June 2010 --
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: *
34: * op( A ) = A or op( A ) = A'.
35: *
36: * A is in Rectangular Full Packed (RFP) Format.
37: *
38: * The matrix X is overwritten on B.
39: *
40: * Arguments
41: * ==========
42: *
43: * TRANSR (input) CHARACTER
44: * = 'N': The Normal Form of RFP A is stored;
45: * = 'T': The Transpose Form of RFP A is stored.
46: *
47: * SIDE (input) CHARACTER
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: *
57: * UPLO (input) CHARACTER
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: *
65: * TRANS (input) CHARACTER
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: *
75: * DIAG (input) CHARACTER
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,
225: + NOTRANS
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' ) )
256: + THEN
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 ) )
273: + RETURN
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,
328: + A, M, B, LDB )
329: ELSE
330: CALL DTRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA,
331: + A( 0 ), M, B, LDB )
332: CALL DGEMM( 'N', 'N', M2, N, M1, -ONE, A( M1 ),
333: + M, B, LDB, ALPHA, B( M1, 0 ), LDB )
334: CALL DTRSM( 'L', 'U', 'T', DIAG, M2, N, ONE,
335: + A( M ), M, B( M1, 0 ), LDB )
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,
345: + A( 0 ), M, B, LDB )
346: ELSE
347: CALL DTRSM( 'L', 'U', 'N', DIAG, M2, N, ALPHA,
348: + A( M ), M, B( M1, 0 ), LDB )
349: CALL DGEMM( 'T', 'N', M1, N, M2, -ONE, A( M1 ),
350: + M, B( M1, 0 ), LDB, ALPHA, B, LDB )
351: CALL DTRSM( 'L', 'L', 'T', DIAG, M1, N, ONE,
352: + A( 0 ), M, B, LDB )
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,
367: + A( M2 ), M, B, LDB )
368: CALL DGEMM( 'T', 'N', M2, N, M1, -ONE, A( 0 ), M,
369: + B, LDB, ALPHA, B( M1, 0 ), LDB )
370: CALL DTRSM( 'L', 'U', 'T', DIAG, M2, N, ONE,
371: + A( M1 ), M, B( M1, 0 ), LDB )
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,
379: + A( M1 ), M, B( M1, 0 ), LDB )
380: CALL DGEMM( 'N', 'N', M1, N, M2, -ONE, A( 0 ), M,
381: + B( M1, 0 ), LDB, ALPHA, B, LDB )
382: CALL DTRSM( 'L', 'L', 'T', DIAG, M1, N, ONE,
383: + A( M2 ), M, B, LDB )
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,
404: + A( 0 ), M1, B, LDB )
405: ELSE
406: CALL DTRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA,
407: + A( 0 ), M1, B, LDB )
408: CALL DGEMM( 'T', 'N', M2, N, M1, -ONE,
409: + A( M1*M1 ), M1, B, LDB, ALPHA,
410: + B( M1, 0 ), LDB )
411: CALL DTRSM( 'L', 'L', 'N', DIAG, M2, N, ONE,
412: + A( 1 ), M1, B( M1, 0 ), LDB )
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,
422: + A( 0 ), M1, B, LDB )
423: ELSE
424: CALL DTRSM( 'L', 'L', 'T', DIAG, M2, N, ALPHA,
425: + A( 1 ), M1, B( M1, 0 ), LDB )
426: CALL DGEMM( 'N', 'N', M1, N, M2, -ONE,
427: + A( M1*M1 ), M1, B( M1, 0 ), LDB,
428: + ALPHA, B, LDB )
429: CALL DTRSM( 'L', 'U', 'N', DIAG, M1, N, ONE,
430: + A( 0 ), M1, B, LDB )
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,
445: + A( M2*M2 ), M2, B, LDB )
446: CALL DGEMM( 'N', 'N', M2, N, M1, -ONE, A( 0 ), M2,
447: + B, LDB, ALPHA, B( M1, 0 ), LDB )
448: CALL DTRSM( 'L', 'L', 'N', DIAG, M2, N, ONE,
449: + A( M1*M2 ), M2, B( M1, 0 ), LDB )
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,
457: + A( M1*M2 ), M2, B( M1, 0 ), LDB )
458: CALL DGEMM( 'T', 'N', M1, N, M2, -ONE, A( 0 ), M2,
459: + B( M1, 0 ), LDB, ALPHA, B, LDB )
460: CALL DTRSM( 'L', 'U', 'N', DIAG, M1, N, ONE,
461: + A( M2*M2 ), M2, B, LDB )
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,
487: + A( 1 ), M+1, B, LDB )
488: CALL DGEMM( 'N', 'N', K, N, K, -ONE, A( K+1 ),
489: + M+1, B, LDB, ALPHA, B( K, 0 ), LDB )
490: CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ONE,
491: + A( 0 ), M+1, B( K, 0 ), LDB )
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,
499: + A( 0 ), M+1, B( K, 0 ), LDB )
500: CALL DGEMM( 'T', 'N', K, N, K, -ONE, A( K+1 ),
501: + M+1, B( K, 0 ), LDB, ALPHA, B, LDB )
502: CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ONE,
503: + A( 1 ), M+1, B, LDB )
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,
517: + A( K+1 ), M+1, B, LDB )
518: CALL DGEMM( 'T', 'N', K, N, K, -ONE, A( 0 ), M+1,
519: + B, LDB, ALPHA, B( K, 0 ), LDB )
520: CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ONE,
521: + A( K ), M+1, B( K, 0 ), LDB )
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,
528: + A( K ), M+1, B( K, 0 ), LDB )
529: CALL DGEMM( 'N', 'N', K, N, K, -ONE, A( 0 ), M+1,
530: + B( K, 0 ), LDB, ALPHA, B, LDB )
531: CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ONE,
532: + A( K+1 ), M+1, B, LDB )
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,
552: + A( K ), K, B, LDB )
553: CALL DGEMM( 'T', 'N', K, N, K, -ONE,
554: + A( K*( K+1 ) ), K, B, LDB, ALPHA,
555: + B( K, 0 ), LDB )
556: CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ONE,
557: + A( 0 ), K, B( K, 0 ), LDB )
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,
565: + A( 0 ), K, B( K, 0 ), LDB )
566: CALL DGEMM( 'N', 'N', K, N, K, -ONE,
567: + A( K*( K+1 ) ), K, B( K, 0 ), LDB,
568: + ALPHA, B, LDB )
569: CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ONE,
570: + A( K ), K, B, LDB )
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,
584: + A( K*( K+1 ) ), K, B, LDB )
585: CALL DGEMM( 'N', 'N', K, N, K, -ONE, A( 0 ), K, B,
586: + LDB, ALPHA, B( K, 0 ), LDB )
587: CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ONE,
588: + A( K*K ), K, B( K, 0 ), LDB )
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,
596: + A( K*K ), K, B( K, 0 ), LDB )
597: CALL DGEMM( 'T', 'N', K, N, K, -ONE, A( 0 ), K,
598: + B( K, 0 ), LDB, ALPHA, B, LDB )
599: CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ONE,
600: + A( K*( K+1 ) ), K, B, LDB )
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,
650: + A( N ), N, B( 0, N1 ), LDB )
651: CALL DGEMM( 'N', 'N', M, N1, N2, -ONE, B( 0, N1 ),
652: + LDB, A( N1 ), N, ALPHA, B( 0, 0 ),
653: + LDB )
654: CALL DTRSM( 'R', 'L', 'N', DIAG, M, N1, ONE,
655: + A( 0 ), N, B( 0, 0 ), LDB )
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,
663: + A( 0 ), N, B( 0, 0 ), LDB )
664: CALL DGEMM( 'N', 'T', M, N2, N1, -ONE, B( 0, 0 ),
665: + LDB, A( N1 ), N, ALPHA, B( 0, N1 ),
666: + LDB )
667: CALL DTRSM( 'R', 'U', 'N', DIAG, M, N2, ONE,
668: + A( N ), N, B( 0, N1 ), LDB )
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,
682: + A( N2 ), N, B( 0, 0 ), LDB )
683: CALL DGEMM( 'N', 'N', M, N2, N1, -ONE, B( 0, 0 ),
684: + LDB, A( 0 ), N, ALPHA, B( 0, N1 ),
685: + LDB )
686: CALL DTRSM( 'R', 'U', 'N', DIAG, M, N2, ONE,
687: + A( N1 ), N, B( 0, N1 ), LDB )
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,
695: + A( N1 ), N, B( 0, N1 ), LDB )
696: CALL DGEMM( 'N', 'T', M, N1, N2, -ONE, B( 0, N1 ),
697: + LDB, A( 0 ), N, ALPHA, B( 0, 0 ), LDB )
698: CALL DTRSM( 'R', 'L', 'N', DIAG, M, N1, ONE,
699: + A( N2 ), N, B( 0, 0 ), LDB )
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,
719: + A( 1 ), N1, B( 0, N1 ), LDB )
720: CALL DGEMM( 'N', 'T', M, N1, N2, -ONE, B( 0, N1 ),
721: + LDB, A( N1*N1 ), N1, ALPHA, B( 0, 0 ),
722: + LDB )
723: CALL DTRSM( 'R', 'U', 'T', DIAG, M, N1, ONE,
724: + A( 0 ), N1, B( 0, 0 ), LDB )
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,
732: + A( 0 ), N1, B( 0, 0 ), LDB )
733: CALL DGEMM( 'N', 'N', M, N2, N1, -ONE, B( 0, 0 ),
734: + LDB, A( N1*N1 ), N1, ALPHA, B( 0, N1 ),
735: + LDB )
736: CALL DTRSM( 'R', 'L', 'T', DIAG, M, N2, ONE,
737: + A( 1 ), N1, B( 0, N1 ), LDB )
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,
751: + A( N2*N2 ), N2, B( 0, 0 ), LDB )
752: CALL DGEMM( 'N', 'T', M, N2, N1, -ONE, B( 0, 0 ),
753: + LDB, A( 0 ), N2, ALPHA, B( 0, N1 ),
754: + LDB )
755: CALL DTRSM( 'R', 'L', 'T', DIAG, M, N2, ONE,
756: + A( N1*N2 ), N2, B( 0, N1 ), LDB )
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,
764: + A( N1*N2 ), N2, B( 0, N1 ), LDB )
765: CALL DGEMM( 'N', 'N', M, N1, N2, -ONE, B( 0, N1 ),
766: + LDB, A( 0 ), N2, ALPHA, B( 0, 0 ),
767: + LDB )
768: CALL DTRSM( 'R', 'U', 'T', DIAG, M, N1, ONE,
769: + A( N2*N2 ), N2, B( 0, 0 ), LDB )
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,
795: + A( 0 ), N+1, B( 0, K ), LDB )
796: CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, K ),
797: + LDB, A( K+1 ), N+1, ALPHA, B( 0, 0 ),
798: + LDB )
799: CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ONE,
800: + A( 1 ), N+1, B( 0, 0 ), LDB )
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,
808: + A( 1 ), N+1, B( 0, 0 ), LDB )
809: CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, 0 ),
810: + LDB, A( K+1 ), N+1, ALPHA, B( 0, K ),
811: + LDB )
812: CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ONE,
813: + A( 0 ), N+1, B( 0, K ), LDB )
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,
827: + A( K+1 ), N+1, B( 0, 0 ), LDB )
828: CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, 0 ),
829: + LDB, A( 0 ), N+1, ALPHA, B( 0, K ),
830: + LDB )
831: CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ONE,
832: + A( K ), N+1, B( 0, K ), LDB )
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,
840: + A( K ), N+1, B( 0, K ), LDB )
841: CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, K ),
842: + LDB, A( 0 ), N+1, ALPHA, B( 0, 0 ),
843: + LDB )
844: CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ONE,
845: + A( K+1 ), N+1, B( 0, 0 ), LDB )
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,
865: + A( 0 ), K, B( 0, K ), LDB )
866: CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, K ),
867: + LDB, A( ( K+1 )*K ), K, ALPHA,
868: + B( 0, 0 ), LDB )
869: CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ONE,
870: + A( K ), K, B( 0, 0 ), LDB )
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,
878: + A( K ), K, B( 0, 0 ), LDB )
879: CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, 0 ),
880: + LDB, A( ( K+1 )*K ), K, ALPHA,
881: + B( 0, K ), LDB )
882: CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ONE,
883: + A( 0 ), K, B( 0, K ), LDB )
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,
897: + A( ( K+1 )*K ), K, B( 0, 0 ), LDB )
898: CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, 0 ),
899: + LDB, A( 0 ), K, ALPHA, B( 0, K ), LDB )
900: CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ONE,
901: + A( K*K ), K, B( 0, K ), LDB )
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,
909: + A( K*K ), K, B( 0, K ), LDB )
910: CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, K ),
911: + LDB, A( 0 ), K, ALPHA, B( 0, 0 ), LDB )
912: CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ONE,
913: + A( ( K+1 )*K ), K, B( 0, 0 ), LDB )
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>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtfsm.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack