dppsvx.c 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456
  1. /* dppsvx.f -- translated by f2c (version 20061008).
  2. You must link the resulting object file with libf2c:
  3. on Microsoft Windows system, link with libf2c.lib;
  4. on Linux or Unix systems, link with .../path/to/libf2c.a -lm
  5. or, if you install libf2c.a in a standard place, with -lf2c -lm
  6. -- in that order, at the end of the command line, as in
  7. cc *.o -lf2c -lm
  8. Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
  9. http://www.netlib.org/f2c/libf2c.zip
  10. */
  11. #include "f2c.h"
  12. #include "blaswrap.h"
  13. /* Table of constant values */
  14. static integer c__1 = 1;
  15. /* Subroutine */ int dppsvx_(char *fact, char *uplo, integer *n, integer *
  16. nrhs, doublereal *ap, doublereal *afp, char *equed, doublereal *s,
  17. doublereal *b, integer *ldb, doublereal *x, integer *ldx, doublereal *
  18. rcond, doublereal *ferr, doublereal *berr, doublereal *work, integer *
  19. iwork, integer *info)
  20. {
  21. /* System generated locals */
  22. integer b_dim1, b_offset, x_dim1, x_offset, i__1, i__2;
  23. doublereal d__1, d__2;
  24. /* Local variables */
  25. integer i__, j;
  26. doublereal amax, smin, smax;
  27. extern logical lsame_(char *, char *);
  28. doublereal scond, anorm;
  29. extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
  30. doublereal *, integer *);
  31. logical equil, rcequ;
  32. extern doublereal dlamch_(char *);
  33. logical nofact;
  34. extern /* Subroutine */ int dlacpy_(char *, integer *, integer *,
  35. doublereal *, integer *, doublereal *, integer *),
  36. xerbla_(char *, integer *);
  37. doublereal bignum;
  38. extern doublereal dlansp_(char *, char *, integer *, doublereal *,
  39. doublereal *);
  40. extern /* Subroutine */ int dppcon_(char *, integer *, doublereal *,
  41. doublereal *, doublereal *, doublereal *, integer *, integer *), dlaqsp_(char *, integer *, doublereal *, doublereal *,
  42. doublereal *, doublereal *, char *);
  43. integer infequ;
  44. extern /* Subroutine */ int dppequ_(char *, integer *, doublereal *,
  45. doublereal *, doublereal *, doublereal *, integer *),
  46. dpprfs_(char *, integer *, integer *, doublereal *, doublereal *,
  47. doublereal *, integer *, doublereal *, integer *, doublereal *,
  48. doublereal *, doublereal *, integer *, integer *),
  49. dpptrf_(char *, integer *, doublereal *, integer *);
  50. doublereal smlnum;
  51. extern /* Subroutine */ int dpptrs_(char *, integer *, integer *,
  52. doublereal *, doublereal *, integer *, integer *);
  53. /* -- LAPACK driver routine (version 3.2) -- */
  54. /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  55. /* November 2006 */
  56. /* .. Scalar Arguments .. */
  57. /* .. */
  58. /* .. Array Arguments .. */
  59. /* .. */
  60. /* Purpose */
  61. /* ======= */
  62. /* DPPSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to */
  63. /* compute the solution to a real system of linear equations */
  64. /* A * X = B, */
  65. /* where A is an N-by-N symmetric positive definite matrix stored in */
  66. /* packed format and X and B are N-by-NRHS matrices. */
  67. /* Error bounds on the solution and a condition estimate are also */
  68. /* provided. */
  69. /* Description */
  70. /* =========== */
  71. /* The following steps are performed: */
  72. /* 1. If FACT = 'E', real scaling factors are computed to equilibrate */
  73. /* the system: */
  74. /* diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B */
  75. /* Whether or not the system will be equilibrated depends on the */
  76. /* scaling of the matrix A, but if equilibration is used, A is */
  77. /* overwritten by diag(S)*A*diag(S) and B by diag(S)*B. */
  78. /* 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to */
  79. /* factor the matrix A (after equilibration if FACT = 'E') as */
  80. /* A = U**T* U, if UPLO = 'U', or */
  81. /* A = L * L**T, if UPLO = 'L', */
  82. /* where U is an upper triangular matrix and L is a lower triangular */
  83. /* matrix. */
  84. /* 3. If the leading i-by-i principal minor is not positive definite, */
  85. /* then the routine returns with INFO = i. Otherwise, the factored */
  86. /* form of A is used to estimate the condition number of the matrix */
  87. /* A. If the reciprocal of the condition number is less than machine */
  88. /* precision, INFO = N+1 is returned as a warning, but the routine */
  89. /* still goes on to solve for X and compute error bounds as */
  90. /* described below. */
  91. /* 4. The system of equations is solved for X using the factored form */
  92. /* of A. */
  93. /* 5. Iterative refinement is applied to improve the computed solution */
  94. /* matrix and calculate error bounds and backward error estimates */
  95. /* for it. */
  96. /* 6. If equilibration was used, the matrix X is premultiplied by */
  97. /* diag(S) so that it solves the original system before */
  98. /* equilibration. */
  99. /* Arguments */
  100. /* ========= */
  101. /* FACT (input) CHARACTER*1 */
  102. /* Specifies whether or not the factored form of the matrix A is */
  103. /* supplied on entry, and if not, whether the matrix A should be */
  104. /* equilibrated before it is factored. */
  105. /* = 'F': On entry, AFP contains the factored form of A. */
  106. /* If EQUED = 'Y', the matrix A has been equilibrated */
  107. /* with scaling factors given by S. AP and AFP will not */
  108. /* be modified. */
  109. /* = 'N': The matrix A will be copied to AFP and factored. */
  110. /* = 'E': The matrix A will be equilibrated if necessary, then */
  111. /* copied to AFP and factored. */
  112. /* UPLO (input) CHARACTER*1 */
  113. /* = 'U': Upper triangle of A is stored; */
  114. /* = 'L': Lower triangle of A is stored. */
  115. /* N (input) INTEGER */
  116. /* The number of linear equations, i.e., the order of the */
  117. /* matrix A. N >= 0. */
  118. /* NRHS (input) INTEGER */
  119. /* The number of right hand sides, i.e., the number of columns */
  120. /* of the matrices B and X. NRHS >= 0. */
  121. /* AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2) */
  122. /* On entry, the upper or lower triangle of the symmetric matrix */
  123. /* A, packed columnwise in a linear array, except if FACT = 'F' */
  124. /* and EQUED = 'Y', then A must contain the equilibrated matrix */
  125. /* diag(S)*A*diag(S). The j-th column of A is stored in the */
  126. /* array AP as follows: */
  127. /* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */
  128. /* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. */
  129. /* See below for further details. A is not modified if */
  130. /* FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit. */
  131. /* On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by */
  132. /* diag(S)*A*diag(S). */
  133. /* AFP (input or output) DOUBLE PRECISION array, dimension */
  134. /* (N*(N+1)/2) */
  135. /* If FACT = 'F', then AFP is an input argument and on entry */
  136. /* contains the triangular factor U or L from the Cholesky */
  137. /* factorization A = U'*U or A = L*L', in the same storage */
  138. /* format as A. If EQUED .ne. 'N', then AFP is the factored */
  139. /* form of the equilibrated matrix A. */
  140. /* If FACT = 'N', then AFP is an output argument and on exit */
  141. /* returns the triangular factor U or L from the Cholesky */
  142. /* factorization A = U'*U or A = L*L' of the original matrix A. */
  143. /* If FACT = 'E', then AFP is an output argument and on exit */
  144. /* returns the triangular factor U or L from the Cholesky */
  145. /* factorization A = U'*U or A = L*L' of the equilibrated */
  146. /* matrix A (see the description of AP for the form of the */
  147. /* equilibrated matrix). */
  148. /* EQUED (input or output) CHARACTER*1 */
  149. /* Specifies the form of equilibration that was done. */
  150. /* = 'N': No equilibration (always true if FACT = 'N'). */
  151. /* = 'Y': Equilibration was done, i.e., A has been replaced by */
  152. /* diag(S) * A * diag(S). */
  153. /* EQUED is an input argument if FACT = 'F'; otherwise, it is an */
  154. /* output argument. */
  155. /* S (input or output) DOUBLE PRECISION array, dimension (N) */
  156. /* The scale factors for A; not accessed if EQUED = 'N'. S is */
  157. /* an input argument if FACT = 'F'; otherwise, S is an output */
  158. /* argument. If FACT = 'F' and EQUED = 'Y', each element of S */
  159. /* must be positive. */
  160. /* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) */
  161. /* On entry, the N-by-NRHS right hand side matrix B. */
  162. /* On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y', */
  163. /* B is overwritten by diag(S) * B. */
  164. /* LDB (input) INTEGER */
  165. /* The leading dimension of the array B. LDB >= max(1,N). */
  166. /* X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) */
  167. /* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to */
  168. /* the original system of equations. Note that if EQUED = 'Y', */
  169. /* A and B are modified on exit, and the solution to the */
  170. /* equilibrated system is inv(diag(S))*X. */
  171. /* LDX (input) INTEGER */
  172. /* The leading dimension of the array X. LDX >= max(1,N). */
  173. /* RCOND (output) DOUBLE PRECISION */
  174. /* The estimate of the reciprocal condition number of the matrix */
  175. /* A after equilibration (if done). If RCOND is less than the */
  176. /* machine precision (in particular, if RCOND = 0), the matrix */
  177. /* is singular to working precision. This condition is */
  178. /* indicated by a return code of INFO > 0. */
  179. /* FERR (output) DOUBLE PRECISION array, dimension (NRHS) */
  180. /* The estimated forward error bound for each solution vector */
  181. /* X(j) (the j-th column of the solution matrix X). */
  182. /* If XTRUE is the true solution corresponding to X(j), FERR(j) */
  183. /* is an estimated upper bound for the magnitude of the largest */
  184. /* element in (X(j) - XTRUE) divided by the magnitude of the */
  185. /* largest element in X(j). The estimate is as reliable as */
  186. /* the estimate for RCOND, and is almost always a slight */
  187. /* overestimate of the true error. */
  188. /* BERR (output) DOUBLE PRECISION array, dimension (NRHS) */
  189. /* The componentwise relative backward error of each solution */
  190. /* vector X(j) (i.e., the smallest relative change in */
  191. /* any element of A or B that makes X(j) an exact solution). */
  192. /* WORK (workspace) DOUBLE PRECISION array, dimension (3*N) */
  193. /* IWORK (workspace) INTEGER array, dimension (N) */
  194. /* INFO (output) INTEGER */
  195. /* = 0: successful exit */
  196. /* < 0: if INFO = -i, the i-th argument had an illegal value */
  197. /* > 0: if INFO = i, and i is */
  198. /* <= N: the leading minor of order i of A is */
  199. /* not positive definite, so the factorization */
  200. /* could not be completed, and the solution has not */
  201. /* been computed. RCOND = 0 is returned. */
  202. /* = N+1: U is nonsingular, but RCOND is less than machine */
  203. /* precision, meaning that the matrix is singular */
  204. /* to working precision. Nevertheless, the */
  205. /* solution and error bounds are computed because */
  206. /* there are a number of situations where the */
  207. /* computed solution can be more accurate than the */
  208. /* value of RCOND would suggest. */
  209. /* Further Details */
  210. /* =============== */
  211. /* The packed storage scheme is illustrated by the following example */
  212. /* when N = 4, UPLO = 'U': */
  213. /* Two-dimensional storage of the symmetric matrix A: */
  214. /* a11 a12 a13 a14 */
  215. /* a22 a23 a24 */
  216. /* a33 a34 (aij = conjg(aji)) */
  217. /* a44 */
  218. /* Packed storage of the upper triangle of A: */
  219. /* AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ] */
  220. /* ===================================================================== */
  221. /* .. Parameters .. */
  222. /* .. */
  223. /* .. Local Scalars .. */
  224. /* .. */
  225. /* .. External Functions .. */
  226. /* .. */
  227. /* .. External Subroutines .. */
  228. /* .. */
  229. /* .. Intrinsic Functions .. */
  230. /* .. */
  231. /* .. Executable Statements .. */
  232. /* Parameter adjustments */
  233. --ap;
  234. --afp;
  235. --s;
  236. b_dim1 = *ldb;
  237. b_offset = 1 + b_dim1;
  238. b -= b_offset;
  239. x_dim1 = *ldx;
  240. x_offset = 1 + x_dim1;
  241. x -= x_offset;
  242. --ferr;
  243. --berr;
  244. --work;
  245. --iwork;
  246. /* Function Body */
  247. *info = 0;
  248. nofact = lsame_(fact, "N");
  249. equil = lsame_(fact, "E");
  250. if (nofact || equil) {
  251. *(unsigned char *)equed = 'N';
  252. rcequ = FALSE_;
  253. } else {
  254. rcequ = lsame_(equed, "Y");
  255. smlnum = dlamch_("Safe minimum");
  256. bignum = 1. / smlnum;
  257. }
  258. /* Test the input parameters. */
  259. if (! nofact && ! equil && ! lsame_(fact, "F")) {
  260. *info = -1;
  261. } else if (! lsame_(uplo, "U") && ! lsame_(uplo,
  262. "L")) {
  263. *info = -2;
  264. } else if (*n < 0) {
  265. *info = -3;
  266. } else if (*nrhs < 0) {
  267. *info = -4;
  268. } else if (lsame_(fact, "F") && ! (rcequ || lsame_(
  269. equed, "N"))) {
  270. *info = -7;
  271. } else {
  272. if (rcequ) {
  273. smin = bignum;
  274. smax = 0.;
  275. i__1 = *n;
  276. for (j = 1; j <= i__1; ++j) {
  277. /* Computing MIN */
  278. d__1 = smin, d__2 = s[j];
  279. smin = min(d__1,d__2);
  280. /* Computing MAX */
  281. d__1 = smax, d__2 = s[j];
  282. smax = max(d__1,d__2);
  283. /* L10: */
  284. }
  285. if (smin <= 0.) {
  286. *info = -8;
  287. } else if (*n > 0) {
  288. scond = max(smin,smlnum) / min(smax,bignum);
  289. } else {
  290. scond = 1.;
  291. }
  292. }
  293. if (*info == 0) {
  294. if (*ldb < max(1,*n)) {
  295. *info = -10;
  296. } else if (*ldx < max(1,*n)) {
  297. *info = -12;
  298. }
  299. }
  300. }
  301. if (*info != 0) {
  302. i__1 = -(*info);
  303. xerbla_("DPPSVX", &i__1);
  304. return 0;
  305. }
  306. if (equil) {
  307. /* Compute row and column scalings to equilibrate the matrix A. */
  308. dppequ_(uplo, n, &ap[1], &s[1], &scond, &amax, &infequ);
  309. if (infequ == 0) {
  310. /* Equilibrate the matrix. */
  311. dlaqsp_(uplo, n, &ap[1], &s[1], &scond, &amax, equed);
  312. rcequ = lsame_(equed, "Y");
  313. }
  314. }
  315. /* Scale the right-hand side. */
  316. if (rcequ) {
  317. i__1 = *nrhs;
  318. for (j = 1; j <= i__1; ++j) {
  319. i__2 = *n;
  320. for (i__ = 1; i__ <= i__2; ++i__) {
  321. b[i__ + j * b_dim1] = s[i__] * b[i__ + j * b_dim1];
  322. /* L20: */
  323. }
  324. /* L30: */
  325. }
  326. }
  327. if (nofact || equil) {
  328. /* Compute the Cholesky factorization A = U'*U or A = L*L'. */
  329. i__1 = *n * (*n + 1) / 2;
  330. dcopy_(&i__1, &ap[1], &c__1, &afp[1], &c__1);
  331. dpptrf_(uplo, n, &afp[1], info);
  332. /* Return if INFO is non-zero. */
  333. if (*info > 0) {
  334. *rcond = 0.;
  335. return 0;
  336. }
  337. }
  338. /* Compute the norm of the matrix A. */
  339. anorm = dlansp_("I", uplo, n, &ap[1], &work[1]);
  340. /* Compute the reciprocal of the condition number of A. */
  341. dppcon_(uplo, n, &afp[1], &anorm, rcond, &work[1], &iwork[1], info);
  342. /* Compute the solution matrix X. */
  343. dlacpy_("Full", n, nrhs, &b[b_offset], ldb, &x[x_offset], ldx);
  344. dpptrs_(uplo, n, nrhs, &afp[1], &x[x_offset], ldx, info);
  345. /* Use iterative refinement to improve the computed solution and */
  346. /* compute error bounds and backward error estimates for it. */
  347. dpprfs_(uplo, n, nrhs, &ap[1], &afp[1], &b[b_offset], ldb, &x[x_offset],
  348. ldx, &ferr[1], &berr[1], &work[1], &iwork[1], info);
  349. /* Transform the solution matrix X to a solution of the original */
  350. /* system. */
  351. if (rcequ) {
  352. i__1 = *nrhs;
  353. for (j = 1; j <= i__1; ++j) {
  354. i__2 = *n;
  355. for (i__ = 1; i__ <= i__2; ++i__) {
  356. x[i__ + j * x_dim1] = s[i__] * x[i__ + j * x_dim1];
  357. /* L40: */
  358. }
  359. /* L50: */
  360. }
  361. i__1 = *nrhs;
  362. for (j = 1; j <= i__1; ++j) {
  363. ferr[j] /= scond;
  364. /* L60: */
  365. }
  366. }
  367. /* Set INFO = N+1 if the matrix is singular to working precision. */
  368. if (*rcond < dlamch_("Epsilon")) {
  369. *info = *n + 1;
  370. }
  371. return 0;
  372. /* End of DPPSVX */
  373. } /* dppsvx_ */