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: *> \ingroup doubleSYcomputational
111: *
112: *> \par Contributors:
113: * ==================
114: *>
115: *> \verbatim
116: *>
117: *> April 2012, Igor Kozachenko,
118: *> Computer Science Division,
119: *> University of California, Berkeley
120: *>
121: *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
122: *> School of Mathematics,
123: *> University of Manchester
124: *>
125: *> \endverbatim
126: *
127: * =====================================================================
128: SUBROUTINE DSYTRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO )
129: *
130: * -- LAPACK computational routine --
131: * -- LAPACK is a software package provided by Univ. of Tennessee, --
132: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
133: *
134: * .. Scalar Arguments ..
135: CHARACTER UPLO
136: INTEGER INFO, LDA, N
137: * ..
138: * .. Array Arguments ..
139: INTEGER IPIV( * )
140: DOUBLE PRECISION A( LDA, * ), WORK( * )
141: * ..
142: *
143: * =====================================================================
144: *
145: * .. Parameters ..
146: DOUBLE PRECISION ONE, ZERO
147: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
148: * ..
149: * .. Local Scalars ..
150: LOGICAL UPPER
151: INTEGER K, KP, KSTEP
152: DOUBLE PRECISION AK, AKKP1, AKP1, D, T, TEMP
153: * ..
154: * .. External Functions ..
155: LOGICAL LSAME
156: DOUBLE PRECISION DDOT
157: EXTERNAL LSAME, DDOT
158: * ..
159: * .. External Subroutines ..
160: EXTERNAL DCOPY, DSWAP, DSYMV, XERBLA
161: * ..
162: * .. Intrinsic Functions ..
163: INTRINSIC ABS, MAX
164: * ..
165: * .. Executable Statements ..
166: *
167: * Test the input parameters.
168: *
169: INFO = 0
170: UPPER = LSAME( UPLO, 'U' )
171: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
172: INFO = -1
173: ELSE IF( N.LT.0 ) THEN
174: INFO = -2
175: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
176: INFO = -4
177: END IF
178: IF( INFO.NE.0 ) THEN
179: CALL XERBLA( 'DSYTRI_ROOK', -INFO )
180: RETURN
181: END IF
182: *
183: * Quick return if possible
184: *
185: IF( N.EQ.0 )
186: $ RETURN
187: *
188: * Check that the diagonal matrix D is nonsingular.
189: *
190: IF( UPPER ) THEN
191: *
192: * Upper triangular storage: examine D from bottom to top
193: *
194: DO 10 INFO = N, 1, -1
195: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
196: $ RETURN
197: 10 CONTINUE
198: ELSE
199: *
200: * Lower triangular storage: examine D from top to bottom.
201: *
202: DO 20 INFO = 1, N
203: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
204: $ RETURN
205: 20 CONTINUE
206: END IF
207: INFO = 0
208: *
209: IF( UPPER ) THEN
210: *
211: * Compute inv(A) from the factorization A = U*D*U**T.
212: *
213: * K is the main loop index, increasing from 1 to N in steps of
214: * 1 or 2, depending on the size of the diagonal blocks.
215: *
216: K = 1
217: 30 CONTINUE
218: *
219: * If K > N, exit from loop.
220: *
221: IF( K.GT.N )
222: $ GO TO 40
223: *
224: IF( IPIV( K ).GT.0 ) THEN
225: *
226: * 1 x 1 diagonal block
227: *
228: * Invert the diagonal block.
229: *
230: A( K, K ) = ONE / A( K, K )
231: *
232: * Compute column K of the inverse.
233: *
234: IF( K.GT.1 ) THEN
235: CALL DCOPY( K-1, A( 1, K ), 1, WORK, 1 )
236: CALL DSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
237: $ A( 1, K ), 1 )
238: A( K, K ) = A( K, K ) - DDOT( K-1, WORK, 1, A( 1, K ),
239: $ 1 )
240: END IF
241: KSTEP = 1
242: ELSE
243: *
244: * 2 x 2 diagonal block
245: *
246: * Invert the diagonal block.
247: *
248: T = ABS( A( K, K+1 ) )
249: AK = A( K, K ) / T
250: AKP1 = A( K+1, K+1 ) / T
251: AKKP1 = A( K, K+1 ) / T
252: D = T*( AK*AKP1-ONE )
253: A( K, K ) = AKP1 / D
254: A( K+1, K+1 ) = AK / D
255: A( K, K+1 ) = -AKKP1 / D
256: *
257: * Compute columns K and K+1 of the inverse.
258: *
259: IF( K.GT.1 ) THEN
260: CALL DCOPY( K-1, A( 1, K ), 1, WORK, 1 )
261: CALL DSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
262: $ A( 1, K ), 1 )
263: A( K, K ) = A( K, K ) - DDOT( K-1, WORK, 1, A( 1, K ),
264: $ 1 )
265: A( K, K+1 ) = A( K, K+1 ) -
266: $ DDOT( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
267: CALL DCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
268: CALL DSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
269: $ A( 1, K+1 ), 1 )
270: A( K+1, K+1 ) = A( K+1, K+1 ) -
271: $ DDOT( K-1, WORK, 1, A( 1, K+1 ), 1 )
272: END IF
273: KSTEP = 2
274: END IF
275: *
276: IF( KSTEP.EQ.1 ) THEN
277: *
278: * Interchange rows and columns K and IPIV(K) in the leading
279: * submatrix A(1:k+1,1:k+1)
280: *
281: KP = IPIV( K )
282: IF( KP.NE.K ) THEN
283: IF( KP.GT.1 )
284: $ CALL DSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
285: CALL DSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
286: TEMP = A( K, K )
287: A( K, K ) = A( KP, KP )
288: A( KP, KP ) = TEMP
289: END IF
290: ELSE
291: *
292: * Interchange rows and columns K and K+1 with -IPIV(K) and
293: * -IPIV(K+1)in the leading submatrix A(1:k+1,1:k+1)
294: *
295: KP = -IPIV( K )
296: IF( KP.NE.K ) THEN
297: IF( KP.GT.1 )
298: $ CALL DSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
299: CALL DSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
300: *
301: TEMP = A( K, K )
302: A( K, K ) = A( KP, KP )
303: A( KP, KP ) = TEMP
304: TEMP = A( K, K+1 )
305: A( K, K+1 ) = A( KP, K+1 )
306: A( KP, K+1 ) = TEMP
307: END IF
308: *
309: K = K + 1
310: KP = -IPIV( K )
311: IF( KP.NE.K ) THEN
312: IF( KP.GT.1 )
313: $ CALL DSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
314: CALL DSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
315: TEMP = A( K, K )
316: A( K, K ) = A( KP, KP )
317: A( KP, KP ) = TEMP
318: END IF
319: END IF
320: *
321: K = K + 1
322: GO TO 30
323: 40 CONTINUE
324: *
325: ELSE
326: *
327: * Compute inv(A) from the factorization A = L*D*L**T.
328: *
329: * K is the main loop index, increasing from 1 to N in steps of
330: * 1 or 2, depending on the size of the diagonal blocks.
331: *
332: K = N
333: 50 CONTINUE
334: *
335: * If K < 1, exit from loop.
336: *
337: IF( K.LT.1 )
338: $ GO TO 60
339: *
340: IF( IPIV( K ).GT.0 ) THEN
341: *
342: * 1 x 1 diagonal block
343: *
344: * Invert the diagonal block.
345: *
346: A( K, K ) = ONE / A( K, K )
347: *
348: * Compute column K of the inverse.
349: *
350: IF( K.LT.N ) THEN
351: CALL DCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
352: CALL DSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
353: $ ZERO, A( K+1, K ), 1 )
354: A( K, K ) = A( K, K ) - DDOT( N-K, WORK, 1, A( K+1, K ),
355: $ 1 )
356: END IF
357: KSTEP = 1
358: ELSE
359: *
360: * 2 x 2 diagonal block
361: *
362: * Invert the diagonal block.
363: *
364: T = ABS( A( K, K-1 ) )
365: AK = A( K-1, K-1 ) / T
366: AKP1 = A( K, K ) / T
367: AKKP1 = A( K, K-1 ) / T
368: D = T*( AK*AKP1-ONE )
369: A( K-1, K-1 ) = AKP1 / D
370: A( K, K ) = AK / D
371: A( K, K-1 ) = -AKKP1 / D
372: *
373: * Compute columns K-1 and K of the inverse.
374: *
375: IF( K.LT.N ) THEN
376: CALL DCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
377: CALL DSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
378: $ ZERO, A( K+1, K ), 1 )
379: A( K, K ) = A( K, K ) - DDOT( N-K, WORK, 1, A( K+1, K ),
380: $ 1 )
381: A( K, K-1 ) = A( K, K-1 ) -
382: $ DDOT( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
383: $ 1 )
384: CALL DCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
385: CALL DSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
386: $ ZERO, A( K+1, K-1 ), 1 )
387: A( K-1, K-1 ) = A( K-1, K-1 ) -
388: $ DDOT( N-K, WORK, 1, A( K+1, K-1 ), 1 )
389: END IF
390: KSTEP = 2
391: END IF
392: *
393: IF( KSTEP.EQ.1 ) THEN
394: *
395: * Interchange rows and columns K and IPIV(K) in the trailing
396: * submatrix A(k-1:n,k-1:n)
397: *
398: KP = IPIV( K )
399: IF( KP.NE.K ) THEN
400: IF( KP.LT.N )
401: $ CALL DSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
402: CALL DSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
403: TEMP = A( K, K )
404: A( K, K ) = A( KP, KP )
405: A( KP, KP ) = TEMP
406: END IF
407: ELSE
408: *
409: * Interchange rows and columns K and K-1 with -IPIV(K) and
410: * -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
411: *
412: KP = -IPIV( K )
413: IF( KP.NE.K ) THEN
414: IF( KP.LT.N )
415: $ CALL DSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
416: CALL DSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
417: *
418: TEMP = A( K, K )
419: A( K, K ) = A( KP, KP )
420: A( KP, KP ) = TEMP
421: TEMP = A( K, K-1 )
422: A( K, K-1 ) = A( KP, K-1 )
423: A( KP, K-1 ) = TEMP
424: END IF
425: *
426: K = K - 1
427: KP = -IPIV( K )
428: IF( KP.NE.K ) THEN
429: IF( KP.LT.N )
430: $ CALL DSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
431: CALL DSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
432: TEMP = A( K, K )
433: A( K, K ) = A( KP, KP )
434: A( KP, KP ) = TEMP
435: END IF
436: END IF
437: *
438: K = K - 1
439: GO TO 50
440: 60 CONTINUE
441: END IF
442: *
443: RETURN
444: *
445: * End of DSYTRI_ROOK
446: *
447: 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