LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zhetri2x ( character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex*16, dimension( n+nb+1,* )  WORK,
integer  NB,
integer  INFO 
)

ZHETRI2X

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

Purpose:
 ZHETRI2X computes the inverse of a COMPLEX*16 Hermitian indefinite matrix
 A using the factorization A = U*D*U**H or A = L*D*L**H computed by
 ZHETRF.
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*16 array, dimension (LDA,N)
          On entry, the NNB diagonal matrix D and the multipliers
          used to obtain the factor U or L as computed by ZHETRF.

          On exit, if INFO = 0, the (symmetric) 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 NNB structure of D
          as determined by ZHETRF.
[out]WORK
          WORK is COMPLEX*16 array, dimension (N+NB+1,NB+3)
[in]NB
          NB is INTEGER
          Block size
[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 2015

Definition at line 122 of file zhetri2x.f.

122 *
123 * -- LAPACK computational routine (version 3.6.0) --
124 * -- LAPACK is a software package provided by Univ. of Tennessee, --
125 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
126 * November 2015
127 *
128 * .. Scalar Arguments ..
129  CHARACTER uplo
130  INTEGER info, lda, n, nb
131 * ..
132 * .. Array Arguments ..
133  INTEGER ipiv( * )
134  COMPLEX*16 a( lda, * ), work( n+nb+1,* )
135 * ..
136 *
137 * =====================================================================
138 *
139 * .. Parameters ..
140  REAL one
141  COMPLEX*16 cone, zero
142  parameter ( one = 1.0d+0,
143  $ cone = ( 1.0d+0, 0.0d+0 ),
144  $ zero = ( 0.0d+0, 0.0d+0 ) )
145 * ..
146 * .. Local Scalars ..
147  LOGICAL upper
148  INTEGER i, iinfo, ip, k, cut, nnb
149  INTEGER count
150  INTEGER j, u11, invd
151 
152  COMPLEX*16 ak, akkp1, akp1, d, t
153  COMPLEX*16 u01_i_j, u01_ip1_j
154  COMPLEX*16 u11_i_j, u11_ip1_j
155 * ..
156 * .. External Functions ..
157  LOGICAL lsame
158  EXTERNAL lsame
159 * ..
160 * .. External Subroutines ..
161  EXTERNAL zsyconv, xerbla, ztrtri
162  EXTERNAL zgemm, ztrmm, zheswapr
163 * ..
164 * .. Intrinsic Functions ..
165  INTRINSIC max
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 *
181 * Quick return if possible
182 *
183 *
184  IF( info.NE.0 ) THEN
185  CALL xerbla( 'ZHETRI2X', -info )
186  RETURN
187  END IF
188  IF( n.EQ.0 )
189  $ RETURN
190 *
191 * Convert A
192 * Workspace got Non-diag elements of D
193 *
194  CALL zsyconv( uplo, 'C', n, a, lda, ipiv, work, iinfo )
195 *
196 * Check that the diagonal matrix D is nonsingular.
197 *
198  IF( upper ) THEN
199 *
200 * Upper triangular storage: examine D from bottom to top
201 *
202  DO info = n, 1, -1
203  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
204  $ RETURN
205  END DO
206  ELSE
207 *
208 * Lower triangular storage: examine D from top to bottom.
209 *
210  DO info = 1, n
211  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
212  $ RETURN
213  END DO
214  END IF
215  info = 0
216 *
217 * Splitting Workspace
218 * U01 is a block (N,NB+1)
219 * The first element of U01 is in WORK(1,1)
220 * U11 is a block (NB+1,NB+1)
221 * The first element of U11 is in WORK(N+1,1)
222  u11 = n
223 * INVD is a block (N,2)
224 * The first element of INVD is in WORK(1,INVD)
225  invd = nb+2
226 
227  IF( upper ) THEN
228 *
229 * invA = P * inv(U**H)*inv(D)*inv(U)*P**H.
230 *
231  CALL ztrtri( uplo, 'U', n, a, lda, info )
232 *
233 * inv(D) and inv(D)*inv(U)
234 *
235  k=1
236  DO WHILE ( k .LE. n )
237  IF( ipiv( k ).GT.0 ) THEN
238 * 1 x 1 diagonal NNB
239  work(k,invd) = one / REAL ( A( K, K ) )
240  work(k,invd+1) = 0
241  k=k+1
242  ELSE
243 * 2 x 2 diagonal NNB
244  t = abs( work(k+1,1) )
245  ak = REAL ( A( K, K ) ) / t
246  akp1 = REAL ( A( K+1, K+1 ) ) / t
247  akkp1 = work(k+1,1) / t
248  d = t*( ak*akp1-one )
249  work(k,invd) = akp1 / d
250  work(k+1,invd+1) = ak / d
251  work(k,invd+1) = -akkp1 / d
252  work(k+1,invd) = dconjg(work(k,invd+1) )
253  k=k+2
254  END IF
255  END DO
256 *
257 * inv(U**H) = (inv(U))**H
258 *
259 * inv(U**H)*inv(D)*inv(U)
260 *
261  cut=n
262  DO WHILE (cut .GT. 0)
263  nnb=nb
264  IF (cut .LE. nnb) THEN
265  nnb=cut
266  ELSE
267  count = 0
268 * count negative elements,
269  DO i=cut+1-nnb,cut
270  IF (ipiv(i) .LT. 0) count=count+1
271  END DO
272 * need a even number for a clear cut
273  IF (mod(count,2) .EQ. 1) nnb=nnb+1
274  END IF
275 
276  cut=cut-nnb
277 *
278 * U01 Block
279 *
280  DO i=1,cut
281  DO j=1,nnb
282  work(i,j)=a(i,cut+j)
283  END DO
284  END DO
285 *
286 * U11 Block
287 *
288  DO i=1,nnb
289  work(u11+i,i)=cone
290  DO j=1,i-1
291  work(u11+i,j)=zero
292  END DO
293  DO j=i+1,nnb
294  work(u11+i,j)=a(cut+i,cut+j)
295  END DO
296  END DO
297 *
298 * invD*U01
299 *
300  i=1
301  DO WHILE (i .LE. cut)
302  IF (ipiv(i) > 0) THEN
303  DO j=1,nnb
304  work(i,j)=work(i,invd)*work(i,j)
305  END DO
306  i=i+1
307  ELSE
308  DO j=1,nnb
309  u01_i_j = work(i,j)
310  u01_ip1_j = work(i+1,j)
311  work(i,j)=work(i,invd)*u01_i_j+
312  $ work(i,invd+1)*u01_ip1_j
313  work(i+1,j)=work(i+1,invd)*u01_i_j+
314  $ work(i+1,invd+1)*u01_ip1_j
315  END DO
316  i=i+2
317  END IF
318  END DO
319 *
320 * invD1*U11
321 *
322  i=1
323  DO WHILE (i .LE. nnb)
324  IF (ipiv(cut+i) > 0) THEN
325  DO j=i,nnb
326  work(u11+i,j)=work(cut+i,invd)*work(u11+i,j)
327  END DO
328  i=i+1
329  ELSE
330  DO j=i,nnb
331  u11_i_j = work(u11+i,j)
332  u11_ip1_j = work(u11+i+1,j)
333  work(u11+i,j)=work(cut+i,invd)*work(u11+i,j) +
334  $ work(cut+i,invd+1)*work(u11+i+1,j)
335  work(u11+i+1,j)=work(cut+i+1,invd)*u11_i_j+
336  $ work(cut+i+1,invd+1)*u11_ip1_j
337  END DO
338  i=i+2
339  END IF
340  END DO
341 *
342 * U11**H*invD1*U11->U11
343 *
344  CALL ztrmm('L','U','C','U',nnb, nnb,
345  $ cone,a(cut+1,cut+1),lda,work(u11+1,1),n+nb+1)
346 *
347  DO i=1,nnb
348  DO j=i,nnb
349  a(cut+i,cut+j)=work(u11+i,j)
350  END DO
351  END DO
352 *
353 * U01**H*invD*U01->A(CUT+I,CUT+J)
354 *
355  CALL zgemm('C','N',nnb,nnb,cut,cone,a(1,cut+1),lda,
356  $ work,n+nb+1, zero, work(u11+1,1), n+nb+1)
357 *
358 * U11 = U11**H*invD1*U11 + U01**H*invD*U01
359 *
360  DO i=1,nnb
361  DO j=i,nnb
362  a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
363  END DO
364  END DO
365 *
366 * U01 = U00**H*invD0*U01
367 *
368  CALL ztrmm('L',uplo,'C','U',cut, nnb,
369  $ cone,a,lda,work,n+nb+1)
370 
371 *
372 * Update U01
373 *
374  DO i=1,cut
375  DO j=1,nnb
376  a(i,cut+j)=work(i,j)
377  END DO
378  END DO
379 *
380 * Next Block
381 *
382  END DO
383 *
384 * Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H
385 *
386  i=1
387  DO WHILE ( i .LE. n )
388  IF( ipiv(i) .GT. 0 ) THEN
389  ip=ipiv(i)
390  IF (i .LT. ip) CALL zheswapr( uplo, n, a, lda, i ,ip )
391  IF (i .GT. ip) CALL zheswapr( uplo, n, a, lda, ip ,i )
392  ELSE
393  ip=-ipiv(i)
394  i=i+1
395  IF ( (i-1) .LT. ip)
396  $ CALL zheswapr( uplo, n, a, lda, i-1 ,ip )
397  IF ( (i-1) .GT. ip)
398  $ CALL zheswapr( uplo, n, a, lda, ip ,i-1 )
399  ENDIF
400  i=i+1
401  END DO
402  ELSE
403 *
404 * LOWER...
405 *
406 * invA = P * inv(U**H)*inv(D)*inv(U)*P**H.
407 *
408  CALL ztrtri( uplo, 'U', n, a, lda, info )
409 *
410 * inv(D) and inv(D)*inv(U)
411 *
412  k=n
413  DO WHILE ( k .GE. 1 )
414  IF( ipiv( k ).GT.0 ) THEN
415 * 1 x 1 diagonal NNB
416  work(k,invd) = one / REAL ( A( K, K ) )
417  work(k,invd+1) = 0
418  k=k-1
419  ELSE
420 * 2 x 2 diagonal NNB
421  t = abs( work(k-1,1) )
422  ak = REAL ( A( K-1, K-1 ) ) / t
423  akp1 = REAL ( A( K, K ) ) / t
424  akkp1 = work(k-1,1) / t
425  d = t*( ak*akp1-one )
426  work(k-1,invd) = akp1 / d
427  work(k,invd) = ak / d
428  work(k,invd+1) = -akkp1 / d
429  work(k-1,invd+1) = dconjg(work(k,invd+1) )
430  k=k-2
431  END IF
432  END DO
433 *
434 * inv(U**H) = (inv(U))**H
435 *
436 * inv(U**H)*inv(D)*inv(U)
437 *
438  cut=0
439  DO WHILE (cut .LT. n)
440  nnb=nb
441  IF (cut + nnb .GE. n) THEN
442  nnb=n-cut
443  ELSE
444  count = 0
445 * count negative elements,
446  DO i=cut+1,cut+nnb
447  IF (ipiv(i) .LT. 0) count=count+1
448  END DO
449 * need a even number for a clear cut
450  IF (mod(count,2) .EQ. 1) nnb=nnb+1
451  END IF
452 * L21 Block
453  DO i=1,n-cut-nnb
454  DO j=1,nnb
455  work(i,j)=a(cut+nnb+i,cut+j)
456  END DO
457  END DO
458 * L11 Block
459  DO i=1,nnb
460  work(u11+i,i)=cone
461  DO j=i+1,nnb
462  work(u11+i,j)=zero
463  END DO
464  DO j=1,i-1
465  work(u11+i,j)=a(cut+i,cut+j)
466  END DO
467  END DO
468 *
469 * invD*L21
470 *
471  i=n-cut-nnb
472  DO WHILE (i .GE. 1)
473  IF (ipiv(cut+nnb+i) > 0) THEN
474  DO j=1,nnb
475  work(i,j)=work(cut+nnb+i,invd)*work(i,j)
476  END DO
477  i=i-1
478  ELSE
479  DO j=1,nnb
480  u01_i_j = work(i,j)
481  u01_ip1_j = work(i-1,j)
482  work(i,j)=work(cut+nnb+i,invd)*u01_i_j+
483  $ work(cut+nnb+i,invd+1)*u01_ip1_j
484  work(i-1,j)=work(cut+nnb+i-1,invd+1)*u01_i_j+
485  $ work(cut+nnb+i-1,invd)*u01_ip1_j
486  END DO
487  i=i-2
488  END IF
489  END DO
490 *
491 * invD1*L11
492 *
493  i=nnb
494  DO WHILE (i .GE. 1)
495  IF (ipiv(cut+i) > 0) THEN
496  DO j=1,nnb
497  work(u11+i,j)=work(cut+i,invd)*work(u11+i,j)
498  END DO
499  i=i-1
500  ELSE
501  DO j=1,nnb
502  u11_i_j = work(u11+i,j)
503  u11_ip1_j = work(u11+i-1,j)
504  work(u11+i,j)=work(cut+i,invd)*work(u11+i,j) +
505  $ work(cut+i,invd+1)*u11_ip1_j
506  work(u11+i-1,j)=work(cut+i-1,invd+1)*u11_i_j+
507  $ work(cut+i-1,invd)*u11_ip1_j
508  END DO
509  i=i-2
510  END IF
511  END DO
512 *
513 * L11**H*invD1*L11->L11
514 *
515  CALL ztrmm('L',uplo,'C','U',nnb, nnb,
516  $ cone,a(cut+1,cut+1),lda,work(u11+1,1),n+nb+1)
517 *
518  DO i=1,nnb
519  DO j=1,i
520  a(cut+i,cut+j)=work(u11+i,j)
521  END DO
522  END DO
523 *
524  IF ( (cut+nnb) .LT. n ) THEN
525 *
526 * L21**H*invD2*L21->A(CUT+I,CUT+J)
527 *
528  CALL zgemm('C','N',nnb,nnb,n-nnb-cut,cone,a(cut+nnb+1,cut+1)
529  $ ,lda,work,n+nb+1, zero, work(u11+1,1), n+nb+1)
530 
531 *
532 * L11 = L11**H*invD1*L11 + U01**H*invD*U01
533 *
534  DO i=1,nnb
535  DO j=1,i
536  a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
537  END DO
538  END DO
539 *
540 * L01 = L22**H*invD2*L21
541 *
542  CALL ztrmm('L',uplo,'C','U', n-nnb-cut, nnb,
543  $ cone,a(cut+nnb+1,cut+nnb+1),lda,work,n+nb+1)
544 
545 * Update L21
546  DO i=1,n-cut-nnb
547  DO j=1,nnb
548  a(cut+nnb+i,cut+j)=work(i,j)
549  END DO
550  END DO
551  ELSE
552 *
553 * L11 = L11**H*invD1*L11
554 *
555  DO i=1,nnb
556  DO j=1,i
557  a(cut+i,cut+j)=work(u11+i,j)
558  END DO
559  END DO
560  END IF
561 *
562 * Next Block
563 *
564  cut=cut+nnb
565  END DO
566 *
567 * Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H
568 *
569  i=n
570  DO WHILE ( i .GE. 1 )
571  IF( ipiv(i) .GT. 0 ) THEN
572  ip=ipiv(i)
573  IF (i .LT. ip) CALL zheswapr( uplo, n, a, lda, i ,ip )
574  IF (i .GT. ip) CALL zheswapr( uplo, n, a, lda, ip ,i )
575  ELSE
576  ip=-ipiv(i)
577  IF ( i .LT. ip) CALL zheswapr( uplo, n, a, lda, i ,ip )
578  IF ( i .GT. ip) CALL zheswapr( uplo, n, a, lda, ip ,i )
579  i=i-1
580  ENDIF
581  i=i-1
582  END DO
583  END IF
584 *
585  RETURN
586 *
587 * End of ZHETRI2X
588 *
subroutine zsyconv(UPLO, WAY, N, A, LDA, IPIV, E, INFO)
ZSYCONV
Definition: zsyconv.f:116
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM
Definition: ztrmm.f:179
subroutine ztrtri(UPLO, DIAG, N, A, LDA, INFO)
ZTRTRI
Definition: ztrtri.f:111
subroutine zheswapr(UPLO, N, A, LDA, I1, I2)
ZHESWAPR applies an elementary permutation on the rows and columns of a Hermitian matrix...
Definition: zheswapr.f:104
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: