zspsvx.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 302 行 · 第 1/2 页
HTML
302 行
</span><span class="comment">*</span><span class="comment"> INFO > 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> FERR (output) DOUBLE PRECISION array, dimension (NRHS)
</span><span class="comment">*</span><span class="comment"> The estimated forward error bound for each solution vector
</span><span class="comment">*</span><span class="comment"> X(j) (the j-th column of the solution matrix X).
</span><span class="comment">*</span><span class="comment"> If XTRUE is the true solution corresponding to X(j), FERR(j)
</span><span class="comment">*</span><span class="comment"> is an estimated upper bound for the magnitude of the largest
</span><span class="comment">*</span><span class="comment"> element in (X(j) - XTRUE) divided by the magnitude of the
</span><span class="comment">*</span><span class="comment"> largest element in X(j). The estimate is as reliable as
</span><span class="comment">*</span><span class="comment"> the estimate for RCOND, and is almost always a slight
</span><span class="comment">*</span><span class="comment"> overestimate of the true error.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> BERR (output) DOUBLE PRECISION array, dimension (NRHS)
</span><span class="comment">*</span><span class="comment"> The componentwise relative backward error of each solution
</span><span class="comment">*</span><span class="comment"> vector X(j) (i.e., the smallest relative change in
</span><span class="comment">*</span><span class="comment"> any element of A or B that makes X(j) an exact solution).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> WORK (workspace) COMPLEX*16 array, dimension (2*N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> RWORK (workspace) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> INFO (output) INTEGER
</span><span class="comment">*</span><span class="comment"> = 0: successful exit
</span><span class="comment">*</span><span class="comment"> < 0: if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment"> > 0: if INFO = i, and i is
</span><span class="comment">*</span><span class="comment"> <= N: D(i,i) is exactly zero. The factorization
</span><span class="comment">*</span><span class="comment"> has been completed but the factor D is exactly
</span><span class="comment">*</span><span class="comment"> singular, so the solution and error bounds could
</span><span class="comment">*</span><span class="comment"> not be computed. RCOND = 0 is returned.
</span><span class="comment">*</span><span class="comment"> = N+1: D is nonsingular, but RCOND is less than machine
</span><span class="comment">*</span><span class="comment"> precision, meaning that the matrix is singular
</span><span class="comment">*</span><span class="comment"> to working precision. Nevertheless, the
</span><span class="comment">*</span><span class="comment"> solution and error bounds are computed because
</span><span class="comment">*</span><span class="comment"> there are a number of situations where the
</span><span class="comment">*</span><span class="comment"> computed solution can be more accurate than the
</span><span class="comment">*</span><span class="comment"> value of RCOND would suggest.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Further Details
</span><span class="comment">*</span><span class="comment"> ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The packed storage scheme is illustrated by the following example
</span><span class="comment">*</span><span class="comment"> when N = 4, UPLO = 'U':
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Two-dimensional storage of the symmetric matrix A:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> a11 a12 a13 a14
</span><span class="comment">*</span><span class="comment"> a22 a23 a24
</span><span class="comment">*</span><span class="comment"> a33 a34 (aij = aji)
</span><span class="comment">*</span><span class="comment"> a44
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Packed storage of the upper triangle of A:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. Parameters ..
</span> DOUBLE PRECISION ZERO
PARAMETER ( ZERO = 0.0D+0 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Scalars ..
</span> LOGICAL NOFACT
DOUBLE PRECISION ANORM
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Functions ..
</span> LOGICAL <a name="LSAME.198"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
DOUBLE PRECISION <a name="DLAMCH.199"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="ZLANSP.199"></a><a href="zlansp.f.html#ZLANSP.1">ZLANSP</a>
EXTERNAL <a name="LSAME.200"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="DLAMCH.200"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="ZLANSP.200"></a><a href="zlansp.f.html#ZLANSP.1">ZLANSP</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Subroutines ..
</span> EXTERNAL <a name="XERBLA.203"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>, ZCOPY, <a name="ZLACPY.203"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>, <a name="ZSPCON.203"></a><a href="zspcon.f.html#ZSPCON.1">ZSPCON</a>, <a name="ZSPRFS.203"></a><a href="zsprfs.f.html#ZSPRFS.1">ZSPRFS</a>, <a name="ZSPTRF.203"></a><a href="zsptrf.f.html#ZSPTRF.1">ZSPTRF</a>,
$ <a name="ZSPTRS.204"></a><a href="zsptrs.f.html#ZSPTRS.1">ZSPTRS</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Intrinsic Functions ..
</span> INTRINSIC MAX
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span> INFO = 0
NOFACT = <a name="LSAME.214"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( FACT, <span class="string">'N'</span> )
IF( .NOT.NOFACT .AND. .NOT.<a name="LSAME.215"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( FACT, <span class="string">'F'</span> ) ) THEN
INFO = -1
ELSE IF( .NOT.<a name="LSAME.217"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> ) .AND. .NOT.<a name="LSAME.217"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'L'</span> ) )
$ THEN
INFO = -2
ELSE IF( N.LT.0 ) THEN
INFO = -3
ELSE IF( NRHS.LT.0 ) THEN
INFO = -4
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
INFO = -9
ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
INFO = -11
END IF
IF( INFO.NE.0 ) THEN
CALL <a name="XERBLA.230"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZSPSVX.230"></a><a href="zspsvx.f.html#ZSPSVX.1">ZSPSVX</a>'</span>, -INFO )
RETURN
END IF
<span class="comment">*</span><span class="comment">
</span> IF( NOFACT ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute the factorization A = U*D*U' or A = L*D*L'.
</span><span class="comment">*</span><span class="comment">
</span> CALL ZCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 )
CALL <a name="ZSPTRF.239"></a><a href="zsptrf.f.html#ZSPTRF.1">ZSPTRF</a>( UPLO, N, AFP, IPIV, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Return if INFO is non-zero.
</span><span class="comment">*</span><span class="comment">
</span> IF( INFO.GT.0 )THEN
RCOND = ZERO
RETURN
END IF
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute the norm of the matrix A.
</span><span class="comment">*</span><span class="comment">
</span> ANORM = <a name="ZLANSP.251"></a><a href="zlansp.f.html#ZLANSP.1">ZLANSP</a>( <span class="string">'I'</span>, UPLO, N, AP, RWORK )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute the reciprocal of the condition number of A.
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="ZSPCON.255"></a><a href="zspcon.f.html#ZSPCON.1">ZSPCON</a>( UPLO, N, AFP, IPIV, ANORM, RCOND, WORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute the solution vectors X.
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="ZLACPY.259"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'Full'</span>, N, NRHS, B, LDB, X, LDX )
CALL <a name="ZSPTRS.260"></a><a href="zsptrs.f.html#ZSPTRS.1">ZSPTRS</a>( UPLO, N, NRHS, AFP, IPIV, X, LDX, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Use iterative refinement to improve the computed solutions and
</span><span class="comment">*</span><span class="comment"> compute error bounds and backward error estimates for them.
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="ZSPRFS.265"></a><a href="zsprfs.f.html#ZSPRFS.1">ZSPRFS</a>( UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, LDX, FERR,
$ BERR, WORK, RWORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Set INFO = N+1 if the matrix is singular to working precision.
</span><span class="comment">*</span><span class="comment">
</span> IF( RCOND.LT.<a name="DLAMCH.270"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Epsilon'</span> ) )
$ INFO = N + 1
<span class="comment">*</span><span class="comment">
</span> RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> End of <a name="ZSPSVX.275"></a><a href="zspsvx.f.html#ZSPSVX.1">ZSPSVX</a>
</span><span class="comment">*</span><span class="comment">
</span> END
</pre>
</body>
</html>
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?