001:       SUBROUTINE DTFSM( TRANSR, SIDE, UPLO, TRANS, DIAG, M, N, ALPHA, A,
002:      +                  B, LDB )
003: *
004: *  -- LAPACK routine (version 3.2.1)                                    --
005: *
006: *  -- Contributed by Fred Gustavson of the IBM Watson Research Center --
007: *  -- April 2009                                                      --
008: *
009: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
010: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
011: *
012: *     ..
013: *     .. Scalar Arguments ..
014:       CHARACTER          TRANSR, DIAG, SIDE, TRANS, UPLO
015:       INTEGER            LDB, M, N
016:       DOUBLE PRECISION   ALPHA
017: *     ..
018: *     .. Array Arguments ..
019:       DOUBLE PRECISION   A( 0: * ), B( 0: LDB-1, 0: * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  Level 3 BLAS like routine for A in RFP Format.
026: *
027: *  DTFSM  solves the matrix equation
028: *
029: *     op( A )*X = alpha*B  or  X*op( A ) = alpha*B
030: *
031: *  where alpha is a scalar, X and B are m by n matrices, A is a unit, or
032: *  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
033: *
034: *     op( A ) = A   or   op( A ) = A'.
035: *
036: *  A is in Rectangular Full Packed (RFP) Format.
037: *
038: *  The matrix X is overwritten on B.
039: *
040: *  Arguments
041: *  ==========
042: *
043: *  TRANSR - (input) CHARACTER
044: *          = 'N':  The Normal Form of RFP A is stored;
045: *          = 'T':  The Transpose Form of RFP A is stored.
046: *
047: *  SIDE   - (input) CHARACTER
048: *           On entry, SIDE specifies whether op( A ) appears on the left
049: *           or right of X as follows:
050: *
051: *              SIDE = 'L' or 'l'   op( A )*X = alpha*B.
052: *
053: *              SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
054: *
055: *           Unchanged on exit.
056: *
057: *  UPLO   - (input) CHARACTER
058: *           On entry, UPLO specifies whether the RFP matrix A came from
059: *           an upper or lower triangular matrix as follows:
060: *           UPLO = 'U' or 'u' RFP A came from an upper triangular matrix
061: *           UPLO = 'L' or 'l' RFP A came from a  lower triangular matrix
062: *
063: *           Unchanged on exit.
064: *
065: *  TRANS  - (input) CHARACTER
066: *           On entry, TRANS  specifies the form of op( A ) to be used
067: *           in the matrix multiplication as follows:
068: *
069: *              TRANS  = 'N' or 'n'   op( A ) = A.
070: *
071: *              TRANS  = 'T' or 't'   op( A ) = A'.
072: *
073: *           Unchanged on exit.
074: *
075: *  DIAG   - (input) CHARACTER
076: *           On entry, DIAG specifies whether or not RFP A is unit
077: *           triangular as follows:
078: *
079: *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
080: *
081: *              DIAG = 'N' or 'n'   A is not assumed to be unit
082: *                                  triangular.
083: *
084: *           Unchanged on exit.
085: *
086: *  M      - (input) INTEGER.
087: *           On entry, M specifies the number of rows of B. M must be at
088: *           least zero.
089: *           Unchanged on exit.
090: *
091: *  N      - (input) INTEGER.
092: *           On entry, N specifies the number of columns of B.  N must be
093: *           at least zero.
094: *           Unchanged on exit.
095: *
096: *  ALPHA  - (input) DOUBLE PRECISION.
097: *           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
098: *           zero then  A is not referenced and  B need not be set before
099: *           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/ouptut) 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 first 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
929: