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