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