LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
zsytri_rook.f
Go to the documentation of this file.
1 *> \brief \b ZSYTRI_ROOK
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZSYTRI_ROOK + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytri_rook.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytri_rook.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytri_rook.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZSYTRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER UPLO
25 * INTEGER INFO, LDA, N
26 * ..
27 * .. Array Arguments ..
28 * INTEGER IPIV( * )
29 * COMPLEX*16 A( LDA, * ), WORK( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> ZSYTRI_ROOK computes the inverse of a complex symmetric
39 *> matrix A using the factorization A = U*D*U**T or A = L*D*L**T
40 *> computed by ZSYTRF_ROOK.
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 block diagonal matrix D and the multipliers
65 *> used to obtain the factor U or L as computed by ZSYTRF_ROOK.
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 block structure of D
85 *> as determined by ZSYTRF_ROOK.
86 *> \endverbatim
87 *>
88 *> \param[out] WORK
89 *> \verbatim
90 *> WORK is COMPLEX*16 array, dimension (N)
91 *> \endverbatim
92 *>
93 *> \param[out] INFO
94 *> \verbatim
95 *> INFO is INTEGER
96 *> = 0: successful exit
97 *> < 0: if INFO = -i, the i-th argument had an illegal value
98 *> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
99 *> inverse could not be computed.
100 *> \endverbatim
101 *
102 * Authors:
103 * ========
104 *
105 *> \author Univ. of Tennessee
106 *> \author Univ. of California Berkeley
107 *> \author Univ. of Colorado Denver
108 *> \author NAG Ltd.
109 *
110 *> \date December 2016
111 *
112 *> \ingroup complex16SYcomputational
113 *
114 *> \par Contributors:
115 * ==================
116 *>
117 *> \verbatim
118 *>
119 *> December 2016, Igor Kozachenko,
120 *> Computer Science Division,
121 *> University of California, Berkeley
122 *>
123 *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
124 *> School of Mathematics,
125 *> University of Manchester
126 *>
127 *> \endverbatim
128 *
129 * =====================================================================
130  SUBROUTINE zsytri_rook( UPLO, N, A, LDA, IPIV, WORK, INFO )
131 *
132 * -- LAPACK computational routine (version 3.7.0) --
133 * -- LAPACK is a software package provided by Univ. of Tennessee, --
134 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
135 * December 2016
136 *
137 * .. Scalar Arguments ..
138  CHARACTER UPLO
139  INTEGER INFO, LDA, N
140 * ..
141 * .. Array Arguments ..
142  INTEGER IPIV( * )
143  COMPLEX*16 A( LDA, * ), WORK( * )
144 * ..
145 *
146 * =====================================================================
147 *
148 * .. Parameters ..
149  COMPLEX*16 CONE, CZERO
150  parameter( cone = ( 1.0d+0, 0.0d+0 ),
151  $ czero = ( 0.0d+0, 0.0d+0 ) )
152 * ..
153 * .. Local Scalars ..
154  LOGICAL UPPER
155  INTEGER K, KP, KSTEP
156  COMPLEX*16 AK, AKKP1, AKP1, D, T, TEMP
157 * ..
158 * .. External Functions ..
159  LOGICAL LSAME
160  COMPLEX*16 ZDOTU
161  EXTERNAL lsame, zdotu
162 * ..
163 * .. External Subroutines ..
164  EXTERNAL zcopy, zswap, zsymv, xerbla
165 * ..
166 * .. Intrinsic Functions ..
167  INTRINSIC max
168 * ..
169 * .. Executable Statements ..
170 *
171 * Test the input parameters.
172 *
173  info = 0
174  upper = lsame( uplo, 'U' )
175  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
176  info = -1
177  ELSE IF( n.LT.0 ) THEN
178  info = -2
179  ELSE IF( lda.LT.max( 1, n ) ) THEN
180  info = -4
181  END IF
182  IF( info.NE.0 ) THEN
183  CALL xerbla( 'ZSYTRI_ROOK', -info )
184  RETURN
185  END IF
186 *
187 * Quick return if possible
188 *
189  IF( n.EQ.0 )
190  $ RETURN
191 *
192 * Check that the diagonal matrix D is nonsingular.
193 *
194  IF( upper ) THEN
195 *
196 * Upper triangular storage: examine D from bottom to top
197 *
198  DO 10 info = n, 1, -1
199  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.czero )
200  $ RETURN
201  10 CONTINUE
202  ELSE
203 *
204 * Lower triangular storage: examine D from top to bottom.
205 *
206  DO 20 info = 1, n
207  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.czero )
208  $ RETURN
209  20 CONTINUE
210  END IF
211  info = 0
212 *
213  IF( upper ) THEN
214 *
215 * Compute inv(A) from the factorization A = U*D*U**T.
216 *
217 * K is the main loop index, increasing from 1 to N in steps of
218 * 1 or 2, depending on the size of the diagonal blocks.
219 *
220  k = 1
221  30 CONTINUE
222 *
223 * If K > N, exit from loop.
224 *
225  IF( k.GT.n )
226  $ GO TO 40
227 *
228  IF( ipiv( k ).GT.0 ) THEN
229 *
230 * 1 x 1 diagonal block
231 *
232 * Invert the diagonal block.
233 *
234  a( k, k ) = cone / a( k, k )
235 *
236 * Compute column K of the inverse.
237 *
238  IF( k.GT.1 ) THEN
239  CALL zcopy( k-1, a( 1, k ), 1, work, 1 )
240  CALL zsymv( uplo, k-1, -cone, a, lda, work, 1, czero,
241  $ a( 1, k ), 1 )
242  a( k, k ) = a( k, k ) - zdotu( k-1, work, 1, a( 1, k ),
243  $ 1 )
244  END IF
245  kstep = 1
246  ELSE
247 *
248 * 2 x 2 diagonal block
249 *
250 * Invert the diagonal block.
251 *
252  t = a( k, k+1 )
253  ak = a( k, k ) / t
254  akp1 = a( k+1, k+1 ) / t
255  akkp1 = a( k, k+1 ) / t
256  d = t*( ak*akp1-cone )
257  a( k, k ) = akp1 / d
258  a( k+1, k+1 ) = ak / d
259  a( k, k+1 ) = -akkp1 / d
260 *
261 * Compute columns K and K+1 of the inverse.
262 *
263  IF( k.GT.1 ) THEN
264  CALL zcopy( k-1, a( 1, k ), 1, work, 1 )
265  CALL zsymv( uplo, k-1, -cone, a, lda, work, 1, czero,
266  $ a( 1, k ), 1 )
267  a( k, k ) = a( k, k ) - zdotu( k-1, work, 1, a( 1, k ),
268  $ 1 )
269  a( k, k+1 ) = a( k, k+1 ) -
270  $ zdotu( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
271  CALL zcopy( k-1, a( 1, k+1 ), 1, work, 1 )
272  CALL zsymv( uplo, k-1, -cone, a, lda, work, 1, czero,
273  $ a( 1, k+1 ), 1 )
274  a( k+1, k+1 ) = a( k+1, k+1 ) -
275  $ zdotu( k-1, work, 1, a( 1, k+1 ), 1 )
276  END IF
277  kstep = 2
278  END IF
279 *
280  IF( kstep.EQ.1 ) THEN
281 *
282 * Interchange rows and columns K and IPIV(K) in the leading
283 * submatrix A(1:k+1,1:k+1)
284 *
285  kp = ipiv( k )
286  IF( kp.NE.k ) THEN
287  IF( kp.GT.1 )
288  $ CALL zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
289  CALL zswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
290  temp = a( k, k )
291  a( k, k ) = a( kp, kp )
292  a( kp, kp ) = temp
293  END IF
294  ELSE
295 *
296 * Interchange rows and columns K and K+1 with -IPIV(K) and
297 * -IPIV(K+1)in the leading submatrix A(1:k+1,1:k+1)
298 *
299  kp = -ipiv( k )
300  IF( kp.NE.k ) THEN
301  IF( kp.GT.1 )
302  $ CALL zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
303  CALL zswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
304 *
305  temp = a( k, k )
306  a( k, k ) = a( kp, kp )
307  a( kp, kp ) = temp
308  temp = a( k, k+1 )
309  a( k, k+1 ) = a( kp, k+1 )
310  a( kp, k+1 ) = temp
311  END IF
312 *
313  k = k + 1
314  kp = -ipiv( k )
315  IF( kp.NE.k ) THEN
316  IF( kp.GT.1 )
317  $ CALL zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
318  CALL zswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
319  temp = a( k, k )
320  a( k, k ) = a( kp, kp )
321  a( kp, kp ) = temp
322  END IF
323  END IF
324 *
325  k = k + 1
326  GO TO 30
327  40 CONTINUE
328 *
329  ELSE
330 *
331 * Compute inv(A) from the factorization A = L*D*L**T.
332 *
333 * K is the main loop index, increasing from 1 to N in steps of
334 * 1 or 2, depending on the size of the diagonal blocks.
335 *
336  k = n
337  50 CONTINUE
338 *
339 * If K < 1, exit from loop.
340 *
341  IF( k.LT.1 )
342  $ GO TO 60
343 *
344  IF( ipiv( k ).GT.0 ) THEN
345 *
346 * 1 x 1 diagonal block
347 *
348 * Invert the diagonal block.
349 *
350  a( k, k ) = cone / a( k, k )
351 *
352 * Compute column K of the inverse.
353 *
354  IF( k.LT.n ) THEN
355  CALL zcopy( n-k, a( k+1, k ), 1, work, 1 )
356  CALL zsymv( uplo, n-k,-cone, a( k+1, k+1 ), lda, work, 1,
357  $ czero, a( k+1, k ), 1 )
358  a( k, k ) = a( k, k ) - zdotu( n-k, work, 1, a( k+1, k ),
359  $ 1 )
360  END IF
361  kstep = 1
362  ELSE
363 *
364 * 2 x 2 diagonal block
365 *
366 * Invert the diagonal block.
367 *
368  t = a( k, k-1 )
369  ak = a( k-1, k-1 ) / t
370  akp1 = a( k, k ) / t
371  akkp1 = a( k, k-1 ) / t
372  d = t*( ak*akp1-cone )
373  a( k-1, k-1 ) = akp1 / d
374  a( k, k ) = ak / d
375  a( k, k-1 ) = -akkp1 / d
376 *
377 * Compute columns K-1 and K of the inverse.
378 *
379  IF( k.LT.n ) THEN
380  CALL zcopy( n-k, a( k+1, k ), 1, work, 1 )
381  CALL zsymv( uplo, n-k,-cone, a( k+1, k+1 ), lda, work, 1,
382  $ czero, a( k+1, k ), 1 )
383  a( k, k ) = a( k, k ) - zdotu( n-k, work, 1, a( k+1, k ),
384  $ 1 )
385  a( k, k-1 ) = a( k, k-1 ) -
386  $ zdotu( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
387  $ 1 )
388  CALL zcopy( n-k, a( k+1, k-1 ), 1, work, 1 )
389  CALL zsymv( uplo, n-k,-cone, a( k+1, k+1 ), lda, work, 1,
390  $ czero, a( k+1, k-1 ), 1 )
391  a( k-1, k-1 ) = a( k-1, k-1 ) -
392  $ zdotu( n-k, work, 1, a( k+1, k-1 ), 1 )
393  END IF
394  kstep = 2
395  END IF
396 *
397  IF( kstep.EQ.1 ) THEN
398 *
399 * Interchange rows and columns K and IPIV(K) in the trailing
400 * submatrix A(k-1:n,k-1:n)
401 *
402  kp = ipiv( k )
403  IF( kp.NE.k ) THEN
404  IF( kp.LT.n )
405  $ CALL zswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
406  CALL zswap( kp-k-1, a( k+1, k ), 1, a( kp, k+1 ), lda )
407  temp = a( k, k )
408  a( k, k ) = a( kp, kp )
409  a( kp, kp ) = temp
410  END IF
411  ELSE
412 *
413 * Interchange rows and columns K and K-1 with -IPIV(K) and
414 * -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
415 *
416  kp = -ipiv( k )
417  IF( kp.NE.k ) THEN
418  IF( kp.LT.n )
419  $ CALL zswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
420  CALL zswap( kp-k-1, a( k+1, k ), 1, a( kp, k+1 ), lda )
421 *
422  temp = a( k, k )
423  a( k, k ) = a( kp, kp )
424  a( kp, kp ) = temp
425  temp = a( k, k-1 )
426  a( k, k-1 ) = a( kp, k-1 )
427  a( kp, k-1 ) = temp
428  END IF
429 *
430  k = k - 1
431  kp = -ipiv( k )
432  IF( kp.NE.k ) THEN
433  IF( kp.LT.n )
434  $ CALL zswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
435  CALL zswap( kp-k-1, a( k+1, k ), 1, a( kp, k+1 ), lda )
436  temp = a( k, k )
437  a( k, k ) = a( kp, kp )
438  a( kp, kp ) = temp
439  END IF
440  END IF
441 *
442  k = k - 1
443  GO TO 50
444  60 CONTINUE
445  END IF
446 *
447  RETURN
448 *
449 * End of ZSYTRI_ROOK
450 *
451  END
zsytri_rook
subroutine zsytri_rook(UPLO, N, A, LDA, IPIV, WORK, INFO)
ZSYTRI_ROOK
Definition: zsytri_rook.f:131
zsymv
subroutine zsymv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZSYMV computes a matrix-vector product for a complex symmetric matrix.
Definition: zsymv.f:159
zcopy
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:83
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
zswap
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:83