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

Functions

subroutine dla_syamv (UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
 DLA_SYAMV computes a matrix-vector product using a symmetric indefinite matrix to calculate error bounds. More...
 
double precision function dla_syrcond (UPLO, N, A, LDA, AF, LDAF, IPIV, CMODE, C, INFO, WORK, IWORK)
 DLA_SYRCOND estimates the Skeel condition number for a symmetric indefinite matrix. More...
 
subroutine dla_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)
 DLA_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 dla_syrpvgrw (UPLO, N, INFO, A, LDA, AF, LDAF, IPIV, WORK)
 DLA_SYRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a symmetric indefinite matrix. More...
 
subroutine dlasyf (UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
 DLASYF computes a partial factorization of a real symmetric matrix using the Bunch-Kaufman diagonal pivoting method. More...
 
subroutine dlasyf_rook (UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
 DLASYF_ROOK *> DLASYF_ROOK computes a partial factorization of a real symmetric matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method. More...
 
subroutine dsycon (UPLO, N, A, LDA, IPIV, ANORM, RCOND, WORK, IWORK, INFO)
 DSYCON More...
 
subroutine dsycon_rook (UPLO, N, A, LDA, IPIV, ANORM, RCOND, WORK, IWORK, INFO)
 DSYCON_ROOK More...
 
subroutine dsyconv (UPLO, WAY, N, A, LDA, IPIV, E, INFO)
 DSYCONV More...
 
subroutine dsyequb (UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO)
 DSYEQUB More...
 
subroutine dsygs2 (ITYPE, UPLO, N, A, LDA, B, LDB, INFO)
 DSYGS2 reduces a symmetric definite generalized eigenproblem to standard form, using the factorization results obtained from spotrf (unblocked algorithm). More...
 
subroutine dsygst (ITYPE, UPLO, N, A, LDA, B, LDB, INFO)
 DSYGST More...
 
subroutine dsyrfs (UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO)
 DSYRFS More...
 
subroutine dsyrfsx (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, IWORK, INFO)
 DSYRFSX More...
 
subroutine dsytd2 (UPLO, N, A, LDA, D, E, TAU, INFO)
 DSYTD2 reduces a symmetric matrix to real symmetric tridiagonal form by an orthogonal similarity transformation (unblocked algorithm). More...
 
subroutine dsytf2 (UPLO, N, A, LDA, IPIV, INFO)
 DSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting method (unblocked algorithm). More...
 
subroutine dsytf2_rook (UPLO, N, A, LDA, IPIV, INFO)
 DSYTF2_ROOK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method (unblocked algorithm). More...
 
subroutine dsytrd (UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
 DSYTRD More...
 
subroutine dsytrf (UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO)
 DSYTRF More...
 
subroutine dsytrf_rook (UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO)
 DSYTRF_ROOK More...
 
subroutine dsytri (UPLO, N, A, LDA, IPIV, WORK, INFO)
 DSYTRI More...
 
subroutine dsytri2 (UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO)
 DSYTRI2 More...
 
subroutine dsytri2x (UPLO, N, A, LDA, IPIV, WORK, NB, INFO)
 DSYTRI2X More...
 
subroutine dsytri_rook (UPLO, N, A, LDA, IPIV, WORK, INFO)
 DSYTRI_ROOK More...
 
subroutine dsytrs (UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
 DSYTRS More...
 
subroutine dsytrs2 (UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, INFO)
 DSYTRS2 More...
 
subroutine dsytrs_rook (UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
 DSYTRS_ROOK More...
 
subroutine dtgsyl (TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)
 DTGSYL More...
 
subroutine dtrsyl (TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, LDC, SCALE, INFO)
 DTRSYL More...
 

Detailed Description

This is the group of double computational functions for SY matrices

Function Documentation

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

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

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

Purpose:
 DLA_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 DOUBLE PRECISION array of 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 DOUBLE PRECISION array, dimension
           ( 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 179 of file dla_syamv.f.

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

DLA_SYRCOND estimates the Skeel condition number for a symmetric indefinite matrix.

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

Purpose:
    DLA_SYRCOND estimates the Skeel condition number of  op(A) * op2(C)
    where op2 is determined by CMODE as follows
    CMODE =  1    op2(C) = C
    CMODE =  0    op2(C) = I
    CMODE = -1    op2(C) = inv(C)
    The Skeel condition number cond(A) = norminf( |inv(A)||A| )
    is computed by computing scaling factors R such that
    diag(R)*A*op2(C) is row equilibrated and computing the standard
    infinity-norm condition number.
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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDAF,N)
     The block diagonal matrix D and the multipliers used to
     obtain the factor U or L 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]CMODE
          CMODE is INTEGER
     Determines op2(C) in the formula op(A) * op2(C) as follows:
     CMODE =  1    op2(C) = C
     CMODE =  0    op2(C) = I
     CMODE = -1    op2(C) = inv(C)
[in]C
          C is DOUBLE PRECISION array, dimension (N)
     The vector C in the formula op(A) * op2(C).
[out]INFO
          INFO is INTEGER
       = 0:  Successful exit.
     i > 0:  The ith argument is invalid.
[in]WORK
          WORK is DOUBLE PRECISION array, dimension (3*N).
     Workspace.
[in]IWORK
          IWORK is INTEGER array, dimension (N).
     Workspace.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 150 of file dla_syrcond.f.

150 *
151 * -- LAPACK computational routine (version 3.4.2) --
152 * -- LAPACK is a software package provided by Univ. of Tennessee, --
153 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
154 * September 2012
155 *
156 * .. Scalar Arguments ..
157  CHARACTER uplo
158  INTEGER n, lda, ldaf, info, cmode
159 * ..
160 * .. Array Arguments
161  INTEGER iwork( * ), ipiv( * )
162  DOUBLE PRECISION a( lda, * ), af( ldaf, * ), work( * ), c( * )
163 * ..
164 *
165 * =====================================================================
166 *
167 * .. Local Scalars ..
168  CHARACTER normin
169  INTEGER kase, i, j
170  DOUBLE PRECISION ainvnm, smlnum, tmp
171  LOGICAL up
172 * ..
173 * .. Local Arrays ..
174  INTEGER isave( 3 )
175 * ..
176 * .. External Functions ..
177  LOGICAL lsame
178  INTEGER idamax
179  DOUBLE PRECISION dlamch
180  EXTERNAL lsame, idamax, dlamch
181 * ..
182 * .. External Subroutines ..
183  EXTERNAL dlacn2, dlatrs, drscl, xerbla, dsytrs
184 * ..
185 * .. Intrinsic Functions ..
186  INTRINSIC abs, max
187 * ..
188 * .. Executable Statements ..
189 *
190  dla_syrcond = 0.0d+0
191 *
192  info = 0
193  IF( n.LT.0 ) THEN
194  info = -2
195  ELSE IF( lda.LT.max( 1, n ) ) THEN
196  info = -4
197  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
198  info = -6
199  END IF
200  IF( info.NE.0 ) THEN
201  CALL xerbla( 'DLA_SYRCOND', -info )
202  RETURN
203  END IF
204  IF( n.EQ.0 ) THEN
205  dla_syrcond = 1.0d+0
206  RETURN
207  END IF
208  up = .false.
209  IF ( lsame( uplo, 'U' ) ) up = .true.
210 *
211 * Compute the equilibration matrix R such that
212 * inv(R)*A*C has unit 1-norm.
213 *
214  IF ( up ) THEN
215  DO i = 1, n
216  tmp = 0.0d+0
217  IF ( cmode .EQ. 1 ) THEN
218  DO j = 1, i
219  tmp = tmp + abs( a( j, i ) * c( j ) )
220  END DO
221  DO j = i+1, n
222  tmp = tmp + abs( a( i, j ) * c( j ) )
223  END DO
224  ELSE IF ( cmode .EQ. 0 ) THEN
225  DO j = 1, i
226  tmp = tmp + abs( a( j, i ) )
227  END DO
228  DO j = i+1, n
229  tmp = tmp + abs( a( i, j ) )
230  END DO
231  ELSE
232  DO j = 1, i
233  tmp = tmp + abs( a( j, i ) / c( j ) )
234  END DO
235  DO j = i+1, n
236  tmp = tmp + abs( a( i, j ) / c( j ) )
237  END DO
238  END IF
239  work( 2*n+i ) = tmp
240  END DO
241  ELSE
242  DO i = 1, n
243  tmp = 0.0d+0
244  IF ( cmode .EQ. 1 ) THEN
245  DO j = 1, i
246  tmp = tmp + abs( a( i, j ) * c( j ) )
247  END DO
248  DO j = i+1, n
249  tmp = tmp + abs( a( j, i ) * c( j ) )
250  END DO
251  ELSE IF ( cmode .EQ. 0 ) THEN
252  DO j = 1, i
253  tmp = tmp + abs( a( i, j ) )
254  END DO
255  DO j = i+1, n
256  tmp = tmp + abs( a( j, i ) )
257  END DO
258  ELSE
259  DO j = 1, i
260  tmp = tmp + abs( a( i, j) / c( j ) )
261  END DO
262  DO j = i+1, n
263  tmp = tmp + abs( a( j, i) / c( j ) )
264  END DO
265  END IF
266  work( 2*n+i ) = tmp
267  END DO
268  ENDIF
269 *
270 * Estimate the norm of inv(op(A)).
271 *
272  smlnum = dlamch( 'Safe minimum' )
273  ainvnm = 0.0d+0
274  normin = 'N'
275 
276  kase = 0
277  10 CONTINUE
278  CALL dlacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
279  IF( kase.NE.0 ) THEN
280  IF( kase.EQ.2 ) THEN
281 *
282 * Multiply by R.
283 *
284  DO i = 1, n
285  work( i ) = work( i ) * work( 2*n+i )
286  END DO
287 
288  IF ( up ) THEN
289  CALL dsytrs( 'U', n, 1, af, ldaf, ipiv, work, n, info )
290  ELSE
291  CALL dsytrs( 'L', n, 1, af, ldaf, ipiv, work, n, info )
292  ENDIF
293 *
294 * Multiply by inv(C).
295 *
296  IF ( cmode .EQ. 1 ) THEN
297  DO i = 1, n
298  work( i ) = work( i ) / c( i )
299  END DO
300  ELSE IF ( cmode .EQ. -1 ) THEN
301  DO i = 1, n
302  work( i ) = work( i ) * c( i )
303  END DO
304  END IF
305  ELSE
306 *
307 * Multiply by inv(C**T).
308 *
309  IF ( cmode .EQ. 1 ) THEN
310  DO i = 1, n
311  work( i ) = work( i ) / c( i )
312  END DO
313  ELSE IF ( cmode .EQ. -1 ) THEN
314  DO i = 1, n
315  work( i ) = work( i ) * c( i )
316  END DO
317  END IF
318 
319  IF ( up ) THEN
320  CALL dsytrs( 'U', n, 1, af, ldaf, ipiv, work, n, info )
321  ELSE
322  CALL dsytrs( 'L', n, 1, af, ldaf, ipiv, work, n, info )
323  ENDIF
324 *
325 * Multiply by R.
326 *
327  DO i = 1, n
328  work( i ) = work( i ) * work( 2*n+i )
329  END DO
330  END IF
331 *
332  GO TO 10
333  END IF
334 *
335 * Compute the estimate of the reciprocal condition number.
336 *
337  IF( ainvnm .NE. 0.0d+0 )
338  $ dla_syrcond = ( 1.0d+0 / ainvnm )
339 *
340  RETURN
341 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
DLATRS solves a triangular system of equations with the scale factor set to prevent overflow...
Definition: dlatrs.f:240
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dsytrs(UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
DSYTRS
Definition: dsytrs.f:122
subroutine dlacn2(N, V, X, ISGN, EST, KASE, ISAVE)
DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: dlacn2.f:138
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:53
subroutine drscl(N, SA, SX, INCX)
DRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: drscl.f:86
double precision function dla_syrcond(UPLO, N, A, LDA, AF, LDAF, IPIV, CMODE, C, INFO, WORK, IWORK)
DLA_SYRCOND estimates the Skeel condition number for a symmetric indefinite matrix.
Definition: dla_syrcond.f:150

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dla_syrfsx_extended ( integer  PREC_TYPE,
character  UPLO,
integer  N,
integer  NRHS,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldaf, * )  AF,
integer  LDAF,
integer, dimension( * )  IPIV,
logical  COLEQU,
double precision, dimension( * )  C,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, 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,
double precision, dimension( * )  RES,
double precision, dimension( * )  AYB,
double precision, dimension( * )  DY,
double precision, dimension( * )  Y_TAIL,
double precision  RCOND,
integer  ITHRESH,
double precision  RTHRESH,
double precision  DZ_UB,
logical  IGNORE_CWISE,
integer  INFO 
)

DLA_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 DLA_SYRFSX_EXTENDED + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DLA_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 DSYRFSX 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDAF,N)
     The block diagonal matrix D and the multipliers used to
     obtain the factor U or L 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]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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension
                    (LDY,NRHS)
     On entry, the solution matrix X, as computed by DSYTRS.
     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 DLA_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 DOUBLE PRECISION array, dimension (N)
     Workspace to hold the intermediate residual.
[in]AYB
          AYB is DOUBLE PRECISION array, dimension (N)
     Workspace. This can be the same workspace passed for Y_TAIL.
[in]DY
          DY is DOUBLE PRECISION array, dimension (N)
     Workspace to hold the intermediate solution.
[in]Y_TAIL
          Y_TAIL is DOUBLE PRECISION 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 DLA_SYRFSX_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 401 of file dla_syrfsx_extended.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

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

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

Purpose:
 DLA_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 DSYTRF, .i.e., the pivot in
     column INFO is exactly 0.
[in]A
          A is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDAF,N)
     The block diagonal matrix D and the multipliers used to
     obtain the factor U or L 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]WORK
          WORK is DOUBLE PRECISION array, dimension (2*N)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 124 of file dla_syrpvgrw.f.

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

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

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

