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

Definition at line 119 of file zsytri2x.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*16 A( LDA, * ), WORK( N+NB+1,* )
132 * ..
133 *
134 * =====================================================================
135 *
136 * .. Parameters ..
137  COMPLEX*16 ONE, ZERO
138  parameter( one = ( 1.0d+0, 0.0d+0 ),
139  $ zero = ( 0.0d+0, 0.0d+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*16 AK, AKKP1, AKP1, D, T
148  COMPLEX*16 U01_I_J, U01_IP1_J
149  COMPLEX*16 U11_I_J, U11_IP1_J
150 * ..
151 * .. External Functions ..
152  LOGICAL LSAME
153  EXTERNAL lsame
154 * ..
155 * .. External Subroutines ..
156  EXTERNAL zsyconv, xerbla, ztrtri
157  EXTERNAL zgemm, ztrmm, zsyswapr
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( 'ZSYTRI2X', -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 zsyconv( 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 ztrtri( 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 ztrmm('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 zgemm('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 ztrmm('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 zsyswapr( uplo, n, a, lda, i ,ip )
386  IF (i .GT. ip) CALL zsyswapr( uplo, n, a, lda, ip ,i )
387  ELSE
388  ip=-ipiv(i)
389  i=i+1
390  IF ( (i-1) .LT. ip)
391  $ CALL zsyswapr( uplo, n, a, lda, i-1 ,ip )
392  IF ( (i-1) .GT. ip)
393  $ CALL zsyswapr( 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 ztrtri( 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 ztrmm('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 
520  IF ( (cut+nnb) .LT. n ) THEN
521 *
522 * L21**T*invD2*L21->A(CUT+I,CUT+J)
523 *
524  CALL zgemm('T','N',nnb,nnb,n-nnb-cut,one,a(cut+nnb+1,cut+1)
525  $ ,lda,work,n+nb+1, zero, work(u11+1,1), n+nb+1)
526 
527 *
528 * L11 = L11**T*invD1*L11 + U01**T*invD*U01
529 *
530  DO i=1,nnb
531  DO j=1,i
532  a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
533  END DO
534  END DO
535 *
536 * U01 = L22**T*invD2*L21
537 *
538  CALL ztrmm('L',uplo,'T','U', n-nnb-cut, nnb,
539  $ one,a(cut+nnb+1,cut+nnb+1),lda,work,n+nb+1)
540 
541 * Update L21
542  DO i=1,n-cut-nnb
543  DO j=1,nnb
544  a(cut+nnb+i,cut+j)=work(i,j)
545  END DO
546  END DO
547  ELSE
548 *
549 * L11 = L11**T*invD1*L11
550 *
551  DO i=1,nnb
552  DO j=1,i
553  a(cut+i,cut+j)=work(u11+i,j)
554  END DO
555  END DO
556  END IF
557 *
558 * Next Block
559 *
560  cut=cut+nnb
561  END DO
562 *
563 * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
564 *
565  i=n
566  DO WHILE ( i .GE. 1 )
567  IF( ipiv(i) .GT. 0 ) THEN
568  ip=ipiv(i)
569  IF (i .LT. ip) CALL zsyswapr( uplo, n, a, lda, i ,ip )
570  IF (i .GT. ip) CALL zsyswapr( uplo, n, a, lda, ip ,i )
571  ELSE
572  ip=-ipiv(i)
573  IF ( i .LT. ip) CALL zsyswapr( uplo, n, a, lda, i ,ip )
574  IF ( i .GT. ip) CALL zsyswapr( uplo, n, a, lda, ip ,i )
575  i=i-1
576  ENDIF
577  i=i-1
578  END DO
579  END IF
580 *
581  RETURN
582 *
583 * End of ZSYTRI2X
584 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM
Definition: ztrmm.f:177
subroutine ztrtri(UPLO, DIAG, N, A, LDA, INFO)
ZTRTRI
Definition: ztrtri.f:109
subroutine zsyswapr(UPLO, N, A, LDA, I1, I2)
ZSYSWAPR
Definition: zsyswapr.f:102
subroutine zsyconv(UPLO, WAY, N, A, LDA, IPIV, E, INFO)
ZSYCONV
Definition: zsyconv.f:114
Here is the call graph for this function:
Here is the caller graph for this function: