LAPACK  3.10.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+NB+1,NB+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 *> \ingroup complex16SYcomputational
117 *
118 * =====================================================================
119  SUBROUTINE zsytri2x( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
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 *
585  END
586 
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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
subroutine zsytri2x(UPLO, N, A, LDA, IPIV, WORK, NB, INFO)
ZSYTRI2X
Definition: zsytri2x.f:120