001:       SUBROUTINE SPFTRS( TRANSR, UPLO, N, NRHS, A, B, LDB, INFO )
002: *
003: *  -- LAPACK routine (version 3.2)                                    --
004: *
005: *  -- Contributed by Fred Gustavson of the IBM Watson Research Center --
006: *  -- November 2008                                                   --
007: *
008: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
009: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
010: *
011: *     .. Scalar Arguments ..
012:       CHARACTER          TRANSR, UPLO
013:       INTEGER            INFO, LDB, N, NRHS
014: *     ..
015: *     .. Array Arguments ..
016:       REAL               A( 0: * ), B( LDB, * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  SPFTRS solves a system of linear equations A*X = B with a symmetric
023: *  positive definite matrix A using the Cholesky factorization
024: *  A = U**T*U or A = L*L**T computed by SPFTRF.
025: *
026: *  Arguments
027: *  =========
028: *
029: *  TRANSR    (input) CHARACTER
030: *          = 'N':  The Normal TRANSR of RFP A is stored;
031: *          = 'T':  The Transpose TRANSR of RFP A is stored.
032: *
033: *  UPLO    (input) CHARACTER
034: *          = 'U':  Upper triangle of RFP A is stored;
035: *          = 'L':  Lower triangle of RFP A is stored.
036: *
037: *  N       (input) INTEGER
038: *          The order of the matrix A.  N >= 0.
039: *
040: *  NRHS    (input) INTEGER
041: *          The number of right hand sides, i.e., the number of columns
042: *          of the matrix B.  NRHS >= 0.
043: *
044: *  A       (input) REAL array, dimension ( N*(N+1)/2 )
045: *          The triangular factor U or L from the Cholesky factorization
046: *          of RFP A = U**H*U or RFP A = L*L**T, as computed by SPFTRF.
047: *          See note below for more details about RFP A.
048: *
049: *  B       (input/output) REAL array, dimension (LDB,NRHS)
050: *          On entry, the right hand side matrix B.
051: *          On exit, the solution matrix X.
052: *
053: *  LDB     (input) INTEGER
054: *          The leading dimension of the array B.  LDB >= max(1,N).
055: *
056: *  INFO    (output) INTEGER
057: *          = 0:  successful exit
058: *          < 0:  if INFO = -i, the i-th argument had an illegal value
059: *
060: *  Notes
061: *  =====
062: *
063: *  We first consider Rectangular Full Packed (RFP) Format when N is
064: *  even. We give an example where N = 6.
065: *
066: *      AP is Upper             AP is Lower
067: *
068: *   00 01 02 03 04 05       00
069: *      11 12 13 14 15       10 11
070: *         22 23 24 25       20 21 22
071: *            33 34 35       30 31 32 33
072: *               44 45       40 41 42 43 44
073: *                  55       50 51 52 53 54 55
074: *
075: *
076: *  Let TRANSR = 'N'. RFP holds AP as follows:
077: *  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
078: *  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
079: *  the transpose of the first three columns of AP upper.
080: *  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
081: *  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
082: *  the transpose of the last three columns of AP lower.
083: *  This covers the case N even and TRANSR = 'N'.
084: *
085: *         RFP A                   RFP A
086: *
087: *        03 04 05                33 43 53
088: *        13 14 15                00 44 54
089: *        23 24 25                10 11 55
090: *        33 34 35                20 21 22
091: *        00 44 45                30 31 32
092: *        01 11 55                40 41 42
093: *        02 12 22                50 51 52
094: *
095: *  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
096: *  transpose of RFP A above. One therefore gets:
097: *
098: *
099: *           RFP A                   RFP A
100: *
101: *     03 13 23 33 00 01 02    33 00 10 20 30 40 50
102: *     04 14 24 34 44 11 12    43 44 11 21 31 41 51
103: *     05 15 25 35 45 55 22    53 54 55 22 32 42 52
104: *
105: *
106: *  We first consider Rectangular Full Packed (RFP) Format when N is
107: *  odd. We give an example where N = 5.
108: *
109: *     AP is Upper                 AP is Lower
110: *
111: *   00 01 02 03 04              00
112: *      11 12 13 14              10 11
113: *         22 23 24              20 21 22
114: *            33 34              30 31 32 33
115: *               44              40 41 42 43 44
116: *
117: *
118: *  Let TRANSR = 'N'. RFP holds AP as follows:
119: *  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
120: *  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
121: *  the transpose of the first two columns of AP upper.
122: *  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
123: *  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
124: *  the transpose of the last two columns of AP lower.
125: *  This covers the case N odd and TRANSR = 'N'.
126: *
127: *         RFP A                   RFP A
128: *
129: *        02 03 04                00 33 43
130: *        12 13 14                10 11 44
131: *        22 23 24                20 21 22
132: *        00 33 34                30 31 32
133: *        01 11 44                40 41 42
134: *
135: *  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
136: *  transpose of RFP A above. One therefore gets:
137: *
138: *           RFP A                   RFP A
139: *
140: *     02 12 22 00 01             00 10 20 30 40 50
141: *     03 13 23 33 11             33 11 21 31 41 51
142: *     04 14 24 34 44             43 44 22 32 42 52
143: *
144: *  =====================================================================
145: *
146: *     .. Parameters ..
147:       REAL               ONE
148:       PARAMETER          ( ONE = 1.0E+0 )
149: *     ..
150: *     .. Local Scalars ..
151:       LOGICAL            LOWER, NORMALTRANSR
152: *     ..
153: *     .. External Functions ..
154:       LOGICAL            LSAME
155:       EXTERNAL           LSAME
156: *     ..
157: *     .. External Subroutines ..
158:       EXTERNAL           XERBLA, STFSM
159: *     ..
160: *     .. Intrinsic Functions ..
161:       INTRINSIC          MAX
162: *     ..
163: *     .. Executable Statements ..
164: *
165: *     Test the input parameters.
166: *
167:       INFO = 0
168:       NORMALTRANSR = LSAME( TRANSR, 'N' )
169:       LOWER = LSAME( UPLO, 'L' )
170:       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
171:          INFO = -1
172:       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
173:          INFO = -2
174:       ELSE IF( N.LT.0 ) THEN
175:          INFO = -3
176:       ELSE IF( NRHS.LT.0 ) THEN
177:          INFO = -4
178:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
179:          INFO = -7
180:       END IF
181:       IF( INFO.NE.0 ) THEN
182:          CALL XERBLA( 'SPFTRS', -INFO )
183:          RETURN
184:       END IF
185: *
186: *     Quick return if possible
187: *
188:       IF( N.EQ.0 .OR. NRHS.EQ.0 )
189:      +   RETURN
190: *
191: *     start execution: there are two triangular solves
192: *
193:       IF( LOWER ) THEN
194:          CALL STFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, ONE, A, B,
195:      +               LDB )
196:          CALL STFSM( TRANSR, 'L', UPLO, 'T', 'N', N, NRHS, ONE, A, B,
197:      +               LDB )
198:       ELSE
199:          CALL STFSM( TRANSR, 'L', UPLO, 'T', 'N', N, NRHS, ONE, A, B,
200:      +               LDB )
201:          CALL STFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, ONE, A, B,
202:      +               LDB )
203:       END IF
204: *
205:       RETURN
206: *
207: *     End of SPFTRS
208: *
209:       END
210: