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