Purpose:
 DLASYF computes a partial factorization of a real 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.

 DLASYF is an auxiliary routine called by DSYTRF. 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 178 of file dlasyf.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dlasyf_rook ( character  UPLO,
integer  N,
integer  NB,
integer  KB,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
double precision, dimension( ldw, * )  W,
integer  LDW,
integer  INFO 
)

DLASYF_ROOK *> DLASYF_ROOK computes a partial factorization of a real symmetric matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method.

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

Purpose:
 DLASYF_ROOK computes a partial factorization of a real 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.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

DSYCON

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

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

 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 DOUBLE PRECISION array, dimension (LDA,N)
          The block diagonal matrix D and the multipliers used to
          obtain the factor U or L as computed by DSYTRF.
[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 DSYTRF.
[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 DOUBLE PRECISION array, dimension (2*N)
[out]IWORK
          IWORK is INTEGER array, dimension (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 132 of file dsycon.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

DSYCON_ROOK

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

Purpose:
 DSYCON_ROOK estimates the reciprocal of the condition number (in the
 1-norm) of a real symmetric matrix A using the factorization
 A = U*D*U**T or A = L*D*L**T computed by DSYTRF_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 DOUBLE PRECISION array, dimension (LDA,N)
          The block diagonal matrix D and the multipliers used to
          obtain the factor U or L as computed by DSYTRF_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 DSYTRF_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 DOUBLE PRECISION array, dimension (2*N)
[out]IWORK
          IWORK is INTEGER array, dimension (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
April 2012
Contributors:

April 2012, 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 146 of file dsycon_rook.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dsyconv ( character  UPLO,
character  WAY,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
double precision, dimension( * )  E,
integer  INFO 
)

DSYCONV

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

Purpose:
 DSYCONV convert A given by TRF into L and D and vice-versa.
 Get Non-diag 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 DOUBLE PRECISION array, dimension (LDA,N)
          The block diagonal matrix D and the multipliers used to
          obtain the factor U or L as computed by DSYTRF.
[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 DSYTRF.
[out]E
          E is DOUBLE PRECISION 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 dsyconv.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  DOUBLE PRECISION a( lda, * ), e( * )
129 * ..
130 *
131 * =====================================================================
132 *
133 * .. Parameters ..
134  DOUBLE PRECISION zero
135  parameter( zero = 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  DOUBLE PRECISION 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( 'DSYCONV', -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 * Convert A (A is upper)
178 *
179 * Convert VALUE
180 *
181  IF ( convert ) THEN
182  i=n
183  e(1)=zero
184  DO WHILE ( i .GT. 1 )
185  IF( ipiv(i) .LT. 0 ) THEN
186  e(i)=a(i-1,i)
187  e(i-1)=zero
188  a(i-1,i)=zero
189  i=i-1
190  ELSE
191  e(i)=zero
192  ENDIF
193  i=i-1
194  END DO
195 *
196 * Convert PERMUTATIONS
197 *
198  i=n
199  DO WHILE ( i .GE. 1 )
200  IF( ipiv(i) .GT. 0) THEN
201  ip=ipiv(i)
202  IF( i .LT. n) THEN
203  DO 12 j= i+1,n
204  temp=a(ip,j)
205  a(ip,j)=a(i,j)
206  a(i,j)=temp
207  12 CONTINUE
208  ENDIF
209  ELSE
210  ip=-ipiv(i)
211  IF( i .LT. n) THEN
212  DO 13 j= i+1,n
213  temp=a(ip,j)
214  a(ip,j)=a(i-1,j)
215  a(i-1,j)=temp
216  13 CONTINUE
217  ENDIF
218  i=i-1
219  ENDIF
220  i=i-1
221  END DO
222 
223  ELSE
224 *
225 * Revert A (A is upper)
226 *
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  ELSE
267 *
268 * A is LOWER
269 *
270  IF ( convert ) THEN
271 *
272 * Convert A (A is lower)
273 *
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  ELSE
318 *
319 * Revert A (A is lower)
320 *
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 DSYCONV
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 dsyequb ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  S,
double precision  SCOND,
double precision  AMAX,
double precision, dimension( * )  WORK,
integer  INFO 
)

DSYEQUB

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

Purpose:
 DSYEQUB 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 137 of file dsyequb.f.

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

DSYGS2 reduces a symmetric definite generalized eigenproblem to standard form, using the factorization results obtained from spotrf (unblocked algorithm).

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

Purpose:
 DSYGS2 reduces a real symmetric-definite generalized eigenproblem
 to standard form.

 If ITYPE = 1, the problem is A*x = lambda*B*x,
 and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)

 If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
 B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T *A*L.

 B must have been previously factorized as U**T *U or L*L**T by DPOTRF.
Parameters
[in]ITYPE
          ITYPE is INTEGER
          = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);
          = 2 or 3: compute U*A*U**T or L**T *A*L.
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the upper or lower triangular part of the
          symmetric matrix A is stored, and how B has been factorized.
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[in]N
          N is INTEGER
          The order of the matrices A and B.  N >= 0.
[in,out]A
          A is DOUBLE PRECISION 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, if INFO = 0, the transformed matrix, stored in the
          same format as A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]B
          B is DOUBLE PRECISION array, dimension (LDB,N)
          The triangular factor from the Cholesky factorization of B,
          as returned by DPOTRF.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,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
September 2012

Definition at line 129 of file dsygs2.f.

129 *
130 * -- LAPACK computational routine (version 3.4.2) --
131 * -- LAPACK is a software package provided by Univ. of Tennessee, --
132 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
133 * September 2012
134 *
135 * .. Scalar Arguments ..
136  CHARACTER uplo
137  INTEGER info, itype, lda, ldb, n
138 * ..
139 * .. Array Arguments ..
140  DOUBLE PRECISION a( lda, * ), b( ldb, * )
141 * ..
142 *
143 * =====================================================================
144 *
145 * .. Parameters ..
146  DOUBLE PRECISION one, half
147  parameter( one = 1.0d0, half = 0.5d0 )
148 * ..
149 * .. Local Scalars ..
150  LOGICAL upper
151  INTEGER k
152  DOUBLE PRECISION akk, bkk, ct
153 * ..
154 * .. External Subroutines ..
155  EXTERNAL daxpy, dscal, dsyr2, dtrmv, dtrsv, xerbla
156 * ..
157 * .. Intrinsic Functions ..
158  INTRINSIC max
159 * ..
160 * .. External Functions ..
161  LOGICAL lsame
162  EXTERNAL lsame
163 * ..
164 * .. Executable Statements ..
165 *
166 * Test the input parameters.
167 *
168  info = 0
169  upper = lsame( uplo, 'U' )
170  IF( itype.LT.1 .OR. itype.GT.3 ) THEN
171  info = -1
172  ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
173  info = -2
174  ELSE IF( n.LT.0 ) THEN
175  info = -3
176  ELSE IF( lda.LT.max( 1, n ) ) THEN
177  info = -5
178  ELSE IF( ldb.LT.max( 1, n ) ) THEN
179  info = -7
180  END IF
181  IF( info.NE.0 ) THEN
182  CALL xerbla( 'DSYGS2', -info )
183  RETURN
184  END IF
185 *
186  IF( itype.EQ.1 ) THEN
187  IF( upper ) THEN
188 *
189 * Compute inv(U**T)*A*inv(U)
190 *
191  DO 10 k = 1, n
192 *
193 * Update the upper triangle of A(k:n,k:n)
194 *
195  akk = a( k, k )
196  bkk = b( k, k )
197  akk = akk / bkk**2
198  a( k, k ) = akk
199  IF( k.LT.n ) THEN
200  CALL dscal( n-k, one / bkk, a( k, k+1 ), lda )
201  ct = -half*akk
202  CALL daxpy( n-k, ct, b( k, k+1 ), ldb, a( k, k+1 ),
203  $ lda )
204  CALL dsyr2( uplo, n-k, -one, a( k, k+1 ), lda,
205  $ b( k, k+1 ), ldb, a( k+1, k+1 ), lda )
206  CALL daxpy( n-k, ct, b( k, k+1 ), ldb, a( k, k+1 ),
207  $ lda )
208  CALL dtrsv( uplo, 'Transpose', 'Non-unit', n-k,
209  $ b( k+1, k+1 ), ldb, a( k, k+1 ), lda )
210  END IF
211  10 CONTINUE
212  ELSE
213 *
214 * Compute inv(L)*A*inv(L**T)
215 *
216  DO 20 k = 1, n
217 *
218 * Update the lower triangle of A(k:n,k:n)
219 *
220  akk = a( k, k )
221  bkk = b( k, k )
222  akk = akk / bkk**2
223  a( k, k ) = akk
224  IF( k.LT.n ) THEN
225  CALL dscal( n-k, one / bkk, a( k+1, k ), 1 )
226  ct = -half*akk
227  CALL daxpy( n-k, ct, b( k+1, k ), 1, a( k+1, k ), 1 )
228  CALL dsyr2( uplo, n-k, -one, a( k+1, k ), 1,
229  $ b( k+1, k ), 1, a( k+1, k+1 ), lda )
230  CALL daxpy( n-k, ct, b( k+1, k ), 1, a( k+1, k ), 1 )
231  CALL dtrsv( uplo, 'No transpose', 'Non-unit', n-k,
232  $ b( k+1, k+1 ), ldb, a( k+1, k ), 1 )
233  END IF
234  20 CONTINUE
235  END IF
236  ELSE
237  IF( upper ) THEN
238 *
239 * Compute U*A*U**T
240 *
241  DO 30 k = 1, n
242 *
243 * Update the upper triangle of A(1:k,1:k)
244 *
245  akk = a( k, k )
246  bkk = b( k, k )
247  CALL dtrmv( uplo, 'No transpose', 'Non-unit', k-1, b,
248  $ ldb, a( 1, k ), 1 )
249  ct = half*akk
250  CALL daxpy( k-1, ct, b( 1, k ), 1, a( 1, k ), 1 )
251  CALL dsyr2( uplo, k-1, one, a( 1, k ), 1, b( 1, k ), 1,
252  $ a, lda )
253  CALL daxpy( k-1, ct, b( 1, k ), 1, a( 1, k ), 1 )
254  CALL dscal( k-1, bkk, a( 1, k ), 1 )
255  a( k, k ) = akk*bkk**2
256  30 CONTINUE
257  ELSE
258 *
259 * Compute L**T *A*L
260 *
261  DO 40 k = 1, n
262 *
263 * Update the lower triangle of A(1:k,1:k)
264 *
265  akk = a( k, k )
266  bkk = b( k, k )
267  CALL dtrmv( uplo, 'Transpose', 'Non-unit', k-1, b, ldb,
268  $ a( k, 1 ), lda )
269  ct = half*akk
270  CALL daxpy( k-1, ct, b( k, 1 ), ldb, a( k, 1 ), lda )
271  CALL dsyr2( uplo, k-1, one, a( k, 1 ), lda, b( k, 1 ),
272  $ ldb, a, lda )
273  CALL daxpy( k-1, ct, b( k, 1 ), ldb, a( k, 1 ), lda )
274  CALL dscal( k-1, bkk, a( k, 1 ), lda )
275  a( k, k ) = akk*bkk**2
276  40 CONTINUE
277  END IF
278  END IF
279  RETURN
280 *
281 * End of DSYGS2
282 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dsyr2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DSYR2
Definition: dsyr2.f:149
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dtrsv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
DTRSV
Definition: dtrsv.f:145
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
subroutine dtrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
DTRMV
Definition: dtrmv.f:149
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:54

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dsygst ( integer  ITYPE,
character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)

DSYGST

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

Purpose:
 DSYGST reduces a real symmetric-definite generalized eigenproblem
 to standard form.

 If ITYPE = 1, the problem is A*x = lambda*B*x,
 and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)

 If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
 B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.

 B must have been previously factorized as U**T*U or L*L**T by DPOTRF.
Parameters
[in]ITYPE
          ITYPE is INTEGER
          = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);
          = 2 or 3: compute U*A*U**T or L**T*A*L.
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored and B is factored as
                  U**T*U;
          = 'L':  Lower triangle of A is stored and B is factored as
                  L*L**T.
[in]N
          N is INTEGER
          The order of the matrices A and B.  N >= 0.
[in,out]A
          A is DOUBLE PRECISION 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, if INFO = 0, the transformed matrix, stored in the
          same format as A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]B
          B is DOUBLE PRECISION array, dimension (LDB,N)
          The triangular factor from the Cholesky factorization of B,
          as returned by DPOTRF.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,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 129 of file dsygst.f.

