1: *> \brief \b DSYTRI_ROOK
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DSYTRI_ROOK + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytri_rook.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytri_rook.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytri_rook.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DSYTRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO )
22: *
23: * .. Scalar Arguments ..
24: * CHARACTER UPLO
25: * INTEGER INFO, LDA, N
26: * ..
27: * .. Array Arguments ..
28: * INTEGER IPIV( * )
29: * DOUBLE PRECISION A( LDA, * ), WORK( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DSYTRI_ROOK computes the inverse of a real symmetric
39: *> matrix A using the factorization A = U*D*U**T or A = L*D*L**T
40: *> computed by DSYTRF_ROOK.
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 block diagonal matrix D and the multipliers
65: *> used to obtain the factor U or L as computed by DSYTRF_ROOK.
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 block structure of D
85: *> as determined by DSYTRF_ROOK.
86: *> \endverbatim
87: *>
88: *> \param[out] WORK
89: *> \verbatim
90: *> WORK is DOUBLE PRECISION array, dimension (N)
91: *> \endverbatim
92: *>
93: *> \param[out] INFO
94: *> \verbatim
95: *> INFO is INTEGER
96: *> = 0: successful exit
97: *> < 0: if INFO = -i, the i-th argument had an illegal value
98: *> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
99: *> inverse could not be computed.
100: *> \endverbatim
101: *
102: * Authors:
103: * ========
104: *
105: *> \author Univ. of Tennessee
106: *> \author Univ. of California Berkeley
107: *> \author Univ. of Colorado Denver
108: *> \author NAG Ltd.
109: *
110: *> \date April 2012
111: *
112: *> \ingroup doubleSYcomputational
113: *
114: *> \par Contributors:
115: * ==================
116: *>
117: *> \verbatim
118: *>
119: *> April 2012, Igor Kozachenko,
120: *> Computer Science Division,
121: *> University of California, Berkeley
122: *>
123: *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
124: *> School of Mathematics,
125: *> University of Manchester
126: *>
127: *> \endverbatim
128: *
129: * =====================================================================
130: SUBROUTINE DSYTRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO )
131: *
132: * -- LAPACK computational routine (version 3.4.1) --
133: * -- LAPACK is a software package provided by Univ. of Tennessee, --
134: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
135: * April 2012
136: *
137: * .. Scalar Arguments ..
138: CHARACTER UPLO
139: INTEGER INFO, LDA, N
140: * ..
141: * .. Array Arguments ..
142: INTEGER IPIV( * )
143: DOUBLE PRECISION A( LDA, * ), WORK( * )
144: * ..
145: *
146: * =====================================================================
147: *
148: * .. Parameters ..
149: DOUBLE PRECISION ONE, ZERO
150: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
151: * ..
152: * .. Local Scalars ..
153: LOGICAL UPPER
154: INTEGER K, KP, KSTEP
155: DOUBLE PRECISION AK, AKKP1, AKP1, D, T, TEMP
156: * ..
157: * .. External Functions ..
158: LOGICAL LSAME
159: DOUBLE PRECISION DDOT
160: EXTERNAL LSAME, DDOT
161: * ..
162: * .. External Subroutines ..
163: EXTERNAL DCOPY, DSWAP, DSYMV, XERBLA
164: * ..
165: * .. Intrinsic Functions ..
166: INTRINSIC ABS, MAX
167: * ..
168: * .. Executable Statements ..
169: *
170: * Test the input parameters.
171: *
172: INFO = 0
173: UPPER = LSAME( UPLO, 'U' )
174: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
175: INFO = -1
176: ELSE IF( N.LT.0 ) THEN
177: INFO = -2
178: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
179: INFO = -4
180: END IF
181: IF( INFO.NE.0 ) THEN
182: CALL XERBLA( 'DSYTRI_ROOK', -INFO )
183: RETURN
184: END IF
185: *
186: * Quick return if possible
187: *
188: IF( N.EQ.0 )
189: $ RETURN
190: *
191: * Check that the diagonal matrix D is nonsingular.
192: *
193: IF( UPPER ) THEN
194: *
195: * Upper triangular storage: examine D from bottom to top
196: *
197: DO 10 INFO = N, 1, -1
198: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
199: $ RETURN
200: 10 CONTINUE
201: ELSE
202: *
203: * Lower triangular storage: examine D from top to bottom.
204: *
205: DO 20 INFO = 1, N
206: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
207: $ RETURN
208: 20 CONTINUE
209: END IF
210: INFO = 0
211: *
212: IF( UPPER ) THEN
213: *
214: * Compute inv(A) from the factorization A = U*D*U**T.
215: *
216: * K is the main loop index, increasing from 1 to N in steps of
217: * 1 or 2, depending on the size of the diagonal blocks.
218: *
219: K = 1
220: 30 CONTINUE
221: *
222: * If K > N, exit from loop.
223: *
224: IF( K.GT.N )
225: $ GO TO 40
226: *
227: IF( IPIV( K ).GT.0 ) THEN
228: *
229: * 1 x 1 diagonal block
230: *
231: * Invert the diagonal block.
232: *
233: A( K, K ) = ONE / A( K, K )
234: *
235: * Compute column K of the inverse.
236: *
237: IF( K.GT.1 ) THEN
238: CALL DCOPY( K-1, A( 1, K ), 1, WORK, 1 )
239: CALL DSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
240: $ A( 1, K ), 1 )
241: A( K, K ) = A( K, K ) - DDOT( K-1, WORK, 1, A( 1, K ),
242: $ 1 )
243: END IF
244: KSTEP = 1
245: ELSE
246: *
247: * 2 x 2 diagonal block
248: *
249: * Invert the diagonal block.
250: *
251: T = ABS( A( K, K+1 ) )
252: AK = A( K, K ) / T
253: AKP1 = A( K+1, K+1 ) / T
254: AKKP1 = A( K, K+1 ) / T
255: D = T*( AK*AKP1-ONE )
256: A( K, K ) = AKP1 / D
257: A( K+1, K+1 ) = AK / D
258: A( K, K+1 ) = -AKKP1 / D
259: *
260: * Compute columns K and K+1 of the inverse.
261: *
262: IF( K.GT.1 ) THEN
263: CALL DCOPY( K-1, A( 1, K ), 1, WORK, 1 )
264: CALL DSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
265: $ A( 1, K ), 1 )
266: A( K, K ) = A( K, K ) - DDOT( K-1, WORK, 1, A( 1, K ),
267: $ 1 )
268: A( K, K+1 ) = A( K, K+1 ) -
269: $ DDOT( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
270: CALL DCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
271: CALL DSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
272: $ A( 1, K+1 ), 1 )
273: A( K+1, K+1 ) = A( K+1, K+1 ) -
274: $ DDOT( K-1, WORK, 1, A( 1, K+1 ), 1 )
275: END IF
276: KSTEP = 2
277: END IF
278: *
279: IF( KSTEP.EQ.1 ) THEN
280: *
281: * Interchange rows and columns K and IPIV(K) in the leading
282: * submatrix A(1:k+1,1:k+1)
283: *
284: KP = IPIV( K )
285: IF( KP.NE.K ) THEN
286: IF( KP.GT.1 )
287: $ CALL DSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
288: CALL DSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
289: TEMP = A( K, K )
290: A( K, K ) = A( KP, KP )
291: A( KP, KP ) = TEMP
292: END IF
293: ELSE
294: *
295: * Interchange rows and columns K and K+1 with -IPIV(K) and
296: * -IPIV(K+1)in the leading submatrix A(1:k+1,1:k+1)
297: *
298: KP = -IPIV( K )
299: IF( KP.NE.K ) THEN
300: IF( KP.GT.1 )
301: $ CALL DSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
302: CALL DSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
303: *
304: TEMP = A( K, K )
305: A( K, K ) = A( KP, KP )
306: A( KP, KP ) = TEMP
307: TEMP = A( K, K+1 )
308: A( K, K+1 ) = A( KP, K+1 )
309: A( KP, K+1 ) = TEMP
310: END IF
311: *
312: K = K + 1
313: KP = -IPIV( K )
314: IF( KP.NE.K ) THEN
315: IF( KP.GT.1 )
316: $ CALL DSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
317: CALL DSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
318: TEMP = A( K, K )
319: A( K, K ) = A( KP, KP )
320: A( KP, KP ) = TEMP
321: END IF
322: END IF
323: *
324: K = K + 1
325: GO TO 30
326: 40 CONTINUE
327: *
328: ELSE
329: *
330: * Compute inv(A) from the factorization A = L*D*L**T.
331: *
332: * K is the main loop index, increasing from 1 to N in steps of
333: * 1 or 2, depending on the size of the diagonal blocks.
334: *
335: K = N
336: 50 CONTINUE
337: *
338: * If K < 1, exit from loop.
339: *
340: IF( K.LT.1 )
341: $ GO TO 60
342: *
343: IF( IPIV( K ).GT.0 ) THEN
344: *
345: * 1 x 1 diagonal block
346: *
347: * Invert the diagonal block.
348: *
349: A( K, K ) = ONE / A( K, K )
350: *
351: * Compute column K of the inverse.
352: *
353: IF( K.LT.N ) THEN
354: CALL DCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
355: CALL DSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
356: $ ZERO, A( K+1, K ), 1 )
357: A( K, K ) = A( K, K ) - DDOT( N-K, WORK, 1, A( K+1, K ),
358: $ 1 )
359: END IF
360: KSTEP = 1
361: ELSE
362: *
363: * 2 x 2 diagonal block
364: *
365: * Invert the diagonal block.
366: *
367: T = ABS( A( K, K-1 ) )
368: AK = A( K-1, K-1 ) / T
369: AKP1 = A( K, K ) / T
370: AKKP1 = A( K, K-1 ) / T
371: D = T*( AK*AKP1-ONE )
372: A( K-1, K-1 ) = AKP1 / D
373: A( K, K ) = AK / D
374: A( K, K-1 ) = -AKKP1 / D
375: *
376: * Compute columns K-1 and K of the inverse.
377: *
378: IF( K.LT.N ) THEN
379: CALL DCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
380: CALL DSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
381: $ ZERO, A( K+1, K ), 1 )
382: A( K, K ) = A( K, K ) - DDOT( N-K, WORK, 1, A( K+1, K ),
383: $ 1 )
384: A( K, K-1 ) = A( K, K-1 ) -
385: $ DDOT( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
386: $ 1 )
387: CALL DCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
388: CALL DSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
389: $ ZERO, A( K+1, K-1 ), 1 )
390: A( K-1, K-1 ) = A( K-1, K-1 ) -
391: $ DDOT( N-K, WORK, 1, A( K+1, K-1 ), 1 )
392: END IF
393: KSTEP = 2
394: END IF
395: *
396: IF( KSTEP.EQ.1 ) THEN
397: *
398: * Interchange rows and columns K and IPIV(K) in the trailing
399: * submatrix A(k-1:n,k-1:n)
400: *
401: KP = IPIV( K )
402: IF( KP.NE.K ) THEN
403: IF( KP.LT.N )
404: $ CALL DSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
405: CALL DSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
406: TEMP = A( K, K )
407: A( K, K ) = A( KP, KP )
408: A( KP, KP ) = TEMP
409: END IF
410: ELSE
411: *
412: * Interchange rows and columns K and K-1 with -IPIV(K) and
413: * -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
414: *
415: KP = -IPIV( K )
416: IF( KP.NE.K ) THEN
417: IF( KP.LT.N )
418: $ CALL DSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
419: CALL DSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
420: *
421: TEMP = A( K, K )
422: A( K, K ) = A( KP, KP )
423: A( KP, KP ) = TEMP
424: TEMP = A( K, K-1 )
425: A( K, K-1 ) = A( KP, K-1 )
426: A( KP, K-1 ) = TEMP
427: END IF
428: *
429: K = K - 1
430: KP = -IPIV( K )
431: IF( KP.NE.K ) THEN
432: IF( KP.LT.N )
433: $ CALL DSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
434: CALL DSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
435: TEMP = A( K, K )
436: A( K, K ) = A( KP, KP )
437: A( KP, KP ) = TEMP
438: END IF
439: END IF
440: *
441: K = K - 1
442: GO TO 50
443: 60 CONTINUE
444: END IF
445: *
446: RETURN
447: *
448: * End of DSYTRI_ROOK
449: *
450: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsytri_rook.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack