1 SUBROUTINE DPFTRS( TRANSR, UPLO, N, NRHS, A, B, LDB, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 *
5 * -- Contributed by Fred Gustavson of the IBM Watson Research Center --
6 * -- April 2011 --
7 *
8 * -- LAPACK is a software package provided by Univ. of Tennessee, --
9 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
10 *
11 * .. Scalar Arguments ..
12 CHARACTER TRANSR, UPLO
13 INTEGER INFO, LDB, N, NRHS
14 * ..
15 * .. Array Arguments ..
16 DOUBLE PRECISION A( 0: * ), B( LDB, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * DPFTRS solves a system of linear equations A*X = B with a symmetric
23 * positive definite matrix A using the Cholesky factorization
24 * A = U**T*U or A = L*L**T computed by DPFTRF.
25 *
26 * Arguments
27 * =========
28 *
29 * TRANSR (input) CHARACTER*1
30 * = 'N': The Normal TRANSR of RFP A is stored;
31 * = 'T': The Transpose TRANSR of RFP A is stored.
32 *
33 * UPLO (input) CHARACTER*1
34 * = 'U': Upper triangle of RFP A is stored;
35 * = 'L': Lower triangle of RFP A is stored.
36 *
37 * N (input) INTEGER
38 * The order of the matrix A. N >= 0.
39 *
40 * NRHS (input) INTEGER
41 * The number of right hand sides, i.e., the number of columns
42 * of the matrix B. NRHS >= 0.
43 *
44 * A (input) DOUBLE PRECISION array, dimension ( N*(N+1)/2 ).
45 * The triangular factor U or L from the Cholesky factorization
46 * of RFP A = U**T*U or RFP A = L*L**T, as computed by DPFTRF.
47 * See note below for more details about RFP A.
48 *
49 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
50 * On entry, the right hand side matrix B.
51 * On exit, the solution matrix X.
52 *
53 * LDB (input) INTEGER
54 * The leading dimension of the array B. LDB >= max(1,N).
55 *
56 * INFO (output) INTEGER
57 * = 0: successful exit
58 * < 0: if INFO = -i, the i-th argument had an illegal value
59 *
60 * Further Details
61 * ===============
62 *
63 * We first consider Rectangular Full Packed (RFP) Format when N is
64 * even. We give an example where N = 6.
65 *
66 * AP is Upper AP is Lower
67 *
68 * 00 01 02 03 04 05 00
69 * 11 12 13 14 15 10 11
70 * 22 23 24 25 20 21 22
71 * 33 34 35 30 31 32 33
72 * 44 45 40 41 42 43 44
73 * 55 50 51 52 53 54 55
74 *
75 *
76 * Let TRANSR = 'N'. RFP holds AP as follows:
77 * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
78 * three columns of AP upper. The lower triangle A(4:6,0:2) consists of
79 * the transpose of the first three columns of AP upper.
80 * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
81 * three columns of AP lower. The upper triangle A(0:2,0:2) consists of
82 * the transpose of the last three columns of AP lower.
83 * This covers the case N even and TRANSR = 'N'.
84 *
85 * RFP A RFP A
86 *
87 * 03 04 05 33 43 53
88 * 13 14 15 00 44 54
89 * 23 24 25 10 11 55
90 * 33 34 35 20 21 22
91 * 00 44 45 30 31 32
92 * 01 11 55 40 41 42
93 * 02 12 22 50 51 52
94 *
95 * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
96 * transpose of RFP A above. One therefore gets:
97 *
98 *
99 * 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 then 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 DOUBLE PRECISION ONE
148 PARAMETER ( ONE = 1.0D+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, DTFSM
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( 'DPFTRS', -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 DTFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, ONE, A, B,
195 $ LDB )
196 CALL DTFSM( TRANSR, 'L', UPLO, 'T', 'N', N, NRHS, ONE, A, B,
197 $ LDB )
198 ELSE
199 CALL DTFSM( TRANSR, 'L', UPLO, 'T', 'N', N, NRHS, ONE, A, B,
200 $ LDB )
201 CALL DTFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, ONE, A, B,
202 $ LDB )
203 END IF
204 *
205 RETURN
206 *
207 * End of DPFTRS
208 *
209 END
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 *
5 * -- Contributed by Fred Gustavson of the IBM Watson Research Center --
6 * -- April 2011 --
7 *
8 * -- LAPACK is a software package provided by Univ. of Tennessee, --
9 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
10 *
11 * .. Scalar Arguments ..
12 CHARACTER TRANSR, UPLO
13 INTEGER INFO, LDB, N, NRHS
14 * ..
15 * .. Array Arguments ..
16 DOUBLE PRECISION A( 0: * ), B( LDB, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * DPFTRS solves a system of linear equations A*X = B with a symmetric
23 * positive definite matrix A using the Cholesky factorization
24 * A = U**T*U or A = L*L**T computed by DPFTRF.
25 *
26 * Arguments
27 * =========
28 *
29 * TRANSR (input) CHARACTER*1
30 * = 'N': The Normal TRANSR of RFP A is stored;
31 * = 'T': The Transpose TRANSR of RFP A is stored.
32 *
33 * UPLO (input) CHARACTER*1
34 * = 'U': Upper triangle of RFP A is stored;
35 * = 'L': Lower triangle of RFP A is stored.
36 *
37 * N (input) INTEGER
38 * The order of the matrix A. N >= 0.
39 *
40 * NRHS (input) INTEGER
41 * The number of right hand sides, i.e., the number of columns
42 * of the matrix B. NRHS >= 0.
43 *
44 * A (input) DOUBLE PRECISION array, dimension ( N*(N+1)/2 ).
45 * The triangular factor U or L from the Cholesky factorization
46 * of RFP A = U**T*U or RFP A = L*L**T, as computed by DPFTRF.
47 * See note below for more details about RFP A.
48 *
49 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
50 * On entry, the right hand side matrix B.
51 * On exit, the solution matrix X.
52 *
53 * LDB (input) INTEGER
54 * The leading dimension of the array B. LDB >= max(1,N).
55 *
56 * INFO (output) INTEGER
57 * = 0: successful exit
58 * < 0: if INFO = -i, the i-th argument had an illegal value
59 *
60 * Further Details
61 * ===============
62 *
63 * We first consider Rectangular Full Packed (RFP) Format when N is
64 * even. We give an example where N = 6.
65 *
66 * AP is Upper AP is Lower
67 *
68 * 00 01 02 03 04 05 00
69 * 11 12 13 14 15 10 11
70 * 22 23 24 25 20 21 22
71 * 33 34 35 30 31 32 33
72 * 44 45 40 41 42 43 44
73 * 55 50 51 52 53 54 55
74 *
75 *
76 * Let TRANSR = 'N'. RFP holds AP as follows:
77 * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
78 * three columns of AP upper. The lower triangle A(4:6,0:2) consists of
79 * the transpose of the first three columns of AP upper.
80 * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
81 * three columns of AP lower. The upper triangle A(0:2,0:2) consists of
82 * the transpose of the last three columns of AP lower.
83 * This covers the case N even and TRANSR = 'N'.
84 *
85 * RFP A RFP A
86 *
87 * 03 04 05 33 43 53
88 * 13 14 15 00 44 54
89 * 23 24 25 10 11 55
90 * 33 34 35 20 21 22
91 * 00 44 45 30 31 32
92 * 01 11 55 40 41 42
93 * 02 12 22 50 51 52
94 *
95 * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
96 * transpose of RFP A above. One therefore gets:
97 *
98 *
99 * 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 then 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 DOUBLE PRECISION ONE
148 PARAMETER ( ONE = 1.0D+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, DTFSM
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( 'DPFTRS', -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 DTFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, ONE, A, B,
195 $ LDB )
196 CALL DTFSM( TRANSR, 'L', UPLO, 'T', 'N', N, NRHS, ONE, A, B,
197 $ LDB )
198 ELSE
199 CALL DTFSM( TRANSR, 'L', UPLO, 'T', 'N', N, NRHS, ONE, A, B,
200 $ LDB )
201 CALL DTFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, ONE, A, B,
202 $ LDB )
203 END IF
204 *
205 RETURN
206 *
207 * End of DPFTRS
208 *
209 END