1 /* ./src_f77/zspsvx.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static integer c__1 = 1;
11
zspsvx_(char * fact,char * uplo,integer * n,integer * nrhs,doublecomplex * ap,doublecomplex * afp,integer * ipiv,doublecomplex * b,integer * ldb,doublecomplex * x,integer * ldx,doublereal * rcond,doublereal * ferr,doublereal * berr,doublecomplex * work,doublereal * rwork,integer * info,ftnlen fact_len,ftnlen uplo_len)12 /* Subroutine */ int zspsvx_(char *fact, char *uplo, integer *n, integer *
13 nrhs, doublecomplex *ap, doublecomplex *afp, integer *ipiv,
14 doublecomplex *b, integer *ldb, doublecomplex *x, integer *ldx,
15 doublereal *rcond, doublereal *ferr, doublereal *berr, doublecomplex *
16 work, doublereal *rwork, integer *info, ftnlen fact_len, ftnlen
17 uplo_len)
18 {
19 /* System generated locals */
20 integer b_dim1, b_offset, x_dim1, x_offset, i__1;
21
22 /* Local variables */
23 extern logical lsame_(char *, char *, ftnlen, ftnlen);
24 static doublereal anorm;
25 extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *,
26 doublecomplex *, integer *);
27 extern doublereal dlamch_(char *, ftnlen);
28 static logical nofact;
29 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen), zlacpy_(
30 char *, integer *, integer *, doublecomplex *, integer *,
31 doublecomplex *, integer *, ftnlen);
32 extern doublereal zlansp_(char *, char *, integer *, doublecomplex *,
33 doublereal *, ftnlen, ftnlen);
34 extern /* Subroutine */ int zspcon_(char *, integer *, doublecomplex *,
35 integer *, doublereal *, doublereal *, doublecomplex *, integer *,
36 ftnlen), zsprfs_(char *, integer *, integer *, doublecomplex *,
37 doublecomplex *, integer *, doublecomplex *, integer *,
38 doublecomplex *, integer *, doublereal *, doublereal *,
39 doublecomplex *, doublereal *, integer *, ftnlen), zsptrf_(char *,
40 integer *, doublecomplex *, integer *, integer *, ftnlen),
41 zsptrs_(char *, integer *, integer *, doublecomplex *, integer *,
42 doublecomplex *, integer *, integer *, ftnlen);
43
44
45 /* -- LAPACK driver routine (version 3.0) -- */
46 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
47 /* Courant Institute, Argonne National Lab, and Rice University */
48 /* June 30, 1999 */
49
50 /* .. Scalar Arguments .. */
51 /* .. */
52 /* .. Array Arguments .. */
53 /* .. */
54
55 /* Purpose */
56 /* ======= */
57
58 /* ZSPSVX uses the diagonal pivoting factorization A = U*D*U**T or */
59 /* A = L*D*L**T to compute the solution to a complex system of linear */
60 /* equations A * X = B, where A is an N-by-N symmetric matrix stored */
61 /* in packed format and X and B are N-by-NRHS matrices. */
62
63 /* Error bounds on the solution and a condition estimate are also */
64 /* provided. */
65
66 /* Description */
67 /* =========== */
68
69 /* The following steps are performed: */
70
71 /* 1. If FACT = 'N', the diagonal pivoting method is used to factor A as */
72 /* A = U * D * U**T, if UPLO = 'U', or */
73 /* A = L * D * L**T, if UPLO = 'L', */
74 /* where U (or L) is a product of permutation and unit upper (lower) */
75 /* triangular matrices and D is symmetric and block diagonal with */
76 /* 1-by-1 and 2-by-2 diagonal blocks. */
77
78 /* 2. If some D(i,i)=0, so that D is exactly singular, then the routine */
79 /* returns with INFO = i. Otherwise, the factored form of A is used */
80 /* to estimate the condition number of the matrix A. If the */
81 /* reciprocal of the condition number is less than machine precision, */
82 /* INFO = N+1 is returned as a warning, but the routine still goes on */
83 /* to solve for X and compute error bounds as described below. */
84
85 /* 3. The system of equations is solved for X using the factored form */
86 /* of A. */
87
88 /* 4. Iterative refinement is applied to improve the computed solution */
89 /* matrix and calculate error bounds and backward error estimates */
90 /* for it. */
91
92 /* Arguments */
93 /* ========= */
94
95 /* FACT (input) CHARACTER*1 */
96 /* Specifies whether or not the factored form of A has been */
97 /* supplied on entry. */
98 /* = 'F': On entry, AFP and IPIV contain the factored form */
99 /* of A. AP, AFP and IPIV will not be modified. */
100 /* = 'N': The matrix A will be copied to AFP and factored. */
101
102 /* UPLO (input) CHARACTER*1 */
103 /* = 'U': Upper triangle of A is stored; */
104 /* = 'L': Lower triangle of A is stored. */
105
106 /* N (input) INTEGER */
107 /* The number of linear equations, i.e., the order of the */
108 /* matrix A. N >= 0. */
109
110 /* NRHS (input) INTEGER */
111 /* The number of right hand sides, i.e., the number of columns */
112 /* of the matrices B and X. NRHS >= 0. */
113
114 /* AP (input) COMPLEX*16 array, dimension (N*(N+1)/2) */
115 /* The upper or lower triangle of the symmetric matrix A, packed */
116 /* columnwise in a linear array. The j-th column of A is stored */
117 /* in the array AP as follows: */
118 /* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */
119 /* if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. */
120 /* See below for further details. */
121
122 /* AFP (input or output) COMPLEX*16 array, dimension (N*(N+1)/2) */
123 /* If FACT = 'F', then AFP is an input argument and on entry */
124 /* contains the block diagonal matrix D and the multipliers used */
125 /* to obtain the factor U or L from the factorization */
126 /* A = U*D*U**T or A = L*D*L**T as computed by ZSPTRF, stored as */
127 /* a packed triangular matrix in the same storage format as A. */
128
129 /* If FACT = 'N', then AFP is an output argument and on exit */
130 /* contains the block diagonal matrix D and the multipliers used */
131 /* to obtain the factor U or L from the factorization */
132 /* A = U*D*U**T or A = L*D*L**T as computed by ZSPTRF, stored as */
133 /* a packed triangular matrix in the same storage format as A. */
134
135 /* IPIV (input or output) INTEGER array, dimension (N) */
136 /* If FACT = 'F', then IPIV is an input argument and on entry */
137 /* contains details of the interchanges and the block structure */
138 /* of D, as determined by ZSPTRF. */
139 /* If IPIV(k) > 0, then rows and columns k and IPIV(k) were */
140 /* interchanged and D(k,k) is a 1-by-1 diagonal block. */
141 /* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and */
142 /* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) */
143 /* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = */
144 /* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were */
145 /* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. */
146
147 /* If FACT = 'N', then IPIV is an output argument and on exit */
148 /* contains details of the interchanges and the block structure */
149 /* of D, as determined by ZSPTRF. */
150
151 /* B (input) COMPLEX*16 array, dimension (LDB,NRHS) */
152 /* The N-by-NRHS right hand side matrix B. */
153
154 /* LDB (input) INTEGER */
155 /* The leading dimension of the array B. LDB >= max(1,N). */
156
157 /* X (output) COMPLEX*16 array, dimension (LDX,NRHS) */
158 /* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. */
159
160 /* LDX (input) INTEGER */
161 /* The leading dimension of the array X. LDX >= max(1,N). */
162
163 /* RCOND (output) DOUBLE PRECISION */
164 /* The estimate of the reciprocal condition number of the matrix */
165 /* A. If RCOND is less than the machine precision (in */
166 /* particular, if RCOND = 0), the matrix is singular to working */
167 /* precision. This condition is indicated by a return code of */
168 /* INFO > 0. */
169
170 /* FERR (output) DOUBLE PRECISION array, dimension (NRHS) */
171 /* The estimated forward error bound for each solution vector */
172 /* X(j) (the j-th column of the solution matrix X). */
173 /* If XTRUE is the true solution corresponding to X(j), FERR(j) */
174 /* is an estimated upper bound for the magnitude of the largest */
175 /* element in (X(j) - XTRUE) divided by the magnitude of the */
176 /* largest element in X(j). The estimate is as reliable as */
177 /* the estimate for RCOND, and is almost always a slight */
178 /* overestimate of the true error. */
179
180 /* BERR (output) DOUBLE PRECISION array, dimension (NRHS) */
181 /* The componentwise relative backward error of each solution */
182 /* vector X(j) (i.e., the smallest relative change in */
183 /* any element of A or B that makes X(j) an exact solution). */
184
185 /* WORK (workspace) COMPLEX*16 array, dimension (2*N) */
186
187 /* RWORK (workspace) DOUBLE PRECISION array, dimension (N) */
188
189 /* INFO (output) INTEGER */
190 /* = 0: successful exit */
191 /* < 0: if INFO = -i, the i-th argument had an illegal value */
192 /* > 0: if INFO = i, and i is */
193 /* <= N: D(i,i) is exactly zero. The factorization */
194 /* has been completed but the factor D is exactly */
195 /* singular, so the solution and error bounds could */
196 /* not be computed. RCOND = 0 is returned. */
197 /* = N+1: D is nonsingular, but RCOND is less than machine */
198 /* precision, meaning that the matrix is singular */
199 /* to working precision. Nevertheless, the */
200 /* solution and error bounds are computed because */
201 /* there are a number of situations where the */
202 /* computed solution can be more accurate than the */
203 /* value of RCOND would suggest. */
204
205 /* Further Details */
206 /* =============== */
207
208 /* The packed storage scheme is illustrated by the following example */
209 /* when N = 4, UPLO = 'U': */
210
211 /* Two-dimensional storage of the symmetric matrix A: */
212
213 /* a11 a12 a13 a14 */
214 /* a22 a23 a24 */
215 /* a33 a34 (aij = aji) */
216 /* a44 */
217
218 /* Packed storage of the upper triangle of A: */
219
220 /* AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ] */
221
222 /* ===================================================================== */
223
224 /* .. Parameters .. */
225 /* .. */
226 /* .. Local Scalars .. */
227 /* .. */
228 /* .. External Functions .. */
229 /* .. */
230 /* .. External Subroutines .. */
231 /* .. */
232 /* .. Intrinsic Functions .. */
233 /* .. */
234 /* .. Executable Statements .. */
235
236 /* Test the input parameters. */
237
238 /* Parameter adjustments */
239 --ap;
240 --afp;
241 --ipiv;
242 b_dim1 = *ldb;
243 b_offset = 1 + b_dim1;
244 b -= b_offset;
245 x_dim1 = *ldx;
246 x_offset = 1 + x_dim1;
247 x -= x_offset;
248 --ferr;
249 --berr;
250 --work;
251 --rwork;
252
253 /* Function Body */
254 *info = 0;
255 nofact = lsame_(fact, "N", (ftnlen)1, (ftnlen)1);
256 if (! nofact && ! lsame_(fact, "F", (ftnlen)1, (ftnlen)1)) {
257 *info = -1;
258 } else if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo,
259 "L", (ftnlen)1, (ftnlen)1)) {
260 *info = -2;
261 } else if (*n < 0) {
262 *info = -3;
263 } else if (*nrhs < 0) {
264 *info = -4;
265 } else if (*ldb < max(1,*n)) {
266 *info = -9;
267 } else if (*ldx < max(1,*n)) {
268 *info = -11;
269 }
270 if (*info != 0) {
271 i__1 = -(*info);
272 xerbla_("ZSPSVX", &i__1, (ftnlen)6);
273 return 0;
274 }
275
276 if (nofact) {
277
278 /* Compute the factorization A = U*D*U' or A = L*D*L'. */
279
280 i__1 = *n * (*n + 1) / 2;
281 zcopy_(&i__1, &ap[1], &c__1, &afp[1], &c__1);
282 zsptrf_(uplo, n, &afp[1], &ipiv[1], info, (ftnlen)1);
283
284 /* Return if INFO is non-zero. */
285
286 if (*info != 0) {
287 if (*info > 0) {
288 *rcond = 0.;
289 }
290 return 0;
291 }
292 }
293
294 /* Compute the norm of the matrix A. */
295
296 anorm = zlansp_("I", uplo, n, &ap[1], &rwork[1], (ftnlen)1, (ftnlen)1);
297
298 /* Compute the reciprocal of the condition number of A. */
299
300 zspcon_(uplo, n, &afp[1], &ipiv[1], &anorm, rcond, &work[1], info, (
301 ftnlen)1);
302
303 /* Set INFO = N+1 if the matrix is singular to working precision. */
304
305 if (*rcond < dlamch_("Epsilon", (ftnlen)7)) {
306 *info = *n + 1;
307 }
308
309 /* Compute the solution vectors X. */
310
311 zlacpy_("Full", n, nrhs, &b[b_offset], ldb, &x[x_offset], ldx, (ftnlen)4);
312 zsptrs_(uplo, n, nrhs, &afp[1], &ipiv[1], &x[x_offset], ldx, info, (
313 ftnlen)1);
314
315 /* Use iterative refinement to improve the computed solutions and */
316 /* compute error bounds and backward error estimates for them. */
317
318 zsprfs_(uplo, n, nrhs, &ap[1], &afp[1], &ipiv[1], &b[b_offset], ldb, &x[
319 x_offset], ldx, &ferr[1], &berr[1], &work[1], &rwork[1], info, (
320 ftnlen)1);
321
322 return 0;
323
324 /* End of ZSPSVX */
325
326 } /* zspsvx_ */
327
328