LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
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

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.```

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)
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 csyswapr(uplo, n, a, lda, i1, i2)
CSYSWAPR
Definition csyswapr.f:100
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: