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:
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsytri2x.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack