LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
chetri2x.f
Go to the documentation of this file.
1*> \brief \b CHETRI2X
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CHETRI2X + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetri2x.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chetri2x.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetri2x.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER UPLO
25* INTEGER INFO, LDA, N, NB
26* ..
27* .. Array Arguments ..
28* INTEGER IPIV( * )
29* COMPLEX A( LDA, * ), WORK( N+NB+1,* )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> CHETRI2X computes the inverse of a complex Hermitian indefinite matrix
39*> A using the factorization A = U*D*U**H or A = L*D*L**H computed by
40*> CHETRF.
41*> \endverbatim
42*
43* Arguments:
44* ==========
45*
46*> \param[in] UPLO
47*> \verbatim
48*> UPLO is CHARACTER*1
49*> Specifies whether the details of the factorization are stored
50*> as an upper or lower triangular matrix.
51*> = 'U': Upper triangular, form is A = U*D*U**H;
52*> = 'L': Lower triangular, form is A = L*D*L**H.
53*> \endverbatim
54*>
55*> \param[in] N
56*> \verbatim
57*> N is INTEGER
58*> The order of the matrix A. N >= 0.
59*> \endverbatim
60*>
61*> \param[in,out] A
62*> \verbatim
63*> A is COMPLEX array, dimension (LDA,N)
64*> On entry, the NNB diagonal matrix D and the multipliers
65*> used to obtain the factor U or L as computed by CHETRF.
66*>
67*> On exit, if INFO = 0, the (symmetric) inverse of the original
68*> matrix. If UPLO = 'U', the upper triangular part of the
69*> inverse is formed and the part of A below the diagonal is not
70*> referenced; if UPLO = 'L' the lower triangular part of the
71*> inverse is formed and the part of A above the diagonal is
72*> not referenced.
73*> \endverbatim
74*>
75*> \param[in] LDA
76*> \verbatim
77*> LDA is INTEGER
78*> The leading dimension of the array A. LDA >= max(1,N).
79*> \endverbatim
80*>
81*> \param[in] IPIV
82*> \verbatim
83*> IPIV is INTEGER array, dimension (N)
84*> Details of the interchanges and the NNB structure of D
85*> as determined by CHETRF.
86*> \endverbatim
87*>
88*> \param[out] WORK
89*> \verbatim
90*> WORK is COMPLEX array, dimension (N+NB+1,NB+3)
91*> \endverbatim
92*>
93*> \param[in] NB
94*> \verbatim
95*> NB is INTEGER
96*> Block size
97*> \endverbatim
98*>
99*> \param[out] INFO
100*> \verbatim
101*> INFO is INTEGER
102*> = 0: successful exit
103*> < 0: if INFO = -i, the i-th argument had an illegal value
104*> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
105*> inverse could not be computed.
106*> \endverbatim
107*
108* Authors:
109* ========
110*
111*> \author Univ. of Tennessee
112*> \author Univ. of California Berkeley
113*> \author Univ. of Colorado Denver
114*> \author NAG Ltd.
115*
116*> \ingroup hetri2x
117*
118* =====================================================================
119 SUBROUTINE chetri2x( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
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*
586 END
587
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
subroutine chetri2x(uplo, n, a, lda, ipiv, work, nb, info)
CHETRI2X
Definition chetri2x.f:120
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