1 SUBROUTINE ZSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
2 $ WORK, INFO )
3 *
4 * -- LAPACK PROTOTYPE routine (version 3.3.1) --
5 *
6 * -- Written by Julie Langou of the Univ. of TN --
7 * -- April 2011 --
8 *
9 * -- LAPACK is a software package provided by Univ. of Tennessee, --
10 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
11 *
12 * .. Scalar Arguments ..
13 CHARACTER UPLO
14 INTEGER INFO, LDA, LDB, N, NRHS
15 * ..
16 * .. Array Arguments ..
17 INTEGER IPIV( * )
18 DOUBLE COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * ZSYTRS2 solves a system of linear equations A*X = B with a real
25 * symmetric matrix A using the factorization A = U*D*U**T or
26 * A = L*D*L**T computed by ZSYTRF and converted by ZSYCONV.
27 *
28 * Arguments
29 * =========
30 *
31 * UPLO (input) CHARACTER*1
32 * Specifies whether the details of the factorization are stored
33 * as an upper or lower triangular matrix.
34 * = 'U': Upper triangular, form is A = U*D*U**T;
35 * = 'L': Lower triangular, form is A = L*D*L**T.
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 COMPLEX array, dimension (LDA,N)
45 * The block diagonal matrix D and the multipliers used to
46 * obtain the factor U or L as computed by ZSYTRF.
47 *
48 * LDA (input) INTEGER
49 * The leading dimension of the array A. LDA >= max(1,N).
50 *
51 * IPIV (input) INTEGER array, dimension (N)
52 * Details of the interchanges and the block structure of D
53 * as determined by ZSYTRF.
54 *
55 * B (input/output) DOUBLE COMPLEX array, dimension (LDB,NRHS)
56 * On entry, the right hand side matrix B.
57 * On exit, the solution matrix X.
58 *
59 * LDB (input) INTEGER
60 * The leading dimension of the array B. LDB >= max(1,N).
61 *
62 * WORK (workspace) REAL array, dimension (N)
63 *
64 * INFO (output) INTEGER
65 * = 0: successful exit
66 * < 0: if INFO = -i, the i-th argument had an illegal value
67 *
68 * =====================================================================
69 *
70 * .. Parameters ..
71 DOUBLE COMPLEX ONE
72 PARAMETER ( ONE = (1.0D+0,0.0D+0) )
73 * ..
74 * .. Local Scalars ..
75 LOGICAL UPPER
76 INTEGER I, IINFO, J, K, KP
77 DOUBLE COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM
78 * ..
79 * .. External Functions ..
80 LOGICAL LSAME
81 EXTERNAL LSAME
82 * ..
83 * .. External Subroutines ..
84 EXTERNAL ZSCAL, ZSYCONV, ZSWAP, ZTRSM, XERBLA
85 * ..
86 * .. Intrinsic Functions ..
87 INTRINSIC MAX
88 * ..
89 * .. Executable Statements ..
90 *
91 INFO = 0
92 UPPER = LSAME( UPLO, 'U' )
93 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
94 INFO = -1
95 ELSE IF( N.LT.0 ) THEN
96 INFO = -2
97 ELSE IF( NRHS.LT.0 ) THEN
98 INFO = -3
99 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
100 INFO = -5
101 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
102 INFO = -8
103 END IF
104 IF( INFO.NE.0 ) THEN
105 CALL XERBLA( 'ZSYTRS2', -INFO )
106 RETURN
107 END IF
108 *
109 * Quick return if possible
110 *
111 IF( N.EQ.0 .OR. NRHS.EQ.0 )
112 $ RETURN
113 *
114 * Convert A
115 *
116 CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
117 *
118 IF( UPPER ) THEN
119 *
120 * Solve A*X = B, where A = U*D*U**T.
121 *
122 * P**T * B
123 K=N
124 DO WHILE ( K .GE. 1 )
125 IF( IPIV( K ).GT.0 ) THEN
126 * 1 x 1 diagonal block
127 * Interchange rows K and IPIV(K).
128 KP = IPIV( K )
129 IF( KP.NE.K )
130 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
131 K=K-1
132 ELSE
133 * 2 x 2 diagonal block
134 * Interchange rows K-1 and -IPIV(K).
135 KP = -IPIV( K )
136 IF( KP.EQ.-IPIV( K-1 ) )
137 $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
138 K=K-2
139 END IF
140 END DO
141 *
142 * Compute (U \P**T * B) -> B [ (U \P**T * B) ]
143 *
144 CALL ZTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB)
145 *
146 * Compute D \ B -> B [ D \ (U \P**T * B) ]
147 *
148 I=N
149 DO WHILE ( I .GE. 1 )
150 IF( IPIV(I) .GT. 0 ) THEN
151 CALL ZSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
152 ELSEIF ( I .GT. 1) THEN
153 IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
154 AKM1K = WORK(I)
155 AKM1 = A( I-1, I-1 ) / AKM1K
156 AK = A( I, I ) / AKM1K
157 DENOM = AKM1*AK - ONE
158 DO 15 J = 1, NRHS
159 BKM1 = B( I-1, J ) / AKM1K
160 BK = B( I, J ) / AKM1K
161 B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
162 B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
163 15 CONTINUE
164 I = I - 1
165 ENDIF
166 ENDIF
167 I = I - 1
168 END DO
169 *
170 * Compute (U**T \ B) -> B [ U**T \ (D \ (U \P**T * B) ) ]
171 *
172 CALL ZTRSM('L','U','T','U',N,NRHS,ONE,A,LDA,B,LDB)
173 *
174 * P * B [ P * (U**T \ (D \ (U \P**T * B) )) ]
175 *
176 K=1
177 DO WHILE ( K .LE. N )
178 IF( IPIV( K ).GT.0 ) THEN
179 * 1 x 1 diagonal block
180 * Interchange rows K and IPIV(K).
181 KP = IPIV( K )
182 IF( KP.NE.K )
183 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
184 K=K+1
185 ELSE
186 * 2 x 2 diagonal block
187 * Interchange rows K-1 and -IPIV(K).
188 KP = -IPIV( K )
189 IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
190 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
191 K=K+2
192 ENDIF
193 END DO
194 *
195 ELSE
196 *
197 * Solve A*X = B, where A = L*D*L**T.
198 *
199 * P**T * B
200 K=1
201 DO WHILE ( K .LE. N )
202 IF( IPIV( K ).GT.0 ) THEN
203 * 1 x 1 diagonal block
204 * Interchange rows K and IPIV(K).
205 KP = IPIV( K )
206 IF( KP.NE.K )
207 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
208 K=K+1
209 ELSE
210 * 2 x 2 diagonal block
211 * Interchange rows K and -IPIV(K+1).
212 KP = -IPIV( K+1 )
213 IF( KP.EQ.-IPIV( K ) )
214 $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
215 K=K+2
216 ENDIF
217 END DO
218 *
219 * Compute (L \P**T * B) -> B [ (L \P**T * B) ]
220 *
221 CALL ZTRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB)
222 *
223 * Compute D \ B -> B [ D \ (L \P**T * B) ]
224 *
225 I=1
226 DO WHILE ( I .LE. N )
227 IF( IPIV(I) .GT. 0 ) THEN
228 CALL ZSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
229 ELSE
230 AKM1K = WORK(I)
231 AKM1 = A( I, I ) / AKM1K
232 AK = A( I+1, I+1 ) / AKM1K
233 DENOM = AKM1*AK - ONE
234 DO 25 J = 1, NRHS
235 BKM1 = B( I, J ) / AKM1K
236 BK = B( I+1, J ) / AKM1K
237 B( I, J ) = ( AK*BKM1-BK ) / DENOM
238 B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
239 25 CONTINUE
240 I = I + 1
241 ENDIF
242 I = I + 1
243 END DO
244 *
245 * Compute (L**T \ B) -> B [ L**T \ (D \ (L \P**T * B) ) ]
246 *
247 CALL ZTRSM('L','L','T','U',N,NRHS,ONE,A,LDA,B,LDB)
248 *
249 * P * B [ P * (L**T \ (D \ (L \P**T * B) )) ]
250 *
251 K=N
252 DO WHILE ( K .GE. 1 )
253 IF( IPIV( K ).GT.0 ) THEN
254 * 1 x 1 diagonal block
255 * Interchange rows K and IPIV(K).
256 KP = IPIV( K )
257 IF( KP.NE.K )
258 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
259 K=K-1
260 ELSE
261 * 2 x 2 diagonal block
262 * Interchange rows K-1 and -IPIV(K).
263 KP = -IPIV( K )
264 IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
265 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
266 K=K-2
267 ENDIF
268 END DO
269 *
270 END IF
271 *
272 * Revert A
273 *
274 CALL ZSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
275 *
276 RETURN
277 *
278 * End of ZSYTRS2
279 *
280 END
2 $ WORK, INFO )
3 *
4 * -- LAPACK PROTOTYPE routine (version 3.3.1) --
5 *
6 * -- Written by Julie Langou of the Univ. of TN --
7 * -- April 2011 --
8 *
9 * -- LAPACK is a software package provided by Univ. of Tennessee, --
10 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
11 *
12 * .. Scalar Arguments ..
13 CHARACTER UPLO
14 INTEGER INFO, LDA, LDB, N, NRHS
15 * ..
16 * .. Array Arguments ..
17 INTEGER IPIV( * )
18 DOUBLE COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * ZSYTRS2 solves a system of linear equations A*X = B with a real
25 * symmetric matrix A using the factorization A = U*D*U**T or
26 * A = L*D*L**T computed by ZSYTRF and converted by ZSYCONV.
27 *
28 * Arguments
29 * =========
30 *
31 * UPLO (input) CHARACTER*1
32 * Specifies whether the details of the factorization are stored
33 * as an upper or lower triangular matrix.
34 * = 'U': Upper triangular, form is A = U*D*U**T;
35 * = 'L': Lower triangular, form is A = L*D*L**T.
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 COMPLEX array, dimension (LDA,N)
45 * The block diagonal matrix D and the multipliers used to
46 * obtain the factor U or L as computed by ZSYTRF.
47 *
48 * LDA (input) INTEGER
49 * The leading dimension of the array A. LDA >= max(1,N).
50 *
51 * IPIV (input) INTEGER array, dimension (N)
52 * Details of the interchanges and the block structure of D
53 * as determined by ZSYTRF.
54 *
55 * B (input/output) DOUBLE COMPLEX array, dimension (LDB,NRHS)
56 * On entry, the right hand side matrix B.
57 * On exit, the solution matrix X.
58 *
59 * LDB (input) INTEGER
60 * The leading dimension of the array B. LDB >= max(1,N).
61 *
62 * WORK (workspace) REAL array, dimension (N)
63 *
64 * INFO (output) INTEGER
65 * = 0: successful exit
66 * < 0: if INFO = -i, the i-th argument had an illegal value
67 *
68 * =====================================================================
69 *
70 * .. Parameters ..
71 DOUBLE COMPLEX ONE
72 PARAMETER ( ONE = (1.0D+0,0.0D+0) )
73 * ..
74 * .. Local Scalars ..
75 LOGICAL UPPER
76 INTEGER I, IINFO, J, K, KP
77 DOUBLE COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM
78 * ..
79 * .. External Functions ..
80 LOGICAL LSAME
81 EXTERNAL LSAME
82 * ..
83 * .. External Subroutines ..
84 EXTERNAL ZSCAL, ZSYCONV, ZSWAP, ZTRSM, XERBLA
85 * ..
86 * .. Intrinsic Functions ..
87 INTRINSIC MAX
88 * ..
89 * .. Executable Statements ..
90 *
91 INFO = 0
92 UPPER = LSAME( UPLO, 'U' )
93 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
94 INFO = -1
95 ELSE IF( N.LT.0 ) THEN
96 INFO = -2
97 ELSE IF( NRHS.LT.0 ) THEN
98 INFO = -3
99 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
100 INFO = -5
101 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
102 INFO = -8
103 END IF
104 IF( INFO.NE.0 ) THEN
105 CALL XERBLA( 'ZSYTRS2', -INFO )
106 RETURN
107 END IF
108 *
109 * Quick return if possible
110 *
111 IF( N.EQ.0 .OR. NRHS.EQ.0 )
112 $ RETURN
113 *
114 * Convert A
115 *
116 CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
117 *
118 IF( UPPER ) THEN
119 *
120 * Solve A*X = B, where A = U*D*U**T.
121 *
122 * P**T * B
123 K=N
124 DO WHILE ( K .GE. 1 )
125 IF( IPIV( K ).GT.0 ) THEN
126 * 1 x 1 diagonal block
127 * Interchange rows K and IPIV(K).
128 KP = IPIV( K )
129 IF( KP.NE.K )
130 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
131 K=K-1
132 ELSE
133 * 2 x 2 diagonal block
134 * Interchange rows K-1 and -IPIV(K).
135 KP = -IPIV( K )
136 IF( KP.EQ.-IPIV( K-1 ) )
137 $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
138 K=K-2
139 END IF
140 END DO
141 *
142 * Compute (U \P**T * B) -> B [ (U \P**T * B) ]
143 *
144 CALL ZTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB)
145 *
146 * Compute D \ B -> B [ D \ (U \P**T * B) ]
147 *
148 I=N
149 DO WHILE ( I .GE. 1 )
150 IF( IPIV(I) .GT. 0 ) THEN
151 CALL ZSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
152 ELSEIF ( I .GT. 1) THEN
153 IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
154 AKM1K = WORK(I)
155 AKM1 = A( I-1, I-1 ) / AKM1K
156 AK = A( I, I ) / AKM1K
157 DENOM = AKM1*AK - ONE
158 DO 15 J = 1, NRHS
159 BKM1 = B( I-1, J ) / AKM1K
160 BK = B( I, J ) / AKM1K
161 B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
162 B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
163 15 CONTINUE
164 I = I - 1
165 ENDIF
166 ENDIF
167 I = I - 1
168 END DO
169 *
170 * Compute (U**T \ B) -> B [ U**T \ (D \ (U \P**T * B) ) ]
171 *
172 CALL ZTRSM('L','U','T','U',N,NRHS,ONE,A,LDA,B,LDB)
173 *
174 * P * B [ P * (U**T \ (D \ (U \P**T * B) )) ]
175 *
176 K=1
177 DO WHILE ( K .LE. N )
178 IF( IPIV( K ).GT.0 ) THEN
179 * 1 x 1 diagonal block
180 * Interchange rows K and IPIV(K).
181 KP = IPIV( K )
182 IF( KP.NE.K )
183 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
184 K=K+1
185 ELSE
186 * 2 x 2 diagonal block
187 * Interchange rows K-1 and -IPIV(K).
188 KP = -IPIV( K )
189 IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
190 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
191 K=K+2
192 ENDIF
193 END DO
194 *
195 ELSE
196 *
197 * Solve A*X = B, where A = L*D*L**T.
198 *
199 * P**T * B
200 K=1
201 DO WHILE ( K .LE. N )
202 IF( IPIV( K ).GT.0 ) THEN
203 * 1 x 1 diagonal block
204 * Interchange rows K and IPIV(K).
205 KP = IPIV( K )
206 IF( KP.NE.K )
207 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
208 K=K+1
209 ELSE
210 * 2 x 2 diagonal block
211 * Interchange rows K and -IPIV(K+1).
212 KP = -IPIV( K+1 )
213 IF( KP.EQ.-IPIV( K ) )
214 $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
215 K=K+2
216 ENDIF
217 END DO
218 *
219 * Compute (L \P**T * B) -> B [ (L \P**T * B) ]
220 *
221 CALL ZTRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB)
222 *
223 * Compute D \ B -> B [ D \ (L \P**T * B) ]
224 *
225 I=1
226 DO WHILE ( I .LE. N )
227 IF( IPIV(I) .GT. 0 ) THEN
228 CALL ZSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
229 ELSE
230 AKM1K = WORK(I)
231 AKM1 = A( I, I ) / AKM1K
232 AK = A( I+1, I+1 ) / AKM1K
233 DENOM = AKM1*AK - ONE
234 DO 25 J = 1, NRHS
235 BKM1 = B( I, J ) / AKM1K
236 BK = B( I+1, J ) / AKM1K
237 B( I, J ) = ( AK*BKM1-BK ) / DENOM
238 B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
239 25 CONTINUE
240 I = I + 1
241 ENDIF
242 I = I + 1
243 END DO
244 *
245 * Compute (L**T \ B) -> B [ L**T \ (D \ (L \P**T * B) ) ]
246 *
247 CALL ZTRSM('L','L','T','U',N,NRHS,ONE,A,LDA,B,LDB)
248 *
249 * P * B [ P * (L**T \ (D \ (L \P**T * B) )) ]
250 *
251 K=N
252 DO WHILE ( K .GE. 1 )
253 IF( IPIV( K ).GT.0 ) THEN
254 * 1 x 1 diagonal block
255 * Interchange rows K and IPIV(K).
256 KP = IPIV( K )
257 IF( KP.NE.K )
258 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
259 K=K-1
260 ELSE
261 * 2 x 2 diagonal block
262 * Interchange rows K-1 and -IPIV(K).
263 KP = -IPIV( K )
264 IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
265 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
266 K=K-2
267 ENDIF
268 END DO
269 *
270 END IF
271 *
272 * Revert A
273 *
274 CALL ZSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
275 *
276 RETURN
277 *
278 * End of ZSYTRS2
279 *
280 END