From 22f703cab05b7cd368f4de9e03991b7664dc5022 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B6rg=20Frings-F=C3=BCrst?= Date: Mon, 1 Sep 2014 13:56:46 +0200 Subject: Initial import of argyll version 1.5.1-8 --- numlib/dnsqtest.c | 192 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 192 insertions(+) create mode 100644 numlib/dnsqtest.c (limited to 'numlib/dnsqtest.c') 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; +} + -- cgit v1.2.3