LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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  ELSE IF( m.LT.0 ) THEN
334  info = -6
335  ELSE IF( n.LT.0 ) THEN
336  info = -7
337  ELSE IF( ldb.LT.max( 1, m ) ) THEN
338  info = -11
339  END IF
340  IF( info.NE.0 ) THEN
341  CALL xerbla( 'DTFSM ', -info )
342  RETURN
343  END IF
344 *
345 * Quick return when ( (N.EQ.0).OR.(M.EQ.0) )
346 *
347  IF( ( m.EQ.0 ) .OR. ( n.EQ.0 ) )
348  $ RETURN
349 *
350 * Quick return when ALPHA.EQ.(0D+0)
351 *
352  IF( alpha.EQ.zero ) THEN
353  DO 20 j = 0, n - 1
354  DO 10 i = 0, m - 1
355  b( i, j ) = zero
356  10 CONTINUE
357  20 CONTINUE
358  RETURN
359  END IF
360 *
361  IF( lside ) THEN
362 *
363 * SIDE = 'L'
364 *
365 * A is M-by-M.
366 * If M is odd, set NISODD = .TRUE., and M1 and M2.
367 * If M is even, NISODD = .FALSE., and M.
368 *
369  IF( mod( m, 2 ).EQ.0 ) THEN
370  misodd = .false.
371  k = m / 2
372  ELSE
373  misodd = .true.
374  IF( lower ) THEN
375  m2 = m / 2
376  m1 = m - m2
377  ELSE
378  m1 = m / 2
379  m2 = m - m1
380  END IF
381  END IF
382 *
383 *
384  IF( misodd ) THEN
385 *
386 * SIDE = 'L' and N is odd
387 *
388  IF( normaltransr ) THEN
389 *
390 * SIDE = 'L', N is odd, and TRANSR = 'N'
391 *
392  IF( lower ) THEN
393 *
394 * SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'L'
395 *
396  IF( notrans ) THEN
397 *
398 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
399 * TRANS = 'N'
400 *
401  IF( m.EQ.1 ) THEN
402  CALL dtrsm( 'L', 'L', 'N', diag, m1, n, alpha,
403  $ a, m, b, ldb )
404  ELSE
405  CALL dtrsm( 'L', 'L', 'N', diag, m1, n, alpha,
406  $ a( 0 ), m, b, ldb )
407  CALL dgemm( 'N', 'N', m2, n, m1, -one, a( m1 ),
408  $ m, b, ldb, alpha, b( m1, 0 ), ldb )
409  CALL dtrsm( 'L', 'U', 'T', diag, m2, n, one,
410  $ a( m ), m, b( m1, 0 ), ldb )
411  END IF
412 *
413  ELSE
414 *
415 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
416 * TRANS = 'T'
417 *
418  IF( m.EQ.1 ) THEN
419  CALL dtrsm( 'L', 'L', 'T', diag, m1, n, alpha,
420  $ a( 0 ), m, b, ldb )
421  ELSE
422  CALL dtrsm( 'L', 'U', 'N', diag, m2, n, alpha,
423  $ a( m ), m, b( m1, 0 ), ldb )
424  CALL dgemm( 'T', 'N', m1, n, m2, -one, a( m1 ),
425  $ m, b( m1, 0 ), ldb, alpha, b, ldb )
426  CALL dtrsm( 'L', 'L', 'T', diag, m1, n, one,
427  $ a( 0 ), m, b, ldb )
428  END IF
429 *
430  END IF
431 *
432  ELSE
433 *
434 * SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'U'
435 *
436  IF( .NOT.notrans ) THEN
437 *
438 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
439 * TRANS = 'N'
440 *
441  CALL dtrsm( 'L', 'L', 'N', diag, m1, n, alpha,
442  $ a( m2 ), m, b, ldb )
443  CALL dgemm( 'T', 'N', m2, n, m1, -one, a( 0 ), m,
444  $ b, ldb, alpha, b( m1, 0 ), ldb )
445  CALL dtrsm( 'L', 'U', 'T', diag, m2, n, one,
446  $ a( m1 ), m, b( m1, 0 ), ldb )
447 *
448  ELSE
449 *
450 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
451 * TRANS = 'T'
452 *
453  CALL dtrsm( 'L', 'U', 'N', diag, m2, n, alpha,
454  $ a( m1 ), m, b( m1, 0 ), ldb )
455  CALL dgemm( 'N', 'N', m1, n, m2, -one, a( 0 ), m,
456  $ b( m1, 0 ), ldb, alpha, b, ldb )
457  CALL dtrsm( 'L', 'L', 'T', diag, m1, n, one,
458  $ a( m2 ), m, b, ldb )
459 *
460  END IF
461 *
462  END IF
463 *
464  ELSE
465 *
466 * SIDE = 'L', N is odd, and TRANSR = 'T'
467 *
468  IF( lower ) THEN
469 *
470 * SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'L'
471 *
472  IF( notrans ) THEN
473 *
474 * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and
475 * TRANS = 'N'
476 *
477  IF( m.EQ.1 ) THEN
478  CALL dtrsm( 'L', 'U', 'T', diag, m1, n, alpha,
479  $ a( 0 ), m1, b, ldb )
480  ELSE
481  CALL dtrsm( 'L', 'U', 'T', diag, m1, n, alpha,
482  $ a( 0 ), m1, b, ldb )
483  CALL dgemm( 'T', 'N', m2, n, m1, -one,
484  $ a( m1*m1 ), m1, b, ldb, alpha,
485  $ b( m1, 0 ), ldb )
486  CALL dtrsm( 'L', 'L', 'N', diag, m2, n, one,
487  $ a( 1 ), m1, b( m1, 0 ), ldb )
488  END IF
489 *
490  ELSE
491 *
492 * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and
493 * TRANS = 'T'
494 *
495  IF( m.EQ.1 ) THEN
496  CALL dtrsm( 'L', 'U', 'N', diag, m1, n, alpha,
497  $ a( 0 ), m1, b, ldb )
498  ELSE
499  CALL dtrsm( 'L', 'L', 'T', diag, m2, n, alpha,
500  $ a( 1 ), m1, b( m1, 0 ), ldb )
501  CALL dgemm( 'N', 'N', m1, n, m2, -one,
502  $ a( m1*m1 ), m1, b( m1, 0 ), ldb,
503  $ alpha, b, ldb )
504  CALL dtrsm( 'L', 'U', 'N', diag, m1, n, one,
505  $ a( 0 ), m1, b, ldb )
506  END IF
507 *
508  END IF
509 *
510  ELSE
511 *
512 * SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'U'
513 *
514  IF( .NOT.notrans ) THEN
515 *
516 * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and
517 * TRANS = 'N'
518 *
519  CALL dtrsm( 'L', 'U', 'T', diag, m1, n, alpha,
520  $ a( m2*m2 ), m2, b, ldb )
521  CALL dgemm( 'N', 'N', m2, n, m1, -one, a( 0 ), m2,
522  $ b, ldb, alpha, b( m1, 0 ), ldb )
523  CALL dtrsm( 'L', 'L', 'N', diag, m2, n, one,
524  $ a( m1*m2 ), m2, b( m1, 0 ), ldb )
525 *
526  ELSE
527 *
528 * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and
529 * TRANS = 'T'
530 *
531  CALL dtrsm( 'L', 'L', 'T', diag, m2, n, alpha,
532  $ a( m1*m2 ), m2, b( m1, 0 ), ldb )
533  CALL dgemm( 'T', 'N', m1, n, m2, -one, a( 0 ), m2,
534  $ b( m1, 0 ), ldb, alpha, b, ldb )
535  CALL dtrsm( 'L', 'U', 'N', diag, m1, n, one,
536  $ a( m2*m2 ), m2, b, ldb )
537 *
538  END IF
539 *
540  END IF
541 *
542  END IF
543 *
544  ELSE
545 *
546 * SIDE = 'L' and N is even
547 *
548  IF( normaltransr ) THEN
549 *
550 * SIDE = 'L', N is even, and TRANSR = 'N'
551 *
552  IF( lower ) THEN
553 *
554 * SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'L'
555 *
556  IF( notrans ) THEN
557 *
558 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L',
559 * and TRANS = 'N'
560 *
561  CALL dtrsm( 'L', 'L', 'N', diag, k, n, alpha,
562  $ a( 1 ), m+1, b, ldb )
563  CALL dgemm( 'N', 'N', k, n, k, -one, a( k+1 ),
564  $ m+1, b, ldb, alpha, b( k, 0 ), ldb )
565  CALL dtrsm( 'L', 'U', 'T', diag, k, n, one,
566  $ a( 0 ), m+1, b( k, 0 ), ldb )
567 *
568  ELSE
569 *
570 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L',
571 * and TRANS = 'T'
572 *
573  CALL dtrsm( 'L', 'U', 'N', diag, k, n, alpha,
574  $ a( 0 ), m+1, b( k, 0 ), ldb )
575  CALL dgemm( 'T', 'N', k, n, k, -one, a( k+1 ),
576  $ m+1, b( k, 0 ), ldb, alpha, b, ldb )
577  CALL dtrsm( 'L', 'L', 'T', diag, k, n, one,
578  $ a( 1 ), m+1, b, ldb )
579 *
580  END IF
581 *
582  ELSE
583 *
584 * SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'U'
585 *
586  IF( .NOT.notrans ) THEN
587 *
588 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U',
589 * and TRANS = 'N'
590 *
591  CALL dtrsm( 'L', 'L', 'N', diag, k, n, alpha,
592  $ a( k+1 ), m+1, b, ldb )
593  CALL dgemm( 'T', 'N', k, n, k, -one, a( 0 ), m+1,
594  $ b, ldb, alpha, b( k, 0 ), ldb )
595  CALL dtrsm( 'L', 'U', 'T', diag, k, n, one,
596  $ a( k ), m+1, b( k, 0 ), ldb )
597 *
598  ELSE
599 *
600 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U',
601 * and TRANS = 'T'
602  CALL dtrsm( 'L', 'U', 'N', diag, k, n, alpha,
603  $ a( k ), m+1, b( k, 0 ), ldb )
604  CALL dgemm( 'N', 'N', k, n, k, -one, a( 0 ), m+1,
605  $ b( k, 0 ), ldb, alpha, b, ldb )
606  CALL dtrsm( 'L', 'L', 'T', diag, k, n, one,
607  $ a( k+1 ), m+1, b, ldb )
608 *
609  END IF
610 *
611  END IF
612 *
613  ELSE
614 *
615 * SIDE = 'L', N is even, and TRANSR = 'T'
616 *
617  IF( lower ) THEN
618 *
619 * SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'L'
620 *
621  IF( notrans ) THEN
622 *
623 * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L',
624 * and TRANS = 'N'
625 *
626  CALL dtrsm( 'L', 'U', 'T', diag, k, n, alpha,
627  $ a( k ), k, b, ldb )
628  CALL dgemm( 'T', 'N', k, n, k, -one,
629  $ a( k*( k+1 ) ), k, b, ldb, alpha,
630  $ b( k, 0 ), ldb )
631  CALL dtrsm( 'L', 'L', 'N', diag, k, n, one,
632  $ a( 0 ), k, b( k, 0 ), ldb )
633 *
634  ELSE
635 *
636 * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L',
637 * and TRANS = 'T'
638 *
639  CALL dtrsm( 'L', 'L', 'T', diag, k, n, alpha,
640  $ a( 0 ), k, b( k, 0 ), ldb )
641  CALL dgemm( 'N', 'N', k, n, k, -one,
642  $ a( k*( k+1 ) ), k, b( k, 0 ), ldb,
643  $ alpha, b, ldb )
644  CALL dtrsm( 'L', 'U', 'N', diag, k, n, one,
645  $ a( k ), k, b, ldb )
646 *
647  END IF
648 *
649  ELSE
650 *
651 * SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'U'
652 *
653  IF( .NOT.notrans ) THEN
654 *
655 * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U',
656 * and TRANS = 'N'
657 *
658  CALL dtrsm( 'L', 'U', 'T', diag, k, n, alpha,
659  $ a( k*( k+1 ) ), k, b, ldb )
660  CALL dgemm( 'N', 'N', k, n, k, -one, a( 0 ), k, b,
661  $ ldb, alpha, b( k, 0 ), ldb )
662  CALL dtrsm( 'L', 'L', 'N', diag, k, n, one,
663  $ a( k*k ), k, b( k, 0 ), ldb )
664 *
665  ELSE
666 *
667 * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U',
668 * and TRANS = 'T'
669 *
670  CALL dtrsm( 'L', 'L', 'T', diag, k, n, alpha,
671  $ a( k*k ), k, b( k, 0 ), ldb )
672  CALL dgemm( 'T', 'N', k, n, k, -one, a( 0 ), k,
673  $ b( k, 0 ), ldb, alpha, b, ldb )
674  CALL dtrsm( 'L', 'U', 'N', diag, k, n, one,
675  $ a( k*( k+1 ) ), k, b, ldb )
676 *
677  END IF
678 *
679  END IF
680 *
681  END IF
682 *
683  END IF
684 *
685  ELSE
686 *
687 * SIDE = 'R'
688 *
689 * A is N-by-N.
690 * If N is odd, set NISODD = .TRUE., and N1 and N2.
691 * If N is even, NISODD = .FALSE., and K.
692 *
693  IF( mod( n, 2 ).EQ.0 ) THEN
694  nisodd = .false.
695  k = n / 2
696  ELSE
697  nisodd = .true.
698  IF( lower ) THEN
699  n2 = n / 2
700  n1 = n - n2
701  ELSE
702  n1 = n / 2
703  n2 = n - n1
704  END IF
705  END IF
706 *
707  IF( nisodd ) THEN
708 *
709 * SIDE = 'R' and N is odd
710 *
711  IF( normaltransr ) THEN
712 *
713 * SIDE = 'R', N is odd, and TRANSR = 'N'
714 *
715  IF( lower ) THEN
716 *
717 * SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'L'
718 *
719  IF( notrans ) THEN
720 *
721 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
722 * TRANS = 'N'
723 *
724  CALL dtrsm( 'R', 'U', 'T', diag, m, n2, alpha,
725  $ a( n ), n, b( 0, n1 ), ldb )
726  CALL dgemm( 'N', 'N', m, n1, n2, -one, b( 0, n1 ),
727  $ ldb, a( n1 ), n, alpha, b( 0, 0 ),
728  $ ldb )
729  CALL dtrsm( 'R', 'L', 'N', diag, m, n1, one,
730  $ a( 0 ), n, b( 0, 0 ), ldb )
731 *
732  ELSE
733 *
734 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
735 * TRANS = 'T'
736 *
737  CALL dtrsm( 'R', 'L', 'T', diag, m, n1, alpha,
738  $ a( 0 ), n, b( 0, 0 ), ldb )
739  CALL dgemm( 'N', 'T', m, n2, n1, -one, b( 0, 0 ),
740  $ ldb, a( n1 ), n, alpha, b( 0, n1 ),
741  $ ldb )
742  CALL dtrsm( 'R', 'U', 'N', diag, m, n2, one,
743  $ a( n ), n, b( 0, n1 ), ldb )
744 *
745  END IF
746 *
747  ELSE
748 *
749 * SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'U'
750 *
751  IF( notrans ) THEN
752 *
753 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
754 * TRANS = 'N'
755 *
756  CALL dtrsm( 'R', 'L', 'T', diag, m, n1, alpha,
757  $ a( n2 ), n, b( 0, 0 ), ldb )
758  CALL dgemm( 'N', 'N', m, n2, n1, -one, b( 0, 0 ),
759  $ ldb, a( 0 ), n, alpha, b( 0, n1 ),
760  $ ldb )
761  CALL dtrsm( 'R', 'U', 'N', diag, m, n2, one,
762  $ a( n1 ), n, b( 0, n1 ), ldb )
763 *
764  ELSE
765 *
766 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
767 * TRANS = 'T'
768 *
769  CALL dtrsm( 'R', 'U', 'T', diag, m, n2, alpha,
770  $ a( n1 ), n, b( 0, n1 ), ldb )
771  CALL dgemm( 'N', 'T', m, n1, n2, -one, b( 0, n1 ),
772  $ ldb, a( 0 ), n, alpha, b( 0, 0 ), ldb )
773  CALL dtrsm( 'R', 'L', 'N', diag, m, n1, one,
774  $ a( n2 ), n, b( 0, 0 ), ldb )
775 *
776  END IF
777 *
778  END IF
779 *
780  ELSE
781 *
782 * SIDE = 'R', N is odd, and TRANSR = 'T'
783 *
784  IF( lower ) THEN
785 *
786 * SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'L'
787 *
788  IF( notrans ) THEN
789 *
790 * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and
791 * TRANS = 'N'
792 *
793  CALL dtrsm( 'R', 'L', 'N', diag, m, n2, alpha,
794  $ a( 1 ), n1, b( 0, n1 ), ldb )
795  CALL dgemm( 'N', 'T', m, n1, n2, -one, b( 0, n1 ),
796  $ ldb, a( n1*n1 ), n1, alpha, b( 0, 0 ),
797  $ ldb )
798  CALL dtrsm( 'R', 'U', 'T', diag, m, n1, one,
799  $ a( 0 ), n1, b( 0, 0 ), ldb )
800 *
801  ELSE
802 *
803 * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and
804 * TRANS = 'T'
805 *
806  CALL dtrsm( 'R', 'U', 'N', diag, m, n1, alpha,
807  $ a( 0 ), n1, b( 0, 0 ), ldb )
808  CALL dgemm( 'N', 'N', m, n2, n1, -one, b( 0, 0 ),
809  $ ldb, a( n1*n1 ), n1, alpha, b( 0, n1 ),
810  $ ldb )
811  CALL dtrsm( 'R', 'L', 'T', diag, m, n2, one,
812  $ a( 1 ), n1, b( 0, n1 ), ldb )
813 *
814  END IF
815 *
816  ELSE
817 *
818 * SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'U'
819 *
820  IF( notrans ) THEN
821 *
822 * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and
823 * TRANS = 'N'
824 *
825  CALL dtrsm( 'R', 'U', 'N', diag, m, n1, alpha,
826  $ a( n2*n2 ), n2, b( 0, 0 ), ldb )
827  CALL dgemm( 'N', 'T', m, n2, n1, -one, b( 0, 0 ),
828  $ ldb, a( 0 ), n2, alpha, b( 0, n1 ),
829  $ ldb )
830  CALL dtrsm( 'R', 'L', 'T', diag, m, n2, one,
831  $ a( n1*n2 ), n2, b( 0, n1 ), ldb )
832 *
833  ELSE
834 *
835 * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and
836 * TRANS = 'T'
837 *
838  CALL dtrsm( 'R', 'L', 'N', diag, m, n2, alpha,
839  $ a( n1*n2 ), n2, b( 0, n1 ), ldb )
840  CALL dgemm( 'N', 'N', m, n1, n2, -one, b( 0, n1 ),
841  $ ldb, a( 0 ), n2, alpha, b( 0, 0 ),
842  $ ldb )
843  CALL dtrsm( 'R', 'U', 'T', diag, m, n1, one,
844  $ a( n2*n2 ), n2, b( 0, 0 ), ldb )
845 *
846  END IF
847 *
848  END IF
849 *
850  END IF
851 *
852  ELSE
853 *
854 * SIDE = 'R' and N is even
855 *
856  IF( normaltransr ) THEN
857 *
858 * SIDE = 'R', N is even, and TRANSR = 'N'
859 *
860  IF( lower ) THEN
861 *
862 * SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'L'
863 *
864  IF( notrans ) THEN
865 *
866 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L',
867 * and TRANS = 'N'
868 *
869  CALL dtrsm( 'R', 'U', 'T', diag, m, k, alpha,
870  $ a( 0 ), n+1, b( 0, k ), ldb )
871  CALL dgemm( 'N', 'N', m, k, k, -one, b( 0, k ),
872  $ ldb, a( k+1 ), n+1, alpha, b( 0, 0 ),
873  $ ldb )
874  CALL dtrsm( 'R', 'L', 'N', diag, m, k, one,
875  $ a( 1 ), n+1, b( 0, 0 ), ldb )
876 *
877  ELSE
878 *
879 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L',
880 * and TRANS = 'T'
881 *
882  CALL dtrsm( 'R', 'L', 'T', diag, m, k, alpha,
883  $ a( 1 ), n+1, b( 0, 0 ), ldb )
884  CALL dgemm( 'N', 'T', m, k, k, -one, b( 0, 0 ),
885  $ ldb, a( k+1 ), n+1, alpha, b( 0, k ),
886  $ ldb )
887  CALL dtrsm( 'R', 'U', 'N', diag, m, k, one,
888  $ a( 0 ), n+1, b( 0, k ), ldb )
889 *
890  END IF
891 *
892  ELSE
893 *
894 * SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'U'
895 *
896  IF( notrans ) THEN
897 *
898 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U',
899 * and TRANS = 'N'
900 *
901  CALL dtrsm( 'R', 'L', 'T', diag, m, k, alpha,
902  $ a( k+1 ), n+1, b( 0, 0 ), ldb )
903  CALL dgemm( 'N', 'N', m, k, k, -one, b( 0, 0 ),
904  $ ldb, a( 0 ), n+1, alpha, b( 0, k ),
905  $ ldb )
906  CALL dtrsm( 'R', 'U', 'N', diag, m, k, one,
907  $ a( k ), n+1, b( 0, k ), ldb )
908 *
909  ELSE
910 *
911 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U',
912 * and TRANS = 'T'
913 *
914  CALL dtrsm( 'R', 'U', 'T', diag, m, k, alpha,
915  $ a( k ), n+1, b( 0, k ), ldb )
916  CALL dgemm( 'N', 'T', m, k, k, -one, b( 0, k ),
917  $ ldb, a( 0 ), n+1, alpha, b( 0, 0 ),
918  $ ldb )
919  CALL dtrsm( 'R', 'L', 'N', diag, m, k, one,
920  $ a( k+1 ), n+1, b( 0, 0 ), ldb )
921 *
922  END IF
923 *
924  END IF
925 *
926  ELSE
927 *
928 * SIDE = 'R', N is even, and TRANSR = 'T'
929 *
930  IF( lower ) THEN
931 *
932 * SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'L'
933 *
934  IF( notrans ) THEN
935 *
936 * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L',
937 * and TRANS = 'N'
938 *
939  CALL dtrsm( 'R', 'L', 'N', diag, m, k, alpha,
940  $ a( 0 ), k, b( 0, k ), ldb )
941  CALL dgemm( 'N', 'T', m, k, k, -one, b( 0, k ),
942  $ ldb, a( ( k+1 )*k ), k, alpha,
943  $ b( 0, 0 ), ldb )
944  CALL dtrsm( 'R', 'U', 'T', diag, m, k, one,
945  $ a( k ), k, b( 0, 0 ), ldb )
946 *
947  ELSE
948 *
949 * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L',
950 * and TRANS = 'T'
951 *
952  CALL dtrsm( 'R', 'U', 'N', diag, m, k, alpha,
953  $ a( k ), k, b( 0, 0 ), ldb )
954  CALL dgemm( 'N', 'N', m, k, k, -one, b( 0, 0 ),
955  $ ldb, a( ( k+1 )*k ), k, alpha,
956  $ b( 0, k ), ldb )
957  CALL dtrsm( 'R', 'L', 'T', diag, m, k, one,
958  $ a( 0 ), k, b( 0, k ), ldb )
959 *
960  END IF
961 *
962  ELSE
963 *
964 * SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'U'
965 *
966  IF( notrans ) THEN
967 *
968 * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U',
969 * and TRANS = 'N'
970 *
971  CALL dtrsm( 'R', 'U', 'N', diag, m, k, alpha,
972  $ a( ( k+1 )*k ), k, b( 0, 0 ), ldb )
973  CALL dgemm( 'N', 'T', m, k, k, -one, b( 0, 0 ),
974  $ ldb, a( 0 ), k, alpha, b( 0, k ), ldb )
975  CALL dtrsm( 'R', 'L', 'T', diag, m, k, one,
976  $ a( k*k ), k, b( 0, k ), ldb )
977 *
978  ELSE
979 *
980 * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U',
981 * and TRANS = 'T'
982 *
983  CALL dtrsm( 'R', 'L', 'N', diag, m, k, alpha,
984  $ a( k*k ), k, b( 0, k ), ldb )
985  CALL dgemm( 'N', 'N', m, k, k, -one, b( 0, k ),
986  $ ldb, a( 0 ), k, alpha, b( 0, 0 ), ldb )
987  CALL dtrsm( 'R', 'U', 'T', diag, m, k, one,
988  $ a( ( k+1 )*k ), k, b( 0, 0 ), ldb )
989 *
990  END IF
991 *
992  END IF
993 *
994  END IF
995 *
996  END IF
997  END IF
998 *
999  RETURN
1000 *
1001 * End of DTFSM
1002 *
1003  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
Definition: dtrsm.f:181
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
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