LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
stfsm.f
Go to the documentation of this file.
1 *> \brief \b STFSM 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 STFSM + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stfsm.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stfsm.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stfsm.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE STFSM( 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 * REAL ALPHA
28 * ..
29 * .. Array Arguments ..
30 * REAL 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 *> STFSM 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 REAL
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 REAL 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 REAL 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 June 2017
186 *
187 *> \ingroup realOTHERcomputational
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 stfsm( TRANSR, SIDE, UPLO, TRANS, DIAG, M, N, ALPHA, A,
278  $ B, LDB )
279 *
280 * -- LAPACK computational routine (version 3.7.1) --
281 * -- LAPACK is a software package provided by Univ. of Tennessee, --
282 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
283 * June 2017
284 *
285 * .. Scalar Arguments ..
286  CHARACTER TRANSR, DIAG, SIDE, TRANS, UPLO
287  INTEGER LDB, M, N
288  REAL ALPHA
289 * ..
290 * .. Array Arguments ..
291  REAL A( 0: * ), B( 0: LDB-1, 0: * )
292 * ..
293 *
294 * =====================================================================
295 *
296 * ..
297 * .. Parameters ..
298  REAL ONE, ZERO
299  parameter( one = 1.0e+0, zero = 0.0e+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 sgemm, strsm, xerbla
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( 'STFSM ', -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  IF( misodd ) THEN
387 *
388 * SIDE = 'L' and N is odd
389 *
390  IF( normaltransr ) THEN
391 *
392 * SIDE = 'L', N is odd, and TRANSR = 'N'
393 *
394  IF( lower ) THEN
395 *
396 * SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'L'
397 *
398  IF( notrans ) THEN
399 *
400 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
401 * TRANS = 'N'
402 *
403  IF( m.EQ.1 ) THEN
404  CALL strsm( 'L', 'L', 'N', diag, m1, n, alpha,
405  $ a, m, b, ldb )
406  ELSE
407  CALL strsm( 'L', 'L', 'N', diag, m1, n, alpha,
408  $ a( 0 ), m, b, ldb )
409  CALL sgemm( 'N', 'N', m2, n, m1, -one, a( m1 ),
410  $ m, b, ldb, alpha, b( m1, 0 ), ldb )
411  CALL strsm( 'L', 'U', 'T', diag, m2, n, one,
412  $ a( m ), m, b( m1, 0 ), ldb )
413  END IF
414 *
415  ELSE
416 *
417 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
418 * TRANS = 'T'
419 *
420  IF( m.EQ.1 ) THEN
421  CALL strsm( 'L', 'L', 'T', diag, m1, n, alpha,
422  $ a( 0 ), m, b, ldb )
423  ELSE
424  CALL strsm( 'L', 'U', 'N', diag, m2, n, alpha,
425  $ a( m ), m, b( m1, 0 ), ldb )
426  CALL sgemm( 'T', 'N', m1, n, m2, -one, a( m1 ),
427  $ m, b( m1, 0 ), ldb, alpha, b, ldb )
428  CALL strsm( 'L', 'L', 'T', diag, m1, n, one,
429  $ a( 0 ), m, b, ldb )
430  END IF
431 *
432  END IF
433 *
434  ELSE
435 *
436 * SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'U'
437 *
438  IF( .NOT.notrans ) THEN
439 *
440 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
441 * TRANS = 'N'
442 *
443  CALL strsm( 'L', 'L', 'N', diag, m1, n, alpha,
444  $ a( m2 ), m, b, ldb )
445  CALL sgemm( 'T', 'N', m2, n, m1, -one, a( 0 ), m,
446  $ b, ldb, alpha, b( m1, 0 ), ldb )
447  CALL strsm( 'L', 'U', 'T', diag, m2, n, one,
448  $ a( m1 ), m, b( m1, 0 ), ldb )
449 *
450  ELSE
451 *
452 * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
453 * TRANS = 'T'
454 *
455  CALL strsm( 'L', 'U', 'N', diag, m2, n, alpha,
456  $ a( m1 ), m, b( m1, 0 ), ldb )
457  CALL sgemm( 'N', 'N', m1, n, m2, -one, a( 0 ), m,
458  $ b( m1, 0 ), ldb, alpha, b, ldb )
459  CALL strsm( 'L', 'L', 'T', diag, m1, n, one,
460  $ a( m2 ), m, b, ldb )
461 *
462  END IF
463 *
464  END IF
465 *
466  ELSE
467 *
468 * SIDE = 'L', N is odd, and TRANSR = 'T'
469 *
470  IF( lower ) THEN
471 *
472 * SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'L'
473 *
474  IF( notrans ) THEN
475 *
476 * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and
477 * TRANS = 'N'
478 *
479  IF( m.EQ.1 ) THEN
480  CALL strsm( 'L', 'U', 'T', diag, m1, n, alpha,
481  $ a( 0 ), m1, b, ldb )
482  ELSE
483  CALL strsm( 'L', 'U', 'T', diag, m1, n, alpha,
484  $ a( 0 ), m1, b, ldb )
485  CALL sgemm( 'T', 'N', m2, n, m1, -one,
486  $ a( m1*m1 ), m1, b, ldb, alpha,
487  $ b( m1, 0 ), ldb )
488  CALL strsm( 'L', 'L', 'N', diag, m2, n, one,
489  $ a( 1 ), m1, b( m1, 0 ), ldb )
490  END IF
491 *
492  ELSE
493 *
494 * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and
495 * TRANS = 'T'
496 *
497  IF( m.EQ.1 ) THEN
498  CALL strsm( 'L', 'U', 'N', diag, m1, n, alpha,
499  $ a( 0 ), m1, b, ldb )
500  ELSE
501  CALL strsm( 'L', 'L', 'T', diag, m2, n, alpha,
502  $ a( 1 ), m1, b( m1, 0 ), ldb )
503  CALL sgemm( 'N', 'N', m1, n, m2, -one,
504  $ a( m1*m1 ), m1, b( m1, 0 ), ldb,
505  $ alpha, b, ldb )
506  CALL strsm( 'L', 'U', 'N', diag, m1, n, one,
507  $ a( 0 ), m1, b, ldb )
508  END IF
509 *
510  END IF
511 *
512  ELSE
513 *
514 * SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'U'
515 *
516  IF( .NOT.notrans ) THEN
517 *
518 * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and
519 * TRANS = 'N'
520 *
521  CALL strsm( 'L', 'U', 'T', diag, m1, n, alpha,
522  $ a( m2*m2 ), m2, b, ldb )
523  CALL sgemm( 'N', 'N', m2, n, m1, -one, a( 0 ), m2,
524  $ b, ldb, alpha, b( m1, 0 ), ldb )
525  CALL strsm( 'L', 'L', 'N', diag, m2, n, one,
526  $ a( m1*m2 ), m2, b( m1, 0 ), ldb )
527 *
528  ELSE
529 *
530 * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and
531 * TRANS = 'T'
532 *
533  CALL strsm( 'L', 'L', 'T', diag, m2, n, alpha,
534  $ a( m1*m2 ), m2, b( m1, 0 ), ldb )
535  CALL sgemm( 'T', 'N', m1, n, m2, -one, a( 0 ), m2,
536  $ b( m1, 0 ), ldb, alpha, b, ldb )
537  CALL strsm( 'L', 'U', 'N', diag, m1, n, one,
538  $ a( m2*m2 ), m2, b, ldb )
539 *
540  END IF
541 *
542  END IF
543 *
544  END IF
545 *
546  ELSE
547 *
548 * SIDE = 'L' and N is even
549 *
550  IF( normaltransr ) THEN
551 *
552 * SIDE = 'L', N is even, and TRANSR = 'N'
553 *
554  IF( lower ) THEN
555 *
556 * SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'L'
557 *
558  IF( notrans ) THEN
559 *
560 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L',
561 * and TRANS = 'N'
562 *
563  CALL strsm( 'L', 'L', 'N', diag, k, n, alpha,
564  $ a( 1 ), m+1, b, ldb )
565  CALL sgemm( 'N', 'N', k, n, k, -one, a( k+1 ),
566  $ m+1, b, ldb, alpha, b( k, 0 ), ldb )
567  CALL strsm( 'L', 'U', 'T', diag, k, n, one,
568  $ a( 0 ), m+1, b( k, 0 ), ldb )
569 *
570  ELSE
571 *
572 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L',
573 * and TRANS = 'T'
574 *
575  CALL strsm( 'L', 'U', 'N', diag, k, n, alpha,
576  $ a( 0 ), m+1, b( k, 0 ), ldb )
577  CALL sgemm( 'T', 'N', k, n, k, -one, a( k+1 ),
578  $ m+1, b( k, 0 ), ldb, alpha, b, ldb )
579  CALL strsm( 'L', 'L', 'T', diag, k, n, one,
580  $ a( 1 ), m+1, b, ldb )
581 *
582  END IF
583 *
584  ELSE
585 *
586 * SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'U'
587 *
588  IF( .NOT.notrans ) THEN
589 *
590 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U',
591 * and TRANS = 'N'
592 *
593  CALL strsm( 'L', 'L', 'N', diag, k, n, alpha,
594  $ a( k+1 ), m+1, b, ldb )
595  CALL sgemm( 'T', 'N', k, n, k, -one, a( 0 ), m+1,
596  $ b, ldb, alpha, b( k, 0 ), ldb )
597  CALL strsm( 'L', 'U', 'T', diag, k, n, one,
598  $ a( k ), m+1, b( k, 0 ), ldb )
599 *
600  ELSE
601 *
602 * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U',
603 * and TRANS = 'T'
604  CALL strsm( 'L', 'U', 'N', diag, k, n, alpha,
605  $ a( k ), m+1, b( k, 0 ), ldb )
606  CALL sgemm( 'N', 'N', k, n, k, -one, a( 0 ), m+1,
607  $ b( k, 0 ), ldb, alpha, b, ldb )
608  CALL strsm( 'L', 'L', 'T', diag, k, n, one,
609  $ a( k+1 ), m+1, b, ldb )
610 *
611  END IF
612 *
613  END IF
614 *
615  ELSE
616 *
617 * SIDE = 'L', N is even, and TRANSR = 'T'
618 *
619  IF( lower ) THEN
620 *
621 * SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'L'
622 *
623  IF( notrans ) THEN
624 *
625 * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L',
626 * and TRANS = 'N'
627 *
628  CALL strsm( 'L', 'U', 'T', diag, k, n, alpha,
629  $ a( k ), k, b, ldb )
630  CALL sgemm( 'T', 'N', k, n, k, -one,
631  $ a( k*( k+1 ) ), k, b, ldb, alpha,
632  $ b( k, 0 ), ldb )
633  CALL strsm( 'L', 'L', 'N', diag, k, n, one,
634  $ a( 0 ), k, b( k, 0 ), ldb )
635 *
636  ELSE
637 *
638 * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L',
639 * and TRANS = 'T'
640 *
641  CALL strsm( 'L', 'L', 'T', diag, k, n, alpha,
642  $ a( 0 ), k, b( k, 0 ), ldb )
643  CALL sgemm( 'N', 'N', k, n, k, -one,
644  $ a( k*( k+1 ) ), k, b( k, 0 ), ldb,
645  $ alpha, b, ldb )
646  CALL strsm( 'L', 'U', 'N', diag, k, n, one,
647  $ a( k ), k, b, ldb )
648 *
649  END IF
650 *
651  ELSE
652 *
653 * SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'U'
654 *
655  IF( .NOT.notrans ) THEN
656 *
657 * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U',
658 * and TRANS = 'N'
659 *
660  CALL strsm( 'L', 'U', 'T', diag, k, n, alpha,
661  $ a( k*( k+1 ) ), k, b, ldb )
662  CALL sgemm( 'N', 'N', k, n, k, -one, a( 0 ), k, b,
663  $ ldb, alpha, b( k, 0 ), ldb )
664  CALL strsm( 'L', 'L', 'N', diag, k, n, one,
665  $ a( k*k ), k, b( k, 0 ), ldb )
666 *
667  ELSE
668 *
669 * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U',
670 * and TRANS = 'T'
671 *
672  CALL strsm( 'L', 'L', 'T', diag, k, n, alpha,
673  $ a( k*k ), k, b( k, 0 ), ldb )
674  CALL sgemm( 'T', 'N', k, n, k, -one, a( 0 ), k,
675  $ b( k, 0 ), ldb, alpha, b, ldb )
676  CALL strsm( 'L', 'U', 'N', diag, k, n, one,
677  $ a( k*( k+1 ) ), k, b, ldb )
678 *
679  END IF
680 *
681  END IF
682 *
683  END IF
684 *
685  END IF
686 *
687  ELSE
688 *
689 * SIDE = 'R'
690 *
691 * A is N-by-N.
692 * If N is odd, set NISODD = .TRUE., and N1 and N2.
693 * If N is even, NISODD = .FALSE., and K.
694 *
695  IF( mod( n, 2 ).EQ.0 ) THEN
696  nisodd = .false.
697  k = n / 2
698  ELSE
699  nisodd = .true.
700  IF( lower ) THEN
701  n2 = n / 2
702  n1 = n - n2
703  ELSE
704  n1 = n / 2
705  n2 = n - n1
706  END IF
707  END IF
708 *
709  IF( nisodd ) THEN
710 *
711 * SIDE = 'R' and N is odd
712 *
713  IF( normaltransr ) THEN
714 *
715 * SIDE = 'R', N is odd, and TRANSR = 'N'
716 *
717  IF( lower ) THEN
718 *
719 * SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'L'
720 *
721  IF( notrans ) THEN
722 *
723 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
724 * TRANS = 'N'
725 *
726  CALL strsm( 'R', 'U', 'T', diag, m, n2, alpha,
727  $ a( n ), n, b( 0, n1 ), ldb )
728  CALL sgemm( 'N', 'N', m, n1, n2, -one, b( 0, n1 ),
729  $ ldb, a( n1 ), n, alpha, b( 0, 0 ),
730  $ ldb )
731  CALL strsm( 'R', 'L', 'N', diag, m, n1, one,
732  $ a( 0 ), n, b( 0, 0 ), ldb )
733 *
734  ELSE
735 *
736 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
737 * TRANS = 'T'
738 *
739  CALL strsm( 'R', 'L', 'T', diag, m, n1, alpha,
740  $ a( 0 ), n, b( 0, 0 ), ldb )
741  CALL sgemm( 'N', 'T', m, n2, n1, -one, b( 0, 0 ),
742  $ ldb, a( n1 ), n, alpha, b( 0, n1 ),
743  $ ldb )
744  CALL strsm( 'R', 'U', 'N', diag, m, n2, one,
745  $ a( n ), n, b( 0, n1 ), ldb )
746 *
747  END IF
748 *
749  ELSE
750 *
751 * SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'U'
752 *
753  IF( notrans ) THEN
754 *
755 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
756 * TRANS = 'N'
757 *
758  CALL strsm( 'R', 'L', 'T', diag, m, n1, alpha,
759  $ a( n2 ), n, b( 0, 0 ), ldb )
760  CALL sgemm( 'N', 'N', m, n2, n1, -one, b( 0, 0 ),
761  $ ldb, a( 0 ), n, alpha, b( 0, n1 ),
762  $ ldb )
763  CALL strsm( 'R', 'U', 'N', diag, m, n2, one,
764  $ a( n1 ), n, b( 0, n1 ), ldb )
765 *
766  ELSE
767 *
768 * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
769 * TRANS = 'T'
770 *
771  CALL strsm( 'R', 'U', 'T', diag, m, n2, alpha,
772  $ a( n1 ), n, b( 0, n1 ), ldb )
773  CALL sgemm( 'N', 'T', m, n1, n2, -one, b( 0, n1 ),
774  $ ldb, a( 0 ), n, alpha, b( 0, 0 ), ldb )
775  CALL strsm( 'R', 'L', 'N', diag, m, n1, one,
776  $ a( n2 ), n, b( 0, 0 ), ldb )
777 *
778  END IF
779 *
780  END IF
781 *
782  ELSE
783 *
784 * SIDE = 'R', N is odd, and TRANSR = 'T'
785 *
786  IF( lower ) THEN
787 *
788 * SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'L'
789 *
790  IF( notrans ) THEN
791 *
792 * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and
793 * TRANS = 'N'
794 *
795  CALL strsm( 'R', 'L', 'N', diag, m, n2, alpha,
796  $ a( 1 ), n1, b( 0, n1 ), ldb )
797  CALL sgemm( 'N', 'T', m, n1, n2, -one, b( 0, n1 ),
798  $ ldb, a( n1*n1 ), n1, alpha, b( 0, 0 ),
799  $ ldb )
800  CALL strsm( 'R', 'U', 'T', diag, m, n1, one,
801  $ a( 0 ), n1, b( 0, 0 ), ldb )
802 *
803  ELSE
804 *
805 * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and
806 * TRANS = 'T'
807 *
808  CALL strsm( 'R', 'U', 'N', diag, m, n1, alpha,
809  $ a( 0 ), n1, b( 0, 0 ), ldb )
810  CALL sgemm( 'N', 'N', m, n2, n1, -one, b( 0, 0 ),
811  $ ldb, a( n1*n1 ), n1, alpha, b( 0, n1 ),
812  $ ldb )
813  CALL strsm( 'R', 'L', 'T', diag, m, n2, one,
814  $ a( 1 ), n1, b( 0, n1 ), ldb )
815 *
816  END IF
817 *
818  ELSE
819 *
820 * SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'U'
821 *
822  IF( notrans ) THEN
823 *
824 * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and
825 * TRANS = 'N'
826 *
827  CALL strsm( 'R', 'U', 'N', diag, m, n1, alpha,
828  $ a( n2*n2 ), n2, b( 0, 0 ), ldb )
829  CALL sgemm( 'N', 'T', m, n2, n1, -one, b( 0, 0 ),
830  $ ldb, a( 0 ), n2, alpha, b( 0, n1 ),
831  $ ldb )
832  CALL strsm( 'R', 'L', 'T', diag, m, n2, one,
833  $ a( n1*n2 ), n2, b( 0, n1 ), ldb )
834 *
835  ELSE
836 *
837 * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and
838 * TRANS = 'T'
839 *
840  CALL strsm( 'R', 'L', 'N', diag, m, n2, alpha,
841  $ a( n1*n2 ), n2, b( 0, n1 ), ldb )
842  CALL sgemm( 'N', 'N', m, n1, n2, -one, b( 0, n1 ),
843  $ ldb, a( 0 ), n2, alpha, b( 0, 0 ),
844  $ ldb )
845  CALL strsm( 'R', 'U', 'T', diag, m, n1, one,
846  $ a( n2*n2 ), n2, b( 0, 0 ), ldb )
847 *
848  END IF
849 *
850  END IF
851 *
852  END IF
853 *
854  ELSE
855 *
856 * SIDE = 'R' and N is even
857 *
858  IF( normaltransr ) THEN
859 *
860 * SIDE = 'R', N is even, and TRANSR = 'N'
861 *
862  IF( lower ) THEN
863 *
864 * SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'L'
865 *
866  IF( notrans ) THEN
867 *
868 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L',
869 * and TRANS = 'N'
870 *
871  CALL strsm( 'R', 'U', 'T', diag, m, k, alpha,
872  $ a( 0 ), n+1, b( 0, k ), ldb )
873  CALL sgemm( 'N', 'N', m, k, k, -one, b( 0, k ),
874  $ ldb, a( k+1 ), n+1, alpha, b( 0, 0 ),
875  $ ldb )
876  CALL strsm( 'R', 'L', 'N', diag, m, k, one,
877  $ a( 1 ), n+1, b( 0, 0 ), ldb )
878 *
879  ELSE
880 *
881 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L',
882 * and TRANS = 'T'
883 *
884  CALL strsm( 'R', 'L', 'T', diag, m, k, alpha,
885  $ a( 1 ), n+1, b( 0, 0 ), ldb )
886  CALL sgemm( 'N', 'T', m, k, k, -one, b( 0, 0 ),
887  $ ldb, a( k+1 ), n+1, alpha, b( 0, k ),
888  $ ldb )
889  CALL strsm( 'R', 'U', 'N', diag, m, k, one,
890  $ a( 0 ), n+1, b( 0, k ), ldb )
891 *
892  END IF
893 *
894  ELSE
895 *
896 * SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'U'
897 *
898  IF( notrans ) THEN
899 *
900 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U',
901 * and TRANS = 'N'
902 *
903  CALL strsm( 'R', 'L', 'T', diag, m, k, alpha,
904  $ a( k+1 ), n+1, b( 0, 0 ), ldb )
905  CALL sgemm( 'N', 'N', m, k, k, -one, b( 0, 0 ),
906  $ ldb, a( 0 ), n+1, alpha, b( 0, k ),
907  $ ldb )
908  CALL strsm( 'R', 'U', 'N', diag, m, k, one,
909  $ a( k ), n+1, b( 0, k ), ldb )
910 *
911  ELSE
912 *
913 * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U',
914 * and TRANS = 'T'
915 *
916  CALL strsm( 'R', 'U', 'T', diag, m, k, alpha,
917  $ a( k ), n+1, b( 0, k ), ldb )
918  CALL sgemm( 'N', 'T', m, k, k, -one, b( 0, k ),
919  $ ldb, a( 0 ), n+1, alpha, b( 0, 0 ),
920  $ ldb )
921  CALL strsm( 'R', 'L', 'N', diag, m, k, one,
922  $ a( k+1 ), n+1, b( 0, 0 ), ldb )
923 *
924  END IF
925 *
926  END IF
927 *
928  ELSE
929 *
930 * SIDE = 'R', N is even, and TRANSR = 'T'
931 *
932  IF( lower ) THEN
933 *
934 * SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'L'
935 *
936  IF( notrans ) THEN
937 *
938 * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L',
939 * and TRANS = 'N'
940 *
941  CALL strsm( 'R', 'L', 'N', diag, m, k, alpha,
942  $ a( 0 ), k, b( 0, k ), ldb )
943  CALL sgemm( 'N', 'T', m, k, k, -one, b( 0, k ),
944  $ ldb, a( ( k+1 )*k ), k, alpha,
945  $ b( 0, 0 ), ldb )
946  CALL strsm( 'R', 'U', 'T', diag, m, k, one,
947  $ a( k ), k, b( 0, 0 ), ldb )
948 *
949  ELSE
950 *
951 * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L',
952 * and TRANS = 'T'
953 *
954  CALL strsm( 'R', 'U', 'N', diag, m, k, alpha,
955  $ a( k ), k, b( 0, 0 ), ldb )
956  CALL sgemm( 'N', 'N', m, k, k, -one, b( 0, 0 ),
957  $ ldb, a( ( k+1 )*k ), k, alpha,
958  $ b( 0, k ), ldb )
959  CALL strsm( 'R', 'L', 'T', diag, m, k, one,
960  $ a( 0 ), k, b( 0, k ), ldb )
961 *
962  END IF
963 *
964  ELSE
965 *
966 * SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'U'
967 *
968  IF( notrans ) THEN
969 *
970 * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U',
971 * and TRANS = 'N'
972 *
973  CALL strsm( 'R', 'U', 'N', diag, m, k, alpha,
974  $ a( ( k+1 )*k ), k, b( 0, 0 ), ldb )
975  CALL sgemm( 'N', 'T', m, k, k, -one, b( 0, 0 ),
976  $ ldb, a( 0 ), k, alpha, b( 0, k ), ldb )
977  CALL strsm( 'R', 'L', 'T', diag, m, k, one,
978  $ a( k*k ), k, b( 0, k ), ldb )
979 *
980  ELSE
981 *
982 * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U',
983 * and TRANS = 'T'
984 *
985  CALL strsm( 'R', 'L', 'N', diag, m, k, alpha,
986  $ a( k*k ), k, b( 0, k ), ldb )
987  CALL sgemm( 'N', 'N', m, k, k, -one, b( 0, k ),
988  $ ldb, a( 0 ), k, alpha, b( 0, 0 ), ldb )
989  CALL strsm( 'R', 'U', 'T', diag, m, k, one,
990  $ a( ( k+1 )*k ), k, b( 0, 0 ), ldb )
991 *
992  END IF
993 *
994  END IF
995 *
996  END IF
997 *
998  END IF
999  END IF
1000 *
1001  RETURN
1002 *
1003 * End of STFSM
1004 *
1005  END
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
stfsm
subroutine stfsm(TRANSR, SIDE, UPLO, TRANS, DIAG, M, N, ALPHA, A, B, LDB)
STFSM solves a matrix equation (one operand is a triangular matrix in RFP format).
Definition: stfsm.f:279
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
strsm
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
Definition: strsm.f:183