129 *
130 * -- LAPACK computational routine (version 3.4.0) --
131 * -- LAPACK is a software package provided by Univ. of Tennessee, --
132 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
133 * November 2011
134 *
135 * .. Scalar Arguments ..
136  CHARACTER uplo
137  INTEGER info, itype, lda, ldb, n
138 * ..
139 * .. Array Arguments ..
140  DOUBLE PRECISION a( lda, * ), b( ldb, * )
141 * ..
142 *
143 * =====================================================================
144 *
145 * .. Parameters ..
146  DOUBLE PRECISION one, half
147  parameter( one = 1.0d0, half = 0.5d0 )
148 * ..
149 * .. Local Scalars ..
150  LOGICAL upper
151  INTEGER k, kb, nb
152 * ..
153 * .. External Subroutines ..
154  EXTERNAL dsygs2, dsymm, dsyr2k, dtrmm, dtrsm, xerbla
155 * ..
156 * .. Intrinsic Functions ..
157  INTRINSIC max, min
158 * ..
159 * .. External Functions ..
160  LOGICAL lsame
161  INTEGER ilaenv
162  EXTERNAL lsame, ilaenv
163 * ..
164 * .. Executable Statements ..
165 *
166 * Test the input parameters.
167 *
168  info = 0
169  upper = lsame( uplo, 'U' )
170  IF( itype.LT.1 .OR. itype.GT.3 ) THEN
171  info = -1
172  ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
173  info = -2
174  ELSE IF( n.LT.0 ) THEN
175  info = -3
176  ELSE IF( lda.LT.max( 1, n ) ) THEN
177  info = -5
178  ELSE IF( ldb.LT.max( 1, n ) ) THEN
179  info = -7
180  END IF
181  IF( info.NE.0 ) THEN
182  CALL xerbla( 'DSYGST', -info )
183  RETURN
184  END IF
185 *
186 * Quick return if possible
187 *
188  IF( n.EQ.0 )
189  $ RETURN
190 *
191 * Determine the block size for this environment.
192 *
193  nb = ilaenv( 1, 'DSYGST', uplo, n, -1, -1, -1 )
194 *
195  IF( nb.LE.1 .OR. nb.GE.n ) THEN
196 *
197 * Use unblocked code
198 *
199  CALL dsygs2( itype, uplo, n, a, lda, b, ldb, info )
200  ELSE
201 *
202 * Use blocked code
203 *
204  IF( itype.EQ.1 ) THEN
205  IF( upper ) THEN
206 *
207 * Compute inv(U**T)*A*inv(U)
208 *
209  DO 10 k = 1, n, nb
210  kb = min( n-k+1, nb )
211 *
212 * Update the upper triangle of A(k:n,k:n)
213 *
214  CALL dsygs2( itype, uplo, kb, a( k, k ), lda,
215  $ b( k, k ), ldb, info )
216  IF( k+kb.LE.n ) THEN
217  CALL dtrsm( 'Left', uplo, 'Transpose', 'Non-unit',
218  $ kb, n-k-kb+1, one, b( k, k ), ldb,
219  $ a( k, k+kb ), lda )
220  CALL dsymm( 'Left', uplo, kb, n-k-kb+1, -half,
221  $ a( k, k ), lda, b( k, k+kb ), ldb, one,
222  $ a( k, k+kb ), lda )
223  CALL dsyr2k( uplo, 'Transpose', n-k-kb+1, kb, -one,
224  $ a( k, k+kb ), lda, b( k, k+kb ), ldb,
225  $ one, a( k+kb, k+kb ), lda )
226  CALL dsymm( 'Left', uplo, kb, n-k-kb+1, -half,
227  $ a( k, k ), lda, b( k, k+kb ), ldb, one,
228  $ a( k, k+kb ), lda )
229  CALL dtrsm( 'Right', uplo, 'No transpose',
230  $ 'Non-unit', kb, n-k-kb+1, one,
231  $ b( k+kb, k+kb ), ldb, a( k, k+kb ),
232  $ lda )
233  END IF
234  10 CONTINUE
235  ELSE
236 *
237 * Compute inv(L)*A*inv(L**T)
238 *
239  DO 20 k = 1, n, nb
240  kb = min( n-k+1, nb )
241 *
242 * Update the lower triangle of A(k:n,k:n)
243 *
244  CALL dsygs2( itype, uplo, kb, a( k, k ), lda,
245  $ b( k, k ), ldb, info )
246  IF( k+kb.LE.n ) THEN
247  CALL dtrsm( 'Right', uplo, 'Transpose', 'Non-unit',
248  $ n-k-kb+1, kb, one, b( k, k ), ldb,
249  $ a( k+kb, k ), lda )
250  CALL dsymm( 'Right', uplo, n-k-kb+1, kb, -half,
251  $ a( k, k ), lda, b( k+kb, k ), ldb, one,
252  $ a( k+kb, k ), lda )
253  CALL dsyr2k( uplo, 'No transpose', n-k-kb+1, kb,
254  $ -one, a( k+kb, k ), lda, b( k+kb, k ),
255  $ ldb, one, a( k+kb, k+kb ), lda )
256  CALL dsymm( 'Right', uplo, n-k-kb+1, kb, -half,
257  $ a( k, k ), lda, b( k+kb, k ), ldb, one,
258  $ a( k+kb, k ), lda )
259  CALL dtrsm( 'Left', uplo, 'No transpose',
260  $ 'Non-unit', n-k-kb+1, kb, one,
261  $ b( k+kb, k+kb ), ldb, a( k+kb, k ),
262  $ lda )
263  END IF
264  20 CONTINUE
265  END IF
266  ELSE
267  IF( upper ) THEN
268 *
269 * Compute U*A*U**T
270 *
271  DO 30 k = 1, n, nb
272  kb = min( n-k+1, nb )
273 *
274 * Update the upper triangle of A(1:k+kb-1,1:k+kb-1)
275 *
276  CALL dtrmm( 'Left', uplo, 'No transpose', 'Non-unit',
277  $ k-1, kb, one, b, ldb, a( 1, k ), lda )
278  CALL dsymm( 'Right', uplo, k-1, kb, half, a( k, k ),
279  $ lda, b( 1, k ), ldb, one, a( 1, k ), lda )
280  CALL dsyr2k( uplo, 'No transpose', k-1, kb, one,
281  $ a( 1, k ), lda, b( 1, k ), ldb, one, a,
282  $ lda )
283  CALL dsymm( 'Right', uplo, k-1, kb, half, a( k, k ),
284  $ lda, b( 1, k ), ldb, one, a( 1, k ), lda )
285  CALL dtrmm( 'Right', uplo, 'Transpose', 'Non-unit',
286  $ k-1, kb, one, b( k, k ), ldb, a( 1, k ),
287  $ lda )
288  CALL dsygs2( itype, uplo, kb, a( k, k ), lda,
289  $ b( k, k ), ldb, info )
290  30 CONTINUE
291  ELSE
292 *
293 * Compute L**T*A*L
294 *
295  DO 40 k = 1, n, nb
296  kb = min( n-k+1, nb )
297 *
298 * Update the lower triangle of A(1:k+kb-1,1:k+kb-1)
299 *
300  CALL dtrmm( 'Right', uplo, 'No transpose', 'Non-unit',
301  $ kb, k-1, one, b, ldb, a( k, 1 ), lda )
302  CALL dsymm( 'Left', uplo, kb, k-1, half, a( k, k ),
303  $ lda, b( k, 1 ), ldb, one, a( k, 1 ), lda )
304  CALL dsyr2k( uplo, 'Transpose', k-1, kb, one,
305  $ a( k, 1 ), lda, b( k, 1 ), ldb, one, a,
306  $ lda )
307  CALL dsymm( 'Left', uplo, kb, k-1, half, a( k, k ),
308  $ lda, b( k, 1 ), ldb, one, a( k, 1 ), lda )
309  CALL dtrmm( 'Left', uplo, 'Transpose', 'Non-unit', kb,
310  $ k-1, one, b( k, k ), ldb, a( k, 1 ), lda )
311  CALL dsygs2( itype, uplo, kb, a( k, k ), lda,
312  $ b( k, k ), ldb, info )
313  40 CONTINUE
314  END IF
315  END IF
316  END IF
317  RETURN
318 *
319 * End of DSYGST
320 *
subroutine dsymm(SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DSYMM
Definition: dsymm.f:191
subroutine dtrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRMM
Definition: dtrmm.f:179
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dsyr2k(UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DSYR2K
Definition: dsyr2k.f:194
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
Definition: dtrsm.f:183
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 dsygs2(ITYPE, UPLO, N, A, LDA, B, LDB, INFO)
DSYGS2 reduces a symmetric definite generalized eigenproblem to standard form, using the factorizatio...
Definition: dsygs2.f:129

Here is the call graph for this function:

Here is the caller graph for this function:

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

DSYRFS

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

Purpose:
 DSYRFS 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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]B
          B is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by DSYTRS.
          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 DOUBLE PRECISION array, dimension (3*N)
[out]IWORK
          IWORK is INTEGER 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 193 of file dsyrfs.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dsyrfsx ( character  UPLO,
character  EQUED,
integer  N,
integer  NRHS,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldaf, * )  AF,
integer  LDAF,
integer, dimension( * )  IPIV,
double precision, dimension( * )  S,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, 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,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

DSYRFSX

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

Purpose:
    DSYRFSX 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (4*N)
[out]IWORK
          IWORK is INTEGER array, dimension (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 dsyrfsx.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( * ), iwork( * )
418  DOUBLE PRECISION a( lda, * ), af( ldaf, * ), b( ldb, * ),
419  $ x( ldx, * ), work( * )
420  DOUBLE PRECISION s( * ), params( * ), berr( * ),
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, n_norms
452  DOUBLE PRECISION anorm, rcond_tmp
453  DOUBLE PRECISION illrcond_thresh, err_lbnd, cwise_wrong
454  LOGICAL ignore_cwise
455  INTEGER ithresh
456  DOUBLE PRECISION rthresh, unstable_thresh
457 * ..
458 * .. External Subroutines ..
460 * ..
461 * .. Intrinsic Functions ..
462  INTRINSIC max, sqrt
463 * ..
464 * .. External Functions ..
465  EXTERNAL lsame, blas_fpinfo_x, ilatrans, ilaprec
466  EXTERNAL dlamch, dlansy, dla_syrcond
467  DOUBLE PRECISION dlamch, dlansy, dla_syrcond
468  LOGICAL lsame
469  INTEGER blas_fpinfo_x
470  INTEGER ilatrans, ilaprec
471 * ..
472 * .. Executable Statements ..
473 *
474 * Check the input parameters.
475 *
476  info = 0
477  ref_type = int( itref_default )
478  IF ( nparams .GE. la_linrx_itref_i ) THEN
479  IF ( params( la_linrx_itref_i ) .LT. 0.0d+0 ) THEN
480  params( la_linrx_itref_i ) = itref_default
481  ELSE
482  ref_type = params( la_linrx_itref_i )
483  END IF
484  END IF
485 *
486 * Set default parameters.
487 *
488  illrcond_thresh = dble( n )*dlamch( 'Epsilon' )
489  ithresh = int( ithresh_default )
490  rthresh = rthresh_default
491  unstable_thresh = dzthresh_default
492  ignore_cwise = componentwise_default .EQ. 0.0d+0
493 *
494  IF ( nparams.GE.la_linrx_ithresh_i ) THEN
495  IF ( params( la_linrx_ithresh_i ).LT.0.0d+0 ) THEN
496  params( la_linrx_ithresh_i ) = ithresh
497  ELSE
498  ithresh = int( params( la_linrx_ithresh_i ) )
499  END IF
500  END IF
501  IF ( nparams.GE.la_linrx_cwise_i ) THEN
502  IF ( params( la_linrx_cwise_i ).LT.0.0d+0 ) THEN
503  IF ( ignore_cwise ) THEN
504  params( la_linrx_cwise_i ) = 0.0d+0
505  ELSE
506  params( la_linrx_cwise_i ) = 1.0d+0
507  END IF
508  ELSE
509  ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.0d+0
510  END IF
511  END IF
512  IF ( ref_type .EQ. 0 .OR. n_err_bnds .EQ. 0 ) THEN
513  n_norms = 0
514  ELSE IF ( ignore_cwise ) THEN
515  n_norms = 1
516  ELSE
517  n_norms = 2
518  END IF
519 *
520  rcequ = lsame( equed, 'Y' )
521 *
522 * Test input parameters.
523 *
524  IF ( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
525  info = -1
526  ELSE IF( .NOT.rcequ .AND. .NOT.lsame( equed, 'N' ) ) THEN
527  info = -2
528  ELSE IF( n.LT.0 ) THEN
529  info = -3
530  ELSE IF( nrhs.LT.0 ) THEN
531  info = -4
532  ELSE IF( lda.LT.max( 1, n ) ) THEN
533  info = -6
534  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
535  info = -8
536  ELSE IF( ldb.LT.max( 1, n ) ) THEN
537  info = -12
538  ELSE IF( ldx.LT.max( 1, n ) ) THEN
539  info = -14
540  END IF
541  IF( info.NE.0 ) THEN
542  CALL xerbla( 'DSYRFSX', -info )
543  RETURN
544  END IF
545 *
546 * Quick return if possible.
547 *
548  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
549  rcond = 1.0d+0
550  DO j = 1, nrhs
551  berr( j ) = 0.0d+0
552  IF ( n_err_bnds .GE. 1 ) THEN
553  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
554  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
555  END IF
556  IF ( n_err_bnds .GE. 2 ) THEN
557  err_bnds_norm( j, la_linrx_err_i ) = 0.0d+0
558  err_bnds_comp( j, la_linrx_err_i ) = 0.0d+0
559  END IF
560  IF ( n_err_bnds .GE. 3 ) THEN
561  err_bnds_norm( j, la_linrx_rcond_i ) = 1.0d+0
562  err_bnds_comp( j, la_linrx_rcond_i ) = 1.0d+0
563  END IF
564  END DO
565  RETURN
566  END IF
567 *
568 * Default to failure.
569 *
570  rcond = 0.0d+0
571  DO j = 1, nrhs
572  berr( j ) = 1.0d+0
573  IF ( n_err_bnds .GE. 1 ) THEN
574  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
575  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
576  END IF
577  IF ( n_err_bnds .GE. 2 ) THEN
578  err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
579  err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
580  END IF
581  IF ( n_err_bnds .GE. 3 ) THEN
582  err_bnds_norm( j, la_linrx_rcond_i ) = 0.0d+0
583  err_bnds_comp( j, la_linrx_rcond_i ) = 0.0d+0
584  END IF
585  END DO
586 *
587 * Compute the norm of A and the reciprocal of the condition
588 * number of A.
589 *
590  norm = 'I'
591  anorm = dlansy( norm, uplo, n, a, lda, work )
592  CALL dsycon( uplo, n, af, ldaf, ipiv, anorm, rcond, work,
593  $ iwork, info )
594 *
595 * Perform refinement on each right-hand side
596 *
597  IF ( ref_type .NE. 0 ) THEN
598 
599  prec_type = ilaprec( 'E' )
600 
601  CALL dla_syrfsx_extended( prec_type, uplo, n,
602  $ nrhs, a, lda, af, ldaf, ipiv, rcequ, s, b,
603  $ ldb, x, ldx, berr, n_norms, err_bnds_norm, err_bnds_comp,
604  $ work( n+1 ), work( 1 ), work( 2*n+1 ), work( 1 ), rcond,
605  $ ithresh, rthresh, unstable_thresh, ignore_cwise,
606  $ info )
607  END IF
608 
609  err_lbnd = max( 10.0d+0, sqrt( dble( n ) ) )*dlamch( 'Epsilon' )
610  IF (n_err_bnds .GE. 1 .AND. n_norms .GE. 1) THEN
611 *
612 * Compute scaled normwise condition number cond(A*C).
613 *
614  IF ( rcequ ) THEN
615  rcond_tmp = dla_syrcond( uplo, n, a, lda, af, ldaf, ipiv,
616  $ -1, s, info, work, iwork )
617  ELSE
618  rcond_tmp = dla_syrcond( uplo, n, a, lda, af, ldaf, ipiv,
619  $ 0, s, info, work, iwork )
620  END IF
621  DO j = 1, nrhs
622 *
623 * Cap the error at 1.0.
624 *
625  IF (n_err_bnds .GE. la_linrx_err_i
626  $ .AND. err_bnds_norm( j, la_linrx_err_i ) .GT. 1.0d+0)
627  $ err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
628 *
629 * Threshold the error (see LAWN).
630 *
631  IF ( rcond_tmp .LT. illrcond_thresh ) THEN
632  err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
633  err_bnds_norm( j, la_linrx_trust_i ) = 0.0d+0
634  IF ( info .LE. n ) info = n + j
635  ELSE IF (err_bnds_norm( j, la_linrx_err_i ) .LT. err_lbnd)
636  $ THEN
637  err_bnds_norm( j, la_linrx_err_i ) = err_lbnd
638  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
639  END IF
640 *
641 * Save the condition number.
642 *
643  IF (n_err_bnds .GE. la_linrx_rcond_i) THEN
644  err_bnds_norm( j, la_linrx_rcond_i ) = rcond_tmp
645  END IF
646  END DO
647  END IF
648 
649  IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 2 ) THEN
650 *
651 * Compute componentwise condition number cond(A*diag(Y(:,J))) for
652 * each right-hand side using the current solution as an estimate of
653 * the true solution. If the componentwise error estimate is too
654 * large, then the solution is a lousy estimate of truth and the
655 * estimated RCOND may be too optimistic. To avoid misleading users,
656 * the inverse condition number is set to 0.0 when the estimated
657 * cwise error is at least CWISE_WRONG.
658 *
659  cwise_wrong = sqrt( dlamch( 'Epsilon' ) )
660  DO j = 1, nrhs
661  IF ( err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
662  $ THEN
663  rcond_tmp = dla_syrcond( uplo, n, a, lda, af, ldaf, ipiv,
664  $ 1, x(1,j), info, work, iwork )
665  ELSE
666  rcond_tmp = 0.0d+0
667  END IF
668 *
669 * Cap the error at 1.0.
670 *
671  IF ( n_err_bnds .GE. la_linrx_err_i
672  $ .AND. err_bnds_comp( j, la_linrx_err_i ) .GT. 1.0d+0 )
673  $ err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
674 *
675 * Threshold the error (see LAWN).
676 *
677  IF ( rcond_tmp .LT. illrcond_thresh ) THEN
678  err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
679  err_bnds_comp( j, la_linrx_trust_i ) = 0.0d+0
680  IF ( .NOT. ignore_cwise
681  $ .AND. info.LT.n + j ) info = n + j
682  ELSE IF ( err_bnds_comp( j, la_linrx_err_i )
683  $ .LT. err_lbnd ) THEN
684  err_bnds_comp( j, la_linrx_err_i ) = err_lbnd
685  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
686  END IF
687 *
688 * Save the condition number.
689 *
690  IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
691  err_bnds_comp( j, la_linrx_rcond_i ) = rcond_tmp
692  END IF
693 
694  END DO
695  END IF
696 *
697  RETURN
698 *
699 * End of DSYRFSX
700 *
integer function ilatrans(TRANS)
ILATRANS
Definition: ilatrans.f:60
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dla_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)
DLA_SYRFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric inde...
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 dsycon(UPLO, N, A, LDA, IPIV, ANORM, RCOND, WORK, IWORK, INFO)
DSYCON
Definition: dsycon.f:132
double precision function dla_syrcond(UPLO, N, A, LDA, AF, LDAF, IPIV, CMODE, C, INFO, WORK, IWORK)
DLA_SYRCOND estimates the Skeel condition number for a symmetric indefinite matrix.
Definition: dla_syrcond.f:150
double precision function dlansy(NORM, UPLO, N, A, LDA, WORK)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric matrix.
Definition: dlansy.f:124

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dsytd2 ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( * )  TAU,
integer  INFO 
)

DSYTD2 reduces a symmetric matrix to real symmetric tridiagonal form by an orthogonal similarity transformation (unblocked algorithm).

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

Purpose:
 DSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal
 form T by an orthogonal similarity transformation: Q**T * A * Q = T.
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 DOUBLE PRECISION 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, if UPLO = 'U', the diagonal and first superdiagonal
          of A are overwritten by the corresponding elements of the
          tridiagonal matrix T, and the elements above the first
          superdiagonal, with the array TAU, represent the orthogonal
          matrix Q as a product of elementary reflectors; if UPLO
          = 'L', the diagonal and first subdiagonal of A are over-
          written by the corresponding elements of the tridiagonal
          matrix T, and the elements below the first subdiagonal, with
          the array TAU, represent the orthogonal matrix Q as a product
          of elementary reflectors. See Further Details.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]D
          D is DOUBLE PRECISION array, dimension (N)
          The diagonal elements of the tridiagonal matrix T:
          D(i) = A(i,i).
