00001 /* zspsvx.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static integer c__1 = 1; 00019 00020 /* Subroutine */ int zspsvx_(char *fact, char *uplo, integer *n, integer * 00021 nrhs, doublecomplex *ap, doublecomplex *afp, integer *ipiv, 00022 doublecomplex *b, integer *ldb, doublecomplex *x, integer *ldx, 00023 doublereal *rcond, doublereal *ferr, doublereal *berr, doublecomplex * 00024 work, doublereal *rwork, integer *info) 00025 { 00026 /* System generated locals */ 00027 integer b_dim1, b_offset, x_dim1, x_offset, i__1; 00028 00029 /* Local variables */ 00030 extern logical lsame_(char *, char *); 00031 doublereal anorm; 00032 extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 00033 doublecomplex *, integer *); 00034 extern doublereal dlamch_(char *); 00035 logical nofact; 00036 extern /* Subroutine */ int xerbla_(char *, integer *), zlacpy_( 00037 char *, integer *, integer *, doublecomplex *, integer *, 00038 doublecomplex *, integer *); 00039 extern doublereal zlansp_(char *, char *, integer *, doublecomplex *, 00040 doublereal *); 00041 extern /* Subroutine */ int zspcon_(char *, integer *, doublecomplex *, 00042 integer *, doublereal *, doublereal *, doublecomplex *, integer *), zsprfs_(char *, integer *, integer *, doublecomplex *, 00043 doublecomplex *, integer *, doublecomplex *, integer *, 00044 doublecomplex *, integer *, doublereal *, doublereal *, 00045 doublecomplex *, doublereal *, integer *), zsptrf_(char *, 00046 integer *, doublecomplex *, integer *, integer *), 00047 zsptrs_(char *, integer *, integer *, doublecomplex *, integer *, 00048 doublecomplex *, integer *, integer *); 00049 00050 00051 /* -- LAPACK driver routine (version 3.2) -- */ 00052 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00053 /* November 2006 */ 00054 00055 /* .. Scalar Arguments .. */ 00056 /* .. */ 00057 /* .. Array Arguments .. */ 00058 /* .. */ 00059 00060 /* Purpose */ 00061 /* ======= */ 00062 00063 /* ZSPSVX uses the diagonal pivoting factorization A = U*D*U**T or */ 00064 /* A = L*D*L**T to compute the solution to a complex system of linear */ 00065 /* equations A * X = B, where A is an N-by-N symmetric matrix stored */ 00066 /* in packed format and X and B are N-by-NRHS matrices. */ 00067 00068 /* Error bounds on the solution and a condition estimate are also */ 00069 /* provided. */ 00070 00071 /* Description */ 00072 /* =========== */ 00073 00074 /* The following steps are performed: */ 00075 00076 /* 1. If FACT = 'N', the diagonal pivoting method is used to factor A as */ 00077 /* A = U * D * U**T, if UPLO = 'U', or */ 00078 /* A = L * D * L**T, if UPLO = 'L', */ 00079 /* where U (or L) is a product of permutation and unit upper (lower) */ 00080 /* triangular matrices and D is symmetric and block diagonal with */ 00081 /* 1-by-1 and 2-by-2 diagonal blocks. */ 00082 00083 /* 2. If some D(i,i)=0, so that D is exactly singular, then the routine */ 00084 /* returns with INFO = i. Otherwise, the factored form of A is used */ 00085 /* to estimate the condition number of the matrix A. If the */ 00086 /* reciprocal of the condition number is less than machine precision, */ 00087 /* INFO = N+1 is returned as a warning, but the routine still goes on */ 00088 /* to solve for X and compute error bounds as described below. */ 00089 00090 /* 3. The system of equations is solved for X using the factored form */ 00091 /* of A. */ 00092 00093 /* 4. Iterative refinement is applied to improve the computed solution */ 00094 /* matrix and calculate error bounds and backward error estimates */ 00095 /* for it. */ 00096 00097 /* Arguments */ 00098 /* ========= */ 00099 00100 /* FACT (input) CHARACTER*1 */ 00101 /* Specifies whether or not the factored form of A has been */ 00102 /* supplied on entry. */ 00103 /* = 'F': On entry, AFP and IPIV contain the factored form */ 00104 /* of A. AP, AFP and IPIV will not be modified. */ 00105 /* = 'N': The matrix A will be copied to AFP and factored. */ 00106 00107 /* UPLO (input) CHARACTER*1 */ 00108 /* = 'U': Upper triangle of A is stored; */ 00109 /* = 'L': Lower triangle of A is stored. */ 00110 00111 /* N (input) INTEGER */ 00112 /* The number of linear equations, i.e., the order of the */ 00113 /* matrix A. N >= 0. */ 00114 00115 /* NRHS (input) INTEGER */ 00116 /* The number of right hand sides, i.e., the number of columns */ 00117 /* of the matrices B and X. NRHS >= 0. */ 00118 00119 /* AP (input) COMPLEX*16 array, dimension (N*(N+1)/2) */ 00120 /* The upper or lower triangle of the symmetric matrix A, packed */ 00121 /* columnwise in a linear array. The j-th column of A is stored */ 00122 /* in the array AP as follows: */ 00123 /* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */ 00124 /* if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. */ 00125 /* See below for further details. */ 00126 00127 /* AFP (input or output) COMPLEX*16 array, dimension (N*(N+1)/2) */ 00128 /* If FACT = 'F', then AFP is an input argument and on entry */ 00129 /* contains the block diagonal matrix D and the multipliers used */ 00130 /* to obtain the factor U or L from the factorization */ 00131 /* A = U*D*U**T or A = L*D*L**T as computed by ZSPTRF, stored as */ 00132 /* a packed triangular matrix in the same storage format as A. */ 00133 00134 /* If FACT = 'N', then AFP is an output argument and on exit */ 00135 /* contains the block diagonal matrix D and the multipliers used */ 00136 /* to obtain the factor U or L from the factorization */ 00137 /* A = U*D*U**T or A = L*D*L**T as computed by ZSPTRF, stored as */ 00138 /* a packed triangular matrix in the same storage format as A. */ 00139 00140 /* IPIV (input or output) INTEGER array, dimension (N) */ 00141 /* If FACT = 'F', then IPIV is an input argument and on entry */ 00142 /* contains details of the interchanges and the block structure */ 00143 /* of D, as determined by ZSPTRF. */ 00144 /* If IPIV(k) > 0, then rows and columns k and IPIV(k) were */ 00145 /* interchanged and D(k,k) is a 1-by-1 diagonal block. */ 00146 /* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and */ 00147 /* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) */ 00148 /* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = */ 00149 /* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were */ 00150 /* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. */ 00151 00152 /* If FACT = 'N', then IPIV is an output argument and on exit */ 00153 /* contains details of the interchanges and the block structure */ 00154 /* of D, as determined by ZSPTRF. */ 00155 00156 /* B (input) COMPLEX*16 array, dimension (LDB,NRHS) */ 00157 /* The N-by-NRHS right hand side matrix B. */ 00158 00159 /* LDB (input) INTEGER */ 00160 /* The leading dimension of the array B. LDB >= max(1,N). */ 00161 00162 /* X (output) COMPLEX*16 array, dimension (LDX,NRHS) */ 00163 /* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. */ 00164 00165 /* LDX (input) INTEGER */ 00166 /* The leading dimension of the array X. LDX >= max(1,N). */ 00167 00168 /* RCOND (output) DOUBLE PRECISION */ 00169 /* The estimate of the reciprocal condition number of the matrix */ 00170 /* A. If RCOND is less than the machine precision (in */ 00171 /* particular, if RCOND = 0), the matrix is singular to working */ 00172 /* precision. This condition is indicated by a return code of */ 00173 /* INFO > 0. */ 00174 00175 /* FERR (output) DOUBLE PRECISION array, dimension (NRHS) */ 00176 /* The estimated forward error bound for each solution vector */ 00177 /* X(j) (the j-th column of the solution matrix X). */ 00178 /* If XTRUE is the true solution corresponding to X(j), FERR(j) */ 00179 /* is an estimated upper bound for the magnitude of the largest */ 00180 /* element in (X(j) - XTRUE) divided by the magnitude of the */ 00181 /* largest element in X(j). The estimate is as reliable as */ 00182 /* the estimate for RCOND, and is almost always a slight */ 00183 /* overestimate of the true error. */ 00184 00185 /* BERR (output) DOUBLE PRECISION array, dimension (NRHS) */ 00186 /* The componentwise relative backward error of each solution */ 00187 /* vector X(j) (i.e., the smallest relative change in */ 00188 /* any element of A or B that makes X(j) an exact solution). */ 00189 00190 /* WORK (workspace) COMPLEX*16 array, dimension (2*N) */ 00191 00192 /* RWORK (workspace) DOUBLE PRECISION array, dimension (N) */ 00193 00194 /* INFO (output) INTEGER */ 00195 /* = 0: successful exit */ 00196 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00197 /* > 0: if INFO = i, and i is */ 00198 /* <= N: D(i,i) is exactly zero. The factorization */ 00199 /* has been completed but the factor D is exactly */ 00200 /* singular, so the solution and error bounds could */ 00201 /* not be computed. RCOND = 0 is returned. */ 00202 /* = N+1: D is nonsingular, but RCOND is less than machine */ 00203 /* precision, meaning that the matrix is singular */ 00204 /* to working precision. Nevertheless, the */ 00205 /* solution and error bounds are computed because */ 00206 /* there are a number of situations where the */ 00207 /* computed solution can be more accurate than the */ 00208 /* value of RCOND would suggest. */ 00209 00210 /* Further Details */ 00211 /* =============== */ 00212 00213 /* The packed storage scheme is illustrated by the following example */ 00214 /* when N = 4, UPLO = 'U': */ 00215 00216 /* Two-dimensional storage of the symmetric matrix A: */ 00217 00218 /* a11 a12 a13 a14 */ 00219 /* a22 a23 a24 */ 00220 /* a33 a34 (aij = aji) */ 00221 /* a44 */ 00222 00223 /* Packed storage of the upper triangle of A: */ 00224 00225 /* AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ] */ 00226 00227 /* ===================================================================== */ 00228 00229 /* .. Parameters .. */ 00230 /* .. */ 00231 /* .. Local Scalars .. */ 00232 /* .. */ 00233 /* .. External Functions .. */ 00234 /* .. */ 00235 /* .. External Subroutines .. */ 00236 /* .. */ 00237 /* .. Intrinsic Functions .. */ 00238 /* .. */ 00239 /* .. Executable Statements .. */ 00240 00241 /* Test the input parameters. */ 00242 00243 /* Parameter adjustments */ 00244 --ap; 00245 --afp; 00246 --ipiv; 00247 b_dim1 = *ldb; 00248 b_offset = 1 + b_dim1; 00249 b -= b_offset; 00250 x_dim1 = *ldx; 00251 x_offset = 1 + x_dim1; 00252 x -= x_offset; 00253 --ferr; 00254 --berr; 00255 --work; 00256 --rwork; 00257 00258 /* Function Body */ 00259 *info = 0; 00260 nofact = lsame_(fact, "N"); 00261 if (! nofact && ! lsame_(fact, "F")) { 00262 *info = -1; 00263 } else if (! lsame_(uplo, "U") && ! lsame_(uplo, 00264 "L")) { 00265 *info = -2; 00266 } else if (*n < 0) { 00267 *info = -3; 00268 } else if (*nrhs < 0) { 00269 *info = -4; 00270 } else if (*ldb < max(1,*n)) { 00271 *info = -9; 00272 } else if (*ldx < max(1,*n)) { 00273 *info = -11; 00274 } 00275 if (*info != 0) { 00276 i__1 = -(*info); 00277 xerbla_("ZSPSVX", &i__1); 00278 return 0; 00279 } 00280 00281 if (nofact) { 00282 00283 /* Compute the factorization A = U*D*U' or A = L*D*L'. */ 00284 00285 i__1 = *n * (*n + 1) / 2; 00286 zcopy_(&i__1, &ap[1], &c__1, &afp[1], &c__1); 00287 zsptrf_(uplo, n, &afp[1], &ipiv[1], info); 00288 00289 /* Return if INFO is non-zero. */ 00290 00291 if (*info > 0) { 00292 *rcond = 0.; 00293 return 0; 00294 } 00295 } 00296 00297 /* Compute the norm of the matrix A. */ 00298 00299 anorm = zlansp_("I", uplo, n, &ap[1], &rwork[1]); 00300 00301 /* Compute the reciprocal of the condition number of A. */ 00302 00303 zspcon_(uplo, n, &afp[1], &ipiv[1], &anorm, rcond, &work[1], info); 00304 00305 /* Compute the solution vectors X. */ 00306 00307 zlacpy_("Full", n, nrhs, &b[b_offset], ldb, &x[x_offset], ldx); 00308 zsptrs_(uplo, n, nrhs, &afp[1], &ipiv[1], &x[x_offset], ldx, info); 00309 00310 /* Use iterative refinement to improve the computed solutions and */ 00311 /* compute error bounds and backward error estimates for them. */ 00312 00313 zsprfs_(uplo, n, nrhs, &ap[1], &afp[1], &ipiv[1], &b[b_offset], ldb, &x[ 00314 x_offset], ldx, &ferr[1], &berr[1], &work[1], &rwork[1], info); 00315 00316 /* Set INFO = N+1 if the matrix is singular to working precision. */ 00317 00318 if (*rcond < dlamch_("Epsilon")) { 00319 *info = *n + 1; 00320 } 00321 00322 return 0; 00323 00324 /* End of ZSPSVX */ 00325 00326 } /* zspsvx_ */