diff options
Diffstat (limited to 'numlib/tconjgrad.c')
-rwxr-xr-x | numlib/tconjgrad.c | 162 |
1 files changed, 162 insertions, 0 deletions
diff --git a/numlib/tconjgrad.c b/numlib/tconjgrad.c new file mode 100755 index 0000000..aa98871 --- /dev/null +++ b/numlib/tconjgrad.c @@ -0,0 +1,162 @@ +/* Code to test the conjgrad minimiser */ +/* + * Copyright 1999, 2018 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. + */ + +#include <stdio.h> +#include "numlib.h" + +/* Final approximate solution: */ + +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 }; + +double fcn( /* Return function value */ + void *fdata, /* Opaque data pointer */ + double tp[] /* Multivriate input value */ +); + +static double dfcn( + void *fdata, + double dp[], + double tp[] +); + +#define N 9 + +static void progress(void *pdata, int perc) { + printf("%c% 3d%%",cr_char,perc); + if (perc == 100) + printf("\n"); + fflush(stdout); +} + +int main(void) +{ + double cp[N]; /* Function input values */ + double s[N]; /* Search area */ + double err; + int j; + int nprint = 0; /* Itteration debugging print = off */ + int rc; + + error_program = "tpowell"; /* Set global error reporting string */ + check_if_not_interactive(); + + /* The following starting values provide a rough solution. */ + for (j = 0; j < N; j++) { + cp[j] = -1.f; + s[j] = 0.9; + } + + nprint = 0; + + /* Set tol to the square root of the machine precision. */ + /* Unless high precision solutions are required, */ + /* this is the recommended setting. */ + + rc = conjgrad( + &err, + N, /* Dimentionality */ + cp, /* Initial starting point */ + s, /* Size of initial search area */ + 0.00000001, /* Tollerance of error change to stop on */ + 1000, /* Maximum iterations allowed */ + fcn, /* Error function to evaluate */ + dfcn, /* Partial derivative function */ + NULL, /* Opaque data needed by function */ + progress, /* Progress callback */ + NULL /* Context for callback */ + ); + + + fprintf(stdout,"Status = %d, final approximate solution err = %f:\n",rc,err); + for (j = 0; j < N; j++) { + fprintf(stdout,"cp[%d] = %e, expect %e\n",j,cp[j],expect[j]); + } + + return 0; +} /* main() */ + +/* Function being minimized */ +double fcn( /* Return function value */ +void *fdata, /* Opaque data pointer */ +double tp[] /* Multivriate input value */ +) { + double err, tt; + double temp, temp1, temp2; + int k; + + /* Function Body */ + err = 0.0; + for (k = 0; k < N; ++k) { + temp = (3.0 - 2.0 * tp[k]) * tp[k]; + temp1 = 0.0; + if (k != 0) { + temp1 = tp[k-1]; + } + temp2 = 0.0; + if (k != ((N)-1)) + temp2 = tp[k+1]; + tt = temp - temp1 - 2.0 * temp2 + 1.0; + err += tt * tt; + } + err = sqrt(err); +//printf("Returning %16.14f\n",err); + return err; +} + + +/* Compute aprox. forward difference */ + +#define JEPS 1.0e-8 /* Aprox. sqrt of machine precision */ + +static double dfcn( +void *fdata, +double dp[], +double tp[] +) { + int i; + double d, h, temp; + + d = fcn(fdata, tp); + + for (i = 0; i < N; i++) { + temp = tp[i]; + + h = JEPS * fabs(temp); + if (h == 0.0) + h = JEPS; + tp[i] = temp + h; /* Add delta */ + h = tp[i] - temp; /* Actual delta with fp precision limits */ + + dp[i] = (fcn(fdata, tp) - d)/h; + + tp[i] = temp; /* Restore value */ + } + + return DFUNC_NRV; +} + + + + + + + + + + + |