LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
Collaboration diagram for complex16:

Functions

subroutine zla_syamv (UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
 ZLA_SYAMV computes a matrix-vector product using a symmetric indefinite matrix to calculate error bounds. More...
 
double precision function zla_syrcond_c (UPLO, N, A, LDA, AF, LDAF, IPIV, C, CAPPLY, INFO, WORK, RWORK)
 ZLA_SYRCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for symmetric indefinite matrices. More...
 
double precision function zla_syrcond_x (UPLO, N, A, LDA, AF, LDAF, IPIV, X, INFO, WORK, RWORK)
 ZLA_SYRCOND_X computes the infinity norm condition number of op(A)*diag(x) for symmetric indefinite matrices. More...
 
subroutine zla_syrfsx_extended (PREC_TYPE, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, COLEQU, C, B, LDB, Y, LDY, BERR_OUT, N_NORMS, ERR_BNDS_NORM, ERR_BNDS_COMP, RES, AYB, DY, Y_TAIL, RCOND, ITHRESH, RTHRESH, DZ_UB, IGNORE_CWISE, INFO)
 ZLA_SYRFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric indefinite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution. More...
 
double precision function zla_syrpvgrw (UPLO, N, INFO, A, LDA, AF, LDAF, IPIV, WORK)
 ZLA_SYRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a symmetric indefinite matrix. More...
 
subroutine zlasyf (UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
 ZLASYF computes a partial factorization of a complex symmetric matrix using the Bunch-Kaufman diagonal pivoting method. More...
 
subroutine zlasyf_rook (UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
 ZLASYF_ROOK computes a partial factorization of a complex symmetric matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method. More...
 
subroutine zsycon (UPLO, N, A, LDA, IPIV, ANORM, RCOND, WORK, INFO)
 ZSYCON More...
 
subroutine zsycon_rook (UPLO, N, A, LDA, IPIV, ANORM, RCOND, WORK, INFO)
 ZSYCON_ROOK More...
 
subroutine zsyconv (UPLO, WAY, N, A, LDA, IPIV, E, INFO)
 ZSYCONV More...
 
subroutine zsyequb (UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO)
 ZSYEQUB More...
 
subroutine zsyrfs (UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
 ZSYRFS More...
 
subroutine zsyrfsx (UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV, S, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK, INFO)
 ZSYRFSX More...
 
subroutine zsytf2 (UPLO, N, A, LDA, IPIV, INFO)
 ZSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting method (unblocked algorithm). More...
 
subroutine zsytf2_rook (UPLO, N, A, LDA, IPIV, INFO)
 ZSYTF2_ROOK computes the factorization of a complex symmetric indefinite matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method (unblocked algorithm). More...
 
subroutine zsytrf (UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO)
 ZSYTRF More...
 
subroutine zsytrf_rook (UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO)
 ZSYTRF_ROOK More...
 
subroutine zsytri (UPLO, N, A, LDA, IPIV, WORK, INFO)
 ZSYTRI More...
 
subroutine zsytri2 (UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO)
 ZSYTRI2 More...
 
subroutine zsytri2x (UPLO, N, A, LDA, IPIV, WORK, NB, INFO)
 ZSYTRI2X More...
 
subroutine zsytri_rook (UPLO, N, A, LDA, IPIV, WORK, INFO)
 ZSYTRI_ROOK More...
 
subroutine zsytrs (UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
 ZSYTRS More...
 
subroutine zsytrs2 (UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, INFO)
 ZSYTRS2 More...
 
subroutine zsytrs_rook (UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
 ZSYTRS_ROOK More...
 
subroutine ztgsyl (TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)
 ZTGSYL More...
 
subroutine ztrsyl (TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, LDC, SCALE, INFO)
 ZTRSYL More...
 

Detailed Description

This is the group of complex16 computational functions for SY matrices

Function Documentation

subroutine zla_syamv ( integer  UPLO,
integer  N,
double precision  ALPHA,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  X,
integer  INCX,
double precision  BETA,
double precision, dimension( * )  Y,
integer  INCY 
)

ZLA_SYAMV computes a matrix-vector product using a symmetric indefinite matrix to calculate error bounds.

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

Purpose:
 ZLA_SYAMV  performs the matrix-vector operation

         y := alpha*abs(A)*abs(x) + beta*abs(y),

 where alpha and beta are scalars, x and y are vectors and A is an
 n by n symmetric matrix.

 This function is primarily used in calculating error bounds.
 To protect against underflow during evaluation, components in
 the resulting vector are perturbed away from zero by (N+1)
 times the underflow threshold.  To prevent unnecessarily large
 errors for block-structure embedded in general matrices,
 "symbolically" zero components are not perturbed.  A zero
 entry is considered "symbolic" if all multiplications involved
 in computing that entry have at least one zero multiplicand.
Parameters
[in]UPLO
          UPLO is INTEGER
           On entry, UPLO specifies whether the upper or lower
           triangular part of the array A is to be referenced as
           follows:

              UPLO = BLAS_UPPER   Only the upper triangular part of A
                                  is to be referenced.

              UPLO = BLAS_LOWER   Only the lower triangular part of A
                                  is to be referenced.

           Unchanged on exit.
[in]N
          N is INTEGER
           On entry, N specifies the number of columns of the matrix A.
           N must be at least zero.
           Unchanged on exit.
[in]ALPHA
          ALPHA is DOUBLE PRECISION .
           On entry, ALPHA specifies the scalar alpha.
           Unchanged on exit.
[in]A
          A is COMPLEX*16 array, DIMENSION ( LDA, n ).
           Before entry, the leading m by n part of the array A must
           contain the matrix of coefficients.
           Unchanged on exit.
[in]LDA
          LDA is INTEGER
           On entry, LDA specifies the first dimension of A as declared
           in the calling (sub) program. LDA must be at least
           max( 1, n ).
           Unchanged on exit.
[in]X
          X is COMPLEX*16 array, DIMENSION at least
           ( 1 + ( n - 1 )*abs( INCX ) )
           Before entry, the incremented array X must contain the
           vector x.
           Unchanged on exit.
[in]INCX
          INCX is INTEGER
           On entry, INCX specifies the increment for the elements of
           X. INCX must not be zero.
           Unchanged on exit.
[in]BETA
          BETA is DOUBLE PRECISION .
           On entry, BETA specifies the scalar beta. When BETA is
           supplied as zero then Y need not be set on input.
           Unchanged on exit.
[in,out]Y
          Y is DOUBLE PRECISION array, dimension
           ( 1 + ( n - 1 )*abs( INCY ) )
           Before entry with BETA non-zero, the incremented array Y
           must contain the vector y. On exit, Y is overwritten by the
           updated vector y.
[in]INCY
          INCY is INTEGER
           On entry, INCY specifies the increment for the elements of
           Y. INCY must not be zero.
           Unchanged on exit.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
  Level 2 Blas routine.

  -- Written on 22-October-1986.
     Jack Dongarra, Argonne National Lab.
     Jeremy Du Croz, Nag Central Office.
     Sven Hammarling, Nag Central Office.
     Richard Hanson, Sandia National Labs.
  -- Modified for the absolute-value product, April 2006
     Jason Riedy, UC Berkeley

Definition at line 181 of file zla_syamv.f.

181 *
182 * -- LAPACK computational routine (version 3.4.2) --
183 * -- LAPACK is a software package provided by Univ. of Tennessee, --
184 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
185 * September 2012
186 *
187 * .. Scalar Arguments ..
188  DOUBLE PRECISION alpha, beta
189  INTEGER incx, incy, lda, n
190  INTEGER uplo
191 * ..
192 * .. Array Arguments ..
193  COMPLEX*16 a( lda, * ), x( * )
194  DOUBLE PRECISION y( * )
195 * ..
196 *
197 * =====================================================================
198 *
199 * .. Parameters ..
200  DOUBLE PRECISION one, zero
201  parameter( one = 1.0d+0, zero = 0.0d+0 )
202 * ..
203 * .. Local Scalars ..
204  LOGICAL symb_zero
205  DOUBLE PRECISION temp, safe1
206  INTEGER i, info, iy, j, jx, kx, ky
207  COMPLEX*16 zdum
208 * ..
209 * .. External Subroutines ..
210  EXTERNAL xerbla, dlamch
211  DOUBLE PRECISION dlamch
212 * ..
213 * .. External Functions ..
214  EXTERNAL ilauplo
215  INTEGER ilauplo
216 * ..
217 * .. Intrinsic Functions ..
218  INTRINSIC max, abs, sign, REAL, dimag
219 * ..
220 * .. Statement Functions ..
221  DOUBLE PRECISION cabs1
222 * ..
223 * .. Statement Function Definitions ..
224  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
225 * ..
226 * .. Executable Statements ..
227 *
228 * Test the input parameters.
229 *
230  info = 0
231  IF ( uplo.NE.ilauplo( 'U' ) .AND.
232  $ uplo.NE.ilauplo( 'L' ) )THEN
233  info = 1
234  ELSE IF( n.LT.0 )THEN
235  info = 2
236  ELSE IF( lda.LT.max( 1, n ) )THEN
237  info = 5
238  ELSE IF( incx.EQ.0 )THEN
239  info = 7
240  ELSE IF( incy.EQ.0 )THEN
241  info = 10
242  END IF
243  IF( info.NE.0 )THEN
244  CALL xerbla( 'DSYMV ', info )
245  RETURN
246  END IF
247 *
248 * Quick return if possible.
249 *
250  IF( ( n.EQ.0 ).OR.( ( alpha.EQ.zero ).AND.( beta.EQ.one ) ) )
251  $ RETURN
252 *
253 * Set up the start points in X and Y.
254 *
255  IF( incx.GT.0 )THEN
256  kx = 1
257  ELSE
258  kx = 1 - ( n - 1 )*incx
259  END IF
260  IF( incy.GT.0 )THEN
261  ky = 1
262  ELSE
263  ky = 1 - ( n - 1 )*incy
264  END IF
265 *
266 * Set SAFE1 essentially to be the underflow threshold times the
267 * number of additions in each row.
268 *
269  safe1 = dlamch( 'Safe minimum' )
270  safe1 = (n+1)*safe1
271 *
272 * Form y := alpha*abs(A)*abs(x) + beta*abs(y).
273 *
274 * The O(N^2) SYMB_ZERO tests could be replaced by O(N) queries to
275 * the inexact flag. Still doesn't help change the iteration order
276 * to per-column.
277 *
278  iy = ky
279  IF ( incx.EQ.1 ) THEN
280  IF ( uplo .EQ. ilauplo( 'U' ) ) THEN
281  DO i = 1, n
282  IF ( beta .EQ. zero ) THEN
283  symb_zero = .true.
284  y( iy ) = 0.0d+0
285  ELSE IF ( y( iy ) .EQ. zero ) THEN
286  symb_zero = .true.
287  ELSE
288  symb_zero = .false.
289  y( iy ) = beta * abs( y( iy ) )
290  END IF
291  IF ( alpha .NE. zero ) THEN
292  DO j = 1, i
293  temp = cabs1( a( j, i ) )
294  symb_zero = symb_zero .AND.
295  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
296 
297  y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
298  END DO
299  DO j = i+1, n
300  temp = cabs1( a( i, j ) )
301  symb_zero = symb_zero .AND.
302  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
303 
304  y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
305  END DO
306  END IF
307 
308  IF ( .NOT.symb_zero )
309  $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
310 
311  iy = iy + incy
312  END DO
313  ELSE
314  DO i = 1, n
315  IF ( beta .EQ. zero ) THEN
316  symb_zero = .true.
317  y( iy ) = 0.0d+0
318  ELSE IF ( y( iy ) .EQ. zero ) THEN
319  symb_zero = .true.
320  ELSE
321  symb_zero = .false.
322  y( iy ) = beta * abs( y( iy ) )
323  END IF
324  IF ( alpha .NE. zero ) THEN
325  DO j = 1, i
326  temp = cabs1( a( i, j ) )
327  symb_zero = symb_zero .AND.
328  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
329 
330  y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
331  END DO
332  DO j = i+1, n
333  temp = cabs1( a( j, i ) )
334  symb_zero = symb_zero .AND.
335  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
336 
337  y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
338  END DO
339  END IF
340 
341  IF ( .NOT.symb_zero )
342  $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
343 
344  iy = iy + incy
345  END DO
346  END IF
347  ELSE
348  IF ( uplo .EQ. ilauplo( 'U' ) ) THEN
349  DO i = 1, n
350  IF ( beta .EQ. zero ) THEN
351  symb_zero = .true.
352  y( iy ) = 0.0d+0
353  ELSE IF ( y( iy ) .EQ. zero ) THEN
354  symb_zero = .true.
355  ELSE
356  symb_zero = .false.
357  y( iy ) = beta * abs( y( iy ) )
358  END IF
359  jx = kx
360  IF ( alpha .NE. zero ) THEN
361  DO j = 1, i
362  temp = cabs1( a( j, i ) )
363  symb_zero = symb_zero .AND.
364  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
365 
366  y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
367  jx = jx + incx
368  END DO
369  DO j = i+1, n
370  temp = cabs1( a( i, j ) )
371  symb_zero = symb_zero .AND.
372  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
373 
374  y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
375  jx = jx + incx
376  END DO
377  END IF
378 
379  IF ( .NOT.symb_zero )
380  $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
381 
382  iy = iy + incy
383  END DO
384  ELSE
385  DO i = 1, n
386  IF ( beta .EQ. zero ) THEN
387  symb_zero = .true.
388  y( iy ) = 0.0d+0
389  ELSE IF ( y( iy ) .EQ. zero ) THEN
390  symb_zero = .true.
391  ELSE
392  symb_zero = .false.
393  y( iy ) = beta * abs( y( iy ) )
394  END IF
395  jx = kx
396  IF ( alpha .NE. zero ) THEN
397  DO j = 1, i
398  temp = cabs1( a( i, j ) )
399  symb_zero = symb_zero .AND.
400  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
401 
402  y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
403  jx = jx + incx
404  END DO
405  DO j = i+1, n
406  temp = cabs1( a( j, i ) )
407  symb_zero = symb_zero .AND.
408  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
409 
410  y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
411  jx = jx + incx
412  END DO
413  END IF
414 
415  IF ( .NOT.symb_zero )
416  $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
417 
418  iy = iy + incy
419  END DO
420  END IF
421 
422  END IF
423 *
424  RETURN
425 *
426 * End of ZLA_SYAMV
427 *
integer function ilauplo(UPLO)
ILAUPLO
Definition: ilauplo.f:60
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65

Here is the call graph for this function:

Here is the caller graph for this function:

double precision function zla_syrcond_c ( character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldaf, * )  AF,
integer  LDAF,
integer, dimension( * )  IPIV,
double precision, dimension( * )  C,
logical  CAPPLY,
integer  INFO,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK 
)

ZLA_SYRCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for symmetric indefinite matrices.

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

Purpose:
    ZLA_SYRCOND_C Computes the infinity norm condition number of
    op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
       = 'U':  Upper triangle of A is stored;
       = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
     The number of linear equations, i.e., the order of the
     matrix A.  N >= 0.
[in]A
          A is COMPLEX*16 array, dimension (LDA,N)
     On entry, the N-by-N matrix A
[in]LDA
          LDA is INTEGER
     The leading dimension of the array A.  LDA >= max(1,N).
[in]AF
          AF is COMPLEX*16 array, dimension (LDAF,N)
     The block diagonal matrix D and the multipliers used to
     obtain the factor U or L as computed by ZSYTRF.
[in]LDAF
          LDAF is INTEGER
     The leading dimension of the array AF.  LDAF >= 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.
[in]C
          C is DOUBLE PRECISION array, dimension (N)
     The vector C in the formula op(A) * inv(diag(C)).
[in]CAPPLY
          CAPPLY is LOGICAL
     If .TRUE. then access the vector C in the formula above.
[out]INFO
          INFO is INTEGER
       = 0:  Successful exit.
     i > 0:  The ith argument is invalid.
[in]WORK
          WORK is COMPLEX*16 array, dimension (2*N).
     Workspace.
[in]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N).
     Workspace.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 142 of file zla_syrcond_c.f.

142 *
143 * -- LAPACK computational routine (version 3.4.2) --
144 * -- LAPACK is a software package provided by Univ. of Tennessee, --
145 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
146 * September 2012
147 *
148 * .. Scalar Arguments ..
149  CHARACTER uplo
150  LOGICAL capply
151  INTEGER n, lda, ldaf, info
152 * ..
153 * .. Array Arguments ..
154  INTEGER ipiv( * )
155  COMPLEX*16 a( lda, * ), af( ldaf, * ), work( * )
156  DOUBLE PRECISION c( * ), rwork( * )
157 * ..
158 *
159 * =====================================================================
160 *
161 * .. Local Scalars ..
162  INTEGER kase
163  DOUBLE PRECISION ainvnm, anorm, tmp
164  INTEGER i, j
165  LOGICAL up, upper
166  COMPLEX*16 zdum
167 * ..
168 * .. Local Arrays ..
169  INTEGER isave( 3 )
170 * ..
171 * .. External Functions ..
172  LOGICAL lsame
173  EXTERNAL lsame
174 * ..
175 * .. External Subroutines ..
176  EXTERNAL zlacn2, zsytrs, xerbla
177 * ..
178 * .. Intrinsic Functions ..
179  INTRINSIC abs, max
180 * ..
181 * .. Statement Functions ..
182  DOUBLE PRECISION cabs1
183 * ..
184 * .. Statement Function Definitions ..
185  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
186 * ..
187 * .. Executable Statements ..
188 *
189  zla_syrcond_c = 0.0d+0
190 *
191  info = 0
192  upper = lsame( uplo, 'U' )
193  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
194  info = -1
195  ELSE IF( n.LT.0 ) THEN
196  info = -2
197  ELSE IF( lda.LT.max( 1, n ) ) THEN
198  info = -4
199  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
200  info = -6
201  END IF
202  IF( info.NE.0 ) THEN
203  CALL xerbla( 'ZLA_SYRCOND_C', -info )
204  RETURN
205  END IF
206  up = .false.
207  IF ( lsame( uplo, 'U' ) ) up = .true.
208 *
209 * Compute norm of op(A)*op2(C).
210 *
211  anorm = 0.0d+0
212  IF ( up ) THEN
213  DO i = 1, n
214  tmp = 0.0d+0
215  IF ( capply ) THEN
216  DO j = 1, i
217  tmp = tmp + cabs1( a( j, i ) ) / c( j )
218  END DO
219  DO j = i+1, n
220  tmp = tmp + cabs1( a( i, j ) ) / c( j )
221  END DO
222  ELSE
223  DO j = 1, i
224  tmp = tmp + cabs1( a( j, i ) )
225  END DO
226  DO j = i+1, n
227  tmp = tmp + cabs1( a( i, j ) )
228  END DO
229  END IF
230  rwork( i ) = tmp
231  anorm = max( anorm, tmp )
232  END DO
233  ELSE
234  DO i = 1, n
235  tmp = 0.0d+0
236  IF ( capply ) THEN
237  DO j = 1, i
238  tmp = tmp + cabs1( a( i, j ) ) / c( j )
239  END DO
240  DO j = i+1, n
241  tmp = tmp + cabs1( a( j, i ) ) / c( j )
242  END DO
243  ELSE
244  DO j = 1, i
245  tmp = tmp + cabs1( a( i, j ) )
246  END DO
247  DO j = i+1, n
248  tmp = tmp + cabs1( a( j, i ) )
249  END DO
250  END IF
251  rwork( i ) = tmp
252  anorm = max( anorm, tmp )
253  END DO
254  END IF
255 *
256 * Quick return if possible.
257 *
258  IF( n.EQ.0 ) THEN
259  zla_syrcond_c = 1.0d+0
260  RETURN
261  ELSE IF( anorm .EQ. 0.0d+0 ) THEN
262  RETURN
263  END IF
264 *
265 * Estimate the norm of inv(op(A)).
266 *
267  ainvnm = 0.0d+0
268 *
269  kase = 0
270  10 CONTINUE
271  CALL zlacn2( n, work( n+1 ), work, ainvnm, kase, isave )
272  IF( kase.NE.0 ) THEN
273  IF( kase.EQ.2 ) THEN
274 *
275 * Multiply by R.
276 *
277  DO i = 1, n
278  work( i ) = work( i ) * rwork( i )
279  END DO
280 *
281  IF ( up ) THEN
282  CALL zsytrs( 'U', n, 1, af, ldaf, ipiv,
283  $ work, n, info )
284  ELSE
285  CALL zsytrs( 'L', n, 1, af, ldaf, ipiv,
286  $ work, n, info )
287  ENDIF
288 *
289 * Multiply by inv(C).
290 *
291  IF ( capply ) THEN
292  DO i = 1, n
293  work( i ) = work( i ) * c( i )
294  END DO
295  END IF
296  ELSE
297 *
298 * Multiply by inv(C**T).
299 *
300  IF ( capply ) THEN
301  DO i = 1, n
302  work( i ) = work( i ) * c( i )
303  END DO
304  END IF
305 *
306  IF ( up ) THEN
307  CALL zsytrs( 'U', n, 1, af, ldaf, ipiv,
308  $ work, n, info )
309  ELSE
310  CALL zsytrs( 'L', n, 1, af, ldaf, ipiv,
311  $ work, n, info )
312  END IF
313 *
314 * Multiply by R.
315 *
316  DO i = 1, n
317  work( i ) = work( i ) * rwork( i )
318  END DO
319  END IF
320  GO TO 10
321  END IF
322 *
323 * Compute the estimate of the reciprocal condition number.
324 *
325  IF( ainvnm .NE. 0.0d+0 )
326  $ zla_syrcond_c = 1.0d+0 / ainvnm
327 *
328  RETURN
329 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlacn2(N, V, X, EST, KASE, ISAVE)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: zlacn2.f:135
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zsytrs(UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZSYTRS
Definition: zsytrs.f:122
double precision function zla_syrcond_c(UPLO, N, A, LDA, AF, LDAF, IPIV, C, CAPPLY, INFO, WORK, RWORK)
ZLA_SYRCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for symmetric indefin...

Here is the call graph for this function:

Here is the caller graph for this function:

double precision function zla_syrcond_x ( character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldaf, * )  AF,
integer  LDAF,
integer, dimension( * )  IPIV,
complex*16, dimension( * )  X,
integer  INFO,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK 
)

ZLA_SYRCOND_X computes the infinity norm condition number of op(A)*diag(x) for symmetric indefinite matrices.

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

Purpose:
    ZLA_SYRCOND_X Computes the infinity norm condition number of
    op(A) * diag(X) where X is a COMPLEX*16 vector.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
       = 'U':  Upper triangle of A is stored;
       = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
     The number of linear equations, i.e., the order of the
     matrix A.  N >= 0.
[in]A
          A is COMPLEX*16 array, dimension (LDA,N)
     On entry, the N-by-N matrix A.
[in]LDA
          LDA is INTEGER
     The leading dimension of the array A.  LDA >= max(1,N).
[in]AF
          AF is COMPLEX*16 array, dimension (LDAF,N)
     The block diagonal matrix D and the multipliers used to
     obtain the factor U or L as computed by ZSYTRF.
[in]LDAF
          LDAF is INTEGER
     The leading dimension of the array AF.  LDAF >= 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.
[in]X
          X is COMPLEX*16 array, dimension (N)
     The vector X in the formula op(A) * diag(X).
[out]INFO
          INFO is INTEGER
       = 0:  Successful exit.
     i > 0:  The ith argument is invalid.
[in]WORK
          WORK is COMPLEX*16 array, dimension (2*N).
     Workspace.
[in]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N).
     Workspace.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 135 of file zla_syrcond_x.f.

135 *
136 * -- LAPACK computational routine (version 3.4.2) --
137 * -- LAPACK is a software package provided by Univ. of Tennessee, --
138 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
139 * September 2012
140 *
141 * .. Scalar Arguments ..
142  CHARACTER uplo
143  INTEGER n, lda, ldaf, info
144 * ..
145 * .. Array Arguments ..
146  INTEGER ipiv( * )
147  COMPLEX*16 a( lda, * ), af( ldaf, * ), work( * ), x( * )
148  DOUBLE PRECISION rwork( * )
149 * ..
150 *
151 * =====================================================================
152 *
153 * .. Local Scalars ..
154  INTEGER kase
155  DOUBLE PRECISION ainvnm, anorm, tmp
156  INTEGER i, j
157  LOGICAL up, upper
158  COMPLEX*16 zdum
159 * ..
160 * .. Local Arrays ..
161  INTEGER isave( 3 )
162 * ..
163 * .. External Functions ..
164  LOGICAL lsame
165  EXTERNAL lsame
166 * ..
167 * .. External Subroutines ..
168  EXTERNAL zlacn2, zsytrs, xerbla
169 * ..
170 * .. Intrinsic Functions ..
171  INTRINSIC abs, max
172 * ..
173 * .. Statement Functions ..
174  DOUBLE PRECISION cabs1
175 * ..
176 * .. Statement Function Definitions ..
177  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
178 * ..
179 * .. Executable Statements ..
180 *
181  zla_syrcond_x = 0.0d+0
182 *
183  info = 0
184  upper = lsame( uplo, 'U' )
185  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
186  info = -1
187  ELSE IF( n.LT.0 ) THEN
188  info = -2
189  ELSE IF( lda.LT.max( 1, n ) ) THEN
190  info = -4
191  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
192  info = -6
193  END IF
194  IF( info.NE.0 ) THEN
195  CALL xerbla( 'ZLA_SYRCOND_X', -info )
196  RETURN
197  END IF
198  up = .false.
199  IF ( lsame( uplo, 'U' ) ) up = .true.
200 *
201 * Compute norm of op(A)*op2(C).
202 *
203  anorm = 0.0d+0
204  IF ( up ) THEN
205  DO i = 1, n
206  tmp = 0.0d+0
207  DO j = 1, i
208  tmp = tmp + cabs1( a( j, i ) * x( j ) )
209  END DO
210  DO j = i+1, n
211  tmp = tmp + cabs1( a( i, j ) * x( j ) )
212  END DO
213  rwork( i ) = tmp
214  anorm = max( anorm, tmp )
215  END DO
216  ELSE
217  DO i = 1, n
218  tmp = 0.0d+0
219  DO j = 1, i
220  tmp = tmp + cabs1( a( i, j ) * x( j ) )
221  END DO
222  DO j = i+1, n
223  tmp = tmp + cabs1( a( j, i ) * x( j ) )
224  END DO
225  rwork( i ) = tmp
226  anorm = max( anorm, tmp )
227  END DO
228  END IF
229 *
230 * Quick return if possible.
231 *
232  IF( n.EQ.0 ) THEN
233  zla_syrcond_x = 1.0d+0
234  RETURN
235  ELSE IF( anorm .EQ. 0.0d+0 ) THEN
236  RETURN
237  END IF
238 *
239 * Estimate the norm of inv(op(A)).
240 *
241  ainvnm = 0.0d+0
242 *
243  kase = 0
244  10 CONTINUE
245  CALL zlacn2( n, work( n+1 ), work, ainvnm, kase, isave )
246  IF( kase.NE.0 ) THEN
247  IF( kase.EQ.2 ) THEN
248 *
249 * Multiply by R.
250 *
251  DO i = 1, n
252  work( i ) = work( i ) * rwork( i )
253  END DO
254 *
255  IF ( up ) THEN
256  CALL zsytrs( 'U', n, 1, af, ldaf, ipiv,
257  $ work, n, info )
258  ELSE
259  CALL zsytrs( 'L', n, 1, af, ldaf, ipiv,
260  $ work, n, info )
261  ENDIF
262 *
263 * Multiply by inv(X).
264 *
265  DO i = 1, n
266  work( i ) = work( i ) / x( i )
267  END DO
268  ELSE
269 *
270 * Multiply by inv(X**T).
271 *
272  DO i = 1, n
273  work( i ) = work( i ) / x( i )
274  END DO
275 *
276  IF ( up ) THEN
277  CALL zsytrs( 'U', n, 1, af, ldaf, ipiv,
278  $ work, n, info )
279  ELSE
280  CALL zsytrs( 'L', n, 1, af, ldaf, ipiv,
281  $ work, n, info )
282  END IF
283 *
284 * Multiply by R.
285 *
286  DO i = 1, n
287  work( i ) = work( i ) * rwork( i )
288  END DO
289  END IF
290  GO TO 10
291  END IF
292 *
293 * Compute the estimate of the reciprocal condition number.
294 *
295  IF( ainvnm .NE. 0.0d+0 )
296  $ zla_syrcond_x = 1.0d+0 / ainvnm
297 *
298  RETURN
299 *
double precision function zla_syrcond_x(UPLO, N, A, LDA, AF, LDAF, IPIV, X, INFO, WORK, RWORK)
ZLA_SYRCOND_X computes the infinity norm condition number of op(A)*diag(x) for symmetric indefinite m...
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlacn2(N, V, X, EST, KASE, ISAVE)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: zlacn2.f:135
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zsytrs(UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZSYTRS
Definition: zsytrs.f:122

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zla_syrfsx_extended ( integer  PREC_TYPE,
character  UPLO,
integer  N,
integer  NRHS,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldaf, * )  AF,
integer  LDAF,
integer, dimension( * )  IPIV,
logical  COLEQU,
double precision, dimension( * )  C,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( ldy, * )  Y,
integer  LDY,
double precision, dimension( * )  BERR_OUT,
integer  N_NORMS,
double precision, dimension( nrhs, * )  ERR_BNDS_NORM,
double precision, dimension( nrhs, * )  ERR_BNDS_COMP,
complex*16, dimension( * )  RES,
double precision, dimension( * )  AYB,
complex*16, dimension( * )  DY,
complex*16, dimension( * )  Y_TAIL,
double precision  RCOND,
integer  ITHRESH,
double precision  RTHRESH,
double precision  DZ_UB,
logical  IGNORE_CWISE,
integer  INFO 
)

ZLA_SYRFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric indefinite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.

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

Purpose:
 ZLA_SYRFSX_EXTENDED improves the computed solution to a system of
 linear equations by performing extra-precise iterative refinement
 and provides error bounds and backward error estimates for the solution.
 This subroutine is called by ZSYRFSX to perform iterative refinement.
 In addition to normwise error bound, the code provides maximum
 componentwise error bound if possible. See comments for ERR_BNDS_NORM
 and ERR_BNDS_COMP for details of the error bounds. Note that this
 subroutine is only resonsible for setting the second fields of
 ERR_BNDS_NORM and ERR_BNDS_COMP.
Parameters
[in]PREC_TYPE
          PREC_TYPE is INTEGER
     Specifies the intermediate precision to be used in refinement.
     The value is defined by ILAPREC(P) where P is a CHARACTER and
     P    = 'S':  Single
          = 'D':  Double
          = 'I':  Indigenous
          = 'X', 'E':  Extra
[in]UPLO
          UPLO is CHARACTER*1
       = 'U':  Upper triangle of A is stored;
       = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
     The number of linear equations, i.e., the order of the
     matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
     The number of right-hand-sides, i.e., the number of columns of the
     matrix B.
[in]A
          A is COMPLEX*16 array, dimension (LDA,N)
     On entry, the N-by-N matrix A.
[in]LDA
          LDA is INTEGER
     The leading dimension of the array A.  LDA >= max(1,N).
[in]AF
          AF is COMPLEX*16 array, dimension (LDAF,N)
     The block diagonal matrix D and the multipliers used to
     obtain the factor U or L as computed by ZSYTRF.
[in]LDAF
          LDAF is INTEGER
     The leading dimension of the array AF.  LDAF >= 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.
[in]COLEQU
          COLEQU is LOGICAL
     If .TRUE. then column equilibration was done to A before calling
     this routine. This is needed to compute the solution and error
     bounds correctly.
[in]C
          C is DOUBLE PRECISION array, dimension (N)
     The column scale factors for A. If COLEQU = .FALSE., C
     is not accessed. If C is input, each element of C should be a power
     of the radix to ensure a reliable solution and error estimates.
     Scaling by powers of the radix does not cause rounding errors unless
     the result underflows or overflows. Rounding errors during scaling
     lead to refining with a matrix that is not equivalent to the
     input matrix, producing error estimates that may not be
     reliable.
[in]B
          B is COMPLEX*16 array, dimension (LDB,NRHS)
     The right-hand-side matrix B.
[in]LDB
          LDB is INTEGER
     The leading dimension of the array B.  LDB >= max(1,N).
[in,out]Y
          Y is COMPLEX*16 array, dimension
                    (LDY,NRHS)
     On entry, the solution matrix X, as computed by ZSYTRS.
     On exit, the improved solution matrix Y.
[in]LDY
          LDY is INTEGER
     The leading dimension of the array Y.  LDY >= max(1,N).
[out]BERR_OUT
          BERR_OUT is DOUBLE PRECISION array, dimension (NRHS)
     On exit, BERR_OUT(j) contains the componentwise relative backward
     error for right-hand-side j from the formula
         max(i) ( abs(RES(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
     where abs(Z) is the componentwise absolute value of the matrix
     or vector Z. This is computed by ZLA_LIN_BERR.
[in]N_NORMS
          N_NORMS is INTEGER
     Determines which error bounds to return (see ERR_BNDS_NORM
     and ERR_BNDS_COMP).
     If N_NORMS >= 1 return normwise error bounds.
     If N_NORMS >= 2 return componentwise error bounds.
[in,out]ERR_BNDS_NORM
          ERR_BNDS_NORM is DOUBLE PRECISION array, dimension
                    (NRHS, N_ERR_BNDS)
     For each right-hand side, this array contains information about
     various error bounds and condition numbers corresponding to the
     normwise relative error, which is defined as follows:

     Normwise relative error in the ith solution vector:
             max_j (abs(XTRUE(j,i) - X(j,i)))
            ------------------------------
                  max_j abs(X(j,i))

     The array is indexed by the type of error information as described
     below. There currently are up to three pieces of information
     returned.

     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
     right-hand side.

     The second index in ERR_BNDS_NORM(:,err) contains the following
     three fields:
     err = 1 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * slamch('Epsilon').

     err = 2 "Guaranteed" error bound: The estimated forward error,
              almost certainly within a factor of 10 of the true error
              so long as the next entry is greater than the threshold
              sqrt(n) * slamch('Epsilon'). This error bound should only
              be trusted if the previous boolean is true.

     err = 3  Reciprocal condition number: Estimated normwise
              reciprocal condition number.  Compared with the threshold
              sqrt(n) * slamch('Epsilon') to determine if the error
              estimate is "guaranteed". These reciprocal condition
              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
              appropriately scaled matrix Z.
              Let Z = S*A, where S scales each row by a power of the
              radix so all absolute row sums of Z are approximately 1.

     This subroutine is only responsible for setting the second field
     above.
     See Lapack Working Note 165 for further details and extra
     cautions.
[in,out]ERR_BNDS_COMP
          ERR_BNDS_COMP is DOUBLE PRECISION array, dimension
                    (NRHS, N_ERR_BNDS)
     For each right-hand side, this array contains information about
     various error bounds and condition numbers corresponding to the
     componentwise relative error, which is defined as follows:

     Componentwise relative error in the ith solution vector:
                    abs(XTRUE(j,i) - X(j,i))
             max_j ----------------------
                         abs(X(j,i))

     The array is indexed by the right-hand side i (on which the
     componentwise relative error depends), and the type of error
     information as described below. There currently are up to three
     pieces of information returned for each right-hand side. If
     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS .LT. 3, then at most
     the first (:,N_ERR_BNDS) entries are returned.

     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
     right-hand side.

     The second index in ERR_BNDS_COMP(:,err) contains the following
     three fields:
     err = 1 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * slamch('Epsilon').

     err = 2 "Guaranteed" error bound: The estimated forward error,
              almost certainly within a factor of 10 of the true error
              so long as the next entry is greater than the threshold
              sqrt(n) * slamch('Epsilon'). This error bound should only
              be trusted if the previous boolean is true.

     err = 3  Reciprocal condition number: Estimated componentwise
              reciprocal condition number.  Compared with the threshold
              sqrt(n) * slamch('Epsilon') to determine if the error
              estimate is "guaranteed". These reciprocal condition
              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
              appropriately scaled matrix Z.
              Let Z = S*(A*diag(x)), where x is the solution for the
              current right-hand side and S scales each row of
              A*diag(x) by a power of the radix so all absolute row
              sums of Z are approximately 1.

     This subroutine is only responsible for setting the second field
     above.
     See Lapack Working Note 165 for further details and extra
     cautions.
[in]RES
          RES is COMPLEX*16 array, dimension (N)
     Workspace to hold the intermediate residual.
[in]AYB
          AYB is DOUBLE PRECISION array, dimension (N)
     Workspace.
[in]DY
          DY is COMPLEX*16 array, dimension (N)
     Workspace to hold the intermediate solution.
[in]Y_TAIL
          Y_TAIL is COMPLEX*16 array, dimension (N)
     Workspace to hold the trailing bits of the intermediate solution.
[in]RCOND
          RCOND is DOUBLE PRECISION
     Reciprocal scaled condition number.  This is an estimate of the
     reciprocal Skeel condition number of the matrix A after
     equilibration (if done).  If this is less than the machine
     precision (in particular, if it is zero), the matrix is singular
     to working precision.  Note that the error may still be small even
     if this number is very small and the matrix appears ill-
     conditioned.
[in]ITHRESH
          ITHRESH is INTEGER
     The maximum number of residual computations allowed for
     refinement. The default is 10. For 'aggressive' set to 100 to
     permit convergence using approximate factorizations or
     factorizations other than LU. If the factorization uses a
     technique other than Gaussian elimination, the guarantees in
     ERR_BNDS_NORM and ERR_BNDS_COMP may no longer be trustworthy.
[in]RTHRESH
          RTHRESH is DOUBLE PRECISION
     Determines when to stop refinement if the error estimate stops
     decreasing. Refinement will stop when the next solution no longer
     satisfies norm(dx_{i+1}) < RTHRESH * norm(dx_i) where norm(Z) is
     the infinity norm of Z. RTHRESH satisfies 0 < RTHRESH <= 1. The
     default value is 0.5. For 'aggressive' set to 0.9 to permit
     convergence on extremely ill-conditioned matrices. See LAWN 165
     for more details.
[in]DZ_UB
          DZ_UB is DOUBLE PRECISION
     Determines when to start considering componentwise convergence.
     Componentwise convergence is only considered after each component
     of the solution Y is stable, which we definte as the relative
     change in each component being less than DZ_UB. The default value
     is 0.25, requiring the first bit to be stable. See LAWN 165 for
     more details.
[in]IGNORE_CWISE
          IGNORE_CWISE is LOGICAL
     If .TRUE. then ignore componentwise convergence. Default value
     is .FALSE..
[out]INFO
          INFO is INTEGER
       = 0:  Successful exit.
       < 0:  if INFO = -i, the ith argument to ZLA_HERFSX_EXTENDED had an illegal
             value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 400 of file zla_syrfsx_extended.f.

400 *
401 * -- LAPACK computational routine (version 3.4.2) --
402 * -- LAPACK is a software package provided by Univ. of Tennessee, --
403 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
404 * September 2012
405 *
406 * .. Scalar Arguments ..
407  INTEGER info, lda, ldaf, ldb, ldy, n, nrhs, prec_type,
408  $ n_norms, ithresh
409  CHARACTER uplo
410  LOGICAL colequ, ignore_cwise
411  DOUBLE PRECISION rthresh, dz_ub
412 * ..
413 * .. Array Arguments ..
414  INTEGER ipiv( * )
415  COMPLEX*16 a( lda, * ), af( ldaf, * ), b( ldb, * ),
416  $ y( ldy, * ), res( * ), dy( * ), y_tail( * )
417  DOUBLE PRECISION c( * ), ayb( * ), rcond, berr_out( * ),
418  $ err_bnds_norm( nrhs, * ),
419  $ err_bnds_comp( nrhs, * )
420 * ..
421 *
422 * =====================================================================
423 *
424 * .. Local Scalars ..
425  INTEGER uplo2, cnt, i, j, x_state, z_state,
426  $ y_prec_state
427  DOUBLE PRECISION yk, dyk, ymin, normy, normx, normdx, dxrat,
428  $ dzrat, prevnormdx, prev_dz_z, dxratmax,
429  $ dzratmax, dx_x, dz_z, final_dx_x, final_dz_z,
430  $ eps, hugeval, incr_thresh
431  LOGICAL incr_prec, upper
432  COMPLEX*16 zdum
433 * ..
434 * .. Parameters ..
435  INTEGER unstable_state, working_state, conv_state,
436  $ noprog_state, base_residual, extra_residual,
437  $ extra_y
438  parameter( unstable_state = 0, working_state = 1,
439  $ conv_state = 2, noprog_state = 3 )
440  parameter( base_residual = 0, extra_residual = 1,
441  $ extra_y = 2 )
442  INTEGER final_nrm_err_i, final_cmp_err_i, berr_i
443  INTEGER rcond_i, nrm_rcond_i, nrm_err_i, cmp_rcond_i
444  INTEGER cmp_err_i, piv_growth_i
445  parameter( final_nrm_err_i = 1, final_cmp_err_i = 2,
446  $ berr_i = 3 )
447  parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
448  parameter( cmp_rcond_i = 7, cmp_err_i = 8,
449  $ piv_growth_i = 9 )
450  INTEGER la_linrx_itref_i, la_linrx_ithresh_i,
451  $ la_linrx_cwise_i
452  parameter( la_linrx_itref_i = 1,
453  $ la_linrx_ithresh_i = 2 )
454  parameter( la_linrx_cwise_i = 3 )
455  INTEGER la_linrx_trust_i, la_linrx_err_i,
456  $ la_linrx_rcond_i
457  parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
458  parameter( la_linrx_rcond_i = 3 )
459 * ..
460 * .. External Functions ..
461  LOGICAL lsame
462  EXTERNAL ilauplo
463  INTEGER ilauplo
464 * ..
465 * .. External Subroutines ..
466  EXTERNAL zaxpy, zcopy, zsytrs, zsymv, blas_zsymv_x,
467  $ blas_zsymv2_x, zla_syamv, zla_wwaddw,
468  $ zla_lin_berr
469  DOUBLE PRECISION dlamch
470 * ..
471 * .. Intrinsic Functions ..
472  INTRINSIC abs, REAL, dimag, max, min
473 * ..
474 * .. Statement Functions ..
475  DOUBLE PRECISION cabs1
476 * ..
477 * .. Statement Function Definitions ..
478  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
479 * ..
480 * .. Executable Statements ..
481 *
482  info = 0
483  upper = lsame( uplo, 'U' )
484  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
485  info = -2
486  ELSE IF( n.LT.0 ) THEN
487  info = -3
488  ELSE IF( nrhs.LT.0 ) THEN
489  info = -4
490  ELSE IF( lda.LT.max( 1, n ) ) THEN
491  info = -6
492  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
493  info = -8
494  ELSE IF( ldb.LT.max( 1, n ) ) THEN
495  info = -13
496  ELSE IF( ldy.LT.max( 1, n ) ) THEN
497  info = -15
498  END IF
499  IF( info.NE.0 ) THEN
500  CALL xerbla( 'ZLA_HERFSX_EXTENDED', -info )
501  RETURN
502  END IF
503  eps = dlamch( 'Epsilon' )
504  hugeval = dlamch( 'Overflow' )
505 * Force HUGEVAL to Inf
506  hugeval = hugeval * hugeval
507 * Using HUGEVAL may lead to spurious underflows.
508  incr_thresh = dble( n ) * eps
509 
510  IF ( lsame( uplo, 'L' ) ) THEN
511  uplo2 = ilauplo( 'L' )
512  ELSE
513  uplo2 = ilauplo( 'U' )
514  ENDIF
515 
516  DO j = 1, nrhs
517  y_prec_state = extra_residual
518  IF ( y_prec_state .EQ. extra_y ) THEN
519  DO i = 1, n
520  y_tail( i ) = 0.0d+0
521  END DO
522  END IF
523 
524  dxrat = 0.0d+0
525  dxratmax = 0.0d+0
526  dzrat = 0.0d+0
527  dzratmax = 0.0d+0
528  final_dx_x = hugeval
529  final_dz_z = hugeval
530  prevnormdx = hugeval
531  prev_dz_z = hugeval
532  dz_z = hugeval
533  dx_x = hugeval
534 
535  x_state = working_state
536  z_state = unstable_state
537  incr_prec = .false.
538 
539  DO cnt = 1, ithresh
540 *
541 * Compute residual RES = B_s - op(A_s) * Y,
542 * op(A) = A, A**T, or A**H depending on TRANS (and type).
543 *
544  CALL zcopy( n, b( 1, j ), 1, res, 1 )
545  IF ( y_prec_state .EQ. base_residual ) THEN
546  CALL zsymv( uplo, n, dcmplx(-1.0d+0), a, lda, y(1,j), 1,
547  $ dcmplx(1.0d+0), res, 1 )
548  ELSE IF ( y_prec_state .EQ. extra_residual ) THEN
549  CALL blas_zsymv_x( uplo2, n, dcmplx(-1.0d+0), a, lda,
550  $ y( 1, j ), 1, dcmplx(1.0d+0), res, 1, prec_type )
551  ELSE
552  CALL blas_zsymv2_x(uplo2, n, dcmplx(-1.0d+0), a, lda,
553  $ y(1, j), y_tail, 1, dcmplx(1.0d+0), res, 1,
554  $ prec_type)
555  END IF
556 
557 ! XXX: RES is no longer needed.
558  CALL zcopy( n, res, 1, dy, 1 )
559  CALL zsytrs( uplo, n, 1, af, ldaf, ipiv, dy, n, info )
560 *
561 * Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
562 *
563  normx = 0.0d+0
564  normy = 0.0d+0
565  normdx = 0.0d+0
566  dz_z = 0.0d+0
567  ymin = hugeval
568 
569  DO i = 1, n
570  yk = cabs1( y( i, j ) )
571  dyk = cabs1( dy( i ) )
572 
573  IF ( yk .NE. 0.0d+0 ) THEN
574  dz_z = max( dz_z, dyk / yk )
575  ELSE IF ( dyk .NE. 0.0d+0 ) THEN
576  dz_z = hugeval
577  END IF
578 
579  ymin = min( ymin, yk )
580 
581  normy = max( normy, yk )
582 
583  IF ( colequ ) THEN
584  normx = max( normx, yk * c( i ) )
585  normdx = max( normdx, dyk * c( i ) )
586  ELSE
587  normx = normy
588  normdx = max( normdx, dyk )
589  END IF
590  END DO
591 
592  IF ( normx .NE. 0.0d+0 ) THEN
593  dx_x = normdx / normx
594  ELSE IF ( normdx .EQ. 0.0d+0 ) THEN
595  dx_x = 0.0d+0
596  ELSE
597  dx_x = hugeval
598  END IF
599 
600  dxrat = normdx / prevnormdx
601  dzrat = dz_z / prev_dz_z
602 *
603 * Check termination criteria.
604 *
605  IF ( ymin*rcond .LT. incr_thresh*normy
606  $ .AND. y_prec_state .LT. extra_y )
607  $ incr_prec = .true.
608 
609  IF ( x_state .EQ. noprog_state .AND. dxrat .LE. rthresh )
610  $ x_state = working_state
611  IF ( x_state .EQ. working_state ) THEN
612  IF ( dx_x .LE. eps ) THEN
613  x_state = conv_state
614  ELSE IF ( dxrat .GT. rthresh ) THEN
615  IF ( y_prec_state .NE. extra_y ) THEN
616  incr_prec = .true.
617  ELSE
618  x_state = noprog_state
619  END IF
620  ELSE
621  IF (dxrat .GT. dxratmax) dxratmax = dxrat
622  END IF
623  IF ( x_state .GT. working_state ) final_dx_x = dx_x
624  END IF
625 
626  IF ( z_state .EQ. unstable_state .AND. dz_z .LE. dz_ub )
627  $ z_state = working_state
628  IF ( z_state .EQ. noprog_state .AND. dzrat .LE. rthresh )
629  $ z_state = working_state
630  IF ( z_state .EQ. working_state ) THEN
631  IF ( dz_z .LE. eps ) THEN
632  z_state = conv_state
633  ELSE IF ( dz_z .GT. dz_ub ) THEN
634  z_state = unstable_state
635  dzratmax = 0.0d+0
636  final_dz_z = hugeval
637  ELSE IF ( dzrat .GT. rthresh ) THEN
638  IF ( y_prec_state .NE. extra_y ) THEN
639  incr_prec = .true.
640  ELSE
641  z_state = noprog_state
642  END IF
643  ELSE
644  IF ( dzrat .GT. dzratmax ) dzratmax = dzrat
645  END IF
646  IF ( z_state .GT. working_state ) final_dz_z = dz_z
647  END IF
648 
649  IF ( x_state.NE.working_state.AND.
650  $ ( ignore_cwise.OR.z_state.NE.working_state ) )
651  $ GOTO 666
652 
653  IF ( incr_prec ) THEN
654  incr_prec = .false.
655  y_prec_state = y_prec_state + 1
656  DO i = 1, n
657  y_tail( i ) = 0.0d+0
658  END DO
659  END IF
660 
661  prevnormdx = normdx
662  prev_dz_z = dz_z
663 *
664 * Update soluton.
665 *
666  IF ( y_prec_state .LT. extra_y ) THEN
667  CALL zaxpy( n, dcmplx(1.0d+0), dy, 1, y(1,j), 1 )
668  ELSE
669  CALL zla_wwaddw( n, y(1,j), y_tail, dy )
670  END IF
671 
672  END DO
673 * Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT.
674  666 CONTINUE
675 *
676 * Set final_* when cnt hits ithresh.
677 *
678  IF ( x_state .EQ. working_state ) final_dx_x = dx_x
679  IF ( z_state .EQ. working_state ) final_dz_z = dz_z
680 *
681 * Compute error bounds.
682 *
683  IF ( n_norms .GE. 1 ) THEN
684  err_bnds_norm( j, la_linrx_err_i ) =
685  $ final_dx_x / (1 - dxratmax)
686  END IF
687  IF ( n_norms .GE. 2 ) THEN
688  err_bnds_comp( j, la_linrx_err_i ) =
689  $ final_dz_z / (1 - dzratmax)
690  END IF
691 *
692 * Compute componentwise relative backward error from formula
693 * max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
694 * where abs(Z) is the componentwise absolute value of the matrix
695 * or vector Z.
696 *
697 * Compute residual RES = B_s - op(A_s) * Y,
698 * op(A) = A, A**T, or A**H depending on TRANS (and type).
699 *
700  CALL zcopy( n, b( 1, j ), 1, res, 1 )
701  CALL zsymv( uplo, n, dcmplx(-1.0d+0), a, lda, y(1,j), 1,
702  $ dcmplx(1.0d+0), res, 1 )
703 
704  DO i = 1, n
705  ayb( i ) = cabs1( b( i, j ) )
706  END DO
707 *
708 * Compute abs(op(A_s))*abs(Y) + abs(B_s).
709 *
710  CALL zla_syamv ( uplo2, n, 1.0d+0,
711  $ a, lda, y(1, j), 1, 1.0d+0, ayb, 1 )
712 
713  CALL zla_lin_berr ( n, n, 1, res, ayb, berr_out( j ) )
714 *
715 * End of loop for each RHS.
716 *
717  END DO
718 *
719  RETURN
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
integer function ilauplo(UPLO)
ILAUPLO
Definition: ilauplo.f:60
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zla_lin_berr(N, NZ, NRHS, RES, AYB, BERR)
ZLA_LIN_BERR computes a component-wise relative backward error.
Definition: zla_lin_berr.f:103
subroutine zla_syamv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZLA_SYAMV computes a matrix-vector product using a symmetric indefinite matrix to calculate error bou...
Definition: zla_syamv.f:181
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
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zla_wwaddw(N, X, Y, W)
ZLA_WWADDW adds a vector into a doubled-single vector.
Definition: zla_wwaddw.f:83
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:53
subroutine zsytrs(UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZSYTRS
Definition: zsytrs.f:122

Here is the call graph for this function:

Here is the caller graph for this function:

double precision function zla_syrpvgrw ( character*1  UPLO,
integer  N,
integer  INFO,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldaf, * )  AF,
integer  LDAF,
integer, dimension( * )  IPIV,
double precision, dimension( * )  WORK 
)

ZLA_SYRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a symmetric indefinite matrix.

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

Purpose:
 ZLA_SYRPVGRW computes the reciprocal pivot growth factor
 norm(A)/norm(U). The "max absolute element" norm is used. If this is
 much less than 1, the stability of the LU factorization of the
 (equilibrated) matrix A could be poor. This also means that the
 solution X, estimated condition numbers, and error bounds could be
 unreliable.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
       = 'U':  Upper triangle of A is stored;
       = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
     The number of linear equations, i.e., the order of the
     matrix A.  N >= 0.
[in]INFO
          INFO is INTEGER
     The value of INFO returned from ZSYTRF, .i.e., the pivot in
     column INFO is exactly 0.
[in]A
          A is COMPLEX*16 array, dimension (LDA,N)
     On entry, the N-by-N matrix A.
[in]LDA
          LDA is INTEGER
     The leading dimension of the array A.  LDA >= max(1,N).
[in]AF
          AF is COMPLEX*16 array, dimension (LDAF,N)
     The block diagonal matrix D and the multipliers used to
     obtain the factor U or L as computed by ZSYTRF.
[in]LDAF
          LDAF is INTEGER
     The leading dimension of the array AF.  LDAF >= 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.
[in]WORK
          WORK is DOUBLE PRECISION array, dimension (2*N)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015

Definition at line 125 of file zla_syrpvgrw.f.

125 *
126 * -- LAPACK computational routine (version 3.6.0) --
127 * -- LAPACK is a software package provided by Univ. of Tennessee, --
128 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
129 * November 2015
130 *
131 * .. Scalar Arguments ..
132  CHARACTER*1 uplo
133  INTEGER n, info, lda, ldaf
134 * ..
135 * .. Array Arguments ..
136  COMPLEX*16 a( lda, * ), af( ldaf, * )
137  DOUBLE PRECISION work( * )
138  INTEGER ipiv( * )
139 * ..
140 *
141 * =====================================================================
142 *
143 * .. Local Scalars ..
144  INTEGER ncols, i, j, k, kp
145  DOUBLE PRECISION amax, umax, rpvgrw, tmp
146  LOGICAL upper
147  COMPLEX*16 zdum
148 * ..
149 * .. Intrinsic Functions ..
150  INTRINSIC abs, REAL, dimag, max, min
151 * ..
152 * .. External Subroutines ..
153  EXTERNAL lsame, zlaset
154  LOGICAL lsame
155 * ..
156 * .. Statement Functions ..
157  DOUBLE PRECISION cabs1
158 * ..
159 * .. Statement Function Definitions ..
160  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
161 * ..
162 * .. Executable Statements ..
163 *
164  upper = lsame( 'Upper', uplo )
165  IF ( info.EQ.0 ) THEN
166  IF ( upper ) THEN
167  ncols = 1
168  ELSE
169  ncols = n
170  END IF
171  ELSE
172  ncols = info
173  END IF
174 
175  rpvgrw = 1.0d+0
176  DO i = 1, 2*n
177  work( i ) = 0.0d+0
178  END DO
179 *
180 * Find the max magnitude entry of each column of A. Compute the max
181 * for all N columns so we can apply the pivot permutation while
182 * looping below. Assume a full factorization is the common case.
183 *
184  IF ( upper ) THEN
185  DO j = 1, n
186  DO i = 1, j
187  work( n+i ) = max( cabs1( a( i, j ) ), work( n+i ) )
188  work( n+j ) = max( cabs1( a( i, j ) ), work( n+j ) )
189  END DO
190  END DO
191  ELSE
192  DO j = 1, n
193  DO i = j, n
194  work( n+i ) = max( cabs1( a( i, j ) ), work( n+i ) )
195  work( n+j ) = max( cabs1( a( i, j ) ), work( n+j ) )
196  END DO
197  END DO
198  END IF
199 *
200 * Now find the max magnitude entry of each column of U or L. Also
201 * permute the magnitudes of A above so they're in the same order as
202 * the factor.
203 *
204 * The iteration orders and permutations were copied from zsytrs.
205 * Calls to SSWAP would be severe overkill.
206 *
207  IF ( upper ) THEN
208  k = n
209  DO WHILE ( k .LT. ncols .AND. k.GT.0 )
210  IF ( ipiv( k ).GT.0 ) THEN
211 ! 1x1 pivot
212  kp = ipiv( k )
213  IF ( kp .NE. k ) THEN
214  tmp = work( n+k )
215  work( n+k ) = work( n+kp )
216  work( n+kp ) = tmp
217  END IF
218  DO i = 1, k
219  work( k ) = max( cabs1( af( i, k ) ), work( k ) )
220  END DO
221  k = k - 1
222  ELSE
223 ! 2x2 pivot
224  kp = -ipiv( k )
225  tmp = work( n+k-1 )
226  work( n+k-1 ) = work( n+kp )
227  work( n+kp ) = tmp
228  DO i = 1, k-1
229  work( k ) = max( cabs1( af( i, k ) ), work( k ) )
230  work( k-1 ) =
231  $ max( cabs1( af( i, k-1 ) ), work( k-1 ) )
232  END DO
233  work( k ) = max( cabs1( af( k, k ) ), work( k ) )
234  k = k - 2
235  END IF
236  END DO
237  k = ncols
238  DO WHILE ( k .LE. n )
239  IF ( ipiv( k ).GT.0 ) THEN
240  kp = ipiv( k )
241  IF ( kp .NE. k ) THEN
242  tmp = work( n+k )
243  work( n+k ) = work( n+kp )
244  work( n+kp ) = tmp
245  END IF
246  k = k + 1
247  ELSE
248  kp = -ipiv( k )
249  tmp = work( n+k )
250  work( n+k ) = work( n+kp )
251  work( n+kp ) = tmp
252  k = k + 2
253  END IF
254  END DO
255  ELSE
256  k = 1
257  DO WHILE ( k .LE. ncols )
258  IF ( ipiv( k ).GT.0 ) THEN
259 ! 1x1 pivot
260  kp = ipiv( k )
261  IF ( kp .NE. k ) THEN
262  tmp = work( n+k )
263  work( n+k ) = work( n+kp )
264  work( n+kp ) = tmp
265  END IF
266  DO i = k, n
267  work( k ) = max( cabs1( af( i, k ) ), work( k ) )
268  END DO
269  k = k + 1
270  ELSE
271 ! 2x2 pivot
272  kp = -ipiv( k )
273  tmp = work( n+k+1 )
274  work( n+k+1 ) = work( n+kp )
275  work( n+kp ) = tmp
276  DO i = k+1, n
277  work( k ) = max( cabs1( af( i, k ) ), work( k ) )
278  work( k+1 ) =
279  $ max( cabs1( af( i, k+1 ) ), work( k+1 ) )
280  END DO
281  work( k ) = max( cabs1( af( k, k ) ), work( k ) )
282  k = k + 2
283  END IF
284  END DO
285  k = ncols
286  DO WHILE ( k .GE. 1 )
287  IF ( ipiv( k ).GT.0 ) THEN
288  kp = ipiv( k )
289  IF ( kp .NE. k ) THEN
290  tmp = work( n+k )
291  work( n+k ) = work( n+kp )
292  work( n+kp ) = tmp
293  END IF
294  k = k - 1
295  ELSE
296  kp = -ipiv( k )
297  tmp = work( n+k )
298  work( n+k ) = work( n+kp )
299  work( n+kp ) = tmp
300  k = k - 2
301  ENDIF
302  END DO
303  END IF
304 *
305 * Compute the *inverse* of the max element growth factor. Dividing
306 * by zero would imply the largest entry of the factor's column is
307 * zero. Than can happen when either the column of A is zero or
308 * massive pivots made the factor underflow to zero. Neither counts
309 * as growth in itself, so simply ignore terms with zero
310 * denominators.
311 *
312  IF ( upper ) THEN
313  DO i = ncols, n
314  umax = work( i )
315  amax = work( n+i )
316  IF ( umax /= 0.0d+0 ) THEN
317  rpvgrw = min( amax / umax, rpvgrw )
318  END IF
319  END DO
320  ELSE
321  DO i = 1, ncols
322  umax = work( i )
323  amax = work( n+i )
324  IF ( umax /= 0.0d+0 ) THEN
325  rpvgrw = min( amax / umax, rpvgrw )
326  END IF
327  END DO
328  END IF
329 
330  zla_syrpvgrw = rpvgrw
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
double precision function zla_syrpvgrw(UPLO, N, INFO, A, LDA, AF, LDAF, IPIV, WORK)
ZLA_SYRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a symmetric indefinite m...
Definition: zla_syrpvgrw.f:125

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zlasyf ( character  UPLO,
integer  N,
integer  NB,
integer  KB,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex*16, dimension( ldw, * )  W,
integer  LDW,
integer  INFO 
)

ZLASYF computes a partial factorization of a complex symmetric matrix using the Bunch-Kaufman diagonal pivoting method.

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

Purpose:
 ZLASYF computes a partial factorization of a complex symmetric matrix
 A using the Bunch-Kaufman diagonal pivoting method. The partial
 factorization has the form:

 A  =  ( I  U12 ) ( A11  0  ) (  I       0    )  if UPLO = 'U', or:
       ( 0  U22 ) (  0   D  ) ( U12**T U22**T )

 A  =  ( L11  0 ) ( D    0  ) ( L11**T L21**T )  if UPLO = 'L'
       ( L21  I ) ( 0   A22 ) (  0       I    )

 where the order of D is at most NB. The actual order is returned in
 the argument KB, and is either NB or NB-1, or N if N <= NB.
 Note that U**T denotes the transpose of U.

 ZLASYF is an auxiliary routine called by ZSYTRF. It uses blocked code
 (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
 A22 (if UPLO = 'L').
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the upper or lower triangular part of the
          symmetric matrix A is stored:
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]NB
          NB is INTEGER
          The maximum number of columns of the matrix A that should be
          factored.  NB should be at least 2 to allow for 2-by-2 pivot
          blocks.
[out]KB
          KB is INTEGER
          The number of columns of A that were actually factored.
          KB is either NB-1 or NB, or N if N <= NB.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
          n-by-n upper triangular part of A contains the upper
          triangular part of the matrix A, and the strictly lower
          triangular part of A is not referenced.  If UPLO = 'L', the
          leading n-by-n lower triangular part of A contains the lower
          triangular part of the matrix A, and the strictly upper
          triangular part of A is not referenced.
          On exit, A contains details of the partial factorization.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D.

          If UPLO = 'U':
             Only the last KB elements of IPIV are set.

             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
             interchanged and D(k,k) is a 1-by-1 diagonal block.

             If IPIV(k) = IPIV(k-1) < 0, then rows and columns
             k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
             is a 2-by-2 diagonal block.

          If UPLO = 'L':
             Only the first KB elements of IPIV are set.

             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
             interchanged and D(k,k) is a 1-by-1 diagonal block.

             If IPIV(k) = IPIV(k+1) < 0, then rows and columns
             k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1)
             is a 2-by-2 diagonal block.
[out]W
          W is COMPLEX*16 array, dimension (LDW,NB)
[in]LDW
          LDW is INTEGER
          The leading dimension of the array W.  LDW >= max(1,N).
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
               has been completed, but the block diagonal matrix D is
               exactly singular.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2013
Contributors:
  November 2013,  Igor Kozachenko,
                  Computer Science Division,
                  University of California, Berkeley

Definition at line 179 of file zlasyf.f.

179 *
180 * -- LAPACK computational routine (version 3.5.0) --
181 * -- LAPACK is a software package provided by Univ. of Tennessee, --
182 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
183 * November 2013
184 *
185 * .. Scalar Arguments ..
186  CHARACTER uplo
187  INTEGER info, kb, lda, ldw, n, nb
188 * ..
189 * .. Array Arguments ..
190  INTEGER ipiv( * )
191  COMPLEX*16 a( lda, * ), w( ldw, * )
192 * ..
193 *
194 * =====================================================================
195 *
196 * .. Parameters ..
197  DOUBLE PRECISION zero, one
198  parameter( zero = 0.0d+0, one = 1.0d+0 )
199  DOUBLE PRECISION eight, sevten
200  parameter( eight = 8.0d+0, sevten = 17.0d+0 )
201  COMPLEX*16 cone
202  parameter( cone = ( 1.0d+0, 0.0d+0 ) )
203 * ..
204 * .. Local Scalars ..
205  INTEGER imax, j, jb, jj, jmax, jp, k, kk, kkw, kp,
206  $ kstep, kw
207  DOUBLE PRECISION absakk, alpha, colmax, rowmax
208  COMPLEX*16 d11, d21, d22, r1, t, z
209 * ..
210 * .. External Functions ..
211  LOGICAL lsame
212  INTEGER izamax
213  EXTERNAL lsame, izamax
214 * ..
215 * .. External Subroutines ..
216  EXTERNAL zcopy, zgemm, zgemv, zscal, zswap
217 * ..
218 * .. Intrinsic Functions ..
219  INTRINSIC abs, dble, dimag, max, min, sqrt
220 * ..
221 * .. Statement Functions ..
222  DOUBLE PRECISION cabs1
223 * ..
224 * .. Statement Function definitions ..
225  cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
226 * ..
227 * .. Executable Statements ..
228 *
229  info = 0
230 *
231 * Initialize ALPHA for use in choosing pivot block size.
232 *
233  alpha = ( one+sqrt( sevten ) ) / eight
234 *
235  IF( lsame( uplo, 'U' ) ) THEN
236 *
237 * Factorize the trailing columns of A using the upper triangle
238 * of A and working backwards, and compute the matrix W = U12*D
239 * for use in updating A11
240 *
241 * K is the main loop index, decreasing from N in steps of 1 or 2
242 *
243 * KW is the column of W which corresponds to column K of A
244 *
245  k = n
246  10 CONTINUE
247  kw = nb + k - n
248 *
249 * Exit from loop
250 *
251  IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
252  $ GO TO 30
253 *
254 * Copy column K of A to column KW of W and update it
255 *
256  CALL zcopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
257  IF( k.LT.n )
258  $ CALL zgemv( 'No transpose', k, n-k, -cone, a( 1, k+1 ), lda,
259  $ w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
260 *
261  kstep = 1
262 *
263 * Determine rows and columns to be interchanged and whether
264 * a 1-by-1 or 2-by-2 pivot block will be used
265 *
266  absakk = cabs1( w( k, kw ) )
267 *
268 * IMAX is the row-index of the largest off-diagonal element in
269 
270 *
271  IF( k.GT.1 ) THEN
272  imax = izamax( k-1, w( 1, kw ), 1 )
273  colmax = cabs1( w( imax, kw ) )
274  ELSE
275  colmax = zero
276  END IF
277 *
278  IF( max( absakk, colmax ).EQ.zero ) THEN
279 *
280 * Column K is zero or underflow: set INFO and continue
281 *
282  IF( info.EQ.0 )
283  $ info = k
284  kp = k
285  ELSE
286  IF( absakk.GE.alpha*colmax ) THEN
287 *
288 * no interchange, use 1-by-1 pivot block
289 *
290  kp = k
291  ELSE
292 *
293 * Copy column IMAX to column KW-1 of W and update it
294 *
295  CALL zcopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
296  CALL zcopy( k-imax, a( imax, imax+1 ), lda,
297  $ w( imax+1, kw-1 ), 1 )
298  IF( k.LT.n )
299  $ CALL zgemv( 'No transpose', k, n-k, -cone,
300  $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
301  $ cone, w( 1, kw-1 ), 1 )
302 *
303 * JMAX is the column-index of the largest off-diagonal
304 * element in row IMAX, and ROWMAX is its absolute value
305 *
306  jmax = imax + izamax( k-imax, w( imax+1, kw-1 ), 1 )
307  rowmax = cabs1( w( jmax, kw-1 ) )
308  IF( imax.GT.1 ) THEN
309  jmax = izamax( imax-1, w( 1, kw-1 ), 1 )
310  rowmax = max( rowmax, cabs1( w( jmax, kw-1 ) ) )
311  END IF
312 *
313  IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
314 *
315 * no interchange, use 1-by-1 pivot block
316 *
317  kp = k
318  ELSE IF( cabs1( w( imax, kw-1 ) ).GE.alpha*rowmax ) THEN
319 *
320 * interchange rows and columns K and IMAX, use 1-by-1
321 * pivot block
322 *
323  kp = imax
324 *
325 * copy column KW-1 of W to column KW of W
326 *
327  CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
328  ELSE
329 *
330 * interchange rows and columns K-1 and IMAX, use 2-by-2
331 * pivot block
332 *
333  kp = imax
334  kstep = 2
335  END IF
336  END IF
337 *
338 * ============================================================
339 *
340 * KK is the column of A where pivoting step stopped
341 *
342  kk = k - kstep + 1
343 *
344 * KKW is the column of W which corresponds to column KK of A
345 *
346  kkw = nb + kk - n
347 *
348 * Interchange rows and columns KP and KK.
349 * Updated column KP is already stored in column KKW of W.
350 *
351  IF( kp.NE.kk ) THEN
352 *
353 * Copy non-updated column KK to column KP of submatrix A
354 * at step K. No need to copy element into column K
355 * (or K and K-1 for 2-by-2 pivot) of A, since these columns
356 * will be later overwritten.
357 *
358  a( kp, kp ) = a( kk, kk )
359  CALL zcopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
360  $ lda )
361  IF( kp.GT.1 )
362  $ CALL zcopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
363 *
364 * Interchange rows KK and KP in last K+1 to N columns of A
365 * (columns K (or K and K-1 for 2-by-2 pivot) of A will be
366 * later overwritten). Interchange rows KK and KP
367 * in last KKW to NB columns of W.
368 *
369  IF( k.LT.n )
370  $ CALL zswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
371  $ lda )
372  CALL zswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
373  $ ldw )
374  END IF
375 *
376  IF( kstep.EQ.1 ) THEN
377 *
378 * 1-by-1 pivot block D(k): column kw of W now holds
379 *
380 * W(kw) = U(k)*D(k),
381 *
382 * where U(k) is the k-th column of U
383 *
384 * Store subdiag. elements of column U(k)
385 * and 1-by-1 block D(k) in column k of A.
386 * NOTE: Diagonal element U(k,k) is a UNIT element
387 * and not stored.
388 * A(k,k) := D(k,k) = W(k,kw)
389 * A(1:k-1,k) := U(1:k-1,k) = W(1:k-1,kw)/D(k,k)
390 *
391  CALL zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
392  r1 = cone / a( k, k )
393  CALL zscal( k-1, r1, a( 1, k ), 1 )
394 *
395  ELSE
396 *
397 * 2-by-2 pivot block D(k): columns kw and kw-1 of W now hold
398 *
399 * ( W(kw-1) W(kw) ) = ( U(k-1) U(k) )*D(k)
400 *
401 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
402 * of U
403 *
404 * Store U(1:k-2,k-1) and U(1:k-2,k) and 2-by-2
405 * block D(k-1:k,k-1:k) in columns k-1 and k of A.
406 * NOTE: 2-by-2 diagonal block U(k-1:k,k-1:k) is a UNIT
407 * block and not stored.
408 * A(k-1:k,k-1:k) := D(k-1:k,k-1:k) = W(k-1:k,kw-1:kw)
409 * A(1:k-2,k-1:k) := U(1:k-2,k:k-1:k) =
410 * = W(1:k-2,kw-1:kw) * ( D(k-1:k,k-1:k)**(-1) )
411 *
412  IF( k.GT.2 ) THEN
413 *
414 * Compose the columns of the inverse of 2-by-2 pivot
415 * block D in the following way to reduce the number
416 * of FLOPS when we myltiply panel ( W(kw-1) W(kw) ) by
417 * this inverse
418 *
419 * D**(-1) = ( d11 d21 )**(-1) =
420 * ( d21 d22 )
421 *
422 * = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
423 * ( (-d21 ) ( d11 ) )
424 *
425 * = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
426 *
427 * * ( ( d22/d21 ) ( -1 ) ) =
428 * ( ( -1 ) ( d11/d21 ) )
429 *
430 * = 1/d21 * 1/(D22*D11-1) * ( ( D11 ) ( -1 ) ) =
431 * ( ( -1 ) ( D22 ) )
432 *
433 * = 1/d21 * T * ( ( D11 ) ( -1 ) )
434 * ( ( -1 ) ( D22 ) )
435 *
436 * = D21 * ( ( D11 ) ( -1 ) )
437 * ( ( -1 ) ( D22 ) )
438 *
439  d21 = w( k-1, kw )
440  d11 = w( k, kw ) / d21
441  d22 = w( k-1, kw-1 ) / d21
442  t = cone / ( d11*d22-cone )
443  d21 = t / d21
444 *
445 * Update elements in columns A(k-1) and A(k) as
446 * dot products of rows of ( W(kw-1) W(kw) ) and columns
447 * of D**(-1)
448 *
449  DO 20 j = 1, k - 2
450  a( j, k-1 ) = d21*( d11*w( j, kw-1 )-w( j, kw ) )
451  a( j, k ) = d21*( d22*w( j, kw )-w( j, kw-1 ) )
452  20 CONTINUE
453  END IF
454 *
455 * Copy D(k) to A
456 *
457  a( k-1, k-1 ) = w( k-1, kw-1 )
458  a( k-1, k ) = w( k-1, kw )
459  a( k, k ) = w( k, kw )
460 *
461  END IF
462 *
463  END IF
464 *
465 * Store details of the interchanges in IPIV
466 *
467  IF( kstep.EQ.1 ) THEN
468  ipiv( k ) = kp
469  ELSE
470  ipiv( k ) = -kp
471  ipiv( k-1 ) = -kp
472  END IF
473 *
474 * Decrease K and return to the start of the main loop
475 *
476  k = k - kstep
477  GO TO 10
478 *
479  30 CONTINUE
480 *
481 * Update the upper triangle of A11 (= A(1:k,1:k)) as
482 *
483 * A11 := A11 - U12*D*U12**T = A11 - U12*W**T
484 *
485 * computing blocks of NB columns at a time
486 *
487  DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
488  jb = min( nb, k-j+1 )
489 *
490 * Update the upper triangle of the diagonal block
491 *
492  DO 40 jj = j, j + jb - 1
493  CALL zgemv( 'No transpose', jj-j+1, n-k, -cone,
494  $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
495  $ a( j, jj ), 1 )
496  40 CONTINUE
497 *
498 * Update the rectangular superdiagonal block
499 *
500  CALL zgemm( 'No transpose', 'Transpose', j-1, jb, n-k,
501  $ -cone, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
502  $ cone, a( 1, j ), lda )
503  50 CONTINUE
504 *
505 * Put U12 in standard form by partially undoing the interchanges
506 * in columns k+1:n looping backwards from k+1 to n
507 *
508  j = k + 1
509  60 CONTINUE
510 *
511 * Undo the interchanges (if any) of rows JJ and JP at each
512 * step J
513 *
514 * (Here, J is a diagonal index)
515  jj = j
516  jp = ipiv( j )
517  IF( jp.LT.0 ) THEN
518  jp = -jp
519 * (Here, J is a diagonal index)
520  j = j + 1
521  END IF
522 * (NOTE: Here, J is used to determine row length. Length N-J+1
523 * of the rows to swap back doesn't include diagonal element)
524  j = j + 1
525  IF( jp.NE.jj .AND. j.LE.n )
526  $ CALL zswap( n-j+1, a( jp, j ), lda, a( jj, j ), lda )
527  IF( j.LT.n )
528  $ GO TO 60
529 *
530 * Set KB to the number of columns factorized
531 *
532  kb = n - k
533 *
534  ELSE
535 *
536 * Factorize the leading columns of A using the lower triangle
537 * of A and working forwards, and compute the matrix W = L21*D
538 * for use in updating A22
539 *
540 * K is the main loop index, increasing from 1 in steps of 1 or 2
541 *
542  k = 1
543  70 CONTINUE
544 *
545 * Exit from loop
546 *
547  IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
548  $ GO TO 90
549 *
550 * Copy column K of A to column K of W and update it
551 *
552  CALL zcopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
553  CALL zgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ), lda,
554  $ w( k, 1 ), ldw, cone, w( k, k ), 1 )
555 *
556  kstep = 1
557 *
558 * Determine rows and columns to be interchanged and whether
559 * a 1-by-1 or 2-by-2 pivot block will be used
560 *
561  absakk = cabs1( w( k, k ) )
562 *
563 * IMAX is the row-index of the largest off-diagonal element in
564 
565 *
566  IF( k.LT.n ) THEN
567  imax = k + izamax( n-k, w( k+1, k ), 1 )
568  colmax = cabs1( w( imax, k ) )
569  ELSE
570  colmax = zero
571  END IF
572 *
573  IF( max( absakk, colmax ).EQ.zero ) THEN
574 *
575 * Column K is zero or underflow: set INFO and continue
576 *
577  IF( info.EQ.0 )
578  $ info = k
579  kp = k
580  ELSE
581  IF( absakk.GE.alpha*colmax ) THEN
582 *
583 * no interchange, use 1-by-1 pivot block
584 *
585  kp = k
586  ELSE
587 *
588 * Copy column IMAX to column K+1 of W and update it
589 *
590  CALL zcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1 )
591  CALL zcopy( n-imax+1, a( imax, imax ), 1, w( imax, k+1 ),
592  $ 1 )
593  CALL zgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
594  $ lda, w( imax, 1 ), ldw, cone, w( k, k+1 ),
595  $ 1 )
596 *
597 * JMAX is the column-index of the largest off-diagonal
598 * element in row IMAX, and ROWMAX is its absolute value
599 *
600  jmax = k - 1 + izamax( imax-k, w( k, k+1 ), 1 )
601  rowmax = cabs1( w( jmax, k+1 ) )
602  IF( imax.LT.n ) THEN
603  jmax = imax + izamax( n-imax, w( imax+1, k+1 ), 1 )
604  rowmax = max( rowmax, cabs1( w( jmax, k+1 ) ) )
605  END IF
606 *
607  IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
608 *
609 * no interchange, use 1-by-1 pivot block
610 *
611  kp = k
612  ELSE IF( cabs1( w( imax, k+1 ) ).GE.alpha*rowmax ) THEN
613 *
614 * interchange rows and columns K and IMAX, use 1-by-1
615 * pivot block
616 *
617  kp = imax
618 *
619 * copy column K+1 of W to column K of W
620 *
621  CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
622  ELSE
623 *
624 * interchange rows and columns K+1 and IMAX, use 2-by-2
625 * pivot block
626 *
627  kp = imax
628  kstep = 2
629  END IF
630  END IF
631 *
632 * ============================================================
633 *
634 * KK is the column of A where pivoting step stopped
635 *
636  kk = k + kstep - 1
637 *
638 * Interchange rows and columns KP and KK.
639 * Updated column KP is already stored in column KK of W.
640 *
641  IF( kp.NE.kk ) THEN
642 *
643 * Copy non-updated column KK to column KP of submatrix A
644 * at step K. No need to copy element into column K
645 * (or K and K+1 for 2-by-2 pivot) of A, since these columns
646 * will be later overwritten.
647 *
648  a( kp, kp ) = a( kk, kk )
649  CALL zcopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
650  $ lda )
651  IF( kp.LT.n )
652  $ CALL zcopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
653 *
654 * Interchange rows KK and KP in first K-1 columns of A
655 * (columns K (or K and K+1 for 2-by-2 pivot) of A will be
656 * later overwritten). Interchange rows KK and KP
657 * in first KK columns of W.
658 *
659  IF( k.GT.1 )
660  $ CALL zswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
661  CALL zswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
662  END IF
663 *
664  IF( kstep.EQ.1 ) THEN
665 *
666 * 1-by-1 pivot block D(k): column k of W now holds
667 *
668 * W(k) = L(k)*D(k),
669 *
670 * where L(k) is the k-th column of L
671 *
672 * Store subdiag. elements of column L(k)
673 * and 1-by-1 block D(k) in column k of A.
674 * (NOTE: Diagonal element L(k,k) is a UNIT element
675 * and not stored)
676 * A(k,k) := D(k,k) = W(k,k)
677 * A(k+1:N,k) := L(k+1:N,k) = W(k+1:N,k)/D(k,k)
678 *
679  CALL zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
680  IF( k.LT.n ) THEN
681  r1 = cone / a( k, k )
682  CALL zscal( n-k, r1, a( k+1, k ), 1 )
683  END IF
684 *
685  ELSE
686 *
687 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
688 *
689 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
690 *
691 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
692 * of L
693 *
694 * Store L(k+2:N,k) and L(k+2:N,k+1) and 2-by-2
695 * block D(k:k+1,k:k+1) in columns k and k+1 of A.
696 * (NOTE: 2-by-2 diagonal block L(k:k+1,k:k+1) is a UNIT
697 * block and not stored)
698 * A(k:k+1,k:k+1) := D(k:k+1,k:k+1) = W(k:k+1,k:k+1)
699 * A(k+2:N,k:k+1) := L(k+2:N,k:k+1) =
700 * = W(k+2:N,k:k+1) * ( D(k:k+1,k:k+1)**(-1) )
701 *
702  IF( k.LT.n-1 ) THEN
703 *
704 * Compose the columns of the inverse of 2-by-2 pivot
705 * block D in the following way to reduce the number
706 * of FLOPS when we myltiply panel ( W(k) W(k+1) ) by
707 * this inverse
708 *
709 * D**(-1) = ( d11 d21 )**(-1) =
710 * ( d21 d22 )
711 *
712 * = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
713 * ( (-d21 ) ( d11 ) )
714 *
715 * = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
716 *
717 * * ( ( d22/d21 ) ( -1 ) ) =
718 * ( ( -1 ) ( d11/d21 ) )
719 *
720 * = 1/d21 * 1/(D22*D11-1) * ( ( D11 ) ( -1 ) ) =
721 * ( ( -1 ) ( D22 ) )
722 *
723 * = 1/d21 * T * ( ( D11 ) ( -1 ) )
724 * ( ( -1 ) ( D22 ) )
725 *
726 * = D21 * ( ( D11 ) ( -1 ) )
727 * ( ( -1 ) ( D22 ) )
728 *
729  d21 = w( k+1, k )
730  d11 = w( k+1, k+1 ) / d21
731  d22 = w( k, k ) / d21
732  t = cone / ( d11*d22-cone )
733  d21 = t / d21
734 *
735 * Update elements in columns A(k) and A(k+1) as
736 * dot products of rows of ( W(k) W(k+1) ) and columns
737 * of D**(-1)
738 *
739  DO 80 j = k + 2, n
740  a( j, k ) = d21*( d11*w( j, k )-w( j, k+1 ) )
741  a( j, k+1 ) = d21*( d22*w( j, k+1 )-w( j, k ) )
742  80 CONTINUE
743  END IF
744 *
745 * Copy D(k) to A
746 *
747  a( k, k ) = w( k, k )
748  a( k+1, k ) = w( k+1, k )
749  a( k+1, k+1 ) = w( k+1, k+1 )
750 *
751  END IF
752 *
753  END IF
754 *
755 * Store details of the interchanges in IPIV
756 *
757  IF( kstep.EQ.1 ) THEN
758  ipiv( k ) = kp
759  ELSE
760  ipiv( k ) = -kp
761  ipiv( k+1 ) = -kp
762  END IF
763 *
764 * Increase K and return to the start of the main loop
765 *
766  k = k + kstep
767  GO TO 70
768 *
769  90 CONTINUE
770 *
771 * Update the lower triangle of A22 (= A(k:n,k:n)) as
772 *
773 * A22 := A22 - L21*D*L21**T = A22 - L21*W**T
774 *
775 * computing blocks of NB columns at a time
776 *
777  DO 110 j = k, n, nb
778  jb = min( nb, n-j+1 )
779 *
780 * Update the lower triangle of the diagonal block
781 *
782  DO 100 jj = j, j + jb - 1
783  CALL zgemv( 'No transpose', j+jb-jj, k-1, -cone,
784  $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
785  $ a( jj, jj ), 1 )
786  100 CONTINUE
787 *
788 * Update the rectangular subdiagonal block
789 *
790  IF( j+jb.LE.n )
791  $ CALL zgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
792  $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
793  $ ldw, cone, a( j+jb, j ), lda )
794  110 CONTINUE
795 *
796 * Put L21 in standard form by partially undoing the interchanges
797 * of rows in columns 1:k-1 looping backwards from k-1 to 1
798 *
799  j = k - 1
800  120 CONTINUE
801 *
802 * Undo the interchanges (if any) of rows JJ and JP at each
803 * step J
804 *
805 * (Here, J is a diagonal index)
806  jj = j
807  jp = ipiv( j )
808  IF( jp.LT.0 ) THEN
809  jp = -jp
810 * (Here, J is a diagonal index)
811  j = j - 1
812  END IF
813 * (NOTE: Here, J is used to determine row length. Length J
814 * of the rows to swap back doesn't include diagonal element)
815  j = j - 1
816  IF( jp.NE.jj .AND. j.GE.1 )
817  $ CALL zswap( j, a( jp, 1 ), lda, a( jj, 1 ), lda )
818  IF( j.GT.1 )
819  $ GO TO 120
820 *
821 * Set KB to the number of columns factorized
822 *
823  kb = k - 1
824 *
825  END IF
826  RETURN
827 *
828 * End of ZLASYF
829 *
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:54
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zlasyf_rook ( character  UPLO,
integer  N,
integer  NB,
integer  KB,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex*16, dimension( ldw, * )  W,
integer  LDW,
integer  INFO 
)

ZLASYF_ROOK computes a partial factorization of a complex symmetric matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method.

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

Purpose:
 ZLASYF_ROOK computes a partial factorization of a complex symmetric
 matrix A using the bounded Bunch-Kaufman ("rook") diagonal
 pivoting method. The partial factorization has the form:

 A  =  ( I  U12 ) ( A11  0  ) (  I       0    )  if UPLO = 'U', or:
       ( 0  U22 ) (  0   D  ) ( U12**T U22**T )

 A  =  ( L11  0 ) (  D   0  ) ( L11**T L21**T )  if UPLO = 'L'
       ( L21  I ) (  0  A22 ) (  0       I    )

 where the order of D is at most NB. The actual order is returned in
 the argument KB, and is either NB or NB-1, or N if N <= NB.

 ZLASYF_ROOK is an auxiliary routine called by ZSYTRF_ROOK. It uses
 blocked code (calling Level 3 BLAS) to update the submatrix
 A11 (if UPLO = 'U') or A22 (if UPLO = 'L').
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the upper or lower triangular part of the
          symmetric matrix A is stored:
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]NB
          NB is INTEGER
          The maximum number of columns of the matrix A that should be
          factored.  NB should be at least 2 to allow for 2-by-2 pivot
          blocks.