[out]E
          E is DOUBLE PRECISION array, dimension (N-1)
          The off-diagonal elements of the tridiagonal matrix T:
          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
[out]TAU
          TAU is DOUBLE PRECISION array, dimension (N-1)
          The scalar factors of the elementary reflectors (see Further
          Details).
[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
September 2012
Further Details:
  If UPLO = 'U', the matrix Q is represented as a product of elementary
  reflectors

     Q = H(n-1) . . . H(2) H(1).

  Each H(i) has the form

     H(i) = I - tau * v * v**T

  where tau is a real scalar, and v is a real vector with
  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
  A(1:i-1,i+1), and tau in TAU(i).

  If UPLO = 'L', the matrix Q is represented as a product of elementary
  reflectors

     Q = H(1) H(2) . . . H(n-1).

  Each H(i) has the form

     H(i) = I - tau * v * v**T

  where tau is a real scalar, and v is a real vector with
  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
  and tau in TAU(i).

  The contents of A on exit are illustrated by the following examples
  with n = 5:

  if UPLO = 'U':                       if UPLO = 'L':

    (  d   e   v2  v3  v4 )              (  d                  )
    (      d   e   v3  v4 )              (  e   d              )
    (          d   e   v4 )              (  v1  e   d          )
    (              d   e  )              (  v1  v2  e   d      )
    (                  d  )              (  v1  v2  v3  e   d  )

  where d and e denote diagonal and off-diagonal elements of T, and vi
  denotes an element of the vector defining H(i).

Definition at line 175 of file dsytd2.f.

175 *
176 * -- LAPACK computational routine (version 3.4.2) --
177 * -- LAPACK is a software package provided by Univ. of Tennessee, --
178 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
179 * September 2012
180 *
181 * .. Scalar Arguments ..
182  CHARACTER uplo
183  INTEGER info, lda, n
184 * ..
185 * .. Array Arguments ..
186  DOUBLE PRECISION a( lda, * ), d( * ), e( * ), tau( * )
187 * ..
188 *
189 * =====================================================================
190 *
191 * .. Parameters ..
192  DOUBLE PRECISION one, zero, half
193  parameter( one = 1.0d0, zero = 0.0d0,
194  $ half = 1.0d0 / 2.0d0 )
195 * ..
196 * .. Local Scalars ..
197  LOGICAL upper
198  INTEGER i
199  DOUBLE PRECISION alpha, taui
200 * ..
201 * .. External Subroutines ..
202  EXTERNAL daxpy, dlarfg, dsymv, dsyr2, xerbla
203 * ..
204 * .. External Functions ..
205  LOGICAL lsame
206  DOUBLE PRECISION ddot
207  EXTERNAL lsame, ddot
208 * ..
209 * .. Intrinsic Functions ..
210  INTRINSIC max, min
211 * ..
212 * .. Executable Statements ..
213 *
214 * Test the input parameters
215 *
216  info = 0
217  upper = lsame( uplo, 'U' )
218  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
219  info = -1
220  ELSE IF( n.LT.0 ) THEN
221  info = -2
222  ELSE IF( lda.LT.max( 1, n ) ) THEN
223  info = -4
224  END IF
225  IF( info.NE.0 ) THEN
226  CALL xerbla( 'DSYTD2', -info )
227  RETURN
228  END IF
229 *
230 * Quick return if possible
231 *
232  IF( n.LE.0 )
233  $ RETURN
234 *
235  IF( upper ) THEN
236 *
237 * Reduce the upper triangle of A
238 *
239  DO 10 i = n - 1, 1, -1
240 *
241 * Generate elementary reflector H(i) = I - tau * v * v**T
242 * to annihilate A(1:i-1,i+1)
243 *
244  CALL dlarfg( i, a( i, i+1 ), a( 1, i+1 ), 1, taui )
245  e( i ) = a( i, i+1 )
246 *
247  IF( taui.NE.zero ) THEN
248 *
249 * Apply H(i) from both sides to A(1:i,1:i)
250 *
251  a( i, i+1 ) = one
252 *
253 * Compute x := tau * A * v storing x in TAU(1:i)
254 *
255  CALL dsymv( uplo, i, taui, a, lda, a( 1, i+1 ), 1, zero,
256  $ tau, 1 )
257 *
258 * Compute w := x - 1/2 * tau * (x**T * v) * v
259 *
260  alpha = -half*taui*ddot( i, tau, 1, a( 1, i+1 ), 1 )
261  CALL daxpy( i, alpha, a( 1, i+1 ), 1, tau, 1 )
262 *
263 * Apply the transformation as a rank-2 update:
264 * A := A - v * w**T - w * v**T
265 *
266  CALL dsyr2( uplo, i, -one, a( 1, i+1 ), 1, tau, 1, a,
267  $ lda )
268 *
269  a( i, i+1 ) = e( i )
270  END IF
271  d( i+1 ) = a( i+1, i+1 )
272  tau( i ) = taui
273  10 CONTINUE
274  d( 1 ) = a( 1, 1 )
275  ELSE
276 *
277 * Reduce the lower triangle of A
278 *
279  DO 20 i = 1, n - 1
280 *
281 * Generate elementary reflector H(i) = I - tau * v * v**T
282 * to annihilate A(i+2:n,i)
283 *
284  CALL dlarfg( n-i, a( i+1, i ), a( min( i+2, n ), i ), 1,
285  $ taui )
286  e( i ) = a( i+1, i )
287 *
288  IF( taui.NE.zero ) THEN
289 *
290 * Apply H(i) from both sides to A(i+1:n,i+1:n)
291 *
292  a( i+1, i ) = one
293 *
294 * Compute x := tau * A * v storing y in TAU(i:n-1)
295 *
296  CALL dsymv( uplo, n-i, taui, a( i+1, i+1 ), lda,
297  $ a( i+1, i ), 1, zero, tau( i ), 1 )
298 *
299 * Compute w := x - 1/2 * tau * (x**T * v) * v
300 *
301  alpha = -half*taui*ddot( n-i, tau( i ), 1, a( i+1, i ),
302  $ 1 )
303  CALL daxpy( n-i, alpha, a( i+1, i ), 1, tau( i ), 1 )
304 *
305 * Apply the transformation as a rank-2 update:
306 * A := A - v * w**T - w * v**T
307 *
308  CALL dsyr2( uplo, n-i, -one, a( i+1, i ), 1, tau( i ), 1,
309  $ a( i+1, i+1 ), lda )
310 *
311  a( i+1, i ) = e( i )
312  END IF
313  d( i ) = a( i, i )
314  tau( i ) = taui
315  20 CONTINUE
316  d( n ) = a( n, n )
317  END IF
318 *
319  RETURN
320 *
321 * End of DSYTD2
322 *
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:108
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dsyr2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DSYR2
Definition: dsyr2.f:149
subroutine dsymv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DSYMV
Definition: dsymv.f:154
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:53
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:54

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dsytf2 ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
integer  INFO 
)

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

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

