--- rpl/lapack/lapack/zspsvx.f 2011/07/22 07:38:20 1.8 +++ rpl/lapack/lapack/zspsvx.f 2011/11/21 20:43:20 1.9 @@ -1,10 +1,286 @@ +*> \brief ZSPSVX computes the solution to system of linear equations A * X = B for OTHER matrices +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZSPSVX + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZSPSVX( FACT, UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, +* LDX, RCOND, FERR, BERR, WORK, RWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER FACT, UPLO +* INTEGER INFO, LDB, LDX, N, NRHS +* DOUBLE PRECISION RCOND +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ) +* COMPLEX*16 AFP( * ), AP( * ), B( LDB, * ), WORK( * ), +* $ X( LDX, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZSPSVX uses the diagonal pivoting factorization A = U*D*U**T or +*> A = L*D*L**T to compute the solution to a complex system of linear +*> equations A * X = B, where A is an N-by-N symmetric matrix stored +*> in packed format and X and B are N-by-NRHS matrices. +*> +*> Error bounds on the solution and a condition estimate are also +*> provided. +*> \endverbatim +* +*> \par Description: +* ================= +*> +*> \verbatim +*> +*> The following steps are performed: +*> +*> 1. If FACT = 'N', the diagonal pivoting method is used to factor A as +*> A = U * D * U**T, if UPLO = 'U', or +*> A = L * D * L**T, if UPLO = 'L', +*> 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. +*> +*> 2. If some D(i,i)=0, so that D is exactly singular, then the routine +*> returns with INFO = i. Otherwise, the factored form of A is used +*> to estimate the condition number of the matrix A. If the +*> reciprocal of the condition number is less than machine precision, +*> INFO = N+1 is returned as a warning, but the routine still goes on +*> to solve for X and compute error bounds as described below. +*> +*> 3. The system of equations is solved for X using the factored form +*> of A. +*> +*> 4. Iterative refinement is applied to improve the computed solution +*> matrix and calculate error bounds and backward error estimates +*> for it. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] FACT +*> \verbatim +*> FACT is CHARACTER*1 +*> Specifies whether or not the factored form of A has been +*> supplied on entry. +*> = 'F': On entry, AFP and IPIV contain the factored form +*> of A. AP, AFP and IPIV will not be modified. +*> = 'N': The matrix A will be copied to AFP and factored. +*> \endverbatim +*> +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> = 'U': Upper triangle of A is stored; +*> = 'L': Lower triangle of A is stored. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of linear equations, i.e., the order of the +*> matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] NRHS +*> \verbatim +*> NRHS is INTEGER +*> The number of right hand sides, i.e., the number of columns +*> of the matrices B and X. NRHS >= 0. +*> \endverbatim +*> +*> \param[in] AP +*> \verbatim +*> AP is COMPLEX*16 array, dimension (N*(N+1)/2) +*> The upper or lower triangle of the symmetric matrix A, packed +*> columnwise in a linear array. The j-th column of A is stored +*> in the array AP as follows: +*> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; +*> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. +*> See below for further details. +*> \endverbatim +*> +*> \param[in,out] AFP +*> \verbatim +*> AFP is or output) COMPLEX*16 array, dimension (N*(N+1)/2) +*> If FACT = 'F', then AFP is an input argument and on entry +*> 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 ZSPTRF, stored as +*> a packed triangular matrix in the same storage format as A. +*> +*> If FACT = 'N', then AFP is an output argument and on exit +*> 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 ZSPTRF, stored as +*> a packed triangular matrix in the same storage format as A. +*> \endverbatim +*> +*> \param[in,out] IPIV +*> \verbatim +*> IPIV is or output) INTEGER array, dimension (N) +*> If FACT = 'F', then IPIV is an input argument and on entry +*> contains details of the interchanges and the block structure +*> of D, as determined by ZSPTRF. +*> 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. +*> +*> If FACT = 'N', then IPIV is an output argument and on exit +*> contains details of the interchanges and the block structure +*> of D, as determined by ZSPTRF. +*> \endverbatim +*> +*> \param[in] B +*> \verbatim +*> B is COMPLEX*16 array, dimension (LDB,NRHS) +*> The N-by-NRHS right hand side matrix B. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,N). +*> \endverbatim +*> +*> \param[out] X +*> \verbatim +*> X is COMPLEX*16 array, dimension (LDX,NRHS) +*> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. +*> \endverbatim +*> +*> \param[in] LDX +*> \verbatim +*> LDX is INTEGER +*> The leading dimension of the array X. LDX >= max(1,N). +*> \endverbatim +*> +*> \param[out] RCOND +*> \verbatim +*> RCOND is DOUBLE PRECISION +*> The estimate of the reciprocal condition number of the matrix +*> A. If RCOND is less than the machine precision (in +*> particular, if RCOND = 0), the matrix is singular to working +*> precision. This condition is indicated by a return code of +*> INFO > 0. +*> \endverbatim +*> +*> \param[out] FERR +*> \verbatim +*> 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. +*> \endverbatim +*> +*> \param[out] BERR +*> \verbatim +*> 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). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX*16 array, dimension (2*N) +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is DOUBLE PRECISION array, dimension (N) +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> > 0: if INFO = i, and i is +*> <= N: D(i,i) is exactly zero. The factorization +*> has been completed but the factor D is exactly +*> singular, so the solution and error bounds could +*> not be computed. RCOND = 0 is returned. +*> = N+1: D is nonsingular, but RCOND is less than machine +*> precision, meaning that the matrix is singular +*> to working precision. Nevertheless, the +*> solution and error bounds are computed because +*> there are a number of situations where the +*> computed solution can be more accurate than the +*> value of RCOND would suggest. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complex16OTHERsolve +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The packed storage scheme is illustrated by the following example +*> when N = 4, UPLO = 'U': +*> +*> Two-dimensional storage of the symmetric matrix A: +*> +*> a11 a12 a13 a14 +*> a22 a23 a24 +*> a33 a34 (aij = aji) +*> a44 +*> +*> Packed storage of the upper triangle of A: +*> +*> AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ] +*> \endverbatim +*> +* ===================================================================== SUBROUTINE ZSPSVX( FACT, UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, $ LDX, RCOND, FERR, BERR, WORK, RWORK, INFO ) * -* -- LAPACK driver routine (version 3.3.1) -- +* -- LAPACK driver routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -- April 2011 -- +* November 2011 * * .. Scalar Arguments .. CHARACTER FACT, UPLO @@ -18,173 +294,6 @@ $ X( LDX, * ) * .. * -* Purpose -* ======= -* -* ZSPSVX uses the diagonal pivoting factorization A = U*D*U**T or -* A = L*D*L**T to compute the solution to a complex system of linear -* equations A * X = B, where A is an N-by-N symmetric matrix stored -* in packed format and X and B are N-by-NRHS matrices. -* -* Error bounds on the solution and a condition estimate are also -* provided. -* -* Description -* =========== -* -* The following steps are performed: -* -* 1. If FACT = 'N', the diagonal pivoting method is used to factor A as -* A = U * D * U**T, if UPLO = 'U', or -* A = L * D * L**T, if UPLO = 'L', -* 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. -* -* 2. If some D(i,i)=0, so that D is exactly singular, then the routine -* returns with INFO = i. Otherwise, the factored form of A is used -* to estimate the condition number of the matrix A. If the -* reciprocal of the condition number is less than machine precision, -* INFO = N+1 is returned as a warning, but the routine still goes on -* to solve for X and compute error bounds as described below. -* -* 3. The system of equations is solved for X using the factored form -* of A. -* -* 4. Iterative refinement is applied to improve the computed solution -* matrix and calculate error bounds and backward error estimates -* for it. -* -* Arguments -* ========= -* -* FACT (input) CHARACTER*1 -* Specifies whether or not the factored form of A has been -* supplied on entry. -* = 'F': On entry, AFP and IPIV contain the factored form -* of A. AP, AFP and IPIV will not be modified. -* = 'N': The matrix A will be copied to AFP and factored. -* -* UPLO (input) CHARACTER*1 -* = 'U': Upper triangle of A is stored; -* = 'L': Lower triangle of A is stored. -* -* N (input) INTEGER -* The number of linear equations, i.e., the order of the -* matrix A. N >= 0. -* -* NRHS (input) INTEGER -* The number of right hand sides, i.e., the number of columns -* of the matrices B and X. NRHS >= 0. -* -* AP (input) COMPLEX*16 array, dimension (N*(N+1)/2) -* The upper or lower triangle of the symmetric matrix A, packed -* columnwise in a linear array. The j-th column of A is stored -* in the array AP as follows: -* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; -* if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. -* See below for further details. -* -* AFP (input or output) COMPLEX*16 array, dimension (N*(N+1)/2) -* If FACT = 'F', then AFP is an input argument and on entry -* 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 ZSPTRF, stored as -* a packed triangular matrix in the same storage format as A. -* -* If FACT = 'N', then AFP is an output argument and on exit -* 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 ZSPTRF, stored as -* a packed triangular matrix in the same storage format as A. -* -* IPIV (input or output) INTEGER array, dimension (N) -* If FACT = 'F', then IPIV is an input argument and on entry -* contains details of the interchanges and the block structure -* of D, as determined by ZSPTRF. -* 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. -* -* If FACT = 'N', then IPIV is an output argument and on exit -* contains details of the interchanges and the block structure -* of D, as determined by ZSPTRF. -* -* B (input) COMPLEX*16 array, dimension (LDB,NRHS) -* The N-by-NRHS right hand side matrix B. -* -* LDB (input) INTEGER -* The leading dimension of the array B. LDB >= max(1,N). -* -* X (output) COMPLEX*16 array, dimension (LDX,NRHS) -* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. -* -* LDX (input) INTEGER -* The leading dimension of the array X. LDX >= max(1,N). -* -* RCOND (output) DOUBLE PRECISION -* The estimate of the reciprocal condition number of the matrix -* A. If RCOND is less than the machine precision (in -* particular, if RCOND = 0), the matrix is singular to working -* precision. This condition is indicated by a return code of -* INFO > 0. -* -* FERR (output) 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. -* -* BERR (output) 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). -* -* WORK (workspace) COMPLEX*16 array, dimension (2*N) -* -* RWORK (workspace) DOUBLE PRECISION array, dimension (N) -* -* INFO (output) INTEGER -* = 0: successful exit -* < 0: if INFO = -i, the i-th argument had an illegal value -* > 0: if INFO = i, and i is -* <= N: D(i,i) is exactly zero. The factorization -* has been completed but the factor D is exactly -* singular, so the solution and error bounds could -* not be computed. RCOND = 0 is returned. -* = N+1: D is nonsingular, but RCOND is less than machine -* precision, meaning that the matrix is singular -* to working precision. Nevertheless, the -* solution and error bounds are computed because -* there are a number of situations where the -* computed solution can be more accurate than the -* value of RCOND would suggest. -* -* Further Details -* =============== -* -* The packed storage scheme is illustrated by the following example -* when N = 4, UPLO = 'U': -* -* Two-dimensional storage of the symmetric matrix A: -* -* a11 a12 a13 a14 -* a22 a23 a24 -* a33 a34 (aij = aji) -* a44 -* -* Packed storage of the upper triangle of A: -* -* AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ] -* * ===================================================================== * * .. Parameters ..