1 SUBROUTINE ZPFTRS( 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 COMPLEX*16 A( 0: * ), B( LDB, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * ZPFTRS solves a system of linear equations A*X = B with a Hermitian
23 * positive definite matrix A using the Cholesky factorization
24 * A = U**H*U or A = L*L**H computed by ZPFTRF.
25 *
26 * Arguments
27 * =========
28 *
29 * TRANSR (input) CHARACTER*1
30 * = 'N': The Normal TRANSR of RFP A is stored;
31 * = 'C': The Conjugate-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) COMPLEX*16 array, dimension ( N*(N+1)/2 );
45 * The triangular factor U or L from the Cholesky factorization
46 * of RFP A = U**H*U or RFP A = L*L**H, as computed by ZPFTRF.
47 * See note below for more details about RFP A.
48 *
49 * B (input/output) COMPLEX*16 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 Standard Packed Format when N is even.
64 * 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 * conjugate-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 * conjugate-transpose of the last three columns of AP lower.
83 * To denote conjugate we place -- above the element. This covers the
84 * case N even and TRANSR = 'N'.
85 *
86 * RFP A RFP A
87 *
88 * -- -- --
89 * 03 04 05 33 43 53
90 * -- --
91 * 13 14 15 00 44 54
92 * --
93 * 23 24 25 10 11 55
94 *
95 * 33 34 35 20 21 22
96 * --
97 * 00 44 45 30 31 32
98 * -- --
99 * 01 11 55 40 41 42
100 * -- -- --
101 * 02 12 22 50 51 52
102 *
103 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
104 * transpose of RFP A above. One therefore gets:
105 *
106 *
107 * RFP A RFP A
108 *
109 * -- -- -- -- -- -- -- -- -- --
110 * 03 13 23 33 00 01 02 33 00 10 20 30 40 50
111 * -- -- -- -- -- -- -- -- -- --
112 * 04 14 24 34 44 11 12 43 44 11 21 31 41 51
113 * -- -- -- -- -- -- -- -- -- --
114 * 05 15 25 35 45 55 22 53 54 55 22 32 42 52
115 *
116 *
117 * We next consider Standard Packed Format when N is odd.
118 * We give an example where N = 5.
119 *
120 * AP is Upper AP is Lower
121 *
122 * 00 01 02 03 04 00
123 * 11 12 13 14 10 11
124 * 22 23 24 20 21 22
125 * 33 34 30 31 32 33
126 * 44 40 41 42 43 44
127 *
128 *
129 * Let TRANSR = 'N'. RFP holds AP as follows:
130 * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
131 * three columns of AP upper. The lower triangle A(3:4,0:1) consists of
132 * conjugate-transpose of the first two columns of AP upper.
133 * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
134 * three columns of AP lower. The upper triangle A(0:1,1:2) consists of
135 * conjugate-transpose of the last two columns of AP lower.
136 * To denote conjugate we place -- above the element. This covers the
137 * case N odd and TRANSR = 'N'.
138 *
139 * RFP A RFP A
140 *
141 * -- --
142 * 02 03 04 00 33 43
143 * --
144 * 12 13 14 10 11 44
145 *
146 * 22 23 24 20 21 22
147 * --
148 * 00 33 34 30 31 32
149 * -- --
150 * 01 11 44 40 41 42
151 *
152 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
153 * transpose of RFP A above. One therefore gets:
154 *
155 *
156 * RFP A RFP A
157 *
158 * -- -- -- -- -- -- -- -- --
159 * 02 12 22 00 01 00 10 20 30 40 50
160 * -- -- -- -- -- -- -- -- --
161 * 03 13 23 33 11 33 11 21 31 41 51
162 * -- -- -- -- -- -- -- -- --
163 * 04 14 24 34 44 43 44 22 32 42 52
164 *
165 * =====================================================================
166 *
167 * .. Parameters ..
168 COMPLEX*16 CONE
169 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
170 * ..
171 * .. Local Scalars ..
172 LOGICAL LOWER, NORMALTRANSR
173 * ..
174 * .. External Functions ..
175 LOGICAL LSAME
176 EXTERNAL LSAME
177 * ..
178 * .. External Subroutines ..
179 EXTERNAL XERBLA, ZTFSM
180 * ..
181 * .. Intrinsic Functions ..
182 INTRINSIC MAX
183 * ..
184 * .. Executable Statements ..
185 *
186 * Test the input parameters.
187 *
188 INFO = 0
189 NORMALTRANSR = LSAME( TRANSR, 'N' )
190 LOWER = LSAME( UPLO, 'L' )
191 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
192 INFO = -1
193 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
194 INFO = -2
195 ELSE IF( N.LT.0 ) THEN
196 INFO = -3
197 ELSE IF( NRHS.LT.0 ) THEN
198 INFO = -4
199 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
200 INFO = -7
201 END IF
202 IF( INFO.NE.0 ) THEN
203 CALL XERBLA( 'ZPFTRS', -INFO )
204 RETURN
205 END IF
206 *
207 * Quick return if possible
208 *
209 IF( N.EQ.0 .OR. NRHS.EQ.0 )
210 $ RETURN
211 *
212 * start execution: there are two triangular solves
213 *
214 IF( LOWER ) THEN
215 CALL ZTFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, CONE, A, B,
216 $ LDB )
217 CALL ZTFSM( TRANSR, 'L', UPLO, 'C', 'N', N, NRHS, CONE, A, B,
218 $ LDB )
219 ELSE
220 CALL ZTFSM( TRANSR, 'L', UPLO, 'C', 'N', N, NRHS, CONE, A, B,
221 $ LDB )
222 CALL ZTFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, CONE, A, B,
223 $ LDB )
224 END IF
225 *
226 RETURN
227 *
228 * End of ZPFTRS
229 *
230 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 COMPLEX*16 A( 0: * ), B( LDB, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * ZPFTRS solves a system of linear equations A*X = B with a Hermitian
23 * positive definite matrix A using the Cholesky factorization
24 * A = U**H*U or A = L*L**H computed by ZPFTRF.
25 *
26 * Arguments
27 * =========
28 *
29 * TRANSR (input) CHARACTER*1
30 * = 'N': The Normal TRANSR of RFP A is stored;
31 * = 'C': The Conjugate-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) COMPLEX*16 array, dimension ( N*(N+1)/2 );
45 * The triangular factor U or L from the Cholesky factorization
46 * of RFP A = U**H*U or RFP A = L*L**H, as computed by ZPFTRF.
47 * See note below for more details about RFP A.
48 *
49 * B (input/output) COMPLEX*16 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 Standard Packed Format when N is even.
64 * 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 * conjugate-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 * conjugate-transpose of the last three columns of AP lower.
83 * To denote conjugate we place -- above the element. This covers the
84 * case N even and TRANSR = 'N'.
85 *
86 * RFP A RFP A
87 *
88 * -- -- --
89 * 03 04 05 33 43 53
90 * -- --
91 * 13 14 15 00 44 54
92 * --
93 * 23 24 25 10 11 55
94 *
95 * 33 34 35 20 21 22
96 * --
97 * 00 44 45 30 31 32
98 * -- --
99 * 01 11 55 40 41 42
100 * -- -- --
101 * 02 12 22 50 51 52
102 *
103 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
104 * transpose of RFP A above. One therefore gets:
105 *
106 *
107 * RFP A RFP A
108 *
109 * -- -- -- -- -- -- -- -- -- --
110 * 03 13 23 33 00 01 02 33 00 10 20 30 40 50
111 * -- -- -- -- -- -- -- -- -- --
112 * 04 14 24 34 44 11 12 43 44 11 21 31 41 51
113 * -- -- -- -- -- -- -- -- -- --
114 * 05 15 25 35 45 55 22 53 54 55 22 32 42 52
115 *
116 *
117 * We next consider Standard Packed Format when N is odd.
118 * We give an example where N = 5.
119 *
120 * AP is Upper AP is Lower
121 *
122 * 00 01 02 03 04 00
123 * 11 12 13 14 10 11
124 * 22 23 24 20 21 22
125 * 33 34 30 31 32 33
126 * 44 40 41 42 43 44
127 *
128 *
129 * Let TRANSR = 'N'. RFP holds AP as follows:
130 * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
131 * three columns of AP upper. The lower triangle A(3:4,0:1) consists of
132 * conjugate-transpose of the first two columns of AP upper.
133 * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
134 * three columns of AP lower. The upper triangle A(0:1,1:2) consists of
135 * conjugate-transpose of the last two columns of AP lower.
136 * To denote conjugate we place -- above the element. This covers the
137 * case N odd and TRANSR = 'N'.
138 *
139 * RFP A RFP A
140 *
141 * -- --
142 * 02 03 04 00 33 43
143 * --
144 * 12 13 14 10 11 44
145 *
146 * 22 23 24 20 21 22
147 * --
148 * 00 33 34 30 31 32
149 * -- --
150 * 01 11 44 40 41 42
151 *
152 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
153 * transpose of RFP A above. One therefore gets:
154 *
155 *
156 * RFP A RFP A
157 *
158 * -- -- -- -- -- -- -- -- --
159 * 02 12 22 00 01 00 10 20 30 40 50
160 * -- -- -- -- -- -- -- -- --
161 * 03 13 23 33 11 33 11 21 31 41 51
162 * -- -- -- -- -- -- -- -- --
163 * 04 14 24 34 44 43 44 22 32 42 52
164 *
165 * =====================================================================
166 *
167 * .. Parameters ..
168 COMPLEX*16 CONE
169 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
170 * ..
171 * .. Local Scalars ..
172 LOGICAL LOWER, NORMALTRANSR
173 * ..
174 * .. External Functions ..
175 LOGICAL LSAME
176 EXTERNAL LSAME
177 * ..
178 * .. External Subroutines ..
179 EXTERNAL XERBLA, ZTFSM
180 * ..
181 * .. Intrinsic Functions ..
182 INTRINSIC MAX
183 * ..
184 * .. Executable Statements ..
185 *
186 * Test the input parameters.
187 *
188 INFO = 0
189 NORMALTRANSR = LSAME( TRANSR, 'N' )
190 LOWER = LSAME( UPLO, 'L' )
191 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
192 INFO = -1
193 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
194 INFO = -2
195 ELSE IF( N.LT.0 ) THEN
196 INFO = -3
197 ELSE IF( NRHS.LT.0 ) THEN
198 INFO = -4
199 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
200 INFO = -7
201 END IF
202 IF( INFO.NE.0 ) THEN
203 CALL XERBLA( 'ZPFTRS', -INFO )
204 RETURN
205 END IF
206 *
207 * Quick return if possible
208 *
209 IF( N.EQ.0 .OR. NRHS.EQ.0 )
210 $ RETURN
211 *
212 * start execution: there are two triangular solves
213 *
214 IF( LOWER ) THEN
215 CALL ZTFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, CONE, A, B,
216 $ LDB )
217 CALL ZTFSM( TRANSR, 'L', UPLO, 'C', 'N', N, NRHS, CONE, A, B,
218 $ LDB )
219 ELSE
220 CALL ZTFSM( TRANSR, 'L', UPLO, 'C', 'N', N, NRHS, CONE, A, B,
221 $ LDB )
222 CALL ZTFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, CONE, A, B,
223 $ LDB )
224 END IF
225 *
226 RETURN
227 *
228 * End of ZPFTRS
229 *
230 END