1 SUBROUTINE DSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * -- Written by Julie Langou of the Univ. of TN --
9 *
10 * .. Scalar Arguments ..
11 CHARACTER UPLO
12 INTEGER INFO, LDA, N, NB
13 * ..
14 * .. Array Arguments ..
15 INTEGER IPIV( * )
16 DOUBLE PRECISION A( LDA, * ), WORK( N+NB+1,* )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * DSYTRI2X computes the inverse of a real symmetric indefinite matrix
23 * A using the factorization A = U*D*U**T or A = L*D*L**T computed by
24 * DSYTRF.
25 *
26 * Arguments
27 * =========
28 *
29 * UPLO (input) CHARACTER*1
30 * Specifies whether the details of the factorization are stored
31 * as an upper or lower triangular matrix.
32 * = 'U': Upper triangular, form is A = U*D*U**T;
33 * = 'L': Lower triangular, form is A = L*D*L**T.
34 *
35 * N (input) INTEGER
36 * The order of the matrix A. N >= 0.
37 *
38 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
39 * On entry, the NNB diagonal matrix D and the multipliers
40 * used to obtain the factor U or L as computed by DSYTRF.
41 *
42 * On exit, if INFO = 0, the (symmetric) inverse of the original
43 * matrix. If UPLO = 'U', the upper triangular part of the
44 * inverse is formed and the part of A below the diagonal is not
45 * referenced; if UPLO = 'L' the lower triangular part of the
46 * inverse is formed and the part of A above the diagonal is
47 * not referenced.
48 *
49 * LDA (input) INTEGER
50 * The leading dimension of the array A. LDA >= max(1,N).
51 *
52 * IPIV (input) INTEGER array, dimension (N)
53 * Details of the interchanges and the NNB structure of D
54 * as determined by DSYTRF.
55 *
56 * WORK (workspace) DOUBLE PRECISION array, dimension (N+NNB+1,NNB+3)
57 *
58 * NB (input) INTEGER
59 * Block size
60 *
61 * INFO (output) INTEGER
62 * = 0: successful exit
63 * < 0: if INFO = -i, the i-th argument had an illegal value
64 * > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
65 * inverse could not be computed.
66 *
67 * =====================================================================
68 *
69 * .. Parameters ..
70 DOUBLE PRECISION ONE, ZERO
71 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
72 * ..
73 * .. Local Scalars ..
74 LOGICAL UPPER
75 INTEGER I, IINFO, IP, K, CUT, NNB
76 INTEGER COUNT
77 INTEGER J, U11, INVD
78
79 DOUBLE PRECISION AK, AKKP1, AKP1, D, T
80 DOUBLE PRECISION U01_I_J, U01_IP1_J
81 DOUBLE PRECISION U11_I_J, U11_IP1_J
82 * ..
83 * .. External Functions ..
84 LOGICAL LSAME
85 EXTERNAL LSAME
86 * ..
87 * .. External Subroutines ..
88 EXTERNAL DSYCONV, XERBLA, DTRTRI
89 EXTERNAL DGEMM, DTRMM, DSYSWAPR
90 * ..
91 * .. Intrinsic Functions ..
92 INTRINSIC MAX
93 * ..
94 * .. Executable Statements ..
95 *
96 * Test the input parameters.
97 *
98 INFO = 0
99 UPPER = LSAME( UPLO, 'U' )
100 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
101 INFO = -1
102 ELSE IF( N.LT.0 ) THEN
103 INFO = -2
104 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
105 INFO = -4
106 END IF
107 *
108 * Quick return if possible
109 *
110 *
111 IF( INFO.NE.0 ) THEN
112 CALL XERBLA( 'DSYTRI2X', -INFO )
113 RETURN
114 END IF
115 IF( N.EQ.0 )
116 $ RETURN
117 *
118 * Convert A
119 * Workspace got Non-diag elements of D
120 *
121 CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
122 *
123 * Check that the diagonal matrix D is nonsingular.
124 *
125 IF( UPPER ) THEN
126 *
127 * Upper triangular storage: examine D from bottom to top
128 *
129 DO INFO = N, 1, -1
130 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
131 $ RETURN
132 END DO
133 ELSE
134 *
135 * Lower triangular storage: examine D from top to bottom.
136 *
137 DO INFO = 1, N
138 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
139 $ RETURN
140 END DO
141 END IF
142 INFO = 0
143 *
144 * Splitting Workspace
145 * U01 is a block (N,NB+1)
146 * The first element of U01 is in WORK(1,1)
147 * U11 is a block (NB+1,NB+1)
148 * The first element of U11 is in WORK(N+1,1)
149 U11 = N
150 * INVD is a block (N,2)
151 * The first element of INVD is in WORK(1,INVD)
152 INVD = NB+2
153
154 IF( UPPER ) THEN
155 *
156 * invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
157 *
158 CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO )
159 *
160 * inv(D) and inv(D)*inv(U)
161 *
162 K=1
163 DO WHILE ( K .LE. N )
164 IF( IPIV( K ).GT.0 ) THEN
165 * 1 x 1 diagonal NNB
166 WORK(K,INVD) = ONE / A( K, K )
167 WORK(K,INVD+1) = 0
168 K=K+1
169 ELSE
170 * 2 x 2 diagonal NNB
171 T = WORK(K+1,1)
172 AK = A( K, K ) / T
173 AKP1 = A( K+1, K+1 ) / T
174 AKKP1 = WORK(K+1,1) / T
175 D = T*( AK*AKP1-ONE )
176 WORK(K,INVD) = AKP1 / D
177 WORK(K+1,INVD+1) = AK / D
178 WORK(K,INVD+1) = -AKKP1 / D
179 WORK(K+1,INVD) = -AKKP1 / D
180 K=K+2
181 END IF
182 END DO
183 *
184 * inv(U**T) = (inv(U))**T
185 *
186 * inv(U**T)*inv(D)*inv(U)
187 *
188 CUT=N
189 DO WHILE (CUT .GT. 0)
190 NNB=NB
191 IF (CUT .LE. NNB) THEN
192 NNB=CUT
193 ELSE
194 COUNT = 0
195 * count negative elements,
196 DO I=CUT+1-NNB,CUT
197 IF (IPIV(I) .LT. 0) COUNT=COUNT+1
198 END DO
199 * need a even number for a clear cut
200 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
201 END IF
202
203 CUT=CUT-NNB
204 *
205 * U01 Block
206 *
207 DO I=1,CUT
208 DO J=1,NNB
209 WORK(I,J)=A(I,CUT+J)
210 END DO
211 END DO
212 *
213 * U11 Block
214 *
215 DO I=1,NNB
216 WORK(U11+I,I)=ONE
217 DO J=1,I-1
218 WORK(U11+I,J)=ZERO
219 END DO
220 DO J=I+1,NNB
221 WORK(U11+I,J)=A(CUT+I,CUT+J)
222 END DO
223 END DO
224 *
225 * invD*U01
226 *
227 I=1
228 DO WHILE (I .LE. CUT)
229 IF (IPIV(I) > 0) THEN
230 DO J=1,NNB
231 WORK(I,J)=WORK(I,INVD)*WORK(I,J)
232 END DO
233 I=I+1
234 ELSE
235 DO J=1,NNB
236 U01_I_J = WORK(I,J)
237 U01_IP1_J = WORK(I+1,J)
238 WORK(I,J)=WORK(I,INVD)*U01_I_J+
239 $ WORK(I,INVD+1)*U01_IP1_J
240 WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+
241 $ WORK(I+1,INVD+1)*U01_IP1_J
242 END DO
243 I=I+2
244 END IF
245 END DO
246 *
247 * invD1*U11
248 *
249 I=1
250 DO WHILE (I .LE. NNB)
251 IF (IPIV(CUT+I) > 0) THEN
252 DO J=I,NNB
253 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
254 END DO
255 I=I+1
256 ELSE
257 DO J=I,NNB
258 U11_I_J = WORK(U11+I,J)
259 U11_IP1_J = WORK(U11+I+1,J)
260 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
261 $ WORK(CUT+I,INVD+1)*WORK(U11+I+1,J)
262 WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+
263 $ WORK(CUT+I+1,INVD+1)*U11_IP1_J
264 END DO
265 I=I+2
266 END IF
267 END DO
268 *
269 * U11**T*invD1*U11->U11
270 *
271 CALL DTRMM('L','U','T','U',NNB, NNB,
272 $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
273 *
274 DO I=1,NNB
275 DO J=I,NNB
276 A(CUT+I,CUT+J)=WORK(U11+I,J)
277 END DO
278 END DO
279 *
280 * U01**T*invD*U01->A(CUT+I,CUT+J)
281 *
282 CALL DGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
283 $ WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
284
285 *
286 * U11 = U11**T*invD1*U11 + U01**T*invD*U01
287 *
288 DO I=1,NNB
289 DO J=I,NNB
290 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
291 END DO
292 END DO
293 *
294 * U01 = U00**T*invD0*U01
295 *
296 CALL DTRMM('L',UPLO,'T','U',CUT, NNB,
297 $ ONE,A,LDA,WORK,N+NB+1)
298
299 *
300 * Update U01
301 *
302 DO I=1,CUT
303 DO J=1,NNB
304 A(I,CUT+J)=WORK(I,J)
305 END DO
306 END DO
307 *
308 * Next Block
309 *
310 END DO
311 *
312 * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
313 *
314 I=1
315 DO WHILE ( I .LE. N )
316 IF( IPIV(I) .GT. 0 ) THEN
317 IP=IPIV(I)
318 IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP )
319 IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I )
320 ELSE
321 IP=-IPIV(I)
322 I=I+1
323 IF ( (I-1) .LT. IP)
324 $ CALL DSYSWAPR( UPLO, N, A, LDA, I-1 ,IP )
325 IF ( (I-1) .GT. IP)
326 $ CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I-1 )
327 ENDIF
328 I=I+1
329 END DO
330 ELSE
331 *
332 * LOWER...
333 *
334 * invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
335 *
336 CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO )
337 *
338 * inv(D) and inv(D)*inv(U)
339 *
340 K=N
341 DO WHILE ( K .GE. 1 )
342 IF( IPIV( K ).GT.0 ) THEN
343 * 1 x 1 diagonal NNB
344 WORK(K,INVD) = ONE / A( K, K )
345 WORK(K,INVD+1) = 0
346 K=K-1
347 ELSE
348 * 2 x 2 diagonal NNB
349 T = WORK(K-1,1)
350 AK = A( K-1, K-1 ) / T
351 AKP1 = A( K, K ) / T
352 AKKP1 = WORK(K-1,1) / T
353 D = T*( AK*AKP1-ONE )
354 WORK(K-1,INVD) = AKP1 / D
355 WORK(K,INVD) = AK / D
356 WORK(K,INVD+1) = -AKKP1 / D
357 WORK(K-1,INVD+1) = -AKKP1 / D
358 K=K-2
359 END IF
360 END DO
361 *
362 * inv(U**T) = (inv(U))**T
363 *
364 * inv(U**T)*inv(D)*inv(U)
365 *
366 CUT=0
367 DO WHILE (CUT .LT. N)
368 NNB=NB
369 IF (CUT + NNB .GT. N) THEN
370 NNB=N-CUT
371 ELSE
372 COUNT = 0
373 * count negative elements,
374 DO I=CUT+1,CUT+NNB
375 IF (IPIV(I) .LT. 0) COUNT=COUNT+1
376 END DO
377 * need a even number for a clear cut
378 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
379 END IF
380 * L21 Block
381 DO I=1,N-CUT-NNB
382 DO J=1,NNB
383 WORK(I,J)=A(CUT+NNB+I,CUT+J)
384 END DO
385 END DO
386 * L11 Block
387 DO I=1,NNB
388 WORK(U11+I,I)=ONE
389 DO J=I+1,NNB
390 WORK(U11+I,J)=ZERO
391 END DO
392 DO J=1,I-1
393 WORK(U11+I,J)=A(CUT+I,CUT+J)
394 END DO
395 END DO
396 *
397 * invD*L21
398 *
399 I=N-CUT-NNB
400 DO WHILE (I .GE. 1)
401 IF (IPIV(CUT+NNB+I) > 0) THEN
402 DO J=1,NNB
403 WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
404 END DO
405 I=I-1
406 ELSE
407 DO J=1,NNB
408 U01_I_J = WORK(I,J)
409 U01_IP1_J = WORK(I-1,J)
410 WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
411 $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
412 WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
413 $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
414 END DO
415 I=I-2
416 END IF
417 END DO
418 *
419 * invD1*L11
420 *
421 I=NNB
422 DO WHILE (I .GE. 1)
423 IF (IPIV(CUT+I) > 0) THEN
424 DO J=1,NNB
425 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
426 END DO
427 I=I-1
428 ELSE
429 DO J=1,NNB
430 U11_I_J = WORK(U11+I,J)
431 U11_IP1_J = WORK(U11+I-1,J)
432 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
433 $ WORK(CUT+I,INVD+1)*U11_IP1_J
434 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
435 $ WORK(CUT+I-1,INVD)*U11_IP1_J
436 END DO
437 I=I-2
438 END IF
439 END DO
440 *
441 * L11**T*invD1*L11->L11
442 *
443 CALL DTRMM('L',UPLO,'T','U',NNB, NNB,
444 $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
445
446 *
447 DO I=1,NNB
448 DO J=1,I
449 A(CUT+I,CUT+J)=WORK(U11+I,J)
450 END DO
451 END DO
452 *
453 IF ( (CUT+NNB) .LT. N ) THEN
454 *
455 * L21**T*invD2*L21->A(CUT+I,CUT+J)
456 *
457 CALL DGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
458 $ ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
459
460 *
461 * L11 = L11**T*invD1*L11 + U01**T*invD*U01
462 *
463 DO I=1,NNB
464 DO J=1,I
465 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
466 END DO
467 END DO
468 *
469 * L01 = L22**T*invD2*L21
470 *
471 CALL DTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB,
472 $ ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
473 *
474 * Update L21
475 *
476 DO I=1,N-CUT-NNB
477 DO J=1,NNB
478 A(CUT+NNB+I,CUT+J)=WORK(I,J)
479 END DO
480 END DO
481
482 ELSE
483 *
484 * L11 = L11**T*invD1*L11
485 *
486 DO I=1,NNB
487 DO J=1,I
488 A(CUT+I,CUT+J)=WORK(U11+I,J)
489 END DO
490 END DO
491 END IF
492 *
493 * Next Block
494 *
495 CUT=CUT+NNB
496 END DO
497 *
498 * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
499 *
500 I=N
501 DO WHILE ( I .GE. 1 )
502 IF( IPIV(I) .GT. 0 ) THEN
503 IP=IPIV(I)
504 IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP )
505 IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I )
506 ELSE
507 IP=-IPIV(I)
508 IF ( I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP )
509 IF ( I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP, I )
510 I=I-1
511 ENDIF
512 I=I-1
513 END DO
514 END IF
515 *
516 RETURN
517 *
518 * End of DSYTRI2X
519 *
520 END
521
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * -- Written by Julie Langou of the Univ. of TN --
9 *
10 * .. Scalar Arguments ..
11 CHARACTER UPLO
12 INTEGER INFO, LDA, N, NB
13 * ..
14 * .. Array Arguments ..
15 INTEGER IPIV( * )
16 DOUBLE PRECISION A( LDA, * ), WORK( N+NB+1,* )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * DSYTRI2X computes the inverse of a real symmetric indefinite matrix
23 * A using the factorization A = U*D*U**T or A = L*D*L**T computed by
24 * DSYTRF.
25 *
26 * Arguments
27 * =========
28 *
29 * UPLO (input) CHARACTER*1
30 * Specifies whether the details of the factorization are stored
31 * as an upper or lower triangular matrix.
32 * = 'U': Upper triangular, form is A = U*D*U**T;
33 * = 'L': Lower triangular, form is A = L*D*L**T.
34 *
35 * N (input) INTEGER
36 * The order of the matrix A. N >= 0.
37 *
38 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
39 * On entry, the NNB diagonal matrix D and the multipliers
40 * used to obtain the factor U or L as computed by DSYTRF.
41 *
42 * On exit, if INFO = 0, the (symmetric) inverse of the original
43 * matrix. If UPLO = 'U', the upper triangular part of the
44 * inverse is formed and the part of A below the diagonal is not
45 * referenced; if UPLO = 'L' the lower triangular part of the
46 * inverse is formed and the part of A above the diagonal is
47 * not referenced.
48 *
49 * LDA (input) INTEGER
50 * The leading dimension of the array A. LDA >= max(1,N).
51 *
52 * IPIV (input) INTEGER array, dimension (N)
53 * Details of the interchanges and the NNB structure of D
54 * as determined by DSYTRF.
55 *
56 * WORK (workspace) DOUBLE PRECISION array, dimension (N+NNB+1,NNB+3)
57 *
58 * NB (input) INTEGER
59 * Block size
60 *
61 * INFO (output) INTEGER
62 * = 0: successful exit
63 * < 0: if INFO = -i, the i-th argument had an illegal value
64 * > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
65 * inverse could not be computed.
66 *
67 * =====================================================================
68 *
69 * .. Parameters ..
70 DOUBLE PRECISION ONE, ZERO
71 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
72 * ..
73 * .. Local Scalars ..
74 LOGICAL UPPER
75 INTEGER I, IINFO, IP, K, CUT, NNB
76 INTEGER COUNT
77 INTEGER J, U11, INVD
78
79 DOUBLE PRECISION AK, AKKP1, AKP1, D, T
80 DOUBLE PRECISION U01_I_J, U01_IP1_J
81 DOUBLE PRECISION U11_I_J, U11_IP1_J
82 * ..
83 * .. External Functions ..
84 LOGICAL LSAME
85 EXTERNAL LSAME
86 * ..
87 * .. External Subroutines ..
88 EXTERNAL DSYCONV, XERBLA, DTRTRI
89 EXTERNAL DGEMM, DTRMM, DSYSWAPR
90 * ..
91 * .. Intrinsic Functions ..
92 INTRINSIC MAX
93 * ..
94 * .. Executable Statements ..
95 *
96 * Test the input parameters.
97 *
98 INFO = 0
99 UPPER = LSAME( UPLO, 'U' )
100 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
101 INFO = -1
102 ELSE IF( N.LT.0 ) THEN
103 INFO = -2
104 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
105 INFO = -4
106 END IF
107 *
108 * Quick return if possible
109 *
110 *
111 IF( INFO.NE.0 ) THEN
112 CALL XERBLA( 'DSYTRI2X', -INFO )
113 RETURN
114 END IF
115 IF( N.EQ.0 )
116 $ RETURN
117 *
118 * Convert A
119 * Workspace got Non-diag elements of D
120 *
121 CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
122 *
123 * Check that the diagonal matrix D is nonsingular.
124 *
125 IF( UPPER ) THEN
126 *
127 * Upper triangular storage: examine D from bottom to top
128 *
129 DO INFO = N, 1, -1
130 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
131 $ RETURN
132 END DO
133 ELSE
134 *
135 * Lower triangular storage: examine D from top to bottom.
136 *
137 DO INFO = 1, N
138 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
139 $ RETURN
140 END DO
141 END IF
142 INFO = 0
143 *
144 * Splitting Workspace
145 * U01 is a block (N,NB+1)
146 * The first element of U01 is in WORK(1,1)
147 * U11 is a block (NB+1,NB+1)
148 * The first element of U11 is in WORK(N+1,1)
149 U11 = N
150 * INVD is a block (N,2)
151 * The first element of INVD is in WORK(1,INVD)
152 INVD = NB+2
153
154 IF( UPPER ) THEN
155 *
156 * invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
157 *
158 CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO )
159 *
160 * inv(D) and inv(D)*inv(U)
161 *
162 K=1
163 DO WHILE ( K .LE. N )
164 IF( IPIV( K ).GT.0 ) THEN
165 * 1 x 1 diagonal NNB
166 WORK(K,INVD) = ONE / A( K, K )
167 WORK(K,INVD+1) = 0
168 K=K+1
169 ELSE
170 * 2 x 2 diagonal NNB
171 T = WORK(K+1,1)
172 AK = A( K, K ) / T
173 AKP1 = A( K+1, K+1 ) / T
174 AKKP1 = WORK(K+1,1) / T
175 D = T*( AK*AKP1-ONE )
176 WORK(K,INVD) = AKP1 / D
177 WORK(K+1,INVD+1) = AK / D
178 WORK(K,INVD+1) = -AKKP1 / D
179 WORK(K+1,INVD) = -AKKP1 / D
180 K=K+2
181 END IF
182 END DO
183 *
184 * inv(U**T) = (inv(U))**T
185 *
186 * inv(U**T)*inv(D)*inv(U)
187 *
188 CUT=N
189 DO WHILE (CUT .GT. 0)
190 NNB=NB
191 IF (CUT .LE. NNB) THEN
192 NNB=CUT
193 ELSE
194 COUNT = 0
195 * count negative elements,
196 DO I=CUT+1-NNB,CUT
197 IF (IPIV(I) .LT. 0) COUNT=COUNT+1
198 END DO
199 * need a even number for a clear cut
200 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
201 END IF
202
203 CUT=CUT-NNB
204 *
205 * U01 Block
206 *
207 DO I=1,CUT
208 DO J=1,NNB
209 WORK(I,J)=A(I,CUT+J)
210 END DO
211 END DO
212 *
213 * U11 Block
214 *
215 DO I=1,NNB
216 WORK(U11+I,I)=ONE
217 DO J=1,I-1
218 WORK(U11+I,J)=ZERO
219 END DO
220 DO J=I+1,NNB
221 WORK(U11+I,J)=A(CUT+I,CUT+J)
222 END DO
223 END DO
224 *
225 * invD*U01
226 *
227 I=1
228 DO WHILE (I .LE. CUT)
229 IF (IPIV(I) > 0) THEN
230 DO J=1,NNB
231 WORK(I,J)=WORK(I,INVD)*WORK(I,J)
232 END DO
233 I=I+1
234 ELSE
235 DO J=1,NNB
236 U01_I_J = WORK(I,J)
237 U01_IP1_J = WORK(I+1,J)
238 WORK(I,J)=WORK(I,INVD)*U01_I_J+
239 $ WORK(I,INVD+1)*U01_IP1_J
240 WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+
241 $ WORK(I+1,INVD+1)*U01_IP1_J
242 END DO
243 I=I+2
244 END IF
245 END DO
246 *
247 * invD1*U11
248 *
249 I=1
250 DO WHILE (I .LE. NNB)
251 IF (IPIV(CUT+I) > 0) THEN
252 DO J=I,NNB
253 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
254 END DO
255 I=I+1
256 ELSE
257 DO J=I,NNB
258 U11_I_J = WORK(U11+I,J)
259 U11_IP1_J = WORK(U11+I+1,J)
260 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
261 $ WORK(CUT+I,INVD+1)*WORK(U11+I+1,J)
262 WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+
263 $ WORK(CUT+I+1,INVD+1)*U11_IP1_J
264 END DO
265 I=I+2
266 END IF
267 END DO
268 *
269 * U11**T*invD1*U11->U11
270 *
271 CALL DTRMM('L','U','T','U',NNB, NNB,
272 $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
273 *
274 DO I=1,NNB
275 DO J=I,NNB
276 A(CUT+I,CUT+J)=WORK(U11+I,J)
277 END DO
278 END DO
279 *
280 * U01**T*invD*U01->A(CUT+I,CUT+J)
281 *
282 CALL DGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
283 $ WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
284
285 *
286 * U11 = U11**T*invD1*U11 + U01**T*invD*U01
287 *
288 DO I=1,NNB
289 DO J=I,NNB
290 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
291 END DO
292 END DO
293 *
294 * U01 = U00**T*invD0*U01
295 *
296 CALL DTRMM('L',UPLO,'T','U',CUT, NNB,
297 $ ONE,A,LDA,WORK,N+NB+1)
298
299 *
300 * Update U01
301 *
302 DO I=1,CUT
303 DO J=1,NNB
304 A(I,CUT+J)=WORK(I,J)
305 END DO
306 END DO
307 *
308 * Next Block
309 *
310 END DO
311 *
312 * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
313 *
314 I=1
315 DO WHILE ( I .LE. N )
316 IF( IPIV(I) .GT. 0 ) THEN
317 IP=IPIV(I)
318 IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP )
319 IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I )
320 ELSE
321 IP=-IPIV(I)
322 I=I+1
323 IF ( (I-1) .LT. IP)
324 $ CALL DSYSWAPR( UPLO, N, A, LDA, I-1 ,IP )
325 IF ( (I-1) .GT. IP)
326 $ CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I-1 )
327 ENDIF
328 I=I+1
329 END DO
330 ELSE
331 *
332 * LOWER...
333 *
334 * invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
335 *
336 CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO )
337 *
338 * inv(D) and inv(D)*inv(U)
339 *
340 K=N
341 DO WHILE ( K .GE. 1 )
342 IF( IPIV( K ).GT.0 ) THEN
343 * 1 x 1 diagonal NNB
344 WORK(K,INVD) = ONE / A( K, K )
345 WORK(K,INVD+1) = 0
346 K=K-1
347 ELSE
348 * 2 x 2 diagonal NNB
349 T = WORK(K-1,1)
350 AK = A( K-1, K-1 ) / T
351 AKP1 = A( K, K ) / T
352 AKKP1 = WORK(K-1,1) / T
353 D = T*( AK*AKP1-ONE )
354 WORK(K-1,INVD) = AKP1 / D
355 WORK(K,INVD) = AK / D
356 WORK(K,INVD+1) = -AKKP1 / D
357 WORK(K-1,INVD+1) = -AKKP1 / D
358 K=K-2
359 END IF
360 END DO
361 *
362 * inv(U**T) = (inv(U))**T
363 *
364 * inv(U**T)*inv(D)*inv(U)
365 *
366 CUT=0
367 DO WHILE (CUT .LT. N)
368 NNB=NB
369 IF (CUT + NNB .GT. N) THEN
370 NNB=N-CUT
371 ELSE
372 COUNT = 0
373 * count negative elements,
374 DO I=CUT+1,CUT+NNB
375 IF (IPIV(I) .LT. 0) COUNT=COUNT+1
376 END DO
377 * need a even number for a clear cut
378 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
379 END IF
380 * L21 Block
381 DO I=1,N-CUT-NNB
382 DO J=1,NNB
383 WORK(I,J)=A(CUT+NNB+I,CUT+J)
384 END DO
385 END DO
386 * L11 Block
387 DO I=1,NNB
388 WORK(U11+I,I)=ONE
389 DO J=I+1,NNB
390 WORK(U11+I,J)=ZERO
391 END DO
392 DO J=1,I-1
393 WORK(U11+I,J)=A(CUT+I,CUT+J)
394 END DO
395 END DO
396 *
397 * invD*L21
398 *
399 I=N-CUT-NNB
400 DO WHILE (I .GE. 1)
401 IF (IPIV(CUT+NNB+I) > 0) THEN
402 DO J=1,NNB
403 WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
404 END DO
405 I=I-1
406 ELSE
407 DO J=1,NNB
408 U01_I_J = WORK(I,J)
409 U01_IP1_J = WORK(I-1,J)
410 WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
411 $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
412 WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
413 $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
414 END DO
415 I=I-2
416 END IF
417 END DO
418 *
419 * invD1*L11
420 *
421 I=NNB
422 DO WHILE (I .GE. 1)
423 IF (IPIV(CUT+I) > 0) THEN
424 DO J=1,NNB
425 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
426 END DO
427 I=I-1
428 ELSE
429 DO J=1,NNB
430 U11_I_J = WORK(U11+I,J)
431 U11_IP1_J = WORK(U11+I-1,J)
432 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
433 $ WORK(CUT+I,INVD+1)*U11_IP1_J
434 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
435 $ WORK(CUT+I-1,INVD)*U11_IP1_J
436 END DO
437 I=I-2
438 END IF
439 END DO
440 *
441 * L11**T*invD1*L11->L11
442 *
443 CALL DTRMM('L',UPLO,'T','U',NNB, NNB,
444 $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
445
446 *
447 DO I=1,NNB
448 DO J=1,I
449 A(CUT+I,CUT+J)=WORK(U11+I,J)
450 END DO
451 END DO
452 *
453 IF ( (CUT+NNB) .LT. N ) THEN
454 *
455 * L21**T*invD2*L21->A(CUT+I,CUT+J)
456 *
457 CALL DGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
458 $ ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
459
460 *
461 * L11 = L11**T*invD1*L11 + U01**T*invD*U01
462 *
463 DO I=1,NNB
464 DO J=1,I
465 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
466 END DO
467 END DO
468 *
469 * L01 = L22**T*invD2*L21
470 *
471 CALL DTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB,
472 $ ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
473 *
474 * Update L21
475 *
476 DO I=1,N-CUT-NNB
477 DO J=1,NNB
478 A(CUT+NNB+I,CUT+J)=WORK(I,J)
479 END DO
480 END DO
481
482 ELSE
483 *
484 * L11 = L11**T*invD1*L11
485 *
486 DO I=1,NNB
487 DO J=1,I
488 A(CUT+I,CUT+J)=WORK(U11+I,J)
489 END DO
490 END DO
491 END IF
492 *
493 * Next Block
494 *
495 CUT=CUT+NNB
496 END DO
497 *
498 * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
499 *
500 I=N
501 DO WHILE ( I .GE. 1 )
502 IF( IPIV(I) .GT. 0 ) THEN
503 IP=IPIV(I)
504 IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP )
505 IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I )
506 ELSE
507 IP=-IPIV(I)
508 IF ( I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP )
509 IF ( I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP, I )
510 I=I-1
511 ENDIF
512 I=I-1
513 END DO
514 END IF
515 *
516 RETURN
517 *
518 * End of DSYTRI2X
519 *
520 END
521