LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zhetri_rook.f
Go to the documentation of this file.
1*> \brief \b ZHETRI_ROOK computes the inverse of HE matrix using the factorization obtained with the bounded Bunch-Kaufman ("rook") diagonal pivoting method.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZHETRI_ROOK + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetri_rook.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetri_rook.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetri_rook.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZHETRI_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* COMPLEX*16 A( LDA, * ), WORK( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> ZHETRI_ROOK computes the inverse of a complex Hermitian indefinite matrix
39*> A using the factorization A = U*D*U**H or A = L*D*L**H computed by
40*> ZHETRF_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**H;
52*> = 'L': Lower triangular, form is A = L*D*L**H.
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 COMPLEX*16 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 ZHETRF_ROOK.
66*>
67*> On exit, if INFO = 0, the (Hermitian) 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 ZHETRF_ROOK.
86*> \endverbatim
87*>
88*> \param[out] WORK
89*> \verbatim
90*> WORK is COMPLEX*16 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 hetri_rook
111*
112*> \par Contributors:
113* ==================
114*>
115*> \verbatim
116*>
117*> November 2013, 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*> \endverbatim
125*
126* =====================================================================
127 SUBROUTINE zhetri_rook( UPLO, N, A, LDA, IPIV, WORK, INFO )
128*
129* -- LAPACK computational routine --
130* -- LAPACK is a software package provided by Univ. of Tennessee, --
131* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
132*
133* .. Scalar Arguments ..
134 CHARACTER UPLO
135 INTEGER INFO, LDA, N
136* ..
137* .. Array Arguments ..
138 INTEGER IPIV( * )
139 COMPLEX*16 A( LDA, * ), WORK( * )
140* ..
141*
142* =====================================================================
143*
144* .. Parameters ..
145 DOUBLE PRECISION ONE
146 COMPLEX*16 CONE, CZERO
147 parameter( one = 1.0d+0, cone = ( 1.0d+0, 0.0d+0 ),
148 $ czero = ( 0.0d+0, 0.0d+0 ) )
149* ..
150* .. Local Scalars ..
151 LOGICAL UPPER
152 INTEGER J, K, KP, KSTEP
153 DOUBLE PRECISION AK, AKP1, D, T
154 COMPLEX*16 AKKP1, TEMP
155* ..
156* .. External Functions ..
157 LOGICAL LSAME
158 COMPLEX*16 ZDOTC
159 EXTERNAL lsame, zdotc
160* ..
161* .. External Subroutines ..
162 EXTERNAL zcopy, zhemv, zswap, xerbla
163* ..
164* .. Intrinsic Functions ..
165 INTRINSIC abs, dconjg, max, dble
166* ..
167* .. Executable Statements ..
168*
169* Test the input parameters.
170*
171 info = 0
172 upper = lsame( uplo, 'U' )
173 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
174 info = -1
175 ELSE IF( n.LT.0 ) THEN
176 info = -2
177 ELSE IF( lda.LT.max( 1, n ) ) THEN
178 info = -4
179 END IF
180 IF( info.NE.0 ) THEN
181 CALL xerbla( 'ZHETRI_ROOK', -info )
182 RETURN
183 END IF
184*
185* Quick return if possible
186*
187 IF( n.EQ.0 )
188 $ RETURN
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 10 info = n, 1, -1
197 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.czero )
198 $ RETURN
199 10 CONTINUE
200 ELSE
201*
202* Lower triangular storage: examine D from top to bottom.
203*
204 DO 20 info = 1, n
205 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.czero )
206 $ RETURN
207 20 CONTINUE
208 END IF
209 info = 0
210*
211 IF( upper ) THEN
212*
213* Compute inv(A) from the factorization A = U*D*U**H.
214*
215* K is the main loop index, increasing from 1 to N in steps of
216* 1 or 2, depending on the size of the diagonal blocks.
217*
218 k = 1
219 30 CONTINUE
220*
221* If K > N, exit from loop.
222*
223 IF( k.GT.n )
224 $ GO TO 70
225*
226 IF( ipiv( k ).GT.0 ) THEN
227*
228* 1 x 1 diagonal block
229*
230* Invert the diagonal block.
231*
232 a( k, k ) = one / dble( a( k, k ) )
233*
234* Compute column K of the inverse.
235*
236 IF( k.GT.1 ) THEN
237 CALL zcopy( k-1, a( 1, k ), 1, work, 1 )
238 CALL zhemv( uplo, k-1, -cone, a, lda, work, 1, czero,
239 $ a( 1, k ), 1 )
240 a( k, k ) = a( k, k ) - dble( zdotc( k-1, work, 1, a( 1,
241 $ k ), 1 ) )
242 END IF
243 kstep = 1
244 ELSE
245*
246* 2 x 2 diagonal block
247*
248* Invert the diagonal block.
249*
250 t = abs( a( k, k+1 ) )
251 ak = dble( a( k, k ) ) / t
252 akp1 = dble( a( k+1, k+1 ) ) / t
253 akkp1 = a( k, k+1 ) / t
254 d = t*( ak*akp1-one )
255 a( k, k ) = akp1 / d
256 a( k+1, k+1 ) = ak / d
257 a( k, k+1 ) = -akkp1 / d
258*
259* Compute columns K and K+1 of the inverse.
260*
261 IF( k.GT.1 ) THEN
262 CALL zcopy( k-1, a( 1, k ), 1, work, 1 )
263 CALL zhemv( uplo, k-1, -cone, a, lda, work, 1, czero,
264 $ a( 1, k ), 1 )
265 a( k, k ) = a( k, k ) - dble( zdotc( k-1, work, 1, a( 1,
266 $ k ), 1 ) )
267 a( k, k+1 ) = a( k, k+1 ) -
268 $ zdotc( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
269 CALL zcopy( k-1, a( 1, k+1 ), 1, work, 1 )
270 CALL zhemv( uplo, k-1, -cone, a, lda, work, 1, czero,
271 $ a( 1, k+1 ), 1 )
272 a( k+1, k+1 ) = a( k+1, k+1 ) -
273 $ dble( zdotc( k-1, work, 1, a( 1, k+1 ),
274 $ 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:k)
283*
284 kp = ipiv( k )
285 IF( kp.NE.k ) THEN
286*
287 IF( kp.GT.1 )
288 $ CALL zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
289*
290 DO 40 j = kp + 1, k - 1
291 temp = dconjg( a( j, k ) )
292 a( j, k ) = dconjg( a( kp, j ) )
293 a( kp, j ) = temp
294 40 CONTINUE
295*
296 a( kp, k ) = dconjg( a( kp, k ) )
297*
298 temp = a( k, k )
299 a( k, k ) = a( kp, kp )
300 a( kp, kp ) = temp
301 END IF
302 ELSE
303*
304* Interchange rows and columns K and K+1 with -IPIV(K) and
305* -IPIV(K+1) in the leading submatrix A(k+1:n,k+1:n)
306*
307* (1) Interchange rows and columns K and -IPIV(K)
308*
309 kp = -ipiv( k )
310 IF( kp.NE.k ) THEN
311*
312 IF( kp.GT.1 )
313 $ CALL zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
314*
315 DO 50 j = kp + 1, k - 1
316 temp = dconjg( a( j, k ) )
317 a( j, k ) = dconjg( a( kp, j ) )
318 a( kp, j ) = temp
319 50 CONTINUE
320*
321 a( kp, k ) = dconjg( a( kp, k ) )
322*
323 temp = a( k, k )
324 a( k, k ) = a( kp, kp )
325 a( kp, kp ) = temp
326*
327 temp = a( k, k+1 )
328 a( k, k+1 ) = a( kp, k+1 )
329 a( kp, k+1 ) = temp
330 END IF
331*
332* (2) Interchange rows and columns K+1 and -IPIV(K+1)
333*
334 k = k + 1
335 kp = -ipiv( k )
336 IF( kp.NE.k ) THEN
337*
338 IF( kp.GT.1 )
339 $ CALL zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
340*
341 DO 60 j = kp + 1, k - 1
342 temp = dconjg( a( j, k ) )
343 a( j, k ) = dconjg( a( kp, j ) )
344 a( kp, j ) = temp
345 60 CONTINUE
346*
347 a( kp, k ) = dconjg( a( kp, k ) )
348*
349 temp = a( k, k )
350 a( k, k ) = a( kp, kp )
351 a( kp, kp ) = temp
352 END IF
353 END IF
354*
355 k = k + 1
356 GO TO 30
357 70 CONTINUE
358*
359 ELSE
360*
361* Compute inv(A) from the factorization A = L*D*L**H.
362*
363* K is the main loop index, decreasing from N to 1 in steps of
364* 1 or 2, depending on the size of the diagonal blocks.
365*
366 k = n
367 80 CONTINUE
368*
369* If K < 1, exit from loop.
370*
371 IF( k.LT.1 )
372 $ GO TO 120
373*
374 IF( ipiv( k ).GT.0 ) THEN
375*
376* 1 x 1 diagonal block
377*
378* Invert the diagonal block.
379*
380 a( k, k ) = one / dble( a( k, k ) )
381*
382* Compute column K of the inverse.
383*
384 IF( k.LT.n ) THEN
385 CALL zcopy( n-k, a( k+1, k ), 1, work, 1 )
386 CALL zhemv( uplo, n-k, -cone, a( k+1, k+1 ), lda, work,
387 $ 1, czero, a( k+1, k ), 1 )
388 a( k, k ) = a( k, k ) - dble( zdotc( n-k, work, 1,
389 $ a( k+1, k ), 1 ) )
390 END IF
391 kstep = 1
392 ELSE
393*
394* 2 x 2 diagonal block
395*
396* Invert the diagonal block.
397*
398 t = abs( a( k, k-1 ) )
399 ak = dble( a( k-1, k-1 ) ) / t
400 akp1 = dble( a( k, k ) ) / t
401 akkp1 = a( k, k-1 ) / t
402 d = t*( ak*akp1-one )
403 a( k-1, k-1 ) = akp1 / d
404 a( k, k ) = ak / d
405 a( k, k-1 ) = -akkp1 / d
406*
407* Compute columns K-1 and K of the inverse.
408*
409 IF( k.LT.n ) THEN
410 CALL zcopy( n-k, a( k+1, k ), 1, work, 1 )
411 CALL zhemv( uplo, n-k, -cone, a( k+1, k+1 ), lda, work,
412 $ 1, czero, a( k+1, k ), 1 )
413 a( k, k ) = a( k, k ) - dble( zdotc( n-k, work, 1,
414 $ a( k+1, k ), 1 ) )
415 a( k, k-1 ) = a( k, k-1 ) -
416 $ zdotc( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
417 $ 1 )
418 CALL zcopy( n-k, a( k+1, k-1 ), 1, work, 1 )
419 CALL zhemv( uplo, n-k, -cone, a( k+1, k+1 ), lda, work,
420 $ 1, czero, a( k+1, k-1 ), 1 )
421 a( k-1, k-1 ) = a( k-1, k-1 ) -
422 $ dble( zdotc( n-k, work, 1, a( k+1, k-1 ),
423 $ 1 ) )
424 END IF
425 kstep = 2
426 END IF
427*
428 IF( kstep.EQ.1 ) THEN
429*
430* Interchange rows and columns K and IPIV(K) in the trailing
431* submatrix A(k:n,k:n)
432*
433 kp = ipiv( k )
434 IF( kp.NE.k ) THEN
435*
436 IF( kp.LT.n )
437 $ CALL zswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
438*
439 DO 90 j = k + 1, kp - 1
440 temp = dconjg( a( j, k ) )
441 a( j, k ) = dconjg( a( kp, j ) )
442 a( kp, j ) = temp
443 90 CONTINUE
444*
445 a( kp, k ) = dconjg( a( kp, k ) )
446*
447 temp = a( k, k )
448 a( k, k ) = a( kp, kp )
449 a( kp, kp ) = temp
450 END IF
451 ELSE
452*
453* Interchange rows and columns K and K-1 with -IPIV(K) and
454* -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
455*
456* (1) Interchange rows and columns K and -IPIV(K)
457*
458 kp = -ipiv( k )
459 IF( kp.NE.k ) THEN
460*
461 IF( kp.LT.n )
462 $ CALL zswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
463*
464 DO 100 j = k + 1, kp - 1
465 temp = dconjg( a( j, k ) )
466 a( j, k ) = dconjg( a( kp, j ) )
467 a( kp, j ) = temp
468 100 CONTINUE
469*
470 a( kp, k ) = dconjg( a( kp, k ) )
471*
472 temp = a( k, k )
473 a( k, k ) = a( kp, kp )
474 a( kp, kp ) = temp
475*
476 temp = a( k, k-1 )
477 a( k, k-1 ) = a( kp, k-1 )
478 a( kp, k-1 ) = temp
479 END IF
480*
481* (2) Interchange rows and columns K-1 and -IPIV(K-1)
482*
483 k = k - 1
484 kp = -ipiv( k )
485 IF( kp.NE.k ) THEN
486*
487 IF( kp.LT.n )
488 $ CALL zswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
489*
490 DO 110 j = k + 1, kp - 1
491 temp = dconjg( a( j, k ) )
492 a( j, k ) = dconjg( a( kp, j ) )
493 a( kp, j ) = temp
494 110 CONTINUE
495*
496 a( kp, k ) = dconjg( a( kp, k ) )
497*
498 temp = a( k, k )
499 a( k, k ) = a( kp, kp )
500 a( kp, kp ) = temp
501 END IF
502 END IF
503*
504 k = k - 1
505 GO TO 80
506 120 CONTINUE
507 END IF
508*
509 RETURN
510*
511* End of ZHETRI_ROOK
512*
513 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
subroutine zhemv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
ZHEMV
Definition zhemv.f:154
subroutine zhetri_rook(uplo, n, a, lda, ipiv, work, info)
ZHETRI_ROOK computes the inverse of HE matrix using the factorization obtained with the bounded Bunch...
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81