Purpose:
 DSYTF2 computes the factorization of a real 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 DOUBLE PRECISION 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.204 and l.372
         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
    by
         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN

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

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dsytf2_rook ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
integer  INFO 
)

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

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

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dsytrd ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DSYTRD

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

Purpose:
 DSYTRD reduces a real symmetric matrix A to real symmetric
 tridiagonal form T by an orthogonal similarity transformation:
 Q**T * A * Q = T.
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 DOUBLE PRECISION 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, if UPLO = 'U', the diagonal and first superdiagonal
          of A are overwritten by the corresponding elements of the
          tridiagonal matrix T, and the elements above the first
          superdiagonal, with the array TAU, represent the orthogonal
          matrix Q as a product of elementary reflectors; if UPLO
          = 'L', the diagonal and first subdiagonal of A are over-
          written by the corresponding elements of the tridiagonal
          matrix T, and the elements below the first subdiagonal, with
          the array TAU, represent the orthogonal matrix Q as a product
          of elementary reflectors. See Further Details.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]D
          D is DOUBLE PRECISION array, dimension (N)
          The diagonal elements of the tridiagonal matrix T:
          D(i) = A(i,i).
[out]E
          E is DOUBLE PRECISION array, dimension (N-1)
          The off-diagonal elements of the tridiagonal matrix T:
          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
