LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zsytri_rook()

subroutine zsytri_rook ( character  uplo,
integer  n,
complex*16, dimension( lda, * )  a,
integer  lda,
integer, dimension( * )  ipiv,
complex*16, dimension( * )  work,
integer  info 
)

ZSYTRI_ROOK

Download ZSYTRI_ROOK + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZSYTRI_ROOK computes the inverse of a complex symmetric
 matrix A using the factorization A = U*D*U**T or A = L*D*L**T
 computed by ZSYTRF_ROOK.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the details of the factorization are stored
          as an upper or lower triangular matrix.
          = 'U':  Upper triangular, form is A = U*D*U**T;
          = 'L':  Lower triangular, form is A = L*D*L**T.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the block diagonal matrix D and the multipliers
          used to obtain the factor U or L as computed by ZSYTRF_ROOK.

          On exit, if INFO = 0, the (symmetric) inverse of the original
          matrix.  If UPLO = 'U', the upper triangular part of the
          inverse is formed and the part of A below the diagonal is not
          referenced; if UPLO = 'L' the lower triangular part of the
          inverse is formed and the part of A above the diagonal is
          not referenced.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D
          as determined by ZSYTRF_ROOK.
[out]WORK
          WORK is COMPLEX*16 array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
               inverse could not be computed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
   December 2016, Igor Kozachenko,
                  Computer Science Division,
                  University of California, Berkeley

  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
                  School of Mathematics,
                  University of Manchester

Definition at line 128 of file zsytri_rook.f.

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*16 A( LDA, * ), WORK( * )
141* ..
142*
143* =====================================================================
144*
145* .. Parameters ..
146 COMPLEX*16 CONE, CZERO
147 parameter( cone = ( 1.0d+0, 0.0d+0 ),
148 $ czero = ( 0.0d+0, 0.0d+0 ) )
149* ..
150* .. Local Scalars ..
151 LOGICAL UPPER
152 INTEGER K, KP, KSTEP
153 COMPLEX*16 AK, AKKP1, AKP1, D, T, TEMP
154* ..
155* .. External Functions ..
156 LOGICAL LSAME
157 COMPLEX*16 ZDOTU
158 EXTERNAL lsame, zdotu
159* ..
160* .. External Subroutines ..
161 EXTERNAL zcopy, zswap, zsymv, 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( 'ZSYTRI_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 zcopy( k-1, a( 1, k ), 1, work, 1 )
237 CALL zsymv( uplo, k-1, -cone, a, lda, work, 1, czero,
238 $ a( 1, k ), 1 )
239 a( k, k ) = a( k, k ) - zdotu( 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 zcopy( k-1, a( 1, k ), 1, work, 1 )
262 CALL zsymv( uplo, k-1, -cone, a, lda, work, 1, czero,
263 $ a( 1, k ), 1 )
264 a( k, k ) = a( k, k ) - zdotu( k-1, work, 1, a( 1, k ),
265 $ 1 )
266 a( k, k+1 ) = a( k, k+1 ) -
267 $ zdotu( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
268 CALL zcopy( k-1, a( 1, k+1 ), 1, work, 1 )
269 CALL zsymv( 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 $ zdotu( 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 zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
286 CALL zswap( 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 zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
300 CALL zswap( 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 zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
315 CALL zswap( 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 zcopy( n-k, a( k+1, k ), 1, work, 1 )
353 CALL zsymv( 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 ) - zdotu( 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 zcopy( n-k, a( k+1, k ), 1, work, 1 )
378 CALL zsymv( 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 ) - zdotu( n-k, work, 1, a( k+1, k ),
381 $ 1 )
382 a( k, k-1 ) = a( k, k-1 ) -
383 $ zdotu( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
384 $ 1 )
385 CALL zcopy( n-k, a( k+1, k-1 ), 1, work, 1 )
386 CALL zsymv( 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 $ zdotu( 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 zswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
403 CALL zswap( 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 zswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
417 CALL zswap( 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 zswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
432 CALL zswap( 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 ZSYTRI_ROOK
447*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
complex *16 function zdotu(n, zx, incx, zy, incy)
ZDOTU
Definition zdotu.f:83
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:157
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: