00001 /* zsysvx.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 static integer c_n1 = -1; 00020 00021 /* Subroutine */ int zsysvx_(char *fact, char *uplo, integer *n, integer * 00022 nrhs, doublecomplex *a, integer *lda, doublecomplex *af, integer * 00023 ldaf, integer *ipiv, doublecomplex *b, integer *ldb, doublecomplex *x, 00024 integer *ldx, doublereal *rcond, doublereal *ferr, doublereal *berr, 00025 doublecomplex *work, integer *lwork, doublereal *rwork, integer *info) 00026 { 00027 /* System generated locals */ 00028 integer a_dim1, a_offset, af_dim1, af_offset, b_dim1, b_offset, x_dim1, 00029 x_offset, i__1, i__2; 00030 00031 /* Local variables */ 00032 integer nb; 00033 extern logical lsame_(char *, char *); 00034 doublereal anorm; 00035 extern doublereal dlamch_(char *); 00036 logical nofact; 00037 extern /* Subroutine */ int xerbla_(char *, integer *); 00038 extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 00039 integer *, integer *); 00040 extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 00041 doublecomplex *, integer *, doublecomplex *, integer *); 00042 integer lwkopt; 00043 logical lquery; 00044 extern doublereal zlansy_(char *, char *, integer *, doublecomplex *, 00045 integer *, doublereal *); 00046 extern /* Subroutine */ int zsycon_(char *, integer *, doublecomplex *, 00047 integer *, integer *, doublereal *, doublereal *, doublecomplex *, 00048 integer *), zsyrfs_(char *, integer *, integer *, 00049 doublecomplex *, integer *, doublecomplex *, integer *, integer *, 00050 doublecomplex *, integer *, doublecomplex *, integer *, 00051 doublereal *, doublereal *, doublecomplex *, doublereal *, 00052 integer *), zsytrf_(char *, integer *, doublecomplex *, 00053 integer *, integer *, doublecomplex *, integer *, integer *), zsytrs_(char *, integer *, integer *, doublecomplex *, 00054 integer *, integer *, doublecomplex *, integer *, integer *); 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 /* ZSYSVX uses the diagonal pivoting factorization to compute the */ 00070 /* solution to a complex system of linear equations A * X = B, */ 00071 /* where A is an N-by-N symmetric matrix and X and B are N-by-NRHS */ 00072 /* matrices. */ 00073 00074 /* Error bounds on the solution and a condition estimate are also */ 00075 /* provided. */ 00076 00077 /* Description */ 00078 /* =========== */ 00079 00080 /* The following steps are performed: */ 00081 00082 /* 1. If FACT = 'N', the diagonal pivoting method is used to factor A. */ 00083 /* The form of the factorization is */ 00084 /* A = U * D * U**T, if UPLO = 'U', or */ 00085 /* A = L * D * L**T, if UPLO = 'L', */ 00086 /* where U (or L) is a product of permutation and unit upper (lower) */ 00087 /* triangular matrices, and D is symmetric and block diagonal with */ 00088 /* 1-by-1 and 2-by-2 diagonal blocks. */ 00089 00090 /* 2. If some D(i,i)=0, so that D is exactly singular, then the routine */ 00091 /* returns with INFO = i. Otherwise, the factored form of A is used */ 00092 /* to estimate the condition number of the matrix A. If the */ 00093 /* reciprocal of the condition number is less than machine precision, */ 00094 /* INFO = N+1 is returned as a warning, but the routine still goes on */ 00095 /* to solve for X and compute error bounds as described below. */ 00096 00097 /* 3. The system of equations is solved for X using the factored form */ 00098 /* of A. */ 00099 00100 /* 4. Iterative refinement is applied to improve the computed solution */ 00101 /* matrix and calculate error bounds and backward error estimates */ 00102 /* for it. */ 00103 00104 /* Arguments */ 00105 /* ========= */ 00106 00107 /* FACT (input) CHARACTER*1 */ 00108 /* Specifies whether or not the factored form of A has been */ 00109 /* supplied on entry. */ 00110 /* = 'F': On entry, AF and IPIV contain the factored form */ 00111 /* of A. A, AF and IPIV will not be modified. */ 00112 /* = 'N': The matrix A will be copied to AF and factored. */ 00113 00114 /* UPLO (input) CHARACTER*1 */ 00115 /* = 'U': Upper triangle of A is stored; */ 00116 /* = 'L': Lower triangle of A is stored. */ 00117 00118 /* N (input) INTEGER */ 00119 /* The number of linear equations, i.e., the order of the */ 00120 /* matrix A. N >= 0. */ 00121 00122 /* NRHS (input) INTEGER */ 00123 /* The number of right hand sides, i.e., the number of columns */ 00124 /* of the matrices B and X. NRHS >= 0. */ 00125 00126 /* A (input) COMPLEX*16 array, dimension (LDA,N) */ 00127 /* The symmetric matrix A. If UPLO = 'U', the leading N-by-N */ 00128 /* upper triangular part of A contains the upper triangular part */ 00129 /* of the matrix A, and the strictly lower triangular part of A */ 00130 /* is not referenced. If UPLO = 'L', the leading N-by-N lower */ 00131 /* triangular part of A contains the lower triangular part of */ 00132 /* the matrix A, and the strictly upper triangular part of A is */ 00133 /* not referenced. */ 00134 00135 /* LDA (input) INTEGER */ 00136 /* The leading dimension of the array A. LDA >= max(1,N). */ 00137 00138 /* AF (input or output) COMPLEX*16 array, dimension (LDAF,N) */ 00139 /* If FACT = 'F', then AF is an input argument and on entry */ 00140 /* contains the block diagonal matrix D and the multipliers used */ 00141 /* to obtain the factor U or L from the factorization */ 00142 /* A = U*D*U**T or A = L*D*L**T as computed by ZSYTRF. */ 00143 00144 /* If FACT = 'N', then AF is an output argument and on exit */ 00145 /* returns the block diagonal matrix D and the multipliers used */ 00146 /* to obtain the factor U or L from the factorization */ 00147 /* A = U*D*U**T or A = L*D*L**T. */ 00148 00149 /* LDAF (input) INTEGER */ 00150 /* The leading dimension of the array AF. LDAF >= max(1,N). */ 00151 00152 /* IPIV (input or output) INTEGER array, dimension (N) */ 00153 /* If FACT = 'F', then IPIV is an input argument and on entry */ 00154 /* contains details of the interchanges and the block structure */ 00155 /* of D, as determined by ZSYTRF. */ 00156 /* If IPIV(k) > 0, then rows and columns k and IPIV(k) were */ 00157 /* interchanged and D(k,k) is a 1-by-1 diagonal block. */ 00158 /* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and */ 00159 /* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) */ 00160 /* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = */ 00161 /* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were */ 00162 /* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. */ 00163 00164 /* If FACT = 'N', then IPIV is an output argument and on exit */ 00165 /* contains details of the interchanges and the block structure */ 00166 /* of D, as determined by ZSYTRF. */ 00167 00168 /* B (input) COMPLEX*16 array, dimension (LDB,NRHS) */ 00169 /* The N-by-NRHS right hand side matrix B. */ 00170 00171 /* LDB (input) INTEGER */ 00172 /* The leading dimension of the array B. LDB >= max(1,N). */ 00173 00174 /* X (output) COMPLEX*16 array, dimension (LDX,NRHS) */ 00175 /* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. */ 00176 00177 /* LDX (input) INTEGER */ 00178 /* The leading dimension of the array X. LDX >= max(1,N). */ 00179 00180 /* RCOND (output) DOUBLE PRECISION */ 00181 /* The estimate of the reciprocal condition number of the matrix */ 00182 /* A. If RCOND is less than the machine precision (in */ 00183 /* particular, if RCOND = 0), the matrix is singular to working */ 00184 /* precision. This condition is indicated by a return code of */ 00185 /* INFO > 0. */ 00186 00187 /* FERR (output) DOUBLE PRECISION array, dimension (NRHS) */ 00188 /* The estimated forward error bound for each solution vector */ 00189 /* X(j) (the j-th column of the solution matrix X). */ 00190 /* If XTRUE is the true solution corresponding to X(j), FERR(j) */ 00191 /* is an estimated upper bound for the magnitude of the largest */ 00192 /* element in (X(j) - XTRUE) divided by the magnitude of the */ 00193 /* largest element in X(j). The estimate is as reliable as */ 00194 /* the estimate for RCOND, and is almost always a slight */ 00195 /* overestimate of the true error. */ 00196 00197 /* BERR (output) DOUBLE PRECISION array, dimension (NRHS) */ 00198 /* The componentwise relative backward error of each solution */ 00199 /* vector X(j) (i.e., the smallest relative change in */ 00200 /* any element of A or B that makes X(j) an exact solution). */ 00201 00202 /* WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) */ 00203 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */ 00204 00205 /* LWORK (input) INTEGER */ 00206 /* The length of WORK. LWORK >= max(1,2*N), and for best */ 00207 /* performance, when FACT = 'N', LWORK >= max(1,2*N,N*NB), where */ 00208 /* NB is the optimal blocksize for ZSYTRF. */ 00209 00210 /* If LWORK = -1, then a workspace query is assumed; the routine */ 00211 /* only calculates the optimal size of the WORK array, returns */ 00212 /* this value as the first entry of the WORK array, and no error */ 00213 /* message related to LWORK is issued by XERBLA. */ 00214 00215 /* RWORK (workspace) DOUBLE PRECISION array, dimension (N) */ 00216 00217 /* INFO (output) INTEGER */ 00218 /* = 0: successful exit */ 00219 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00220 /* > 0: if INFO = i, and i is */ 00221 /* <= N: D(i,i) is exactly zero. The factorization */ 00222 /* has been completed but the factor D is exactly */ 00223 /* singular, so the solution and error bounds could */ 00224 /* not be computed. RCOND = 0 is returned. */ 00225 /* = N+1: D is nonsingular, but RCOND is less than machine */ 00226 /* precision, meaning that the matrix is singular */ 00227 /* to working precision. Nevertheless, the */ 00228 /* solution and error bounds are computed because */ 00229 /* there are a number of situations where the */ 00230 /* computed solution can be more accurate than the */ 00231 /* value of RCOND would suggest. */ 00232 00233 /* ===================================================================== */ 00234 00235 /* .. Parameters .. */ 00236 /* .. */ 00237 /* .. Local Scalars .. */ 00238 /* .. */ 00239 /* .. External Functions .. */ 00240 /* .. */ 00241 /* .. External Subroutines .. */ 00242 /* .. */ 00243 /* .. Intrinsic Functions .. */ 00244 /* .. */ 00245 /* .. Executable Statements .. */ 00246 00247 /* Test the input parameters. */ 00248 00249 /* Parameter adjustments */ 00250 a_dim1 = *lda; 00251 a_offset = 1 + a_dim1; 00252 a -= a_offset; 00253 af_dim1 = *ldaf; 00254 af_offset = 1 + af_dim1; 00255 af -= af_offset; 00256 --ipiv; 00257 b_dim1 = *ldb; 00258 b_offset = 1 + b_dim1; 00259 b -= b_offset; 00260 x_dim1 = *ldx; 00261 x_offset = 1 + x_dim1; 00262 x -= x_offset; 00263 --ferr; 00264 --berr; 00265 --work; 00266 --rwork; 00267 00268 /* Function Body */ 00269 *info = 0; 00270 nofact = lsame_(fact, "N"); 00271 lquery = *lwork == -1; 00272 if (! nofact && ! lsame_(fact, "F")) { 00273 *info = -1; 00274 } else if (! lsame_(uplo, "U") && ! lsame_(uplo, 00275 "L")) { 00276 *info = -2; 00277 } else if (*n < 0) { 00278 *info = -3; 00279 } else if (*nrhs < 0) { 00280 *info = -4; 00281 } else if (*lda < max(1,*n)) { 00282 *info = -6; 00283 } else if (*ldaf < max(1,*n)) { 00284 *info = -8; 00285 } else if (*ldb < max(1,*n)) { 00286 *info = -11; 00287 } else if (*ldx < max(1,*n)) { 00288 *info = -13; 00289 } else /* if(complicated condition) */ { 00290 /* Computing MAX */ 00291 i__1 = 1, i__2 = *n << 1; 00292 if (*lwork < max(i__1,i__2) && ! lquery) { 00293 *info = -18; 00294 } 00295 } 00296 00297 if (*info == 0) { 00298 /* Computing MAX */ 00299 i__1 = 1, i__2 = *n << 1; 00300 lwkopt = max(i__1,i__2); 00301 if (nofact) { 00302 nb = ilaenv_(&c__1, "ZSYTRF", uplo, n, &c_n1, &c_n1, &c_n1); 00303 /* Computing MAX */ 00304 i__1 = lwkopt, i__2 = *n * nb; 00305 lwkopt = max(i__1,i__2); 00306 } 00307 work[1].r = (doublereal) lwkopt, work[1].i = 0.; 00308 } 00309 00310 if (*info != 0) { 00311 i__1 = -(*info); 00312 xerbla_("ZSYSVX", &i__1); 00313 return 0; 00314 } else if (lquery) { 00315 return 0; 00316 } 00317 00318 if (nofact) { 00319 00320 /* Compute the factorization A = U*D*U' or A = L*D*L'. */ 00321 00322 zlacpy_(uplo, n, n, &a[a_offset], lda, &af[af_offset], ldaf); 00323 zsytrf_(uplo, n, &af[af_offset], ldaf, &ipiv[1], &work[1], lwork, 00324 info); 00325 00326 /* Return if INFO is non-zero. */ 00327 00328 if (*info > 0) { 00329 *rcond = 0.; 00330 return 0; 00331 } 00332 } 00333 00334 /* Compute the norm of the matrix A. */ 00335 00336 anorm = zlansy_("I", uplo, n, &a[a_offset], lda, &rwork[1]); 00337 00338 /* Compute the reciprocal of the condition number of A. */ 00339 00340 zsycon_(uplo, n, &af[af_offset], ldaf, &ipiv[1], &anorm, rcond, &work[1], 00341 info); 00342 00343 /* Compute the solution vectors X. */ 00344 00345 zlacpy_("Full", n, nrhs, &b[b_offset], ldb, &x[x_offset], ldx); 00346 zsytrs_(uplo, n, nrhs, &af[af_offset], ldaf, &ipiv[1], &x[x_offset], ldx, 00347 info); 00348 00349 /* Use iterative refinement to improve the computed solutions and */ 00350 /* compute error bounds and backward error estimates for them. */ 00351 00352 zsyrfs_(uplo, n, nrhs, &a[a_offset], lda, &af[af_offset], ldaf, &ipiv[1], 00353 &b[b_offset], ldb, &x[x_offset], ldx, &ferr[1], &berr[1], &work[1] 00354 , &rwork[1], info); 00355 00356 /* Set INFO = N+1 if the matrix is singular to working precision. */ 00357 00358 if (*rcond < dlamch_("Epsilon")) { 00359 *info = *n + 1; 00360 } 00361 00362 work[1].r = (doublereal) lwkopt, work[1].i = 0.; 00363 00364 return 0; 00365 00366 /* End of ZSYSVX */ 00367 00368 } /* zsysvx_ */