Skip to content
Snippets Groups Projects
Commit a3b27e37 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

Reorganize the Nuemric lib

parent 1f482884
No related branches found
No related tags found
No related merge requests found
# $Id: Makefile,v 1.3 2002-05-18 16:31:16 geuzaine Exp $
# $Id: Makefile,v 1.4 2002-05-18 21:34:29 geuzaine Exp $
#
# Makefile for "libGmshNumeric.a"
#
......@@ -19,16 +19,10 @@ RMFLAGS = -f
CFLAGS = $(OPT_FLAGS) $(OS_FLAGS) $(VERSION_FLAGS) $(INCLUDE)
SRC = Numeric.cpp\
mnbrak.cpp \
brent.cpp \
nrutil.cpp \
dsvdcmp.cpp \
newt.cpp \
fmin.cpp \
fdjac.cpp \
lnsrch.cpp \
lubksb.cpp \
ludcmp.cpp
Newton.cpp \
Opti.cpp \
SVD.cpp \
NRUtil.cpp \
OBJ = $(SRC:.cpp=.o)
......@@ -53,17 +47,14 @@ depend:
$(RM) $(RMFLAGS) Makefile.new
# DO NOT DELETE THIS LINE
mnbrak.o: mnbrak.cpp nrutil.h Numeric.h
brent.o: brent.cpp nrutil.h Numeric.h
nrutil.o: nrutil.cpp ../Common/Gmsh.h ../Common/Message.h \
Numeric.o: Numeric.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h
dsvdcmp.o: dsvdcmp.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/avl.h ../DataStr/Tools.h Numeric.h
Newton.o: Newton.cpp NRUtil.h Numeric.h
Opti.o: Opti.cpp NRUtil.h Numeric.h
SVD.o: SVD.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
NRUtil.h Numeric.h
NRUtil.o: NRUtil.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h nrutil.h Numeric.h
newt.o: newt.cpp nrutil.h Numeric.h
fmin.o: fmin.cpp nrutil.h Numeric.h
fdjac.o: fdjac.cpp nrutil.h Numeric.h
lnsrch.o: lnsrch.cpp nrutil.h Numeric.h
lubksb.o: lubksb.cpp
ludcmp.o: ludcmp.cpp nrutil.h Numeric.h
../DataStr/avl.h ../DataStr/Tools.h
// $Id: nrutil.cpp,v 1.1 2002-05-18 07:18:03 geuzaine Exp $
// $Id: NRUtil.cpp,v 1.1 2002-05-18 21:34:29 geuzaine Exp $
/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
......
File moved
// $Id: Newton.cpp,v 1.1 2002-05-18 21:34:29 geuzaine Exp $
/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
#include <math.h>
#define NRANSI
#include "NRUtil.h"
#define SQR(a) ((a)*(a))
// global variables
int nn;
double *fvec;
void (*nrfuncv) (int n, double v[], double f[]);
// fmin.cpp
double
fmin (double x[])
{
int i;
double sum;
(*nrfuncv) (nn, x, fvec);
for (sum = 0.0, i = 1; i <= nn; i++)
sum += SQR (fvec[i]);
return 0.5 * sum;
}
// newt.cpp
#define MAXITS 200
#define TOLF 1.0e-10 /* 1.0e-4 */
#define TOLMIN 1.0e-12 /* 1.0e-6 */
#define TOLX 1.0e-14 /* 1.0e-7 */
#define STPMX 100.0
#define FREERETURN {free_dvector(fvec,1,n);free_dvector(xold,1,n);\
free_dvector(p,1,n);free_dvector(g,1,n);free_dmatrix(fjac,1,n,1,n);\
free_ivector(indx,1,n);return;}
void
newt (double x[], int n, int *check,
void (*vecfunc) (int, double[], double[]))
{
void fdjac (int n, double x[], double fvec[], double **df,
void (*vecfunc) (int, double[], double[]));
double fmin (double x[]);
void lnsrch (int n, double xold[], double fold, double g[], double p[],
double x[], double *f, double stpmax, int *check,
double (*func) (double[]));
void lubksb (double **a, int n, int *indx, double b[]);
void ludcmp (double **a, int n, int *indx, double *d);
int i, its, j, *indx;
double d, den, f, fold, stpmax, sum, temp, test, **fjac, *g, *p, *xold;
indx = ivector (1, n);
fjac = dmatrix (1, n, 1, n);
g = dvector (1, n);
p = dvector (1, n);
xold = dvector (1, n);
fvec = dvector (1, n);
nn = n;
nrfuncv = vecfunc;
f = fmin (x);
test = 0.0;
for (i = 1; i <= n; i++)
if (fabs (fvec[i]) > test)
test = fabs (fvec[i]);
if (test < 0.01 * TOLF)
FREERETURN;
for (sum = 0.0, i = 1; i <= n; i++)
sum += SQR (x[i]);
stpmax = STPMX * MAX (sqrt (sum), (double) n);
for (its = 1; its <= MAXITS; its++) {
fdjac (n, x, fvec, fjac, vecfunc);
//vecjac (n, x, fvec, fjac, vecfunc);
for (i = 1; i <= n; i++) {
for (sum = 0.0, j = 1; j <= n; j++)
sum += fjac[j][i] * fvec[j];
g[i] = sum;
}
for (i = 1; i <= n; i++)
xold[i] = x[i];
fold = f;
for (i = 1; i <= n; i++)
p[i] = -fvec[i];
ludcmp (fjac, n, indx, &d);
lubksb (fjac, n, indx, p);
lnsrch (n, xold, fold, g, p, x, &f, stpmax, check, fmin);
test = 0.0;
for (i = 1; i <= n; i++)
if (fabs (fvec[i]) > test)
test = fabs (fvec[i]);
if (test < TOLF) {
*check = 0;
FREERETURN;
}
if (*check) {
test = 0.0;
den = MAX (f, 0.5 * n);
for (i = 1; i <= n; i++) {
temp = fabs (g[i]) * MAX (fabs (x[i]), 1.0) / den;
if (temp > test)
test = temp;
}
*check = (test < TOLMIN ? 1 : 0);
FREERETURN;
}
test = 0.0;
for (i = 1; i <= n; i++) {
temp = (fabs (x[i] - xold[i])) / MAX (fabs (x[i]), 1.0);
if (temp > test)
test = temp;
}
if (test < TOLX)
FREERETURN;
}
//Msg(WARNING, "MAXITS exceeded in newt");
*check=1;
}
#undef MAXITS
#undef TOLF
#undef TOLMIN
#undef TOLX
#undef STPMX
#undef FREERETURN
// lnsrch.cpp
#define ALF 1.0e-8 /* 1.0e-4 */
#define TOLX 1.0e-14 /* 1.0e-7 */
void
lnsrch (int n, double xold[], double fold, double g[], double p[], double x[],
double *f, double stpmax, int *check, double (*func) (double[]))
{
int i;
double a, alam, alam2, alamin, b, disc, f2, fold2, rhs1, rhs2, slope,
sum, temp, test, tmplam;
*check = 0;
for (sum = 0.0, i = 1; i <= n; i++)
sum += p[i] * p[i];
sum = sqrt (sum);
if (sum > stpmax)
for (i = 1; i <= n; i++)
p[i] *= stpmax / sum;
for (slope = 0.0, i = 1; i <= n; i++)
slope += g[i] * p[i];
test = 0.0;
for (i = 1; i <= n; i++) {
temp = fabs (p[i]) / MAX (fabs (xold[i]), 1.0);
if (temp > test)
test = temp;
}
alamin = TOLX / test;
alam = 1.0;
for (;;) {
for (i = 1; i <= n; i++)
x[i] = xold[i] + alam * p[i];
*f = (*func) (x);
if (alam < alamin) {
for (i = 1; i <= n; i++)
x[i] = xold[i];
*check = 1;
return;
}
else if (*f <= fold + ALF * alam * slope)
return;
else {
if (alam == 1.0)
tmplam = -slope / (2.0 * (*f - fold - slope));
else {
rhs1 = *f - fold - alam * slope;
rhs2 = f2 - fold2 - alam2 * slope;
a = (rhs1 / (alam * alam) - rhs2 / (alam2 * alam2)) / (alam - alam2);
b =
(-alam2 * rhs1 / (alam * alam) +
alam * rhs2 / (alam2 * alam2)) / (alam - alam2);
if (a == 0.0)
tmplam = -slope / (2.0 * b);
else {
disc = b * b - 3.0 * a * slope;
if (disc < 0.0)
nrerror ("Roundoff problem in lnsrch.");
else
tmplam = (-b + sqrt (disc)) / (3.0 * a);
}
if (tmplam > 0.5 * alam)
tmplam = 0.5 * alam;
}
}
alam2 = alam;
f2 = *f;
fold2 = fold;
alam = MAX (tmplam, 0.1 * alam);
}
}
#undef ALF
#undef TOLX
// fdjac.cpp
#define EPS 1.0e-4 /* 1.0e-4 */
void
fdjac (int n, double x[], double fvec[], double **df,
void (*vecfunc) (int, double[], double[]))
{
int i, j;
double h, temp, *f;
f = dvector (1, n);
for (j = 1; j <= n; j++) {
temp = x[j];
h = EPS * fabs (temp);
if (h == 0.0)
h = EPS;
x[j] = temp + h;
h = x[j] - temp;
(*vecfunc) (n, x, f);
x[j] = temp;
for (i = 1; i <= n; i++)
df[i][j] = (f[i] - fvec[i]) / h;
}
free_dvector (f, 1, n);
}
#undef EPS
// ludcmp.cpp
#define TINY 1.0e-20;
void
ludcmp (double **a, int n, int *indx, double *d)
{
int i, imax, j, k;
double big, dum, sum, temp;
double *vv;
vv = dvector (1, n);
*d = 1.0;
for (i = 1; i <= n; i++) {
big = 0.0;
for (j = 1; j <= n; j++)
if ((temp = fabs (a[i][j])) > big)
big = temp;
if (big == 0.0)
nrerror ("Singular matrix in routine ludcmp");
vv[i] = 1.0 / big;
}
for (j = 1; j <= n; j++) {
for (i = 1; i < j; i++) {
sum = a[i][j];
for (k = 1; k < i; k++)
sum -= a[i][k] * a[k][j];
a[i][j] = sum;
}
big = 0.0;
for (i = j; i <= n; i++) {
sum = a[i][j];
for (k = 1; k < j; k++)
sum -= a[i][k] * a[k][j];
a[i][j] = sum;
if ((dum = vv[i] * fabs (sum)) >= big) {
big = dum;
imax = i;
}
}
if (j != imax) {
for (k = 1; k <= n; k++) {
dum = a[imax][k];
a[imax][k] = a[j][k];
a[j][k] = dum;
}
*d = -(*d);
vv[imax] = vv[j];
}
indx[j] = imax;
if (a[j][j] == 0.0)
a[j][j] = TINY;
if (j != n) {
dum = 1.0 / (a[j][j]);
for (i = j + 1; i <= n; i++)
a[i][j] *= dum;
}
}
free_dvector (vv, 1, n);
}
#undef TINY
// lubksb.cpp
void
lubksb (double **a, int n, int *indx, double b[])
{
int i, ii = 0, ip, j;
double sum;
for (i = 1; i <= n; i++) {
ip = indx[i];
sum = b[ip];
b[ip] = b[i];
if (ii)
for (j = ii; j <= i - 1; j++)
sum -= a[i][j] * b[j];
else if (sum)
ii = i;
b[i] = sum;
}
for (i = n; i >= 1; i--) {
sum = b[i];
for (j = i + 1; j <= n; j++)
sum -= a[i][j] * b[j];
b[i] = sum / a[i][i];
}
}
// $Id: brent.cpp,v 1.1 2002-05-18 07:18:03 geuzaine Exp $
// $Id: Opti.cpp,v 1.1 2002-05-18 21:34:29 geuzaine Exp $
/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
#include <math.h>
#define NRANSI
#include "nrutil.h"
#include "NRUtil.h"
// mnbrak
#define GOLD 1.618034
#define GLIMIT 100.0
#define TINY 1.0e-20
#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
void
mnbrak (double *ax, double *bx, double *cx, double *fa, double *fb, double *fc,
double (*func) (double))
{
double ulim, u, r, q, fu, dum;
*fa = (*func) (*ax);
*fb = (*func) (*bx);
if (*fb > *fa)
{
SHFT (dum, *ax, *bx, dum)
SHFT (dum, *fb, *fa, dum)
}
*cx = (*bx) + GOLD * (*bx - *ax);
*fc = (*func) (*cx);
while (*fb > *fc)
{
r = (*bx - *ax) * (*fb - *fc);
q = (*bx - *cx) * (*fb - *fa);
u = (*bx) - ((*bx - *cx) * q - (*bx - *ax) * r) /
(2.0 * SIGN (FMAX (fabs (q - r), TINY), q - r));
ulim = (*bx) + GLIMIT * (*cx - *bx);
if ((*bx - u) * (u - *cx) > 0.0)
{
fu = (*func) (u);
if (fu < *fc)
{
*ax = (*bx);
*bx = u;
*fa = (*fb);
*fb = fu;
return;
}
else if (fu > *fb)
{
*cx = u;
*fc = fu;
return;
}
u = (*cx) + GOLD * (*cx - *bx);
fu = (*func) (u);
}
else if ((*cx - u) * (u - ulim) > 0.0)
{
fu = (*func) (u);
if (fu < *fc)
{
SHFT (*bx, *cx, u, *cx + GOLD * (*cx - *bx))
SHFT (*fb, *fc, fu, (*func) (u))
}
}
else if ((u - ulim) * (ulim - *cx) >= 0.0)
{
u = ulim;
fu = (*func) (u);
}
else
{
u = (*cx) + GOLD * (*cx - *bx);
fu = (*func) (u);
}
SHFT (*ax, *bx, *cx, u)
SHFT (*fa, *fb, *fc, fu)
}
}
#undef GOLD
#undef GLIMIT
#undef TINY
#undef SHFT
#undef NRANSI
// brent
#define ITMAX 100
#define CGOLD 0.3819660
#define ZEPS 1.0e-10
......@@ -96,4 +179,3 @@ brent (double ax, double bx, double cx, double (*f) (double), double tol,
#undef ZEPS
#undef SHFT
#undef NRANSI
/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
// $Id: dsvdcmp.cpp,v 1.1 2002-05-18 07:18:03 geuzaine Exp $
// $Id: SVD.cpp,v 1.1 2002-05-18 21:34:29 geuzaine Exp $
/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
#include "Gmsh.h"
#include "nrutil.h"
#include "NRUtil.h"
double dpythag(double a, double b)
{
......
// $Id: fdjac.cpp,v 1.1 2002-05-18 07:18:03 geuzaine Exp $
#include <math.h>
#define NRANSI
#include "nrutil.h"
#define EPS 1.0e-4
void
fdjac (int n, float x[], float fvec[], float **df,
void (*vecfunc) (int, float[], float[]))
{
int i, j;
float h, temp, *f;
f = vector (1, n);
for (j = 1; j <= n; j++)
{
temp = x[j];
h = EPS * fabs (temp);
if (h == 0.0)
h = EPS;
x[j] = temp + h;
h = x[j] - temp;
(*vecfunc) (n, x, f);
x[j] = temp;
for (i = 1; i <= n; i++)
df[i][j] = (f[i] - fvec[i]) / h;
}
free_vector (f, 1, n);
}
#undef EPS
#undef NRANSI
// $Id: fmin.cpp,v 1.1 2002-05-18 07:18:03 geuzaine Exp $
#define NRANSI
#include "nrutil.h"
extern int nn;
extern float *fvec;
extern void (*nrfuncv) (int n, float v[], float f[]);
float
fmin (float x[])
{
int i;
float sum;
(*nrfuncv) (nn, x, fvec);
for (sum = 0.0, i = 1; i <= nn; i++)
sum += SQR (fvec[i]);
return 0.5 * sum;
}
#undef NRANSI
// $Id: lnsrch.cpp,v 1.1 2002-05-18 07:18:03 geuzaine Exp $
#include <math.h>
#define NRANSI
#include "nrutil.h"
#define ALF 1.0e-4
#define TOLX 1.0e-7
void
lnsrch (int n, float xold[], float fold, float g[], float p[], float x[],
float *f, float stpmax, int *check, float (*func) (float[]))
{
int i;
float a, alam, alam2, alamin, b, disc, f2, fold2, rhs1, rhs2, slope,
sum, temp, test, tmplam;
*check = 0;
for (sum = 0.0, i = 1; i <= n; i++)
sum += p[i] * p[i];
sum = sqrt (sum);
if (sum > stpmax)
for (i = 1; i <= n; i++)
p[i] *= stpmax / sum;
for (slope = 0.0, i = 1; i <= n; i++)
slope += g[i] * p[i];
test = 0.0;
for (i = 1; i <= n; i++)
{
temp = fabs (p[i]) / FMAX (fabs (xold[i]), 1.0);
if (temp > test)
test = temp;
}
alamin = TOLX / test;
alam = 1.0;
for (;;)
{
for (i = 1; i <= n; i++)
x[i] = xold[i] + alam * p[i];
*f = (*func) (x);
if (alam < alamin)
{
for (i = 1; i <= n; i++)
x[i] = xold[i];
*check = 1;
return;
}
else if (*f <= fold + ALF * alam * slope)
return;
else
{
if (alam == 1.0)
tmplam = -slope / (2.0 * (*f - fold - slope));
else
{
rhs1 = *f - fold - alam * slope;
rhs2 = f2 - fold2 - alam2 * slope;
a = (rhs1 / (alam * alam) - rhs2 / (alam2 * alam2)) / (alam - alam2);
b = (-alam2 * rhs1 / (alam * alam) + alam * rhs2 / (alam2 * alam2)) / (alam - alam2);
if (a == 0.0)
tmplam = -slope / (2.0 * b);
else
{
disc = b * b - 3.0 * a * slope;
if (disc < 0.0)
nrerror ("Roundoff problem in lnsrch.");
else
tmplam = (-b + sqrt (disc)) / (3.0 * a);
}
if (tmplam > 0.5 * alam)
tmplam = 0.5 * alam;
}
}
alam2 = alam;
f2 = *f;
fold2 = fold;
alam = FMAX (tmplam, 0.1 * alam);
}
}
#undef ALF
#undef TOLX
#undef NRANSI
// $Id: lubksb.cpp,v 1.1 2002-05-18 07:18:03 geuzaine Exp $
void
lubksb (float **a, int n, int *indx, float b[])
{
int i, ii = 0, ip, j;
float sum;
for (i = 1; i <= n; i++)
{
ip = indx[i];
sum = b[ip];
b[ip] = b[i];
if (ii)
for (j = ii; j <= i - 1; j++)
sum -= a[i][j] * b[j];
else if (sum)
ii = i;
b[i] = sum;
}
for (i = n; i >= 1; i--)
{
sum = b[i];
for (j = i + 1; j <= n; j++)
sum -= a[i][j] * b[j];
b[i] = sum / a[i][i];
}
}
// $Id: ludcmp.cpp,v 1.1 2002-05-18 07:18:03 geuzaine Exp $
#include <math.h>
#define NRANSI
#include "nrutil.h"
#define TINY 1.0e-20;
void
ludcmp (float **a, int n, int *indx, float *d)
{
int i, imax, j, k;
float big, dum, sum, temp;
float *vv;
vv = vector (1, n);
*d = 1.0;
for (i = 1; i <= n; i++)
{
big = 0.0;
for (j = 1; j <= n; j++)
if ((temp = fabs (a[i][j])) > big)
big = temp;
if (big == 0.0)
nrerror ("Singular matrix in routine ludcmp");
vv[i] = 1.0 / big;
}
for (j = 1; j <= n; j++)
{
for (i = 1; i < j; i++)
{
sum = a[i][j];
for (k = 1; k < i; k++)
sum -= a[i][k] * a[k][j];
a[i][j] = sum;
}
big = 0.0;
for (i = j; i <= n; i++)
{
sum = a[i][j];
for (k = 1; k < j; k++)
sum -= a[i][k] * a[k][j];
a[i][j] = sum;
if ((dum = vv[i] * fabs (sum)) >= big)
{
big = dum;
imax = i;
}
}
if (j != imax)
{
for (k = 1; k <= n; k++)
{
dum = a[imax][k];
a[imax][k] = a[j][k];
a[j][k] = dum;
}
*d = -(*d);
vv[imax] = vv[j];
}
indx[j] = imax;
if (a[j][j] == 0.0)
a[j][j] = TINY;
if (j != n)
{
dum = 1.0 / (a[j][j]);
for (i = j + 1; i <= n; i++)
a[i][j] *= dum;
}
}
free_vector (vv, 1, n);
}
#undef TINY
#undef NRANSI
// $Id: mnbrak.cpp,v 1.1 2002-05-18 07:18:03 geuzaine Exp $
#include <math.h>
#define NRANSI
#include "nrutil.h"
#define GOLD 1.618034
#define GLIMIT 100.0
#define TINY 1.0e-20
#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
void
mnbrak (double *ax, double *bx, double *cx, double *fa, double *fb, double *fc,
double (*func) (double))
{
double ulim, u, r, q, fu, dum;
*fa = (*func) (*ax);
*fb = (*func) (*bx);
if (*fb > *fa)
{
SHFT (dum, *ax, *bx, dum)
SHFT (dum, *fb, *fa, dum)
}
*cx = (*bx) + GOLD * (*bx - *ax);
*fc = (*func) (*cx);
while (*fb > *fc)
{
r = (*bx - *ax) * (*fb - *fc);
q = (*bx - *cx) * (*fb - *fa);
u = (*bx) - ((*bx - *cx) * q - (*bx - *ax) * r) /
(2.0 * SIGN (FMAX (fabs (q - r), TINY), q - r));
ulim = (*bx) + GLIMIT * (*cx - *bx);
if ((*bx - u) * (u - *cx) > 0.0)
{
fu = (*func) (u);
if (fu < *fc)
{
*ax = (*bx);
*bx = u;
*fa = (*fb);
*fb = fu;
return;
}
else if (fu > *fb)
{
*cx = u;
*fc = fu;
return;
}
u = (*cx) + GOLD * (*cx - *bx);
fu = (*func) (u);
}
else if ((*cx - u) * (u - ulim) > 0.0)
{
fu = (*func) (u);
if (fu < *fc)
{
SHFT (*bx, *cx, u, *cx + GOLD * (*cx - *bx))
SHFT (*fb, *fc, fu, (*func) (u))
}
}
else if ((u - ulim) * (ulim - *cx) >= 0.0)
{
u = ulim;
fu = (*func) (u);
}
else
{
u = (*cx) + GOLD * (*cx - *bx);
fu = (*func) (u);
}
SHFT (*ax, *bx, *cx, u)
SHFT (*fa, *fb, *fc, fu)
}
}
#undef GOLD
#undef GLIMIT
#undef TINY
#undef SHFT
#undef NRANSI
// $Id: newt.cpp,v 1.1 2002-05-18 07:18:03 geuzaine Exp $
#include <math.h>
#define NRANSI
#include "nrutil.h"
#define MAXITS 200
#define TOLF 1.0e-4
#define TOLMIN 1.0e-6
#define TOLX 1.0e-7
#define STPMX 100.0
int nn;
float *fvec;
void (*nrfuncv) (int n, float v[], float f[]);
#define FREERETURN {free_vector(fvec,1,n);free_vector(xold,1,n);\
free_vector(p,1,n);free_vector(g,1,n);free_matrix(fjac,1,n,1,n);\
free_ivector(indx,1,n);return;}
void
newt (float x[], int n, int *check,
void (*vecfunc) (int, float[], float[]))
{
void fdjac (int n, float x[], float fvec[], float **df,
void (*vecfunc) (int, float[], float[]));
float fmin (float x[]);
void lnsrch (int n, float xold[], float fold, float g[], float p[], float x[],
float *f, float stpmax, int *check, float (*func) (float[]));
void lubksb (float **a, int n, int *indx, float b[]);
void ludcmp (float **a, int n, int *indx, float *d);
int i, its, j, *indx;
float d, den, f, fold, stpmax, sum, temp, test, **fjac, *g, *p, *xold;
indx = ivector (1, n);
fjac = matrix (1, n, 1, n);
g = vector (1, n);
p = vector (1, n);
xold = vector (1, n);
fvec = vector (1, n);
nn = n;
nrfuncv = vecfunc;
f = fmin (x);
test = 0.0;
for (i = 1; i <= n; i++)
if (fabs (fvec[i]) > test)
test = fabs (fvec[i]);
if (test < 0.01 * TOLF)
FREERETURN
for (sum = 0.0, i = 1; i <= n; i++)
sum += SQR (x[i]);
stpmax = STPMX * FMAX (sqrt (sum), (float) n);
for (its = 1; its <= MAXITS; its++)
{
fdjac (n, x, fvec, fjac, vecfunc);
for (i = 1; i <= n; i++)
{
for (sum = 0.0, j = 1; j <= n; j++)
sum += fjac[j][i] * fvec[j];
g[i] = sum;
}
for (i = 1; i <= n; i++)
xold[i] = x[i];
fold = f;
for (i = 1; i <= n; i++)
p[i] = -fvec[i];
ludcmp (fjac, n, indx, &d);
lubksb (fjac, n, indx, p);
lnsrch (n, xold, fold, g, p, x, &f, stpmax, check, fmin);
test = 0.0;
for (i = 1; i <= n; i++)
if (fabs (fvec[i]) > test)
test = fabs (fvec[i]);
if (test < TOLF)
{
*check = 0;
FREERETURN
}
if (*check)
{
test = 0.0;
den = FMAX (f, 0.5 * n);
for (i = 1; i <= n; i++)
{
temp = fabs (g[i]) * FMAX (fabs (x[i]), 1.0) / den;
if (temp > test)
test = temp;
}
*check = (test < TOLMIN ? 1 : 0);
FREERETURN
}
test = 0.0;
for (i = 1; i <= n; i++)
{
temp = (fabs (x[i] - xold[i])) / FMAX (fabs (x[i]), 1.0);
if (temp > test)
test = temp;
}
if (test < TOLX)
FREERETURN
}
nrerror ("MAXITS exceeded in newt");
}
#undef MAXITS
#undef TOLF
#undef TOLMIN
#undef TOLX
#undef STPMX
#undef FREERETURN
#undef NRANSI
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment