LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
zsytri2x.f
Go to the documentation of this file.
1 *> \brief \b ZSYTRI2X
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZSYTRI2X + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytri2x.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytri2x.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytri2x.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER UPLO
25 * INTEGER INFO, LDA, N, NB
26 * ..
27 * .. Array Arguments ..
28 * INTEGER IPIV( * )
29 * COMPLEX*16 A( LDA, * ), WORK( N+NB+1,* )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> ZSYTRI2X computes the inverse of a complex symmetric indefinite matrix
39 *> A using the factorization A = U*D*U**T or A = L*D*L**T computed by
40 *> ZSYTRF.
41 *> \endverbatim
42 *
43 * Arguments:
44 * ==========
45 *
46 *> \param[in] UPLO
47 *> \verbatim
48 *> UPLO is CHARACTER*1
49 *> Specifies whether the details of the factorization are stored
50 *> as an upper or lower triangular matrix.
51 *> = 'U': Upper triangular, form is A = U*D*U**T;
52 *> = 'L': Lower triangular, form is A = L*D*L**T.
53 *> \endverbatim
54 *>
55 *> \param[in] N
56 *> \verbatim
57 *> N is INTEGER
58 *> The order of the matrix A. N >= 0.
59 *> \endverbatim
60 *>
61 *> \param[in,out] A
62 *> \verbatim
63 *> A is COMPLEX*16 array, dimension (LDA,N)
64 *> On entry, the NNB diagonal matrix D and the multipliers
65 *> used to obtain the factor U or L as computed by ZSYTRF.
66 *>
67 *> On exit, if INFO = 0, the (symmetric) inverse of the original
68 *> matrix. If UPLO = 'U', the upper triangular part of the
69 *> inverse is formed and the part of A below the diagonal is not
70 *> referenced; if UPLO = 'L' the lower triangular part of the
71 *> inverse is formed and the part of A above the diagonal is
72 *> not referenced.
73 *> \endverbatim
74 *>
75 *> \param[in] LDA
76 *> \verbatim
77 *> LDA is INTEGER
78 *> The leading dimension of the array A. LDA >= max(1,N).
79 *> \endverbatim
80 *>
81 *> \param[in] IPIV
82 *> \verbatim
83 *> IPIV is INTEGER array, dimension (N)
84 *> Details of the interchanges and the NNB structure of D
85 *> as determined by ZSYTRF.
86 *> \endverbatim
87 *>
88 *> \param[out] WORK
89 *> \verbatim
90 *> WORK is COMPLEX*16 array, dimension (N+NNB+1,NNB+3)
91 *> \endverbatim
92 *>
93 *> \param[in] NB
94 *> \verbatim
95 *> NB is INTEGER
96 *> Block size
97 *> \endverbatim
98 *>
99 *> \param[out] INFO
100 *> \verbatim
101 *> INFO is INTEGER
102 *> = 0: successful exit
103 *> < 0: if INFO = -i, the i-th argument had an illegal value
104 *> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
105 *> inverse could not be computed.
106 *> \endverbatim
107 *
108 * Authors:
109 * ========
110 *
111 *> \author Univ. of Tennessee
112 *> \author Univ. of California Berkeley
113 *> \author Univ. of Colorado Denver
114 *> \author NAG Ltd.
115 *
116 *> \date November 2011
117 *
118 *> \ingroup complex16SYcomputational
119 *
120 * =====================================================================
121  SUBROUTINE zsytri2x( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
122 *
123 * -- LAPACK computational routine (version 3.4.0) --
124 * -- LAPACK is a software package provided by Univ. of Tennessee, --
125 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
126 * November 2011
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 *
588  END
589 
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 zsytri2x(UPLO, N, A, LDA, IPIV, WORK, NB, INFO)
ZSYTRI2X
Definition: zsytri2x.f:122
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM
Definition: ztrmm.f:179
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