00001 /* sspsvx.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 sspsvx_(char *fact, char *uplo, integer *n, integer * 00021 nrhs, real *ap, real *afp, integer *ipiv, real *b, integer *ldb, real 00022 *x, integer *ldx, real *rcond, real *ferr, real *berr, real *work, 00023 integer *iwork, 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 scopy_(integer *, real *, integer *, real *, 00032 integer *); 00033 extern doublereal slamch_(char *); 00034 logical nofact; 00035 extern /* Subroutine */ int xerbla_(char *, integer *), slacpy_( 00036 char *, integer *, integer *, real *, integer *, real *, integer * 00037 ); 00038 extern doublereal slansp_(char *, char *, integer *, real *, real *); 00039 extern /* Subroutine */ int sspcon_(char *, integer *, real *, integer *, 00040 real *, real *, real *, integer *, integer *), ssprfs_( 00041 char *, integer *, integer *, real *, real *, integer *, real *, 00042 integer *, real *, integer *, real *, real *, real *, integer *, 00043 integer *), ssptrf_(char *, integer *, real *, integer *, 00044 integer *), ssptrs_(char *, integer *, integer *, real *, 00045 integer *, real *, 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 /* SSPSVX uses the diagonal pivoting factorization A = U*D*U**T or */ 00061 /* A = L*D*L**T to compute the solution to a real 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 of */ 00101 /* 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) REAL 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) REAL array, dimension */ 00125 /* (N*(N+1)/2) */ 00126 /* If FACT = 'F', then AFP is an input argument and on entry */ 00127 /* contains the block diagonal matrix D and the multipliers used */ 00128 /* to obtain the factor U or L from the factorization */ 00129 /* A = U*D*U**T or A = L*D*L**T as computed by SSPTRF, stored as */ 00130 /* a packed triangular matrix in the same storage format as A. */ 00131 00132 /* If FACT = 'N', then AFP is an output argument and on exit */ 00133 /* contains the block diagonal matrix D and the multipliers used */ 00134 /* to obtain the factor U or L from the factorization */ 00135 /* A = U*D*U**T or A = L*D*L**T as computed by SSPTRF, stored as */ 00136 /* a packed triangular matrix in the same storage format as A. */ 00137 00138 /* IPIV (input or output) INTEGER array, dimension (N) */ 00139 /* If FACT = 'F', then IPIV is an input argument and on entry */ 00140 /* contains details of the interchanges and the block structure */ 00141 /* of D, as determined by SSPTRF. */ 00142 /* If IPIV(k) > 0, then rows and columns k and IPIV(k) were */ 00143 /* interchanged and D(k,k) is a 1-by-1 diagonal block. */ 00144 /* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and */ 00145 /* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) */ 00146 /* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = */ 00147 /* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were */ 00148 /* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. */ 00149 00150 /* If FACT = 'N', then IPIV is an output argument and on exit */ 00151 /* contains details of the interchanges and the block structure */ 00152 /* of D, as determined by SSPTRF. */ 00153 00154 /* B (input) REAL array, dimension (LDB,NRHS) */ 00155 /* The N-by-NRHS right hand side matrix B. */ 00156 00157 /* LDB (input) INTEGER */ 00158 /* The leading dimension of the array B. LDB >= max(1,N). */ 00159 00160 /* X (output) REAL array, dimension (LDX,NRHS) */ 00161 /* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. */ 00162 00163 /* LDX (input) INTEGER */ 00164 /* The leading dimension of the array X. LDX >= max(1,N). */ 00165 00166 /* RCOND (output) REAL */ 00167 /* The estimate of the reciprocal condition number of the matrix */ 00168 /* A. If RCOND is less than the machine precision (in */ 00169 /* particular, if RCOND = 0), the matrix is singular to working */ 00170 /* precision. This condition is indicated by a return code of */ 00171 /* INFO > 0. */ 00172 00173 /* FERR (output) REAL array, dimension (NRHS) */ 00174 /* The estimated forward error bound for each solution vector */ 00175 /* X(j) (the j-th column of the solution matrix X). */ 00176 /* If XTRUE is the true solution corresponding to X(j), FERR(j) */ 00177 /* is an estimated upper bound for the magnitude of the largest */ 00178 /* element in (X(j) - XTRUE) divided by the magnitude of the */ 00179 /* largest element in X(j). The estimate is as reliable as */ 00180 /* the estimate for RCOND, and is almost always a slight */ 00181 /* overestimate of the true error. */ 00182 00183 /* BERR (output) REAL array, dimension (NRHS) */ 00184 /* The componentwise relative backward error of each solution */ 00185 /* vector X(j) (i.e., the smallest relative change in */ 00186 /* any element of A or B that makes X(j) an exact solution). */ 00187 00188 /* WORK (workspace) REAL array, dimension (3*N) */ 00189 00190 /* IWORK (workspace) INTEGER array, dimension (N) */ 00191 00192 /* INFO (output) INTEGER */ 00193 /* = 0: successful exit */ 00194 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00195 /* > 0: if INFO = i, and i is */ 00196 /* <= N: D(i,i) is exactly zero. The factorization */ 00197 /* has been completed but the factor D is exactly */ 00198 /* singular, so the solution and error bounds could */ 00199 /* not be computed. RCOND = 0 is returned. */ 00200 /* = N+1: D is nonsingular, but RCOND is less than machine */ 00201 /* precision, meaning that the matrix is singular */ 00202 /* to working precision. Nevertheless, the */ 00203 /* solution and error bounds are computed because */ 00204 /* there are a number of situations where the */ 00205 /* computed solution can be more accurate than the */ 00206 /* value of RCOND would suggest. */ 00207 00208 /* Further Details */ 00209 /* =============== */ 00210 00211 /* The packed storage scheme is illustrated by the following example */ 00212 /* when N = 4, UPLO = 'U': */ 00213 00214 /* Two-dimensional storage of the symmetric matrix A: */ 00215 00216 /* a11 a12 a13 a14 */ 00217 /* a22 a23 a24 */ 00218 /* a33 a34 (aij = aji) */ 00219 /* a44 */ 00220 00221 /* Packed storage of the upper triangle of A: */ 00222 00223 /* AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ] */ 00224 00225 /* ===================================================================== */ 00226 00227 /* .. Parameters .. */ 00228 /* .. */ 00229 /* .. Local Scalars .. */ 00230 /* .. */ 00231 /* .. External Functions .. */ 00232 /* .. */ 00233 /* .. External Subroutines .. */ 00234 /* .. */ 00235 /* .. Intrinsic Functions .. */ 00236 /* .. */ 00237 /* .. Executable Statements .. */ 00238 00239 /* Test the input parameters. */ 00240 00241 /* Parameter adjustments */ 00242 --ap; 00243 --afp; 00244 --ipiv; 00245 b_dim1 = *ldb; 00246 b_offset = 1 + b_dim1; 00247 b -= b_offset; 00248 x_dim1 = *ldx; 00249 x_offset = 1 + x_dim1; 00250 x -= x_offset; 00251 --ferr; 00252 --berr; 00253 --work; 00254 --iwork; 00255 00256 /* Function Body */ 00257 *info = 0; 00258 nofact = lsame_(fact, "N"); 00259 if (! nofact && ! lsame_(fact, "F")) { 00260 *info = -1; 00261 } else if (! lsame_(uplo, "U") && ! lsame_(uplo, 00262 "L")) { 00263 *info = -2; 00264 } else if (*n < 0) { 00265 *info = -3; 00266 } else if (*nrhs < 0) { 00267 *info = -4; 00268 } else if (*ldb < max(1,*n)) { 00269 *info = -9; 00270 } else if (*ldx < max(1,*n)) { 00271 *info = -11; 00272 } 00273 if (*info != 0) { 00274 i__1 = -(*info); 00275 xerbla_("SSPSVX", &i__1); 00276 return 0; 00277 } 00278 00279 if (nofact) { 00280 00281 /* Compute the factorization A = U*D*U' or A = L*D*L'. */ 00282 00283 i__1 = *n * (*n + 1) / 2; 00284 scopy_(&i__1, &ap[1], &c__1, &afp[1], &c__1); 00285 ssptrf_(uplo, n, &afp[1], &ipiv[1], info); 00286 00287 /* Return if INFO is non-zero. */ 00288 00289 if (*info > 0) { 00290 *rcond = 0.f; 00291 return 0; 00292 } 00293 } 00294 00295 /* Compute the norm of the matrix A. */ 00296 00297 anorm = slansp_("I", uplo, n, &ap[1], &work[1]); 00298 00299 /* Compute the reciprocal of the condition number of A. */ 00300 00301 sspcon_(uplo, n, &afp[1], &ipiv[1], &anorm, rcond, &work[1], &iwork[1], 00302 info); 00303 00304 /* Compute the solution vectors X. */ 00305 00306 slacpy_("Full", n, nrhs, &b[b_offset], ldb, &x[x_offset], ldx); 00307 ssptrs_(uplo, n, nrhs, &afp[1], &ipiv[1], &x[x_offset], ldx, info); 00308 00309 /* Use iterative refinement to improve the computed solutions and */ 00310 /* compute error bounds and backward error estimates for them. */ 00311 00312 ssprfs_(uplo, n, nrhs, &ap[1], &afp[1], &ipiv[1], &b[b_offset], ldb, &x[ 00313 x_offset], ldx, &ferr[1], &berr[1], &work[1], &iwork[1], info); 00314 00315 /* Set INFO = N+1 if the matrix is singular to working precision. */ 00316 00317 if (*rcond < slamch_("Epsilon")) { 00318 *info = *n + 1; 00319 } 00320 00321 return 0; 00322 00323 /* End of SSPSVX */ 00324 00325 } /* sspsvx_ */