[out]TAU
          TAU is DOUBLE PRECISION array, dimension (N-1)
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= 1.
          For optimum performance LWORK >= N*NB, where NB is the
          optimal blocksize.

          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
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Further Details:
  If UPLO = 'U', the matrix Q is represented as a product of elementary
  reflectors

     Q = H(n-1) . . . H(2) H(1).

  Each H(i) has the form

     H(i) = I - tau * v * v**T

  where tau is a real scalar, and v is a real vector with
  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
  A(1:i-1,i+1), and tau in TAU(i).

  If UPLO = 'L', the matrix Q is represented as a product of elementary
  reflectors

     Q = H(1) H(2) . . . H(n-1).

  Each H(i) has the form

     H(i) = I - tau * v * v**T

  where tau is a real scalar, and v is a real vector with
  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
  and tau in TAU(i).

  The contents of A on exit are illustrated by the following examples
  with n = 5:

  if UPLO = 'U':                       if UPLO = 'L':

    (  d   e   v2  v3  v4 )              (  d                  )
    (      d   e   v3  v4 )              (  e   d              )
    (          d   e   v4 )              (  v1  e   d          )
    (              d   e  )              (  v1  v2  e   d      )
    (                  d  )              (  v1  v2  v3  e   d  )

  where d and e denote diagonal and off-diagonal elements of T, and vi
  denotes an element of the vector defining H(i).

