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