![[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|
Return to dsytri2x.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
2: *
3: * -- LAPACK routine (version 3.3.0) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2010
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')*inv(D)*inv(U)*P'.
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') = (inv(U))'
185: *
186: * inv(U')*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: * U11T*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: * U01'invD*U01->A(CUT+I,CUT+J)
275: *
276: CALL DGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
277: $ WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
278: *
279: * U11 = U11T*invD1*U11 + U01'invD*U01
280: *
281: DO I=1,NNB
282: DO J=I,NNB
283: A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
284: END DO
285: END DO
286: *
287: * U01 = U00T*invD0*U01
288: *
289: CALL DTRMM('L',UPLO,'T','U',CUT, NNB,
290: $ ONE,A,LDA,WORK,N+NB+1)
291:
292: *
293: * Update U01
294: *
295: DO I=1,CUT
296: DO J=1,NNB
297: A(I,CUT+J)=WORK(I,J)
298: END DO
299: END DO
300: *
301: * Next Block
302: *
303: END DO
304: *
305: * Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P'
306: *
307: I=1
308: DO WHILE ( I .LE. N )
309: IF( IPIV(I) .GT. 0 ) THEN
310: IP=IPIV(I)
311: IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, I ,IP )
312: IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, IP ,I )
313: ELSE
314: IP=-IPIV(I)
315: I=I+1
316: IF ( (I-1) .LT. IP)
317: $ CALL DSYSWAPR( UPLO, N, A, I-1 ,IP )
318: IF ( (I-1) .GT. IP)
319: $ CALL DSYSWAPR( UPLO, N, A, IP ,I-1 )
320: ENDIF
321: I=I+1
322: END DO
323: ELSE
324: *
325: * LOWER...
326: *
327: * invA = P * inv(U')*inv(D)*inv(U)*P'.
328: *
329: CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO )
330: *
331: * inv(D) and inv(D)*inv(U)
332: *
333: K=N
334: DO WHILE ( K .GE. 1 )
335: IF( IPIV( K ).GT.0 ) THEN
336: * 1 x 1 diagonal NNB
337: WORK(K,INVD) = ONE / A( K, K )
338: WORK(K,INVD+1) = 0
339: K=K-1
340: ELSE
341: * 2 x 2 diagonal NNB
342: T = WORK(K-1,1)
343: AK = A( K-1, K-1 ) / T
344: AKP1 = A( K, K ) / T
345: AKKP1 = WORK(K-1,1) / T
346: D = T*( AK*AKP1-ONE )
347: WORK(K-1,INVD) = AKP1 / D
348: WORK(K,INVD) = AK / D
349: WORK(K,INVD+1) = -AKKP1 / D
350: WORK(K-1,INVD+1) = -AKKP1 / D
351: K=K-2
352: END IF
353: END DO
354: *
355: * inv(U') = (inv(U))'
356: *
357: * inv(U')*inv(D)*inv(U)
358: *
359: CUT=0
360: DO WHILE (CUT .LT. N)
361: NNB=NB
362: IF (CUT + NNB .GT. N) THEN
363: NNB=N-CUT
364: ELSE
365: COUNT = 0
366: * count negative elements,
367: DO I=CUT+1,CUT+NNB
368: IF (IPIV(I) .LT. 0) COUNT=COUNT+1
369: END DO
370: * need a even number for a clear cut
371: IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
372: END IF
373: * L21 Block
374: DO I=1,N-CUT-NNB
375: DO J=1,NNB
376: WORK(I,J)=A(CUT+NNB+I,CUT+J)
377: END DO
378: END DO
379: * L11 Block
380: DO I=1,NNB
381: WORK(U11+I,I)=ONE
382: DO J=I+1,NNB
383: WORK(U11+I,J)=ZERO
384: END DO
385: DO J=1,I-1
386: WORK(U11+I,J)=A(CUT+I,CUT+J)
387: END DO
388: END DO
389: *
390: * invD*L21
391: *
392: I=N-CUT-NNB
393: DO WHILE (I .GE. 1)
394: IF (IPIV(CUT+NNB+I) > 0) THEN
395: DO J=1,NNB
396: WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
397: END DO
398: I=I-1
399: ELSE
400: DO J=1,NNB
401: U01_I_J = WORK(I,J)
402: U01_IP1_J = WORK(I-1,J)
403: WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
404: $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
405: WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
406: $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
407: END DO
408: I=I-2
409: END IF
410: END DO
411: *
412: * invD1*L11
413: *
414: I=NNB
415: DO WHILE (I .GE. 1)
416: IF (IPIV(CUT+I) > 0) THEN
417: DO J=1,NNB
418: WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
419: END DO
420: I=I-1
421: ELSE
422: DO J=1,NNB
423: U11_I_J = WORK(U11+I,J)
424: U11_IP1_J = WORK(U11+I-1,J)
425: WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
426: $ WORK(CUT+I,INVD+1)*U11_IP1_J
427: WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
428: $ WORK(CUT+I-1,INVD)*U11_IP1_J
429: END DO
430: I=I-2
431: END IF
432: END DO
433: *
434: * L11T*invD1*L11->L11
435: *
436: CALL DTRMM('L',UPLO,'T','U',NNB, NNB,
437: $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
438:
439: IF ( (CUT+NNB) .LT. N ) THEN
440: *
441: * L21T*invD2*L21->A(CUT+I,CUT+J)
442: *
443: CALL DGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
444: $ ,LDA,WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
445:
446: *
447: * L11 = L11T*invD1*L11 + U01'invD*U01
448: *
449: DO I=1,NNB
450: DO J=1,I
451: A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
452: END DO
453: END DO
454: *
455: * L01 = L22T*invD2*L21
456: *
457: CALL DTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB,
458: $ ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
459: *
460: * Update L21
461: *
462: DO I=1,N-CUT-NNB
463: DO J=1,NNB
464: A(CUT+NNB+I,CUT+J)=WORK(I,J)
465: END DO
466: END DO
467:
468: ELSE
469: *
470: * L11 = L11T*invD1*L11
471: *
472: DO I=1,NNB
473: DO J=1,I
474: A(CUT+I,CUT+J)=WORK(U11+I,J)
475: END DO
476: END DO
477: END IF
478: *
479: * Next Block
480: *
481: CUT=CUT+NNB
482: END DO
483: *
484: * Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P'
485: *
486: I=N
487: DO WHILE ( I .GE. 1 )
488: IF( IPIV(I) .GT. 0 ) THEN
489: IP=IPIV(I)
490: IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, I ,IP )
491: IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, IP ,I )
492: ELSE
493: IP=-IPIV(I)
494: IF ( I .LT. IP) CALL DSYSWAPR( UPLO, N, A, I ,IP )
495: IF ( I .GT. IP) CALL DSYSWAPR( UPLO, N, A, IP ,I )
496: I=I-1
497: ENDIF
498: I=I-1
499: END DO
500: END IF
501: *
502: RETURN
503: *
504: * End of DSYTRI2X
505: *
506: END
507: