LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ 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.
Date
June 2017

Definition at line 122 of file csytri2x.f.

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