diff --git a/Adapt/Adapt.cpp b/Adapt/Adapt.cpp deleted file mode 100644 index bf7224fec9d6783a1a32fc2fe662b05113e57f7c..0000000000000000000000000000000000000000 --- a/Adapt/Adapt.cpp +++ /dev/null @@ -1,191 +0,0 @@ -// $Id: Adapt.cpp,v 1.7 2001-04-26 17:58:00 remacle Exp $ - -#include "Gmsh.h" -#include "Adapt.h" -#include "nrutil.h" - -#define TOL 1.e-08 -#define MAXDEG 999 - -static int NN ; -static double MINH , *ERR , *HH , *PP , E0, DIM ; - -/* ------------------------------------------------------------------------ */ -/* f XXX */ -/* ------------------------------------------------------------------------ */ - -/* h-type version 1 : minimize the number of elements while keeping a - given global error */ - -double fH1 (double l){ - int i; - double val1 = 0.0, val2 = 0.0; - - for(i = 1 ; i <= NN ; i++){ - val1 += pow(2.*l*DSQR(ERR[i])*PP[i]/DIM, DIM/(2.*PP[i]+DIM)); - val2 += DSQR(ERR[i]) * pow(2.*l*DSQR(ERR[i])*PP[i]/DIM, -2.*PP[i]/(2.*PP[i]+DIM)); - } - - return( -(val1 + l * (val2 - DSQR(E0))) ); -} - -/* h-type version 2 : minimize the error while keeping a given number - of elements */ - -double fH2 (double l){ - int i; - double val1 = 0.0, val2 = 0.0, qi; - - for(i = 1 ; i <= NN ; i++){ - qi = pow(DIM*l/(2.*PP[i]*DSQR(ERR[i])), -DIM/(DIM+2.*PP[i])); - val1 += DSQR(ERR[i]) * pow(qi, -2.*PP[i]/DIM); - val2 += qi; - } - - return( -(val1 + l * (val2 - E0)) ); -} - -/* p-type : minimize error by modifying the interpolation order vector */ - -double fP1 (double l){ - int i; - double val1 = 0.0, val2 = 0.0, qi, e; - - for(i = 1 ; i <= NN ; i++){ - e = ERR[i]; - if(e == 0.0) e=1.e-12; - qi = - log(2.*l*log(HH[i]/MINH)*DSQR(e)) / log(HH[i]/MINH); - val1 -= 0.5 * qi; - val2 += pow(HH[i]/MINH, qi) * DSQR(e); - } - - return( -(val1 + l * (val2 - DSQR(E0))) ); -} - - -/* ------------------------------------------------------------------------ */ -/* A d a p t */ -/* ------------------------------------------------------------------------ */ - -double min1d (int method, double (*funct)(double), double *xmin){ - double xx, fx, fb, fa, bx, ax; - double brent(double ax, double bx, double cx, - double (*f)(double), double tol, double *xmin); - void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, - double *fc, double (*func)(double)); - - switch(method){ - case ADAPT_H1: - case ADAPT_P1: ax=1.e-12; xx=1.e2; break; - default: ax=1.e-15; xx=1.e-12; break; - } - mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,funct); - - return( brent(ax,xx,bx,funct,TOL,xmin) ); -} - -/* Adapt return the constraint (N0 ou e0) for the optimzed problem */ - -double AdaptMesh (int N, /* Number of elements */ - int method, /* ADAPT_H1, ADAPT_H2, ADAPT_P1 or ADAPT_P2 */ - int dim, /* 2 or 3 */ - double *err, /* elementary errors */ - double *h, /* elementary mesh sizes */ - double *p, /* elementary exponents */ - double e0 /* prescribed error or number of elements */){ - int i; - double contr=0.0, pivrai, lambda, minf, qi, ri, pi, obj, obj2, minri, maxri; - double errmin, errmax; - - h[N+1] = 1.0; - p[N+1] = 1.0; - - NN = N; - ERR = err; - HH = h; - PP = p; - E0 = e0; - DIM = (double)dim; - - for(i = 1 ; i <= N ; i++){ - if(i == 1) - errmin = errmax = err[i]; - else{ - errmin = DMIN(errmin, err[i]); - errmax = DMAX(errmax, err[i]); - } - } - - switch (method) { - - case ADAPT_H1 : - minf = min1d (method, fH1, &lambda); - obj = 0.0; - for(i = 1 ; i <= N ; i++){ - qi = pow(2.*lambda*DSQR(err[i])*p[i]/DIM, DIM/(2.*p[i]+DIM)); - ri = pow(qi,1./DIM); - if(i==1) minri = maxri = ri; - if(err[i] == 0.0) ri = .5; - minri = DMIN(minri, ri); - maxri = DMAX(maxri, ri); - obj += DSQR(err[i]) * pow(ri, -2.*p[i]) ; - h[i] = sqrt(2.) * h[i]/ri; - p[i] = ri; - } - contr = fabs(minf); - - printf("H-Refinement 1, Error %g=>%g, Objective %g, Reduction Factor %g->%g", - e0, sqrt(obj), -minf, minri, maxri); - break; - - case ADAPT_H2 : - minf = min1d (method, fH2, &lambda); - obj = 0.0; - for(i = 1 ; i <= N ; i++){ - qi = pow((DIM*lambda)/(2.*DSQR(err[i])*p[i]), -DIM/(DIM+2.*p[i])); - ri = pow(qi, 1./DIM); - if(i == 1) minri = maxri = ri; - minri = DMIN(minri, ri); - maxri = DMAX(maxri, ri); - obj += pow(ri,DIM) ; - h[i] = h[i]/ri; - p[i] = p[i]; - } - contr = sqrt(fabs(minf)); - - printf( "H-Refinement 2, Elements %g=>%g, Objective %g, Reduction Factor %g->%g", - e0, obj, 100. * sqrt(fabs(minf)), minri, maxri); - break; - - case ADAPT_P1 : - MINH = h[1]; - for(i = 1 ; i <= N ; i++) MINH =DMIN(h[i], MINH); - MINH /= 2.; - - minf = min1d (method, fP1, &lambda); - obj = obj2 = 0.0; - for(i = 1 ; i <= N ; i++){ - qi = -log(2.*lambda*DSQR(err[i])*log(h[i]/MINH)) / log(h[i]/MINH); - pi = p[i] - .5 * qi; - pivrai = DMIN(DMAX(1., (double)(int)(pi+.99)), MAXDEG); - obj2 += pow(h[i]/MINH, 2.*(p[i]-pivrai)) * DSQR(err[i]); - obj += DSQR(err[i]) * pow(h[i]/MINH, qi); - h[i] = h[i]; - p[i] = pi; - } - contr = fabs(minf); - - printf("P-Refinement, Error %g=%g=>%g, Objective %g", - e0, sqrt(obj), sqrt(obj2), minf); - break; - - case ADAPT_P2 : - minf = min1d (method, fH1, &lambda); - break; - - default : - printf("Unknown Adaption Method"); - } - - return(contr) ; -} diff --git a/Adapt/Adapt.h b/Adapt/Adapt.h deleted file mode 100644 index d4885075d44088c640e7d4d589ed26e919aaa702..0000000000000000000000000000000000000000 --- a/Adapt/Adapt.h +++ /dev/null @@ -1,17 +0,0 @@ -#ifndef _ADAPT_H_ -#define _ADAPT_H_ - -#define ADAPT_P1 1 -#define ADAPT_P2 2 -#define ADAPT_H1 3 -#define ADAPT_H2 4 - -double AdaptMesh (int N, /* Number of elements */ - int method, /* ADAPT_H1, ADAPT_H2, ADAPT_P1 or ADAPT_P2 */ - int dim, /* 2 or 3 */ - double *err, /* elementary errors */ - double *h, /* elementary mesh sizes */ - double *p, /* elementary exponents */ - double e0 /* prescribed error or number of elements */); - -#endif diff --git a/Adapt/Makefile b/Adapt/Makefile deleted file mode 100644 index 798e5fd0b1932fe3bff62efad32c3b79e5be0535..0000000000000000000000000000000000000000 --- a/Adapt/Makefile +++ /dev/null @@ -1,78 +0,0 @@ -# $Id: Makefile,v 1.25 2002-05-18 00:55:14 geuzaine Exp $ -# -# Makefile for "libGmshAdapt.a" -# - -.IGNORE: - -CXX = c++ -AR = ar ruvs -RM = rm -RANLIB = ranlib - -LIB = ../lib/libGmshAdapt.a -INCLUDE = -I../Common -I../DataStr - -OPT_FLAGS = -g -Wall -OS_FLAGS = -VERSION_FLAGS = - -RMFLAGS = -f -CFLAGS = $(OPT_FLAGS) $(OS_FLAGS) $(VERSION_FLAGS) $(INCLUDE) - -SRC = Adapt.cpp \ - mnbrak.cpp \ - brent.cpp \ - nrutil.cpp \ - dsvdcmp.cpp \ - newt.cpp \ - fmin.cpp \ - fdjac.cpp \ - lnsrch.cpp \ - lubksb.cpp \ - ludcmp.cpp - -OBJ = $(SRC:.cpp=.o) - -.SUFFIXES: .o .cpp - -$(LIB): $(OBJ) - $(AR) $(LIB) $(OBJ) - $(RANLIB) $(LIB) - -.cpp.o: - $(CXX) $(CFLAGS) -c $< - -clean: - $(RM) $(RMFLAGS) *.o - -lint: - $(LINT) $(CFLAGS) $(SRC) - -depend: - (sed '/^# DO NOT DELETE THIS LINE/q' Makefile && \ - $(CXX) -MM $(CFLAGS) ${SRC} \ - ) >Makefile.new - cp Makefile Makefile.bak - cp Makefile.new Makefile - $(RM) $(RMFLAGS) Makefile.new - -# DO NOT DELETE THIS LINE -Adapt.o: Adapt.cpp ../Common/Gmsh.h ../Common/Message.h \ - ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \ - ../DataStr/avl.h ../DataStr/Tools.h Adapt.h nrutil.h \ - ../Common/Numeric.h -mnbrak.o: mnbrak.cpp nrutil.h ../Common/Numeric.h -brent.o: brent.cpp nrutil.h ../Common/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 -dsvdcmp.o: dsvdcmp.cpp ../Common/Gmsh.h ../Common/Message.h \ - ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \ - ../DataStr/avl.h ../DataStr/Tools.h nrutil.h ../Common/Numeric.h -newt.o: newt.cpp nrutil.h ../Common/Numeric.h -fmin.o: fmin.cpp nrutil.h ../Common/Numeric.h -fdjac.o: fdjac.cpp nrutil.h ../Common/Numeric.h -lnsrch.o: lnsrch.cpp nrutil.h ../Common/Numeric.h -lubksb.o: lubksb.cpp -ludcmp.o: ludcmp.cpp nrutil.h ../Common/Numeric.h diff --git a/Adapt/brent.cpp b/Adapt/brent.cpp deleted file mode 100644 index 706d32781079a58085f47aa224e23000fe1eae3d..0000000000000000000000000000000000000000 --- a/Adapt/brent.cpp +++ /dev/null @@ -1,99 +0,0 @@ -// $Id: brent.cpp,v 1.4 2001-01-08 08:05:39 geuzaine Exp $ - - -#include <math.h> -#define NRANSI -#include "nrutil.h" -#define ITMAX 100 -#define CGOLD 0.3819660 -#define ZEPS 1.0e-10 -#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d); - -double -brent (double ax, double bx, double cx, double (*f) (double), double tol, - double *xmin) -{ - int iter; - double a, b, d=0.0, etemp, fu, fv, fw, fx, p, q, r, tol1, tol2, u, v, w, - x, xm; - double e = 0.0; - - a = (ax < cx ? ax : cx); - b = (ax > cx ? ax : cx); - x = w = v = bx; - fw = fv = fx = (*f) (x); - for (iter = 1; iter <= ITMAX; iter++) - { - xm = 0.5 * (a + b); - tol2 = 2.0 * (tol1 = tol * fabs (x) + ZEPS); - if (fabs (x - xm) <= (tol2 - 0.5 * (b - a))) - { - *xmin = x; - return fx; - } - if (fabs (e) > tol1) - { - r = (x - w) * (fx - fv); - q = (x - v) * (fx - fw); - p = (x - v) * q - (x - w) * r; - q = 2.0 * (q - r); - if (q > 0.0) - p = -p; - q = fabs (q); - etemp = e; - e = d; - if (fabs (p) >= fabs (0.5 * q * etemp) || p <= q * (a - x) || p >= q * (b - x)) - d = CGOLD * (e = (x >= xm ? a - x : b - x)); - else - { - d = p / q; - u = x + d; - if (u - a < tol2 || b - u < tol2) - d = SIGN (tol1, xm - x); - } - } - else - { - d = CGOLD * (e = (x >= xm ? a - x : b - x)); - } - u = (fabs (d) >= tol1 ? x + d : x + SIGN (tol1, d)); - fu = (*f) (u); - if (fu <= fx) - { - if (u >= x) - a = x; - else - b = x; - SHFT (v, w, x, u) - SHFT (fv, fw, fx, fu) - } - else - { - if (u < x) - a = u; - else - b = u; - if (fu <= fw || w == x) - { - v = w; - w = u; - fv = fw; - fw = fu; - } - else if (fu <= fv || v == x || v == w) - { - v = u; - fv = fu; - } - } - } - nrerror ("Too many iterations in brent"); - *xmin = x; - return fx; -} -#undef ITMAX -#undef CGOLD -#undef ZEPS -#undef SHFT -#undef NRANSI -/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */ diff --git a/Adapt/dsvdcmp.cpp b/Adapt/dsvdcmp.cpp deleted file mode 100644 index 3a6c849a3611307efce2d59ff49851e7994b2f0b..0000000000000000000000000000000000000000 --- a/Adapt/dsvdcmp.cpp +++ /dev/null @@ -1,215 +0,0 @@ -// $Id: dsvdcmp.cpp,v 1.5 2001-11-01 09:36:59 geuzaine Exp $ - -#include "Gmsh.h" -#include "nrutil.h" - -double dpythag(double a, double b) -{ - double absa,absb; - absa=fabs(a); - absb=fabs(b); - if (absa > absb) return absa*sqrt(1.0+SQR(absb/absa)); - else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+SQR(absa/absb))); -} - -void dsvdcmp(double **a, int m, int n, double w[], double **v) -{ - double dpythag(double a, double b); - int flag,i,its,j,jj,k,l,nm; - double anorm,c,f,g,h,s,scale,x,y,z,*rv1; - - rv1=dvector(1,n); - g=scale=anorm=0.0; - for (i=1;i<=n;i++) { - l=i+1; - rv1[i]=scale*g; - g=s=scale=0.0; - if (i <= m) { - for (k=i;k<=m;k++) scale += fabs(a[k][i]); - if (scale) { - for (k=i;k<=m;k++) { - a[k][i] /= scale; - s += a[k][i]*a[k][i]; - } - f=a[i][i]; - g = -SIGN(sqrt(s),f); - h=f*g-s; - a[i][i]=f-g; - for (j=l;j<=n;j++) { - for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j]; - f=s/h; - for (k=i;k<=m;k++) a[k][j] += f*a[k][i]; - } - for (k=i;k<=m;k++) a[k][i] *= scale; - } - } - w[i]=scale *g; - g=s=scale=0.0; - if (i <= m && i != n) { - for (k=l;k<=n;k++) scale += fabs(a[i][k]); - if (scale) { - for (k=l;k<=n;k++) { - a[i][k] /= scale; - s += a[i][k]*a[i][k]; - } - f=a[i][l]; - g = -SIGN(sqrt(s),f); - h=f*g-s; - a[i][l]=f-g; - for (k=l;k<=n;k++) rv1[k]=a[i][k]/h; - for (j=l;j<=m;j++) { - for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k]; - for (k=l;k<=n;k++) a[j][k] += s*rv1[k]; - } - for (k=l;k<=n;k++) a[i][k] *= scale; - } - } - anorm=DMAX(anorm,(fabs(w[i])+fabs(rv1[i]))); - } - for (i=n;i>=1;i--) { - if (i < n) { - if (g) { - for (j=l;j<=n;j++) v[j][i]=(a[i][j]/a[i][l])/g; - for (j=l;j<=n;j++) { - for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j]; - for (k=l;k<=n;k++) v[k][j] += s*v[k][i]; - } - } - for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0; - } - v[i][i]=1.0; - g=rv1[i]; - l=i; - } - for (i=IMIN(m,n);i>=1;i--) { - l=i+1; - g=w[i]; - for (j=l;j<=n;j++) a[i][j]=0.0; - if (g) { - g=1.0/g; - for (j=l;j<=n;j++) { - for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j]; - f=(s/a[i][i])*g; - for (k=i;k<=m;k++) a[k][j] += f*a[k][i]; - } - for (j=i;j<=m;j++) a[j][i] *= g; - } else for (j=i;j<=m;j++) a[j][i]=0.0; - ++a[i][i]; - } - for (k=n;k>=1;k--) { - for (its=1;its<=30;its++) { - flag=1; - for (l=k;l>=1;l--) { - nm=l-1; - if ((double)(fabs(rv1[l])+anorm) == anorm) { - flag=0; - break; - } - if ((double)(fabs(w[nm])+anorm) == anorm) break; - } - if (flag) { - c=0.0; - s=1.0; - for (i=l;i<=k;i++) { - f=s*rv1[i]; - rv1[i]=c*rv1[i]; - if ((double)(fabs(f)+anorm) == anorm) break; - g=w[i]; - h=dpythag(f,g); - w[i]=h; - h=1.0/h; - c=g*h; - s = -f*h; - for (j=1;j<=m;j++) { - y=a[j][nm]; - z=a[j][i]; - a[j][nm]=y*c+z*s; - a[j][i]=z*c-y*s; - } - } - } - z=w[k]; - if (l == k) { - if (z < 0.0) { - w[k] = -z; - for (j=1;j<=n;j++) v[j][k] = -v[j][k]; - } - break; - } - if (its == 1000){ - Msg(GERROR, "SVD decomposition: no convergence in 1000 iterations"); - } - x=w[l]; - nm=k-1; - y=w[nm]; - g=rv1[nm]; - h=rv1[k]; - f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y); - g=dpythag(f,1.0); - f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x; - c=s=1.0; - for (j=l;j<=nm;j++) { - i=j+1; - g=rv1[i]; - y=w[i]; - h=s*g; - g=c*g; - z=dpythag(f,h); - rv1[j]=z; - c=f/z; - s=h/z; - f=x*c+g*s; - g = g*c-x*s; - h=y*s; - y *= c; - for (jj=1;jj<=n;jj++) { - x=v[jj][j]; - z=v[jj][i]; - v[jj][j]=x*c+z*s; - v[jj][i]=z*c-x*s; - } - z=dpythag(f,h); - w[j]=z; - if (z) { - z=1.0/z; - c=f*z; - s=h*z; - } - f=c*g+s*y; - x=c*y-s*g; - for (jj=1;jj<=m;jj++) { - y=a[jj][j]; - z=a[jj][i]; - a[jj][j]=y*c+z*s; - a[jj][i]=z*c-y*s; - } - } - rv1[l]=0.0; - rv1[k]=f; - w[k]=x; - } - } - free_dvector(rv1,1,n); -} - -void dsvbksb(double **u, double w[], double **v, int m, int n, double b[], double x[]) -{ - int jj,j,i; - double s,*tmp; - - tmp=dvector(1,n); - for (j=1;j<=n;j++) { - s=0.0; - if (w[j]) { - for (i=1;i<=m;i++) s += u[i][j]*b[i]; - s /= w[j]; - } - tmp[j]=s; - } - for (j=1;j<=n;j++) { - s=0.0; - for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj]; - x[j]=s; - } - free_dvector(tmp,1,n); -} diff --git a/Adapt/fdjac.cpp b/Adapt/fdjac.cpp deleted file mode 100644 index 7a02e764987ec196c4b3f5426701a7e874c56571..0000000000000000000000000000000000000000 --- a/Adapt/fdjac.cpp +++ /dev/null @@ -1,31 +0,0 @@ -// $Id: fdjac.cpp,v 1.4 2001-01-08 08:05:39 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 diff --git a/Adapt/fmin.cpp b/Adapt/fmin.cpp deleted file mode 100644 index a81d482a0829034d223bc56a7702730c087c2e4e..0000000000000000000000000000000000000000 --- a/Adapt/fmin.cpp +++ /dev/null @@ -1,20 +0,0 @@ -// $Id: fmin.cpp,v 1.3 2001-01-08 08:05:40 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 diff --git a/Adapt/lnsrch.cpp b/Adapt/lnsrch.cpp deleted file mode 100644 index e45de2cdf0023706491657e6c9ffe64312815c43..0000000000000000000000000000000000000000 --- a/Adapt/lnsrch.cpp +++ /dev/null @@ -1,80 +0,0 @@ -// $Id: lnsrch.cpp,v 1.4 2001-01-08 08:05:40 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 diff --git a/Adapt/lubksb.cpp b/Adapt/lubksb.cpp deleted file mode 100644 index 65d5dce895e4caf582bb38f675550fc280d7be60..0000000000000000000000000000000000000000 --- a/Adapt/lubksb.cpp +++ /dev/null @@ -1,27 +0,0 @@ -// $Id: lubksb.cpp,v 1.4 2001-01-08 08:05:40 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]; - } -} diff --git a/Adapt/ludcmp.cpp b/Adapt/ludcmp.cpp deleted file mode 100644 index c8a00772375dbc2fe1dee875351e367a8851b694..0000000000000000000000000000000000000000 --- a/Adapt/ludcmp.cpp +++ /dev/null @@ -1,72 +0,0 @@ -// $Id: ludcmp.cpp,v 1.4 2001-01-08 08:05:40 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 diff --git a/Adapt/mnbrak.cpp b/Adapt/mnbrak.cpp deleted file mode 100644 index 5c3edaf06310c8aa45cb06b0b449f24d1049e0d0..0000000000000000000000000000000000000000 --- a/Adapt/mnbrak.cpp +++ /dev/null @@ -1,79 +0,0 @@ -// $Id: mnbrak.cpp,v 1.4 2001-01-08 08:05:40 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 diff --git a/Adapt/newt.cpp b/Adapt/newt.cpp deleted file mode 100644 index f01cb9aed6cc40f82b23287dc7c3a19b63d2a5a3..0000000000000000000000000000000000000000 --- a/Adapt/newt.cpp +++ /dev/null @@ -1,107 +0,0 @@ -// $Id: newt.cpp,v 1.4 2001-01-08 08:05:40 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 diff --git a/Adapt/nrutil.cpp b/Adapt/nrutil.cpp deleted file mode 100644 index 39b174257845459cd298d037d18c55e9b64f874c..0000000000000000000000000000000000000000 --- a/Adapt/nrutil.cpp +++ /dev/null @@ -1,296 +0,0 @@ -// $Id: nrutil.cpp,v 1.5 2001-11-19 13:43:33 geuzaine Exp $ -#include <stdio.h> -#include <stddef.h> -#include <stdlib.h> -#include <malloc.h> - -#include "Gmsh.h" - - -#define NR_END 1 -#define FREE_ARG char* - -void nrerror(char error_text[]) -/* Numerical Recipes standard error handler */ -{ - Msg(GERROR, "%s", error_text); - /* - fprintf(stderr,"Numerical Recipes run-time error...\n"); - fprintf(stderr,"%s\n",error_text); - fprintf(stderr,"...now exiting to system...\n"); - exit(1); - */ -} - -float *vector(long nl, long nh) -/* allocate a float vector with subscript range v[nl..nh] */ -{ - float *v; - - v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float))); - if (!v) nrerror("allocation failure in vector()"); - return v-nl+NR_END; -} - -int *ivector(long nl, long nh) -/* allocate an int vector with subscript range v[nl..nh] */ -{ - int *v; - - v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int))); - if (!v) nrerror("allocation failure in ivector()"); - return v-nl+NR_END; -} - -unsigned char *cvector(long nl, long nh) -/* allocate an unsigned char vector with subscript range v[nl..nh] */ -{ - unsigned char *v; - - v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char))); - if (!v) nrerror("allocation failure in cvector()"); - return v-nl+NR_END; -} - -unsigned long *lvector(long nl, long nh) -/* allocate an unsigned long vector with subscript range v[nl..nh] */ -{ - unsigned long *v; - - v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long))); - if (!v) nrerror("allocation failure in lvector()"); - return v-nl+NR_END; -} - -double *dvector(long nl, long nh) -/* allocate a double vector with subscript range v[nl..nh] */ -{ - double *v; - - v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double))); - if (!v) nrerror("allocation failure in dvector()"); - return v-nl+NR_END; -} - -float **matrix(long nrl, long nrh, long ncl, long nch) -/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */ -{ - long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; - float **m; - - /* allocate pointers to rows */ - m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*))); - if (!m) nrerror("allocation failure 1 in matrix()"); - m += NR_END; - m -= nrl; - - /* allocate rows and set pointers to them */ - m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float))); - if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); - m[nrl] += NR_END; - m[nrl] -= ncl; - - for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; - - /* return pointer to array of pointers to rows */ - return m; -} - -double **dmatrix(long nrl, long nrh, long ncl, long nch) -/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */ -{ - long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; - double **m; - - /* allocate pointers to rows */ - m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*))); - if (!m) nrerror("allocation failure 1 in matrix()"); - m += NR_END; - m -= nrl; - - /* allocate rows and set pointers to them */ - m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double))); - if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); - m[nrl] += NR_END; - m[nrl] -= ncl; - - for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; - - /* return pointer to array of pointers to rows */ - return m; -} - -int **imatrix(long nrl, long nrh, long ncl, long nch) -/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */ -{ - long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; - int **m; - - /* allocate pointers to rows */ - m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*))); - if (!m) nrerror("allocation failure 1 in matrix()"); - m += NR_END; - m -= nrl; - - - /* allocate rows and set pointers to them */ - m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int))); - if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); - m[nrl] += NR_END; - m[nrl] -= ncl; - - for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; - - /* return pointer to array of pointers to rows */ - return m; -} - -float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch, - long newrl, long newcl) -/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */ -{ - long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl; - float **m; - - /* allocate array of pointers to rows */ - m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*))); - if (!m) nrerror("allocation failure in submatrix()"); - m += NR_END; - m -= newrl; - - /* set pointers to rows */ - for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol; - - /* return pointer to array of pointers to rows */ - return m; -} - -float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch) -/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix -declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1 -and ncol=nch-ncl+1. The routine should be called with the address -&a[0][0] as the first argument. */ -{ - long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1; - float **m; - - /* allocate pointers to rows */ - m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*))); - if (!m) nrerror("allocation failure in convert_matrix()"); - m += NR_END; - m -= nrl; - - /* set pointers to rows */ - m[nrl]=a-ncl; - for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol; - /* return pointer to array of pointers to rows */ - return m; -} - -float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh) -/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */ -{ - long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1; - float ***t; - - /* allocate pointers to pointers to rows */ - t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**))); - if (!t) nrerror("allocation failure 1 in f3tensor()"); - t += NR_END; - t -= nrl; - - /* allocate pointers to rows and set pointers to them */ - t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*))); - if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()"); - t[nrl] += NR_END; - t[nrl] -= ncl; - - /* allocate rows and set pointers to them */ - t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float))); - if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()"); - t[nrl][ncl] += NR_END; - t[nrl][ncl] -= ndl; - - for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep; - for(i=nrl+1;i<=nrh;i++) { - t[i]=t[i-1]+ncol; - t[i][ncl]=t[i-1][ncl]+ncol*ndep; - for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep; - } - - /* return pointer to array of pointers to rows */ - return t; -} - -void free_vector(float *v, long nl, long nh) -/* free a float vector allocated with vector() */ -{ - free((FREE_ARG) (v+nl-NR_END)); -} - -void free_ivector(int *v, long nl, long nh) -/* free an int vector allocated with ivector() */ -{ - free((FREE_ARG) (v+nl-NR_END)); -} - -void free_cvector(unsigned char *v, long nl, long nh) -/* free an unsigned char vector allocated with cvector() */ -{ - free((FREE_ARG) (v+nl-NR_END)); -} - -void free_lvector(unsigned long *v, long nl, long nh) -/* free an unsigned long vector allocated with lvector() */ -{ - free((FREE_ARG) (v+nl-NR_END)); -} - -void free_dvector(double *v, long nl, long nh) -/* free a double vector allocated with dvector() */ -{ - free((FREE_ARG) (v+nl-NR_END)); -} - -void free_matrix(float **m, long nrl, long nrh, long ncl, long nch) -/* free a float matrix allocated by matrix() */ -{ - free((FREE_ARG) (m[nrl]+ncl-NR_END)); - free((FREE_ARG) (m+nrl-NR_END)); -} - -void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch) -/* free a double matrix allocated by dmatrix() */ -{ - free((FREE_ARG) (m[nrl]+ncl-NR_END)); - free((FREE_ARG) (m+nrl-NR_END)); -} - -void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch) -/* free an int matrix allocated by imatrix() */ -{ - free((FREE_ARG) (m[nrl]+ncl-NR_END)); - free((FREE_ARG) (m+nrl-NR_END)); -} - -void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch) -/* free a submatrix allocated by submatrix() */ -{ - free((FREE_ARG) (b+nrl-NR_END)); -} - -void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch) -/* free a matrix allocated by convert_matrix() */ -{ - free((FREE_ARG) (b+nrl-NR_END)); -} - -void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch, - long ndl, long ndh) -/* free a float f3tensor allocated by f3tensor() */ -{ - free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END)); - free((FREE_ARG) (t[nrl]+ncl-NR_END)); - free((FREE_ARG) (t+nrl-NR_END)); -} - diff --git a/Adapt/nrutil.h b/Adapt/nrutil.h deleted file mode 100644 index 191fc3528e8f13540eedf4f2cb2c706ec1262d69..0000000000000000000000000000000000000000 --- a/Adapt/nrutil.h +++ /dev/null @@ -1,35 +0,0 @@ -#ifndef _NR_UTILS_H_ -#define _NR_UTILS_H_ - -#include "Numeric.h" - -#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) - -void nrerror (char error_text[]); -float *vector (long nl, long nh); -int *ivector (long nl, long nh); -unsigned char *cvector (long nl, long nh); -unsigned long *lvector (long nl, long nh); -double *dvector (long nl, long nh); -float **matrix (long nrl, long nrh, long ncl, long nch); -double **dmatrix (long nrl, long nrh, long ncl, long nch); -int **imatrix (long nrl, long nrh, long ncl, long nch); -float **submatrix (float **a, long oldrl, long oldrh, long oldcl, long oldch, - long newrl, long newcl); -float **convert_matrix (float *a, long nrl, long nrh, long ncl, long nch); -float ***f3tensor (long nrl, long nrh, long ncl, long nch, long ndl, long ndh); -void free_vector (float *v, long nl, long nh); -void free_ivector (int *v, long nl, long nh); -void free_cvector (unsigned char *v, long nl, long nh); -void free_lvector (unsigned long *v, long nl, long nh); -void free_dvector (double *v, long nl, long nh); -void free_matrix (float **m, long nrl, long nrh, long ncl, long nch); -void free_dmatrix (double **m, long nrl, long nrh, long ncl, long nch); -void free_imatrix (int **m, long nrl, long nrh, long ncl, long nch); -void free_submatrix (float **b, long nrl, long nrh, long ncl, long nch); -void free_convert_matrix (float **b, long nrl, long nrh, long ncl, long nch); -void free_f3tensor (float ***t, long nrl, long nrh, long ncl, long nch, - long ndl, long ndh); - - -#endif /* _NR_UTILS_H_ */