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

◆ csytri2x()

subroutine csytri2x ( character  UPLO,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex, dimension( n+nb+1,* )  WORK,
integer  NB,
integer  INFO 
)

CSYTRI2X

Download CSYTRI2X + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 CSYTRI2X computes the inverse of a real symmetric indefinite matrix
 A using the factorization A = U*D*U**T or A = L*D*L**T computed by
 CSYTRF.
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**T;
          = 'L':  Lower triangular, form is A = L*D*L**T.
[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 CSYTRF.

          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 CSYTRF.
[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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 119 of file csytri2x.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 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*
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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
Here is the call graph for this function:
Here is the caller graph for this function: