LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zhetri2x.f
Go to the documentation of this file.
1 *> \brief \b ZHETRI2X
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZHETRI2X + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetri2x.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetri2x.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetri2x.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHETRI2X( 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 *> ZHETRI2X computes the inverse of a COMPLEX*16 Hermitian indefinite matrix
39 *> A using the factorization A = U*D*U**H or A = L*D*L**H computed by
40 *> ZHETRF.
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**H;
52 *> = 'L': Lower triangular, form is A = L*D*L**H.
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 ZHETRF.
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 ZHETRF.
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 complex16HEcomputational
119 *
120 * =====================================================================
121  SUBROUTINE zhetri2x( 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  REAL one
141  COMPLEX*16 cone, zero
142  parameter( one = 1.0d+0,
143  $ cone = ( 1.0d+0, 0.0d+0 ),
144  $ zero = ( 0.0d+0, 0.0d+0 ) )
145 * ..
146 * .. Local Scalars ..
147  LOGICAL upper
148  INTEGER i, iinfo, ip, k, cut, nnb
149  INTEGER count
150  INTEGER j, u11, invd
151 
152  COMPLEX*16 ak, akkp1, akp1, d, t
153  COMPLEX*16 u01_i_j, u01_ip1_j
154  COMPLEX*16 u11_i_j, u11_ip1_j
155 * ..
156 * .. External Functions ..
157  LOGICAL lsame
158  EXTERNAL lsame
159 * ..
160 * .. External Subroutines ..
161  EXTERNAL zsyconv, xerbla, ztrtri
162  EXTERNAL zgemm, ztrmm, zheswapr
163 * ..
164 * .. Intrinsic Functions ..
165  INTRINSIC max
166 * ..
167 * .. Executable Statements ..
168 *
169 * Test the input parameters.
170 *
171  info = 0
172  upper = lsame( uplo, 'U' )
173  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
174  info = -1
175  ELSE IF( n.LT.0 ) THEN
176  info = -2
177  ELSE IF( lda.LT.max( 1, n ) ) THEN
178  info = -4
179  END IF
180 *
181 * Quick return if possible
182 *
183 *
184  IF( info.NE.0 ) THEN
185  CALL xerbla( 'ZHETRI2X', -info )
186  return
187  END IF
188  IF( n.EQ.0 )
189  $ return
190 *
191 * Convert A
192 * Workspace got Non-diag elements of D
193 *
194  CALL zsyconv( uplo, 'C', n, a, lda, ipiv, work, iinfo )
195 *
196 * Check that the diagonal matrix D is nonsingular.
197 *
198  IF( upper ) THEN
199 *
200 * Upper triangular storage: examine D from bottom to top
201 *
202  DO info = n, 1, -1
203  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
204  $ return
205  END DO
206  ELSE
207 *
208 * Lower triangular storage: examine D from top to bottom.
209 *
210  DO info = 1, n
211  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
212  $ return
213  END DO
214  END IF
215  info = 0
216 *
217 * Splitting Workspace
218 * U01 is a block (N,NB+1)
219 * The first element of U01 is in WORK(1,1)
220 * U11 is a block (NB+1,NB+1)
221 * The first element of U11 is in WORK(N+1,1)
222  u11 = n
223 * INVD is a block (N,2)
224 * The first element of INVD is in WORK(1,INVD)
225  invd = nb+2
226 
227  IF( upper ) THEN
228 *
229 * invA = P * inv(U**H)*inv(D)*inv(U)*P**H.
230 *
231  CALL ztrtri( uplo, 'U', n, a, lda, info )
232 *
233 * inv(D) and inv(D)*inv(U)
234 *
235  k=1
236  DO WHILE ( k .LE. n )
237  IF( ipiv( k ).GT.0 ) THEN
238 * 1 x 1 diagonal NNB
239  work(k,invd) = one / REAL ( A( K, K ) )
240  work(k,invd+1) = 0
241  k=k+1
242  ELSE
243 * 2 x 2 diagonal NNB
244  t = abs( work(k+1,1) )
245  ak = REAL ( A( K, K ) ) / t
246  akp1 = REAL ( A( K+1, K+1 ) ) / t
247  akkp1 = work(k+1,1) / t
248  d = t*( ak*akp1-one )
249  work(k,invd) = akp1 / d
250  work(k+1,invd+1) = ak / d
251  work(k,invd+1) = -akkp1 / d
252  work(k+1,invd) = dconjg(work(k,invd+1) )
253  k=k+2
254  END IF
255  END DO
256 *
257 * inv(U**H) = (inv(U))**H
258 *
259 * inv(U**H)*inv(D)*inv(U)
260 *
261  cut=n
262  DO WHILE (cut .GT. 0)
263  nnb=nb
264  IF (cut .LE. nnb) THEN
265  nnb=cut
266  ELSE
267  count = 0
268 * count negative elements,
269  DO i=cut+1-nnb,cut
270  IF (ipiv(i) .LT. 0) count=count+1
271  END DO
272 * need a even number for a clear cut
273  IF (mod(count,2) .EQ. 1) nnb=nnb+1
274  END IF
275 
276  cut=cut-nnb
277 *
278 * U01 Block
279 *
280  DO i=1,cut
281  DO j=1,nnb
282  work(i,j)=a(i,cut+j)
283  END DO
284  END DO
285 *
286 * U11 Block
287 *
288  DO i=1,nnb
289  work(u11+i,i)=cone
290  DO j=1,i-1
291  work(u11+i,j)=zero
292  END DO
293  DO j=i+1,nnb
294  work(u11+i,j)=a(cut+i,cut+j)
295  END DO
296  END DO
297 *
298 * invD*U01
299 *
300  i=1
301  DO WHILE (i .LE. cut)
302  IF (ipiv(i) > 0) THEN
303  DO j=1,nnb
304  work(i,j)=work(i,invd)*work(i,j)
305  END DO
306  i=i+1
307  ELSE
308  DO j=1,nnb
309  u01_i_j = work(i,j)
310  u01_ip1_j = work(i+1,j)
311  work(i,j)=work(i,invd)*u01_i_j+
312  $ work(i,invd+1)*u01_ip1_j
313  work(i+1,j)=work(i+1,invd)*u01_i_j+
314  $ work(i+1,invd+1)*u01_ip1_j
315  END DO
316  i=i+2
317  END IF
318  END DO
319 *
320 * invD1*U11
321 *
322  i=1
323  DO WHILE (i .LE. nnb)
324  IF (ipiv(cut+i) > 0) THEN
325  DO j=i,nnb
326  work(u11+i,j)=work(cut+i,invd)*work(u11+i,j)
327  END DO
328  i=i+1
329  ELSE
330  DO j=i,nnb
331  u11_i_j = work(u11+i,j)
332  u11_ip1_j = work(u11+i+1,j)
333  work(u11+i,j)=work(cut+i,invd)*work(u11+i,j) +
334  $ work(cut+i,invd+1)*work(u11+i+1,j)
335  work(u11+i+1,j)=work(cut+i+1,invd)*u11_i_j+
336  $ work(cut+i+1,invd+1)*u11_ip1_j
337  END DO
338  i=i+2
339  END IF
340  END DO
341 *
342 * U11**H*invD1*U11->U11
343 *
344  CALL ztrmm('L','U','C','U',nnb, nnb,
345  $ cone,a(cut+1,cut+1),lda,work(u11+1,1),n+nb+1)
346 *
347  DO i=1,nnb
348  DO j=i,nnb
349  a(cut+i,cut+j)=work(u11+i,j)
350  END DO
351  END DO
352 *
353 * U01**H*invD*U01->A(CUT+I,CUT+J)
354 *
355  CALL zgemm('C','N',nnb,nnb,cut,cone,a(1,cut+1),lda,
356  $ work,n+nb+1, zero, work(u11+1,1), n+nb+1)
357 *
358 * U11 = U11**H*invD1*U11 + U01**H*invD*U01
359 *
360  DO i=1,nnb
361  DO j=i,nnb
362  a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
363  END DO
364  END DO
365 *
366 * U01 = U00**H*invD0*U01
367 *
368  CALL ztrmm('L',uplo,'C','U',cut, nnb,
369  $ cone,a,lda,work,n+nb+1)
370 
371 *
372 * Update U01
373 *
374  DO i=1,cut
375  DO j=1,nnb
376  a(i,cut+j)=work(i,j)
377  END DO
378  END DO
379 *
380 * Next Block
381 *
382  END DO
383 *
384 * Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H
385 *
386  i=1
387  DO WHILE ( i .LE. n )
388  IF( ipiv(i) .GT. 0 ) THEN
389  ip=ipiv(i)
390  IF (i .LT. ip) CALL zheswapr( uplo, n, a, lda, i ,ip )
391  IF (i .GT. ip) CALL zheswapr( uplo, n, a, lda, ip ,i )
392  ELSE
393  ip=-ipiv(i)
394  i=i+1
395  IF ( (i-1) .LT. ip)
396  $ CALL zheswapr( uplo, n, a, lda, i-1 ,ip )
397  IF ( (i-1) .GT. ip)
398  $ CALL zheswapr( uplo, n, a, lda, ip ,i-1 )
399  ENDIF
400  i=i+1
401  END DO
402  ELSE
403 *
404 * LOWER...
405 *
406 * invA = P * inv(U**H)*inv(D)*inv(U)*P**H.
407 *
408  CALL ztrtri( uplo, 'U', n, a, lda, info )
409 *
410 * inv(D) and inv(D)*inv(U)
411 *
412  k=n
413  DO WHILE ( k .GE. 1 )
414  IF( ipiv( k ).GT.0 ) THEN
415 * 1 x 1 diagonal NNB
416  work(k,invd) = one / REAL ( A( K, K ) )
417  work(k,invd+1) = 0
418  k=k-1
419  ELSE
420 * 2 x 2 diagonal NNB
421  t = abs( work(k-1,1) )
422  ak = REAL ( A( K-1, K-1 ) ) / t
423  akp1 = REAL ( A( K, K ) ) / t
424  akkp1 = work(k-1,1) / t
425  d = t*( ak*akp1-one )
426  work(k-1,invd) = akp1 / d
427  work(k,invd) = ak / d
428  work(k,invd+1) = -akkp1 / d
429  work(k-1,invd+1) = dconjg(work(k,invd+1) )
430  k=k-2
431  END IF
432  END DO
433 *
434 * inv(U**H) = (inv(U))**H
435 *
436 * inv(U**H)*inv(D)*inv(U)
437 *
438  cut=0
439  DO WHILE (cut .LT. n)
440  nnb=nb
441  IF (cut + nnb .GE. n) THEN
442  nnb=n-cut
443  ELSE
444  count = 0
445 * count negative elements,
446  DO i=cut+1,cut+nnb
447  IF (ipiv(i) .LT. 0) count=count+1
448  END DO
449 * need a even number for a clear cut
450  IF (mod(count,2) .EQ. 1) nnb=nnb+1
451  END IF
452 * L21 Block
453  DO i=1,n-cut-nnb
454  DO j=1,nnb
455  work(i,j)=a(cut+nnb+i,cut+j)
456  END DO
457  END DO
458 * L11 Block
459  DO i=1,nnb
460  work(u11+i,i)=cone
461  DO j=i+1,nnb
462  work(u11+i,j)=zero
463  END DO
464  DO j=1,i-1
465  work(u11+i,j)=a(cut+i,cut+j)
466  END DO
467  END DO
468 *
469 * invD*L21
470 *
471  i=n-cut-nnb
472  DO WHILE (i .GE. 1)
473  IF (ipiv(cut+nnb+i) > 0) THEN
474  DO j=1,nnb
475  work(i,j)=work(cut+nnb+i,invd)*work(i,j)
476  END DO
477  i=i-1
478  ELSE
479  DO j=1,nnb
480  u01_i_j = work(i,j)
481  u01_ip1_j = work(i-1,j)
482  work(i,j)=work(cut+nnb+i,invd)*u01_i_j+
483  $ work(cut+nnb+i,invd+1)*u01_ip1_j
484  work(i-1,j)=work(cut+nnb+i-1,invd+1)*u01_i_j+
485  $ work(cut+nnb+i-1,invd)*u01_ip1_j
486  END DO
487  i=i-2
488  END IF
489  END DO
490 *
491 * invD1*L11
492 *
493  i=nnb
494  DO WHILE (i .GE. 1)
495  IF (ipiv(cut+i) > 0) THEN
496  DO j=1,nnb
497  work(u11+i,j)=work(cut+i,invd)*work(u11+i,j)
498  END DO
499  i=i-1
500  ELSE
501  DO j=1,nnb
502  u11_i_j = work(u11+i,j)
503  u11_ip1_j = work(u11+i-1,j)
504  work(u11+i,j)=work(cut+i,invd)*work(u11+i,j) +
505  $ work(cut+i,invd+1)*u11_ip1_j
506  work(u11+i-1,j)=work(cut+i-1,invd+1)*u11_i_j+
507  $ work(cut+i-1,invd)*u11_ip1_j
508  END DO
509  i=i-2
510  END IF
511  END DO
512 *
513 * L11**H*invD1*L11->L11
514 *
515  CALL ztrmm('L',uplo,'C','U',nnb, nnb,
516  $ cone,a(cut+1,cut+1),lda,work(u11+1,1),n+nb+1)
517 *
518  DO i=1,nnb
519  DO j=1,i
520  a(cut+i,cut+j)=work(u11+i,j)
521  END DO
522  END DO
523 *
524  IF ( (cut+nnb) .LT. n ) THEN
525 *
526 * L21**H*invD2*L21->A(CUT+I,CUT+J)
527 *
528  CALL zgemm('C','N',nnb,nnb,n-nnb-cut,cone,a(cut+nnb+1,cut+1)
529  $ ,lda,work,n+nb+1, zero, work(u11+1,1), n+nb+1)
530 
531 *
532 * L11 = L11**H*invD1*L11 + U01**H*invD*U01
533 *
534  DO i=1,nnb
535  DO j=1,i
536  a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
537  END DO
538  END DO
539 *
540 * L01 = L22**H*invD2*L21
541 *
542  CALL ztrmm('L',uplo,'C','U', n-nnb-cut, nnb,
543  $ cone,a(cut+nnb+1,cut+nnb+1),lda,work,n+nb+1)
544 
545 * Update L21
546  DO i=1,n-cut-nnb
547  DO j=1,nnb
548  a(cut+nnb+i,cut+j)=work(i,j)
549  END DO
550  END DO
551  ELSE
552 *
553 * L11 = L11**H*invD1*L11
554 *
555  DO i=1,nnb
556  DO j=1,i
557  a(cut+i,cut+j)=work(u11+i,j)
558  END DO
559  END DO
560  END IF
561 *
562 * Next Block
563 *
564  cut=cut+nnb
565  END DO
566 *
567 * Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H
568 *
569  i=n
570  DO WHILE ( i .GE. 1 )
571  IF( ipiv(i) .GT. 0 ) THEN
572  ip=ipiv(i)
573  IF (i .LT. ip) CALL zheswapr( uplo, n, a, lda, i ,ip )
574  IF (i .GT. ip) CALL zheswapr( uplo, n, a, lda, ip ,i )
575  ELSE
576  ip=-ipiv(i)
577  IF ( i .LT. ip) CALL zheswapr( uplo, n, a, lda, i ,ip )
578  IF ( i .GT. ip) CALL zheswapr( uplo, n, a, lda, ip ,i )
579  i=i-1
580  ENDIF
581  i=i-1
582  END DO
583  END IF
584 *
585  return
586 *
587 * End of ZHETRI2X
588 *
589  END
590