summaryrefslogtreecommitdiff
path: root/numlib/dnsqtest.c
diff options
context:
space:
mode:
Diffstat (limited to 'numlib/dnsqtest.c')
-rw-r--r--numlib/dnsqtest.c192
1 files changed, 192 insertions, 0 deletions
diff --git a/numlib/dnsqtest.c b/numlib/dnsqtest.c
new file mode 100644
index 0000000..3f02819
--- /dev/null
+++ b/numlib/dnsqtest.c
@@ -0,0 +1,192 @@
+
+/*
+ * Copyright 1999 Graeme W. Gill
+ * All rights reserved.
+ *
+ * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
+ * see the License.txt file for licencing details.
+ */
+
+/* Example use of dnsqe() */
+/* */
+/* The problem is to determine the values of X(1), X(2), ..., X(9), */
+/* which solve the system of tridiagonal equations */
+/* */
+/* (3-2*X(1))*X(1) -2*X(2) = -1 */
+/* -X(I-1) + (3-2*X(I))*X(I) -2*X(I+1) = -1, I=2-8 */
+/* -X(8) + (3-2*X(9))*X(9) = -1 */
+/* */
+/* Final approximate solution: */
+/* */
+/* -0.5706545E+00 */
+/* -0.6816283E+00 */
+/* -0.7017325E+00 */
+/* -0.7042129E+00 */
+/* -0.7013690E+00 */
+/* -0.6918656E+00 */
+/* -0.6657920E+00 */
+/* -0.5960342E+00 */
+/* -0.4164121E+00 */
+
+#include "numlib.h"
+
+/* Compute norm of a vector */
+static double denorm(int n, double *x);
+
+int fcn(void *fdata, int n, double *x, double *fvec, int iflag);
+
+double expect[9] = {
+ -0.5706545E+00,
+ -0.6816283E+00,
+ -0.7017325E+00,
+ -0.7042129E+00,
+ -0.7013690E+00,
+ -0.6918656E+00,
+ -0.6657920E+00,
+ -0.5960342E+00,
+ -0.4164121E+00 };
+
+int main(void)
+{
+ int n = 9 /* 9 */; /* Problem vector size */
+ double x[9]; /* Function input values */
+ double fvec[9]; /* Function output values */
+ double ss; /* Search area */
+ int info, j;
+ double fnorm;
+ int nprint = 0; /* Itteration debugging print = off */
+ double tol;
+
+ error_program = "dnsqtest"; /* Set global error reporting string */
+
+ /* Driver for dnsqe example. */
+ /* Not supplying Jacobian, use approximation */
+
+ /* The following starting values provide a rough solution. */
+ for (j = 1; j <= 9; ++j) {
+ x[j - 1] = -1.f;
+ }
+ ss = 0.1;
+
+ nprint = 0;
+
+ /* Set tol to the square root of the machine precision. */
+ /* Unless high precision solutions are required, */
+ /* this is the recommended setting. */
+
+ tol = sqrt(M_DIVER);
+
+ info = dnsqe(NULL, fcn, NULL, n, x, ss, fvec, tol, tol, 0, nprint);
+ fnorm = denorm(n, fvec);
+ fprintf(stdout,"Final L2 norm of the residuals = %e\n",fnorm);
+ fprintf(stdout,"Exit return value = %d (1 = sucess)\n",info);
+ fprintf(stdout,"Final approximate solution:\n");
+ for (j = 0; j < n; j++) {
+ fprintf(stdout,"x[%d] = %f, expect %f\n",j,x[j], expect[j]);
+ }
+
+ return 0;
+} /* main() */
+
+/* Function being solved */
+int fcn(void *fdata, int n, double *x, double *fvec, int iflag)
+{
+ double temp, temp1, temp2;
+ int k;
+
+ /* Function Body */
+ for (k = 0; k < n; ++k) {
+ temp = (3.0 - 2.0 * x[k]) * x[k];
+ temp1 = 0.0;
+ if (k != 0) {
+ temp1 = x[k-1];
+ }
+ temp2 = 0.0;
+ if (k != ((n)-1))
+ temp2 = x[k+1];
+ fvec[k] = temp - temp1 - 2.0 * temp2 + 1.0;
+ if (iflag == 0)
+ printf("x[%d] = %f, fvec[%d] + %f\n",k,x[k],k,fvec[k]);
+
+#ifdef DEBUG
+printf("~~ x[%d] = %f, fvec[%d] + %f\n",k,x[k],k,fvec[k]);
+#endif /* DEBUG */
+ }
+ /* Return < 0 to abort */
+ return 0;
+}
+
+/* - - - - - - - - - - - - - - - - - - - */
+
+static double denorm(
+ int n, /* Size of x[] */
+ double x[]) /* Input vector */
+{
+ /* Initialized data */
+ static double rdwarf = 3.834e-20;
+ static double rgiant = 1.304e19;
+
+ /* Local variables */
+ static double xabs, x1max, x3max;
+ static int i;
+ static double s1, s2, s3, agiant, floatn;
+ double ret_val, td;
+
+ s1 = 0.0; /* Large component */
+ s2 = 0.0; /* Intermedate component */
+ s3 = 0.0; /* Small component */
+ x1max = 0.0;
+ x3max = 0.0;
+ floatn = (double) (n + 1);
+ agiant = rgiant / floatn;
+ for (i = 0; i < n; i++) {
+ xabs = (td = x[i], fabs(td));
+
+ /* Sum for intermediate components. */
+ if (xabs > rdwarf && xabs < agiant) {
+ td = xabs; /* Computing 2nd power */
+ s2 += td * td;
+
+ /* Sum for small components. */
+ } else if (xabs <= rdwarf) {
+ if (xabs <= x3max) {
+ if (xabs != 0.0) { /* Computing 2nd power */
+ td = xabs / x3max;
+ s3 += td * td;
+ }
+ } else { /* Computing 2nd power */
+ td = x3max / xabs;
+ s3 = 1.0 + s3 * (td * td);
+ x3max = xabs;
+ }
+
+ /* Sum for large components. */
+ } else {
+ if (xabs <= x1max) { /* Computing 2nd power */
+ td = xabs / x1max;
+ s1 += td * td;
+ } else { /* Computing 2nd power */
+ td = x1max / xabs;
+ s1 = 1.0 + s1 * (td * td);
+ x1max = xabs;
+ }
+ }
+ }
+
+ /* Calculation of norm. */
+ if (s1 != 0.0) { /* Large is present */
+ ret_val = x1max * sqrt(s1 + s2 / x1max / x1max);
+ } else { /* Medium and small are present */
+ if (s2 == 0.0) {
+ ret_val = x3max * sqrt(s3); /* Small only */
+ } else {
+ if (s2 >= x3max) { /* Medium larger than small */
+ ret_val = sqrt(s2 * (1.0 + x3max / s2 * (x3max * s3)));
+ } else { /* Small large than medium */
+ ret_val = sqrt(x3max * (s2 / x3max + x3max * s3));
+ }
+ }
+ }
+ return ret_val;
+}
+