LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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: