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