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