Definition at line 194 of file dsytrd.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, lwork, n
203 * ..
204 * .. Array Arguments ..
205  DOUBLE PRECISION a( lda, * ), d( * ), e( * ), tau( * ),
206  $ work( * )
207 * ..
208 *
209 * =====================================================================
210 *
211 * .. Parameters ..
212  DOUBLE PRECISION one
213  parameter( one = 1.0d+0 )
214 * ..
215 * .. Local Scalars ..
216  LOGICAL lquery, upper
217  INTEGER i, iinfo, iws, j, kk, ldwork, lwkopt, nb,
218  $ nbmin, nx
219 * ..
220 * .. External Subroutines ..
221  EXTERNAL dlatrd, dsyr2k, dsytd2, xerbla
222 * ..
223 * .. Intrinsic Functions ..
224  INTRINSIC max
225 * ..
226 * .. External Functions ..
227  LOGICAL lsame
228  INTEGER ilaenv
229  EXTERNAL lsame, ilaenv
230 * ..
231 * .. Executable Statements ..
232 *
233 * Test the input parameters
234 *
235  info = 0
236  upper = lsame( uplo, 'U' )
237  lquery = ( lwork.EQ.-1 )
238  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
239  info = -1
240  ELSE IF( n.LT.0 ) THEN
241  info = -2
242  ELSE IF( lda.LT.max( 1, n ) ) THEN
243  info = -4
244  ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
245  info = -9
246  END IF
247 *
248  IF( info.EQ.0 ) THEN
249 *
250 * Determine the block size.
251 *
252  nb = ilaenv( 1, 'DSYTRD', uplo, n, -1, -1, -1 )
253  lwkopt = n*nb
254  work( 1 ) = lwkopt
255  END IF
256 *
257  IF( info.NE.0 ) THEN
258  CALL xerbla( 'DSYTRD', -info )
259  RETURN
260  ELSE IF( lquery ) THEN
261  RETURN
262  END IF
263 *
264 * Quick return if possible
265 *
266  IF( n.EQ.0 ) THEN
267  work( 1 ) = 1
268  RETURN
269  END IF
270 *
271  nx = n
272  iws = 1
273  IF( nb.GT.1 .AND. nb.LT.n ) THEN
274 *
275 * Determine when to cross over from blocked to unblocked code
276 * (last block is always handled by unblocked code).
277 *
278  nx = max( nb, ilaenv( 3, 'DSYTRD', uplo, n, -1, -1, -1 ) )
279  IF( nx.LT.n ) THEN
280 *
281 * Determine if workspace is large enough for blocked code.
282 *
283  ldwork = n
284  iws = ldwork*nb
285  IF( lwork.LT.iws ) THEN
286 *
287 * Not enough workspace to use optimal NB: determine the
288 * minimum value of NB, and reduce NB or force use of
289 * unblocked code by setting NX = N.
290 *
291  nb = max( lwork / ldwork, 1 )
292  nbmin = ilaenv( 2, 'DSYTRD', uplo, n, -1, -1, -1 )
293  IF( nb.LT.nbmin )
294  $ nx = n
295  END IF
296  ELSE
297  nx = n
298  END IF
299  ELSE
300  nb = 1
301  END IF
302 *
303  IF( upper ) THEN
304 *
305 * Reduce the upper triangle of A.
306 * Columns 1:kk are handled by the unblocked method.
307 *
308  kk = n - ( ( n-nx+nb-1 ) / nb )*nb
309  DO 20 i = n - nb + 1, kk + 1, -nb
310 *
311 * Reduce columns i:i+nb-1 to tridiagonal form and form the
312 * matrix W which is needed to update the unreduced part of
313 * the matrix
314 *
315  CALL dlatrd( uplo, i+nb-1, nb, a, lda, e, tau, work,
316  $ ldwork )
317 *
318 * Update the unreduced submatrix A(1:i-1,1:i-1), using an
319 * update of the form: A := A - V*W**T - W*V**T
320 *
321  CALL dsyr2k( uplo, 'No transpose', i-1, nb, -one, a( 1, i ),
322  $ lda, work, ldwork, one, a, lda )
323 *
324 * Copy superdiagonal elements back into A, and diagonal
325 * elements into D
326 *
327  DO 10 j = i, i + nb - 1
328  a( j-1, j ) = e( j-1 )
329  d( j ) = a( j, j )
330  10 CONTINUE
331  20 CONTINUE
332 *
333 * Use unblocked code to reduce the last or only block
334 *
335  CALL dsytd2( uplo, kk, a, lda, d, e, tau, iinfo )
336  ELSE
337 *
338 * Reduce the lower triangle of A
339 *
340  DO 40 i = 1, n - nx, nb
341 *
342 * Reduce columns i:i+nb-1 to tridiagonal form and form the
343 * matrix W which is needed to update the unreduced part of
344 * the matrix
345 *
346  CALL dlatrd( uplo, n-i+1, nb, a( i, i ), lda, e( i ),
347  $ tau( i ), work, ldwork )
348 *
349 * Update the unreduced submatrix A(i+ib:n,i+ib:n), using
350 * an update of the form: A := A - V*W**T - W*V**T
351 *
352  CALL dsyr2k( uplo, 'No transpose', n-i-nb+1, nb, -one,
353  $ a( i+nb, i ), lda, work( nb+1 ), ldwork, one,
354  $ a( i+nb, i+nb ), lda )
355 *
356 * Copy subdiagonal elements back into A, and diagonal
357 * elements into D
358 *
359  DO 30 j = i, i + nb - 1
360  a( j+1, j ) = e( j )
361  d( j ) = a( j, j )
362  30 CONTINUE
363  40 CONTINUE
364 *
365 * Use unblocked code to reduce the last or only block
366 *
367  CALL dsytd2( uplo, n-i+1, a( i, i ), lda, d( i ), e( i ),
368  $ tau( i ), iinfo )
369  END IF
370 *
371  work( 1 ) = lwkopt
372  RETURN
373 *
374 * End of DSYTRD
375 *
subroutine dsytd2(UPLO, N, A, LDA, D, E, TAU, INFO)
DSYTD2 reduces a symmetric matrix to real symmetric tridiagonal form by an orthogonal similarity tran...
Definition: dsytd2.f:175
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dsyr2k(UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DSYR2K
Definition: dsyr2k.f:194
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dlatrd(UPLO, N, NB, A, LDA, E, TAU, W, LDW)
DLATRD reduces the first nb rows and columns of a symmetric/Hermitian matrix A to real tridiagonal fo...
Definition: dlatrd.f:200
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 dsytrf ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DSYTRF

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

Purpose:
 DSYTRF computes the factorization of a real 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
 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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), whe