Annotation of rpl/lapack/lapack/dsysvx.f, revision 1.9
1.9 ! bertrand 1: *> \brief <b> DSYSVX computes the solution to system of linear equations A * X = B for SY matrices</b>
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DSYSVX + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsysvx.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsysvx.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsysvx.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
! 22: * LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
! 23: * IWORK, INFO )
! 24: *
! 25: * .. Scalar Arguments ..
! 26: * CHARACTER FACT, UPLO
! 27: * INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
! 28: * DOUBLE PRECISION RCOND
! 29: * ..
! 30: * .. Array Arguments ..
! 31: * INTEGER IPIV( * ), IWORK( * )
! 32: * DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
! 33: * $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
! 34: * ..
! 35: *
! 36: *
! 37: *> \par Purpose:
! 38: * =============
! 39: *>
! 40: *> \verbatim
! 41: *>
! 42: *> DSYSVX uses the diagonal pivoting factorization to compute the
! 43: *> solution to a real system of linear equations A * X = B,
! 44: *> where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
! 45: *> matrices.
! 46: *>
! 47: *> Error bounds on the solution and a condition estimate are also
! 48: *> provided.
! 49: *> \endverbatim
! 50: *
! 51: *> \par Description:
! 52: * =================
! 53: *>
! 54: *> \verbatim
! 55: *>
! 56: *> The following steps are performed:
! 57: *>
! 58: *> 1. If FACT = 'N', the diagonal pivoting method is used to factor A.
! 59: *> The form of the factorization is
! 60: *> A = U * D * U**T, if UPLO = 'U', or
! 61: *> A = L * D * L**T, if UPLO = 'L',
! 62: *> where U (or L) is a product of permutation and unit upper (lower)
! 63: *> triangular matrices, and D is symmetric and block diagonal with
! 64: *> 1-by-1 and 2-by-2 diagonal blocks.
! 65: *>
! 66: *> 2. If some D(i,i)=0, so that D is exactly singular, then the routine
! 67: *> returns with INFO = i. Otherwise, the factored form of A is used
! 68: *> to estimate the condition number of the matrix A. If the
! 69: *> reciprocal of the condition number is less than machine precision,
! 70: *> INFO = N+1 is returned as a warning, but the routine still goes on
! 71: *> to solve for X and compute error bounds as described below.
! 72: *>
! 73: *> 3. The system of equations is solved for X using the factored form
! 74: *> of A.
! 75: *>
! 76: *> 4. Iterative refinement is applied to improve the computed solution
! 77: *> matrix and calculate error bounds and backward error estimates
! 78: *> for it.
! 79: *> \endverbatim
! 80: *
! 81: * Arguments:
! 82: * ==========
! 83: *
! 84: *> \param[in] FACT
! 85: *> \verbatim
! 86: *> FACT is CHARACTER*1
! 87: *> Specifies whether or not the factored form of A has been
! 88: *> supplied on entry.
! 89: *> = 'F': On entry, AF and IPIV contain the factored form of
! 90: *> A. AF and IPIV will not be modified.
! 91: *> = 'N': The matrix A will be copied to AF and factored.
! 92: *> \endverbatim
! 93: *>
! 94: *> \param[in] UPLO
! 95: *> \verbatim
! 96: *> UPLO is CHARACTER*1
! 97: *> = 'U': Upper triangle of A is stored;
! 98: *> = 'L': Lower triangle of A is stored.
! 99: *> \endverbatim
! 100: *>
! 101: *> \param[in] N
! 102: *> \verbatim
! 103: *> N is INTEGER
! 104: *> The number of linear equations, i.e., the order of the
! 105: *> matrix A. N >= 0.
! 106: *> \endverbatim
! 107: *>
! 108: *> \param[in] NRHS
! 109: *> \verbatim
! 110: *> NRHS is INTEGER
! 111: *> The number of right hand sides, i.e., the number of columns
! 112: *> of the matrices B and X. NRHS >= 0.
! 113: *> \endverbatim
! 114: *>
! 115: *> \param[in] A
! 116: *> \verbatim
! 117: *> A is DOUBLE PRECISION array, dimension (LDA,N)
! 118: *> The symmetric matrix A. If UPLO = 'U', the leading N-by-N
! 119: *> upper triangular part of A contains the upper triangular part
! 120: *> of the matrix A, and the strictly lower triangular part of A
! 121: *> is not referenced. If UPLO = 'L', the leading N-by-N lower
! 122: *> triangular part of A contains the lower triangular part of
! 123: *> the matrix A, and the strictly upper triangular part of A is
! 124: *> not referenced.
! 125: *> \endverbatim
! 126: *>
! 127: *> \param[in] LDA
! 128: *> \verbatim
! 129: *> LDA is INTEGER
! 130: *> The leading dimension of the array A. LDA >= max(1,N).
! 131: *> \endverbatim
! 132: *>
! 133: *> \param[in,out] AF
! 134: *> \verbatim
! 135: *> AF is or output) DOUBLE PRECISION array, dimension (LDAF,N)
! 136: *> If FACT = 'F', then AF is an input argument and on entry
! 137: *> contains the block diagonal matrix D and the multipliers used
! 138: *> to obtain the factor U or L from the factorization
! 139: *> A = U*D*U**T or A = L*D*L**T as computed by DSYTRF.
! 140: *>
! 141: *> If FACT = 'N', then AF is an output argument and on exit
! 142: *> returns the block diagonal matrix D and the multipliers used
! 143: *> to obtain the factor U or L from the factorization
! 144: *> A = U*D*U**T or A = L*D*L**T.
! 145: *> \endverbatim
! 146: *>
! 147: *> \param[in] LDAF
! 148: *> \verbatim
! 149: *> LDAF is INTEGER
! 150: *> The leading dimension of the array AF. LDAF >= max(1,N).
! 151: *> \endverbatim
! 152: *>
! 153: *> \param[in,out] IPIV
! 154: *> \verbatim
! 155: *> IPIV is or output) INTEGER array, dimension (N)
! 156: *> If FACT = 'F', then IPIV is an input argument and on entry
! 157: *> contains details of the interchanges and the block structure
! 158: *> of D, as determined by DSYTRF.
! 159: *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
! 160: *> interchanged and D(k,k) is a 1-by-1 diagonal block.
! 161: *> If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
! 162: *> columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
! 163: *> is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
! 164: *> IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
! 165: *> interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
! 166: *>
! 167: *> If FACT = 'N', then IPIV is an output argument and on exit
! 168: *> contains details of the interchanges and the block structure
! 169: *> of D, as determined by DSYTRF.
! 170: *> \endverbatim
! 171: *>
! 172: *> \param[in] B
! 173: *> \verbatim
! 174: *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
! 175: *> The N-by-NRHS right hand side matrix B.
! 176: *> \endverbatim
! 177: *>
! 178: *> \param[in] LDB
! 179: *> \verbatim
! 180: *> LDB is INTEGER
! 181: *> The leading dimension of the array B. LDB >= max(1,N).
! 182: *> \endverbatim
! 183: *>
! 184: *> \param[out] X
! 185: *> \verbatim
! 186: *> X is DOUBLE PRECISION array, dimension (LDX,NRHS)
! 187: *> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
! 188: *> \endverbatim
! 189: *>
! 190: *> \param[in] LDX
! 191: *> \verbatim
! 192: *> LDX is INTEGER
! 193: *> The leading dimension of the array X. LDX >= max(1,N).
! 194: *> \endverbatim
! 195: *>
! 196: *> \param[out] RCOND
! 197: *> \verbatim
! 198: *> RCOND is DOUBLE PRECISION
! 199: *> The estimate of the reciprocal condition number of the matrix
! 200: *> A. If RCOND is less than the machine precision (in
! 201: *> particular, if RCOND = 0), the matrix is singular to working
! 202: *> precision. This condition is indicated by a return code of
! 203: *> INFO > 0.
! 204: *> \endverbatim
! 205: *>
! 206: *> \param[out] FERR
! 207: *> \verbatim
! 208: *> FERR is DOUBLE PRECISION array, dimension (NRHS)
! 209: *> The estimated forward error bound for each solution vector
! 210: *> X(j) (the j-th column of the solution matrix X).
! 211: *> If XTRUE is the true solution corresponding to X(j), FERR(j)
! 212: *> is an estimated upper bound for the magnitude of the largest
! 213: *> element in (X(j) - XTRUE) divided by the magnitude of the
! 214: *> largest element in X(j). The estimate is as reliable as
! 215: *> the estimate for RCOND, and is almost always a slight
! 216: *> overestimate of the true error.
! 217: *> \endverbatim
! 218: *>
! 219: *> \param[out] BERR
! 220: *> \verbatim
! 221: *> BERR is DOUBLE PRECISION array, dimension (NRHS)
! 222: *> The componentwise relative backward error of each solution
! 223: *> vector X(j) (i.e., the smallest relative change in
! 224: *> any element of A or B that makes X(j) an exact solution).
! 225: *> \endverbatim
! 226: *>
! 227: *> \param[out] WORK
! 228: *> \verbatim
! 229: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 230: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 231: *> \endverbatim
! 232: *>
! 233: *> \param[in] LWORK
! 234: *> \verbatim
! 235: *> LWORK is INTEGER
! 236: *> The length of WORK. LWORK >= max(1,3*N), and for best
! 237: *> performance, when FACT = 'N', LWORK >= max(1,3*N,N*NB), where
! 238: *> NB is the optimal blocksize for DSYTRF.
! 239: *>
! 240: *> If LWORK = -1, then a workspace query is assumed; the routine
! 241: *> only calculates the optimal size of the WORK array, returns
! 242: *> this value as the first entry of the WORK array, and no error
! 243: *> message related to LWORK is issued by XERBLA.
! 244: *> \endverbatim
! 245: *>
! 246: *> \param[out] IWORK
! 247: *> \verbatim
! 248: *> IWORK is INTEGER array, dimension (N)
! 249: *> \endverbatim
! 250: *>
! 251: *> \param[out] INFO
! 252: *> \verbatim
! 253: *> INFO is INTEGER
! 254: *> = 0: successful exit
! 255: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 256: *> > 0: if INFO = i, and i is
! 257: *> <= N: D(i,i) is exactly zero. The factorization
! 258: *> has been completed but the factor D is exactly
! 259: *> singular, so the solution and error bounds could
! 260: *> not be computed. RCOND = 0 is returned.
! 261: *> = N+1: D is nonsingular, but RCOND is less than machine
! 262: *> precision, meaning that the matrix is singular
! 263: *> to working precision. Nevertheless, the
! 264: *> solution and error bounds are computed because
! 265: *> there are a number of situations where the
! 266: *> computed solution can be more accurate than the
! 267: *> value of RCOND would suggest.
! 268: *> \endverbatim
! 269: *
! 270: * Authors:
! 271: * ========
! 272: *
! 273: *> \author Univ. of Tennessee
! 274: *> \author Univ. of California Berkeley
! 275: *> \author Univ. of Colorado Denver
! 276: *> \author NAG Ltd.
! 277: *
! 278: *> \date November 2011
! 279: *
! 280: *> \ingroup doubleSYsolve
! 281: *
! 282: * =====================================================================
1.1 bertrand 283: SUBROUTINE DSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
284: $ LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
285: $ IWORK, INFO )
286: *
1.9 ! bertrand 287: * -- LAPACK driver routine (version 3.4.0) --
1.1 bertrand 288: * -- LAPACK is a software package provided by Univ. of Tennessee, --
289: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 290: * November 2011
1.1 bertrand 291: *
292: * .. Scalar Arguments ..
293: CHARACTER FACT, UPLO
294: INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
295: DOUBLE PRECISION RCOND
296: * ..
297: * .. Array Arguments ..
298: INTEGER IPIV( * ), IWORK( * )
299: DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
300: $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
301: * ..
302: *
303: * =====================================================================
304: *
305: * .. Parameters ..
306: DOUBLE PRECISION ZERO
307: PARAMETER ( ZERO = 0.0D+0 )
308: * ..
309: * .. Local Scalars ..
310: LOGICAL LQUERY, NOFACT
311: INTEGER LWKOPT, NB
312: DOUBLE PRECISION ANORM
313: * ..
314: * .. External Functions ..
315: LOGICAL LSAME
316: INTEGER ILAENV
317: DOUBLE PRECISION DLAMCH, DLANSY
318: EXTERNAL LSAME, ILAENV, DLAMCH, DLANSY
319: * ..
320: * .. External Subroutines ..
321: EXTERNAL DLACPY, DSYCON, DSYRFS, DSYTRF, DSYTRS, XERBLA
322: * ..
323: * .. Intrinsic Functions ..
324: INTRINSIC MAX
325: * ..
326: * .. Executable Statements ..
327: *
328: * Test the input parameters.
329: *
330: INFO = 0
331: NOFACT = LSAME( FACT, 'N' )
332: LQUERY = ( LWORK.EQ.-1 )
333: IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
334: INFO = -1
335: ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
336: $ THEN
337: INFO = -2
338: ELSE IF( N.LT.0 ) THEN
339: INFO = -3
340: ELSE IF( NRHS.LT.0 ) THEN
341: INFO = -4
342: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
343: INFO = -6
344: ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
345: INFO = -8
346: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
347: INFO = -11
348: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
349: INFO = -13
350: ELSE IF( LWORK.LT.MAX( 1, 3*N ) .AND. .NOT.LQUERY ) THEN
351: INFO = -18
352: END IF
353: *
354: IF( INFO.EQ.0 ) THEN
355: LWKOPT = MAX( 1, 3*N )
356: IF( NOFACT ) THEN
357: NB = ILAENV( 1, 'DSYTRF', UPLO, N, -1, -1, -1 )
358: LWKOPT = MAX( LWKOPT, N*NB )
359: END IF
360: WORK( 1 ) = LWKOPT
361: END IF
362: *
363: IF( INFO.NE.0 ) THEN
364: CALL XERBLA( 'DSYSVX', -INFO )
365: RETURN
366: ELSE IF( LQUERY ) THEN
367: RETURN
368: END IF
369: *
370: IF( NOFACT ) THEN
371: *
1.8 bertrand 372: * Compute the factorization A = U*D*U**T or A = L*D*L**T.
1.1 bertrand 373: *
374: CALL DLACPY( UPLO, N, N, A, LDA, AF, LDAF )
375: CALL DSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO )
376: *
377: * Return if INFO is non-zero.
378: *
379: IF( INFO.GT.0 )THEN
380: RCOND = ZERO
381: RETURN
382: END IF
383: END IF
384: *
385: * Compute the norm of the matrix A.
386: *
387: ANORM = DLANSY( 'I', UPLO, N, A, LDA, WORK )
388: *
389: * Compute the reciprocal of the condition number of A.
390: *
391: CALL DSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, IWORK,
392: $ INFO )
393: *
394: * Compute the solution vectors X.
395: *
396: CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
397: CALL DSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
398: *
399: * Use iterative refinement to improve the computed solutions and
400: * compute error bounds and backward error estimates for them.
401: *
402: CALL DSYRFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X,
403: $ LDX, FERR, BERR, WORK, IWORK, INFO )
404: *
405: * Set INFO = N+1 if the matrix is singular to working precision.
406: *
407: IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
408: $ INFO = N + 1
409: *
410: WORK( 1 ) = LWKOPT
411: *
412: RETURN
413: *
414: * End of DSYSVX
415: *
416: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsysvx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack