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