LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ chetri_rook()

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.
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 127 of file chetri_rook.f.

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 A( LDA, * ), WORK( * )
140 * ..
141 *
142 * =====================================================================
143 *
144 * .. Parameters ..
145  REAL ONE
146  COMPLEX CONE, CZERO
147  parameter( one = 1.0e+0, cone = ( 1.0e+0, 0.0e+0 ),
148  $ czero = ( 0.0e+0, 0.0e+0 ) )
149 * ..
150 * .. Local Scalars ..
151  LOGICAL UPPER
152  INTEGER J, K, KP, KSTEP
153  REAL AK, AKP1, D, T
154  COMPLEX AKKP1, TEMP
155 * ..
156 * .. External Functions ..
157  LOGICAL LSAME
158  COMPLEX CDOTC
159  EXTERNAL lsame, cdotc
160 * ..
161 * .. External Subroutines ..
162  EXTERNAL ccopy, chemv, cswap, xerbla
163 * ..
164 * .. Intrinsic Functions ..
165  INTRINSIC abs, conjg, max, real
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( 'CHETRI_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 / real( a( k, k ) )
233 *
234 * Compute column K of the inverse.
235 *
236  IF( k.GT.1 ) THEN
237  CALL ccopy( k-1, a( 1, k ), 1, work, 1 )
238  CALL chemv( uplo, k-1, -cone, a, lda, work, 1, czero,
239  $ a( 1, k ), 1 )
240  a( k, k ) = a( k, k ) - real( cdotc( 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 = real( a( k, k ) ) / t
252  akp1 = real( 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 ccopy( k-1, a( 1, k ), 1, work, 1 )
263  CALL chemv( uplo, k-1, -cone, a, lda, work, 1, czero,
264  $ a( 1, k ), 1 )
265  a( k, k ) = a( k, k ) - real( cdotc( k-1, work, 1, a( 1,
266  $ k ), 1 ) )
267  a( k, k+1 ) = a( k, k+1 ) -
268  $ cdotc( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
269  CALL ccopy( k-1, a( 1, k+1 ), 1, work, 1 )
270  CALL chemv( 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  $ real( cdotc( 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 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,
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 = real( a( k-1, k-1 ) ) / t
400  akp1 = real( 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 ccopy( n-k, a( k+1, k ), 1, work, 1 )
411  CALL chemv( 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 ) - real( cdotc( n-k, work, 1,
414  $ a( k+1, k ), 1 ) )
415  a( k, k-1 ) = a( k, k-1 ) -
416  $ cdotc( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
417  $ 1 )
418  CALL ccopy( n-k, a( k+1, k-1 ), 1, work, 1 )
419  CALL chemv( 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  $ real( cdotc( 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 cswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
438 *
439  DO 90 j = k + 1, kp - 1
440  temp = conjg( a( j, k ) )
441  a( j, k ) = conjg( a( kp, j ) )
442  a( kp, j ) = temp
443  90 CONTINUE
444 *
445  a( kp, k ) = conjg( 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 cswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
463 *
464  DO 100 j = k + 1, kp - 1
465  temp = conjg( a( j, k ) )
466  a( j, k ) = conjg( a( kp, j ) )
467  a( kp, j ) = temp
468  100 CONTINUE
469 *
470  a( kp, k ) = conjg( 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 cswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
489 *
490  DO 110 j = k + 1, kp - 1
491  temp = conjg( a( j, k ) )
492  a( j, k ) = conjg( a( kp, j ) )
493  a( kp, j ) = temp
494  110 CONTINUE
495 *
496  a( kp, k ) = conjg( 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 CHETRI_ROOK
512 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
complex function cdotc(N, CX, INCX, CY, INCY)
CDOTC
Definition: cdotc.f:83
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine chemv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CHEMV
Definition: chemv.f:154
Here is the call graph for this function:
Here is the caller graph for this function: