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