LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
csytri2x.f
Go to the documentation of this file.
1*> \brief \b CSYTRI2X
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytri2x.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytri2x.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytri2x.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CSYTRI2X( 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*> CSYTRI2X computes the inverse of a real symmetric indefinite matrix
39*> A using the factorization A = U*D*U**T or A = L*D*L**T computed by
40*> CSYTRF.
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**T;
52*> = 'L': Lower triangular, form is A = L*D*L**T.
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 CSYTRF.
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 CSYTRF.
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 complexSYcomputational
117*
118* =====================================================================
119 SUBROUTINE csytri2x( 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 COMPLEX ONE, ZERO
138 parameter( one = ( 1.0e+0, 0.0e+0 ),
139 \$ zero = ( 0.0e+0, 0.0e+0 ) )
140* ..
141* .. Local Scalars ..
142 LOGICAL UPPER
143 INTEGER I, IINFO, IP, K, CUT, NNB
144 INTEGER COUNT
145 INTEGER J, U11, INVD
146
147 COMPLEX AK, AKKP1, AKP1, D, T
148 COMPLEX U01_I_J, U01_IP1_J
149 COMPLEX U11_I_J, U11_IP1_J
150* ..
151* .. External Functions ..
152 LOGICAL LSAME
153 EXTERNAL lsame
154* ..
155* .. External Subroutines ..
156 EXTERNAL csyconv, xerbla, ctrtri
157 EXTERNAL cgemm, ctrmm, csyswapr
158* ..
159* .. Intrinsic Functions ..
160 INTRINSIC max
161* ..
162* .. Executable Statements ..
163*
164* Test the input parameters.
165*
166 info = 0
167 upper = lsame( uplo, 'U' )
168 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
169 info = -1
170 ELSE IF( n.LT.0 ) THEN
171 info = -2
172 ELSE IF( lda.LT.max( 1, n ) ) THEN
173 info = -4
174 END IF
175*
176* Quick return if possible
177*
178*
179 IF( info.NE.0 ) THEN
180 CALL xerbla( 'CSYTRI2X', -info )
181 RETURN
182 END IF
183 IF( n.EQ.0 )
184 \$ RETURN
185*
186* Convert A
187* Workspace got Non-diag elements of D
188*
189 CALL csyconv( uplo, 'C', n, a, lda, ipiv, work, iinfo )
190*
191* Check that the diagonal matrix D is nonsingular.
192*
193 IF( upper ) THEN
194*
195* Upper triangular storage: examine D from bottom to top
196*
197 DO info = n, 1, -1
198 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
199 \$ RETURN
200 END DO
201 ELSE
202*
203* Lower triangular storage: examine D from top to bottom.
204*
205 DO info = 1, n
206 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
207 \$ RETURN
208 END DO
209 END IF
210 info = 0
211*
212* Splitting Workspace
213* U01 is a block (N,NB+1)
214* The first element of U01 is in WORK(1,1)
215* U11 is a block (NB+1,NB+1)
216* The first element of U11 is in WORK(N+1,1)
217 u11 = n
218* INVD is a block (N,2)
219* The first element of INVD is in WORK(1,INVD)
220 invd = nb+2
221
222 IF( upper ) THEN
223*
224* invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
225*
226 CALL ctrtri( uplo, 'U', n, a, lda, info )
227*
228* inv(D) and inv(D)*inv(U)
229*
230 k=1
231 DO WHILE ( k .LE. n )
232 IF( ipiv( k ).GT.0 ) THEN
233* 1 x 1 diagonal NNB
234 work(k,invd) = one / a( k, k )
235 work(k,invd+1) = 0
236 k=k+1
237 ELSE
238* 2 x 2 diagonal NNB
239 t = work(k+1,1)
240 ak = a( k, k ) / t
241 akp1 = a( k+1, k+1 ) / t
242 akkp1 = work(k+1,1) / t
243 d = t*( ak*akp1-one )
244 work(k,invd) = akp1 / d
245 work(k+1,invd+1) = ak / d
246 work(k,invd+1) = -akkp1 / d
247 work(k+1,invd) = -akkp1 / d
248 k=k+2
249 END IF
250 END DO
251*
252* inv(U**T) = (inv(U))**T
253*
254* inv(U**T)*inv(D)*inv(U)
255*
256 cut=n
257 DO WHILE (cut .GT. 0)
258 nnb=nb
259 IF (cut .LE. nnb) THEN
260 nnb=cut
261 ELSE
262 count = 0
263* count negative elements,
264 DO i=cut+1-nnb,cut
265 IF (ipiv(i) .LT. 0) count=count+1
266 END DO
267* need a even number for a clear cut
268 IF (mod(count,2) .EQ. 1) nnb=nnb+1
269 END IF
270
271 cut=cut-nnb
272*
273* U01 Block
274*
275 DO i=1,cut
276 DO j=1,nnb
277 work(i,j)=a(i,cut+j)
278 END DO
279 END DO
280*
281* U11 Block
282*
283 DO i=1,nnb
284 work(u11+i,i)=one
285 DO j=1,i-1
286 work(u11+i,j)=zero
287 END DO
288 DO j=i+1,nnb
289 work(u11+i,j)=a(cut+i,cut+j)
290 END DO
291 END DO
292*
293* invD*U01
294*
295 i=1
296 DO WHILE (i .LE. cut)
297 IF (ipiv(i) > 0) THEN
298 DO j=1,nnb
299 work(i,j)=work(i,invd)*work(i,j)
300 END DO
301 i=i+1
302 ELSE
303 DO j=1,nnb
304 u01_i_j = work(i,j)
305 u01_ip1_j = work(i+1,j)
306 work(i,j)=work(i,invd)*u01_i_j+
307 \$ work(i,invd+1)*u01_ip1_j
308 work(i+1,j)=work(i+1,invd)*u01_i_j+
309 \$ work(i+1,invd+1)*u01_ip1_j
310 END DO
311 i=i+2
312 END IF
313 END DO
314*
315* invD1*U11
316*
317 i=1
318 DO WHILE (i .LE. nnb)
319 IF (ipiv(cut+i) > 0) THEN
320 DO j=i,nnb
321 work(u11+i,j)=work(cut+i,invd)*work(u11+i,j)
322 END DO
323 i=i+1
324 ELSE
325 DO j=i,nnb
326 u11_i_j = work(u11+i,j)
327 u11_ip1_j = work(u11+i+1,j)
328 work(u11+i,j)=work(cut+i,invd)*work(u11+i,j) +
329 \$ work(cut+i,invd+1)*work(u11+i+1,j)
330 work(u11+i+1,j)=work(cut+i+1,invd)*u11_i_j+
331 \$ work(cut+i+1,invd+1)*u11_ip1_j
332 END DO
333 i=i+2
334 END IF
335 END DO
336*
337* U11**T*invD1*U11->U11
338*
339 CALL ctrmm('L','U','T','U',nnb, nnb,
340 \$ one,a(cut+1,cut+1),lda,work(u11+1,1),n+nb+1)
341*
342 DO i=1,nnb
343 DO j=i,nnb
344 a(cut+i,cut+j)=work(u11+i,j)
345 END DO
346 END DO
347*
348* U01**T*invD*U01->A(CUT+I,CUT+J)
349*
350 CALL cgemm('T','N',nnb,nnb,cut,one,a(1,cut+1),lda,
351 \$ work,n+nb+1, zero, work(u11+1,1), n+nb+1)
352*
353* U11 = U11**T*invD1*U11 + U01**T*invD*U01
354*
355 DO i=1,nnb
356 DO j=i,nnb
357 a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
358 END DO
359 END DO
360*
361* U01 = U00**T*invD0*U01
362*
363 CALL ctrmm('L',uplo,'T','U',cut, nnb,
364 \$ one,a,lda,work,n+nb+1)
365
366*
367* Update U01
368*
369 DO i=1,cut
370 DO j=1,nnb
371 a(i,cut+j)=work(i,j)
372 END DO
373 END DO
374*
375* Next Block
376*
377 END DO
378*
379* Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
380*
381 i=1
382 DO WHILE ( i .LE. n )
383 IF( ipiv(i) .GT. 0 ) THEN
384 ip=ipiv(i)
385 IF (i .LT. ip) CALL csyswapr( uplo, n, a, lda, i ,ip )
386 IF (i .GT. ip) CALL csyswapr( uplo, n, a, lda, ip ,i )
387 ELSE
388 ip=-ipiv(i)
389 i=i+1
390 IF ( (i-1) .LT. ip)
391 \$ CALL csyswapr( uplo, n, a, lda, i-1 ,ip )
392 IF ( (i-1) .GT. ip)
393 \$ CALL csyswapr( uplo, n, a, lda, ip ,i-1 )
394 ENDIF
395 i=i+1
396 END DO
397 ELSE
398*
399* LOWER...
400*
401* invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
402*
403 CALL ctrtri( uplo, 'U', n, a, lda, info )
404*
405* inv(D) and inv(D)*inv(U)
406*
407 k=n
408 DO WHILE ( k .GE. 1 )
409 IF( ipiv( k ).GT.0 ) THEN
410* 1 x 1 diagonal NNB
411 work(k,invd) = one / a( k, k )
412 work(k,invd+1) = 0
413 k=k-1
414 ELSE
415* 2 x 2 diagonal NNB
416 t = work(k-1,1)
417 ak = a( k-1, k-1 ) / t
418 akp1 = a( k, k ) / t
419 akkp1 = work(k-1,1) / t
420 d = t*( ak*akp1-one )
421 work(k-1,invd) = akp1 / d
422 work(k,invd) = ak / d
423 work(k,invd+1) = -akkp1 / d
424 work(k-1,invd+1) = -akkp1 / d
425 k=k-2
426 END IF
427 END DO
428*
429* inv(U**T) = (inv(U))**T
430*
431* inv(U**T)*inv(D)*inv(U)
432*
433 cut=0
434 DO WHILE (cut .LT. n)
435 nnb=nb
436 IF (cut + nnb .GE. n) THEN
437 nnb=n-cut
438 ELSE
439 count = 0
440* count negative elements,
441 DO i=cut+1,cut+nnb
442 IF (ipiv(i) .LT. 0) count=count+1
443 END DO
444* need a even number for a clear cut
445 IF (mod(count,2) .EQ. 1) nnb=nnb+1
446 END IF
447* L21 Block
448 DO i=1,n-cut-nnb
449 DO j=1,nnb
450 work(i,j)=a(cut+nnb+i,cut+j)
451 END DO
452 END DO
453* L11 Block
454 DO i=1,nnb
455 work(u11+i,i)=one
456 DO j=i+1,nnb
457 work(u11+i,j)=zero
458 END DO
459 DO j=1,i-1
460 work(u11+i,j)=a(cut+i,cut+j)
461 END DO
462 END DO
463*
464* invD*L21
465*
466 i=n-cut-nnb
467 DO WHILE (i .GE. 1)
468 IF (ipiv(cut+nnb+i) > 0) THEN
469 DO j=1,nnb
470 work(i,j)=work(cut+nnb+i,invd)*work(i,j)
471 END DO
472 i=i-1
473 ELSE
474 DO j=1,nnb
475 u01_i_j = work(i,j)
476 u01_ip1_j = work(i-1,j)
477 work(i,j)=work(cut+nnb+i,invd)*u01_i_j+
478 \$ work(cut+nnb+i,invd+1)*u01_ip1_j
479 work(i-1,j)=work(cut+nnb+i-1,invd+1)*u01_i_j+
480 \$ work(cut+nnb+i-1,invd)*u01_ip1_j
481 END DO
482 i=i-2
483 END IF
484 END DO
485*
486* invD1*L11
487*
488 i=nnb
489 DO WHILE (i .GE. 1)
490 IF (ipiv(cut+i) > 0) THEN
491 DO j=1,nnb
492 work(u11+i,j)=work(cut+i,invd)*work(u11+i,j)
493 END DO
494 i=i-1
495 ELSE
496 DO j=1,nnb
497 u11_i_j = work(u11+i,j)
498 u11_ip1_j = work(u11+i-1,j)
499 work(u11+i,j)=work(cut+i,invd)*work(u11+i,j) +
500 \$ work(cut+i,invd+1)*u11_ip1_j
501 work(u11+i-1,j)=work(cut+i-1,invd+1)*u11_i_j+
502 \$ work(cut+i-1,invd)*u11_ip1_j
503 END DO
504 i=i-2
505 END IF
506 END DO
507*
508* L11**T*invD1*L11->L11
509*
510 CALL ctrmm('L',uplo,'T','U',nnb, nnb,
511 \$ one,a(cut+1,cut+1),lda,work(u11+1,1),n+nb+1)
512*
513 DO i=1,nnb
514 DO j=1,i
515 a(cut+i,cut+j)=work(u11+i,j)
516 END DO
517 END DO
518*
519 IF ( (cut+nnb) .LT. n ) THEN
520*
521* L21**T*invD2*L21->A(CUT+I,CUT+J)
522*
523 CALL cgemm('T','N',nnb,nnb,n-nnb-cut,one,a(cut+nnb+1,cut+1)
524 \$ ,lda,work,n+nb+1, zero, work(u11+1,1), n+nb+1)
525
526*
527* L11 = L11**T*invD1*L11 + U01**T*invD*U01
528*
529 DO i=1,nnb
530 DO j=1,i
531 a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
532 END DO
533 END DO
534*
535* L01 = L22**T*invD2*L21
536*
537 CALL ctrmm('L',uplo,'T','U', n-nnb-cut, nnb,
538 \$ one,a(cut+nnb+1,cut+nnb+1),lda,work,n+nb+1)
539
540* Update L21
541 DO i=1,n-cut-nnb
542 DO j=1,nnb
543 a(cut+nnb+i,cut+j)=work(i,j)
544 END DO
545 END DO
546 ELSE
547*
548* L11 = L11**T*invD1*L11
549*
550 DO i=1,nnb
551 DO j=1,i
552 a(cut+i,cut+j)=work(u11+i,j)
553 END DO
554 END DO
555 END IF
556*
557* Next Block
558*
559 cut=cut+nnb
560 END DO
561*
562* Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
563*
564 i=n
565 DO WHILE ( i .GE. 1 )
566 IF( ipiv(i) .GT. 0 ) THEN
567 ip=ipiv(i)
568 IF (i .LT. ip) CALL csyswapr( uplo, n, a, lda, i ,ip )
569 IF (i .GT. ip) CALL csyswapr( uplo, n, a, lda, ip ,i )
570 ELSE
571 ip=-ipiv(i)
572 IF ( i .LT. ip) CALL csyswapr( uplo, n, a, lda, i ,ip )
573 IF ( i .GT. ip) CALL csyswapr( uplo, n, a, lda, ip ,i )
574 i=i-1
575 ENDIF
576 i=i-1
577 END DO
578 END IF
579*
580 RETURN
581*
582* End of CSYTRI2X
583*
584 END
585
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
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
subroutine csyswapr(UPLO, N, A, LDA, I1, I2)
CSYSWAPR
Definition: csyswapr.f:100
subroutine csyconv(UPLO, WAY, N, A, LDA, IPIV, E, INFO)
CSYCONV
Definition: csyconv.f:114
subroutine csytri2x(UPLO, N, A, LDA, IPIV, WORK, NB, INFO)
CSYTRI2X
Definition: csytri2x.f:120