LAPACK  3.10.1 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

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.```

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: