Remove Adapt directory

Adapt/Adapt.cpp

-// $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) ;
-}

Adapt/Adapt.h

-#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

Adapt/Makefile

-# $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

Adapt/brent.cpp

-// $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. */

Adapt/dsvdcmp.cpp

-// $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);
-}

Adapt/fdjac.cpp

-// $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

Adapt/fmin.cpp

-// $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

Adapt/lnsrch.cpp

-// $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

Adapt/lubksb.cpp

-// $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];
-    }
-}

Adapt/ludcmp.cpp

-// $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

Adapt/mnbrak.cpp

-// $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

Adapt/newt.cpp

-// $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

Adapt/nrutil.cpp

-// $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));
-}
-

Adapt/nrutil.h

-#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_ */