00001 /* zppsvx.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 zppsvx_(char *fact, char *uplo, integer *n, integer * 00021 nrhs, doublecomplex *ap, doublecomplex *afp, char *equed, doublereal * 00022 s, 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, i__2, i__3, i__4, i__5; 00028 doublereal d__1, d__2; 00029 doublecomplex z__1; 00030 00031 /* Local variables */ 00032 integer i__, j; 00033 doublereal amax, smin, smax; 00034 extern logical lsame_(char *, char *); 00035 doublereal scond, anorm; 00036 logical equil, rcequ; 00037 extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 00038 doublecomplex *, integer *); 00039 extern doublereal dlamch_(char *); 00040 logical nofact; 00041 extern /* Subroutine */ int xerbla_(char *, integer *); 00042 doublereal bignum; 00043 integer infequ; 00044 extern doublereal zlanhp_(char *, char *, integer *, doublecomplex *, 00045 doublereal *); 00046 extern /* Subroutine */ int zlaqhp_(char *, integer *, doublecomplex *, 00047 doublereal *, doublereal *, doublereal *, char *), 00048 zlacpy_(char *, integer *, integer *, doublecomplex *, integer *, 00049 doublecomplex *, integer *), zppcon_(char *, integer *, 00050 doublecomplex *, doublereal *, doublereal *, doublecomplex *, 00051 doublereal *, integer *); 00052 doublereal smlnum; 00053 extern /* Subroutine */ int zppequ_(char *, integer *, doublecomplex *, 00054 doublereal *, doublereal *, doublereal *, integer *), 00055 zpprfs_(char *, integer *, integer *, doublecomplex *, 00056 doublecomplex *, doublecomplex *, integer *, doublecomplex *, 00057 integer *, doublereal *, doublereal *, doublecomplex *, 00058 doublereal *, integer *), zpptrf_(char *, integer *, 00059 doublecomplex *, integer *), zpptrs_(char *, integer *, 00060 integer *, doublecomplex *, doublecomplex *, integer *, integer *); 00061 00062 00063 /* -- LAPACK driver routine (version 3.2) -- */ 00064 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00065 /* November 2006 */ 00066 00067 /* .. Scalar Arguments .. */ 00068 /* .. */ 00069 /* .. Array Arguments .. */ 00070 /* .. */ 00071 00072 /* Purpose */ 00073 /* ======= */ 00074 00075 /* ZPPSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to */ 00076 /* compute the solution to a complex system of linear equations */ 00077 /* A * X = B, */ 00078 /* where A is an N-by-N Hermitian positive definite matrix stored in */ 00079 /* packed format and X and B are N-by-NRHS matrices. */ 00080 00081 /* Error bounds on the solution and a condition estimate are also */ 00082 /* provided. */ 00083 00084 /* Description */ 00085 /* =========== */ 00086 00087 /* The following steps are performed: */ 00088 00089 /* 1. If FACT = 'E', real scaling factors are computed to equilibrate */ 00090 /* the system: */ 00091 /* diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B */ 00092 /* Whether or not the system will be equilibrated depends on the */ 00093 /* scaling of the matrix A, but if equilibration is used, A is */ 00094 /* overwritten by diag(S)*A*diag(S) and B by diag(S)*B. */ 00095 00096 /* 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to */ 00097 /* factor the matrix A (after equilibration if FACT = 'E') as */ 00098 /* A = U'* U , if UPLO = 'U', or */ 00099 /* A = L * L', if UPLO = 'L', */ 00100 /* where U is an upper triangular matrix, L is a lower triangular */ 00101 /* matrix, and ' indicates conjugate transpose. */ 00102 00103 /* 3. If the leading i-by-i principal minor is not positive definite, */ 00104 /* then the routine returns with INFO = i. Otherwise, the factored */ 00105 /* form of A is used to estimate the condition number of the matrix */ 00106 /* A. If the reciprocal of the condition number is less than machine */ 00107 /* precision, INFO = N+1 is returned as a warning, but the routine */ 00108 /* still goes on to solve for X and compute error bounds as */ 00109 /* described below. */ 00110 00111 /* 4. The system of equations is solved for X using the factored form */ 00112 /* of A. */ 00113 00114 /* 5. Iterative refinement is applied to improve the computed solution */ 00115 /* matrix and calculate error bounds and backward error estimates */ 00116 /* for it. */ 00117 00118 /* 6. If equilibration was used, the matrix X is premultiplied by */ 00119 /* diag(S) so that it solves the original system before */ 00120 /* equilibration. */ 00121 00122 /* Arguments */ 00123 /* ========= */ 00124 00125 /* FACT (input) CHARACTER*1 */ 00126 /* Specifies whether or not the factored form of the matrix A is */ 00127 /* supplied on entry, and if not, whether the matrix A should be */ 00128 /* equilibrated before it is factored. */ 00129 /* = 'F': On entry, AFP contains the factored form of A. */ 00130 /* If EQUED = 'Y', the matrix A has been equilibrated */ 00131 /* with scaling factors given by S. AP and AFP will not */ 00132 /* be modified. */ 00133 /* = 'N': The matrix A will be copied to AFP and factored. */ 00134 /* = 'E': The matrix A will be equilibrated if necessary, then */ 00135 /* copied to AFP and factored. */ 00136 00137 /* UPLO (input) CHARACTER*1 */ 00138 /* = 'U': Upper triangle of A is stored; */ 00139 /* = 'L': Lower triangle of A is stored. */ 00140 00141 /* N (input) INTEGER */ 00142 /* The number of linear equations, i.e., the order of the */ 00143 /* matrix A. N >= 0. */ 00144 00145 /* NRHS (input) INTEGER */ 00146 /* The number of right hand sides, i.e., the number of columns */ 00147 /* of the matrices B and X. NRHS >= 0. */ 00148 00149 /* AP (input/output) COMPLEX*16 array, dimension (N*(N+1)/2) */ 00150 /* On entry, the upper or lower triangle of the Hermitian matrix */ 00151 /* A, packed columnwise in a linear array, except if FACT = 'F' */ 00152 /* and EQUED = 'Y', then A must contain the equilibrated matrix */ 00153 /* diag(S)*A*diag(S). The j-th column of A is stored in the */ 00154 /* array AP as follows: */ 00155 /* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */ 00156 /* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. */ 00157 /* See below for further details. A is not modified if */ 00158 /* FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit. */ 00159 00160 /* On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by */ 00161 /* diag(S)*A*diag(S). */ 00162 00163 /* AFP (input or output) COMPLEX*16 array, dimension (N*(N+1)/2) */ 00164 /* If FACT = 'F', then AFP is an input argument and on entry */ 00165 /* contains the triangular factor U or L from the Cholesky */ 00166 /* factorization A = U**H*U or A = L*L**H, in the same storage */ 00167 /* format as A. If EQUED .ne. 'N', then AFP is the factored */ 00168 /* form of the equilibrated matrix A. */ 00169 00170 /* If FACT = 'N', then AFP is an output argument and on exit */ 00171 /* returns the triangular factor U or L from the Cholesky */ 00172 /* factorization A = U**H*U or A = L*L**H of the original */ 00173 /* matrix A. */ 00174 00175 /* If FACT = 'E', then AFP is an output argument and on exit */ 00176 /* returns the triangular factor U or L from the Cholesky */ 00177 /* factorization A = U**H*U or A = L*L**H of the equilibrated */ 00178 /* matrix A (see the description of AP for the form of the */ 00179 /* equilibrated matrix). */ 00180 00181 /* EQUED (input or output) CHARACTER*1 */ 00182 /* Specifies the form of equilibration that was done. */ 00183 /* = 'N': No equilibration (always true if FACT = 'N'). */ 00184 /* = 'Y': Equilibration was done, i.e., A has been replaced by */ 00185 /* diag(S) * A * diag(S). */ 00186 /* EQUED is an input argument if FACT = 'F'; otherwise, it is an */ 00187 /* output argument. */ 00188 00189 /* S (input or output) DOUBLE PRECISION array, dimension (N) */ 00190 /* The scale factors for A; not accessed if EQUED = 'N'. S is */ 00191 /* an input argument if FACT = 'F'; otherwise, S is an output */ 00192 /* argument. If FACT = 'F' and EQUED = 'Y', each element of S */ 00193 /* must be positive. */ 00194 00195 /* B (input/output) COMPLEX*16 array, dimension (LDB,NRHS) */ 00196 /* On entry, the N-by-NRHS right hand side matrix B. */ 00197 /* On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y', */ 00198 /* B is overwritten by diag(S) * B. */ 00199 00200 /* LDB (input) INTEGER */ 00201 /* The leading dimension of the array B. LDB >= max(1,N). */ 00202 00203 /* X (output) COMPLEX*16 array, dimension (LDX,NRHS) */ 00204 /* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to */ 00205 /* the original system of equations. Note that if EQUED = 'Y', */ 00206 /* A and B are modified on exit, and the solution to the */ 00207 /* equilibrated system is inv(diag(S))*X. */ 00208 00209 /* LDX (input) INTEGER */ 00210 /* The leading dimension of the array X. LDX >= max(1,N). */ 00211 00212 /* RCOND (output) DOUBLE PRECISION */ 00213 /* The estimate of the reciprocal condition number of the matrix */ 00214 /* A after equilibration (if done). If RCOND is less than the */ 00215 /* machine precision (in particular, if RCOND = 0), the matrix */ 00216 /* is singular to working precision. This condition is */ 00217 /* indicated by a return code of INFO > 0. */ 00218 00219 /* FERR (output) DOUBLE PRECISION array, dimension (NRHS) */ 00220 /* The estimated forward error bound for each solution vector */ 00221 /* X(j) (the j-th column of the solution matrix X). */ 00222 /* If XTRUE is the true solution corresponding to X(j), FERR(j) */ 00223 /* is an estimated upper bound for the magnitude of the largest */ 00224 /* element in (X(j) - XTRUE) divided by the magnitude of the */ 00225 /* largest element in X(j). The estimate is as reliable as */ 00226 /* the estimate for RCOND, and is almost always a slight */ 00227 /* overestimate of the true error. */ 00228 00229 /* BERR (output) DOUBLE PRECISION array, dimension (NRHS) */ 00230 /* The componentwise relative backward error of each solution */ 00231 /* vector X(j) (i.e., the smallest relative change in */ 00232 /* any element of A or B that makes X(j) an exact solution). */ 00233 00234 /* WORK (workspace) COMPLEX*16 array, dimension (2*N) */ 00235 00236 /* RWORK (workspace) DOUBLE PRECISION array, dimension (N) */ 00237 00238 /* INFO (output) INTEGER */ 00239 /* = 0: successful exit */ 00240 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00241 /* > 0: if INFO = i, and i is */ 00242 /* <= N: the leading minor of order i of A is */ 00243 /* not positive definite, so the factorization */ 00244 /* could not be completed, and the solution has not */ 00245 /* been computed. RCOND = 0 is returned. */ 00246 /* = N+1: U is nonsingular, but RCOND is less than machine */ 00247 /* precision, meaning that the matrix is singular */ 00248 /* to working precision. Nevertheless, the */ 00249 /* solution and error bounds are computed because */ 00250 /* there are a number of situations where the */ 00251 /* computed solution can be more accurate than the */ 00252 /* value of RCOND would suggest. */ 00253 00254 /* Further Details */ 00255 /* =============== */ 00256 00257 /* The packed storage scheme is illustrated by the following example */ 00258 /* when N = 4, UPLO = 'U': */ 00259 00260 /* Two-dimensional storage of the Hermitian matrix A: */ 00261 00262 /* a11 a12 a13 a14 */ 00263 /* a22 a23 a24 */ 00264 /* a33 a34 (aij = conjg(aji)) */ 00265 /* a44 */ 00266 00267 /* Packed storage of the upper triangle of A: */ 00268 00269 /* AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ] */ 00270 00271 /* ===================================================================== */ 00272 00273 /* .. Parameters .. */ 00274 /* .. */ 00275 /* .. Local Scalars .. */ 00276 /* .. */ 00277 /* .. External Functions .. */ 00278 /* .. */ 00279 /* .. External Subroutines .. */ 00280 /* .. */ 00281 /* .. Intrinsic Functions .. */ 00282 /* .. */ 00283 /* .. Executable Statements .. */ 00284 00285 /* Parameter adjustments */ 00286 --ap; 00287 --afp; 00288 --s; 00289 b_dim1 = *ldb; 00290 b_offset = 1 + b_dim1; 00291 b -= b_offset; 00292 x_dim1 = *ldx; 00293 x_offset = 1 + x_dim1; 00294 x -= x_offset; 00295 --ferr; 00296 --berr; 00297 --work; 00298 --rwork; 00299 00300 /* Function Body */ 00301 *info = 0; 00302 nofact = lsame_(fact, "N"); 00303 equil = lsame_(fact, "E"); 00304 if (nofact || equil) { 00305 *(unsigned char *)equed = 'N'; 00306 rcequ = FALSE_; 00307 } else { 00308 rcequ = lsame_(equed, "Y"); 00309 smlnum = dlamch_("Safe minimum"); 00310 bignum = 1. / smlnum; 00311 } 00312 00313 /* Test the input parameters. */ 00314 00315 if (! nofact && ! equil && ! lsame_(fact, "F")) { 00316 *info = -1; 00317 } else if (! lsame_(uplo, "U") && ! lsame_(uplo, 00318 "L")) { 00319 *info = -2; 00320 } else if (*n < 0) { 00321 *info = -3; 00322 } else if (*nrhs < 0) { 00323 *info = -4; 00324 } else if (lsame_(fact, "F") && ! (rcequ || lsame_( 00325 equed, "N"))) { 00326 *info = -7; 00327 } else { 00328 if (rcequ) { 00329 smin = bignum; 00330 smax = 0.; 00331 i__1 = *n; 00332 for (j = 1; j <= i__1; ++j) { 00333 /* Computing MIN */ 00334 d__1 = smin, d__2 = s[j]; 00335 smin = min(d__1,d__2); 00336 /* Computing MAX */ 00337 d__1 = smax, d__2 = s[j]; 00338 smax = max(d__1,d__2); 00339 /* L10: */ 00340 } 00341 if (smin <= 0.) { 00342 *info = -8; 00343 } else if (*n > 0) { 00344 scond = max(smin,smlnum) / min(smax,bignum); 00345 } else { 00346 scond = 1.; 00347 } 00348 } 00349 if (*info == 0) { 00350 if (*ldb < max(1,*n)) { 00351 *info = -10; 00352 } else if (*ldx < max(1,*n)) { 00353 *info = -12; 00354 } 00355 } 00356 } 00357 00358 if (*info != 0) { 00359 i__1 = -(*info); 00360 xerbla_("ZPPSVX", &i__1); 00361 return 0; 00362 } 00363 00364 if (equil) { 00365 00366 /* Compute row and column scalings to equilibrate the matrix A. */ 00367 00368 zppequ_(uplo, n, &ap[1], &s[1], &scond, &amax, &infequ); 00369 if (infequ == 0) { 00370 00371 /* Equilibrate the matrix. */ 00372 00373 zlaqhp_(uplo, n, &ap[1], &s[1], &scond, &amax, equed); 00374 rcequ = lsame_(equed, "Y"); 00375 } 00376 } 00377 00378 /* Scale the right-hand side. */ 00379 00380 if (rcequ) { 00381 i__1 = *nrhs; 00382 for (j = 1; j <= i__1; ++j) { 00383 i__2 = *n; 00384 for (i__ = 1; i__ <= i__2; ++i__) { 00385 i__3 = i__ + j * b_dim1; 00386 i__4 = i__; 00387 i__5 = i__ + j * b_dim1; 00388 z__1.r = s[i__4] * b[i__5].r, z__1.i = s[i__4] * b[i__5].i; 00389 b[i__3].r = z__1.r, b[i__3].i = z__1.i; 00390 /* L20: */ 00391 } 00392 /* L30: */ 00393 } 00394 } 00395 00396 if (nofact || equil) { 00397 00398 /* Compute the Cholesky factorization A = U'*U or A = L*L'. */ 00399 00400 i__1 = *n * (*n + 1) / 2; 00401 zcopy_(&i__1, &ap[1], &c__1, &afp[1], &c__1); 00402 zpptrf_(uplo, n, &afp[1], info); 00403 00404 /* Return if INFO is non-zero. */ 00405 00406 if (*info > 0) { 00407 *rcond = 0.; 00408 return 0; 00409 } 00410 } 00411 00412 /* Compute the norm of the matrix A. */ 00413 00414 anorm = zlanhp_("I", uplo, n, &ap[1], &rwork[1]); 00415 00416 /* Compute the reciprocal of the condition number of A. */ 00417 00418 zppcon_(uplo, n, &afp[1], &anorm, rcond, &work[1], &rwork[1], info); 00419 00420 /* Compute the solution matrix X. */ 00421 00422 zlacpy_("Full", n, nrhs, &b[b_offset], ldb, &x[x_offset], ldx); 00423 zpptrs_(uplo, n, nrhs, &afp[1], &x[x_offset], ldx, info); 00424 00425 /* Use iterative refinement to improve the computed solution and */ 00426 /* compute error bounds and backward error estimates for it. */ 00427 00428 zpprfs_(uplo, n, nrhs, &ap[1], &afp[1], &b[b_offset], ldb, &x[x_offset], 00429 ldx, &ferr[1], &berr[1], &work[1], &rwork[1], info); 00430 00431 /* Transform the solution matrix X to a solution of the original */ 00432 /* system. */ 00433 00434 if (rcequ) { 00435 i__1 = *nrhs; 00436 for (j = 1; j <= i__1; ++j) { 00437 i__2 = *n; 00438 for (i__ = 1; i__ <= i__2; ++i__) { 00439 i__3 = i__ + j * x_dim1; 00440 i__4 = i__; 00441 i__5 = i__ + j * x_dim1; 00442 z__1.r = s[i__4] * x[i__5].r, z__1.i = s[i__4] * x[i__5].i; 00443 x[i__3].r = z__1.r, x[i__3].i = z__1.i; 00444 /* L40: */ 00445 } 00446 /* L50: */ 00447 } 00448 i__1 = *nrhs; 00449 for (j = 1; j <= i__1; ++j) { 00450 ferr[j] /= scond; 00451 /* L60: */ 00452 } 00453 } 00454 00455 /* Set INFO = N+1 if the matrix is singular to working precision. */ 00456 00457 if (*rcond < dlamch_("Epsilon")) { 00458 *info = *n + 1; 00459 } 00460 00461 return 0; 00462 00463 /* End of ZPPSVX */ 00464 00465 } /* zppsvx_ */