[out]KB
          KB is INTEGER
          The number of columns of A that were actually factored.
          KB is either NB-1 or NB, or N if N <= NB.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
          n-by-n upper triangular part of A contains the upper
          triangular part of the matrix A, and the strictly lower
          triangular part of A is not referenced.  If UPLO = 'L', the
          leading n-by-n lower triangular part of A contains the lower
          triangular part of the matrix A, and the strictly upper
          triangular part of A is not referenced.
          On exit, A contains details of the partial factorization.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D.

          If UPLO = 'U':
             Only the last KB elements of IPIV are set.

             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
             interchanged and D(k,k) is a 1-by-1 diagonal block.

             If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
             columns k and -IPIV(k) were interchanged and rows and
             columns k-1 and -IPIV(k-1) were inerchaged,
             D(k-1:k,k-1:k) is a 2-by-2 diagonal block.

          If UPLO = 'L':
             Only the first KB elements of IPIV are set.

             If IPIV(k) > 0, then rows and columns k and IPIV(k)
             were interchanged and D(k,k) is a 1-by-1 diagonal block.

             If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
             columns k and -IPIV(k) were interchanged and rows and
             columns k+1 and -IPIV(k+1) were inerchaged,
             D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
[out]W
          W is COMPLEX*16 array, dimension (LDW,NB)
