Annotation of rpl/lapack/lapack/dsytri2x.f, revision 1.4
1.4 ! bertrand 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+NNB+1,NNB+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: *> \date November 2011
! 117: *
! 118: *> \ingroup doubleSYcomputational
! 119: *
! 120: * =====================================================================
1.1 bertrand 121: SUBROUTINE DSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
122: *
1.4 ! bertrand 123: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 124: * -- LAPACK is a software package provided by Univ. of Tennessee, --
125: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.4 ! bertrand 126: * November 2011
1.1 bertrand 127: *
128: * .. Scalar Arguments ..
129: CHARACTER UPLO
130: INTEGER INFO, LDA, N, NB
131: * ..
132: * .. Array Arguments ..
133: INTEGER IPIV( * )
134: DOUBLE PRECISION A( LDA, * ), WORK( N+NB+1,* )
135: * ..
136: *
137: * =====================================================================
138: *
139: * .. Parameters ..
140: DOUBLE PRECISION ONE, ZERO
141: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
142: * ..
143: * .. Local Scalars ..
144: LOGICAL UPPER
145: INTEGER I, IINFO, IP, K, CUT, NNB
146: INTEGER COUNT
147: INTEGER J, U11, INVD
148:
149: DOUBLE PRECISION AK, AKKP1, AKP1, D, T
150: DOUBLE PRECISION U01_I_J, U01_IP1_J
151: DOUBLE PRECISION U11_I_J, U11_IP1_J
152: * ..
153: * .. External Functions ..
154: LOGICAL LSAME
155: EXTERNAL LSAME
156: * ..
157: * .. External Subroutines ..
158: EXTERNAL DSYCONV, XERBLA, DTRTRI
159: EXTERNAL DGEMM, DTRMM, DSYSWAPR
160: * ..
161: * .. Intrinsic Functions ..
162: INTRINSIC MAX
163: * ..
164: * .. Executable Statements ..
165: *
166: * Test the input parameters.
167: *
168: INFO = 0
169: UPPER = LSAME( UPLO, 'U' )
170: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
171: INFO = -1
172: ELSE IF( N.LT.0 ) THEN
173: INFO = -2
174: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
175: INFO = -4
176: END IF
177: *
178: * Quick return if possible
179: *
180: *
181: IF( INFO.NE.0 ) THEN
182: CALL XERBLA( 'DSYTRI2X', -INFO )
183: RETURN
184: END IF
185: IF( N.EQ.0 )
186: $ RETURN
187: *
188: * Convert A
189: * Workspace got Non-diag elements of D
190: *
191: CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
192: *
193: * Check that the diagonal matrix D is nonsingular.
194: *
195: IF( UPPER ) THEN
196: *
197: * Upper triangular storage: examine D from bottom to top
198: *
199: DO INFO = N, 1, -1
200: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
201: $ RETURN
202: END DO
203: ELSE
204: *
205: * Lower triangular storage: examine D from top to bottom.
206: *
207: DO INFO = 1, N
208: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
209: $ RETURN
210: END DO
211: END IF
212: INFO = 0
213: *
214: * Splitting Workspace
215: * U01 is a block (N,NB+1)
216: * The first element of U01 is in WORK(1,1)
217: * U11 is a block (NB+1,NB+1)
218: * The first element of U11 is in WORK(N+1,1)
219: U11 = N
220: * INVD is a block (N,2)
221: * The first element of INVD is in WORK(1,INVD)
222: INVD = NB+2
223:
224: IF( UPPER ) THEN
225: *
1.3 bertrand 226: * invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
1.1 bertrand 227: *
228: CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO )
229: *
230: * inv(D) and inv(D)*inv(U)
231: *
232: K=1
233: DO WHILE ( K .LE. N )
234: IF( IPIV( K ).GT.0 ) THEN
235: * 1 x 1 diagonal NNB
236: WORK(K,INVD) = ONE / A( K, K )
237: WORK(K,INVD+1) = 0
238: K=K+1
239: ELSE
240: * 2 x 2 diagonal NNB
241: T = WORK(K+1,1)
242: AK = A( K, K ) / T
243: AKP1 = A( K+1, K+1 ) / T
244: AKKP1 = WORK(K+1,1) / T
245: D = T*( AK*AKP1-ONE )
246: WORK(K,INVD) = AKP1 / D
247: WORK(K+1,INVD+1) = AK / D
248: WORK(K,INVD+1) = -AKKP1 / D
249: WORK(K+1,INVD) = -AKKP1 / D
250: K=K+2
251: END IF
252: END DO
253: *
1.3 bertrand 254: * inv(U**T) = (inv(U))**T
1.1 bertrand 255: *
1.3 bertrand 256: * inv(U**T)*inv(D)*inv(U)
1.1 bertrand 257: *
258: CUT=N
259: DO WHILE (CUT .GT. 0)
260: NNB=NB
261: IF (CUT .LE. NNB) THEN
262: NNB=CUT
263: ELSE
264: COUNT = 0
265: * count negative elements,
266: DO I=CUT+1-NNB,CUT
267: IF (IPIV(I) .LT. 0) COUNT=COUNT+1
268: END DO
269: * need a even number for a clear cut
270: IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
271: END IF
272:
273: CUT=CUT-NNB
274: *
275: * U01 Block
276: *
277: DO I=1,CUT
278: DO J=1,NNB
279: WORK(I,J)=A(I,CUT+J)
280: END DO
281: END DO
282: *
283: * U11 Block
284: *
285: DO I=1,NNB
286: WORK(U11+I,I)=ONE
287: DO J=1,I-1
288: WORK(U11+I,J)=ZERO
289: END DO
290: DO J=I+1,NNB
291: WORK(U11+I,J)=A(CUT+I,CUT+J)
292: END DO
293: END DO
294: *
295: * invD*U01
296: *
297: I=1
298: DO WHILE (I .LE. CUT)
299: IF (IPIV(I) > 0) THEN
300: DO J=1,NNB
301: WORK(I,J)=WORK(I,INVD)*WORK(I,J)
302: END DO
303: I=I+1
304: ELSE
305: DO J=1,NNB
306: U01_I_J = WORK(I,J)
307: U01_IP1_J = WORK(I+1,J)
308: WORK(I,J)=WORK(I,INVD)*U01_I_J+
309: $ WORK(I,INVD+1)*U01_IP1_J
310: WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+
311: $ WORK(I+1,INVD+1)*U01_IP1_J
312: END DO
313: I=I+2
314: END IF
315: END DO
316: *
317: * invD1*U11
318: *
319: I=1
320: DO WHILE (I .LE. NNB)
321: IF (IPIV(CUT+I) > 0) THEN
322: DO J=I,NNB
323: WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
324: END DO
325: I=I+1
326: ELSE
327: DO J=I,NNB
328: U11_I_J = WORK(U11+I,J)
329: U11_IP1_J = WORK(U11+I+1,J)
330: WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
331: $ WORK(CUT+I,INVD+1)*WORK(U11+I+1,J)
332: WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+
333: $ WORK(CUT+I+1,INVD+1)*U11_IP1_J
334: END DO
335: I=I+2
336: END IF
337: END DO
338: *
1.3 bertrand 339: * U11**T*invD1*U11->U11
1.1 bertrand 340: *
341: CALL DTRMM('L','U','T','U',NNB, NNB,
342: $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
343: *
1.3 bertrand 344: DO I=1,NNB
345: DO J=I,NNB
346: A(CUT+I,CUT+J)=WORK(U11+I,J)
347: END DO
348: END DO
349: *
350: * U01**T*invD*U01->A(CUT+I,CUT+J)
1.1 bertrand 351: *
352: CALL DGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
1.3 bertrand 353: $ WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
354:
1.1 bertrand 355: *
1.3 bertrand 356: * U11 = U11**T*invD1*U11 + U01**T*invD*U01
1.1 bertrand 357: *
358: DO I=1,NNB
359: DO J=I,NNB
360: A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
361: END DO
362: END DO
363: *
1.3 bertrand 364: * U01 = U00**T*invD0*U01
1.1 bertrand 365: *
366: CALL DTRMM('L',UPLO,'T','U',CUT, NNB,
367: $ ONE,A,LDA,WORK,N+NB+1)
368:
369: *
370: * Update U01
371: *
372: DO I=1,CUT
373: DO J=1,NNB
374: A(I,CUT+J)=WORK(I,J)
375: END DO
376: END DO
377: *
378: * Next Block
379: *
380: END DO
381: *
1.3 bertrand 382: * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
1.1 bertrand 383: *
384: I=1
385: DO WHILE ( I .LE. N )
386: IF( IPIV(I) .GT. 0 ) THEN
387: IP=IPIV(I)
1.3 bertrand 388: IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP )
389: IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I )
1.1 bertrand 390: ELSE
391: IP=-IPIV(I)
392: I=I+1
393: IF ( (I-1) .LT. IP)
1.3 bertrand 394: $ CALL DSYSWAPR( UPLO, N, A, LDA, I-1 ,IP )
1.1 bertrand 395: IF ( (I-1) .GT. IP)
1.3 bertrand 396: $ CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I-1 )
1.1 bertrand 397: ENDIF
398: I=I+1
399: END DO
400: ELSE
401: *
402: * LOWER...
403: *
1.3 bertrand 404: * invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
1.1 bertrand 405: *
406: CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO )
407: *
408: * inv(D) and inv(D)*inv(U)
409: *
410: K=N
411: DO WHILE ( K .GE. 1 )
412: IF( IPIV( K ).GT.0 ) THEN
413: * 1 x 1 diagonal NNB
414: WORK(K,INVD) = ONE / A( K, K )
415: WORK(K,INVD+1) = 0
416: K=K-1
417: ELSE
418: * 2 x 2 diagonal NNB
419: T = WORK(K-1,1)
420: AK = A( K-1, K-1 ) / T
421: AKP1 = A( K, K ) / T
422: AKKP1 = WORK(K-1,1) / T
423: D = T*( AK*AKP1-ONE )
424: WORK(K-1,INVD) = AKP1 / D
425: WORK(K,INVD) = AK / D
426: WORK(K,INVD+1) = -AKKP1 / D
427: WORK(K-1,INVD+1) = -AKKP1 / D
428: K=K-2
429: END IF
430: END DO
431: *
1.3 bertrand 432: * inv(U**T) = (inv(U))**T
1.1 bertrand 433: *
1.3 bertrand 434: * inv(U**T)*inv(D)*inv(U)
1.1 bertrand 435: *
436: CUT=0
437: DO WHILE (CUT .LT. N)
438: NNB=NB
439: IF (CUT + NNB .GT. N) THEN
440: NNB=N-CUT
441: ELSE
442: COUNT = 0
443: * count negative elements,
444: DO I=CUT+1,CUT+NNB
445: IF (IPIV(I) .LT. 0) COUNT=COUNT+1
446: END DO
447: * need a even number for a clear cut
448: IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
449: END IF
450: * L21 Block
451: DO I=1,N-CUT-NNB
452: DO J=1,NNB
453: WORK(I,J)=A(CUT+NNB+I,CUT+J)
454: END DO
455: END DO
456: * L11 Block
457: DO I=1,NNB
458: WORK(U11+I,I)=ONE
459: DO J=I+1,NNB
460: WORK(U11+I,J)=ZERO
461: END DO
462: DO J=1,I-1
463: WORK(U11+I,J)=A(CUT+I,CUT+J)
464: END DO
465: END DO
466: *
467: * invD*L21
468: *
469: I=N-CUT-NNB
470: DO WHILE (I .GE. 1)
471: IF (IPIV(CUT+NNB+I) > 0) THEN
472: DO J=1,NNB
473: WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
474: END DO
475: I=I-1
476: ELSE
477: DO J=1,NNB
478: U01_I_J = WORK(I,J)
479: U01_IP1_J = WORK(I-1,J)
480: WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
481: $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
482: WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
483: $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
484: END DO
485: I=I-2
486: END IF
487: END DO
488: *
489: * invD1*L11
490: *
491: I=NNB
492: DO WHILE (I .GE. 1)
493: IF (IPIV(CUT+I) > 0) THEN
494: DO J=1,NNB
495: WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
496: END DO
497: I=I-1
498: ELSE
499: DO J=1,NNB
500: U11_I_J = WORK(U11+I,J)
501: U11_IP1_J = WORK(U11+I-1,J)
502: WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
503: $ WORK(CUT+I,INVD+1)*U11_IP1_J
504: WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
505: $ WORK(CUT+I-1,INVD)*U11_IP1_J
506: END DO
507: I=I-2
508: END IF
509: END DO
510: *
1.3 bertrand 511: * L11**T*invD1*L11->L11
1.1 bertrand 512: *
513: CALL DTRMM('L',UPLO,'T','U',NNB, NNB,
514: $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
515:
1.3 bertrand 516: *
517: DO I=1,NNB
518: DO J=1,I
519: A(CUT+I,CUT+J)=WORK(U11+I,J)
520: END DO
521: END DO
522: *
1.1 bertrand 523: IF ( (CUT+NNB) .LT. N ) THEN
524: *
1.3 bertrand 525: * L21**T*invD2*L21->A(CUT+I,CUT+J)
1.1 bertrand 526: *
527: CALL DGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
1.3 bertrand 528: $ ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
1.1 bertrand 529:
530: *
1.3 bertrand 531: * L11 = L11**T*invD1*L11 + U01**T*invD*U01
1.1 bertrand 532: *
533: DO I=1,NNB
534: DO J=1,I
535: A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
536: END DO
537: END DO
538: *
1.3 bertrand 539: * L01 = L22**T*invD2*L21
1.1 bertrand 540: *
541: CALL DTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB,
542: $ ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
543: *
544: * Update L21
545: *
546: DO I=1,N-CUT-NNB
547: DO J=1,NNB
548: A(CUT+NNB+I,CUT+J)=WORK(I,J)
549: END DO
550: END DO
551:
552: ELSE
553: *
1.3 bertrand 554: * L11 = L11**T*invD1*L11
1.1 bertrand 555: *
556: DO I=1,NNB
557: DO J=1,I
558: A(CUT+I,CUT+J)=WORK(U11+I,J)
559: END DO
560: END DO
561: END IF
562: *
563: * Next Block
564: *
565: CUT=CUT+NNB
566: END DO
567: *
1.3 bertrand 568: * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
1.1 bertrand 569: *
570: I=N
571: DO WHILE ( I .GE. 1 )
572: IF( IPIV(I) .GT. 0 ) THEN
573: IP=IPIV(I)
1.3 bertrand 574: IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP )
575: IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I )
1.1 bertrand 576: ELSE
577: IP=-IPIV(I)
1.3 bertrand 578: IF ( I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP )
579: IF ( I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP, I )
1.1 bertrand 580: I=I-1
581: ENDIF
582: I=I-1
583: END DO
584: END IF
585: *
586: RETURN
587: *
588: * End of DSYTRI2X
589: *
590: END
591:
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