LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zhetri2x()

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.

Definition at line 119 of file zhetri2x.f.

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