From 7b03ad3e96658e29ecd214587280f7f4270c58c1 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Mon, 22 Dec 2008 22:34:07 +0000
Subject: [PATCH] More "numeric" work: removed all traces of Numerical Recipes,
as the NR license is incompatible with our use. Replaced with Lapack for
LU/SVD. Still need to recode newt & brent.
---
Geo/GFace.cpp | 81 +--
Geo/Geo.cpp | 9 +-
Numeric/GmshMatrix.h | 233 ++++----
...gmsh_predicates.cpp => GmshPredicates.cpp} | 0
Numeric/Makefile | 17 +-
Numeric/Numeric.cpp | 545 +++++++++++++++---
Numeric/Numeric.h | 35 +-
Numeric/NumericEmbedded.cpp | 38 ++
Numeric/NumericEmbedded.h | 1 +
Numeric/gsl_brent.cpp | 132 -----
Numeric/gsl_min.cpp | 265 ---------
Numeric/gsl_newt.cpp | 108 ----
Post/Makefile | 2 +-
configure | 56 +-
configure.in | 38 +-
contrib/NR/Makefile | 61 --
contrib/NR/brent.cpp | 77 ---
contrib/NR/dpythag.cpp | 14 -
contrib/NR/dsvdcmp.cpp | 183 ------
contrib/NR/fdjac.cpp | 29 -
contrib/NR/fmin.cpp | 20 -
contrib/NR/lnsrch.cpp | 65 ---
contrib/NR/lubksb.cpp | 23 -
contrib/NR/ludcmp.cpp | 60 --
contrib/NR/mnbrak.cpp | 67 ---
contrib/NR/newt.cpp | 92 ---
contrib/NR/nrutil.cpp | 293 ----------
contrib/NR/nrutil.h | 43 --
28 files changed, 680 insertions(+), 1907 deletions(-)
rename Numeric/{gmsh_predicates.cpp => GmshPredicates.cpp} (100%)
delete mode 100644 Numeric/gsl_brent.cpp
delete mode 100644 Numeric/gsl_min.cpp
delete mode 100644 Numeric/gsl_newt.cpp
delete mode 100644 contrib/NR/Makefile
delete mode 100644 contrib/NR/brent.cpp
delete mode 100644 contrib/NR/dpythag.cpp
delete mode 100644 contrib/NR/dsvdcmp.cpp
delete mode 100644 contrib/NR/fdjac.cpp
delete mode 100644 contrib/NR/fmin.cpp
delete mode 100644 contrib/NR/lnsrch.cpp
delete mode 100644 contrib/NR/lubksb.cpp
delete mode 100644 contrib/NR/ludcmp.cpp
delete mode 100644 contrib/NR/mnbrak.cpp
delete mode 100644 contrib/NR/newt.cpp
delete mode 100644 contrib/NR/nrutil.cpp
delete mode 100644 contrib/NR/nrutil.h
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index e001d89afe..85b787e0d5 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -10,6 +10,7 @@
#include "MElement.h"
#include "GmshMessage.h"
#include "VertexArray.h"
+#include "GmshMatrix.h"
#if defined(HAVE_GMSH_EMBEDDED)
#include "GmshEmbedded.h"
@@ -17,14 +18,6 @@
#include "Numeric.h"
#include "GaussLegendre1D.h"
#include "Context.h"
-#if defined(HAVE_GSL)
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_linalg.h>
-#else
-#define NRANSI
-#include "nrutil.h"
-void dsvdcmp(double **a, int m, int n, double w[], double **v);
-#endif
#endif
extern Context_T CTX;
@@ -241,7 +234,6 @@ void GFace::computeMeanPlane(const std::vector<MVertex*> &points)
void GFace::computeMeanPlane(const std::vector<SPoint3> &points)
{
-#if !defined(HAVE_GMSH_EMBEDDED)
// The concept of a mean plane computed in the sense of least
// squares is fine for plane surfaces(!), but not really the best
// one for non-plane surfaces. Indeed, imagine a quarter of a circle
@@ -263,63 +255,29 @@ void GFace::computeMeanPlane(const std::vector<SPoint3> &points)
ym /= (double)ndata;
zm /= (double)ndata;
- int min;
- double res[4], svd[3];
-#if defined(HAVE_GSL)
- gsl_matrix *U = gsl_matrix_alloc(ndata, na);
- gsl_matrix *V = gsl_matrix_alloc(na, na);
- gsl_vector *W = gsl_vector_alloc(na);
- gsl_vector *TMPVEC = gsl_vector_alloc(na);
+ Double_Matrix U(ndata, na), V(na, na);
+ Double_Vector sigma(na);
for(int i = 0; i < ndata; i++) {
- gsl_matrix_set(U, i, 0, points[i].x() - xm);
- gsl_matrix_set(U, i, 1, points[i].y() - ym);
- gsl_matrix_set(U, i, 2, points[i].z() - zm);
+ U(i, 0) = points[i].x() - xm;
+ U(i, 1) = points[i].y() - ym;
+ U(i, 2) = points[i].z() - zm;
}
- gsl_linalg_SV_decomp(U, V, W, TMPVEC);
- svd[0] = gsl_vector_get(W, 0);
- svd[1] = gsl_vector_get(W, 1);
- svd[2] = gsl_vector_get(W, 2);
+ U.svd(V, sigma);
+ double res[4], svd[3];
+ svd[0] = sigma(0);
+ svd[1] = sigma(1);
+ svd[2] = sigma(2);
+ int min;
if(fabs(svd[0]) < fabs(svd[1]) && fabs(svd[0]) < fabs(svd[2]))
min = 0;
else if(fabs(svd[1]) < fabs(svd[0]) && fabs(svd[1]) < fabs(svd[2]))
min = 1;
else
min = 2;
- res[0] = gsl_matrix_get(V, 0, min);
- res[1] = gsl_matrix_get(V, 1, min);
- res[2] = gsl_matrix_get(V, 2, min);
+ res[0] = V(0, min);
+ res[1] = V(1, min);
+ res[2] = V(2, min);
norme(res);
- gsl_matrix_free(U);
- gsl_matrix_free(V);
- gsl_vector_free(W);
- gsl_vector_free(TMPVEC);
-#else
- double **U = dmatrix(1, ndata, 1, na);
- double **V = dmatrix(1, na, 1, na);
- double *W = dvector(1, na);
- for(int i = 0; i < ndata; i++) {
- U[i + 1][1] = points[i].x() - xm;
- U[i + 1][2] = points[i].y() - ym;
- U[i + 1][3] = points[i].z() - zm;
- }
- dsvdcmp(U, ndata, na, W, V);
- if(fabs(W[1]) < fabs(W[2]) && fabs(W[1]) < fabs(W[3]))
- min = 1;
- else if(fabs(W[2]) < fabs(W[1]) && fabs(W[2]) < fabs(W[3]))
- min = 2;
- else
- min = 3;
- svd[0] = W[1];
- svd[1] = W[2];
- svd[2] = W[3];
- res[0] = V[1][min];
- res[1] = V[2][min];
- res[2] = V[3][min];
- norme(res);
- free_dmatrix(U, 1, ndata, 1, na);
- free_dmatrix(V, 1, na, 1, na);
- free_dvector(W, 1, na);
-#endif
double ex[3], t1[3], t2[3];
@@ -428,7 +386,6 @@ end:
}
}
}
-#endif
}
void GFace::getMeanPlaneData(double VX[3], double VY[3],
@@ -481,17 +438,11 @@ double GFace::curvature(const SPoint2 ¶m) const
SVector3 dndu = 500 * (n2 - n1);
SVector3 dndv = 500 * (n4 - n3);
- // double c = fabs(dot(dndu, du) + dot(dndv, dv)) / detJ;
-
-
double ddu = dot(dndu,du);
double ddv = dot(dndv,dv);
double c = std::max(fabs(ddu),fabs(ddv))/detJ;
-
-
// Msg::Info("c = %g detJ %g", c, detJ);
-
return c;
}
@@ -505,7 +456,6 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
double &U, double &V, const double relax,
const bool onSurface) const
{
-#if !defined(HAVE_GMSH_EMBEDDED)
const double Precision = 1.e-8;
const int MaxIter = 25;
const int NumInitGuess = 11;
@@ -584,7 +534,6 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
Msg::Info("point %g %g %g : Relaxation factor = %g", X, Y, Z, 0.75 * relax);
XYZtoUV(X, Y, Z, U, V, 0.75 * relax);
}
-#endif
}
SPoint2 GFace::parFromPoint(const SPoint3 &p) const
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 302b465c00..fe4353a52d 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -2886,9 +2886,9 @@ static Vertex *VERTEX;
static double min1d(double (*funct) (double), double *xmin)
{
- // 0. for tolerance allows for maximum as code in gsl_brent
- return (brent(CURVE->ubeg, 0.5*(CURVE->ubeg + CURVE->uend), CURVE->uend,
- funct, 0., xmin));
+ // tolerance==0. allows for maximum as code in gsl_brent
+ return brent(CURVE->ubeg, 0.5 * (CURVE->ubeg + CURVE->uend), CURVE->uend,
+ funct, 0., xmin);
}
static void projectPS(int N, double x[], double res[])
@@ -2943,8 +2943,8 @@ bool ProjectPointOnSurface(Surface *s, Vertex &p, double u[2])
int check;
SURFACE = s;
VERTEX = &p;
-
newt(x, 2, &check, projectPS);
+
Vertex vv = InterpolateSurface(s, x[1], x[2], 0, 0);
double res[3];
projectPS(2, x, res);
@@ -3091,6 +3091,7 @@ static bool IntersectCurveSurface(Curve *c, Surface *s, double x[4])
SURFACE = s;
CURVE = c;
newt(x, 3, &check, intersectCS);
+
if(check) return false;
return true;
}
diff --git a/Numeric/GmshMatrix.h b/Numeric/GmshMatrix.h
index 0df4528384..dc72f50663 100644
--- a/Numeric/GmshMatrix.h
+++ b/Numeric/GmshMatrix.h
@@ -7,8 +7,34 @@
#define _GMSH_MATRIX_H_
#include <math.h>
+#include <algorithm>
#include "GmshMessage.h"
+#if defined(HAVE_BLAS)
+extern "C" {
+ void dgemm_(const char *transa, const char *transb,
+ int *l, int *n, int *m, double *alpha,
+ const void *a, int *lda, void *b, int *ldb,
+ double *beta, void *c, int *ldc);
+ void dgemv_(const char *trans, int *m, int *n, double *alpha,
+ void *a, int *lda, void *x, int *incx,
+ double *beta, void *y, int *incy);
+}
+#endif
+
+#if defined(HAVE_LAPACK)
+extern "C" {
+ void dgesv_(const int *N, const int *nrhs, double *A, const int *lda,
+ int *ipiv, double *b, const int *ldb, int *info);
+ void dgetrf_(const int *M, const int *N, double *A, const int *lda,
+ int *ipiv, int *info);
+ void dgesvd_(const char* jobu, const char *jobvt, const int *M, const int *N,
+ double *A, const int *lda, double *S, double* U, const int *ldu,
+ double *VT, const int *ldvt, double *work, const int *lwork,
+ const int *info);
+}
+#endif
+
template <class SCALAR> class Gmsh_Matrix;
template <class SCALAR>
@@ -60,79 +86,6 @@ class Gmsh_Matrix
private:
int _r, _c;
SCALAR *_data;
- void _back_substitution(int *indx, SCALAR *b)
- {
- int i, ii = -1, ip, j;
- SCALAR sum;
- for(i = 0; i < _c; i++){
- ip = indx[i];
- sum = b[ip];
- b[ip] = b[i];
- if(ii != -1)
- for(j = ii; j <= i - 1; j++) sum -= (*this)(i, j) * b[j];
- else if(sum) ii = i;
- b[i] = sum;
- }
- for(i = _c - 1; i >= 0; i--){
- sum = b[i];
- for(j = i + 1; j < _c; j++) sum -= (*this)(i, j) * b[j];
- b[i] = sum / (*this)(i, i);
- }
- }
- bool _lu_decomposition(int *indx, SCALAR &determinant)
- {
- if(_r != _c)
- Msg::Fatal("Cannot LU factorize non-square matrix");
- int i, imax, j, k;
- SCALAR big, dum, sum, temp;
- SCALAR *vv = new SCALAR[_c];
- determinant = 1.;
- for(i = 0; i < _c; i++){
- big = 0.;
- for(j = 0; j < _c; j++)
- if((temp = fabs((*this)(i, j))) > big) big = temp;
- if(big == 0.) {
- delete [] vv;
- return false;
- }
- vv[i] = 1. / big;
- }
- for(j = 0; j < _c; j++){
- for(i = 0; i < j; i++){
- sum = (*this)(i, j);
- for(k = 0; k < i; k++) sum -= (*this)(i, k)*(*this)(k, j);
- (*this)(i, j) = sum;
- }
- big = 0.;
- for(i = j; i < _c; i++){
- sum = (*this)(i, j);
- for(k = 0; k < j; k++)
- sum -= (*this)(i, k) * (*this)(k, j);
- (*this)(i, j) = sum;
- if((dum = vv[i] * fabs(sum)) >= big){
- big = dum;
- imax = i;
- }
- }
- if(j != imax){
- for(k = 0; k < _c; k++){
- dum = (*this)(imax, k);
- (*this)(imax, k) = (*this)(j, k);
- (*this)(j, k) = dum;
- }
- determinant = -(determinant);
- vv[imax] = vv[j];
- }
- indx[j] = imax;
- if((*this)(j, j) == 0.) (*this)(j, j) = 1.e-20;
- if(j != _c){
- dum = 1. / ((*this)(j, j));
- for(i = j + 1; i < _c; i++) (*this)(i, j) *= dum;
- }
- }
- delete [] vv;
- return true;
- }
public:
Gmsh_Matrix(int r, int c) : _r(r), _c(c)
{
@@ -179,20 +132,34 @@ class Gmsh_Matrix
}
inline void mult(const Gmsh_Matrix<SCALAR> &b, Gmsh_Matrix<SCALAR> &c)
{
+#if 0 // defined(HAVE_BLAS)
+ void dgemm_(const char *transa, const char *transb,
+ int *l, int *n, int *m, double *alpha,
+ const void *a, int *lda, void *b, int *ldb,
+ double *beta, void *c, int *ldc);
+#else
c.scale(0.);
for(int i = 0; i < _r; i++)
for(int j = 0; j < b.size2(); j++)
for(int k = 0; k < _c; k++)
c._data[i + _r * j] += (*this)(i, k) * b(k, j);
+#endif
}
inline void blas_dgemm(const Gmsh_Matrix<SCALAR> &b, Gmsh_Matrix<SCALAR> &c,
const SCALAR alpha=1., const SCALAR beta=1.)
{
+#if 0 // defined(HAVE_BLAS)
+ void dgemm_(const char *transa, const char *transb,
+ int *l, int *n, int *m, double *alpha,
+ const void *a, int *lda, void *b, int *ldb,
+ double *beta, void *c, int *ldc);
+#else
Gmsh_Matrix<SCALAR> temp(b.size1(), c.size2());
temp.mult(b, c);
scale(beta);
temp.scale(alpha);
add(temp);
+#endif
}
inline void set_all(const SCALAR &m)
{
@@ -217,45 +184,39 @@ class Gmsh_Matrix
}
inline void mult(const Gmsh_Vector<SCALAR> &x, Gmsh_Vector<SCALAR> &y)
{
+#if 0 //defined(HAVE_BLAS)
+ void dgemv_(const char *trans, int *m, int *n, double *alpha,
+ void *a, int *lda, void *x, int *incx,
+ double *beta, void *y, int *incy);
+#else
y.scale(0.);
for(int i = 0; i < _r; i++)
for(int j = 0; j < _c; j++)
y._data[i] += (*this)(i, j) * x(j);
+#endif
}
- inline bool lu_solve(const Gmsh_Vector<SCALAR> &rhs, Gmsh_Vector<SCALAR> &result)
+ inline Gmsh_Matrix<SCALAR> transpose()
{
- int *indx = new int[_c];
- SCALAR d;
- if(!_lu_decomposition(indx, d)){
- delete [] indx;
- return false;
- }
- for(int i = 0; i < _c; i++) result(i) = rhs(i);
- _back_substitution(indx, result._data);
- delete [] indx;
- return true;
+ Gmsh_Matrix<SCALAR> T(size2(), size1());
+ for(int i = 0; i < size1(); i++)
+ for(int j = 0; j < size2(); j++)
+ T(j, i) = (*this)(i, j);
+ return T;
}
- inline bool invert()
+ inline bool lu_solve(const Gmsh_Vector<SCALAR> &rhs, Gmsh_Vector<SCALAR> &result)
{
- Gmsh_Matrix y(_r, _c);
- SCALAR *col = new SCALAR[_c];
- int *indx = new int[_c];
- SCALAR d;
- if (!_lu_decomposition(indx, d)){
- delete [] col;
- delete [] indx;
- return false;
- }
- for(int j = 0; j < _c; j++){
- for(int i = 0; i < _c; i++) col[i] = 0.0;
- col[j] = 1.0;
- _back_substitution(indx, col);
- for(int i = 0; i < _c; i++) y(i, j) = col[i];
- }
- (*this) = y;
- delete [] col;
- delete [] indx;
- return true;
+#if defined(HAVE_LAPACK)
+ int N = size1(), nrhs = 1, lda = N, ldb = N, info;
+ int *ipiv = new int[N];
+ for(int i = 0; i < N; i++) result(i) = rhs(i);
+ dgesv_(&N, &nrhs, _data, &lda, ipiv, result._data, &ldb, &info);
+ delete [] ipiv;
+ if(info == 0) return true;
+ Msg::Error("Problem in LAPACK LU (info=%d)", info);
+#else
+ Msg::Error("LU factorization requires LAPACK");
+#endif
+ return false;
}
Gmsh_Matrix<SCALAR> cofactor(int i, int j) const
{
@@ -272,14 +233,38 @@ class Gmsh_Matrix
}
SCALAR determinant() const
{
- Gmsh_Matrix<SCALAR> tmp = *this;
- SCALAR factor = 1.;
- int *indx = new int[_c];
- if(!tmp._lu_decomposition(indx, factor)) return 0.;
- SCALAR det = factor;
- for(int i = 0; i < _c; i++) det *= tmp(i, i);
- delete [] indx;
+ Gmsh_Matrix<SCALAR> tmp(*this);
+#if defined(HAVE_LAPACK)
+ int M = size1(), N = size2(), lda = size1(), info;
+ int *ipiv = new int[std::min(M, N)];
+ dgetrf_(&M, &N, tmp._data, &lda, ipiv, &info);
+ SCALAR det = 1.;
+ for(int i = 0; i < size1(); i++){
+ det *= tmp(i, i);
+ if(ipiv[i] != i + 1) det = -det;
+ }
return det;
+#else
+ Msg::Error("Determinant computation requires LAPACK");
+ return 0.;
+#endif
+ }
+ bool svd(Gmsh_Matrix<SCALAR> &V, Gmsh_Vector<SCALAR> &S)
+ {
+#if defined(HAVE_LAPACK)
+ Gmsh_Matrix<SCALAR> VT(V.size2(), V.size1());
+ int M = size1(), N = size2(), LDA = size1(), LDVT = VT.size1(), info;
+ int LWORK = std::max(3 * std::min(M, N) + std::max(M, N), 5 * std::min(M, N));
+ Gmsh_Vector<SCALAR> WORK(LWORK);
+ dgesvd_("O", "A", &M, &N, _data, &LDA, S._data, _data, &LDA,
+ VT._data, &LDVT, WORK._data, &LWORK, &info);
+ V = VT.transpose();
+ if(info == 0) return true;
+ Msg::Error("Problem in LAPACK SVD (info=%d)", info);
+#else
+ Msg::Error("Singular value decomposition requires LAPACK");
+#endif
+ return false;
}
};
@@ -403,24 +388,20 @@ class GSL_Matrix
{
gsl_blas_dgemv(CblasNoTrans, 1., _data, x._data, 0., b._data);
}
- inline bool lu_solve(const GSL_Vector &rhs, GSL_Vector &result)
+ inline GSL_Matrix transpose()
{
- int s;
- gsl_permutation * p = gsl_permutation_alloc(size1());
- gsl_linalg_LU_decomp(_data, p, &s);
- gsl_linalg_LU_solve(_data, p, rhs._data, result._data);
- gsl_permutation_free(p);
- return true;
+ GSL_Matrix T(size2(), size1());
+ for(int i = 0; i < size1(); i++)
+ for(int j = 0; j < size2(); j++)
+ T(j, i) = (*this)(i, j);
+ return T;
}
- inline bool invert()
+ inline bool lu_solve(const GSL_Vector &rhs, GSL_Vector &result)
{
int s;
gsl_permutation *p = gsl_permutation_alloc(size1());
gsl_linalg_LU_decomp(_data, p, &s);
- gsl_matrix *inv = gsl_matrix_calloc(size1(), size2());
- gsl_linalg_LU_invert(_data, p, inv) ;
- gsl_matrix_memcpy(_data, inv);
- gsl_matrix_free(inv);
+ gsl_linalg_LU_solve(_data, p, rhs._data, result._data);
gsl_permutation_free(p);
return true;
}
@@ -446,6 +427,12 @@ class GSL_Matrix
gsl_permutation_free(p);
return gsl_linalg_LU_det(tmp._data, s);
}
+ bool svd(GSL_Matrix &V, GSL_Vector &S)
+ {
+ GSL_Vector tmp(S.size());
+ gsl_linalg_SV_decomp(_data, V._data, S._data, tmp._data);
+ return false;
+ }
};
typedef GSL_Matrix Double_Matrix;
diff --git a/Numeric/gmsh_predicates.cpp b/Numeric/GmshPredicates.cpp
similarity index 100%
rename from Numeric/gmsh_predicates.cpp
rename to Numeric/GmshPredicates.cpp
diff --git a/Numeric/Makefile b/Numeric/Makefile
index f07d08cb20..b0eea45326 100644
--- a/Numeric/Makefile
+++ b/Numeric/Makefile
@@ -24,10 +24,7 @@ SRC = Numeric.cpp\
GaussQuadratureTet.cpp\
GaussQuadratureHex.cpp\
GaussLegendreSimplex.cpp\
- gmsh_predicates.cpp\
- gsl_newt.cpp\
- gsl_min.cpp\
- gsl_brent.cpp
+ GmshPredicates.cpp\
OBJ = ${SRC:.cpp=${OBJEXT}}
@@ -56,8 +53,8 @@ depend:
# DO NOT DELETE THIS LINE
Numeric${OBJEXT}: Numeric.cpp ../Common/GmshMessage.h Numeric.h \
- NumericEmbedded.h
-NumericEmbedded${OBJEXT}: NumericEmbedded.cpp NumericEmbedded.h \
+ NumericEmbedded.h GmshMatrix.h
+NumericEmbedded${OBJEXT}: NumericEmbedded.cpp NumericEmbedded.h GmshMatrix.h \
../Common/GmshMessage.h
gmshAssembler${OBJEXT}: gmshAssembler.cpp ../Geo/MVertex.h ../Geo/SPoint2.h \
../Geo/SPoint3.h gmshAssembler.h gmshLinearSystem.h
@@ -122,10 +119,4 @@ GaussQuadratureHex${OBJEXT}: GaussQuadratureHex.cpp Gauss.h Numeric.h \
NumericEmbedded.h GaussLegendre1D.h
GaussLegendreSimplex${OBJEXT}: GaussLegendreSimplex.cpp Gauss.h Numeric.h \
NumericEmbedded.h GaussLegendre1D.h
-gmsh_predicates${OBJEXT}: gmsh_predicates.cpp
-gsl_newt${OBJEXT}: gsl_newt.cpp ../Common/GmshMessage.h Numeric.h \
- NumericEmbedded.h
-gsl_min${OBJEXT}: gsl_min.cpp ../Common/GmshMessage.h Numeric.h \
- NumericEmbedded.h
-gsl_brent${OBJEXT}: gsl_brent.cpp ../Common/GmshMessage.h Numeric.h \
- NumericEmbedded.h
+GmshPredicates${OBJEXT}: GmshPredicates.cpp
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 4c7d0233a2..35b4256978 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -5,16 +5,19 @@
#include "GmshMessage.h"
#include "Numeric.h"
+#include "GmshMatrix.h"
-// Check GSL version. We need at least 1.2, since all versions <=
-// 1.1.1 have a buggy SVD routine (infinite loop fixed on Sun Jun 16
-// 11:45:29 2002 in GSL's cvs tree). We check this at run time since
-// Gmsh could be distributed with the gsl dynamically linked.
+// This file contains the routines that depend on the GSL, some of
+// which need to be reimplemented (e.g. in terms of Gmsh_Matrix)
+// before we can get rid of the GSL completely.
#if defined(HAVE_GSL)
#include <gsl/gsl_errno.h>
-#include <gsl/gsl_vector.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_min.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_multiroots.h>
#include <gsl/gsl_linalg.h>
void my_gsl_msg(const char *reason, const char *file, int line,
@@ -33,90 +36,482 @@ int check_gsl()
return 1;
}
-#else
-
-#define NRANSI
-#include "nrutil.h"
-void dsvdcmp(double **a, int m, int n, double w[], double **v);
+static double (*nrfuncbrent) (double);
-int check_gsl()
+static double fn1(double x, void *params)
{
- // initilize robust geometric predicates
- gmsh::exactinit() ;
- return 1;
+ double val = nrfuncbrent(x);
+ return val;
}
-#endif
+static int gsl_min_fminimizer_set_with_values_MOD(gsl_min_fminimizer *s,
+ gsl_function *f,
+ double x_minimum, double f_minimum,
+ double x_lower, double f_lower,
+ double x_upper, double f_upper)
+{
+ s->function = f;
+ s->x_minimum = x_minimum;
+ s->x_lower = x_lower;
+ s->x_upper = x_upper;
+ if (x_lower > x_upper){
+ GSL_ERROR ("invalid interval (lower > upper)", GSL_EINVAL);
+ }
+ s->f_lower = f_lower;
+ s->f_upper = f_upper;
+ s->f_minimum = f_minimum;
+ return (s->type->set) (s->state, s->function,
+ x_minimum, f_minimum,
+ x_lower, f_lower,
+ x_upper, f_upper);
+}
-void invert_singular_matrix3x3(double MM[3][3], double II[3][3])
+// Returns the minimum betwen ax and cx to a given tolerance tol using
+// brent's method.
+double brent(double ax, double bx, double cx,
+ double (*f) (double), double tol, double *xmin)
{
- int i, j, k, n = 3;
- double TT[3][3];
+#define MAXITER 100
+ // The solver can stall at the following internal tolerance in brent - so
+ // this is about the lowest possible tolerance.
+ if(tol == 0.) tol = 10. * GSL_SQRT_DBL_EPSILON;
+ const double tolsq = tol*tol;
+ int status;
+ int iter = 0;
+ double a, b, c; // a < b < c
+ double fa, fb, fc;
+ const gsl_min_fminimizer_type *T;
+ gsl_min_fminimizer *s;
+ gsl_function F;
- for(i = 1; i <= n; i++) {
- for(j = 1; j <= n; j++) {
- II[i - 1][j - 1] = 0.0;
- TT[i - 1][j - 1] = 0.0;
- }
+ // gsl wants a<b
+ b = bx;
+ if(ax < cx) {
+ a = ax;
+ c = cx;
+ }
+ else {
+ a = cx;
+ c = ax;
}
-#if defined(HAVE_GSL)
- gsl_matrix *M = gsl_matrix_alloc(3, 3);
- gsl_matrix *V = gsl_matrix_alloc(3, 3);
- gsl_vector *W = gsl_vector_alloc(3);
- gsl_vector *TMPVEC = gsl_vector_alloc(3);
- for(i = 1; i <= n; i++) {
- for(j = 1; j <= n; j++) {
- gsl_matrix_set(M, i - 1, j - 1, MM[i - 1][j - 1]);
+ nrfuncbrent = f;
+ F.function = &fn1;
+ F.params = 0;
+
+ fa = f(a);
+ fb = f(b);
+ fc = f(c);
+ if(fb > fa || fb > fc) {
+ if(a < c) {
+ b = a;
+ fb = fa;
}
- }
- gsl_linalg_SV_decomp(M, V, W, TMPVEC);
- for(i = 1; i <= n; i++) {
- for(j = 1; j <= n; j++) {
- double ww = gsl_vector_get(W, i - 1);
- if(fabs(ww) > 1.e-16) { //singular value precision
- TT[i - 1][j - 1] += gsl_matrix_get(M, j - 1, i - 1) / ww;
- }
+ else {
+ b = c;
+ fb = fc;
}
}
- for(i = 1; i <= n; i++) {
- for(j = 1; j <= n; j++) {
- for(k = 1; k <= n; k++) {
- II[i - 1][j - 1] +=
- gsl_matrix_get(V, i - 1, k - 1) * TT[k - 1][j - 1];
- }
- }
+
+ // if a-c < tol, return func(b)
+ if(gsl_min_test_interval(a, c, tolsq, tol) == GSL_SUCCESS) {
+ *xmin = b;
+ return (fb);
}
- gsl_matrix_free(M);
- gsl_matrix_free(V);
- gsl_vector_free(W);
- gsl_vector_free(TMPVEC);
-#else
- double **M = dmatrix(1, 3, 1, 3);
- double **V = dmatrix(1, 3, 1, 3);
- double *W = dvector(1, 3);
- for(i = 1; i <= n; i++) {
- for(j = 1; j <= n; j++) {
- M[i][j] = MM[i - 1][j - 1];
- }
+
+ T = gsl_min_fminimizer_brent;
+ s = gsl_min_fminimizer_alloc(T);
+ // gsl_min_fminimizer_set(s, &F, b, a, c); // Standard
+ // This mod avoids some error checks so we can have a minimum at an extent.
+ gsl_min_fminimizer_set_with_values_MOD(s, &F, b, fb, a, fa, c, fc);
+
+ do {
+ iter++;
+ status = gsl_min_fminimizer_iterate(s);
+ if(status)
+ break; // solver problem
+ a = gsl_min_fminimizer_x_lower(s);
+ c = gsl_min_fminimizer_x_upper(s);
+ // fb = gsl_min_fminimizer_f_minimum(s);
+ // fa = gsl_min_fminimizer_f_lower(s);
+ // fc = gsl_min_fminimizer_f_upper(s);
+ status = gsl_min_test_interval(a, c, tolsq, tol);
+ // ... in case curve becomes too flat - not required?
+ // std::fabs(fa - fc) < (std::fabs(fb) + tol)*tol
}
- dsvdcmp(M, n, n, W, V);
- for(i = 1; i <= n; i++) {
- for(j = 1; j <= n; j++) {
- if(fabs(W[i]) > 1.e-16) { //singular value precision
- TT[i - 1][j - 1] += M[j][i] / W[i];
- }
- }
+ while(status == GSL_CONTINUE && iter < MAXITER);
+ b = gsl_min_fminimizer_x_minimum(s);
+
+ if(status != GSL_SUCCESS)
+ Msg::Error("minimization did not converge: f(%g) = %g", b, fn1(b, NULL));
+
+ *xmin = b;
+ gsl_min_fminimizer_free(s);
+ return fn1(b, NULL);
+}
+
+
+#define MAX_DIM_NEWT 10
+#define MAXITER 100
+#define PREC 1.e-8
+
+static int nrdim;
+static double nru[MAX_DIM_NEWT], nrv[MAX_DIM_NEWT];
+static void (*nrfunc) (int n, double x[], double y[]);
+struct gsl_dummy{ int i; };
+
+static void convert_vector_from_gsl(const gsl_vector * gv, double *v)
+{
+ int i, m;
+ m = gv->size;
+ for(i = 0; i < m; i++) {
+ v[i + 1] = gsl_vector_get(gv, i);
}
- for(i = 1; i <= n; i++) {
- for(j = 1; j <= n; j++) {
- for(k = 1; k <= n; k++) {
- II[i - 1][j - 1] += V[i][k] * TT[k - 1][j - 1];
- }
- }
+}
+
+static void convert_vector_to_gsl(double *v, int n, gsl_vector * gv)
+{
+ int i;
+ for(i = 0; i < n; i++) {
+ gsl_vector_set(gv, i, v[i + 1]);
}
- free_dmatrix(M, 1, n, 1, n);
- free_dmatrix(V, 1, n, 1, n);
- free_dvector(W, 1, n);
-#endif
}
+
+static int gslfunc(const gsl_vector * xx, void *params, gsl_vector * f)
+{
+ convert_vector_from_gsl(xx, nru);
+ (*nrfunc) (nrdim, nru, nrv);
+ // Msg::Info("f(%lf,%lf) = %lf %lf\n",nru[1],nru[2],nrv[1],nrv[2]);
+ convert_vector_to_gsl(nrv, nrdim, f);
+ return GSL_SUCCESS;
+}
+
+// Warning: for compatibility with the old newt from NR, x[] and the
+// arguments of func() are indexed from 1 to n...
+void newt(double x[], int n, int *check,
+ void (*func) (int, double[], double[]))
+{
+ const gsl_multiroot_fsolver_type *T;
+ gsl_multiroot_fsolver *s;
+ int status;
+ size_t iter = 0;
+ struct gsl_dummy p = { 1 };
+ gsl_multiroot_function f = { &gslfunc, n, &p };
+ gsl_vector *xx = gsl_vector_alloc(n);
+
+ if(n > MAX_DIM_NEWT - 1)
+ Msg::Fatal("Maximum Newton dimension exceeded\n");
+ nrdim = n;
+
+ nrfunc = func;
+ convert_vector_to_gsl(x, n, xx);
+
+ T = gsl_multiroot_fsolver_hybrid;
+ s = gsl_multiroot_fsolver_alloc(T, n);
+ gsl_multiroot_fsolver_set(s, &f, xx);
+
+ do {
+ iter++;
+ status = gsl_multiroot_fsolver_iterate(s);
+ // Msg::Info("status %d %d %d %lf %lf\n",
+ // status,n,iter,gsl_vector_get(s->x,0),gsl_vector_get(s->x,1));
+ if(status)
+ break; // solver problem
+ status = gsl_multiroot_test_residual(s->f, n * PREC);
+ }
+ while(status == GSL_CONTINUE && iter < MAXITER);
+
+ if(status == GSL_CONTINUE) {
+ *check = 1; // problem !!!
+ }
+ else {
+ // Msg::Info("status %d %d %d %lf %lf\n",
+ // status,n,iter,gsl_vector_get(s->x,0),gsl_vector_get(s->x,1));
+ convert_vector_from_gsl(s->x, x);
+ *check = 0; // converged
+ }
+
+ gsl_multiroot_fsolver_free(s);
+ gsl_vector_free(xx);
+}
+
+static double (*f_stat) (double, double, void *data);
+static void (*df_stat) (double, double, double &, double &, double &, void *data);
+static double (*f_stat3) (double, double, double, void *data);
+static void (*df_stat3) (double, double, double, double &, double &,
+ double &, double &,void *data);
+static double (*f_statN) (double*, void *data);
+static void (*df_statN) (double*, double *, double &,void *data);
+
+static double fobj(const gsl_vector * x, void * data)
+{
+ double u, v;
+ u = gsl_vector_get(x, 0);
+ v = gsl_vector_get(x, 1);
+ return f_stat(u, v, data);
+}
+
+static double fobj3(const gsl_vector * x, void * data)
+{
+ double u, v,w;
+ u = gsl_vector_get(x, 0);
+ v = gsl_vector_get(x, 1);
+ w = gsl_vector_get(x, 2);
+ return f_stat3(u, v, w, data);
+}
+
+static double fobjN (const gsl_vector * x, void * data)
+{
+ return f_statN(x->data, data);
+}
+
+static void dfobj (const gsl_vector * x, void * params, gsl_vector * g)
+{
+ double u, v, f, duf, dvf;
+ u = gsl_vector_get(x, 0);
+ v = gsl_vector_get(x, 1);
+ df_stat(u, v, f, duf, dvf, params);
+ gsl_vector_set(g, 0, duf);
+ gsl_vector_set(g, 1, dvf);
+}
+
+static void dfobj3 (const gsl_vector * x, void * params, gsl_vector * g)
+{
+ double u, v, w,f, duf, dvf,dwf;
+ u = gsl_vector_get(x, 0);
+ v = gsl_vector_get(x, 1);
+ w = gsl_vector_get(x, 2);
+ df_stat3(u, v, w, f, duf, dvf, dwf, params);
+ gsl_vector_set(g, 0, duf);
+ gsl_vector_set(g, 1, dvf);
+ gsl_vector_set(g, 2, dwf);
+}
+
+static void dfobjN(const gsl_vector * x, void * params, gsl_vector * g)
+{
+ double f;
+ df_statN(x->data, g->data,f,params);
+}
+
+static void fdfobj(const gsl_vector * x, void * params, double * f, gsl_vector * g)
+{
+ double u, v, duf, dvf;
+ u = gsl_vector_get(x, 0);
+ v = gsl_vector_get(x, 1);
+ df_stat(u, v, *f, duf, dvf, params);
+ gsl_vector_set(g, 0, duf);
+ gsl_vector_set(g, 1, dvf);
+}
+
+static void fdfobj3(const gsl_vector * x, void * params, double * f, gsl_vector * g)
+{
+ double u, v, w,duf, dvf,dwf;
+ u = gsl_vector_get(x, 0);
+ v = gsl_vector_get(x, 1);
+ w = gsl_vector_get(x, 2);
+ df_stat3(u, v, w, *f, duf, dvf, dwf, params);
+ gsl_vector_set(g, 0, duf);
+ gsl_vector_set(g, 1, dvf);
+ gsl_vector_set(g, 2, dwf);
+}
+
+static void fdfobjN(const gsl_vector * x, void * params, double * f, gsl_vector * g)
+{
+ df_statN(x->data,g->data,*f,params);
+}
+
+void minimize_2(double (*f) (double, double, void *data),
+ void (*df) (double, double, double &, double &, double &, void *data) ,
+ void *data,int niter,
+ double &U, double &V, double &res)
+{
+ f_stat = f;
+ df_stat = df;
+
+ int iter = 0;
+ int status;
+
+ const gsl_multimin_fdfminimizer_type *T;
+ gsl_multimin_fdfminimizer *s;
+ gsl_vector *x;
+ gsl_multimin_function_fdf my_func;
+
+ my_func.f = &fobj;
+ my_func.df = &dfobj;
+ my_func.fdf = &fdfobj;
+ my_func.n = 2;
+ my_func.params = data;
+
+ x = gsl_vector_alloc (2);
+ gsl_vector_set (x, 0, U);
+ gsl_vector_set (x, 1, V);
+
+ T = gsl_multimin_fdfminimizer_conjugate_fr;
+ s = gsl_multimin_fdfminimizer_alloc (T, 2);
+
+ gsl_multimin_fdfminimizer_set(s, &my_func, x, 0.01, 1e-4);
+
+ do{
+ iter++;
+ status = gsl_multimin_fdfminimizer_iterate(s);
+ if(status)
+ break;
+ status = gsl_multimin_test_gradient(s->gradient, 1e-3);
+ }
+ while (status == GSL_CONTINUE && iter < niter);
+
+ U = gsl_vector_get(s->x, 0);
+ V = gsl_vector_get(s->x, 1);
+ res = s->f;
+ gsl_multimin_fdfminimizer_free (s);
+ gsl_vector_free (x);
+}
+
+void minimize_3(double (*f) (double, double, double,void *data),
+ void (*df) (double , double , double ,
+ double &, double &, double &, double &,
+ void *data) ,
+ void *data,int niter,
+ double &U, double &V, double &W, double &res)
+{
+ f_stat3 = f;
+ df_stat3 = df;
+
+ int iter = 0;
+ int status;
+
+ const gsl_multimin_fdfminimizer_type *T;
+ gsl_multimin_fdfminimizer *s;
+ gsl_vector *x;
+ gsl_multimin_function_fdf my_func;
+
+ my_func.f = &fobj3;
+ my_func.df = &dfobj3;
+ my_func.fdf = &fdfobj3;
+ my_func.n = 3;
+ my_func.params = data;
+
+ x = gsl_vector_alloc (3);
+ gsl_vector_set(x, 0, U);
+ gsl_vector_set(x, 1, V);
+ gsl_vector_set(x, 2, W);
+
+ T = gsl_multimin_fdfminimizer_conjugate_fr;
+ s = gsl_multimin_fdfminimizer_alloc(T, 3);
+
+ gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.01, 1e-4);
+
+ do{
+ iter++;
+ status = gsl_multimin_fdfminimizer_iterate (s);
+ if (status)
+ break;
+ status = gsl_multimin_test_gradient (s->gradient, 1e-3);
+ }
+ while (status == GSL_CONTINUE && iter < niter);
+
+ U = gsl_vector_get (s->x, 0);
+ V = gsl_vector_get (s->x, 1);
+ W = gsl_vector_get (s->x, 2);
+ res = s->f;
+ gsl_multimin_fdfminimizer_free (s);
+ gsl_vector_free (x);
+}
+
+void minimize_N(int N,
+ double (*f) (double*,void *data),
+ void (*df) (double*,double*,double &,void *data) ,
+ void *data,int niter,
+ double *U, double &res)
+{
+ f_statN = f;
+ df_statN = df;
+
+ int iter = 0;
+ int status;
+
+ const gsl_multimin_fdfminimizer_type *T;
+ gsl_multimin_fdfminimizer *s;
+ gsl_vector *x;
+ gsl_multimin_function_fdf my_func;
+
+ my_func.f = &fobjN;
+ my_func.df = &dfobjN;
+ my_func.fdf = &fdfobjN;
+ my_func.n = N;
+ my_func.params = data;
+
+ x = gsl_vector_alloc (N);
+ for (int i=0;i<N;i++)
+ gsl_vector_set (x, i, U[i]);
+
+ T = gsl_multimin_fdfminimizer_conjugate_fr;
+ s = gsl_multimin_fdfminimizer_alloc (T, N);
+
+ gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.01, 1e-4);
+
+ do{
+ iter++;
+ status = gsl_multimin_fdfminimizer_iterate (s);
+
+ if (status)
+ break;
+
+ status = gsl_multimin_test_gradient (s->gradient, 1e-3);
+ }
+ while (status == GSL_CONTINUE && iter < niter);
+
+ for (int i=0;i<N;i++)
+ U[i] = gsl_vector_get (s->x, i);
+ res = s->f;
+ gsl_multimin_fdfminimizer_free (s);
+ gsl_vector_free (x);
+}
+
+#else
+
+int check_gsl()
+{
+ // initilize robust geometric predicates
+ gmsh::exactinit() ;
+ return 1;
+}
+
+double brent(double ax, double bx, double cx,
+ double (*f) (double), double tol, double *xmin)
+{
+ Msg::Error("Gmsh must be compiled with GSL support for brent");
+}
+
+void newt(double x[], int n, int *check,
+ void (*func) (int, double[], double[]))
+{
+ Msg::Error("Gmsh must be compiled with GSL support for newt");
+}
+
+void minimize_2(double (*f) (double, double, void *data),
+ void (*df) (double, double, double &, double &, double &, void *data) ,
+ void *data,int niter,
+ double &U, double &V, double &res)
+{
+ Msg::Error("Gmsh must be compiled with GSL support for minimize_2");
+}
+
+void minimize_3(double (*f) (double, double, double, void *data),
+ void (*df) (double, double, double , double &, double &,
+ double &, double &, void *data) ,
+ void *data,int niter,
+ double &U, double &V, double &W, double &res)
+{
+ Msg::Error("Gmsh must be compiled with GSL support for minimize_3");
+}
+
+void minimize_N(int N, double (*f) (double*, void *data),
+ void (*df) (double*, double*, double &, void *data) ,
+ void *data,int niter,
+ double *, double &res)
+{
+ Msg::Error("Gmsh must be compiled with GSL support for minimize_N");
+}
+
+#endif
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 341abd4728..4dc536ba36 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -9,28 +9,23 @@
#include "NumericEmbedded.h"
int check_gsl();
-void invert_singular_matrix3x3(double MM[3][3], double II[3][3]);
-
-// Numerical routines implemented using either Numerical Recipes or
-// the GSL
double brent(double ax, double bx, double cx,
- double (*f)(double), double tol, double *xmin);
+ double (*f) (double), double tol, double *xmin);
void newt(double x[], int n, int *check,
- void (*vecfunc)(int, double [], double []));
-void minimize_2 (double (*f) (double, double, void *data),
- void (*df) (double, double, double &, double &, double &, void *data) ,
- void *data,int niter,
- double &U, double &V, double &res);
-void minimize_3 (double (*f) (double, double, double, void *data),
- void (*df) (double, double, double , double &, double &,
- double &, double &, void *data),
- void *data,int niter,
- double &U, double &V, double &W, double &res);
-void minimize_N (int N,
- double (*f) (double*, void *data),
- void (*df) (double*, double*, double &, void *data) ,
- void *data,int niter,
- double *, double &res);
+ void (*func) (int, double[], double[]));
+void minimize_2(double (*f) (double, double, void *data),
+ void (*df) (double, double, double &, double &, double &, void *data),
+ void *data,int niter,
+ double &U, double &V, double &res);
+void minimize_3(double (*f) (double, double, double, void *data),
+ void (*df) (double, double, double , double &, double &,
+ double &, double &, void *data),
+ void *data,int niter,
+ double &U, double &V, double &W, double &res);
+void minimize_N(int N, double (*f) (double*, void *data),
+ void (*df) (double*, double*, double &, void *data),
+ void *data,int niter,
+ double *, double &res);
// Robust geometrical predicates
namespace gmsh{
diff --git a/Numeric/NumericEmbedded.cpp b/Numeric/NumericEmbedded.cpp
index c1963fefd1..c625836a55 100644
--- a/Numeric/NumericEmbedded.cpp
+++ b/Numeric/NumericEmbedded.cpp
@@ -8,6 +8,7 @@
// Msg)
#include "NumericEmbedded.h"
+#include "GmshMatrix.h"
#include "GmshMessage.h"
#define SQU(a) ((a)*(a))
@@ -449,3 +450,40 @@ void eigsort(double d[3])
}
}
}
+
+void invert_singular_matrix3x3(double MM[3][3], double II[3][3])
+{
+ int i, j, k, n = 3;
+ double TT[3][3];
+
+ for(i = 1; i <= n; i++) {
+ for(j = 1; j <= n; j++) {
+ II[i - 1][j - 1] = 0.0;
+ TT[i - 1][j - 1] = 0.0;
+ }
+ }
+
+ Double_Matrix M(3, 3), V(3, 3);
+ Double_Vector W(3);
+ for(i = 1; i <= n; i++){
+ for(j = 1; j <= n; j++){
+ M(i - 1, j - 1) = MM[i - 1][j - 1];
+ }
+ }
+ M.svd(V, W);
+ for(i = 1; i <= n; i++) {
+ for(j = 1; j <= n; j++) {
+ double ww = W(i - 1);
+ if(fabs(ww) > 1.e-16) { // singular value precision
+ TT[i - 1][j - 1] += M(j - 1, i - 1) / ww;
+ }
+ }
+ }
+ for(i = 1; i <= n; i++) {
+ for(j = 1; j <= n; j++) {
+ for(k = 1; k <= n; k++) {
+ II[i - 1][j - 1] += V(i - 1, k - 1) * TT[k - 1][j - 1];
+ }
+ }
+ }
+}
diff --git a/Numeric/NumericEmbedded.h b/Numeric/NumericEmbedded.h
index 08523ad8df..1256319121 100644
--- a/Numeric/NumericEmbedded.h
+++ b/Numeric/NumericEmbedded.h
@@ -73,5 +73,6 @@ double InterpolateIso(double *X, double *Y, double *Z,
void gradSimplex(double *x, double *y, double *z, double *v, double *grad);
double ComputeVonMises(double *val);
double ComputeScalarRep(int numComp, double *val);
+void invert_singular_matrix3x3(double MM[3][3], double II[3][3]);
#endif
diff --git a/Numeric/gsl_brent.cpp b/Numeric/gsl_brent.cpp
deleted file mode 100644
index 59cb3be607..0000000000
--- a/Numeric/gsl_brent.cpp
+++ /dev/null
@@ -1,132 +0,0 @@
-// Gmsh - Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-
-#if defined(HAVE_GSL)
-
-#include "GmshMessage.h"
-#include "Numeric.h"
-
-#include <gsl/gsl_errno.h>
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_min.h>
-
-static double (*nrfunc) (double);
-
-double fn1(double x, void *params)
-{
- double val = nrfunc(x);
- return val;
-}
-
-int gsl_min_fminimizer_set_with_values_MOD(gsl_min_fminimizer *s,
- gsl_function *f,
- double x_minimum, double f_minimum,
- double x_lower, double f_lower,
- double x_upper, double f_upper)
-{
- s->function = f;
- s->x_minimum = x_minimum;
- s->x_lower = x_lower;
- s->x_upper = x_upper;
- if (x_lower > x_upper)
- {
- GSL_ERROR ("invalid interval (lower > upper)", GSL_EINVAL);
- }
- s->f_lower = f_lower;
- s->f_upper = f_upper;
- s->f_minimum = f_minimum;
- return (s->type->set) (s->state, s->function,
- x_minimum, f_minimum,
- x_lower, f_lower,
- x_upper, f_upper);
-}
-
-#define MAXITER 100
-
-// Returns the minimum betwen ax and cx to a given tolerance tol using
-// brent's method.
-
-double brent(double ax, double bx, double cx,
- double (*f) (double), double tol, double *xmin)
-{
- // The solver can stall at the following internal tolerance in brent - so
- // this is about the lowest possible tolerance.
- if(tol == 0.) tol = 10.*GSL_SQRT_DBL_EPSILON;
- const double tolsq = tol*tol;
- int status;
- int iter = 0;
- double a, b, c; // a < b < c
- double fa, fb, fc;
- const gsl_min_fminimizer_type *T;
- gsl_min_fminimizer *s;
- gsl_function F;
-
- // gsl wants a<b
- b = bx;
- if(ax < cx) {
- a = ax;
- c = cx;
- }
- else {
- a = cx;
- c = ax;
- }
-
- nrfunc = f;
- F.function = &fn1;
- F.params = 0;
-
- fa = f(a);
- fb = f(b);
- fc = f(c);
- if(fb > fa || fb > fc) {
- if(a < c) {
- b = a;
- fb = fa;
- }
- else {
- b = c;
- fb = fc;
- }
- }
-
- // if a-c < tol, return func(b)
- if(gsl_min_test_interval(a, c, tolsq, tol) == GSL_SUCCESS) {
- *xmin = b;
- return (fb);
- }
-
- T = gsl_min_fminimizer_brent;
- s = gsl_min_fminimizer_alloc(T);
-// gsl_min_fminimizer_set(s, &F, b, a, c); // Standard
- // This mod avoids some error checks so we can have a minimum at an extent.
- gsl_min_fminimizer_set_with_values_MOD(s, &F, b, fb, a, fa, c, fc);
-
- do {
- iter++;
- status = gsl_min_fminimizer_iterate(s);
- if(status)
- break; // solver problem
- a = gsl_min_fminimizer_x_lower(s);
- c = gsl_min_fminimizer_x_upper(s);
-// fb = gsl_min_fminimizer_f_minimum(s);
-// fa = gsl_min_fminimizer_f_lower(s);
-// fc = gsl_min_fminimizer_f_upper(s);
- status = gsl_min_test_interval(a, c, tolsq, tol);
- // ... in case curve becomes too flat - not required?
-// std::fabs(fa - fc) < (std::fabs(fb) + tol)*tol
- }
- while(status == GSL_CONTINUE && iter < MAXITER);
- b = gsl_min_fminimizer_x_minimum(s);
-
- if(status != GSL_SUCCESS)
- Msg::Error("MIN1D not converged: f(%g) = %g", b, fn1(b, NULL));
-
- *xmin = b;
- gsl_min_fminimizer_free(s);
- return fn1(b, NULL);
-}
-
-#endif
diff --git a/Numeric/gsl_min.cpp b/Numeric/gsl_min.cpp
deleted file mode 100644
index 4a425567f5..0000000000
--- a/Numeric/gsl_min.cpp
+++ /dev/null
@@ -1,265 +0,0 @@
-// Gmsh - Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-
-#include "GmshMessage.h"
-
-#if defined(HAVE_GSL)
-
-#include "Numeric.h"
-#include <gsl/gsl_errno.h>
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_multimin.h>
-
-static double (*f_stat) (double, double, void *data);
-static void (*df_stat) (double, double, double &, double &, double &, void *data);
-static double (*f_stat3) (double, double, double, void *data);
-static void (*df_stat3) (double, double, double, double &, double &, double &, double &,void *data);
-static double (*f_statN) (double*, void *data);
-static void (*df_statN) (double*, double *, double &,void *data);
-
-double fobj (const gsl_vector * x, void * data){
- double u, v;
- u = gsl_vector_get(x, 0);
- v = gsl_vector_get(x, 1);
- return f_stat(u,v,data);
-}
-
-
-double fobj3 (const gsl_vector * x, void * data){
- double u, v,w;
- u = gsl_vector_get(x, 0);
- v = gsl_vector_get(x, 1);
- w = gsl_vector_get(x, 2);
- return f_stat3(u,v,w,data);
-}
-
-double fobjN (const gsl_vector * x, void * data){
- return f_statN(x->data,data);
-}
-
-void dfobj (const gsl_vector * x, void * params, gsl_vector * g){
- double u, v, f, duf, dvf;
- u = gsl_vector_get(x, 0);
- v = gsl_vector_get(x, 1);
- df_stat(u,v,f, duf,dvf,params);
- gsl_vector_set(g, 0, duf);
- gsl_vector_set(g, 1, dvf);
-}
-
-void dfobj3 (const gsl_vector * x, void * params, gsl_vector * g){
- double u, v, w,f, duf, dvf,dwf;
- u = gsl_vector_get(x, 0);
- v = gsl_vector_get(x, 1);
- w = gsl_vector_get(x, 2);
- df_stat3(u,v,w,f, duf,dvf,dwf,params);
- gsl_vector_set(g, 0, duf);
- gsl_vector_set(g, 1, dvf);
- gsl_vector_set(g, 2, dwf);
-}
-
-void dfobjN (const gsl_vector * x, void * params, gsl_vector * g){
- double f;
- df_statN(x->data, g->data,f,params);
-}
-
-void fdfobj (const gsl_vector * x, void * params, double * f, gsl_vector * g){
- double u, v, duf, dvf;
- u = gsl_vector_get(x, 0);
- v = gsl_vector_get(x, 1);
- df_stat(u,v,*f, duf,dvf,params);
- gsl_vector_set(g, 0, duf);
- gsl_vector_set(g, 1, dvf);
-}
-
-void fdfobj3 (const gsl_vector * x, void * params, double * f, gsl_vector * g){
- double u, v, w,duf, dvf,dwf;
- u = gsl_vector_get(x, 0);
- v = gsl_vector_get(x, 1);
- w = gsl_vector_get(x, 2);
- df_stat3(u,v,w,*f, duf,dvf,dwf,params);
- gsl_vector_set(g, 0, duf);
- gsl_vector_set(g, 1, dvf);
- gsl_vector_set(g, 2, dwf);
-}
-
-void fdfobjN (const gsl_vector * x, void * params, double * f, gsl_vector * g){
- df_statN(x->data,g->data,*f,params);
-}
-
-void minimize_2 ( double (*f) (double, double, void *data),
- void (*df) (double, double, double &, double &, double &, void *data) ,
- void *data,int niter,
- double &U, double &V, double &res){
- f_stat = f;
- df_stat = df;
-
- int iter = 0;
- int status;
-
- const gsl_multimin_fdfminimizer_type *T;
- gsl_multimin_fdfminimizer *s;
- gsl_vector *x;
- gsl_multimin_function_fdf my_func;
-
- my_func.f = &fobj;
- my_func.df = &dfobj;
- my_func.fdf = &fdfobj;
- my_func.n = 2;
- my_func.params = data;
-
- x = gsl_vector_alloc (2);
- gsl_vector_set (x, 0, U);
- gsl_vector_set (x, 1, V);
-
- T = gsl_multimin_fdfminimizer_conjugate_fr;
- s = gsl_multimin_fdfminimizer_alloc (T, 2);
-
- gsl_multimin_fdfminimizer_set(s, &my_func, x, 0.01, 1e-4);
-
- do{
- iter++;
- status = gsl_multimin_fdfminimizer_iterate(s);
- if(status)
- break;
- status = gsl_multimin_test_gradient(s->gradient, 1e-3);
- }
- while (status == GSL_CONTINUE && iter < niter);
-
- U = gsl_vector_get(s->x, 0);
- V = gsl_vector_get(s->x, 1);
- res = s->f;
- gsl_multimin_fdfminimizer_free (s);
- gsl_vector_free (x);
-}
-
-void minimize_3(double (*f) (double, double, double,void *data),
- void (*df) (double , double , double ,
- double &, double &, double &, double &,
- void *data) ,
- void *data,int niter,
- double &U, double &V, double &W, double &res)
-{
- f_stat3 = f;
- df_stat3 = df;
-
- int iter = 0;
- int status;
-
- const gsl_multimin_fdfminimizer_type *T;
- gsl_multimin_fdfminimizer *s;
- gsl_vector *x;
- gsl_multimin_function_fdf my_func;
-
- my_func.f = &fobj3;
- my_func.df = &dfobj3;
- my_func.fdf = &fdfobj3;
- my_func.n = 3;
- my_func.params = data;
-
- x = gsl_vector_alloc (3);
- gsl_vector_set(x, 0, U);
- gsl_vector_set(x, 1, V);
- gsl_vector_set(x, 2, W);
-
- T = gsl_multimin_fdfminimizer_conjugate_fr;
- s = gsl_multimin_fdfminimizer_alloc(T, 3);
-
- gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.01, 1e-4);
-
- do{
- iter++;
- status = gsl_multimin_fdfminimizer_iterate (s);
- if (status)
- break;
- status = gsl_multimin_test_gradient (s->gradient, 1e-3);
- }
- while (status == GSL_CONTINUE && iter < niter);
-
- U = gsl_vector_get (s->x, 0);
- V = gsl_vector_get (s->x, 1);
- W = gsl_vector_get (s->x, 2);
- res = s->f;
- gsl_multimin_fdfminimizer_free (s);
- gsl_vector_free (x);
-}
-
-void minimize_N(int N,
- double (*f) (double*,void *data),
- void (*df) (double*,double*,double &,void *data) ,
- void *data,int niter,
- double *U, double &res){
- f_statN = f;
- df_statN = df;
-
- int iter = 0;
- int status;
-
- const gsl_multimin_fdfminimizer_type *T;
- gsl_multimin_fdfminimizer *s;
- gsl_vector *x;
- gsl_multimin_function_fdf my_func;
-
- my_func.f = &fobjN;
- my_func.df = &dfobjN;
- my_func.fdf = &fdfobjN;
- my_func.n = N;
- my_func.params = data;
-
- x = gsl_vector_alloc (N);
- for (int i=0;i<N;i++)
- gsl_vector_set (x, i, U[i]);
-
- T = gsl_multimin_fdfminimizer_conjugate_fr;
- s = gsl_multimin_fdfminimizer_alloc (T, N);
-
- gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.01, 1e-4);
-
- do{
- iter++;
- status = gsl_multimin_fdfminimizer_iterate (s);
-
- if (status)
- break;
-
- status = gsl_multimin_test_gradient (s->gradient, 1e-3);
- }
- while (status == GSL_CONTINUE && iter < niter);
-
- for (int i=0;i<N;i++)
- U[i] = gsl_vector_get (s->x, i);
- res = s->f;
- gsl_multimin_fdfminimizer_free (s);
- gsl_vector_free (x);
-
-}
-
-#else
-
-void minimize_2 ( double (*f) (double, double, void *data),
- void (*df) (double, double, double &, double &, double &, void *data) ,
- void *data,int niter,
- double &U, double &V, double &res)
-{
- Msg::Error("Gmsh must be compiled with GSL support for minimize_2");
-}
-
-void minimize_3 ( double (*f) (double, double, double, void *data),
- void (*df) (double, double, double , double &, double &, double &, double &, void *data) ,
- void *data,int niter,
- double &U, double &V, double &W, double &res)
-{
- Msg::Error("Gmsh must be compiled with GSL support for minimize_3");
-}
-
-void minimize_N (int N,
- double (*f) (double*, void *data),
- void (*df) (double*, double*, double &, void *data) ,
- void *data,int niter,
- double *, double &res)
-{
- Msg::Error("Gmsh must be compiled with GSL support for minimize_N");
-}
-
-#endif
diff --git a/Numeric/gsl_newt.cpp b/Numeric/gsl_newt.cpp
deleted file mode 100644
index f8dfc1e3ea..0000000000
--- a/Numeric/gsl_newt.cpp
+++ /dev/null
@@ -1,108 +0,0 @@
-// Gmsh - Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-//
-// Contributor(s):
-// Nicolas Tardieu
-//
-
-// This implements a Newton method using the GSL.
-
-#include "GmshMessage.h"
-
-#if defined(HAVE_GSL)
-
-#include "Numeric.h"
-
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_errno.h>
-#include <gsl/gsl_multiroots.h>
-#include <gsl/gsl_linalg.h>
-
-#define MAX_DIM_NEWT 10
-#define MAXITER 100
-#define PREC 1.e-8
-
-int nrdim;
-double nru[MAX_DIM_NEWT], nrv[MAX_DIM_NEWT];
-static void (*nrfunc) (int n, double x[], double y[]);
-struct gsl_dummy{ int i; };
-
-void convert_vector_from_gsl(const gsl_vector * gv, double *v)
-{
- int i, m;
- m = gv->size;
- for(i = 0; i < m; i++) {
- v[i + 1] = gsl_vector_get(gv, i);
- }
-}
-
-void convert_vector_to_gsl(double *v, int n, gsl_vector * gv)
-{
- int i;
- for(i = 0; i < n; i++) {
- gsl_vector_set(gv, i, v[i + 1]);
- }
-}
-
-int gslfunc(const gsl_vector * xx, void *params, gsl_vector * f)
-{
- convert_vector_from_gsl(xx, nru);
- (*nrfunc) (nrdim, nru, nrv);
- // Msg::Info("f(%lf,%lf) = %lf %lf\n",nru[1],nru[2],nrv[1],nrv[2]);
- convert_vector_to_gsl(nrv, nrdim, f);
- return GSL_SUCCESS;
-}
-
-// Warning: for compatibility with the old newt from NR, x[] and the
-// arguments of func() are indexed from 1 to n...
-
-void newt(double x[], int n, int *check,
- void (*func) (int, double[], double[]))
-{
- const gsl_multiroot_fsolver_type *T;
- gsl_multiroot_fsolver *s;
- int status;
- size_t iter = 0;
- struct gsl_dummy p = { 1 };
- gsl_multiroot_function f = { &gslfunc, n, &p };
- gsl_vector *xx = gsl_vector_alloc(n);
-
- if(n > MAX_DIM_NEWT - 1)
- Msg::Fatal("Maximum Newton dimension exceeded\n");
- nrdim = n;
-
- nrfunc = func;
- convert_vector_to_gsl(x, n, xx);
-
- T = gsl_multiroot_fsolver_hybrid;
- s = gsl_multiroot_fsolver_alloc(T, n);
- gsl_multiroot_fsolver_set(s, &f, xx);
-
- do {
- iter++;
- status = gsl_multiroot_fsolver_iterate(s);
- // Msg::Info("status %d %d %d %lf %lf\n",
- // status,n,iter,gsl_vector_get(s->x,0),gsl_vector_get(s->x,1));
- if(status)
- break; // solver problem
- status = gsl_multiroot_test_residual(s->f, n * PREC);
- }
- while(status == GSL_CONTINUE && iter < MAXITER);
-
- if(status == GSL_CONTINUE) {
- *check = 1; // problem !!!
- }
- else {
- // Msg::Info("status %d %d %d %lf %lf\n",
- // status,n,iter,gsl_vector_get(s->x,0),gsl_vector_get(s->x,1));
- convert_vector_from_gsl(s->x, x);
- *check = 0; // converged
- }
-
- gsl_multiroot_fsolver_free(s);
- gsl_vector_free(xx);
-}
-
-#endif
diff --git a/Post/Makefile b/Post/Makefile
index 36d3fd23bb..8f1f9363a8 100644
--- a/Post/Makefile
+++ b/Post/Makefile
@@ -126,7 +126,7 @@ adaptiveData${OBJEXT}: adaptiveData.cpp adaptiveData.h ../Numeric/GmshMatrix.h \
../Common/GmshMessage.h ../Plugin/Plugin.h ../Common/Options.h \
../Post/ColorTable.h ../Post/PView.h ../Geo/SPoint3.h \
../Post/PViewDataList.h ../Post/PViewData.h ../Geo/SBoundingBox3d.h \
- ../Geo/SPoint3.h ../Common/ListUtils.h
+ ../Geo/SPoint3.h ../Common/ListUtils.h ../Common/OS.h
OctreePost${OBJEXT}: OctreePost.cpp ../Common/Octree.h \
../Common/OctreeInternals.h OctreePost.h ../Common/ListUtils.h PView.h \
../Geo/SPoint3.h PViewDataList.h PViewData.h ../Geo/SBoundingBox3d.h \
diff --git a/configure b/configure
index cb2a8f1277..6b661cfdcc 100755
--- a/configure
+++ b/configure
@@ -1316,7 +1316,8 @@ Optional Packages:
--with-mpi-prefix=PFX prefix where MPI is installed
--with-fftw3-prefix=PFX prefix where FFTW3 is installed
--with-fm-prefix=PFX prefix where FourierModel is installed
- --with-blas-prefix=PFX prefix where BLAS (and LAPACK) are installed
+ --with-blas-lapack-prefix=PFX
+ prefix where BLAS and LAPACK are installed
Some influential environment variables:
CC C compiler command
@@ -1841,9 +1842,9 @@ if test "${with_fm_prefix+set}" = set; then
fi
-# Check whether --with-blas-prefix was given.
-if test "${with_blas_prefix+set}" = set; then
- withval=$with_blas_prefix; BLAS_PREFIX=$withval
+# Check whether --with-blas-lapack-prefix was given.
+if test "${with_blas_lapack_prefix+set}" = set; then
+ withval=$with_blas_lapack_prefix; BLAS_LAPACK_PREFIX=$withval
fi
@@ -5128,43 +5129,6 @@ fi
fi
fi
fi
-if test "x${GSL}" != "xyes"; then
- { echo "$as_me:$LINENO: checking for ./contrib/NR/dsvdcmp.cpp" >&5
-echo $ECHO_N "checking for ./contrib/NR/dsvdcmp.cpp... $ECHO_C" >&6; }
-if test "${ac_cv_file___contrib_NR_dsvdcmp_cpp+set}" = set; then
- echo $ECHO_N "(cached) $ECHO_C" >&6
-else
- test "$cross_compiling" = yes &&
- { { echo "$as_me:$LINENO: error: cannot check for file existence when cross compiling" >&5
-echo "$as_me: error: cannot check for file existence when cross compiling" >&2;}
- { (exit 1); exit 1; }; }
-if test -r "./contrib/NR/dsvdcmp.cpp"; then
- ac_cv_file___contrib_NR_dsvdcmp_cpp=yes
-else
- ac_cv_file___contrib_NR_dsvdcmp_cpp=no
-fi
-fi
-{ echo "$as_me:$LINENO: result: $ac_cv_file___contrib_NR_dsvdcmp_cpp" >&5
-echo "${ECHO_T}$ac_cv_file___contrib_NR_dsvdcmp_cpp" >&6; }
-if test $ac_cv_file___contrib_NR_dsvdcmp_cpp = yes; then
- NR="yes"
-fi
-
- if test "x${NR}" = "xyes"; then
- GMSH_DIRS="${GMSH_DIRS} contrib/NR"
- GMSH_LIBS="${GMSH_LIBS} -lGmshNR"
- echo "********************************************************************"
- echo " You are building a non-free version of Gmsh, using code copyright"
- echo " (C) 1986-92 Numerical Recipes Software."
- echo " To use the GSL instead, run configure again with the --enable-gsl"
- echo " option."
- echo "********************************************************************"
- else
- { { echo "$as_me:$LINENO: error: Could not find GSL, aborting." >&5
-echo "$as_me: error: Could not find GSL, aborting." >&2;}
- { (exit 1); exit 1; }; }
- fi
-fi
if test "x$enable_fm" != "xno"; then
if test "x${FM_PREFIX}" != "x"; then
@@ -5316,8 +5280,8 @@ echo "$as_me: WARNING: Could not find FFTW3: disabling FourierModel." >&2;}
fi
fi
-if test "x${BLAS_PREFIX}" != "x"; then
- LDFLAGS="-L${BLAS_PREFIX} -L${BLAS_PREFIX}/lib ${LDFLAGS}"
+if test "x${BLAS_LAPACK_PREFIX}" != "x"; then
+ LDFLAGS="-L${BLAS_LAPACK_PREFIX} -L${BLAS_LAPACK_PREFIX}/lib ${LDFLAGS}"
fi
{ echo "$as_me:$LINENO: checking for cblas_dgemm in -lcblas" >&5
echo $ECHO_N "checking for cblas_dgemm in -lcblas... $ECHO_C" >&6; }
@@ -5452,8 +5416,8 @@ fi
fi
if test "x${CBLAS}" = "xyes"; then
- if test "x${BLAS_PREFIX}" != "x"; then
- BLAS_PATH="-L${BLAS_PREFIX} -L${BLAS_PREFIX}/lib"
+ if test "x${BLAS_LAPACK_PREFIX}" != "x"; then
+ BLAS_PATH="-L${BLAS_LAPACK_PREFIX} -L${BLAS_LAPACK_PREFIX}/lib"
fi
FLAGS="${FLAGS} -DHAVE_CBLAS"
else
@@ -5463,7 +5427,7 @@ else
fi
fi
-if test "x${FM}" = "xyes"; then
+if test "x${FM}" = "xyes" -o "x${GSL}" != "xyes"; then
ac_ext=f
ac_compile='$F77 -c $FFLAGS conftest.$ac_ext >&5'
ac_link='$F77 -o conftest$ac_exeext $FFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5'
diff --git a/configure.in b/configure.in
index 296cb4a024..3f7eb290f4 100644
--- a/configure.in
+++ b/configure.in
@@ -69,10 +69,10 @@ AC_ARG_WITH(fm-prefix,
AC_HELP_STRING([--with-fm-prefix=PFX],
[prefix where FourierModel is installed]),
[FM_PREFIX=$withval])
-AC_ARG_WITH(blas-prefix,
- AC_HELP_STRING([--with-blas-prefix=PFX],
- [prefix where BLAS (and LAPACK) are installed]),
- [BLAS_PREFIX=$withval])
+AC_ARG_WITH(blas-lapack-prefix,
+ AC_HELP_STRING([--with-blas-lapack-prefix=PFX],
+ [prefix where BLAS and LAPACK are installed]),
+ [BLAS_LAPACK_PREFIX=$withval])
dnl Parse '--enable' command line options
AC_ARG_ENABLE(gsl,
@@ -667,22 +667,6 @@ if test "x$enable_gsl" != "xno"; then
fi
fi
fi
-if test "x${GSL}" != "xyes"; then
- dnl Check if non-free numerical recipes routines are in the tree
- AC_CHECK_FILE(./contrib/NR/dsvdcmp.cpp,NR="yes")
- if test "x${NR}" = "xyes"; then
- GMSH_DIRS="${GMSH_DIRS} contrib/NR"
- GMSH_LIBS="${GMSH_LIBS} -lGmshNR"
- echo "********************************************************************"
- echo " You are building a non-free version of Gmsh, using code copyright"
- echo " (C) 1986-92 Numerical Recipes Software."
- echo " To use the GSL instead, run configure again with the --enable-gsl"
- echo " option."
- echo "********************************************************************"
- else
- AC_MSG_ERROR([Could not find GSL, aborting.])
- fi
-fi
dnl Check for FourierModel
if test "x$enable_fm" != "xno"; then
@@ -718,16 +702,16 @@ if test "x$enable_fm" != "xno"; then
fi
dnl Check for C version of BLAS
-if test "x${BLAS_PREFIX}" != "x"; then
- LDFLAGS="-L${BLAS_PREFIX} -L${BLAS_PREFIX}/lib ${LDFLAGS}"
+if test "x${BLAS_LAPACK_PREFIX}" != "x"; then
+ LDFLAGS="-L${BLAS_LAPACK_PREFIX} -L${BLAS_LAPACK_PREFIX}/lib ${LDFLAGS}"
fi
AC_CHECK_LIB(cblas,cblas_dgemm,CBLAS="yes" BLAS_LIBS="-lcblas")
if test "x${CBLAS}" != "xyes"; then
AC_CHECK_LIB(cblas,cblas_dgemm,CBLAS="yes" BLAS_LIBS="-lcblas -latlas",[],-latlas)
fi
if test "x${CBLAS}" = "xyes"; then
- if test "x${BLAS_PREFIX}" != "x"; then
- BLAS_PATH="-L${BLAS_PREFIX} -L${BLAS_PREFIX}/lib"
+ if test "x${BLAS_LAPACK_PREFIX}" != "x"; then
+ BLAS_PATH="-L${BLAS_LAPACK_PREFIX} -L${BLAS_LAPACK_PREFIX}/lib"
fi
FLAGS="${FLAGS} -DHAVE_CBLAS"
else
@@ -738,9 +722,9 @@ else
fi
fi
-dnl Check for Fortran version of blas and lapack (only used in
-dnl FourierModel at the moment)
-if test "x${FM}" = "xyes"; then
+dnl Check for Fortran version of blas and lapack (only used when not
+dnl using GSL, or of FourierModel is linked in)
+if test "x${FM}" = "xyes" -o "x${GSL}" != "xyes"; then
AC_PROG_F77
case "${F77}" in
*gfortran*)
diff --git a/contrib/NR/Makefile b/contrib/NR/Makefile
deleted file mode 100644
index a56f425626..0000000000
--- a/contrib/NR/Makefile
+++ /dev/null
@@ -1,61 +0,0 @@
-# Gmsh - Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
-#
-# See the LICENSE.txt file for license information. Please report all
-# bugs and problems to <gmsh@geuz.org>.
-
-include ../../variables
-
-LIB = ../../lib/libGmshNR${LIBEXT}
-
-INC = ${DASH}I../../Common ${DASH}I.
-
-# don't optimize this library: there are some problems with gcc...
-CFLAGS = ${FLAGS} ${INC}
-
-SRC = brent.cpp\
- dpythag.cpp\
- dsvdcmp.cpp\
- fdjac.cpp\
- fmin.cpp\
- lnsrch.cpp\
- lubksb.cpp\
- ludcmp.cpp\
- newt.cpp\
- nrutil.cpp
-
-OBJ = ${SRC:.cpp=${OBJEXT}}
-
-.SUFFIXES: ${OBJEXT} .cpp
-
-${LIB}: ${OBJ}
- ${AR} ${ARFLAGS}${LIB} ${OBJ}
- ${RANLIB} ${LIB}
-
-cpobj: ${OBJ}
- cp -f ${OBJ} ../../lib/
-
-.cpp${OBJEXT}:
- ${CXX} ${CFLAGS} ${DASH}c $<
-
-clean:
- ${RM} *.o *.obj
-
-depend:
- (sed '/^# DO NOT DELETE THIS LINE/q' Makefile && \
- ${CXX} -MM ${CFLAGS} ${SRC} | sed 's/.o:/$${OBJEXT}:/g' \
- ) > Makefile.new
- cp Makefile Makefile.bak
- cp Makefile.new Makefile
- rm -f Makefile.new
-
-# DO NOT DELETE THIS LINE
-brent${OBJEXT}: brent.cpp nrutil.h
-dpythag${OBJEXT}: dpythag.cpp nrutil.h
-dsvdcmp${OBJEXT}: dsvdcmp.cpp nrutil.h
-fdjac${OBJEXT}: fdjac.cpp nrutil.h
-fmin${OBJEXT}: fmin.cpp nrutil.h
-lnsrch${OBJEXT}: lnsrch.cpp nrutil.h
-lubksb${OBJEXT}: lubksb.cpp
-ludcmp${OBJEXT}: ludcmp.cpp nrutil.h
-newt${OBJEXT}: newt.cpp nrutil.h
-nrutil${OBJEXT}: nrutil.cpp ../../Common/GmshMessage.h
diff --git a/contrib/NR/brent.cpp b/contrib/NR/brent.cpp
deleted file mode 100644
index 8dc82e8813..0000000000
--- a/contrib/NR/brent.cpp
+++ /dev/null
@@ -1,77 +0,0 @@
-#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);
-
-/* This file has been modified for inclusion in Gmsh (float->double) */
-
-double brent(double ax, double bx, double cx, double (*f)(double), double tol,
- double *xmin)
-{
- int iter;
- double a,b,d,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/contrib/NR/dpythag.cpp b/contrib/NR/dpythag.cpp
deleted file mode 100644
index 4585295c43..0000000000
--- a/contrib/NR/dpythag.cpp
+++ /dev/null
@@ -1,14 +0,0 @@
-#include <math.h>
-#define NRANSI
-#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+DSQR(absb/absa));
- else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+DSQR(absa/absb)));
-}
-#undef NRANSI
-/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
diff --git a/contrib/NR/dsvdcmp.cpp b/contrib/NR/dsvdcmp.cpp
deleted file mode 100644
index 81317ec7db..0000000000
--- a/contrib/NR/dsvdcmp.cpp
+++ /dev/null
@@ -1,183 +0,0 @@
-#include <math.h>
-#define NRANSI
-#include "nrutil.h"
-
-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 == 30) nrerror("no convergence in 30 dsvdcmp 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);
-}
-#undef NRANSI
-/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
diff --git a/contrib/NR/fdjac.cpp b/contrib/NR/fdjac.cpp
deleted file mode 100644
index 7767ff1bc4..0000000000
--- a/contrib/NR/fdjac.cpp
+++ /dev/null
@@ -1,29 +0,0 @@
-#include <math.h>
-#define NRANSI
-#include "nrutil.h"
-#define EPS 1.0e-4
-
-/* This file has been modified for inclusion in Gmsh (float->double) */
-
-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
-#undef NRANSI
-/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
diff --git a/contrib/NR/fmin.cpp b/contrib/NR/fmin.cpp
deleted file mode 100644
index f461aa7a05..0000000000
--- a/contrib/NR/fmin.cpp
+++ /dev/null
@@ -1,20 +0,0 @@
-#define NRANSI
-#include "nrutil.h"
-
-/* This file has been modified for inclusion in Gmsh (float->double) */
-
-extern int nn;
-extern double *fvec;
-extern void (*nrfuncv)(int n, double v[], double f[]);
-
-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;
-}
-#undef NRANSI
-/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
diff --git a/contrib/NR/lnsrch.cpp b/contrib/NR/lnsrch.cpp
deleted file mode 100644
index e4a5346a5f..0000000000
--- a/contrib/NR/lnsrch.cpp
+++ /dev/null
@@ -1,65 +0,0 @@
-#include <math.h>
-#define NRANSI
-#include "nrutil.h"
-#define ALF 1.0e-4
-#define TOLX 1.0e-7
-
-/* This file has been modified for inclusion in Gmsh (float->double) */
-
-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])/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
-/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
diff --git a/contrib/NR/lubksb.cpp b/contrib/NR/lubksb.cpp
deleted file mode 100644
index 202c826e2f..0000000000
--- a/contrib/NR/lubksb.cpp
+++ /dev/null
@@ -1,23 +0,0 @@
-/* This file has been modified for inclusion in Gmsh (float->double) */
-
-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];
- }
-}
-/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
diff --git a/contrib/NR/ludcmp.cpp b/contrib/NR/ludcmp.cpp
deleted file mode 100644
index 8f4e0c17fc..0000000000
--- a/contrib/NR/ludcmp.cpp
+++ /dev/null
@@ -1,60 +0,0 @@
-#include <math.h>
-#define NRANSI
-#include "nrutil.h"
-#define TINY 1.0e-20;
-
-/* This file has been modified for inclusion in Gmsh (float->double) */
-
-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
-#undef NRANSI
-/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
diff --git a/contrib/NR/mnbrak.cpp b/contrib/NR/mnbrak.cpp
deleted file mode 100644
index 9aedccc886..0000000000
--- a/contrib/NR/mnbrak.cpp
+++ /dev/null
@@ -1,67 +0,0 @@
-#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);
-
-/* This file has been modified for inclusion in Gmsh (float->double) */
-
-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
-/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
diff --git a/contrib/NR/newt.cpp b/contrib/NR/newt.cpp
deleted file mode 100644
index b09482177e..0000000000
--- a/contrib/NR/newt.cpp
+++ /dev/null
@@ -1,92 +0,0 @@
-#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
-
-/* This file has been modified for inclusion in Gmsh (float->double) */
-
-int nn;
-double *fvec;
-void (*nrfuncv)(int n, double v[], double f[]);
-#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*FMAX(sqrt(sum),(double)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
-/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
diff --git a/contrib/NR/nrutil.cpp b/contrib/NR/nrutil.cpp
deleted file mode 100644
index 0666bc287a..0000000000
--- a/contrib/NR/nrutil.cpp
+++ /dev/null
@@ -1,293 +0,0 @@
-// This file has been modified for inclusion in Gmsh
-
-#include <stdio.h>
-#include <stddef.h>
-#include <stdlib.h>
-#include "GmshMessage.h"
-
-#define NR_END 1
-#define FREE_ARG char*
-
-void nrerror(char error_text[])
-/* Numerical Recipes standard error handler */
-{
- Msg::Error("%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/contrib/NR/nrutil.h b/contrib/NR/nrutil.h
deleted file mode 100644
index 25aa402060..0000000000
--- a/contrib/NR/nrutil.h
+++ /dev/null
@@ -1,43 +0,0 @@
-// This file has been modified for inclusion in GetDP
-
-#ifndef _NR_UTILS_H_
-#define _NR_UTILS_H_
-
-#include <math.h>
-
-#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
-#define DSQR(a) ((a)*(a))
-#define SQR(a) ((a)*(a))
-#define DMIN(a,b) ((a)<(b) ? (a) : (b))
-#define DMAX(a,b) ((a)>(b) ? (a) : (b))
-#define FMIN(a,b) ((a)<(b) ? (a) : (b))
-#define FMAX(a,b) ((a)>(b) ? (a) : (b))
-#define IMIN(a,b) ((a)<(b) ? (a) : (b))
-#define IMAX(a,b) ((a)>(b) ? (a) : (b))
-
-void nrerror(char error_text[]);
-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
--
GitLab