LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ zsytri2x()

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

ZSYTRI2X

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

Purpose:
 ZSYTRI2X computes the inverse of a complex symmetric indefinite matrix
 A using the factorization A = U*D*U**T or A = L*D*L**T computed by
 ZSYTRF.
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*16 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 ZSYTRF.

          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 ZSYTRF.
[out]WORK
          WORK is COMPLEX*16 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 zsytri2x.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*16 a( lda, * ), work( n+nb+1,* )
135 * ..
136 *
137 * =====================================================================
138 *
139 * .. Parameters ..
140  COMPLEX*16 one, zero
141  parameter( one = ( 1.0d+0, 0.0d+0 ),
142  $ zero = ( 0.0d+0, 0.0d+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*16 ak, akkp1, akp1, d, t
151  COMPLEX*16 u01_i_j, u01_ip1_j
152  COMPLEX*16 u11_i_j, u11_ip1_j
153 * ..
154 * .. External Functions ..
155  LOGICAL lsame
156  EXTERNAL lsame
157 * ..
158 * .. External Subroutines ..
159  EXTERNAL zsyconv, xerbla, ztrtri
160  EXTERNAL zgemm, ztrmm, zsyswapr
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( 'ZSYTRI2X', -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 zsyconv( 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 ztrtri( 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) = 1/ 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 ztrmm('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 zgemm('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 ztrmm('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 zsyswapr( uplo, n, a, lda, i ,ip )
389  IF (i .GT. ip) CALL zsyswapr( uplo, n, a, lda, ip ,i )
390  ELSE
391  ip=-ipiv(i)
392  i=i+1
393  IF ( (i-1) .LT. ip)
394  $ CALL zsyswapr( uplo, n, a, lda, i-1 ,ip )
395  IF ( (i-1) .GT. ip)
396  $ CALL zsyswapr( 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 ztrtri( 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) = 1/ 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 ztrmm('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 
523  IF ( (cut+nnb) .LT. n ) THEN
524 *
525 * L21**T*invD2*L21->A(CUT+I,CUT+J)
526 *
527  CALL zgemm('T','N',nnb,nnb,n-nnb-cut,one,a(cut+nnb+1,cut+1)
528  $ ,lda,work,n+nb+1, zero, work(u11+1,1), n+nb+1)
529 
530 *
531 * L11 = L11**T*invD1*L11 + U01**T*invD*U01
532 *
533  DO i=1,nnb
534  DO j=1,i
535  a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
536  END DO
537  END DO
538 *
539 * U01 = L22**T*invD2*L21
540 *
541  CALL ztrmm('L',uplo,'T','U', n-nnb-cut, nnb,
542  $ one,a(cut+nnb+1,cut+nnb+1),lda,work,n+nb+1)
543 
544 * Update L21
545  DO i=1,n-cut-nnb
546  DO j=1,nnb
547  a(cut+nnb+i,cut+j)=work(i,j)
548  END DO
549  END DO
550  ELSE
551 *
552 * L11 = L11**T*invD1*L11
553 *
554  DO i=1,nnb
555  DO j=1,i
556  a(cut+i,cut+j)=work(u11+i,j)
557  END DO
558  END DO
559  END IF
560 *
561 * Next Block
562 *
563  cut=cut+nnb
564  END DO
565 *
566 * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
567 *
568  i=n
569  DO WHILE ( i .GE. 1 )
570  IF( ipiv(i) .GT. 0 ) THEN
571  ip=ipiv(i)
572  IF (i .LT. ip) CALL zsyswapr( uplo, n, a, lda, i ,ip )
573  IF (i .GT. ip) CALL zsyswapr( uplo, n, a, lda, ip ,i )
574  ELSE
575  ip=-ipiv(i)
576  IF ( i .LT. ip) CALL zsyswapr( uplo, n, a, lda, i ,ip )
577  IF ( i .GT. ip) CALL zsyswapr( uplo, n, a, lda, ip ,i )
578  i=i-1
579  ENDIF
580  i=i-1
581  END DO
582  END IF
583 *
584  RETURN
585 *
586 * End of ZSYTRI2X
587 *
subroutine zsyconv(UPLO, WAY, N, A, LDA, IPIV, E, INFO)
ZSYCONV
Definition: zsyconv.f:116
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM
Definition: ztrmm.f:179
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine ztrtri(UPLO, DIAG, N, A, LDA, INFO)
ZTRTRI
Definition: ztrtri.f:111
subroutine zsyswapr(UPLO, N, A, LDA, I1, I2)
ZSYSWAPR
Definition: zsyswapr.f:104
Here is the call graph for this function:
Here is the caller graph for this function: