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