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: *> \date September 2012
186: *
187: *> \ingroup doubleOTHERcomputational
188: *
189: *> \par Further Details:
190: * =====================
191: *>
192: *> \verbatim
193: *>
194: *> We first consider Rectangular Full Packed (RFP) Format when N is
195: *> even. We give an example where N = 6.
196: *>
197: *> AP is Upper AP is Lower
198: *>
199: *> 00 01 02 03 04 05 00
200: *> 11 12 13 14 15 10 11
201: *> 22 23 24 25 20 21 22
202: *> 33 34 35 30 31 32 33
203: *> 44 45 40 41 42 43 44
204: *> 55 50 51 52 53 54 55
205: *>
206: *>
207: *> Let TRANSR = 'N'. RFP holds AP as follows:
208: *> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
209: *> three columns of AP upper. The lower triangle A(4:6,0:2) consists of
210: *> the transpose of the first three columns of AP upper.
211: *> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
212: *> three columns of AP lower. The upper triangle A(0:2,0:2) consists of
213: *> the transpose of the last three columns of AP lower.
214: *> This covers the case N even and TRANSR = 'N'.
215: *>
216: *> RFP A RFP A
217: *>
218: *> 03 04 05 33 43 53
219: *> 13 14 15 00 44 54
220: *> 23 24 25 10 11 55
221: *> 33 34 35 20 21 22
222: *> 00 44 45 30 31 32
223: *> 01 11 55 40 41 42
224: *> 02 12 22 50 51 52
225: *>
226: *> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
227: *> transpose of RFP A above. One therefore gets:
228: *>
229: *>
230: *> RFP A RFP A
231: *>
232: *> 03 13 23 33 00 01 02 33 00 10 20 30 40 50
233: *> 04 14 24 34 44 11 12 43 44 11 21 31 41 51
234: *> 05 15 25 35 45 55 22 53 54 55 22 32 42 52
235: *>
236: *>
237: *> We then consider Rectangular Full Packed (RFP) Format when N is
238: *> odd. We give an example where N = 5.
239: *>
240: *> AP is Upper AP is Lower
241: *>
242: *> 00 01 02 03 04 00
243: *> 11 12 13 14 10 11
244: *> 22 23 24 20 21 22
245: *> 33 34 30 31 32 33
246: *> 44 40 41 42 43 44
247: *>
248: *>
249: *> Let TRANSR = 'N'. RFP holds AP as follows:
250: *> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
251: *> three columns of AP upper. The lower triangle A(3:4,0:1) consists of
252: *> the transpose of the first two columns of AP upper.
253: *> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
254: *> three columns of AP lower. The upper triangle A(0:1,1:2) consists of
255: *> the transpose of the last two columns of AP lower.
256: *> This covers the case N odd and TRANSR = 'N'.
257: *>
258: *> RFP A RFP A
259: *>
260: *> 02 03 04 00 33 43
261: *> 12 13 14 10 11 44
262: *> 22 23 24 20 21 22
263: *> 00 33 34 30 31 32
264: *> 01 11 44 40 41 42
265: *>
266: *> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
267: *> transpose of RFP A above. One therefore gets:
268: *>
269: *> RFP A RFP A
270: *>
271: *> 02 12 22 00 01 00 10 20 30 40 50
272: *> 03 13 23 33 11 33 11 21 31 41 51
273: *> 04 14 24 34 44 43 44 22 32 42 52
274: *> \endverbatim
275: *
276: * =====================================================================
277: SUBROUTINE DTFSM( TRANSR, SIDE, UPLO, TRANS, DIAG, M, N, ALPHA, A,
278: $ B, LDB )
279: *
280: * -- LAPACK computational routine (version 3.4.2) --
281: * -- LAPACK is a software package provided by Univ. of Tennessee, --
282: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
283: * September 2012
284: *
285: * .. Scalar Arguments ..
286: CHARACTER TRANSR, DIAG, SIDE, TRANS, UPLO
287: INTEGER LDB, M, N
288: DOUBLE PRECISION ALPHA
289: * ..
290: * .. Array Arguments ..
291: DOUBLE PRECISION A( 0: * ), B( 0: LDB-1, 0: * )
292: * ..
293: *
294: * =====================================================================
295: *
296: * ..
297: * .. Parameters ..
298: DOUBLE PRECISION ONE, ZERO
299: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
300: * ..
301: * .. Local Scalars ..
302: LOGICAL LOWER, LSIDE, MISODD, NISODD, NORMALTRANSR,
303: $ NOTRANS
304: INTEGER M1, M2, N1, N2, K, INFO, I, J
305: * ..
306: * .. External Functions ..
307: LOGICAL LSAME
308: EXTERNAL LSAME
309: * ..
310: * .. External Subroutines ..
311: EXTERNAL XERBLA, DGEMM, DTRSM
312: * ..
313: * .. Intrinsic Functions ..
314: INTRINSIC MAX, MOD
315: * ..
316: * .. Executable Statements ..
317: *
318: * Test the input parameters.
319: *
320: INFO = 0
321: NORMALTRANSR = LSAME( TRANSR, 'N' )
322: LSIDE = LSAME( SIDE, 'L' )
323: LOWER = LSAME( UPLO, 'L' )
324: NOTRANS = LSAME( TRANS, 'N' )
325: IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
326: INFO = -1
327: ELSE IF( .NOT.LSIDE .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
328: INFO = -2
329: ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
330: INFO = -3
331: ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
332: INFO = -4
333: ELSE IF( .NOT.LSAME( DIAG, 'N' ) .AND. .NOT.LSAME( DIAG, 'U' ) )
334: $ THEN
335: INFO = -5
336: ELSE IF( M.LT.0 ) THEN
337: INFO = -6
338: ELSE IF( N.LT.0 ) THEN
339: INFO = -7
340: ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
341: INFO = -11
342: END IF
343: IF( INFO.NE.0 ) THEN
344: CALL XERBLA( 'DTFSM ', -INFO )
345: RETURN
346: END IF
347: *
348: * Quick return when ( (N.EQ.0).OR.(M.EQ.0) )
349: *
350: IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )
351: $ RETURN
352: *
353: * Quick return when ALPHA.EQ.(0D+0)
354: *
355: IF( ALPHA.EQ.ZERO ) THEN
356: DO 20 J = 0, N - 1
357: DO 10 I = 0, M - 1
358: B( I, J ) = ZERO
359: 10 CONTINUE
360: 20 CONTINUE
361: RETURN
362: END IF
363: *
364: IF( LSIDE ) THEN
365: *
366: * SIDE = 'L'
367: *
368: * A is M-by-M.
369: * If M is odd, set NISODD = .TRUE., and M1 and M2.
370: * If M is even, NISODD = .FALSE., and M.
371: *
372: IF( MOD( M, 2 ).EQ.0 ) THEN
373: MISODD = .FALSE.
374: K = M / 2
375: ELSE
376: MISODD = .TRUE.
377: IF( LOWER ) THEN
378: M2 = M / 2
379: M1 = M - M2
380: ELSE
381: M1 = M / 2
382: M2 = M - M1
383: END IF
384: END IF
385: *
386: *
387: IF( MISODD ) THEN
388: *
389: * SIDE = 'L' and N is odd
390: *
391: IF( NORMALTRANSR ) THEN
392: *
393: * SIDE = 'L', N is odd, and TRANSR = 'N'
394: *
395: IF( LOWER ) THEN
396: *
397: * SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'L'
398: *
399: IF( NOTRANS ) THEN
400: *
401: * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
402: * TRANS = 'N'
403: *
404: IF( M.EQ.1 ) THEN
405: CALL DTRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA,
406: $ A, M, B, LDB )
407: ELSE
408: CALL DTRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA,
409: $ A( 0 ), M, B, LDB )
410: CALL DGEMM( 'N', 'N', M2, N, M1, -ONE, A( M1 ),
411: $ M, B, LDB, ALPHA, B( M1, 0 ), LDB )
412: CALL DTRSM( 'L', 'U', 'T', DIAG, M2, N, ONE,
413: $ A( M ), M, B( M1, 0 ), LDB )
414: END IF
415: *
416: ELSE
417: *
418: * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
419: * TRANS = 'T'
420: *
421: IF( M.EQ.1 ) THEN
422: CALL DTRSM( 'L', 'L', 'T', DIAG, M1, N, ALPHA,
423: $ A( 0 ), M, B, LDB )
424: ELSE
425: CALL DTRSM( 'L', 'U', 'N', DIAG, M2, N, ALPHA,
426: $ A( M ), M, B( M1, 0 ), LDB )
427: CALL DGEMM( 'T', 'N', M1, N, M2, -ONE, A( M1 ),
428: $ M, B( M1, 0 ), LDB, ALPHA, B, LDB )
429: CALL DTRSM( 'L', 'L', 'T', DIAG, M1, N, ONE,
430: $ A( 0 ), M, B, LDB )
431: END IF
432: *
433: END IF
434: *
435: ELSE
436: *
437: * SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'U'
438: *
439: IF( .NOT.NOTRANS ) THEN
440: *
441: * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
442: * TRANS = 'N'
443: *
444: CALL DTRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA,
445: $ A( M2 ), M, B, LDB )
446: CALL DGEMM( 'T', 'N', M2, N, M1, -ONE, A( 0 ), M,
447: $ B, LDB, ALPHA, B( M1, 0 ), LDB )
448: CALL DTRSM( 'L', 'U', 'T', DIAG, M2, N, ONE,
449: $ A( M1 ), M, B( M1, 0 ), LDB )
450: *
451: ELSE
452: *
453: * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
454: * TRANS = 'T'
455: *
456: CALL DTRSM( 'L', 'U', 'N', DIAG, M2, N, ALPHA,
457: $ A( M1 ), M, B( M1, 0 ), LDB )
458: CALL DGEMM( 'N', 'N', M1, N, M2, -ONE, A( 0 ), M,
459: $ B( M1, 0 ), LDB, ALPHA, B, LDB )
460: CALL DTRSM( 'L', 'L', 'T', DIAG, M1, N, ONE,
461: $ A( M2 ), M, B, LDB )
462: *
463: END IF
464: *
465: END IF
466: *
467: ELSE
468: *
469: * SIDE = 'L', N is odd, and TRANSR = 'T'
470: *
471: IF( LOWER ) THEN
472: *
473: * SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'L'
474: *
475: IF( NOTRANS ) THEN
476: *
477: * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and
478: * TRANS = 'N'
479: *
480: IF( M.EQ.1 ) THEN
481: CALL DTRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA,
482: $ A( 0 ), M1, B, LDB )
483: ELSE
484: CALL DTRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA,
485: $ A( 0 ), M1, B, LDB )
486: CALL DGEMM( 'T', 'N', M2, N, M1, -ONE,
487: $ A( M1*M1 ), M1, B, LDB, ALPHA,
488: $ B( M1, 0 ), LDB )
489: CALL DTRSM( 'L', 'L', 'N', DIAG, M2, N, ONE,
490: $ A( 1 ), M1, B( M1, 0 ), LDB )
491: END IF
492: *
493: ELSE
494: *
495: * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and
496: * TRANS = 'T'
497: *
498: IF( M.EQ.1 ) THEN
499: CALL DTRSM( 'L', 'U', 'N', DIAG, M1, N, ALPHA,
500: $ A( 0 ), M1, B, LDB )
501: ELSE
502: CALL DTRSM( 'L', 'L', 'T', DIAG, M2, N, ALPHA,
503: $ A( 1 ), M1, B( M1, 0 ), LDB )
504: CALL DGEMM( 'N', 'N', M1, N, M2, -ONE,
505: $ A( M1*M1 ), M1, B( M1, 0 ), LDB,
506: $ ALPHA, B, LDB )
507: CALL DTRSM( 'L', 'U', 'N', DIAG, M1, N, ONE,
508: $ A( 0 ), M1, B, LDB )
509: END IF
510: *
511: END IF
512: *
513: ELSE
514: *
515: * SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'U'
516: *
517: IF( .NOT.NOTRANS ) THEN
518: *
519: * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and
520: * TRANS = 'N'
521: *
522: CALL DTRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA,
523: $ A( M2*M2 ), M2, B, LDB )
524: CALL DGEMM( 'N', 'N', M2, N, M1, -ONE, A( 0 ), M2,
525: $ B, LDB, ALPHA, B( M1, 0 ), LDB )
526: CALL DTRSM( 'L', 'L', 'N', DIAG, M2, N, ONE,
527: $ A( M1*M2 ), M2, B( M1, 0 ), LDB )
528: *
529: ELSE
530: *
531: * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and
532: * TRANS = 'T'
533: *
534: CALL DTRSM( 'L', 'L', 'T', DIAG, M2, N, ALPHA,
535: $ A( M1*M2 ), M2, B( M1, 0 ), LDB )
536: CALL DGEMM( 'T', 'N', M1, N, M2, -ONE, A( 0 ), M2,
537: $ B( M1, 0 ), LDB, ALPHA, B, LDB )
538: CALL DTRSM( 'L', 'U', 'N', DIAG, M1, N, ONE,
539: $ A( M2*M2 ), M2, B, LDB )
540: *
541: END IF
542: *
543: END IF
544: *
545: END IF
546: *
547: ELSE
548: *
549: * SIDE = 'L' and N is even
550: *
551: IF( NORMALTRANSR ) THEN
552: *
553: * SIDE = 'L', N is even, and TRANSR = 'N'
554: *
555: IF( LOWER ) THEN
556: *
557: * SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'L'
558: *
559: IF( NOTRANS ) THEN
560: *
561: * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L',
562: * and TRANS = 'N'
563: *
564: CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ALPHA,
565: $ A( 1 ), M+1, B, LDB )
566: CALL DGEMM( 'N', 'N', K, N, K, -ONE, A( K+1 ),
567: $ M+1, B, LDB, ALPHA, B( K, 0 ), LDB )
568: CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ONE,
569: $ A( 0 ), M+1, B( K, 0 ), LDB )
570: *
571: ELSE
572: *
573: * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L',
574: * and TRANS = 'T'
575: *
576: CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ALPHA,
577: $ A( 0 ), M+1, B( K, 0 ), LDB )
578: CALL DGEMM( 'T', 'N', K, N, K, -ONE, A( K+1 ),
579: $ M+1, B( K, 0 ), LDB, ALPHA, B, LDB )
580: CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ONE,
581: $ A( 1 ), M+1, B, LDB )
582: *
583: END IF
584: *
585: ELSE
586: *
587: * SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'U'
588: *
589: IF( .NOT.NOTRANS ) THEN
590: *
591: * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U',
592: * and TRANS = 'N'
593: *
594: CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ALPHA,
595: $ A( K+1 ), M+1, B, LDB )
596: CALL DGEMM( 'T', 'N', K, N, K, -ONE, A( 0 ), M+1,
597: $ B, LDB, ALPHA, B( K, 0 ), LDB )
598: CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ONE,
599: $ A( K ), M+1, B( K, 0 ), LDB )
600: *
601: ELSE
602: *
603: * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U',
604: * and TRANS = 'T'
605: CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ALPHA,
606: $ A( K ), M+1, B( K, 0 ), LDB )
607: CALL DGEMM( 'N', 'N', K, N, K, -ONE, A( 0 ), M+1,
608: $ B( K, 0 ), LDB, ALPHA, B, LDB )
609: CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ONE,
610: $ A( K+1 ), M+1, B, LDB )
611: *
612: END IF
613: *
614: END IF
615: *
616: ELSE
617: *
618: * SIDE = 'L', N is even, and TRANSR = 'T'
619: *
620: IF( LOWER ) THEN
621: *
622: * SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'L'
623: *
624: IF( NOTRANS ) THEN
625: *
626: * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L',
627: * and TRANS = 'N'
628: *
629: CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ALPHA,
630: $ A( K ), K, B, LDB )
631: CALL DGEMM( 'T', 'N', K, N, K, -ONE,
632: $ A( K*( K+1 ) ), K, B, LDB, ALPHA,
633: $ B( K, 0 ), LDB )
634: CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ONE,
635: $ A( 0 ), K, B( K, 0 ), LDB )
636: *
637: ELSE
638: *
639: * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L',
640: * and TRANS = 'T'
641: *
642: CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ALPHA,
643: $ A( 0 ), K, B( K, 0 ), LDB )
644: CALL DGEMM( 'N', 'N', K, N, K, -ONE,
645: $ A( K*( K+1 ) ), K, B( K, 0 ), LDB,
646: $ ALPHA, B, LDB )
647: CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ONE,
648: $ A( K ), K, B, LDB )
649: *
650: END IF
651: *
652: ELSE
653: *
654: * SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'U'
655: *
656: IF( .NOT.NOTRANS ) THEN
657: *
658: * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U',
659: * and TRANS = 'N'
660: *
661: CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ALPHA,
662: $ A( K*( K+1 ) ), K, B, LDB )
663: CALL DGEMM( 'N', 'N', K, N, K, -ONE, A( 0 ), K, B,
664: $ LDB, ALPHA, B( K, 0 ), LDB )
665: CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ONE,
666: $ A( K*K ), K, B( K, 0 ), LDB )
667: *
668: ELSE
669: *
670: * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U',
671: * and TRANS = 'T'
672: *
673: CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ALPHA,
674: $ A( K*K ), K, B( K, 0 ), LDB )
675: CALL DGEMM( 'T', 'N', K, N, K, -ONE, A( 0 ), K,
676: $ B( K, 0 ), LDB, ALPHA, B, LDB )
677: CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ONE,
678: $ A( K*( K+1 ) ), K, B, LDB )
679: *
680: END IF
681: *
682: END IF
683: *
684: END IF
685: *
686: END IF
687: *
688: ELSE
689: *
690: * SIDE = 'R'
691: *
692: * A is N-by-N.
693: * If N is odd, set NISODD = .TRUE., and N1 and N2.
694: * If N is even, NISODD = .FALSE., and K.
695: *
696: IF( MOD( N, 2 ).EQ.0 ) THEN
697: NISODD = .FALSE.
698: K = N / 2
699: ELSE
700: NISODD = .TRUE.
701: IF( LOWER ) THEN
702: N2 = N / 2
703: N1 = N - N2
704: ELSE
705: N1 = N / 2
706: N2 = N - N1
707: END IF
708: END IF
709: *
710: IF( NISODD ) THEN
711: *
712: * SIDE = 'R' and N is odd
713: *
714: IF( NORMALTRANSR ) THEN
715: *
716: * SIDE = 'R', N is odd, and TRANSR = 'N'
717: *
718: IF( LOWER ) THEN
719: *
720: * SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'L'
721: *
722: IF( NOTRANS ) THEN
723: *
724: * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
725: * TRANS = 'N'
726: *
727: CALL DTRSM( 'R', 'U', 'T', DIAG, M, N2, ALPHA,
728: $ A( N ), N, B( 0, N1 ), LDB )
729: CALL DGEMM( 'N', 'N', M, N1, N2, -ONE, B( 0, N1 ),
730: $ LDB, A( N1 ), N, ALPHA, B( 0, 0 ),
731: $ LDB )
732: CALL DTRSM( 'R', 'L', 'N', DIAG, M, N1, ONE,
733: $ A( 0 ), N, B( 0, 0 ), LDB )
734: *
735: ELSE
736: *
737: * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
738: * TRANS = 'T'
739: *
740: CALL DTRSM( 'R', 'L', 'T', DIAG, M, N1, ALPHA,
741: $ A( 0 ), N, B( 0, 0 ), LDB )
742: CALL DGEMM( 'N', 'T', M, N2, N1, -ONE, B( 0, 0 ),
743: $ LDB, A( N1 ), N, ALPHA, B( 0, N1 ),
744: $ LDB )
745: CALL DTRSM( 'R', 'U', 'N', DIAG, M, N2, ONE,
746: $ A( N ), N, B( 0, N1 ), LDB )
747: *
748: END IF
749: *
750: ELSE
751: *
752: * SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'U'
753: *
754: IF( NOTRANS ) THEN
755: *
756: * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
757: * TRANS = 'N'
758: *
759: CALL DTRSM( 'R', 'L', 'T', DIAG, M, N1, ALPHA,
760: $ A( N2 ), N, B( 0, 0 ), LDB )
761: CALL DGEMM( 'N', 'N', M, N2, N1, -ONE, B( 0, 0 ),
762: $ LDB, A( 0 ), N, ALPHA, B( 0, N1 ),
763: $ LDB )
764: CALL DTRSM( 'R', 'U', 'N', DIAG, M, N2, ONE,
765: $ A( N1 ), N, B( 0, N1 ), LDB )
766: *
767: ELSE
768: *
769: * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
770: * TRANS = 'T'
771: *
772: CALL DTRSM( 'R', 'U', 'T', DIAG, M, N2, ALPHA,
773: $ A( N1 ), N, B( 0, N1 ), LDB )
774: CALL DGEMM( 'N', 'T', M, N1, N2, -ONE, B( 0, N1 ),
775: $ LDB, A( 0 ), N, ALPHA, B( 0, 0 ), LDB )
776: CALL DTRSM( 'R', 'L', 'N', DIAG, M, N1, ONE,
777: $ A( N2 ), N, B( 0, 0 ), LDB )
778: *
779: END IF
780: *
781: END IF
782: *
783: ELSE
784: *
785: * SIDE = 'R', N is odd, and TRANSR = 'T'
786: *
787: IF( LOWER ) THEN
788: *
789: * SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'L'
790: *
791: IF( NOTRANS ) THEN
792: *
793: * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and
794: * TRANS = 'N'
795: *
796: CALL DTRSM( 'R', 'L', 'N', DIAG, M, N2, ALPHA,
797: $ A( 1 ), N1, B( 0, N1 ), LDB )
798: CALL DGEMM( 'N', 'T', M, N1, N2, -ONE, B( 0, N1 ),
799: $ LDB, A( N1*N1 ), N1, ALPHA, B( 0, 0 ),
800: $ LDB )
801: CALL DTRSM( 'R', 'U', 'T', DIAG, M, N1, ONE,
802: $ A( 0 ), N1, B( 0, 0 ), LDB )
803: *
804: ELSE
805: *
806: * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and
807: * TRANS = 'T'
808: *
809: CALL DTRSM( 'R', 'U', 'N', DIAG, M, N1, ALPHA,
810: $ A( 0 ), N1, B( 0, 0 ), LDB )
811: CALL DGEMM( 'N', 'N', M, N2, N1, -ONE, B( 0, 0 ),
812: $ LDB, A( N1*N1 ), N1, ALPHA, B( 0, N1 ),
813: $ LDB )
814: CALL DTRSM( 'R', 'L', 'T', DIAG, M, N2, ONE,
815: $ A( 1 ), N1, B( 0, N1 ), LDB )
816: *
817: END IF
818: *
819: ELSE
820: *
821: * SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'U'
822: *
823: IF( NOTRANS ) THEN
824: *
825: * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and
826: * TRANS = 'N'
827: *
828: CALL DTRSM( 'R', 'U', 'N', DIAG, M, N1, ALPHA,
829: $ A( N2*N2 ), N2, B( 0, 0 ), LDB )
830: CALL DGEMM( 'N', 'T', M, N2, N1, -ONE, B( 0, 0 ),
831: $ LDB, A( 0 ), N2, ALPHA, B( 0, N1 ),
832: $ LDB )
833: CALL DTRSM( 'R', 'L', 'T', DIAG, M, N2, ONE,
834: $ A( N1*N2 ), N2, B( 0, N1 ), LDB )
835: *
836: ELSE
837: *
838: * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and
839: * TRANS = 'T'
840: *
841: CALL DTRSM( 'R', 'L', 'N', DIAG, M, N2, ALPHA,
842: $ A( N1*N2 ), N2, B( 0, N1 ), LDB )
843: CALL DGEMM( 'N', 'N', M, N1, N2, -ONE, B( 0, N1 ),
844: $ LDB, A( 0 ), N2, ALPHA, B( 0, 0 ),
845: $ LDB )
846: CALL DTRSM( 'R', 'U', 'T', DIAG, M, N1, ONE,
847: $ A( N2*N2 ), N2, B( 0, 0 ), LDB )
848: *
849: END IF
850: *
851: END IF
852: *
853: END IF
854: *
855: ELSE
856: *
857: * SIDE = 'R' and N is even
858: *
859: IF( NORMALTRANSR ) THEN
860: *
861: * SIDE = 'R', N is even, and TRANSR = 'N'
862: *
863: IF( LOWER ) THEN
864: *
865: * SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'L'
866: *
867: IF( NOTRANS ) THEN
868: *
869: * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L',
870: * and TRANS = 'N'
871: *
872: CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ALPHA,
873: $ A( 0 ), N+1, B( 0, K ), LDB )
874: CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, K ),
875: $ LDB, A( K+1 ), N+1, ALPHA, B( 0, 0 ),
876: $ LDB )
877: CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ONE,
878: $ A( 1 ), N+1, B( 0, 0 ), LDB )
879: *
880: ELSE
881: *
882: * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L',
883: * and TRANS = 'T'
884: *
885: CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ALPHA,
886: $ A( 1 ), N+1, B( 0, 0 ), LDB )
887: CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, 0 ),
888: $ LDB, A( K+1 ), N+1, ALPHA, B( 0, K ),
889: $ LDB )
890: CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ONE,
891: $ A( 0 ), N+1, B( 0, K ), LDB )
892: *
893: END IF
894: *
895: ELSE
896: *
897: * SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'U'
898: *
899: IF( NOTRANS ) THEN
900: *
901: * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U',
902: * and TRANS = 'N'
903: *
904: CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ALPHA,
905: $ A( K+1 ), N+1, B( 0, 0 ), LDB )
906: CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, 0 ),
907: $ LDB, A( 0 ), N+1, ALPHA, B( 0, K ),
908: $ LDB )
909: CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ONE,
910: $ A( K ), N+1, B( 0, K ), LDB )
911: *
912: ELSE
913: *
914: * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U',
915: * and TRANS = 'T'
916: *
917: CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ALPHA,
918: $ A( K ), N+1, B( 0, K ), LDB )
919: CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, K ),
920: $ LDB, A( 0 ), N+1, ALPHA, B( 0, 0 ),
921: $ LDB )
922: CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ONE,
923: $ A( K+1 ), N+1, B( 0, 0 ), LDB )
924: *
925: END IF
926: *
927: END IF
928: *
929: ELSE
930: *
931: * SIDE = 'R', N is even, and TRANSR = 'T'
932: *
933: IF( LOWER ) THEN
934: *
935: * SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'L'
936: *
937: IF( NOTRANS ) THEN
938: *
939: * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L',
940: * and TRANS = 'N'
941: *
942: CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ALPHA,
943: $ A( 0 ), K, B( 0, K ), LDB )
944: CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, K ),
945: $ LDB, A( ( K+1 )*K ), K, ALPHA,
946: $ B( 0, 0 ), LDB )
947: CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ONE,
948: $ A( K ), K, B( 0, 0 ), LDB )
949: *
950: ELSE
951: *
952: * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L',
953: * and TRANS = 'T'
954: *
955: CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ALPHA,
956: $ A( K ), K, B( 0, 0 ), LDB )
957: CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, 0 ),
958: $ LDB, A( ( K+1 )*K ), K, ALPHA,
959: $ B( 0, K ), LDB )
960: CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ONE,
961: $ A( 0 ), K, B( 0, K ), LDB )
962: *
963: END IF
964: *
965: ELSE
966: *
967: * SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'U'
968: *
969: IF( NOTRANS ) THEN
970: *
971: * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U',
972: * and TRANS = 'N'
973: *
974: CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ALPHA,
975: $ A( ( K+1 )*K ), K, B( 0, 0 ), LDB )
976: CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, 0 ),
977: $ LDB, A( 0 ), K, ALPHA, B( 0, K ), LDB )
978: CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ONE,
979: $ A( K*K ), K, B( 0, K ), LDB )
980: *
981: ELSE
982: *
983: * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U',
984: * and TRANS = 'T'
985: *
986: CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ALPHA,
987: $ A( K*K ), K, B( 0, K ), LDB )
988: CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, K ),
989: $ LDB, A( 0 ), K, ALPHA, B( 0, 0 ), LDB )
990: CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ONE,
991: $ A( ( K+1 )*K ), K, B( 0, 0 ), LDB )
992: *
993: END IF
994: *
995: END IF
996: *
997: END IF
998: *
999: END IF
1000: END IF
1001: *
1002: RETURN
1003: *
1004: * End of DTFSM
1005: *
1006: 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