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