LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ chetri2x()

 subroutine chetri2x ( character uplo, integer n, complex, dimension( lda, * ) a, integer lda, integer, dimension( * ) ipiv, complex, dimension( n+nb+1,* ) work, integer nb, integer info )

CHETRI2X

Purpose:
``` CHETRI2X 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.```
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 NNB diagonal matrix D and the multipliers used to obtain the factor U or L as computed by CHETRF. 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 CHETRF.``` [out] WORK ` WORK is COMPLEX 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 chetri2x.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 A( LDA, * ), WORK( N+NB+1,* )
132* ..
133*
134* =====================================================================
135*
136* .. Parameters ..
137 REAL ONE
138 COMPLEX CONE, ZERO
139 parameter( one = 1.0e+0,
140 \$ cone = ( 1.0e+0, 0.0e+0 ),
141 \$ zero = ( 0.0e+0, 0.0e+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 AK, AKKP1, AKP1, D, T
150 COMPLEX U01_I_J, U01_IP1_J
151 COMPLEX U11_I_J, U11_IP1_J
152* ..
153* .. External Functions ..
154 LOGICAL LSAME
155 EXTERNAL lsame
156* ..
157* .. External Subroutines ..
158 EXTERNAL csyconv, xerbla, ctrtri
159 EXTERNAL cgemm, ctrmm, cheswapr
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( 'CHETRI2X', -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 csyconv( 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 ctrtri( 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 = real( a( k, k ) ) / t
243 akp1 = real( 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) = conjg(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 ctrmm('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 cgemm('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 ctrmm('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 cheswapr( uplo, n, a, lda, i ,ip )
388 IF (i .GT. ip) CALL cheswapr( uplo, n, a, lda, ip ,i )
389 ELSE
390 ip=-ipiv(i)
391 i=i+1
392 IF ( (i-1) .LT. ip)
393 \$ CALL cheswapr( uplo, n, a, lda, i-1 ,ip )
394 IF ( (i-1) .GT. ip)
395 \$ CALL cheswapr( 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 ctrtri( 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 = real( a( k-1, k-1 ) ) / t
420 akp1 = real( 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) = conjg(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 ctrmm('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 cgemm('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 ctrmm('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 cheswapr( uplo, n, a, lda, i ,ip )
571 IF (i .GT. ip) CALL cheswapr( uplo, n, a, lda, ip ,i )
572 ELSE
573 ip=-ipiv(i)
574 IF ( i .LT. ip) CALL cheswapr( uplo, n, a, lda, i ,ip )
575 IF ( i .GT. ip) CALL cheswapr( 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 CHETRI2X
585*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
subroutine cheswapr(uplo, n, a, lda, i1, i2)
CHESWAPR applies an elementary permutation on the rows and columns of a Hermitian matrix.
Definition cheswapr.f:102
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine csyconv(uplo, way, n, a, lda, ipiv, e, info)
CSYCONV
Definition csyconv.f:114
subroutine ctrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRMM
Definition ctrmm.f:177
subroutine ctrtri(uplo, diag, n, a, lda, info)
CTRTRI
Definition ctrtri.f:109
Here is the call graph for this function:
Here is the caller graph for this function: