LAPACK  3.10.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.

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:102
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: