LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine chetri_rook ( character  UPLO,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex, dimension( * )  WORK,
integer  INFO 
)

CHETRI_ROOK computes the inverse of HE matrix using the factorization obtained with the bounded Bunch-Kaufman ("rook") diagonal pivoting method.

Download CHETRI_ROOK + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
 CHETRI_ROOK computes the inverse of a complex Hermitian indefinite matrix
 A using the factorization A = U*D*U**H or A = L*D*L**H computed by
 CHETRF_ROOK.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the details of the factorization are stored
          as an upper or lower triangular matrix.
          = 'U':  Upper triangular, form is A = U*D*U**H;
          = 'L':  Lower triangular, form is A = L*D*L**H.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX array, dimension (LDA,N)
          On entry, the block diagonal matrix D and the multipliers
          used to obtain the factor U or L as computed by CHETRF_ROOK.

          On exit, if INFO = 0, the (Hermitian) inverse of the original
          matrix.  If UPLO = 'U', the upper triangular part of the
          inverse is formed and the part of A below the diagonal is not
          referenced; if UPLO = 'L' the lower triangular part of the
          inverse is formed and the part of A above the diagonal is
          not referenced.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D
          as determined by CHETRF_ROOK.
[out]WORK
          WORK is COMPLEX array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
               inverse could not be computed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2013
Contributors: 
  November 2013,  Igor Kozachenko,
                  Computer Science Division,
                  University of California, Berkeley

  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
                  School of Mathematics,
                  University of Manchester

