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