[in]LDW
          LDW is INTEGER
          The leading dimension of the array W.  LDW >= max(1,N).
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
               has been completed, but the block diagonal matrix D is
               exactly singular.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2013
Contributors:
  November 2013,  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 186 of file zlasyf_rook.f.

186 *
187 * -- LAPACK computational routine (version 3.5.0) --
188 * -- LAPACK is a software package provided by Univ. of Tennessee, --
189 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
190 * November 2013
191 *
192 * .. Scalar Arguments ..
193  CHARACTER uplo
194  INTEGER info, kb, lda, ldw, n, nb
195 * ..
196 * .. Array Arguments ..
197  INTEGER ipiv( * )
198  COMPLEX*16 a( lda, * ), w( ldw, * )
199 * ..
200 *
201 * =====================================================================
202 *
203 * .. Parameters ..
204  DOUBLE PRECISION zero, one
205  parameter( zero = 0.0d+0, one = 1.0d+0 )
206  DOUBLE PRECISION eight, sevten
207  parameter( eight = 8.0d+0, sevten = 17.0d+0 )
208  COMPLEX*16 cone, czero
209  parameter( cone = ( 1.0d+0, 0.0d+0 ),
210  $ czero = ( 0.0d+0, 0.0d+0 ) )
211 * ..
212 * .. Local Scalars ..
213  LOGICAL done
214  INTEGER imax, itemp, j, jb, jj, jmax, jp1, jp2, k, kk,
215  $ kw, kkw, kp, kstep, p, ii
216  DOUBLE PRECISION absakk, alpha, colmax, rowmax, dtemp, sfmin
217  COMPLEX*16 d11, d12, d21, d22, r1, t, z
218 * ..
219 * .. External Functions ..
220  LOGICAL lsame
221  INTEGER izamax
222  DOUBLE PRECISION dlamch
223  EXTERNAL lsame, izamax, dlamch
224 * ..
225 * .. External Subroutines ..
226  EXTERNAL zcopy, zgemm, zgemv, zscal, zswap
227 * ..
228 * .. Intrinsic Functions ..
229  INTRINSIC abs, max, min, sqrt, dimag, dble
230 * ..
231 * .. Statement Functions ..
232  DOUBLE PRECISION cabs1
233 * ..
234 * .. Statement Function definitions ..
235  cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
236 * ..
237 * .. Executable Statements ..
238 *
239  info = 0
240 *
241 * Initialize ALPHA for use in choosing pivot block size.
242 *
243  alpha = ( one+sqrt( sevten ) ) / eight
244 *
245 * Compute machine safe minimum
246 *
247  sfmin = dlamch( 'S' )
248 *
249  IF( lsame( uplo, 'U' ) ) THEN
250 *
251 * Factorize the trailing columns of A using the upper triangle
252 * of A and working backwards, and compute the matrix W = U12*D
253 * for use in updating A11
254 *
255 * K is the main loop index, decreasing from N in steps of 1 or 2
256 *
257  k = n
258  10 CONTINUE
259 *
260 * KW is the column of W which corresponds to column K of A
261 *
262  kw = nb + k - n
263 *
264 * Exit from loop
265 *
266  IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
267  $ GO TO 30
268 *
269  kstep = 1
270  p = k
271 *
272 * Copy column K of A to column KW of W and update it
273 *
274  CALL zcopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
275  IF( k.LT.n )
276  $ CALL zgemv( 'No transpose', k, n-k, -cone, a( 1, k+1 ),
277  $ lda, w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
278 *
279 * Determine rows and columns to be interchanged and whether
280 * a 1-by-1 or 2-by-2 pivot block will be used
281 *
282  absakk = cabs1( w( k, kw ) )
283 *
284 * IMAX is the row-index of the largest off-diagonal element in
285 * column K, and COLMAX is its absolute value.
286 * Determine both COLMAX and IMAX.
287 *
288  IF( k.GT.1 ) THEN
289  imax = izamax( k-1, w( 1, kw ), 1 )
290  colmax = cabs1( w( imax, kw ) )
291  ELSE
292  colmax = zero
293  END IF
294 *
295  IF( max( absakk, colmax ).EQ.zero ) THEN
296 *
297 * Column K is zero or underflow: set INFO and continue
298 *
299  IF( info.EQ.0 )
300  $ info = k
301  kp = k
302  CALL zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
303  ELSE
304 *
305 * ============================================================
306 *
307 * Test for interchange
308 *
309 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
310 * (used to handle NaN and Inf)
311 *
312  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
313 *
314 * no interchange, use 1-by-1 pivot block
315 *
316  kp = k
317 *
318  ELSE
319 *
320  done = .false.
321 *
322 * Loop until pivot found
323 *
324  12 CONTINUE
325 *
326 * Begin pivot search loop body
327 *
328 *
329 * Copy column IMAX to column KW-1 of W and update it
330 *
331  CALL zcopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
332  CALL zcopy( k-imax, a( imax, imax+1 ), lda,
333  $ w( imax+1, kw-1 ), 1 )
334 *
335  IF( k.LT.n )
336  $ CALL zgemv( 'No transpose', k, n-k, -cone,
337  $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
338  $ cone, w( 1, kw-1 ), 1 )
339 *
340 * JMAX is the column-index of the largest off-diagonal
341 * element in row IMAX, and ROWMAX is its absolute value.
342 * Determine both ROWMAX and JMAX.
343 *
344  IF( imax.NE.k ) THEN
345  jmax = imax + izamax( k-imax, w( imax+1, kw-1 ),
346  $ 1 )
347  rowmax = cabs1( w( jmax, kw-1 ) )
348  ELSE
349  rowmax = zero
350  END IF
351 *
352  IF( imax.GT.1 ) THEN
353  itemp = izamax( imax-1, w( 1, kw-1 ), 1 )
354  dtemp = cabs1( w( itemp, kw-1 ) )
355  IF( dtemp.GT.rowmax ) THEN
356  rowmax = dtemp
357  jmax = itemp
358  END IF
359  END IF
360 *
361 * Equivalent to testing for
362 * CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX
363 * (used to handle NaN and Inf)
364 *
365  IF( .NOT.(cabs1( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
366  $ THEN
367 *
368 * interchange rows and columns K and IMAX,
369 * use 1-by-1 pivot block
370 *
371  kp = imax
372 *
373 * copy column KW-1 of W to column KW of W
374 *
375  CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
376 *
377  done = .true.
378 *
379 * Equivalent to testing for ROWMAX.EQ.COLMAX,
380 * (used to handle NaN and Inf)
381 *
382  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
383  $ THEN
384 *
385 * interchange rows and columns K-1 and IMAX,
386 * use 2-by-2 pivot block
387 *
388  kp = imax
389  kstep = 2
390  done = .true.
391  ELSE
392 *
393 * Pivot not found: set params and repeat
394 *
395  p = imax
396  colmax = rowmax
397  imax = jmax
398 *
399 * Copy updated JMAXth (next IMAXth) column to Kth of W
400 *
401  CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
402 *
403  END IF
404 *
405 * End pivot search loop body
406 *
407  IF( .NOT. done ) GOTO 12
408 *
409  END IF
410 *
411 * ============================================================
412 *
413  kk = k - kstep + 1
414 *
415 * KKW is the column of W which corresponds to column KK of A
416 *
417  kkw = nb + kk - n
418 *
419  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
420 *
421 * Copy non-updated column K to column P
422 *
423  CALL zcopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
424  CALL zcopy( p, a( 1, k ), 1, a( 1, p ), 1 )
425 *
426 * Interchange rows K and P in last N-K+1 columns of A
427 * and last N-K+2 columns of W
428 *
429  CALL zswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
430  CALL zswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
431  END IF
432 *
433 * Updated column KP is already stored in column KKW of W
434 *
435  IF( kp.NE.kk ) THEN
436 *
437 * Copy non-updated column KK to column KP
438 *
439  a( kp, k ) = a( kk, k )
440  CALL zcopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
441  $ lda )
442  CALL zcopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
443 *
444 * Interchange rows KK and KP in last N-KK+1 columns
445 * of A and W
446 *
447  CALL zswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
448  CALL zswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
449  $ ldw )
450  END IF
451 *
452  IF( kstep.EQ.1 ) THEN
453 *
454 * 1-by-1 pivot block D(k): column KW of W now holds
455 *
456 * W(k) = U(k)*D(k)
457 *
458 * where U(k) is the k-th column of U
459 *
460 * Store U(k) in column k of A
461 *
462  CALL zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
463  IF( k.GT.1 ) THEN
464  IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
465  r1 = cone / a( k, k )
466  CALL zscal( k-1, r1, a( 1, k ), 1 )
467  ELSE IF( a( k, k ).NE.czero ) THEN
468  DO 14 ii = 1, k - 1
469  a( ii, k ) = a( ii, k ) / a( k, k )
470  14 CONTINUE
471  END IF
472  END IF
473 *
474  ELSE
475 *
476 * 2-by-2 pivot block D(k): columns KW and KW-1 of W now
477 * hold
478 *
479 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
480 *
481 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
482 * of U
483 *
484  IF( k.GT.2 ) THEN
485 *
486 * Store U(k) and U(k-1) in columns k and k-1 of A
487 *
488  d12 = w( k-1, kw )
489  d11 = w( k, kw ) / d12
490  d22 = w( k-1, kw-1 ) / d12
491  t = cone / ( d11*d22-cone )
492  DO 20 j = 1, k - 2
493  a( j, k-1 ) = t*( (d11*w( j, kw-1 )-w( j, kw ) ) /
494  $ d12 )
495  a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
496  $ d12 )
497  20 CONTINUE
498  END IF
499 *
500 * Copy D(k) to A
501 *
502  a( k-1, k-1 ) = w( k-1, kw-1 )
503  a( k-1, k ) = w( k-1, kw )
504  a( k, k ) = w( k, kw )
505  END IF
506  END IF
507 *
508 * Store details of the interchanges in IPIV
509 *
510  IF( kstep.EQ.1 ) THEN
511  ipiv( k ) = kp
512  ELSE
513  ipiv( k ) = -p
514  ipiv( k-1 ) = -kp
515  END IF
516 *
517 * Decrease K and return to the start of the main loop
518 *
519  k = k - kstep
520  GO TO 10
521 *
522  30 CONTINUE
523 *
524 * Update the upper triangle of A11 (= A(1:k,1:k)) as
525 *
526 * A11 := A11 - U12*D*U12**T = A11 - U12*W**T
527 *
528 * computing blocks of NB columns at a time
529 *
530  DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
531  jb = min( nb, k-j+1 )
532 *
533 * Update the upper triangle of the diagonal block
534 *
535  DO 40 jj = j, j + jb - 1
536  CALL zgemv( 'No transpose', jj-j+1, n-k, -cone,
537  $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
538  $ a( j, jj ), 1 )
539  40 CONTINUE
540 *
541 * Update the rectangular superdiagonal block
542 *
543  IF( j.GE.2 )
544  $ CALL zgemm( 'No transpose', 'Transpose', j-1, jb,
545  $ n-k, -cone, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
546  $ cone, a( 1, j ), lda )
547  50 CONTINUE
548 *
549 * Put U12 in standard form by partially undoing the interchanges
550 * in columns k+1:n
551 *
552  j = k + 1
553  60 CONTINUE
554 *
555  kstep = 1
556  jp1 = 1
557  jj = j
558  jp2 = ipiv( j )
559  IF( jp2.LT.0 ) THEN
560  jp2 = -jp2
561  j = j + 1
562  jp1 = -ipiv( j )
563  kstep = 2
564  END IF
565 *
566  j = j + 1
567  IF( jp2.NE.jj .AND. j.LE.n )
568  $ CALL zswap( n-j+1, a( jp2, j ), lda, a( jj, j ), lda )
569  jj = j - 1
570  IF( jp1.NE.jj .AND. kstep.EQ.2 )
571  $ CALL zswap( n-j+1, a( jp1, j ), lda, a( jj, j ), lda )
572  IF( j.LE.n )
573  $ GO TO 60
574 *
575 * Set KB to the number of columns factorized
576 *
577  kb = n - k
578 *
579  ELSE
580 *
581 * Factorize the leading columns of A using the lower triangle
582 * of A and working forwards, and compute the matrix W = L21*D
583 * for use in updating A22
584 *
585 * K is the main loop index, increasing from 1 in steps of 1 or 2
586 *
587  k = 1
588  70 CONTINUE
589 *
590 * Exit from loop
591 *
592  IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
593  $ GO TO 90
594 *
595  kstep = 1
596  p = k
597 *
598 * Copy column K of A to column K of W and update it
599 *
600  CALL zcopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
601  IF( k.GT.1 )
602  $ CALL zgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
603  $ lda, w( k, 1 ), ldw, cone, w( k, k ), 1 )
604 *
605 * Determine rows and columns to be interchanged and whether
606 * a 1-by-1 or 2-by-2 pivot block will be used
607 *
608  absakk = cabs1( w( k, k ) )
609 *
610 * IMAX is the row-index of the largest off-diagonal element in
611 * column K, and COLMAX is its absolute value.
612 * Determine both COLMAX and IMAX.
613 *
614  IF( k.LT.n ) THEN
615  imax = k + izamax( n-k, w( k+1, k ), 1 )
616  colmax = cabs1( w( imax, k ) )
617  ELSE
618  colmax = zero
619  END IF
620 *
621  IF( max( absakk, colmax ).EQ.zero ) THEN
622 *
623 * Column K is zero or underflow: set INFO and continue
624 *
625  IF( info.EQ.0 )
626  $ info = k
627  kp = k
628  CALL zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
629  ELSE
630 *
631 * ============================================================
632 *
633 * Test for interchange
634 *
635 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
636 * (used to handle NaN and Inf)
637 *
638  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
639 *
640 * no interchange, use 1-by-1 pivot block
641 *
642  kp = k
643 *
644  ELSE
645 *
646  done = .false.
647 *
648 * Loop until pivot found
649 *
650  72 CONTINUE
651 *
652 * Begin pivot search loop body
653 *
654 *
655 * Copy column IMAX to column K+1 of W and update it
656 *
657  CALL zcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
658  CALL zcopy( n-imax+1, a( imax, imax ), 1,
659  $ w( imax, k+1 ), 1 )
660  IF( k.GT.1 )
661  $ CALL zgemv( 'No transpose', n-k+1, k-1, -cone,
662  $ a( k, 1 ), lda, w( imax, 1 ), ldw,
663  $ cone, w( k, k+1 ), 1 )
664 *
665 * JMAX is the column-index of the largest off-diagonal
666 * element in row IMAX, and ROWMAX is its absolute value.
667 * Determine both ROWMAX and JMAX.
668 *
669  IF( imax.NE.k ) THEN
670  jmax = k - 1 + izamax( imax-k, w( k, k+1 ), 1 )
671  rowmax = cabs1( w( jmax, k+1 ) )
672  ELSE
673  rowmax = zero
674  END IF
675 *
676  IF( imax.LT.n ) THEN
677  itemp = imax + izamax( n-imax, w( imax+1, k+1 ), 1)
678  dtemp = cabs1( w( itemp, k+1 ) )
679  IF( dtemp.GT.rowmax ) THEN
680  rowmax = dtemp
681  jmax = itemp
682  END IF
683  END IF
684 *
685 * Equivalent to testing for
686 * CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX
687 * (used to handle NaN and Inf)
688 *
689  IF( .NOT.( cabs1( w( imax, k+1 ) ).LT.alpha*rowmax ) )
690  $ THEN
691 *
692 * interchange rows and columns K and IMAX,
693 * use 1-by-1 pivot block
694 *
695  kp = imax
696 *
697 * copy column K+1 of W to column K of W
698 *
699  CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
700 *
701  done = .true.
702 *
703 * Equivalent to testing for ROWMAX.EQ.COLMAX,
704 * (used to handle NaN and Inf)
705 *
706  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
707  $ THEN
708 *
709 * interchange rows and columns K+1 and IMAX,
710 * use 2-by-2 pivot block
711 *
712  kp = imax
713  kstep = 2
714  done = .true.
715  ELSE
716 *
717 * Pivot not found: set params and repeat
718 *
719  p = imax
720  colmax = rowmax
721  imax = jmax
722 *
723 * Copy updated JMAXth (next IMAXth) column to Kth of W
724 *
725  CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
726 *
727  END IF
728 *
729 * End pivot search loop body
730 *
731  IF( .NOT. done ) GOTO 72
732 *
733  END IF
734 *
735 * ============================================================
736 *
737  kk = k + kstep - 1
738 *
739  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
740 *
741 * Copy non-updated column K to column P
742 *
743  CALL zcopy( p-k, a( k, k ), 1, a( p, k ), lda )
744  CALL zcopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
745 *
746 * Interchange rows K and P in first K columns of A
747 * and first K+1 columns of W
748 *
749  CALL zswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
750  CALL zswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
751  END IF
752 *
753 * Updated column KP is already stored in column KK of W
754 *
755  IF( kp.NE.kk ) THEN
756 *
757 * Copy non-updated column KK to column KP
758 *
759  a( kp, k ) = a( kk, k )
760  CALL zcopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
761  CALL zcopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
762 *
763 * Interchange rows KK and KP in first KK columns of A and W
764 *
765  CALL zswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
766  CALL zswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
767  END IF
768 *
769  IF( kstep.EQ.1 ) THEN
770 *
771 * 1-by-1 pivot block D(k): column k of W now holds
772 *
773 * W(k) = L(k)*D(k)
774 *
775 * where L(k) is the k-th column of L
776 *
777 * Store L(k) in column k of A
778 *
779  CALL zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
780  IF( k.LT.n ) THEN
781  IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
782  r1 = cone / a( k, k )
783  CALL zscal( n-k, r1, a( k+1, k ), 1 )
784  ELSE IF( a( k, k ).NE.czero ) THEN
785  DO 74 ii = k + 1, n
786  a( ii, k ) = a( ii, k ) / a( k, k )
787  74 CONTINUE
788  END IF
789  END IF
790 *
791  ELSE
792 *
793 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
794 *
795 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
796 *
797 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
798 * of L
799 *
800  IF( k.LT.n-1 ) THEN
801 *
802 * Store L(k) and L(k+1) in columns k and k+1 of A
803 *
804  d21 = w( k+1, k )
805  d11 = w( k+1, k+1 ) / d21
806  d22 = w( k, k ) / d21
807  t = cone / ( d11*d22-cone )
808  DO 80 j = k + 2, n
809  a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
810  $ d21 )
811  a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
812  $ d21 )
813  80 CONTINUE
814  END IF
815 *
816 * Copy D(k) to A
817 *
818  a( k, k ) = w( k, k )
819  a( k+1, k ) = w( k+1, k )
820  a( k+1, k+1 ) = w( k+1, k+1 )
821  END IF
822  END IF
823 *
824 * Store details of the interchanges in IPIV
825 *
826  IF( kstep.EQ.1 ) THEN
827  ipiv( k ) = kp
828  ELSE
829  ipiv( k ) = -p
830  ipiv( k+1 ) = -kp
831  END IF
832 *
833 * Increase K and return to the start of the main loop
834 *
835  k = k + kstep
836  GO TO 70
837 *
838  90 CONTINUE
839 *
840 * Update the lower triangle of A22 (= A(k:n,k:n)) as
841 *
842 * A22 := A22 - L21*D*L21**T = A22 - L21*W**T
843 *
844 * computing blocks of NB columns at a time
845 *
846  DO 110 j = k, n, nb
847  jb = min( nb, n-j+1 )
848 *
849 * Update the lower triangle of the diagonal block
850 *
851  DO 100 jj = j, j + jb - 1
852  CALL zgemv( 'No transpose', j+jb-jj, k-1, -cone,
853  $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
854  $ a( jj, jj ), 1 )
855  100 CONTINUE
856 *
857 * Update the rectangular subdiagonal block
858 *
859  IF( j+jb.LE.n )
860  $ CALL zgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
861  $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ), ldw,
862  $ cone, a( j+jb, j ), lda )
863  110 CONTINUE
864 *
865 * Put L21 in standard form by partially undoing the interchanges
866 * in columns 1:k-1
867 *
868  j = k - 1
869  120 CONTINUE
870 *
871  kstep = 1
872  jp1 = 1
873  jj = j
874  jp2 = ipiv( j )
875  IF( jp2.LT.0 ) THEN
876  jp2 = -jp2
877  j = j - 1
878  jp1 = -ipiv( j )
879  kstep = 2
880  END IF
881 *
882  j = j - 1
883  IF( jp2.NE.jj .AND. j.GE.1 )
884  $ CALL zswap( j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
885  jj = j + 1
886  IF( jp1.NE.jj .AND. kstep.EQ.2 )
887  $ CALL zswap( j, a( jp1, 1 ), lda, a( jj, 1 ), lda )
888  IF( j.GE.1 )
889  $ GO TO 120
890 *
891 * Set KB to the number of columns factorized
892 *
893  kb = k - 1
894 *
895  END IF
896  RETURN
897 *
898 * End of ZLASYF_ROOK
899 *
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:54
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zsycon ( character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
double precision  ANORM,
double precision  RCOND,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZSYCON

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

Purpose:
 ZSYCON estimates the reciprocal of the condition number (in the
 1-norm) of a complex symmetric matrix A using the factorization
 A = U*D*U**T or A = L*D*L**T computed by ZSYTRF.

 An estimate is obtained for norm(inv(A)), and the reciprocal of the
 condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
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]A
          A is COMPLEX*16 array, dimension (LDA,N)
          The block diagonal matrix D and the multipliers used to
          obtain the factor U or L as computed by ZSYTRF.
[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.
[in]ANORM
          ANORM is DOUBLE PRECISION
          The 1-norm of the original matrix A.
[out]RCOND
          RCOND is DOUBLE PRECISION
          The reciprocal of the condition number of the matrix A,
          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
          estimate of the 1-norm of inv(A) computed in this routine.
[out]WORK
          WORK is COMPLEX*16 array, dimension (2*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 127 of file zsycon.f.

127 *
128 * -- LAPACK computational routine (version 3.4.0) --
129 * -- LAPACK is a software package provided by Univ. of Tennessee, --
130 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
131 * November 2011
132 *
133 * .. Scalar Arguments ..
134  CHARACTER uplo
135  INTEGER info, lda, n
136  DOUBLE PRECISION anorm, rcond
137 * ..
138 * .. Array Arguments ..
139  INTEGER ipiv( * )
140  COMPLEX*16 a( lda, * ), work( * )
141 * ..
142 *
143 * =====================================================================
144 *
145 * .. Parameters ..
146  DOUBLE PRECISION one, zero
147  parameter( one = 1.0d+0, zero = 0.0d+0 )
148 * ..
149 * .. Local Scalars ..
150  LOGICAL upper
151  INTEGER i, kase
152  DOUBLE PRECISION ainvnm
153 * ..
154 * .. Local Arrays ..
155  INTEGER isave( 3 )
156 * ..
157 * .. External Functions ..
158  LOGICAL lsame
159  EXTERNAL lsame
160 * ..
161 * .. External Subroutines ..
162  EXTERNAL xerbla, zlacn2, zsytrs
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  ELSE IF( anorm.LT.zero ) THEN
180  info = -6
181  END IF
182  IF( info.NE.0 ) THEN
183  CALL xerbla( 'ZSYCON', -info )
184  RETURN
185  END IF
186 *
187 * Quick return if possible
188 *
189  rcond = zero
190  IF( n.EQ.0 ) THEN
191  rcond = one
192  RETURN
193  ELSE IF( anorm.LE.zero ) THEN
194  RETURN
195  END IF
196 *
197 * Check that the diagonal matrix D is nonsingular.
198 *
199  IF( upper ) THEN
200 *
201 * Upper triangular storage: examine D from bottom to top
202 *
203  DO 10 i = n, 1, -1
204  IF( ipiv( i ).GT.0 .AND. a( i, i ).EQ.zero )
205  $ RETURN
206  10 CONTINUE
207  ELSE
208 *
209 * Lower triangular storage: examine D from top to bottom.
210 *
211  DO 20 i = 1, n
212  IF( ipiv( i ).GT.0 .AND. a( i, i ).EQ.zero )
213  $ RETURN
214  20 CONTINUE
215  END IF
216 *
217 * Estimate the 1-norm of the inverse.
218 *
219  kase = 0
220  30 CONTINUE
221  CALL zlacn2( n, work( n+1 ), work, ainvnm, kase, isave )
222  IF( kase.NE.0 ) THEN
223 *
224 * Multiply by inv(L*D*L**T) or inv(U*D*U**T).
225 *
226  CALL zsytrs( uplo, n, 1, a, lda, ipiv, work, n, info )
227  GO TO 30
228  END IF
229 *
230 * Compute the estimate of the reciprocal condition number.
231 *
232  IF( ainvnm.NE.zero )
233  $ rcond = ( one / ainvnm ) / anorm
234 *
235  RETURN
236 *
237 * End of ZSYCON
238 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlacn2(N, V, X, EST, KASE, ISAVE)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: zlacn2.f:135
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zsytrs(UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZSYTRS
Definition: zsytrs.f:122

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zsycon_rook ( character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
double precision  ANORM,
double precision  RCOND,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZSYCON_ROOK

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

Purpose:
 ZSYCON_ROOK estimates the reciprocal of the condition number (in the
 1-norm) of a complex symmetric matrix A using the factorization
 A = U*D*U**T or A = L*D*L**T computed by ZSYTRF_ROOK.

 An estimate is obtained for norm(inv(A)), and the reciprocal of the
 condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
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]A
          A is COMPLEX*16 array, dimension (LDA,N)
          The block diagonal matrix D and the multipliers used to
          obtain the factor U or L as computed by ZSYTRF_ROOK.
[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.
[in]ANORM
          ANORM is DOUBLE PRECISION
          The 1-norm of the original matrix A.
[out]RCOND
          RCOND is DOUBLE PRECISION
          The reciprocal of the condition number of the matrix A,
          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
          estimate of the 1-norm of inv(A) computed in this routine.
[out]WORK
          WORK is COMPLEX*16 array, dimension (2*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015
Contributors:

November 2015, 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 141 of file zsycon_rook.f.

141 *
142 * -- LAPACK computational routine (version 3.6.0) --
143 * -- LAPACK is a software package provided by Univ. of Tennessee, --
144 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
145 * November 2015
146 *
147 * .. Scalar Arguments ..
148  CHARACTER uplo
149  INTEGER info, lda, n
150  DOUBLE PRECISION anorm, rcond
151 * ..
152 * .. Array Arguments ..
153  INTEGER ipiv( * )
154  COMPLEX*16 a( lda, * ), work( * )
155 * ..
156 *
157 * =====================================================================
158 *
159 * .. Parameters ..
160  DOUBLE PRECISION one, zero
161  parameter( one = 1.0d+0, zero = 0.0d+0 )
162  COMPLEX*16 czero
163  parameter( czero = ( 0.0d+0, 0.0d+0 ) )
164 * ..
165 * .. Local Scalars ..
166  LOGICAL upper
167  INTEGER i, kase
168  DOUBLE PRECISION ainvnm
169 * ..
170 * .. Local Arrays ..
171  INTEGER isave( 3 )
172 * ..
173 * .. External Functions ..
174  LOGICAL lsame
175  EXTERNAL lsame
176 * ..
177 * .. External Subroutines ..
178  EXTERNAL zlacn2, zsytrs_rook, xerbla
179 * ..
180 * .. Intrinsic Functions ..
181  INTRINSIC max
182 * ..
183 * .. Executable Statements ..
184 *
185 * Test the input parameters.
186 *
187  info = 0
188  upper = lsame( uplo, 'U' )
189  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
190  info = -1
191  ELSE IF( n.LT.0 ) THEN
192  info = -2
193  ELSE IF( lda.LT.max( 1, n ) ) THEN
194  info = -4
195  ELSE IF( anorm.LT.zero ) THEN
196  info = -6
197  END IF
198  IF( info.NE.0 ) THEN
199  CALL xerbla( 'ZSYCON_ROOK', -info )
200  RETURN
201  END IF
202 *
203 * Quick return if possible
204 *
205  rcond = zero
206  IF( n.EQ.0 ) THEN
207  rcond = one
208  RETURN
209  ELSE IF( anorm.LE.zero ) THEN
210  RETURN
211  END IF
212 *
213 * Check that the diagonal matrix D is nonsingular.
214 *
215  IF( upper ) THEN
216 *
217 * Upper triangular storage: examine D from bottom to top
218 *
219  DO 10 i = n, 1, -1
220  IF( ipiv( i ).GT.0 .AND. a( i, i ).EQ.czero )
221  $ RETURN
222  10 CONTINUE
223  ELSE
224 *
225 * Lower triangular storage: examine D from top to bottom.
226 *
227  DO 20 i = 1, n
228  IF( ipiv( i ).GT.0 .AND. a( i, i ).EQ.czero )
229  $ RETURN
230  20 CONTINUE
231  END IF
232 *
233 * Estimate the 1-norm of the inverse.
234 *
235  kase = 0
236  30 CONTINUE
237  CALL zlacn2( n, work( n+1 ), work, ainvnm, kase, isave )
238  IF( kase.NE.0 ) THEN
239 *
240 * Multiply by inv(L*D*L**T) or inv(U*D*U**T).
241 *
242  CALL zsytrs_rook( uplo, n, 1, a, lda, ipiv, work, n, info )
243  GO TO 30
244  END IF
245 *
246 * Compute the estimate of the reciprocal condition number.
247 *
248  IF( ainvnm.NE.zero )
249  $ rcond = ( one / ainvnm ) / anorm
250 *
251  RETURN
252 *
253 * End of ZSYCON_ROOK
254 *
subroutine zsytrs_rook(UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZSYTRS_ROOK
Definition: zsytrs_rook.f:138
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlacn2(N, V, X, EST, KASE, ISAVE)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: zlacn2.f:135
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zsyconv ( character  UPLO,
character  WAY,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex*16, dimension( * )  E,
integer  INFO 
)

ZSYCONV

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

Purpose:
 ZSYCONV converts A given by ZHETRF into L and D or vice-versa.
 Get nondiagonal elements of D (returned in workspace) and 
 apply or reverse permutation done in TRF.
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]WAY
          WAY is CHARACTER*1
          = 'C': Convert 
          = 'R': Revert
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          The block diagonal matrix D and the multipliers used to
          obtain the factor U or L as computed by ZSYTRF.
[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.
[out]E
          E is COMPLEX*16 array, dimension (N)
          E stores the supdiagonal/subdiagonal of the symmetric 1-by-1
          or 2-by-2 block diagonal matrix D in LDLT.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015

Definition at line 116 of file zsyconv.f.

116 *
117 * -- LAPACK computational routine (version 3.6.0) --
118 * -- LAPACK is a software package provided by Univ. of Tennessee, --
119 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
120 * November 2015
121 *
122 * .. Scalar Arguments ..
123  CHARACTER uplo, way
124  INTEGER info, lda, n
125 * ..
126 * .. Array Arguments ..
127  INTEGER ipiv( * )
128  COMPLEX*16 a( lda, * ), e( * )
129 * ..
130 *
131 * =====================================================================
132 *
133 * .. Parameters ..
134  COMPLEX*16 zero
135  parameter( zero = (0.0d+0,0.0d+0) )
136 * ..
137 * .. External Functions ..
138  LOGICAL lsame
139  EXTERNAL lsame
140 *
141 * .. External Subroutines ..
142  EXTERNAL xerbla
143 * .. Local Scalars ..
144  LOGICAL upper, convert
145  INTEGER i, ip, j
146  COMPLEX*16 temp
147 * ..
148 * .. Executable Statements ..
149 *
150  info = 0
151  upper = lsame( uplo, 'U' )
152  convert = lsame( way, 'C' )
153  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
154  info = -1
155  ELSE IF( .NOT.convert .AND. .NOT.lsame( way, 'R' ) ) THEN
156  info = -2
157  ELSE IF( n.LT.0 ) THEN
158  info = -3
159  ELSE IF( lda.LT.max( 1, n ) ) THEN
160  info = -5
161 
162  END IF
163  IF( info.NE.0 ) THEN
164  CALL xerbla( 'ZSYCONV', -info )
165  RETURN
166  END IF
167 *
168 * Quick return if possible
169 *
170  IF( n.EQ.0 )
171  $ RETURN
172 *
173  IF( upper ) THEN
174 *
175 * A is UPPER
176 *
177  IF ( convert ) THEN
178 *
179 * Convert A (A is upper)
180 *
181 * Convert VALUE
182 *
183  i=n
184  e(1)=zero
185  DO WHILE ( i .GT. 1 )
186  IF( ipiv(i) .LT. 0 ) THEN
187  e(i)=a(i-1,i)
188  e(i-1)=zero
189  a(i-1,i)=zero
190  i=i-1
191  ELSE
192  e(i)=zero
193  ENDIF
194  i=i-1
195  END DO
196 *
197 * Convert PERMUTATIONS
198 *
199  i=n
200  DO WHILE ( i .GE. 1 )
201  IF( ipiv(i) .GT. 0) THEN
202  ip=ipiv(i)
203  IF( i .LT. n) THEN
204  DO 12 j= i+1,n
205  temp=a(ip,j)
206  a(ip,j)=a(i,j)
207  a(i,j)=temp
208  12 CONTINUE
209  ENDIF
210  ELSE
211  ip=-ipiv(i)
212  IF( i .LT. n) THEN
213  DO 13 j= i+1,n
214  temp=a(ip,j)
215  a(ip,j)=a(i-1,j)
216  a(i-1,j)=temp
217  13 CONTINUE
218  ENDIF
219  i=i-1
220  ENDIF
221  i=i-1
222  END DO
223 *
224  ELSE
225 *
226 * Revert A (A is upper)
227 *
228 * Revert PERMUTATIONS
229 *
230  i=1
231  DO WHILE ( i .LE. n )
232  IF( ipiv(i) .GT. 0 ) THEN
233  ip=ipiv(i)
234  IF( i .LT. n) THEN
235  DO j= i+1,n
236  temp=a(ip,j)
237  a(ip,j)=a(i,j)
238  a(i,j)=temp
239  END DO
240  ENDIF
241  ELSE
242  ip=-ipiv(i)
243  i=i+1
244  IF( i .LT. n) THEN
245  DO j= i+1,n
246  temp=a(ip,j)
247  a(ip,j)=a(i-1,j)
248  a(i-1,j)=temp
249  END DO
250  ENDIF
251  ENDIF
252  i=i+1
253  END DO
254 *
255 * Revert VALUE
256 *
257  i=n
258  DO WHILE ( i .GT. 1 )
259  IF( ipiv(i) .LT. 0 ) THEN
260  a(i-1,i)=e(i)
261  i=i-1
262  ENDIF
263  i=i-1
264  END DO
265  END IF
266 *
267  ELSE
268 *
269 * A is LOWER
270 *
271  IF ( convert ) THEN
272 *
273 * Convert A (A is lower)
274 *
275 * Convert VALUE
276 *
277  i=1
278  e(n)=zero
279  DO WHILE ( i .LE. n )
280  IF( i.LT.n .AND. ipiv(i) .LT. 0 ) THEN
281  e(i)=a(i+1,i)
282  e(i+1)=zero
283  a(i+1,i)=zero
284  i=i+1
285  ELSE
286  e(i)=zero
287  ENDIF
288  i=i+1
289  END DO
290 *
291 * Convert PERMUTATIONS
292 *
293  i=1
294  DO WHILE ( i .LE. n )
295  IF( ipiv(i) .GT. 0 ) THEN
296  ip=ipiv(i)
297  IF (i .GT. 1) THEN
298  DO 22 j= 1,i-1
299  temp=a(ip,j)
300  a(ip,j)=a(i,j)
301  a(i,j)=temp
302  22 CONTINUE
303  ENDIF
304  ELSE
305  ip=-ipiv(i)
306  IF (i .GT. 1) THEN
307  DO 23 j= 1,i-1
308  temp=a(ip,j)
309  a(ip,j)=a(i+1,j)
310  a(i+1,j)=temp
311  23 CONTINUE
312  ENDIF
313  i=i+1
314  ENDIF
315  i=i+1
316  END DO
317 *
318  ELSE
319 *
320 * Revert A (A is lower)
321 *
322 * Revert PERMUTATIONS
323 *
324  i=n
325  DO WHILE ( i .GE. 1 )
326  IF( ipiv(i) .GT. 0 ) THEN
327  ip=ipiv(i)
328  IF (i .GT. 1) THEN
329  DO j= 1,i-1
330  temp=a(i,j)
331  a(i,j)=a(ip,j)
332  a(ip,j)=temp
333  END DO
334  ENDIF
335  ELSE
336  ip=-ipiv(i)
337  i=i-1
338  IF (i .GT. 1) THEN
339  DO j= 1,i-1
340  temp=a(i+1,j)
341  a(i+1,j)=a(ip,j)
342  a(ip,j)=temp
343  END DO
344  ENDIF
345  ENDIF
346  i=i-1
347  END DO
348 *
349 * Revert VALUE
350 *
351  i=1
352  DO WHILE ( i .LE. n-1 )
353  IF( ipiv(i) .LT. 0 ) THEN
354  a(i+1,i)=e(i)
355  i=i+1
356  ENDIF
357  i=i+1
358  END DO
359  END IF
360  END IF
361 *
362  RETURN
363 *
364 * End of ZSYCONV
365 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zsyequb ( character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  S,
double precision  SCOND,
double precision  AMAX,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZSYEQUB

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

Purpose:
 ZSYEQUB computes row and column scalings intended to equilibrate a
 symmetric matrix A and reduce its condition number
 (with respect to the two-norm).  S contains the scale factors,
 S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with
 elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal.  This
 choice of S puts the condition number of B within a factor N of the
 smallest possible condition number over all possible diagonal
 scalings.
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]A
          A is COMPLEX*16 array, dimension (LDA,N)
          The N-by-N symmetric matrix whose scaling
          factors are to be computed.  Only the diagonal elements of A
          are referenced.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]S
          S is DOUBLE PRECISION array, dimension (N)
          If INFO = 0, S contains the scale factors for A.
[out]SCOND
          SCOND is DOUBLE PRECISION
          If INFO = 0, S contains the ratio of the smallest S(i) to
          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
          large nor too small, it is not worth scaling by S.
[out]AMAX
          AMAX is DOUBLE PRECISION
          Absolute value of largest matrix element.  If AMAX is very
          close to overflow or very close to underflow, the matrix
          should be scaled.
[out]WORK
          WORK is COMPLEX*16 array, dimension (3*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, the i-th diagonal element is nonpositive.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
References:
Livne, O.E. and Golub, G.H., "Scaling by Binormalization",
Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004.
DOI 10.1023/B:NUMA.0000016606.32820.69
Tech report version: http://ruready.utah.edu/archive/papers/bin.pdf

Definition at line 138 of file zsyequb.f.

138 *
139 * -- LAPACK computational routine (version 3.4.0) --
140 * -- LAPACK is a software package provided by Univ. of Tennessee, --
141 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
142 * November 2011
143 *
144 * .. Scalar Arguments ..
145  INTEGER info, lda, n
146  DOUBLE PRECISION amax, scond
147  CHARACTER uplo
148 * ..
149 * .. Array Arguments ..
150  COMPLEX*16 a( lda, * ), work( * )
151  DOUBLE PRECISION s( * )
152 * ..
153 *
154 * =====================================================================
155 *
156 * .. Parameters ..
157  DOUBLE PRECISION one, zero
158  parameter( one = 1.0d0, zero = 0.0d0 )
159  INTEGER max_iter
160  parameter( max_iter = 100 )
161 * ..
162 * .. Local Scalars ..
163  INTEGER i, j, iter
164  DOUBLE PRECISION avg, std, tol, c0, c1, c2, t, u, si, d, base,
165  $ smin, smax, smlnum, bignum, scale, sumsq
166  LOGICAL up
167  COMPLEX*16 zdum
168 * ..
169 * .. External Functions ..
170  DOUBLE PRECISION dlamch
171  LOGICAL lsame
172  EXTERNAL dlamch, lsame
173 * ..
174 * .. External Subroutines ..
175  EXTERNAL zlassq
176 * ..
177 * .. Intrinsic Functions ..
178  INTRINSIC abs, dble, dimag, int, log, max, min, sqrt
179 * ..
180 * .. Statement Functions ..
181  DOUBLE PRECISION cabs1
182 * ..
183 * Statement Function Definitions
184  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
185 * ..
186 * .. Executable Statements ..
187 *
188 * Test the input parameters.
189 *
190  info = 0
191  IF ( .NOT. ( lsame( uplo, 'U' ) .OR. lsame( uplo, 'L' ) ) ) THEN
192  info = -1
193  ELSE IF ( n .LT. 0 ) THEN
194  info = -2
195  ELSE IF ( lda .LT. max( 1, n ) ) THEN
196  info = -4
197  END IF
198  IF ( info .NE. 0 ) THEN
199  CALL xerbla( 'ZSYEQUB', -info )
200  RETURN
201  END IF
202 
203  up = lsame( uplo, 'U' )
204  amax = zero
205 *
206 * Quick return if possible.
207 *
208  IF ( n .EQ. 0 ) THEN
209  scond = one
210  RETURN
211  END IF
212 
213  DO i = 1, n
214  s( i ) = zero
215  END DO
216 
217  amax = zero
218  IF ( up ) THEN
219  DO j = 1, n
220  DO i = 1, j-1
221  s( i ) = max( s( i ), cabs1( a( i, j ) ) )
222  s( j ) = max( s( j ), cabs1( a( i, j ) ) )
223  amax = max( amax, cabs1( a( i, j ) ) )
224  END DO
225  s( j ) = max( s( j ), cabs1( a( j, j) ) )
226  amax = max( amax, cabs1( a( j, j ) ) )
227  END DO
228  ELSE
229  DO j = 1, n
230  s( j ) = max( s( j ), cabs1( a( j, j ) ) )
231  amax = max( amax, cabs1( a( j, j ) ) )
232  DO i = j+1, n
233  s( i ) = max( s( i ), cabs1( a( i, j ) ) )
234  s( j ) = max( s( j ), cabs1(a( i, j ) ) )
235  amax = max( amax, cabs1( a( i, j ) ) )
236  END DO
237  END DO
238  END IF
239  DO j = 1, n
240  s( j ) = 1.0d+0 / s( j )
241  END DO
242 
243  tol = one / sqrt( 2.0d0 * n )
244 
245  DO iter = 1, max_iter
246  scale = 0.0d+0
247  sumsq = 0.0d+0
248 * beta = |A|s
249  DO i = 1, n
250  work( i ) = zero
251  END DO
252  IF ( up ) THEN
253  DO j = 1, n
254  DO i = 1, j-1
255  t = cabs1( a( i, j ) )
256  work( i ) = work( i ) + cabs1( a( i, j ) ) * s( j )
257  work( j ) = work( j ) + cabs1( a( i, j ) ) * s( i )
258  END DO
259  work( j ) = work( j ) + cabs1( a( j, j ) ) * s( j )
260  END DO
261  ELSE
262  DO j = 1, n
263  work( j ) = work( j ) + cabs1( a( j, j ) ) * s( j )
264  DO i = j+1, n
265  t = cabs1( a( i, j ) )
266  work( i ) = work( i ) + cabs1( a( i, j ) ) * s( j )
267  work( j ) = work( j ) + cabs1( a( i, j ) ) * s( i )
268  END DO
269  END DO
270  END IF
271 
272 * avg = s^T beta / n
273  avg = 0.0d+0
274  DO i = 1, n
275  avg = avg + s( i )*work( i )
276  END DO
277  avg = avg / n
278 
279  std = 0.0d+0
280  DO i = n+1, 2*n
281  work( i ) = s( i-n ) * work( i-n ) - avg
282  END DO
283  CALL zlassq( n, work( n+1 ), 1, scale, sumsq )
284  std = scale * sqrt( sumsq / n )
285 
286  IF ( std .LT. tol * avg ) GOTO 999
287 
288  DO i = 1, n
289  t = cabs1( a( i, i ) )
290  si = s( i )
291  c2 = ( n-1 ) * t
292  c1 = ( n-2 ) * ( work( i ) - t*si )
293  c0 = -(t*si)*si + 2*work( i )*si - n*avg
294  d = c1*c1 - 4*c0*c2
295 
296  IF ( d .LE. 0 ) THEN
297  info = -1
298  RETURN
299  END IF
300  si = -2*c0 / ( c1 + sqrt( d ) )
301 
302  d = si - s( i )
303  u = zero
304  IF ( up ) THEN
305  DO j = 1, i
306  t = cabs1( a( j, i ) )
307  u = u + s( j )*t
308  work( j ) = work( j ) + d*t
309  END DO
310  DO j = i+1,n
311  t = cabs1( a( i, j ) )
312  u = u + s( j )*t
313  work( j ) = work( j ) + d*t
314  END DO
315  ELSE
316  DO j = 1, i
317  t = cabs1( a( i, j ) )
318  u = u + s( j )*t
319  work( j ) = work( j ) + d*t
320  END DO
321  DO j = i+1,n
322  t = cabs1( a( j, i ) )
323  u = u + s( j )*t
324  work( j ) = work( j ) + d*t
325  END DO
326  END IF
327  avg = avg + ( u + work( i ) ) * d / n
328  s( i ) = si
329  END DO
330  END DO
331 
332  999 CONTINUE
333 
334  smlnum = dlamch( 'SAFEMIN' )
335  bignum = one / smlnum
336  smin = bignum
337  smax = zero
338  t = one / sqrt( avg )
339  base = dlamch( 'B' )
340  u = one / log( base )
341  DO i = 1, n
342  s( i ) = base ** int( u * log( s( i ) * t ) )
343  smin = min( smin, s( i ) )
344  smax = max( smax, s( i ) )
345  END DO
346  scond = max( smin, smlnum ) / min( smax, bignum )
347 *
subroutine zlassq(N, X, INCX, SCALE, SUMSQ)
ZLASSQ updates a sum of squares represented in scaled form.
Definition: zlassq.f:108
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zsyrfs ( character  UPLO,
integer  N,
integer  NRHS,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldaf, * )  AF,
integer  LDAF,
integer, dimension( * )  IPIV,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( ldx, * )  X,
integer  LDX,
double precision, dimension( * )  FERR,
double precision, dimension( * )  BERR,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZSYRFS

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

Purpose:
 ZSYRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is symmetric indefinite, and
 provides error bounds and backward error estimates for the solution.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrices B and X.  NRHS >= 0.
[in]A
          A is COMPLEX*16 array, dimension (LDA,N)
          The symmetric matrix A.  If UPLO = 'U', the leading N-by-N
          upper triangular part of A contains the upper triangular part
          of the matrix A, and the strictly lower triangular part of A
          is not referenced.  If UPLO = 'L', the leading N-by-N lower
          triangular part of A contains the lower triangular part of
          the matrix A, and the strictly upper triangular part of A is
          not referenced.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]AF
          AF is COMPLEX*16 array, dimension (LDAF,N)
          The factored form of the matrix A.  AF contains the block
          diagonal matrix D and the multipliers used to obtain the
          factor U or L from the factorization A = U*D*U**T or
          A = L*D*L**T as computed by ZSYTRF.
[in]LDAF
          LDAF is INTEGER
          The leading dimension of the array AF.  LDAF >= 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.
[in]B
          B is COMPLEX*16 array, dimension (LDB,NRHS)
          The right hand side matrix B.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,N).
[in,out]X
          X is COMPLEX*16 array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by ZSYTRS.
          On exit, the improved solution matrix X.
[in]LDX
          LDX is INTEGER
          The leading dimension of the array X.  LDX >= max(1,N).
[out]FERR
          FERR is DOUBLE PRECISION array, dimension (NRHS)
          The estimated forward error bound for each solution vector
          X(j) (the j-th column of the solution matrix X).
          If XTRUE is the true solution corresponding to X(j), FERR(j)
          is an estimated upper bound for the magnitude of the largest
          element in (X(j) - XTRUE) divided by the magnitude of the
          largest element in X(j).  The estimate is as reliable as
          the estimate for RCOND, and is almost always a slight
          overestimate of the true error.
[out]BERR
          BERR is DOUBLE PRECISION array, dimension (NRHS)
          The componentwise relative backward error of each solution
          vector X(j) (i.e., the smallest relative change in
          any element of A or B that makes X(j) an exact solution).
[out]WORK
          WORK is COMPLEX*16 array, dimension (2*N)
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Internal Parameters:
  ITMAX is the maximum number of steps of iterative refinement.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 194 of file zsyrfs.f.

194 *
195 * -- LAPACK computational routine (version 3.4.0) --
196 * -- LAPACK is a software package provided by Univ. of Tennessee, --
197 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
198 * November 2011
199 *
200 * .. Scalar Arguments ..
201  CHARACTER uplo
202  INTEGER info, lda, ldaf, ldb, ldx, n, nrhs
203 * ..
204 * .. Array Arguments ..
205  INTEGER ipiv( * )
206  DOUBLE PRECISION berr( * ), ferr( * ), rwork( * )
207  COMPLEX*16 a( lda, * ), af( ldaf, * ), b( ldb, * ),
208  $ work( * ), x( ldx, * )
209 * ..
210 *
211 * =====================================================================
212 *
213 * .. Parameters ..
214  INTEGER itmax
215  parameter( itmax = 5 )
216  DOUBLE PRECISION zero
217  parameter( zero = 0.0d+0 )
218  COMPLEX*16 one
219  parameter( one = ( 1.0d+0, 0.0d+0 ) )
220  DOUBLE PRECISION two
221  parameter( two = 2.0d+0 )
222  DOUBLE PRECISION three
223  parameter( three = 3.0d+0 )
224 * ..
225 * .. Local Scalars ..
226  LOGICAL upper
227  INTEGER count, i, j, k, kase, nz
228  DOUBLE PRECISION eps, lstres, s, safe1, safe2, safmin, xk
229  COMPLEX*16 zdum
230 * ..
231 * .. Local Arrays ..
232  INTEGER isave( 3 )
233 * ..
234 * .. External Subroutines ..
235  EXTERNAL xerbla, zaxpy, zcopy, zlacn2, zsymv, zsytrs
236 * ..
237 * .. Intrinsic Functions ..
238  INTRINSIC abs, dble, dimag, max
239 * ..
240 * .. External Functions ..
241  LOGICAL lsame
242  DOUBLE PRECISION dlamch
243  EXTERNAL lsame, dlamch
244 * ..
245 * .. Statement Functions ..
246  DOUBLE PRECISION cabs1
247 * ..
248 * .. Statement Function definitions ..
249  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
250 * ..
251 * .. Executable Statements ..
252 *
253 * Test the input parameters.
254 *
255  info = 0
256  upper = lsame( uplo, 'U' )
257  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
258  info = -1
259  ELSE IF( n.LT.0 ) THEN
260  info = -2
261  ELSE IF( nrhs.LT.0 ) THEN
262  info = -3
263  ELSE IF( lda.LT.max( 1, n ) ) THEN
264  info = -5
265  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
266  info = -7
267  ELSE IF( ldb.LT.max( 1, n ) ) THEN
268  info = -10
269  ELSE IF( ldx.LT.max( 1, n ) ) THEN
270  info = -12
271  END IF
272  IF( info.NE.0 ) THEN
273  CALL xerbla( 'ZSYRFS', -info )
274  RETURN
275  END IF
276 *
277 * Quick return if possible
278 *
279  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
280  DO 10 j = 1, nrhs
281  ferr( j ) = zero
282  berr( j ) = zero
283  10 CONTINUE
284  RETURN
285  END IF
286 *
287 * NZ = maximum number of nonzero elements in each row of A, plus 1
288 *
289  nz = n + 1
290  eps = dlamch( 'Epsilon' )
291  safmin = dlamch( 'Safe minimum' )
292  safe1 = nz*safmin
293  safe2 = safe1 / eps
294 *
295 * Do for each right hand side
296 *
297  DO 140 j = 1, nrhs
298 *
299  count = 1
300  lstres = three
301  20 CONTINUE
302 *
303 * Loop until stopping criterion is satisfied.
304 *
305 * Compute residual R = B - A * X
306 *
307  CALL zcopy( n, b( 1, j ), 1, work, 1 )
308  CALL zsymv( uplo, n, -one, a, lda, x( 1, j ), 1, one, work, 1 )
309 *
310 * Compute componentwise relative backward error from formula
311 *
312 * max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
313 *
314 * where abs(Z) is the componentwise absolute value of the matrix
315 * or vector Z. If the i-th component of the denominator is less
316 * than SAFE2, then SAFE1 is added to the i-th components of the
317 * numerator and denominator before dividing.
318 *
319  DO 30 i = 1, n
320  rwork( i ) = cabs1( b( i, j ) )
321  30 CONTINUE
322 *
323 * Compute abs(A)*abs(X) + abs(B).
324 *
325  IF( upper ) THEN
326  DO 50 k = 1, n
327  s = zero
328  xk = cabs1( x( k, j ) )
329  DO 40 i = 1, k - 1
330  rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
331  s = s + cabs1( a( i, k ) )*cabs1( x( i, j ) )
332  40 CONTINUE
333  rwork( k ) = rwork( k ) + cabs1( a( k, k ) )*xk + s
334  50 CONTINUE
335  ELSE
336  DO 70 k = 1, n
337  s = zero
338  xk = cabs1( x( k, j ) )
339  rwork( k ) = rwork( k ) + cabs1( a( k, k ) )*xk
340  DO 60 i = k + 1, n
341  rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
342  s = s + cabs1( a( i, k ) )*cabs1( x( i, j ) )
343  60 CONTINUE
344  rwork( k ) = rwork( k ) + s
345  70 CONTINUE
346  END IF
347  s = zero
348  DO 80 i = 1, n
349  IF( rwork( i ).GT.safe2 ) THEN
350  s = max( s, cabs1( work( i ) ) / rwork( i ) )
351  ELSE
352  s = max( s, ( cabs1( work( i ) )+safe1 ) /
353  $ ( rwork( i )+safe1 ) )
354  END IF
355  80 CONTINUE
356  berr( j ) = s
357 *
358 * Test stopping criterion. Continue iterating if
359 * 1) The residual BERR(J) is larger than machine epsilon, and
360 * 2) BERR(J) decreased by at least a factor of 2 during the
361 * last iteration, and
362 * 3) At most ITMAX iterations tried.
363 *
364  IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
365  $ count.LE.itmax ) THEN
366 *
367 * Update solution and try again.
368 *
369  CALL zsytrs( uplo, n, 1, af, ldaf, ipiv, work, n, info )
370  CALL zaxpy( n, one, work, 1, x( 1, j ), 1 )
371  lstres = berr( j )
372  count = count + 1
373  GO TO 20
374  END IF
375 *
376 * Bound error from formula
377 *
378 * norm(X - XTRUE) / norm(X) .le. FERR =
379 * norm( abs(inv(A))*
380 * ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
381 *
382 * where
383 * norm(Z) is the magnitude of the largest component of Z
384 * inv(A) is the inverse of A
385 * abs(Z) is the componentwise absolute value of the matrix or
386 * vector Z
387 * NZ is the maximum number of nonzeros in any row of A, plus 1
388 * EPS is machine epsilon
389 *
390 * The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
391 * is incremented by SAFE1 if the i-th component of
392 * abs(A)*abs(X) + abs(B) is less than SAFE2.
393 *
394 * Use ZLACN2 to estimate the infinity-norm of the matrix
395 * inv(A) * diag(W),
396 * where W = abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) )))
397 *
398  DO 90 i = 1, n
399  IF( rwork( i ).GT.safe2 ) THEN
400  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
401  ELSE
402  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
403  $ safe1
404  END IF
405  90 CONTINUE
406 *
407  kase = 0
408  100 CONTINUE
409  CALL zlacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
410  IF( kase.NE.0 ) THEN
411  IF( kase.EQ.1 ) THEN
412 *
413 * Multiply by diag(W)*inv(A**T).
414 *
415  CALL zsytrs( uplo, n, 1, af, ldaf, ipiv, work, n, info )
416  DO 110 i = 1, n
417  work( i ) = rwork( i )*work( i )
418  110 CONTINUE
419  ELSE IF( kase.EQ.2 ) THEN
420 *
421 * Multiply by inv(A)*diag(W).
422 *
423  DO 120 i = 1, n
424  work( i ) = rwork( i )*work( i )
425  120 CONTINUE
426  CALL zsytrs( uplo, n, 1, af, ldaf, ipiv, work, n, info )
427  END IF
428  GO TO 100
429  END IF
430 *
431 * Normalize error.
432 *
433  lstres = zero
434  DO 130 i = 1, n
435  lstres = max( lstres, cabs1( x( i, j ) ) )
436  130 CONTINUE
437  IF( lstres.NE.zero )
438  $ ferr( j ) = ferr( j ) / lstres
439 *
440  140 CONTINUE
441 *
442  RETURN
443 *
444 * End of ZSYRFS
445 *
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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
subroutine zlacn2(N, V, X, EST, KASE, ISAVE)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: zlacn2.f:135
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:53
subroutine zsytrs(UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZSYTRS
Definition: zsytrs.f:122

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zsyrfsx ( character  UPLO,
character  EQUED,
integer  N,
integer  NRHS,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldaf, * )  AF,
integer  LDAF,
integer, dimension( * )  IPIV,
double precision, dimension( * )  S,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( ldx, * )  X,
integer  LDX,
double precision  RCOND,
double precision, dimension( * )  BERR,
integer  N_ERR_BNDS,
double precision, dimension( nrhs, * )  ERR_BNDS_NORM,
double precision, dimension( nrhs, * )  ERR_BNDS_COMP,
integer  NPARAMS,
double precision, dimension( * )  PARAMS,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZSYRFSX

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

Purpose:
    ZSYRFSX improves the computed solution to a system of linear
    equations when the coefficient matrix is symmetric indefinite, and
    provides error bounds and backward error estimates for the
    solution.  In addition to normwise error bound, the code provides
    maximum componentwise error bound if possible.  See comments for
    ERR_BNDS_NORM and ERR_BNDS_COMP for details of the error bounds.

    The original system of linear equations may have been equilibrated
    before calling this routine, as described by arguments EQUED and S
    below. In this case, the solution and error bounds returned are
    for the original unequilibrated system.
     Some optional parameters are bundled in the PARAMS array.  These
     settings determine how refinement is performed, but often the
     defaults are acceptable.  If the defaults are acceptable, users
     can pass NPARAMS = 0 which prevents the source code from accessing
     the PARAMS argument.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
       = 'U':  Upper triangle of A is stored;
       = 'L':  Lower triangle of A is stored.
[in]EQUED
          EQUED is CHARACTER*1
     Specifies the form of equilibration that was done to A
     before calling this routine. This is needed to compute
     the solution and error bounds correctly.
       = 'N':  No equilibration
       = 'Y':  Both row and column equilibration, i.e., A has been
               replaced by diag(S) * A * diag(S).
               The right hand side B has been changed accordingly.
[in]N
          N is INTEGER
     The order of the matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
     The number of right hand sides, i.e., the number of columns
     of the matrices B and X.  NRHS >= 0.
[in]A
          A is COMPLEX*16 array, dimension (LDA,N)
     The symmetric matrix A.  If UPLO = 'U', the leading N-by-N
     upper triangular part of A contains the upper triangular
     part of the matrix A, and the strictly lower triangular
     part of A is not referenced.  If UPLO = 'L', the leading
     N-by-N lower triangular part of A contains the lower
     triangular part of the matrix A, and the strictly upper
     triangular part of A is not referenced.
[in]LDA
          LDA is INTEGER
     The leading dimension of the array A.  LDA >= max(1,N).
[in]AF
          AF is COMPLEX*16 array, dimension (LDAF,N)
     The factored form of the matrix A.  AF contains the block
     diagonal matrix D and the multipliers used to obtain the
     factor U or L from the factorization A = U*D*U**T or A =
     L*D*L**T as computed by DSYTRF.
[in]LDAF
          LDAF is INTEGER
     The leading dimension of the array AF.  LDAF >= max(1,N).
[in]IPIV
          IPIV is INTEGER array, dimension (N)
     Details of the interchanges and the block structure of D
     as determined by DSYTRF.
[in,out]S
          S is DOUBLE PRECISION array, dimension (N)
     The scale factors for A.  If EQUED = 'Y', A is multiplied on
     the left and right by diag(S).  S is an input argument if FACT =
     'F'; otherwise, S is an output argument.  If FACT = 'F' and EQUED
     = 'Y', each element of S must be positive.  If S is output, each
     element of S is a power of the radix. If S is input, each element
     of S should be a power of the radix to ensure a reliable solution
     and error estimates. Scaling by powers of the radix does not cause
     rounding errors unless the result underflows or overflows.
     Rounding errors during scaling lead to refining with a matrix that
     is not equivalent to the input matrix, producing error estimates
     that may not be reliable.
[in]B
          B is COMPLEX*16 array, dimension (LDB,NRHS)
     The right hand side matrix B.
[in]LDB
          LDB is INTEGER
     The leading dimension of the array B.  LDB >= max(1,N).
[in,out]X
          X is COMPLEX*16 array, dimension (LDX,NRHS)
     On entry, the solution matrix X, as computed by DGETRS.
     On exit, the improved solution matrix X.
[in]LDX
          LDX is INTEGER
     The leading dimension of the array X.  LDX >= max(1,N).
[out]RCOND
          RCOND is DOUBLE PRECISION
     Reciprocal scaled condition number.  This is an estimate of the
     reciprocal Skeel condition number of the matrix A after
     equilibration (if done).  If this is less than the machine
     precision (in particular, if it is zero), the matrix is singular
     to working precision.  Note that the error may still be small even
     if this number is very small and the matrix appears ill-
     conditioned.
[out]BERR
          BERR is DOUBLE PRECISION array, dimension (NRHS)
     Componentwise relative backward error.  This is the
     componentwise relative backward error of each solution vector X(j)
     (i.e., the smallest relative change in any element of A or B that
     makes X(j) an exact solution).
[in]N_ERR_BNDS
          N_ERR_BNDS is INTEGER
     Number of error bounds to return for each right hand side
     and each type (normwise or componentwise).  See ERR_BNDS_NORM and
     ERR_BNDS_COMP below.
[out]ERR_BNDS_NORM
          ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
     For each right-hand side, this array contains information about
     various error bounds and condition numbers corresponding to the
     normwise relative error, which is defined as follows:

     Normwise relative error in the ith solution vector:
             max_j (abs(XTRUE(j,i) - X(j,i)))
            ------------------------------
                  max_j abs(X(j,i))

     The array is indexed by the type of error information as described
     below. There currently are up to three pieces of information
     returned.

     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
     right-hand side.

     The second index in ERR_BNDS_NORM(:,err) contains the following
     three fields:
     err = 1 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * dlamch('Epsilon').

     err = 2 "Guaranteed" error bound: The estimated forward error,
              almost certainly within a factor of 10 of the true error
              so long as the next entry is greater than the threshold
              sqrt(n) * dlamch('Epsilon'). This error bound should only
              be trusted if the previous boolean is true.

     err = 3  Reciprocal condition number: Estimated normwise
              reciprocal condition number.  Compared with the threshold
              sqrt(n) * dlamch('Epsilon') to determine if the error
              estimate is "guaranteed". These reciprocal condition
              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
              appropriately scaled matrix Z.
              Let Z = S*A, where S scales each row by a power of the
              radix so all absolute row sums of Z are approximately 1.

     See Lapack Working Note 165 for further details and extra
     cautions.
[out]ERR_BNDS_COMP
          ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
     For each right-hand side, this array contains information about
     various error bounds and condition numbers corresponding to the
     componentwise relative error, which is defined as follows:

     Componentwise relative error in the ith solution vector:
                    abs(XTRUE(j,i) - X(j,i))
             max_j ----------------------
                         abs(X(j,i))

     The array is indexed by the right-hand side i (on which the
     componentwise relative error depends), and the type of error
     information as described below. There currently are up to three
     pieces of information returned for each right-hand side. If
     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS .LT. 3, then at most
     the first (:,N_ERR_BNDS) entries are returned.

     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
     right-hand side.

     The second index in ERR_BNDS_COMP(:,err) contains the following
     three fields:
     err = 1 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * dlamch('Epsilon').

     err = 2 "Guaranteed" error bound: The estimated forward error,
              almost certainly within a factor of 10 of the true error
              so long as the next entry is greater than the threshold
              sqrt(n) * dlamch('Epsilon'). This error bound should only
              be trusted if the previous boolean is true.

     err = 3  Reciprocal condition number: Estimated componentwise
              reciprocal condition number.  Compared with the threshold
              sqrt(n) * dlamch('Epsilon') to determine if the error
              estimate is "guaranteed". These reciprocal condition
              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
              appropriately scaled matrix Z.
              Let Z = S*(A*diag(x)), where x is the solution for the
              current right-hand side and S scales each row of
              A*diag(x) by a power of the radix so all absolute row
              sums of Z are approximately 1.

     See Lapack Working Note 165 for further details and extra
     cautions.
[in]NPARAMS
          NPARAMS is INTEGER
     Specifies the number of parameters set in PARAMS.  If .LE. 0, the
     PARAMS array is never referenced and default values are used.
[in,out]PARAMS
          PARAMS is DOUBLE PRECISION array, dimension NPARAMS
     Specifies algorithm parameters.  If an entry is .LT. 0.0, then
     that entry will be filled with default value used for that
     parameter.  Only positions up to NPARAMS are accessed; defaults
     are used for higher-numbered parameters.

       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
            refinement or not.
         Default: 1.0D+0
            = 0.0 : No refinement is performed, and no error bounds are
                    computed.
            = 1.0 : Use the double-precision refinement algorithm,
                    possibly with doubled-single computations if the
                    compilation environment does not support DOUBLE
                    PRECISION.
              (other values are reserved for future use)

       PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
            computations allowed for refinement.
         Default: 10
         Aggressive: Set to 100 to permit convergence using approximate
                     factorizations or factorizations other than LU. If
                     the factorization uses a technique other than
                     Gaussian elimination, the guarantees in
                     err_bnds_norm and err_bnds_comp may no longer be
                     trustworthy.

       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
            will attempt to find a solution with small componentwise
            relative error in the double-precision algorithm.  Positive
            is true, 0.0 is false.
         Default: 1.0 (attempt componentwise convergence)
[out]WORK
          WORK is COMPLEX*16 array, dimension (2*N)
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (2*N)
[out]INFO
          INFO is INTEGER
       = 0:  Successful exit. The solution to every right-hand side is
         guaranteed.
       < 0:  If INFO = -i, the i-th argument had an illegal value
       > 0 and <= N:  U(INFO,INFO) is exactly zero.  The factorization
         has been completed, but the factor U is exactly singular, so
         the solution and error bounds could not be computed. RCOND = 0
         is returned.
       = N+J: The solution corresponding to the Jth right-hand side is
         not guaranteed. The solutions corresponding to other right-
         hand sides K with K > J may not be guaranteed as well, but
         only the first such right-hand side is reported. If a small
         componentwise error is not requested (PARAMS(3) = 0.0) then
         the Jth right-hand side is the first with a normwise error
         bound that is not guaranteed (the smallest J such
         that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
         the Jth right-hand side is the first with either a normwise or
         componentwise error bound that is not guaranteed (the smallest
         J such that either ERR_BNDS_NORM(J,1) = 0.0 or
         ERR_BNDS_COMP(J,1) = 0.0). See the definition of
         ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
         about all of the right-hand sides check ERR_BNDS_NORM or
         ERR_BNDS_COMP.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
April 2012

Definition at line 404 of file zsyrfsx.f.

404 *
405 * -- LAPACK computational routine (version 3.4.1) --
406 * -- LAPACK is a software package provided by Univ. of Tennessee, --
407 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
408 * April 2012
409 *
410 * .. Scalar Arguments ..
411  CHARACTER uplo, equed
412  INTEGER info, lda, ldaf, ldb, ldx, n, nrhs, nparams,
413  $ n_err_bnds
414  DOUBLE PRECISION rcond
415 * ..
416 * .. Array Arguments ..
417  INTEGER ipiv( * )
418  COMPLEX*16 a( lda, * ), af( ldaf, * ), b( ldb, * ),
419  $ x( ldx, * ), work( * )
420  DOUBLE PRECISION s( * ), params( * ), berr( * ), rwork( * ),
421  $ err_bnds_norm( nrhs, * ),
422  $ err_bnds_comp( nrhs, * )
423 * ..
424 *
425 * ==================================================================
426 *
427 * .. Parameters ..
428  DOUBLE PRECISION zero, one
429  parameter( zero = 0.0d+0, one = 1.0d+0 )
430  DOUBLE PRECISION itref_default, ithresh_default
431  DOUBLE PRECISION componentwise_default, rthresh_default
432  DOUBLE PRECISION dzthresh_default
433  parameter( itref_default = 1.0d+0 )
434  parameter( ithresh_default = 10.0d+0 )
435  parameter( componentwise_default = 1.0d+0 )
436  parameter( rthresh_default = 0.5d+0 )
437  parameter( dzthresh_default = 0.25d+0 )
438  INTEGER la_linrx_itref_i, la_linrx_ithresh_i,
439  $ la_linrx_cwise_i
440  parameter( la_linrx_itref_i = 1,
441  $ la_linrx_ithresh_i = 2 )
442  parameter( la_linrx_cwise_i = 3 )
443  INTEGER la_linrx_trust_i, la_linrx_err_i,
444  $ la_linrx_rcond_i
445  parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
446  parameter( la_linrx_rcond_i = 3 )
447 * ..
448 * .. Local Scalars ..
449  CHARACTER(1) norm
450  LOGICAL rcequ
451  INTEGER j, prec_type, ref_type
452  INTEGER n_norms
453  DOUBLE PRECISION anorm, rcond_tmp
454  DOUBLE PRECISION illrcond_thresh, err_lbnd, cwise_wrong
455  LOGICAL ignore_cwise
456  INTEGER ithresh
457  DOUBLE PRECISION rthresh, unstable_thresh
458 * ..
459 * .. External Subroutines ..
461 * ..
462 * .. Intrinsic Functions ..
463  INTRINSIC max, sqrt, transfer
464 * ..
465 * .. External Functions ..
466  EXTERNAL lsame, blas_fpinfo_x, ilatrans, ilaprec
468  DOUBLE PRECISION dlamch, zlansy, zla_syrcond_x, zla_syrcond_c
469  LOGICAL lsame
470  INTEGER blas_fpinfo_x
471  INTEGER ilatrans, ilaprec
472 * ..
473 * .. Executable Statements ..
474 *
475 * Check the input parameters.
476 *
477  info = 0
478  ref_type = int( itref_default )
479  IF ( nparams .GE. la_linrx_itref_i ) THEN
480  IF ( params( la_linrx_itref_i ) .LT. 0.0d+0 ) THEN
481  params( la_linrx_itref_i ) = itref_default
482  ELSE
483  ref_type = params( la_linrx_itref_i )
484  END IF
485  END IF
486 *
487 * Set default parameters.
488 *
489  illrcond_thresh = dble( n ) * dlamch( 'Epsilon' )
490  ithresh = int( ithresh_default )
491  rthresh = rthresh_default
492  unstable_thresh = dzthresh_default
493  ignore_cwise = componentwise_default .EQ. 0.0d+0
494 *
495  IF ( nparams.GE.la_linrx_ithresh_i ) THEN
496  IF ( params( la_linrx_ithresh_i ).LT.0.0d+0 ) THEN
497  params( la_linrx_ithresh_i ) = ithresh
498  ELSE
499  ithresh = int( params( la_linrx_ithresh_i ) )
500  END IF
501  END IF
502  IF ( nparams.GE.la_linrx_cwise_i ) THEN
503  IF ( params( la_linrx_cwise_i ).LT.0.0d+0 ) THEN
504  IF ( ignore_cwise ) THEN
505  params( la_linrx_cwise_i ) = 0.0d+0
506  ELSE
507  params( la_linrx_cwise_i ) = 1.0d+0
508  END IF
509  ELSE
510  ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.0d+0
511  END IF
512  END IF
513  IF ( ref_type .EQ. 0 .OR. n_err_bnds .EQ. 0 ) THEN
514  n_norms = 0
515  ELSE IF ( ignore_cwise ) THEN
516  n_norms = 1
517  ELSE
518  n_norms = 2
519  END IF
520 *
521  rcequ = lsame( equed, 'Y' )
522 *
523 * Test input parameters.
524 *
525  IF ( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
526  info = -1
527  ELSE IF( .NOT.rcequ .AND. .NOT.lsame( equed, 'N' ) ) THEN
528  info = -2
529  ELSE IF( n.LT.0 ) THEN
530  info = -3
531  ELSE IF( nrhs.LT.0 ) THEN
532  info = -4
533  ELSE IF( lda.LT.max( 1, n ) ) THEN
534  info = -6
535  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
536  info = -8
537  ELSE IF( ldb.LT.max( 1, n ) ) THEN
538  info = -12
539  ELSE IF( ldx.LT.max( 1, n ) ) THEN
540  info = -14
541  END IF
542  IF( info.NE.0 ) THEN
543  CALL xerbla( 'ZSYRFSX', -info )
544  RETURN
545  END IF
546 *
547 * Quick return if possible.
548 *
549  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
550  rcond = 1.0d+0
551  DO j = 1, nrhs
552  berr( j ) = 0.0d+0
553  IF ( n_err_bnds .GE. 1 ) THEN
554  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
555  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
556  END IF
557  IF ( n_err_bnds .GE. 2 ) THEN
558  err_bnds_norm( j, la_linrx_err_i ) = 0.0d+0
559  err_bnds_comp( j, la_linrx_err_i ) = 0.0d+0
560  END IF
561  IF ( n_err_bnds .GE. 3 ) THEN
562  err_bnds_norm( j, la_linrx_rcond_i ) = 1.0d+0
563  err_bnds_comp( j, la_linrx_rcond_i ) = 1.0d+0
564  END IF
565  END DO
566  RETURN
567  END IF
568 *
569 * Default to failure.
570 *
571  rcond = 0.0d+0
572  DO j = 1, nrhs
573  berr( j ) = 1.0d+0
574  IF ( n_err_bnds .GE. 1 ) THEN
575  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
576  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
577  END IF
578  IF ( n_err_bnds .GE. 2 ) THEN
579  err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
580  err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
581  END IF
582  IF ( n_err_bnds .GE. 3 ) THEN
583  err_bnds_norm( j, la_linrx_rcond_i ) = 0.0d+0
584  err_bnds_comp( j, la_linrx_rcond_i ) = 0.0d+0
585  END IF
586  END DO
587 *
588 * Compute the norm of A and the reciprocal of the condition
589 * number of A.
590 *
591  norm = 'I'
592  anorm = zlansy( norm, uplo, n, a, lda, rwork )
593  CALL zsycon( uplo, n, af, ldaf, ipiv, anorm, rcond, work,
594  $ info )
595 *
596 * Perform refinement on each right-hand side
597 *
598  IF ( ref_type .NE. 0 ) THEN
599 
600  prec_type = ilaprec( 'E' )
601 
602  CALL zla_syrfsx_extended( prec_type, uplo, n,
603  $ nrhs, a, lda, af, ldaf, ipiv, rcequ, s, b,
604  $ ldb, x, ldx, berr, n_norms, err_bnds_norm, err_bnds_comp,
605  $ work, rwork, work(n+1),
606  $ transfer(rwork(1:2*n), (/ (zero, zero) /), n), rcond,
607  $ ithresh, rthresh, unstable_thresh, ignore_cwise,
608  $ info )
609  END IF
610 
611  err_lbnd = max( 10.0d+0, sqrt( dble( n ) ) ) * dlamch( 'Epsilon' )
612  IF (n_err_bnds .GE. 1 .AND. n_norms .GE. 1) THEN
613 *
614 * Compute scaled normwise condition number cond(A*C).
615 *
616  IF ( rcequ ) THEN
617  rcond_tmp = zla_syrcond_c( uplo, n, a, lda, af, ldaf, ipiv,
618  $ s, .true., info, work, rwork )
619  ELSE
620  rcond_tmp = zla_syrcond_c( uplo, n, a, lda, af, ldaf, ipiv,
621  $ s, .false., info, work, rwork )
622  END IF
623  DO j = 1, nrhs
624 *
625 * Cap the error at 1.0.
626 *
627  IF ( n_err_bnds .GE. la_linrx_err_i
628  $ .AND. err_bnds_norm( j, la_linrx_err_i ) .GT. 1.0d+0 )
629  $ err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
630 *
631 * Threshold the error (see LAWN).
632 *
633  IF ( rcond_tmp .LT. illrcond_thresh ) THEN
634  err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
635  err_bnds_norm( j, la_linrx_trust_i ) = 0.0d+0
636  IF ( info .LE. n ) info = n + j
637  ELSE IF ( err_bnds_norm( j, la_linrx_err_i ) .LT. err_lbnd )
638  $ THEN
639  err_bnds_norm( j, la_linrx_err_i ) = err_lbnd
640  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
641  END IF
642 *
643 * Save the condition number.
644 *
645  IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
646  err_bnds_norm( j, la_linrx_rcond_i ) = rcond_tmp
647  END IF
648  END DO
649  END IF
650 
651  IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 2 ) THEN
652 *
653 * Compute componentwise condition number cond(A*diag(Y(:,J))) for
654 * each right-hand side using the current solution as an estimate of
655 * the true solution. If the componentwise error estimate is too
656 * large, then the solution is a lousy estimate of truth and the
657 * estimated RCOND may be too optimistic. To avoid misleading users,
658 * the inverse condition number is set to 0.0 when the estimated
659 * cwise error is at least CWISE_WRONG.
660 *
661  cwise_wrong = sqrt( dlamch( 'Epsilon' ) )
662  DO j = 1, nrhs
663  IF ( err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
664  $ THEN
665  rcond_tmp = zla_syrcond_x( uplo, n, a, lda, af, ldaf,
666  $ ipiv, x(1,j), info, work, rwork )
667  ELSE
668  rcond_tmp = 0.0d+0
669  END IF
670 *
671 * Cap the error at 1.0.
672 *
673  IF ( n_err_bnds .GE. la_linrx_err_i
674  $ .AND. err_bnds_comp( j, la_linrx_err_i ) .GT. 1.0d+0 )
675  $ err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
676 
677 *
678 * Threshold the error (see LAWN).
679 *
680  IF ( rcond_tmp .LT. illrcond_thresh ) THEN
681  err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
682  err_bnds_comp( j, la_linrx_trust_i ) = 0.0d+0
683  IF (.NOT. ignore_cwise
684  $ .AND. info.LT.n + j ) info = n + j
685  ELSE IF ( err_bnds_comp( j, la_linrx_err_i )
686  $ .LT. err_lbnd ) THEN
687  err_bnds_comp( j, la_linrx_err_i ) = err_lbnd
688  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
689  END IF
690 *
691 * Save the condition number.
692 *
693  IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
694  err_bnds_comp( j, la_linrx_rcond_i ) = rcond_tmp
695  END IF
696 
697  END DO
698  END IF
699 *
700  RETURN
701 *
702 * End of ZSYRFSX
703 *
integer function ilatrans(TRANS)
ILATRANS
Definition: ilatrans.f:60
double precision function zlansy(NORM, UPLO, N, A, LDA, WORK)
ZLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex symmetric matrix.
Definition: zlansy.f:125
double precision function zla_syrcond_x(UPLO, N, A, LDA, AF, LDAF, IPIV, X, INFO, WORK, RWORK)
ZLA_SYRCOND_X computes the infinity norm condition number of op(A)*diag(x) for symmetric indefinite m...
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
integer function ilaprec(PREC)
ILAPREC
Definition: ilaprec.f:60
subroutine zla_syrfsx_extended(PREC_TYPE, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, COLEQU, C, B, LDB, Y, LDY, BERR_OUT, N_NORMS, ERR_BNDS_NORM, ERR_BNDS_COMP, RES, AYB, DY, Y_TAIL, RCOND, ITHRESH, RTHRESH, DZ_UB, IGNORE_CWISE, INFO)
ZLA_SYRFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric inde...
subroutine zsycon(UPLO, N, A, LDA, IPIV, ANORM, RCOND, WORK, INFO)
ZSYCON
Definition: zsycon.f:127
double precision function zla_syrcond_c(UPLO, N, A, LDA, AF, LDAF, IPIV, C, CAPPLY, INFO, WORK, RWORK)
ZLA_SYRCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for symmetric indefin...

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zsytf2 ( character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
integer  INFO 
)

ZSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting method (unblocked algorithm).

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

Purpose:
 ZSYTF2 computes the factorization of a complex symmetric matrix A
 using the Bunch-Kaufman diagonal pivoting method:

    A = U*D*U**T  or  A = L*D*L**T

 where U (or L) is a product of permutation and unit upper (lower)
 triangular matrices, U**T is the transpose of U, and D is symmetric and
 block diagonal with 1-by-1 and 2-by-2 diagonal blocks.

 This is the unblocked version of the algorithm, calling Level 2 BLAS.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the upper or lower triangular part of the
          symmetric matrix A is stored:
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[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 symmetric matrix A.  If UPLO = 'U', the leading
          n-by-n upper triangular part of A contains the upper
          triangular part of the matrix A, and the strictly lower
          triangular part of A is not referenced.  If UPLO = 'L', the
          leading n-by-n lower triangular part of A contains the lower
          triangular part of the matrix A, and the strictly upper
          triangular part of A is not referenced.

          On exit, the block diagonal matrix D and the multipliers used
          to obtain the factor U or L (see below for further details).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D.

          If UPLO = 'U':
             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
             interchanged and D(k,k) is a 1-by-1 diagonal block.

             If IPIV(k) = IPIV(k-1) < 0, then rows and columns
             k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
             is a 2-by-2 diagonal block.

          If UPLO = 'L':
             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
             interchanged and D(k,k) is a 1-by-1 diagonal block.

             If IPIV(k) = IPIV(k+1) < 0, then rows and columns
             k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1)
             is a 2-by-2 diagonal block.
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -k, the k-th argument had an illegal value
          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
               has been completed, but the block diagonal matrix D is
               exactly singular, and division by zero will occur if it
               is used to solve a system of equations.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2013
Further Details:
  If UPLO = 'U', then A = U*D*U**T, where
     U = P(n)*U(n)* ... *P(k)U(k)* ...,
  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
  that if the diagonal block D(k) is of order s (s = 1 or 2), then

             (   I    v    0   )   k-s
     U(k) =  (   0    I    0   )   s
             (   0    0    I   )   n-k
                k-s   s   n-k

  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
  and A(k,k), and v overwrites A(1:k-2,k-1:k).

  If UPLO = 'L', then A = L*D*L**T, where
     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
  that if the diagonal block D(k) is of order s (s = 1 or 2), then

             (   I    0     0   )  k-1
     L(k) =  (   0    I     0   )  s
             (   0    v     I   )  n-k-s+1
                k-1   s  n-k-s+1

  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
Contributors:
  09-29-06 - patch from
    Bobby Cheng, MathWorks

    Replace l.209 and l.377
         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
    by
         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN

  1-96 - Based on modifications by J. Lewis, Boeing Computer Services
         Company

Definition at line 193 of file zsytf2.f.

193 *
194 * -- LAPACK computational routine (version 3.5.0) --
195 * -- LAPACK is a software package provided by Univ. of Tennessee, --
196 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
197 * November 2013
198 *
199 * .. Scalar Arguments ..
200  CHARACTER uplo
201  INTEGER info, lda, n
202 * ..
203 * .. Array Arguments ..
204  INTEGER ipiv( * )
205  COMPLEX*16 a( lda, * )
206 * ..
207 *
208 * =====================================================================
209 *
210 * .. Parameters ..
211  DOUBLE PRECISION zero, one
212  parameter( zero = 0.0d+0, one = 1.0d+0 )
213  DOUBLE PRECISION eight, sevten
214  parameter( eight = 8.0d+0, sevten = 17.0d+0 )
215  COMPLEX*16 cone
216  parameter( cone = ( 1.0d+0, 0.0d+0 ) )
217 * ..
218 * .. Local Scalars ..
219  LOGICAL upper
220  INTEGER i, imax, j, jmax, k, kk, kp, kstep
221  DOUBLE PRECISION absakk, alpha, colmax, rowmax
222  COMPLEX*16 d11, d12, d21, d22, r1, t, wk, wkm1, wkp1, z
223 * ..
224 * .. External Functions ..
225  LOGICAL disnan, lsame
226  INTEGER izamax
227  EXTERNAL disnan, lsame, izamax
228 * ..
229 * .. External Subroutines ..
230  EXTERNAL xerbla, zscal, zswap, zsyr
231 * ..
232 * .. Intrinsic Functions ..
233  INTRINSIC abs, dble, dimag, max, sqrt
234 * ..
235 * .. Statement Functions ..
236  DOUBLE PRECISION cabs1
237 * ..
238 * .. Statement Function definitions ..
239  cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
240 * ..
241 * .. Executable Statements ..
242 *
243 * Test the input parameters.
244 *
245  info = 0
246  upper = lsame( uplo, 'U' )
247  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
248  info = -1
249  ELSE IF( n.LT.0 ) THEN
250  info = -2
251  ELSE IF( lda.LT.max( 1, n ) ) THEN
252  info = -4
253  END IF
254  IF( info.NE.0 ) THEN
255  CALL xerbla( 'ZSYTF2', -info )
256  RETURN
257  END IF
258 *
259 * Initialize ALPHA for use in choosing pivot block size.
260 *
261  alpha = ( one+sqrt( sevten ) ) / eight
262 *
263  IF( upper ) THEN
264 *
265 * Factorize A as U*D*U**T using the upper triangle of A
266 *
267 * K is the main loop index, decreasing from N to 1 in steps of
268 * 1 or 2
269 *
270  k = n
271  10 CONTINUE
272 *
273 * If K < 1, exit from loop
274 *
275  IF( k.LT.1 )
276  $ GO TO 70
277  kstep = 1
278 *
279 * Determine rows and columns to be interchanged and whether
280 * a 1-by-1 or 2-by-2 pivot block will be used
281 *
282  absakk = cabs1( a( k, k ) )
283 *
284 * IMAX is the row-index of the largest off-diagonal element in
285 * column K, and COLMAX is its absolute value.
286 * Determine both COLMAX and IMAX.
287 *
288  IF( k.GT.1 ) THEN
289  imax = izamax( k-1, a( 1, k ), 1 )
290  colmax = cabs1( a( imax, k ) )
291  ELSE
292  colmax = zero
293  END IF
294 *
295  IF( max( absakk, colmax ).EQ.zero .OR. disnan(absakk) ) THEN
296 *
297 * Column K is zero or underflow, or contains a NaN:
298 * set INFO and continue
299 *
300  IF( info.EQ.0 )
301  $ info = k
302  kp = k
303  ELSE
304  IF( absakk.GE.alpha*colmax ) THEN
305 *
306 * no interchange, use 1-by-1 pivot block
307 *
308  kp = k
309  ELSE
310 *
311 * JMAX is the column-index of the largest off-diagonal
312 * element in row IMAX, and ROWMAX is its absolute value
313 *
314  jmax = imax + izamax( k-imax, a( imax, imax+1 ), lda )
315  rowmax = cabs1( a( imax, jmax ) )
316  IF( imax.GT.1 ) THEN
317  jmax = izamax( imax-1, a( 1, imax ), 1 )
318  rowmax = max( rowmax, cabs1( a( jmax, imax ) ) )
319  END IF
320 *
321  IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
322 *
323 * no interchange, use 1-by-1 pivot block
324 *
325  kp = k
326  ELSE IF( cabs1( a( imax, imax ) ).GE.alpha*rowmax ) THEN
327 *
328 * interchange rows and columns K and IMAX, use 1-by-1
329 * pivot block
330 *
331  kp = imax
332  ELSE
333 *
334 * interchange rows and columns K-1 and IMAX, use 2-by-2
335 * pivot block
336 *
337  kp = imax
338  kstep = 2
339  END IF
340  END IF
341 *
342  kk = k - kstep + 1
343  IF( kp.NE.kk ) THEN
344 *
345 * Interchange rows and columns KK and KP in the leading
346 * submatrix A(1:k,1:k)
347 *
348  CALL zswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
349  CALL zswap( kk-kp-1, a( kp+1, kk ), 1, a( kp, kp+1 ),
350  $ lda )
351  t = a( kk, kk )
352  a( kk, kk ) = a( kp, kp )
353  a( kp, kp ) = t
354  IF( kstep.EQ.2 ) THEN
355  t = a( k-1, k )
356  a( k-1, k ) = a( kp, k )
357  a( kp, k ) = t
358  END IF
359  END IF
360 *
361 * Update the leading submatrix
362 *
363  IF( kstep.EQ.1 ) THEN
364 *
365 * 1-by-1 pivot block D(k): column k now holds
366 *
367 * W(k) = U(k)*D(k)
368 *
369 * where U(k) is the k-th column of U
370 *
371 * Perform a rank-1 update of A(1:k-1,1:k-1) as
372 *
373 * A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T
374 *
375  r1 = cone / a( k, k )
376  CALL zsyr( uplo, k-1, -r1, a( 1, k ), 1, a, lda )
377 *
378 * Store U(k) in column k
379 *
380  CALL zscal( k-1, r1, a( 1, k ), 1 )
381  ELSE
382 *
383 * 2-by-2 pivot block D(k): columns k and k-1 now hold
384 *
385 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
386 *
387 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
388 * of U
389 *
390 * Perform a rank-2 update of A(1:k-2,1:k-2) as
391 *
392 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
393 * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**T
394 *
395  IF( k.GT.2 ) THEN
396 *
397  d12 = a( k-1, k )
398  d22 = a( k-1, k-1 ) / d12
399  d11 = a( k, k ) / d12
400  t = cone / ( d11*d22-cone )
401  d12 = t / d12
402 *
403  DO 30 j = k - 2, 1, -1
404  wkm1 = d12*( d11*a( j, k-1 )-a( j, k ) )
405  wk = d12*( d22*a( j, k )-a( j, k-1 ) )
406  DO 20 i = j, 1, -1
407  a( i, j ) = a( i, j ) - a( i, k )*wk -
408  $ a( i, k-1 )*wkm1
409  20 CONTINUE
410  a( j, k ) = wk
411  a( j, k-1 ) = wkm1
412  30 CONTINUE
413 *
414  END IF
415 *
416  END IF
417  END IF
418 *
419 * Store details of the interchanges in IPIV
420 *
421  IF( kstep.EQ.1 ) THEN
422  ipiv( k ) = kp
423  ELSE
424  ipiv( k ) = -kp
425  ipiv( k-1 ) = -kp
426  END IF
427 *
428 * Decrease K and return to the start of the main loop
429 *
430  k = k - kstep
431  GO TO 10
432 *
433  ELSE
434 *
435 * Factorize A as L*D*L**T using the lower triangle of A
436 *
437 * K is the main loop index, increasing from 1 to N in steps of
438 * 1 or 2
439 *
440  k = 1
441  40 CONTINUE
442 *
443 * If K > N, exit from loop
444 *
445  IF( k.GT.n )
446  $ GO TO 70
447  kstep = 1
448 *
449 * Determine rows and columns to be interchanged and whether
450 * a 1-by-1 or 2-by-2 pivot block will be used
451 *
452  absakk = cabs1( a( k, k ) )
453 *
454 * IMAX is the row-index of the largest off-diagonal element in
455 * column K, and COLMAX is its absolute value.
456 * Determine both COLMAX and IMAX.
457 *
458  IF( k.LT.n ) THEN
459  imax = k + izamax( n-k, a( k+1, k ), 1 )
460  colmax = cabs1( a( imax, k ) )
461  ELSE
462  colmax = zero
463  END IF
464 *
465  IF( max( absakk, colmax ).EQ.zero .OR. disnan(absakk) ) THEN
466 *
467 * Column K is zero or underflow, or contains a NaN:
468 * set INFO and continue
469 *
470  IF( info.EQ.0 )
471  $ info = k
472  kp = k
473  ELSE
474  IF( absakk.GE.alpha*colmax ) THEN
475 *
476 * no interchange, use 1-by-1 pivot block
477 *
478  kp = k
479  ELSE
480 *
481 * JMAX is the column-index of the largest off-diagonal
482 * element in row IMAX, and ROWMAX is its absolute value
483 *
484  jmax = k - 1 + izamax( imax-k, a( imax, k ), lda )
485  rowmax = cabs1( a( imax, jmax ) )
486  IF( imax.LT.n ) THEN
487  jmax = imax + izamax( n-imax, a( imax+1, imax ), 1 )
488  rowmax = max( rowmax, cabs1( a( jmax, imax ) ) )
489  END IF
490 *
491  IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
492 *
493 * no interchange, use 1-by-1 pivot block
494 *
495  kp = k
496  ELSE IF( cabs1( a( imax, imax ) ).GE.alpha*rowmax ) THEN
497 *
498 * interchange rows and columns K and IMAX, use 1-by-1
499 * pivot block
500 *
501  kp = imax
502  ELSE
503 *
504 * interchange rows and columns K+1 and IMAX, use 2-by-2
505 * pivot block
506 *
507  kp = imax
508  kstep = 2
509  END IF
510  END IF
511 *
512  kk = k + kstep - 1
513  IF( kp.NE.kk ) THEN
514 *
515 * Interchange rows and columns KK and KP in the trailing
516 * submatrix A(k:n,k:n)
517 *
518  IF( kp.LT.n )
519  $ CALL zswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
520  CALL zswap( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
521  $ lda )
522  t = a( kk, kk )
523  a( kk, kk ) = a( kp, kp )
524  a( kp, kp ) = t
525  IF( kstep.EQ.2 ) THEN
526  t = a( k+1, k )
527  a( k+1, k ) = a( kp, k )
528  a( kp, k ) = t
529  END IF
530  END IF
531 *
532 * Update the trailing submatrix
533 *
534  IF( kstep.EQ.1 ) THEN
535 *
536 * 1-by-1 pivot block D(k): column k now holds
537 *
538 * W(k) = L(k)*D(k)
539 *
540 * where L(k) is the k-th column of L
541 *
542  IF( k.LT.n ) THEN
543 *
544 * Perform a rank-1 update of A(k+1:n,k+1:n) as
545 *
546 * A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T
547 *
548  r1 = cone / a( k, k )
549  CALL zsyr( uplo, n-k, -r1, a( k+1, k ), 1,
550  $ a( k+1, k+1 ), lda )
551 *
552 * Store L(k) in column K
553 *
554  CALL zscal( n-k, r1, a( k+1, k ), 1 )
555  END IF
556  ELSE
557 *
558 * 2-by-2 pivot block D(k)
559 *
560  IF( k.LT.n-1 ) THEN
561 *
562 * Perform a rank-2 update of A(k+2:n,k+2:n) as
563 *
564 * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**T
565 * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**T
566 *
567 * where L(k) and L(k+1) are the k-th and (k+1)-th
568 * columns of L
569 *
570  d21 = a( k+1, k )
571  d11 = a( k+1, k+1 ) / d21
572  d22 = a( k, k ) / d21
573  t = cone / ( d11*d22-cone )
574  d21 = t / d21
575 *
576  DO 60 j = k + 2, n
577  wk = d21*( d11*a( j, k )-a( j, k+1 ) )
578  wkp1 = d21*( d22*a( j, k+1 )-a( j, k ) )
579  DO 50 i = j, n
580  a( i, j ) = a( i, j ) - a( i, k )*wk -
581  $ a( i, k+1 )*wkp1
582  50 CONTINUE
583  a( j, k ) = wk
584  a( j, k+1 ) = wkp1
585  60 CONTINUE
586  END IF
587  END IF
588  END IF
589 *
590 * Store details of the interchanges in IPIV
591 *
592  IF( kstep.EQ.1 ) THEN
593  ipiv( k ) = kp
594  ELSE
595  ipiv( k ) = -kp
596  ipiv( k+1 ) = -kp
597  END IF
598 *
599 * Increase K and return to the start of the main loop
600 *
601  k = k + kstep
602  GO TO 40
603 *
604  END IF
605 *
606  70 CONTINUE
607  RETURN
608 *
609 * End of ZSYTF2
610 *
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:61
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zsyr(UPLO, N, ALPHA, X, INCX, A, LDA)
ZSYR performs the symmetric rank-1 update of a complex symmetric matrix.
Definition: zsyr.f:137
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:54
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zsytf2_rook ( character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
integer  INFO 
)

ZSYTF2_ROOK computes the factorization of a complex symmetric indefinite matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method (unblocked algorithm).

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

Purpose:
 ZSYTF2_ROOK computes the factorization of a complex symmetric matrix A
 using the bounded Bunch-Kaufman ("rook") diagonal pivoting method:

    A = U*D*U**T  or  A = L*D*L**T

 where U (or L) is a product of permutation and unit upper (lower)
 triangular matrices, U**T is the transpose of U, and D is symmetric and
 block diagonal with 1-by-1 and 2-by-2 diagonal blocks.

 This is the unblocked version of the algorithm, calling Level 2 BLAS.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the upper or lower triangular part of the
          symmetric matrix A is stored:
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[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 symmetric matrix A.  If UPLO = 'U', the leading
          n-by-n upper triangular part of A contains the upper
          triangular part of the matrix A, and the strictly lower
          triangular part of A is not referenced.  If UPLO = 'L', the
          leading n-by-n lower triangular part of A contains the lower
          triangular part of the matrix A, and the strictly upper
          triangular part of A is not referenced.

          On exit, the block diagonal matrix D and the multipliers used
          to obtain the factor U or L (see below for further details).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D.

          If UPLO = 'U':
             If IPIV(k) > 0, then rows and columns k and IPIV(k)
             were interchanged and D(k,k) is a 1-by-1 diagonal block.

             If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
             columns k and -IPIV(k) were interchanged and rows and
             columns k-1 and -IPIV(k-1) were inerchaged,
             D(k-1:k,k-1:k) is a 2-by-2 diagonal block.

          If UPLO = 'L':
             If IPIV(k) > 0, then rows and columns k and IPIV(k)
             were interchanged and D(k,k) is a 1-by-1 diagonal block.

             If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
             columns k and -IPIV(k) were interchanged and rows and
             columns k+1 and -IPIV(k+1) were inerchaged,
             D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -k, the k-th argument had an illegal value
          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
               has been completed, but the block diagonal matrix D is
               exactly singular, and division by zero will occur if it
               is used to solve a system of equations.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2013
Further Details:
  If UPLO = 'U', then A = U*D*U**T, where
     U = P(n)*U(n)* ... *P(k)U(k)* ...,
  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
  that if the diagonal block D(k) is of order s (s = 1 or 2), then

             (   I    v    0   )   k-s
     U(k) =  (   0    I    0   )   s
             (   0    0    I   )   n-k
                k-s   s   n-k

  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
  and A(k,k), and v overwrites A(1:k-2,k-1:k).

  If UPLO = 'L', then A = L*D*L**T, where
     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
  that if the diagonal block D(k) is of order s (s = 1 or 2), then

             (   I    0     0   )  k-1
     L(k) =  (   0    I     0   )  s
             (   0    v     I   )  n-k-s+1
                k-1   s  n-k-s+1

  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
Contributors:
  November 2013,     Igor Kozachenko,
                  Computer Science Division,
                  University of California, Berkeley

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

  01-01-96 - Based on modifications by
    J. Lewis, Boeing Computer Services Company
    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville abd , USA

Definition at line 196 of file zsytf2_rook.f.

196 *
197 * -- LAPACK computational routine (version 3.5.0) --
198 * -- LAPACK is a software package provided by Univ. of Tennessee, --
199 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
200 * November 2013
201 *
202 * .. Scalar Arguments ..
203  CHARACTER uplo
204  INTEGER info, lda, n
205 * ..
206 * .. Array Arguments ..
207  INTEGER ipiv( * )
208  COMPLEX*16 a( lda, * )
209 * ..
210 *
211 * =====================================================================
212 *
213 * .. Parameters ..
214  DOUBLE PRECISION zero, one
215  parameter( zero = 0.0d+0, one = 1.0d+0 )
216  DOUBLE PRECISION eight, sevten
217  parameter( eight = 8.0d+0, sevten = 17.0d+0 )
218  COMPLEX*16 cone
219  parameter( cone = ( 1.0d+0, 0.0d+0 ) )
220 * ..
221 * .. Local Scalars ..
222  LOGICAL upper, done
223  INTEGER i, imax, j, jmax, itemp, k, kk, kp, kstep,
224  $ p, ii
225  DOUBLE PRECISION absakk, alpha, colmax, rowmax, dtemp, sfmin
226  COMPLEX*16 d11, d12, d21, d22, t, wk, wkm1, wkp1, z
227 * ..
228 * .. External Functions ..
229  LOGICAL lsame
230  INTEGER izamax
231  DOUBLE PRECISION dlamch
232  EXTERNAL lsame, izamax, dlamch
233 * ..
234 * .. External Subroutines ..
235  EXTERNAL zscal, zswap, zsyr, xerbla
236 * ..
237 * .. Intrinsic Functions ..
238  INTRINSIC abs, max, sqrt, dimag, dble
239 * ..
240 * .. Statement Functions ..
241  DOUBLE PRECISION cabs1
242 * ..
243 * .. Statement Function definitions ..
244  cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
245 * ..
246 * .. Executable Statements ..
247 *
248 * Test the input parameters.
249 *
250  info = 0
251  upper = lsame( uplo, 'U' )
252  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
253  info = -1
254  ELSE IF( n.LT.0 ) THEN
255  info = -2
256  ELSE IF( lda.LT.max( 1, n ) ) THEN
257  info = -4
258  END IF
259  IF( info.NE.0 ) THEN
260  CALL xerbla( 'ZSYTF2_ROOK', -info )
261  RETURN
262  END IF
263 *
264 * Initialize ALPHA for use in choosing pivot block size.
265 *
266  alpha = ( one+sqrt( sevten ) ) / eight
267 *
268 * Compute machine safe minimum
269 *
270  sfmin = dlamch( 'S' )
271 *
272  IF( upper ) THEN
273 *
274 * Factorize A as U*D*U**T using the upper triangle of A
275 *
276 * K is the main loop index, decreasing from N to 1 in steps of
277 * 1 or 2
278 *
279  k = n
280  10 CONTINUE
281 *
282 * If K < 1, exit from loop
283 *
284  IF( k.LT.1 )
285  $ GO TO 70
286  kstep = 1
287  p = k
288 *
289 * Determine rows and columns to be interchanged and whether
290 * a 1-by-1 or 2-by-2 pivot block will be used
291 *
292  absakk = cabs1( a( k, k ) )
293 *
294 * IMAX is the row-index of the largest off-diagonal element in
295 * column K, and COLMAX is its absolute value.
296 * Determine both COLMAX and IMAX.
297 *
298  IF( k.GT.1 ) THEN
299  imax = izamax( k-1, a( 1, k ), 1 )
300  colmax = cabs1( a( imax, k ) )
301  ELSE
302  colmax = zero
303  END IF
304 *
305  IF( (max( absakk, colmax ).EQ.zero) ) THEN
306 *
307 * Column K is zero or underflow: set INFO and continue
308 *
309  IF( info.EQ.0 )
310  $ info = k
311  kp = k
312  ELSE
313 *
314 * Test for interchange
315 *
316 * Equivalent to testing for (used to handle NaN and Inf)
317 * ABSAKK.GE.ALPHA*COLMAX
318 *
319  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
320 *
321 * no interchange,
322 * use 1-by-1 pivot block
323 *
324  kp = k
325  ELSE
326 *
327  done = .false.
328 *
329 * Loop until pivot found
330 *
331  12 CONTINUE
332 *
333 * Begin pivot search loop body
334 *
335 * JMAX is the column-index of the largest off-diagonal
336 * element in row IMAX, and ROWMAX is its absolute value.
337 * Determine both ROWMAX and JMAX.
338 *
339  IF( imax.NE.k ) THEN
340  jmax = imax + izamax( k-imax, a( imax, imax+1 ),
341  $ lda )
342  rowmax = cabs1( a( imax, jmax ) )
343  ELSE
344  rowmax = zero
345  END IF
346 *
347  IF( imax.GT.1 ) THEN
348  itemp = izamax( imax-1, a( 1, imax ), 1 )
349  dtemp = cabs1( a( itemp, imax ) )
350  IF( dtemp.GT.rowmax ) THEN
351  rowmax = dtemp
352  jmax = itemp
353  END IF
354  END IF
355 *
356 * Equivalent to testing for (used to handle NaN and Inf)
357 * CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
358 *
359  IF( .NOT.( cabs1(a( imax, imax )).LT.alpha*rowmax ) )
360  $ THEN
361 *
362 * interchange rows and columns K and IMAX,
363 * use 1-by-1 pivot block
364 *
365  kp = imax
366  done = .true.
367 *
368 * Equivalent to testing for ROWMAX .EQ. COLMAX,
369 * used to handle NaN and Inf
370 *
371  ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
372 *
373 * interchange rows and columns K+1 and IMAX,
374 * use 2-by-2 pivot block
375 *
376  kp = imax
377  kstep = 2
378  done = .true.
379  ELSE
380 *
381 * Pivot NOT found, set variables and repeat
382 *
383  p = imax
384  colmax = rowmax
385  imax = jmax
386  END IF
387 *
388 * End pivot search loop body
389 *
390  IF( .NOT. done ) GOTO 12
391 *
392  END IF
393 *
394 * Swap TWO rows and TWO columns
395 *
396 * First swap
397 *
398  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
399 *
400 * Interchange rows and column K and P in the leading
401 * submatrix A(1:k,1:k) if we have a 2-by-2 pivot
402 *
403  IF( p.GT.1 )
404  $ CALL zswap( p-1, a( 1, k ), 1, a( 1, p ), 1 )
405  IF( p.LT.(k-1) )
406  $ CALL zswap( k-p-1, a( p+1, k ), 1, a( p, p+1 ),
407  $ lda )
408  t = a( k, k )
409  a( k, k ) = a( p, p )
410  a( p, p ) = t
411  END IF
412 *
413 * Second swap
414 *
415  kk = k - kstep + 1
416  IF( kp.NE.kk ) THEN
417 *
418 * Interchange rows and columns KK and KP in the leading
419 * submatrix A(1:k,1:k)
420 *
421  IF( kp.GT.1 )
422  $ CALL zswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
423  IF( ( kk.GT.1 ) .AND. ( kp.LT.(kk-1) ) )
424  $ CALL zswap( kk-kp-1, a( kp+1, kk ), 1, a( kp, kp+1 ),
425  $ lda )
426  t = a( kk, kk )
427  a( kk, kk ) = a( kp, kp )
428  a( kp, kp ) = t
429  IF( kstep.EQ.2 ) THEN
430  t = a( k-1, k )
431  a( k-1, k ) = a( kp, k )
432  a( kp, k ) = t
433  END IF
434  END IF
435 *
436 * Update the leading submatrix
437 *
438  IF( kstep.EQ.1 ) THEN
439 *
440 * 1-by-1 pivot block D(k): column k now holds
441 *
442 * W(k) = U(k)*D(k)
443 *
444 * where U(k) is the k-th column of U
445 *
446  IF( k.GT.1 ) THEN
447 *
448 * Perform a rank-1 update of A(1:k-1,1:k-1) and
449 * store U(k) in column k
450 *
451  IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
452 *
453 * Perform a rank-1 update of A(1:k-1,1:k-1) as
454 * A := A - U(k)*D(k)*U(k)**T
455 * = A - W(k)*1/D(k)*W(k)**T
456 *
457  d11 = cone / a( k, k )
458  CALL zsyr( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
459 *
460 * Store U(k) in column k
461 *
462  CALL zscal( k-1, d11, a( 1, k ), 1 )
463  ELSE
464 *
465 * Store L(k) in column K
466 *
467  d11 = a( k, k )
468  DO 16 ii = 1, k - 1
469  a( ii, k ) = a( ii, k ) / d11
470  16 CONTINUE
471 *
472 * Perform a rank-1 update of A(k+1:n,k+1:n) as
473 * A := A - U(k)*D(k)*U(k)**T
474 * = A - W(k)*(1/D(k))*W(k)**T
475 * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
476 *
477  CALL zsyr( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
478  END IF
479  END IF
480 *
481  ELSE
482 *
483 * 2-by-2 pivot block D(k): columns k and k-1 now hold
484 *
485 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
486 *
487 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
488 * of U
489 *
490 * Perform a rank-2 update of A(1:k-2,1:k-2) as
491 *
492 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
493 * = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
494 *
495 * and store L(k) and L(k+1) in columns k and k+1
496 *
497  IF( k.GT.2 ) THEN
498 *
499  d12 = a( k-1, k )
500  d22 = a( k-1, k-1 ) / d12
501  d11 = a( k, k ) / d12
502  t = cone / ( d11*d22-cone )
503 *
504  DO 30 j = k - 2, 1, -1
505 *
506  wkm1 = t*( d11*a( j, k-1 )-a( j, k ) )
507  wk = t*( d22*a( j, k )-a( j, k-1 ) )
508 *
509  DO 20 i = j, 1, -1
510  a( i, j ) = a( i, j ) - (a( i, k ) / d12 )*wk -
511  $ ( a( i, k-1 ) / d12 )*wkm1
512  20 CONTINUE
513 *
514 * Store U(k) and U(k-1) in cols k and k-1 for row J
515 *
516  a( j, k ) = wk / d12
517  a( j, k-1 ) = wkm1 / d12
518 *
519  30 CONTINUE
520 *
521  END IF
522 *
523  END IF
524  END IF
525 *
526 * Store details of the interchanges in IPIV
527 *
528  IF( kstep.EQ.1 ) THEN
529  ipiv( k ) = kp
530  ELSE
531  ipiv( k ) = -p
532  ipiv( k-1 ) = -kp
533  END IF
534 *
535 * Decrease K and return to the start of the main loop
536 *
537  k = k - kstep
538  GO TO 10
539 *
540  ELSE
541 *
542 * Factorize A as L*D*L**T using the lower triangle of A
543 *
544 * K is the main loop index, increasing from 1 to N in steps of
545 * 1 or 2
546 *
547  k = 1
548  40 CONTINUE
549 *
550 * If K > N, exit from loop
551 *
552  IF( k.GT.n )
553  $ GO TO 70
554  kstep = 1
555  p = k
556 *
557 * Determine rows and columns to be interchanged and whether
558 * a 1-by-1 or 2-by-2 pivot block will be used
559 *
560  absakk = cabs1( a( k, k ) )
561 *
562 * IMAX is the row-index of the largest off-diagonal element in
563 * column K, and COLMAX is its absolute value.
564 * Determine both COLMAX and IMAX.
565 *
566  IF( k.LT.n ) THEN
567  imax = k + izamax( n-k, a( k+1, k ), 1 )
568  colmax = cabs1( a( imax, k ) )
569  ELSE
570  colmax = zero
571  END IF
572 *
573  IF( ( max( absakk, colmax ).EQ.zero ) ) THEN
574 *
575 * Column K is zero or underflow: set INFO and continue
576 *
577  IF( info.EQ.0 )
578  $ info = k
579  kp = k
580  ELSE
581 *
582 * Test for interchange
583 *
584 * Equivalent to testing for (used to handle NaN and Inf)
585 * ABSAKK.GE.ALPHA*COLMAX
586 *
587  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
588 *
589 * no interchange, use 1-by-1 pivot block
590 *
591  kp = k
592  ELSE
593 *
594  done = .false.
595 *
596 * Loop until pivot found
597 *
598  42 CONTINUE
599 *
600 * Begin pivot search loop body
601 *
602 * JMAX is the column-index of the largest off-diagonal
603 * element in row IMAX, and ROWMAX is its absolute value.
604 * Determine both ROWMAX and JMAX.
605 *
606  IF( imax.NE.k ) THEN
607  jmax = k - 1 + izamax( imax-k, a( imax, k ), lda )
608  rowmax = cabs1( a( imax, jmax ) )
609  ELSE
610  rowmax = zero
611  END IF
612 *
613  IF( imax.LT.n ) THEN
614  itemp = imax + izamax( n-imax, a( imax+1, imax ),
615  $ 1 )
616  dtemp = cabs1( a( itemp, imax ) )
617  IF( dtemp.GT.rowmax ) THEN
618  rowmax = dtemp
619  jmax = itemp
620  END IF
621  END IF
622 *
623 * Equivalent to testing for (used to handle NaN and Inf)
624 * CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
625 *
626  IF( .NOT.( cabs1(a( imax, imax )).LT.alpha*rowmax ) )
627  $ THEN
628 *
629 * interchange rows and columns K and IMAX,
630 * use 1-by-1 pivot block
631 *
632  kp = imax
633  done = .true.
634 *
635 * Equivalent to testing for ROWMAX .EQ. COLMAX,
636 * used to handle NaN and Inf
637 *
638  ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
639 *
640 * interchange rows and columns K+1 and IMAX,
641 * use 2-by-2 pivot block
642 *
643  kp = imax
644  kstep = 2
645  done = .true.
646  ELSE
647 *
648 * Pivot NOT found, set variables and repeat
649 *
650  p = imax
651  colmax = rowmax
652  imax = jmax
653  END IF
654 *
655 * End pivot search loop body
656 *
657  IF( .NOT. done ) GOTO 42
658 *
659  END IF
660 *
661 * Swap TWO rows and TWO columns
662 *
663 * First swap
664 *
665  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
666 *
667 * Interchange rows and column K and P in the trailing
668 * submatrix A(k:n,k:n) if we have a 2-by-2 pivot
669 *
670  IF( p.LT.n )
671  $ CALL zswap( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
672  IF( p.GT.(k+1) )
673  $ CALL zswap( p-k-1, a( k+1, k ), 1, a( p, k+1 ), lda )
674  t = a( k, k )
675  a( k, k ) = a( p, p )
676  a( p, p ) = t
677  END IF
678 *
679 * Second swap
680 *
681  kk = k + kstep - 1
682  IF( kp.NE.kk ) THEN
683 *
684 * Interchange rows and columns KK and KP in the trailing
685 * submatrix A(k:n,k:n)
686 *
687  IF( kp.LT.n )
688  $ CALL zswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
689  IF( ( kk.LT.n ) .AND. ( kp.GT.(kk+1) ) )
690  $ CALL zswap( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
691  $ lda )
692  t = a( kk, kk )
693  a( kk, kk ) = a( kp, kp )
694  a( kp, kp ) = t
695  IF( kstep.EQ.2 ) THEN
696  t = a( k+1, k )
697  a( k+1, k ) = a( kp, k )
698  a( kp, k ) = t
699  END IF
700  END IF
701 *
702 * Update the trailing submatrix
703 *
704  IF( kstep.EQ.1 ) THEN
705 *
706 * 1-by-1 pivot block D(k): column k now holds
707 *
708 * W(k) = L(k)*D(k)
709 *
710 * where L(k) is the k-th column of L
711 *
712  IF( k.LT.n ) THEN
713 *
714 * Perform a rank-1 update of A(k+1:n,k+1:n) and
715 * store L(k) in column k
716 *
717  IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
718 *
719 * Perform a rank-1 update of A(k+1:n,k+1:n) as
720 * A := A - L(k)*D(k)*L(k)**T
721 * = A - W(k)*(1/D(k))*W(k)**T
722 *
723  d11 = cone / a( k, k )
724  CALL zsyr( uplo, n-k, -d11, a( k+1, k ), 1,
725  $ a( k+1, k+1 ), lda )
726 *
727 * Store L(k) in column k
728 *
729  CALL zscal( n-k, d11, a( k+1, k ), 1 )
730  ELSE
731 *
732 * Store L(k) in column k
733 *
734  d11 = a( k, k )
735  DO 46 ii = k + 1, n
736  a( ii, k ) = a( ii, k ) / d11
737  46 CONTINUE
738 *
739 * Perform a rank-1 update of A(k+1:n,k+1:n) as
740 * A := A - L(k)*D(k)*L(k)**T
741 * = A - W(k)*(1/D(k))*W(k)**T
742 * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
743 *
744  CALL zsyr( uplo, n-k, -d11, a( k+1, k ), 1,
745  $ a( k+1, k+1 ), lda )
746  END IF
747  END IF
748 *
749  ELSE
750 *
751 * 2-by-2 pivot block D(k): columns k and k+1 now hold
752 *
753 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
754 *
755 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
756 * of L
757 *
758 *
759 * Perform a rank-2 update of A(k+2:n,k+2:n) as
760 *
761 * A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
762 * = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
763 *
764 * and store L(k) and L(k+1) in columns k and k+1
765 *
766  IF( k.LT.n-1 ) THEN
767 *
768  d21 = a( k+1, k )
769  d11 = a( k+1, k+1 ) / d21
770  d22 = a( k, k ) / d21
771  t = cone / ( d11*d22-cone )
772 *
773  DO 60 j = k + 2, n
774 *
775 * Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
776 *
777  wk = t*( d11*a( j, k )-a( j, k+1 ) )
778  wkp1 = t*( d22*a( j, k+1 )-a( j, k ) )
779 *
780 * Perform a rank-2 update of A(k+2:n,k+2:n)
781 *
782  DO 50 i = j, n
783  a( i, j ) = a( i, j ) - ( a( i, k ) / d21 )*wk -
784  $ ( a( i, k+1 ) / d21 )*wkp1
785  50 CONTINUE
786 *
787 * Store L(k) and L(k+1) in cols k and k+1 for row J
788 *
789  a( j, k ) = wk / d21
790  a( j, k+1 ) = wkp1 / d21
791 *
792  60 CONTINUE
793 *
794  END IF
795 *
796  END IF
797  END IF
798 *
799 * Store details of the interchanges in IPIV
800 *
801  IF( kstep.EQ.1 ) THEN
802  ipiv( k ) = kp
803  ELSE
804  ipiv( k ) = -p
805  ipiv( k+1 ) = -kp
806  END IF
807 *
808 * Increase K and return to the start of the main loop
809 *
810  k = k + kstep
811  GO TO 40
812 *
813  END IF
814 *
815  70 CONTINUE
816 *
817  RETURN
818 *
819 * End of ZSYTF2_ROOK
820 *
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zsyr(UPLO, N, ALPHA, X, INCX, A, LDA)
ZSYR performs the symmetric rank-1 update of a complex symmetric matrix.
Definition: zsyr.f:137
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:54
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zsytrf ( character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

ZSYTRF

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

Purpose:
 ZSYTRF computes the factorization of a complex symmetric matrix A
 using the Bunch-Kaufman diagonal pivoting method.  The form of the
 factorization is

    A = U*D*U**T  or  A = L*D*L**T

 where U (or L) is a product of permutation and unit upper (lower)
 triangular matrices, and D is symmetric and block diagonal with
 with 1-by-1 and 2-by-2 diagonal blocks.

 This is the blocked version of the algorithm, calling Level 3 BLAS.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[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 symmetric matrix A.  If UPLO = 'U', the leading
          N-by-N upper triangular part of A contains the upper
          triangular part of the matrix A, and the strictly lower
          triangular part of A is not referenced.  If UPLO = 'L', the
          leading N-by-N lower triangular part of A contains the lower
          triangular part of the matrix A, and the strictly upper
          triangular part of A is not referenced.

          On exit, the block diagonal matrix D and the multipliers used
          to obtain the factor U or L (see below for further details).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D.
          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
          interchanged and D(k,k) is a 1-by-1 diagonal block.
          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
[out]WORK
          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The length of WORK.  LWORK >=1.  For best performance
          LWORK >= N*NB, where NB is the block size returned by ILAENV.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[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) is exactly zero.  The factorization
                has been completed, but the block diagonal matrix D is
                exactly singular, and division by zero will occur if it
                is used to solve a system of equations.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Further Details:
  If UPLO = 'U', then A = U*D*U**T, where
     U = P(n)*U(n)* ... *P(k)U(k)* ...,
  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
  that if the diagonal block D(k) is of order s (s = 1 or 2), then

             (   I    v    0   )   k-s
     U(k) =  (   0    I    0   )   s
             (   0    0    I   )   n-k
                k-s   s   n-k

  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
  and A(k,k), and v overwrites A(1:k-2,k-1:k).

  If UPLO = 'L', then A = L*D*L**T, where
     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
  that if the diagonal block D(k) is of order s (s = 1 or 2), then

             (   I    0     0   )  k-1
     L(k) =  (   0    I     0   )  s
             (   0    v     I   )  n-k-s+1
                k-1   s  n-k-s+1

  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).

Definition at line 184 of file zsytrf.f.

184 *
185 * -- LAPACK computational routine (version 3.4.0) --
186 * -- LAPACK is a software package provided by Univ. of Tennessee, --
187 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188 * November 2011
189 *
190 * .. Scalar Arguments ..
191  CHARACTER uplo
192  INTEGER info, lda, lwork, n
193 * ..
194 * .. Array Arguments ..
195  INTEGER ipiv( * )
196  COMPLEX*16 a( lda, * ), work( * )
197 * ..
198 *
199 * =====================================================================
200 *
201 * .. Local Scalars ..
202  LOGICAL lquery, upper
203  INTEGER iinfo, iws, j, k, kb, ldwork, lwkopt, nb, nbmin
204 * ..
205 * .. External Functions ..
206  LOGICAL lsame
207  INTEGER ilaenv
208  EXTERNAL lsame, ilaenv
209 * ..
210 * .. External Subroutines ..
211  EXTERNAL xerbla, zlasyf, zsytf2
212 * ..
213 * .. Intrinsic Functions ..
214  INTRINSIC max
215 * ..
216 * .. Executable Statements ..
217 *
218 * Test the input parameters.
219 *
220  info = 0
221  upper = lsame( uplo, 'U' )
222  lquery = ( lwork.EQ.-1 )
223  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
224  info = -1
225  ELSE IF( n.LT.0 ) THEN
226  info = -2
227  ELSE IF( lda.LT.max( 1, n ) ) THEN
228  info = -4
229  ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
230  info = -7
231  END IF
232 *
233  IF( info.EQ.0 ) THEN
234 *
235 * Determine the block size
236 *
237  nb = ilaenv( 1, 'ZSYTRF', uplo, n, -1, -1, -1 )
238  lwkopt = n*nb
239  work( 1 ) = lwkopt
240  END IF
241 *
242  IF( info.NE.0 ) THEN
243  CALL xerbla( 'ZSYTRF', -info )
244  RETURN
245  ELSE IF( lquery ) THEN
246  RETURN
247  END IF
248 *
249  nbmin = 2
250  ldwork = n
251  IF( nb.GT.1 .AND. nb.LT.n ) THEN
252  iws = ldwork*nb
253  IF( lwork.LT.iws ) THEN
254  nb = max( lwork / ldwork, 1 )
255  nbmin = max( 2, ilaenv( 2, 'ZSYTRF', uplo, n, -1, -1, -1 ) )
256  END IF
257  ELSE
258  iws = 1
259  END IF
260  IF( nb.LT.nbmin )
261  $ nb = n
262 *
263  IF( upper ) THEN
264 *
265 * Factorize A as U*D*U**T using the upper triangle of A
266 *
267 * K is the main loop index, decreasing from N to 1 in steps of
268 * KB, where KB is the number of columns factorized by ZLASYF;
269 * KB is either NB or NB-1, or K for the last block
270 *
271  k = n
272  10 CONTINUE
273 *
274 * If K < 1, exit from loop
275 *
276  IF( k.LT.1 )
277  $ GO TO 40
278 *
279  IF( k.GT.nb ) THEN
280 *
281 * Factorize columns k-kb+1:k of A and use blocked code to
282 * update columns 1:k-kb
283 *
284  CALL zlasyf( uplo, k, nb, kb, a, lda, ipiv, work, n, iinfo )
285  ELSE
286 *
287 * Use unblocked code to factorize columns 1:k of A
288 *
289  CALL zsytf2( uplo, k, a, lda, ipiv, iinfo )
290  kb = k
291  END IF
292 *
293 * Set INFO on the first occurrence of a zero pivot
294 *
295  IF( info.EQ.0 .AND. iinfo.GT.0 )
296  $ info = iinfo
297 *
298 * Decrease K and return to the start of the main loop
299 *
300  k = k - kb
301  GO TO 10
302 *
303  ELSE
304 *
305 * Factorize A as L*D*L**T using the lower triangle of A
306 *
307 * K is the main loop index, increasing from 1 to N in steps of
308 * KB, where KB is the number of columns factorized by ZLASYF;
309 * KB is either NB or NB-1, or N-K+1 for the last block
310 *
311  k = 1
312  20 CONTINUE
313 *
314 * If K > N, exit from loop
315 *
316  IF( k.GT.n )
317  $ GO TO 40
318 *
319  IF( k.LE.n-nb ) THEN
320 *
321 * Factorize columns k:k+kb-1 of A and use blocked code to
322 * update columns k+kb:n
323 *
324  CALL zlasyf( uplo, n-k+1, nb, kb, a( k, k ), lda, ipiv( k ),
325  $ work, n, iinfo )
326  ELSE
327 *
328 * Use unblocked code to factorize columns k:n of A
329 *
330  CALL zsytf2( uplo, n-k+1, a( k, k ), lda, ipiv( k ), iinfo )
331  kb = n - k + 1
332  END IF
333 *
334 * Set INFO on the first occurrence of a zero pivot
335 *
336  IF( info.EQ.0 .AND. iinfo.GT.0 )
337  $ info = iinfo + k - 1
338 *
339 * Adjust IPIV
340 *
341  DO 30 j = k, k + kb - 1
342  IF( ipiv( j ).GT.0 ) THEN
343  ipiv( j ) = ipiv( j ) + k - 1
344  ELSE
345  ipiv( j ) = ipiv( j ) - k + 1
346  END IF
347  30 CONTINUE
348 *
349 * Increase K and return to the start of the main loop
350 *
351  k = k + kb
352  GO TO 20
353 *
354  END IF
355 *
356  40 CONTINUE
357  work( 1 ) = lwkopt
358  RETURN
359 *
360 * End of ZSYTRF
361 *
subroutine zlasyf(UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
ZLASYF computes a partial factorization of a complex symmetric matrix using the Bunch-Kaufman diagona...
Definition: zlasyf.f:179
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zsytf2(UPLO, N, A, LDA, IPIV, INFO)
ZSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting ...
Definition: zsytf2.f:193
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zsytrf_rook ( character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

ZSYTRF_ROOK

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

Purpose:
 ZSYTRF_ROOK computes the factorization of a complex symmetric matrix A
 using the bounded Bunch-Kaufman ("rook") diagonal pivoting method.
 The form of the factorization is

    A = U*D*U**T  or  A = L*D*L**T

 where U (or L) is a product of permutation and unit upper (lower)
 triangular matrices, and D is symmetric and block diagonal with
 1-by-1 and 2-by-2 diagonal blocks.

 This is the blocked version of the algorithm, calling Level 3 BLAS.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[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 symmetric matrix A.  If UPLO = 'U', the leading
          N-by-N upper triangular part of A contains the upper
          triangular part of the matrix A, and the strictly lower
          triangular part of A is not referenced.  If UPLO = 'L', the
          leading N-by-N lower triangular part of A contains the lower
          triangular part of the matrix A, and the strictly upper
          triangular part of A is not referenced.

          On exit, the block diagonal matrix D and the multipliers used
          to obtain the factor U or L (see below for further details).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D.

          If UPLO = 'U':
               If IPIV(k) > 0, then rows and columns k and IPIV(k)
               were interchanged and D(k,k) is a 1-by-1 diagonal block.

               If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
               columns k and -IPIV(k) were interchanged and rows and
               columns k-1 and -IPIV(k-1) were inerchaged,
               D(k-1:k,k-1:k) is a 2-by-2 diagonal block.

          If UPLO = 'L':
               If IPIV(k) > 0, then rows and columns k and IPIV(k)
               were interchanged and D(k,k) is a 1-by-1 diagonal block.

               If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
               columns k and -IPIV(k) were interchanged and rows and
               columns k+1 and -IPIV(k+1) were inerchaged,
               D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
[out]WORK
          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)).
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The length of WORK.  LWORK >=1.  For best performance
          LWORK >= N*NB, where NB is the block size returned by ILAENV.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[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) is exactly zero.  The factorization
                has been completed, but the block diagonal matrix D is
                exactly singular, and division by zero will occur if it
                is used to solve a system of equations.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015
Further Details:
  If UPLO = 'U', then A = U*D*U**T, where
     U = P(n)*U(n)* ... *P(k)U(k)* ...,
  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
  that if the diagonal block D(k) is of order s (s = 1 or 2), then

             (   I    v    0   )   k-s
     U(k) =  (   0    I    0   )   s
             (   0    0    I   )   n-k
                k-s   s   n-k

  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
  and A(k,k), and v overwrites A(1:k-2,k-1:k).

  If UPLO = 'L', then A = L*D*L**T, where
     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
  that if the diagonal block D(k) is of order s (s = 1 or 2), then

             (   I    0     0   )  k-1
     L(k) =  (   0    I     0   )  s
             (   0    v     I   )  n-k-s+1
                k-1   s  n-k-s+1

  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
Contributors:
   November 2015, 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 210 of file zsytrf_rook.f.

210 *
211 * -- LAPACK computational routine (version 3.6.0) --
212 * -- LAPACK is a software package provided by Univ. of Tennessee, --
213 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
214 * November 2015
215 *
216 * .. Scalar Arguments ..
217  CHARACTER uplo
218  INTEGER info, lda, lwork, n
219 * ..
220 * .. Array Arguments ..
221  INTEGER ipiv( * )
222  COMPLEX*16 a( lda, * ), work( * )
223 * ..
224 *
225 * =====================================================================
226 *
227 * .. Local Scalars ..
228  LOGICAL lquery, upper
229  INTEGER iinfo, iws, j, k, kb, ldwork, lwkopt, nb, nbmin
230 * ..
231 * .. External Functions ..
232  LOGICAL lsame
233  INTEGER ilaenv
234  EXTERNAL lsame, ilaenv
235 * ..
236 * .. External Subroutines ..
237  EXTERNAL zlasyf_rook, zsytf2_rook, xerbla
238 * ..
239 * .. Intrinsic Functions ..
240  INTRINSIC max
241 * ..
242 * .. Executable Statements ..
243 *
244 * Test the input parameters.
245 *
246  info = 0
247  upper = lsame( uplo, 'U' )
248  lquery = ( lwork.EQ.-1 )
249  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
250  info = -1
251  ELSE IF( n.LT.0 ) THEN
252  info = -2
253  ELSE IF( lda.LT.max( 1, n ) ) THEN
254  info = -4
255  ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
256  info = -7
257  END IF
258 *
259  IF( info.EQ.0 ) THEN
260 *
261 * Determine the block size
262 *
263  nb = ilaenv( 1, 'ZSYTRF_ROOK', uplo, n, -1, -1, -1 )
264  lwkopt = n*nb
265  work( 1 ) = lwkopt
266  END IF
267 *
268  IF( info.NE.0 ) THEN
269  CALL xerbla( 'ZSYTRF_ROOK', -info )
270  RETURN
271  ELSE IF( lquery ) THEN
272  RETURN
273  END IF
274 *
275  nbmin = 2
276  ldwork = n
277  IF( nb.GT.1 .AND. nb.LT.n ) THEN
278  iws = ldwork*nb
279  IF( lwork.LT.iws ) THEN
280  nb = max( lwork / ldwork, 1 )
281  nbmin = max( 2, ilaenv( 2, 'ZSYTRF_ROOK',
282  $ uplo, n, -1, -1, -1 ) )
283  END IF
284  ELSE
285  iws = 1
286  END IF
287  IF( nb.LT.nbmin )
288  $ nb = n
289 *
290  IF( upper ) THEN
291 *
292 * Factorize A as U*D*U**T using the upper triangle of A
293 *
294 * K is the main loop index, decreasing from N to 1 in steps of
295 * KB, where KB is the number of columns factorized by ZLASYF_ROOK;
296 * KB is either NB or NB-1, or K for the last block
297 *
298  k = n
299  10 CONTINUE
300 *
301 * If K < 1, exit from loop
302 *
303  IF( k.LT.1 )
304  $ GO TO 40
305 *
306  IF( k.GT.nb ) THEN
307 *
308 * Factorize columns k-kb+1:k of A and use blocked code to
309 * update columns 1:k-kb
310 *
311  CALL zlasyf_rook( uplo, k, nb, kb, a, lda,
312  $ ipiv, work, ldwork, iinfo )
313  ELSE
314 *
315 * Use unblocked code to factorize columns 1:k of A
316 *
317  CALL zsytf2_rook( uplo, k, a, lda, ipiv, iinfo )
318  kb = k
319  END IF
320 *
321 * Set INFO on the first occurrence of a zero pivot
322 *
323  IF( info.EQ.0 .AND. iinfo.GT.0 )
324  $ info = iinfo
325 *
326 * No need to adjust IPIV
327 *
328 * Decrease K and return to the start of the main loop
329 *
330  k = k - kb
331  GO TO 10
332 *
333  ELSE
334 *
335 * Factorize A as L*D*L**T using the lower triangle of A
336 *
337 * K is the main loop index, increasing from 1 to N in steps of
338 * KB, where KB is the number of columns factorized by ZLASYF_ROOK;
339 * KB is either NB or NB-1, or N-K+1 for the last block
340 *
341  k = 1
342  20 CONTINUE
343 *
344 * If K > N, exit from loop
345 *
346  IF( k.GT.n )
347  $ GO TO 40
348 *
349  IF( k.LE.n-nb ) THEN
350 *
351 * Factorize columns k:k+kb-1 of A and use blocked code to
352 * update columns k+kb:n
353 *
354  CALL zlasyf_rook( uplo, n-k+1, nb, kb, a( k, k ), lda,
355  $ ipiv( k ), work, ldwork, iinfo )
356  ELSE
357 *
358 * Use unblocked code to factorize columns k:n of A
359 *
360  CALL zsytf2_rook( uplo, n-k+1, a( k, k ), lda, ipiv( k ),
361  $ iinfo )
362  kb = n - k + 1
363  END IF
364 *
365 * Set INFO on the first occurrence of a zero pivot
366 *
367  IF( info.EQ.0 .AND. iinfo.GT.0 )
368  $ info = iinfo + k - 1
369 *
370 * Adjust IPIV
371 *
372  DO 30 j = k, k + kb - 1
373  IF( ipiv( j ).GT.0 ) THEN
374  ipiv( j ) = ipiv( j ) + k - 1
375  ELSE
376  ipiv( j ) = ipiv( j ) - k + 1
377  END IF
378  30 CONTINUE
379 *
380 * Increase K and return to the start of the main loop
381 *
382  k = k + kb
383  GO TO 20
384 *
385  END IF
386 *
387  40 CONTINUE
388  work( 1 ) = lwkopt
389  RETURN
390 *
391 * End of ZSYTRF_ROOK
392 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine zlasyf_rook(UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
ZLASYF_ROOK computes a partial factorization of a complex symmetric matrix using the bounded Bunch-Ka...
Definition: zlasyf_rook.f:186
subroutine zsytf2_rook(UPLO, N, A, LDA, IPIV, INFO)
ZSYTF2_ROOK computes the factorization of a complex symmetric indefinite matrix using the bounded Bun...
Definition: zsytf2_rook.f:196

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zsytri ( character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZSYTRI

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

Purpose:
 ZSYTRI computes the inverse of a complex symmetric indefinite matrix
 A using the factorization A = U*D*U**T or A = L*D*L**T computed by
 ZSYTRF.
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.

          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.
[out]WORK
          WORK is COMPLEX*16 array, dimension (2*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.
Date
November 2011

Definition at line 116 of file zsytri.f.

116 *
117 * -- LAPACK computational routine (version 3.4.0) --
118 * -- LAPACK is a software package provided by Univ. of Tennessee, --
119 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
120 * November 2011
121 *
122 * .. Scalar Arguments ..
123  CHARACTER uplo
124  INTEGER info, lda, n
125 * ..
126 * .. Array Arguments ..
127  INTEGER ipiv( * )
128  COMPLEX*16 a( lda, * ), work( * )
129 * ..
130 *
131 * =====================================================================
132 *
133 * .. Parameters ..
134  COMPLEX*16 one, zero
135  parameter( one = ( 1.0d+0, 0.0d+0 ),
136  $ zero = ( 0.0d+0, 0.0d+0 ) )
137 * ..
138 * .. Local Scalars ..
139  LOGICAL upper
140  INTEGER k, kp, kstep
141  COMPLEX*16 ak, akkp1, akp1, d, t, temp
142 * ..
143 * .. External Functions ..
144  LOGICAL lsame
145  COMPLEX*16 zdotu
146  EXTERNAL lsame, zdotu
147 * ..
148 * .. External Subroutines ..
149  EXTERNAL xerbla, zcopy, zswap, zsymv
150 * ..
151 * .. Intrinsic Functions ..
152  INTRINSIC abs, max
153 * ..
154 * .. Executable Statements ..
155 *
156 * Test the input parameters.
157 *
158  info = 0
159  upper = lsame( uplo, 'U' )
160  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
161  info = -1
162  ELSE IF( n.LT.0 ) THEN
163  info = -2
164  ELSE IF( lda.LT.max( 1, n ) ) THEN
165  info = -4
166  END IF
167  IF( info.NE.0 ) THEN
168  CALL xerbla( 'ZSYTRI', -info )
169  RETURN
170  END IF
171 *
172 * Quick return if possible
173 *
174  IF( n.EQ.0 )
175  $ RETURN
176 *
177 * Check that the diagonal matrix D is nonsingular.
178 *
179  IF( upper ) THEN
180 *
181 * Upper triangular storage: examine D from bottom to top
182 *
183  DO 10 info = n, 1, -1
184  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
185  $ RETURN
186  10 CONTINUE
187  ELSE
188 *
189 * Lower triangular storage: examine D from top to bottom.
190 *
191  DO 20 info = 1, n
192  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
193  $ RETURN
194  20 CONTINUE
195  END IF
196  info = 0
197 *
198  IF( upper ) THEN
199 *
200 * Compute inv(A) from the factorization A = U*D*U**T.
201 *
202 * K is the main loop index, increasing from 1 to N in steps of
203 * 1 or 2, depending on the size of the diagonal blocks.
204 *
205  k = 1
206  30 CONTINUE
207 *
208 * If K > N, exit from loop.
209 *
210  IF( k.GT.n )
211  $ GO TO 40
212 *
213  IF( ipiv( k ).GT.0 ) THEN
214 *
215 * 1 x 1 diagonal block
216 *
217 * Invert the diagonal block.
218 *
219  a( k, k ) = one / a( k, k )
220 *
221 * Compute column K of the inverse.
222 *
223  IF( k.GT.1 ) THEN
224  CALL zcopy( k-1, a( 1, k ), 1, work, 1 )
225  CALL zsymv( uplo, k-1, -one, a, lda, work, 1, zero,
226  $ a( 1, k ), 1 )
227  a( k, k ) = a( k, k ) - zdotu( k-1, work, 1, a( 1, k ),
228  $ 1 )
229  END IF
230  kstep = 1
231  ELSE
232 *
233 * 2 x 2 diagonal block
234 *
235 * Invert the diagonal block.
236 *
237  t = a( k, k+1 )
238  ak = a( k, k ) / t
239  akp1 = a( k+1, k+1 ) / t
240  akkp1 = a( k, k+1 ) / t
241  d = t*( ak*akp1-one )
242  a( k, k ) = akp1 / d
243  a( k+1, k+1 ) = ak / d
244  a( k, k+1 ) = -akkp1 / d
245 *
246 * Compute columns K and K+1 of the inverse.
247 *
248  IF( k.GT.1 ) THEN
249  CALL zcopy( k-1, a( 1, k ), 1, work, 1 )
250  CALL zsymv( uplo, k-1, -one, a, lda, work, 1, zero,
251  $ a( 1, k ), 1 )
252  a( k, k ) = a( k, k ) - zdotu( k-1, work, 1, a( 1, k ),
253  $ 1 )
254  a( k, k+1 ) = a( k, k+1 ) -
255  $ zdotu( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
256  CALL zcopy( k-1, a( 1, k+1 ), 1, work, 1 )
257  CALL zsymv( uplo, k-1, -one, a, lda, work, 1, zero,
258  $ a( 1, k+1 ), 1 )
259  a( k+1, k+1 ) = a( k+1, k+1 ) -
260  $ zdotu( k-1, work, 1, a( 1, k+1 ), 1 )
261  END IF
262  kstep = 2
263  END IF
264 *
265  kp = abs( ipiv( k ) )
266  IF( kp.NE.k ) THEN
267 *
268 * Interchange rows and columns K and KP in the leading
269 * submatrix A(1:k+1,1:k+1)
270 *
271  CALL zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
272  CALL zswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
273  temp = a( k, k )
274  a( k, k ) = a( kp, kp )
275  a( kp, kp ) = temp
276  IF( kstep.EQ.2 ) THEN
277  temp = a( k, k+1 )
278  a( k, k+1 ) = a( kp, k+1 )
279  a( kp, k+1 ) = temp
280  END IF
281  END IF
282 *
283  k = k + kstep
284  GO TO 30
285  40 CONTINUE
286 *
287  ELSE
288 *
289 * Compute inv(A) from the factorization A = L*D*L**T.
290 *
291 * K is the main loop index, increasing from 1 to N in steps of
292 * 1 or 2, depending on the size of the diagonal blocks.
293 *
294  k = n
295  50 CONTINUE
296 *
297 * If K < 1, exit from loop.
298 *
299  IF( k.LT.1 )
300  $ GO TO 60
301 *
302  IF( ipiv( k ).GT.0 ) THEN
303 *
304 * 1 x 1 diagonal block
305 *
306 * Invert the diagonal block.
307 *
308  a( k, k ) = one / a( k, k )
309 *
310 * Compute column K of the inverse.
311 *
312  IF( k.LT.n ) THEN
313  CALL zcopy( n-k, a( k+1, k ), 1, work, 1 )
314  CALL zsymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
315  $ zero, a( k+1, k ), 1 )
316  a( k, k ) = a( k, k ) - zdotu( n-k, work, 1, a( k+1, k ),
317  $ 1 )
318  END IF
319  kstep = 1
320  ELSE
321 *
322 * 2 x 2 diagonal block
323 *
324 * Invert the diagonal block.
325 *
326  t = a( k, k-1 )
327  ak = a( k-1, k-1 ) / t
328  akp1 = a( k, k ) / t
329  akkp1 = a( k, k-1 ) / t
330  d = t*( ak*akp1-one )
331  a( k-1, k-1 ) = akp1 / d
332  a( k, k ) = ak / d
333  a( k, k-1 ) = -akkp1 / d
334 *
335 * Compute columns K-1 and K of the inverse.
336 *
337  IF( k.LT.n ) THEN
338  CALL zcopy( n-k, a( k+1, k ), 1, work, 1 )
339  CALL zsymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
340  $ zero, a( k+1, k ), 1 )
341  a( k, k ) = a( k, k ) - zdotu( n-k, work, 1, a( k+1, k ),
342  $ 1 )
343  a( k, k-1 ) = a( k, k-1 ) -
344  $ zdotu( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
345  $ 1 )
346  CALL zcopy( n-k, a( k+1, k-1 ), 1, work, 1 )
347  CALL zsymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
348  $ zero, a( k+1, k-1 ), 1 )
349  a( k-1, k-1 ) = a( k-1, k-1 ) -
350  $ zdotu( n-k, work, 1, a( k+1, k-1 ), 1 )
351  END IF
352  kstep = 2
353  END IF
354 *
355  kp = abs( ipiv( k ) )
356  IF( kp.NE.k ) THEN
357 *
358 * Interchange rows and columns K and KP