Definition at line 130 of file chetri_rook.f.

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, $ K ), 1 ) )
244  END IF
245  kstep = 1
246  ELSE
247 *
248 * 2 x 2 diagonal block
249 *
250 * Invert the diagonal block.
251 *
252  t = abs( a( k, k+1 ) )
253  ak = REAL( A( K, K ) ) / t
254  akp1 = REAL( A( K+1, K+1 ) ) / t
255  akkp1 = a( k, k+1 ) / t
256  d = t*( ak*akp1-one )
257  a( k, k ) = akp1 / d
258  a( k+1, k+1 ) = ak / d
259  a( k, k+1 ) = -akkp1 / d
260 *
261 * Compute columns K and K+1 of the inverse.
262 *
263  IF( k.GT.1 ) THEN
264  CALL ccopy( k-1, a( 1, k ), 1, work, 1 )
265  CALL chemv( uplo, k-1, -cone, a, lda, work, 1, czero,
266  $ a( 1, k ), 1 )
267  a( k, k ) = a( k, k ) - REAL( CDOTC( K-1, WORK, 1, A( 1, $ K ), 1 ) )
268  a( k, k+1 ) = a( k, k+1 ) -
269  $ cdotc( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
270  CALL ccopy( k-1, a( 1, k+1 ), 1, work, 1 )
271  CALL chemv( uplo, k-1, -cone, a, lda, work, 1, czero,
272  $ a( 1, k+1 ), 1 )
273  a( k+1, k+1 ) = a( k+1, k+1 ) -
274  $ REAL( CDOTC( 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:k)
283 *
284  kp = ipiv( k )
285  IF( kp.NE.k ) THEN
286 *
287  IF( kp.GT.1 )
288  $ CALL cswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
289 *
290  DO 40 j = kp + 1, k - 1
291  temp = conjg( a( j, k ) )
292  a( j, k ) = conjg( a( kp, j ) )
293  a( kp, j ) = temp
294  40 CONTINUE
295 *
296  a( kp, k ) = conjg( 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 cswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
314 *
315  DO 50 j = kp + 1, k - 1
316  temp = conjg( a( j, k ) )
317  a( j, k ) = conjg( a( kp, j ) )
318  a( kp, j ) = temp
319  50 CONTINUE
320 *
321  a( kp, k ) = conjg( 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 cswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
340 *
341  DO 60 j = kp + 1, k - 1
342  temp = conjg( a( j, k ) )
343  a( j, k ) = conjg( a( kp, j ) )
344  a( kp, j ) = temp
345  60 CONTINUE
346 *
347  a( kp, k ) = conjg( 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 / REAL( A( K, K ) )
381 *
382 * Compute column K of the inverse.
383 *
384  IF( k.LT.n ) THEN
385  CALL ccopy( n-k, a( k+1, k ), 1, work, 1 )
386  CALL chemv( 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 ) - REAL( CDOTC( N-K, WORK, 1, $ A( K+1, K ), 1 ) )
389  END IF
390  kstep = 1
391  ELSE
392 *
393 * 2 x 2 diagonal block
394 *
395 * Invert the diagonal block.
396 *
397  t = abs( a( k, k-1 ) )
398  ak = REAL( A( K-1, K-1 ) ) / t
399  akp1 = REAL( A( K, K ) ) / t
400  akkp1 = a( k, k-1 ) / t
401  d = t*( ak*akp1-one )
402  a( k-1, k-1 ) = akp1 / d
403  a( k, k ) = ak / d
404  a( k, k-1 ) = -akkp1 / d
405 *
406 * Compute columns K-1 and K of the inverse.
407 *
408  IF( k.LT.n ) THEN
409  CALL ccopy( n-k, a( k+1, k ), 1, work, 1 )
410  CALL chemv( uplo, n-k, -cone, a( k+1, k+1 ), lda, work,
411  $ 1, czero, a( k+1, k ), 1 )
412  a( k, k ) = a( k, k ) - REAL( CDOTC( N-K, WORK, 1, $ A( K+1, K ), 1 ) )
413  a( k, k-1 ) = a( k, k-1 ) -
414  $ cdotc( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
415  $ 1 )
416  CALL ccopy( n-k, a( k+1, k-1 ), 1, work, 1 )
417  CALL chemv( uplo, n-k, -cone, a( k+1, k+1 ), lda, work,
418  $ 1, czero, a( k+1, k-1 ), 1 )
419  a( k-1, k-1 ) = a( k-1, k-1 ) -
420  $ REAL( CDOTC( N-K, WORK, 1, A( K+1, K-1 ), $ 1 ) )
421  END IF
422  kstep = 2
423  END IF
424 *
425  IF( kstep.EQ.1 ) THEN
426 *
427 * Interchange rows and columns K and IPIV(K) in the trailing
428 * submatrix A(k:n,k:n)
429 *
430  kp = ipiv( k )
431  IF( kp.NE.k ) THEN
432 *
433  IF( kp.LT.n )
434  $ CALL cswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
435 *
436  DO 90 j = k + 1, kp - 1
437  temp = conjg( a( j, k ) )
438  a( j, k ) = conjg( a( kp, j ) )
439  a( kp, j ) = temp
440  90 CONTINUE
441 *
442  a( kp, k ) = conjg( a( kp, k ) )
443 *
444  temp = a( k, k )
445  a( k, k ) = a( kp, kp )
446  a( kp, kp ) = temp
447  END IF
448  ELSE
449 *
450 * Interchange rows and columns K and K-1 with -IPIV(K) and
451 * -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
452 *
453 * (1) Interchange rows and columns K and -IPIV(K)
454 *
455  kp = -ipiv( k )
456  IF( kp.NE.k ) THEN
457 *
458  IF( kp.LT.n )
459  $ CALL cswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
460 *
461  DO 100 j = k + 1, kp - 1
462  temp = conjg( a( j, k ) )
463  a( j, k ) = conjg( a( kp, j ) )
464  a( kp, j ) = temp
465  100 CONTINUE
466 *
467  a( kp, k ) = conjg( a( kp, k ) )
468 *
469  temp = a( k, k )
470  a( k, k ) = a( kp, kp )
471  a( kp, kp ) = temp
472 *
473  temp = a( k, k-1 )
474  a( k, k-1 ) = a( kp, k-1 )
475  a( kp, k-1 ) = temp
476  END IF
477 *
478 * (2) Interchange rows and columns K-1 and -IPIV(K-1)
479 *
480  k = k - 1
481  kp = -ipiv( k )
482  IF( kp.NE.k ) THEN
483 *
484  IF( kp.LT.n )
485  $ CALL cswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
486 *
487  DO 110 j = k + 1, kp - 1
488  temp = conjg( a( j, k ) )
489  a( j, k ) = conjg( a( kp, j ) )
490  a( kp, j ) = temp
491  110 CONTINUE
492 *
493  a( kp, k ) = conjg( a( kp, k ) )
494 *
495  temp = a( k, k )
496  a( k, k ) = a( kp, kp )
497  a( kp, kp ) = temp
498  END IF
499  END IF
500 *
501  k = k - 1
502  GO TO 80
503  120 CONTINUE
504  END IF
505 *
506  RETURN
507 *
508 * End of CHETRI_ROOK
509 *
510 
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cdotc
complex function cdotc(N, CX, INCX, CY, INCY)
CDOTC
Definition: cdotc.f:54
chemv
subroutine chemv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CHEMV
Definition: chemv.f:156
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: