From bf682b25567401fd2e9c4c3687643a496baffd94 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 25 Oct 2009 16:28:33 +0000
Subject: [PATCH] - use lapack for all eigenvector computations - more work on
slepc integration
---
CMakeLists.txt | 9 +-
Common/GmshConfig.h.in | 1 +
Geo/GFace.cpp | 43 +-
Geo/SOrientedBoundingBox.cpp | 2 +-
Geo/STensor3.cpp | 2 -
Geo/STensor3.h | 103 ++--
Mesh/HighOrder.cpp | 4 +-
Mesh/highOrderSmoother.cpp | 2 +-
Numeric/CMakeLists.txt | 1 -
Numeric/EigSolve.cpp | 871 ---------------------------------
Numeric/EigSolve.h | 17 -
Numeric/Numeric.cpp | 2 +-
Numeric/fullMatrix.cpp | 58 +--
Numeric/fullMatrix.h | 12 +-
Numeric/functionSpace.cpp | 2 +-
Plugin/Eigenvectors.cpp | 38 +-
Solver/convexCombinationTerm.h | 2 +-
Solver/crossConfTerm.h | 2 +-
Solver/eigenSolve.cpp | 104 ++++
Solver/eigenSolve.h | 47 ++
Solver/elasticityTerm.cpp | 8 +-
Solver/helmholtzTerm.h | 2 +-
Solver/linearSystemFull.h | 4 +-
23 files changed, 296 insertions(+), 1040 deletions(-)
delete mode 100644 Numeric/EigSolve.cpp
delete mode 100644 Numeric/EigSolve.h
create mode 100644 Solver/eigenSolve.cpp
create mode 100644 Solver/eigenSolve.h
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 81dd8af9af..f2346500de 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -562,7 +562,14 @@ if(ENABLE_PETSC)
endif(ENABLE_PETSC)
if(HAVE_PETSC AND ENABLE_SLEPC)
- message(STATUS "Checking for SLEPc... TODO")
+ set(ENV_SLEPC_DIR $ENV{SLEPC_DIR})
+ find_library(SLEPC_LIB slepc PATHS ${ENV_SLEPC_DIR}/${ENV_PETSC_ARCH}/lib)
+ if(SLEPC_LIB)
+ set(HAVE_SLEPC TRUE)
+ list(APPEND CONFIG_OPTIONS "SLEPc")
+ list(APPEND EXTERNAL_LIBRARIES ${SLEPC_LIB})
+ list(APPEND EXTERNAL_INCLUDES ${ENV_SLEPC_DIR}/include)
+ endif(SLEPC_LIB)
endif(HAVE_PETSC AND ENABLE_SLEPC)
if(ENABLE_OCC)
diff --git a/Common/GmshConfig.h.in b/Common/GmshConfig.h.in
index c13e1f8a84..2fadf56577 100644
--- a/Common/GmshConfig.h.in
+++ b/Common/GmshConfig.h.in
@@ -38,6 +38,7 @@
#cmakedefine HAVE_OSMESA
#cmakedefine HAVE_PETSC
#cmakedefine HAVE_QT
+#cmakedefine HAVE_SLEPSC
#cmakedefine HAVE_SOLVER
#cmakedefine HAVE_TAUCS
#cmakedefine HAVE_TETGEN
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 97628868ae..e257f2edd7 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -15,7 +15,6 @@
#include "VertexArray.h"
#include "fullMatrix.h"
#include "Numeric.h"
-#include "EigSolve.h"
#include "GaussLegendre1D.h"
#include "Context.h"
@@ -641,33 +640,33 @@ void GFace::getMetricEigenVectors(const SPoint2 ¶m,
inv_form1[1][0] = inv_form1[0][1] = -1 * inv_det_form1 * form1[0][1];
// N = (inverse of form1) X (form2)
- double N[4]; // { N00 N01 N10 N11 }
- N[0] = inv_form1[0][0] * form2[0][0] + inv_form1[0][1] * form2[1][0];
- N[1] = inv_form1[0][0] * form2[0][1] + inv_form1[0][1] * form2[1][1];
- N[2] = inv_form1[1][0] * form2[0][0] + inv_form1[1][1] * form2[1][0];
- N[3] = inv_form1[1][0] * form2[0][1] + inv_form1[1][1] * form2[1][1];
+ fullMatrix<double> N(2, 2);
+ N(0, 0) = inv_form1[0][0] * form2[0][0] + inv_form1[0][1] * form2[1][0];
+ N(0, 1) = inv_form1[0][0] * form2[0][1] + inv_form1[0][1] * form2[1][1];
+ N(1, 0) = inv_form1[1][0] * form2[0][0] + inv_form1[1][1] * form2[1][0];
+ N(1, 1) = inv_form1[1][0] * form2[0][1] + inv_form1[1][1] * form2[1][1];
// eigen values and vectors of N
- int work1[2];
- double work2[2];
- double eigValI[2];
- if (EigSolve(2, 2, N, eigVal, eigValI, eigVec, work1, work2) != 1) {
+ fullMatrix<double> vl(2, 2), vr(2, 2);
+ fullVector<double> dr(2), di(2);
+ if(N.eig(dr, di, vl, vr, true)){
+ eigVal[0] = fabs(dr(0));
+ eigVal[1] = fabs(dr(1));
+ eigVec[0] = vr(0, 0);
+ eigVec[1] = vr(1, 0);
+ eigVec[2] = vr(0, 1);
+ eigVec[3] = vr(1, 1);
+ }
+ else{
Msg::Error("Problem in eigen vectors computation");
- Msg::Error(" N: %f %f %f %f", N[0], N[1], N[2], N[3]);
- Msg::Error(" * Eigen values:");
- Msg::Error(" %f + i * %f, %f + i * %f",
- eigVal[0], eigValI[0], eigVal[1], eigValI[1]);
- Msg::Error(" * Eigen vectors (trust it only if eigen values are real):");
- Msg::Error(" ( %f, %f ), ( %f, %f )",
- eigVec[0], eigVec[2], eigVec[1], eigVec[3]);
+ Msg::Error(" N = [ %f %f ]", N(0, 0), N(0, 1));
+ Msg::Error(" [ %f %f ]", N(1, 0), N(1, 1));
+ for(int i = 0; i < 2; i++) eigVal[i] = 0.;
+ for(int i = 0; i < 4; i++) eigVec[i] = 0.;
}
- if (fabs(eigValI[0]) > 1.e-12 || fabs(eigValI[1]) > 1.e-12) {
+ if (fabs(di(0)) > 1.e-12 || fabs(di(1)) > 1.e-12) {
Msg::Error("Found imaginary eigenvalues");
}
-
- eigVal[0] = fabs(eigVal[0]);
- eigVal[1] = fabs(eigVal[1]);
- EigSort(2, eigVal, eigValI, eigVec);
}
void GFace::XYZtoUV(const double X, const double Y, const double Z,
diff --git a/Geo/SOrientedBoundingBox.cpp b/Geo/SOrientedBoundingBox.cpp
index cd1431d782..7a9dad595e 100644
--- a/Geo/SOrientedBoundingBox.cpp
+++ b/Geo/SOrientedBoundingBox.cpp
@@ -258,7 +258,7 @@ SOrientedBoundingBox* SOrientedBoundingBox::buildOBB(std::vector<SPoint3> vertic
fullMatrix<double> right_eigv(3,3);
fullVector<double> real_eig(3);
fullVector<double> img_eig(3);
- covariance.eig(left_eigv, real_eig, img_eig, right_eigv,true);
+ covariance.eig(real_eig, img_eig, left_eigv, right_eigv,true);
// Now, project the data in the new basis.
fullMatrix<double> projected(3,num_vertices);
diff --git a/Geo/STensor3.cpp b/Geo/STensor3.cpp
index 979cf9e562..253339d65a 100644
--- a/Geo/STensor3.cpp
+++ b/Geo/STensor3.cpp
@@ -8,7 +8,6 @@ void SMetric3::print (const char *s) const
(*this)(0,1),(*this)(0,2),(*this)(1,2));
}
-
SMetric3 intersection (const SMetric3 &m1, const SMetric3 &m2)
{
SMetric3 im1 = m1.invert();
@@ -77,4 +76,3 @@ SMetric3 interpolation (const SMetric3 &m1,
im1 += im4;
return im1.invert();
}
-
diff --git a/Geo/STensor3.h b/Geo/STensor3.h
index af6c16471b..8ef6c06d33 100644
--- a/Geo/STensor3.h
+++ b/Geo/STensor3.h
@@ -17,49 +17,35 @@ class SMetric3 {
// 00 10 11 20 21 22
double _val[6];
public:
- inline int getIndex(int i, int j) const{
+ inline int getIndex(int i, int j) const
+ {
static int _index[3][3] = {{0,1,3},{1,2,4},{3,4,5}};
return _index[i][j];
}
- void getMat (fullMatrix<double> & mat) const{
- for (int i=0;i<3;i++){
- for (int j=0;j<3;j++){
- mat(i,j) = _val[getIndex(i,j)];
+ void getMat(fullMatrix<double> &mat) const
+ {
+ for (int i = 0; i < 3; i++){
+ for (int j = 0; j < 3; j++){
+ mat(i,j) = _val[getIndex(i, j)];
}
}
}
- void setMat (const fullMatrix<double> & mat){
- for (int i=0;i<3;i++)
- for (int j=i;j<3;j++)
- _val[getIndex(i,j)] = mat(i,j);
+ void setMat (const fullMatrix<double> & mat)
+ {
+ for (int i = 0; i < 3; i++)
+ for (int j = i; j < 3; j++)
+ _val[getIndex(i, j)] = mat(i, j);
}
- SMetric3(const SMetric3& m){
- for (int i=0; i<6; i++) _val[i]=m._val[i];
+ SMetric3(const SMetric3& m)
+ {
+ for (int i = 0; i < 6; i++) _val[i] = m._val[i];
}
// default constructor, identity
- SMetric3(const double v = 1.0){
+ SMetric3(const double v = 1.0)
+ {
_val[0] = _val[2] = _val[5] = v;
_val[1] = _val[3] = _val[4] = 0.0;
}
-/* SMetric3(const double l1, // (h_1)^-2 */
-/* const double l2, */
-/* const double l3, */
-/* const SVector3 &t1, */
-/* const SVector3 &t2, */
-/* const SVector3 &t3){ */
-/* SVector3 t1b = t1; */
-/* SVector3 t2b = t2; */
-/* SVector3 t3b = t3; */
-/* t1b[0] *= l1; t2b[0] *= l1; t3b[0] *= l1; */
-/* t1b[1] *= l2; t2b[1] *= l2; t3b[1] *= l2; */
-/* t1b[2] *= l3; t2b[2] *= l3; t3b[2] *= l3; */
-/* _val[0] = dot (t1b,t1); */
-/* _val[1] = dot (t2b,t1); */
-/* _val[2] = dot (t2b,t2); */
-/* _val[3] = dot (t3b,t1); */
-/* _val[4] = dot (t3b,t2); */
-/* _val[5] = dot (t3b,t3); */
-/* } */
SMetric3(const double l1, // l1 = h1^-2
const double l2,
const double l3,
@@ -97,42 +83,50 @@ class SMetric3 {
_val[4] = tmp(2,0) * e(0,1) + tmp(2,1) * e(1,1) + tmp(2,2) * e(2,1);
_val[5] = tmp(2,0) * e(0,2) + tmp(2,1) * e(1,2) + tmp(2,2) * e(2,2);
}
- inline double &operator()(int i, int j){
- return _val[getIndex(i,j)];
+ inline double &operator()(int i, int j)
+ {
+ return _val[getIndex(i, j)];
}
- inline double operator()(int i, int j) const{
- return _val[getIndex(i,j)];
+ inline double operator()(int i, int j) const
+ {
+ return _val[getIndex(i, j)];
}
- SMetric3 invert () const {
- fullMatrix<double> m(3,3);
+ SMetric3 invert () const
+ {
+ fullMatrix<double> m(3, 3);
getMat(m);
m.invertInPlace();
SMetric3 ithis;
ithis.setMat(m);
return ithis;
}
- SMetric3 operator + (const SMetric3 &other) const {
+ SMetric3 operator + (const SMetric3 &other) const
+ {
SMetric3 res(*this);
- for (int i=0;i<6;i++)res._val[i]+=other._val[i];
+ for (int i = 0; i < 6; i++) res._val[i] += other._val[i];
return res;
}
- SMetric3& operator += (const SMetric3 &other) {
- for (int i=0;i<6;i++)_val[i]+=other._val[i];
+ SMetric3& operator += (const SMetric3 &other)
+ {
+ for (int i = 0; i < 6; i++) _val[i] += other._val[i];
return *this;
}
- SMetric3& operator *= (const double &other) {
- for (int i=0;i<6;i++)_val[i]*=other;
+ SMetric3& operator *= (const double &other)
+ {
+ for (int i = 0; i < 6; i++) _val[i] *= other;
return *this;
}
- SMetric3& operator *= (const SMetric3 &other) {
- fullMatrix<double> m1(3,3),m2(3,3),m3(3,3);
+ SMetric3& operator *= (const SMetric3 &other)
+ {
+ fullMatrix<double> m1(3, 3), m2(3, 3), m3(3, 3);
getMat(m1);
other.getMat(m2);
- m1.mult(m2,m3);
+ m1.mult(m2, m3);
setMat(m3);
return *this;
}
- SMetric3 transform (fullMatrix<double> &V){
+ SMetric3 transform(fullMatrix<double> &V)
+ {
fullMatrix<double> m(3,3);
getMat(m);
fullMatrix<double> result(3,3),temp(3,3);
@@ -142,25 +136,28 @@ class SMetric3 {
return a;
}
// s: true if eigenvalues are sorted (from min to max of the REAL part)
- void eig (fullMatrix<double> &V, fullVector<double> &S, bool s=false) const {
- fullMatrix<double> me(3,3),right(3,3);
+ void eig(fullMatrix<double> &V, fullVector<double> &S, bool s=false) const
+ {
+ fullMatrix<double> me(3, 3), right(3, 3);
fullVector<double> im(3);
getMat(me);
- me.eig(V,S,im,right,s);
+ me.eig(S, im, V, right, s);
}
void print(const char *) const;
};
// scalar product with respect to the metric
inline double dot(const SVector3 &a, const SMetric3 &m, const SVector3 &b)
-{ return
+{
+ return
b.x() * ( m(0,0) * a.x() + m(1,0) * a.y() + m(2,0) * a.z() ) +
b.y() * ( m(0,1) * a.x() + m(1,1) * a.y() + m(2,1) * a.z() ) +
- b.z() * ( m(0,2) * a.x() + m(1,2) * a.y() + m(2,2) * a.z() ) ;}
+ b.z() * ( m(0,2) * a.x() + m(1,2) * a.y() + m(2,2) * a.z() ) ;
+}
// compute the largest inscribed ellipsoid...
SMetric3 intersection (const SMetric3 &m1,
- const SMetric3 &m2);
+ const SMetric3 &m2);
SMetric3 interpolation (const SMetric3 &m1,
const SMetric3 &m2,
const double t);
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index cc2af1b450..f8e1790908 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -118,7 +118,7 @@ static bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N,
if (M == 1)
DU(0) = R(0) / J(0, 0);
else
- J.lu_solve(R, DU);
+ J.luSolve(R, DU);
for (int i = 0; i < M; i++){
u[i+1] -= underRelax*DU(i);
@@ -204,7 +204,7 @@ static bool computeEquidistantParameters(GFace *gf, double u0, double uN,
if (M == 1)
DU(0) = R(0) / J(0, 0);
else
- J.lu_solve(R, DU);
+ J.luSolve(R, DU);
for (int i = 0; i < M; i++){
t[i + 1] -= DU(i);
diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index 79c5cf0578..d72fd166b0 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -566,7 +566,7 @@ double highOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
fullVector<double> D3(n3);
fullVector<double> R2(n2);
fullMatrix<double> J23K33(n2, n3);
- K33.set_all(0.0);
+ K33.setAll(0.0);
SElement se(e);
El.elementMatrix(&se, K33);
computeMetricVector(gf, e, J32, J23, D3);
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index e87dd53eaf..bd8066edd2 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -6,7 +6,6 @@
set(SRC
Numeric.cpp
fullMatrix.cpp
- EigSolve.cpp
functionSpace.cpp
GaussQuadratureLin.cpp
GaussQuadratureTri.cpp
diff --git a/Numeric/EigSolve.cpp b/Numeric/EigSolve.cpp
deleted file mode 100644
index 0b51712afb..0000000000
--- a/Numeric/EigSolve.cpp
+++ /dev/null
@@ -1,871 +0,0 @@
-// Gmsh - Copyright (C) 1997-2009 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):
-// Laurent Stainier
-
-#include <math.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-
-#define PREC 1.e-16
-
-// Solve eigenproblem for a real general matrix (based on a C++
-// reimplementation of lapack routines from Laurent Stainier)
-//
-// IMPORTANT! How to use the output of EigSolve():
-//
-// 1) If the i-th eigenvalue is real, the i-th COLUMN of the
-// eigenvector Matrix contains the corresponding eigenvector.
-//
-// 2) If the i-th eigenvalue is complex with positive imaginary part,
-// COLUMNS i and (i + 1) contain the real and imaginary parts of the
-// corresponding eigenvector. The conjugate of this vector is the
-// eigenvector for the conjugate eigenvalue.
-//
-// Note the return value. If it does not equal 1, some eigenvalues and
-// all eigenvectors are meaningless.
-
-static void swap(double *a,int inca,double *b,int incb,int n)
-{
- double tmp;
- for (int i=0; i < n; i++, a+=inca, b+=incb) {
- tmp = (*a);
- (*a) = (*b);
- (*b) = tmp;
- }
-}
-
-static void cdiv(double ar,double ai,double br,double bi,double& cr,double& ci)
-{
- double s = fabs(br)+fabs(bi);
- double si = 1.e0/s;
- double ars,ais,brs,bis;
- ars = ar*si;
- ais = ai*si;
- brs = br*si;
- bis = bi*si;
- s = brs*brs+bis*bis;
- cr = (ars*brs+ais*bis)/s;
- ci = (ais*brs-ars*bis)/s;
-}
-
-// Balancing operations
-
-static void balanc(int nm,int n,double *A,int& low,int& high,double *scale)
-{
- static const double RADIX = 16.e0;
-
- double b2 = RADIX*RADIX;
- int i,j,ij,ji,k=0,l=n-1;
- bool last = false;
-
- // search for rows isolating an eigenvalue and push them down
- ROWS:
- // check rows
- for (i=l; i >= 0; i--) {
- for (j=0, ij=i*nm; j <= l; j++, ij++) {
- if (i == j) continue;
- if (A[ij] != 0.e0) break;
- }
- if (j > l) {
- scale[l] = i;
- if (i != l) { // exchange columns and rows
- swap(A+i,nm,A+l,nm,l+1);
- swap(A+i*nm+k,1,A+l*nm+k,1,n-k);
- }
- if (l == 0) goto END;
- l--;
- goto ROWS;
- }
- }
-
- // search for columns isolating an eigenvalue and push them left
- COLUMNS:
- // check column
- for (j=k; j <= l; j++) {
- for (i=k, ij=k*nm+j; i <= l; i++, ij+=nm) {
- if (i == j) continue;
- if (A[ij] != 0.e0) break;
- }
- if (i > l) { // exchange columns and rows
- scale[k] = j;
- if (j != k) {
- swap(A+j,nm,A+k,nm,l+1);
- swap(A+j*nm+k,1,A+k*nm+k,1,n-k);
- }
- k++;
- goto COLUMNS;
- }
- }
-
- // balance the submatrix in rows/columns k to l
- for (i=k; i <= l; i++) scale[i] = 1.e0;
-
- // iterative loop for norm reduction
- while (!last) {
- last = true;
-
- for (i=k; i <= l; i++) {
-
- // calculate row and column norm
- double c=0.e0,r=0.e0;
- for (j=k, ij=i*nm+k, ji=k*nm+i; j <=l ; j++, ij++, ji+=nm) {
- if (i == j) continue;
- c += fabs(A[ji]);
- r += fabs(A[ij]);
- }
-
- // if both are non-zero
- double g,f,s;
- if ((c > PREC) && (r > PREC)) {
-
- g = r/RADIX;
- f = 1.e0;
- s = c+r;
- while (c < g) {
- f *= RADIX;
- c *= b2;
- }
-
- g = r*RADIX;
- while (c > g) {
- f /= RADIX;
- c /= b2;
- }
-
- // balance
- if ((c+r)/f < 0.95*s) {
- g = 1.e0/f;
- scale[i] *= f;
- last = false;
-
- // apply similarity transformation
- for (j=k, ij=i*nm+k; j < n; j++, ij++) A[ij] *= g;
- for (j=0, ji=i ; j <= l; j++, ji+=nm) A[ji] *= f;
- }
- }
- }
- }
-
- END:
- low = k;
- high = l;
-}
-
-static void balbak(int nm,int n,int low,int high,double scale[],int m,double *v)
-{
- int i,j,k,ii,ij;
-
- // quick return
- if (m == 0) return;
-
- // step 1
- if (low != high) {
- for (i=low; i <= high; i++) {
- double s = scale[i];
- for (j=0,ij=i; j < m; j++, ij+=nm) v[ij] *= s;
- }
- }
-
- // step 2
- for (ii=0; ii < n; ii++) {
- if (ii >= low && ii <= high) continue;
- if (ii < low)
- i = low-ii;
- else
- i = ii;
-
- k = (int) scale[i];
- if (k == i) continue;
- swap(v+i,nm,v+k,nm,m);
- }
-}
-
-// Transformation of a general matrix to upper Hessenberg form
-
-static void elmhes(int nm,int n,int low,int high,double *A,int work[])
-{
- int i,j,k,l,m,ij,ji,im1,jm1,jm,mj;
-
- k = low+1;
- l = high-1;
- if (l < k) return;
-
- for (m=k; m <= l; m++) {
- int m1 = m-1;
- double x = 0.e0;
-
- // find the pivot
- i = m;
- for (j=m, jm1=m*nm+m-1; j <= high; j++, jm1+=nm) {
- if (fabs(A[jm1]) < fabs(x)) continue;
- x = A[jm1];
- i = j;
- }
- work[m] = i;
-
- // interchange rows and columns
- if (i != m) {
- swap(A+i*nm+m1,1,A+m*nm+m1,1,n-m1);
- swap(A+i,nm,A+m,nm,high+1);
- }
-
- // carry out elimination
- if (fabs(x) > PREC) {
- for (i=m+1, im1=(m+1)*nm+m-1; i <= high; i++, im1+=nm) {
- double y = A[im1];
- if (y != 0.e0) {
- y /= x;
- A[im1] = y;
- for (j=m, ij=i*nm+m, mj=m*nm+m; j < n; j++, ij++, mj++) A[ij] -= y*A[mj];
- for (j=0, ji=i, jm=m; j <= high; j++, ji+=nm, jm+=nm) A[jm] += y*A[ji];
- }
- }
- }
- }
-}
-
-static void eltran(int nm,int n,int low,int high,double *A,int work[],double *v)
-{
- int i,j,m,ij,im,im1,mj;
-
- // initialize v to identity matrix
- for (j=0; j < n; j++) {
- for (i=0, ij=j*nm; i < n; i++, ij++) v[ij] = 0.e0;
- v[j+j*nm] = 1.e0;
- }
-
- // loop from high-1 to low+1
- for (m=high-1; m > low; m--) {
- int m1 = m+1;
- for (i=m1, im=m1*nm+m, im1=m1+m*nm; i <= high; i++, im+=nm, im1++)
- v[im1] = A[im-1];
-
- i = work[m];
- if (i == m) continue;
-
- for (j=m, mj=m+m*nm, ij=i+m*nm; j <= high; j++, mj+=nm, ij+=nm) {
- v[mj] = v[ij];
- v[ij] = 0.e0;
- }
-
- v[i+m*nm] = 1.e0;
- }
-}
-
-// Find eigenvalues of a real upper Hessenberg matrix by the QR method
-
-static int hqr(int nm,int n,int low,int high,double *H,double wr[],double wi[],double *vv)
-{
- int i,j,k,l,m,ij,ij1,ik,ik1,ik2,in,in1,jn,jn1,kj,kj1,kj2,nj,nj1,mmin;
- double p=0.0e0,q=0.0e0,r=0.0e0,s=0.0e0,x,y,z=0.0e0,u,v,w;
-
- // store roots isolated by balanc, and compute matrix norm
- k=0;
- double norm = 0.e0;
- for (i=0; i < n; i++) {
- for (j=k, ij=i*nm+k; j < n; j++, ij++) norm += fabs(H[ij]);
-
- k = i;
- if (i >= low && i <= high) continue;
- wr[i] = H[i*nm+i];
- wi[i] = 0.e0;
- }
-
- // search for next eigenvalues
- int nn=high,nn1;
- double t=0.e0;
- while (nn >= low) {
-
- int iter=0;
- do {
-
- // look for single small subdiagonal element
- for (l=nn; l > low; l--) {
- s = fabs(H[(l-1)*(nm+1)])+fabs(H[l*(nm+1)]);
- if (s == 0.e0) s = norm;
- if ((fabs(H[l*nm+l-1])+s) == s) break;
- }
-
- // form shift
- x = H[nn*nm+nn];
- if (l == nn) { // ... one root found
- wr[nn] = x+t;
- wi[nn] = 0.e0;
- H[nn*nm+nn] = wr[nn];
- nn--;
- }
- else {
- nn1 = nn-1;
- y = H[nn1*nm+nn1];
- w = H[nn*nm+nn1]*H[nn1*nm+nn];
- if (l == nn1) { // ... two roots found
- p = 0.5*(y-x);
- q = p*p+w;
- z = sqrt(fabs(q));
- H[nn*nm+nn] = x+t;
- x += t;
- H[nn1*nm+nn1] = y+t;
- if (q >= 0.e0) { // ... a real pair
- z = p+((p > 0.e0)?(z):(-z));
- wr[nn1] = wr[nn] = x+z;
- if (fabs(z) > PREC) wr[nn] = x-w/z;
- wi[nn1] = wi[nn] = 0.e0;
-
- x = H[nn*nm+nn1];
- s = fabs(x)+fabs(z);
- p = x/s;
- q = z/s;
- r = sqrt(p*p+q*q);
- p /= r;
- q /= r;
- // row modification
- for (j=nn1, nj=nn*nm+nn1, nj1=nn1*nm+nn1; j < n; j++, nj++, nj1++) {
- z = H[nj1];
- H[nj1] = q*z+p*H[nj];
- H[nj] = q*H[nj]-p*z;
- }
- // column modification
- for (i=0, in=nn, in1=nn1; i <= nn; i++, in+=nm, in1+=nm) {
- z = H[in1];
- H[in1] = q*z+p*H[in];
- H[in] = q*H[in]-p*z;
- }
- // accumulate transformations
- for (i=low, in=low+nn*nm, in1=low+nn1*nm; i <= high; i++, in++, in1++) {
- z = vv[in1];
- vv[in1] = q*z+p*vv[in];
- vv[in] = q*vv[in]-p*z;
- }
- }
- else { // ... a complex pair
- wr[nn1] = wr[nn] = x+p;
- wi[nn1] = z;
- wi[nn] = -z;
- }
- nn -= 2;
- }
- else { // ... no roots found. Continue iterations.
- if (iter == 30) return nn+1;
- if (iter == 10 || iter == 20) { // form exceptional shift
- t += x;
- for (i=low; i <=nn; i++) H[i*nm+i] -= x;
- s = fabs(H[nn*nm+nn1])+fabs(H[nn1*nm+nn1-1]);
- y = x = 0.75*s;
- w = -0.4375*s*s;
- }
- iter++;
-
- // look for two consecutive small sub-diagonal elements
- for (m=nn-2; m >= l; m--) {
- z = H[m*nm+m];
- r = x-z;
- s = y-z;
- p = (r*s-w)/H[(m+1)*nm+m]+H[m*nm+m+1];
- q = H[(m+1)*(nm+1)]-z-r-s;
- r = H[(m+2)*nm+m+1];
- s = fabs(p)+fabs(q)+fabs(r);
- p /= s;
- q /= s;
- r /= s;
- if (m == l) break;
- u = fabs(p)*(fabs(H[(m-1)*(nm+1)])
- +fabs(z)+fabs(H[(m+1)*(nm+1)]));
- v = u+fabs(H[m*nm+m-1])*(fabs(q)+fabs(r));
- if (u == v) break;
- }
- for (i=m+2; i <= nn; i++) {
- H[i*nm+i-2] = 0.e0;
- if (i != (m+2)) H[i*nm+i-3] = 0.e0;
- }
-
- // double QR step on rows l to nn and columns m to nn
- for (k=m; k <= nn1; k++) {
-
- if (k != m) { // begin setup of Householder vector
- p = H[k*nm+k-1];
- q = H[(k+1)*nm+k-1];
- r = 0.e0;
- if (k != nn1) r = H[(k+2)*nm+k-1];
- x = fabs(p)+fabs(q)+fabs(r);
- if (x < PREC) continue;
- // scale to prevent underflow or overflow
- p /= x;
- q /= x;
- r /= x;
- }
-
- s = sqrt(p*p+q*q+r*r);
- if (p < 0.e0) s = -s;
- if (k == m) {
- if (l != m) H[k*nm+k-1] = -H[k*nm+k-1];
- }
- else
- H[k*nm+k-1] = -s*x;
- p += s;
- x = p/s;
- y = q/s;
- z = r/s;
- q /= p;
- r /= p;
-
- // row modification
- kj = k*nm+k; kj1 = kj+nm; kj2 = kj1+nm;
- for (j=k; j <= nn; j++, kj++, kj1++,kj2++) {
- p = H[kj]+q*H[kj1];
- if (k != nn1) {
- p += r*H[kj2];
- H[kj2] -= p*z;
- }
- H[kj1] -= p*y;
- H[kj] -= p*x;
- }
-
- // column modification
- mmin = (nn < k+3) ? nn : k+3;
- ik = l*nm+k; ik1 = ik+1; ik2 = ik1+1;
- for (i=l; i <= mmin; i++, ik+=nm, ik1+=nm, ik2+=nm) {
- p = x*H[ik]+y*H[ik1];
- if (k != nn1) {
- p += z*H[ik2];
- H[ik2] -= p*r;
- }
- H[ik1] -= p*q;
- H[ik] -= p;
- }
-
- // accumulate transformations
- ik = low+k*nm; ik1 = ik+nm; ik2 = ik1+nm;
- for (i=low; i <= high; i++, ik++, ik1++, ik2++) {
- p = x*vv[ik]+y*vv[ik1];
- if (k != nn1) {
- p += z*vv[ik2];
- vv[ik2] -= p*r;
- }
- vv[ik1] -= p*q;
- vv[ik] -= p;
- }
- }
- }
- }
-
- } while (l < nn-1);
- }
-
- // Backsubstitute to find vectors of upper triangular form
- for (nn=n-1; nn >= 0; nn--) {
- p = wr[nn];
- q = wi[nn];
- nn1 = nn-1;
- if (q == 0.e0) { // real vector
- m = nn;
- H[nn*nm+nn] = 1.e0;
- for (i=nn1; i >= 0; i--) {
- in = i*nm+nn;
- in1 = i*nm+nn1;
- w = H[i*nm+i]-p;
-
- r = 0.e0;
- for (j=m, ij=i*nm+m, jn=m*nm+nn; j <= nn; j++, ij++, jn+=nm)
- r += H[ij]*H[jn];
-
- if (wi[i] < 0.e0) {
- z = w;
- s = r;
- continue;
- }
-
- m = i;
- if (wi[i] == 0.e0) {
- t = w;
- if (t == 0.e0) {
- u = norm;
- t = u;
- do {
- t *= 0.01;
- v = norm+t;
- } while (v > u);
- }
- H[in] = -r/t;
- }
- else { // solve real equation
- x = H[i*nm+i+1];
- y = H[(i+1)*nm+i];
- q = (wr[i]-p)*(wr[i]-p)+wi[i]*wi[i];
- t = (x*s-z*r)/q;
- H[i*nm+nn] = t;
- if (fabs(x) <= fabs(z))
- H[in+nm] = (-s-y*t)/z;
- else
- H[in+nm] = (-r-w*t)/x;
- }
-
- // overflow control
- t = fabs(H[in]);
- if (t != 0.e0) {
- u = t;
- v = t+1.e0/t;
- if (u >= v)
- for (j=i, jn=in; j <= nn; j++, jn+=nm) H[jn] /= t;
- }
- }
- }
- else if (q < 0.e0) { // complex vector
- m = nn1;
-
- // last vector component chosen imaginary, so that eigenvector matrix is triangular
- if (fabs(H[nn*nm+nn1]) <= fabs(H[nn1*nm+nn]))
- cdiv(0.e0,-H[nn1*nm+nn],H[nn1*nm+nn1]-p,q,H[nn1*nm+nn1],H[nn1*nm+nn]);
- else {
- H[nn1*nm+nn1] = q/H[nn*nm+nn1];
- H[nn1*nm+nn] = -(H[nn*nm+nn]-p)/H[nn*nm+nn1];
- }
-
- for (i=nn-2; i >= 0; i--) {
- in = i*nm+nn;
- in1 = i*nm+nn1;
- w = H[i*nm+i]-p;
-
- double ra = 0.e0;
- double sa = 0.e0;
- for (j=m, ij=i*nm+m, jn=m*nm+nn, jn1=m*nm+nn1; j <= nn;
- j++, ij++, jn+=nm, jn1+=nm) {
- ra += H[ij]*H[jn1];
- sa += H[ij]*H[jn];
- }
-
- if (wi[i] < 0.e0) {
- z = w;
- r = ra;
- s = sa;
- continue;
- }
-
- m = i;
- if (wi[i] == 0.e0) {
- cdiv(-ra,-sa,w,q,H[in1],H[in]);
- }
- else { // solve complex equations
- x = H[i*nm+i+1];
- y = H[(i+1)*nm+i];
- double vr = (wr[i]-p)*(wr[i]-p)+wi[i]*wi[i]-q*q;
- double vi = 2*(wr[i]-p)*q;
- if (vr == 0.e0 && vi == 0.e0) {
- u = norm*(fabs(w)+fabs(q)+fabs(x)
- +fabs(y)+fabs(z));
- vr = u;
- do {
- vr *= 0.01;
- v = u+vr;
- } while (v > u);
- }
- cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi,H[in1],H[in]);
- if (fabs(x) <= (fabs(z)+fabs(q)))
- cdiv(-r-y*H[in1],-s-y*H[in],z,q,
- H[in1+nm],H[in+nm]);
- else {
- H[in1+nm] = (-ra-w*H[in1]+q*H[in])/x;
- H[in+nm] = (-sa-w*H[in]-q*H[in1])/x;
- }
- }
-
- // overflow control
- u = fabs(H[in1]); v = fabs(H[in]);
- t = (u > v) ? u:v;
- if (t != 0.e0) {
- u = t;
- v = t+1.e0/t;
- if (u >= v) {
- for (j=i, jn=in, jn1=in1; j < nn; j++, jn+=nm, jn1+=nm) {
- H[jn1] /= t;
- H[jn] /= t;
- }
- }
- }
- }
- }
- }
-
- // vectors of isolated roots
- for (i=0; i < n; i++) {
- if (i >= low && i <= high) continue;
- for (j=0, ij=i, ij1=i*nm; j < n; j++, ij+=nm, ij1++)
- vv[ij] = H[ij1];
- }
-
- // multiply by transformation matrix to give evctors of original full matrix
- for (j=n-1; j >= 0; j--) {
- m = (j < high) ? j:high;
- for (i=low, ij=low+j*nm; i <= high; i++, ij++) {
- z = 0.e0;
- for (k=low, ik=i+low*nm, kj=low*nm+j; k <= m; k++, ik+=nm, kj+=nm)
- z += vv[ik]*H[kj];
- vv[ij] = z;
- }
- }
-
- return 0;
-}
-
-// Normalizes the eigenvectors. This subroutine is based on the
-// LINPACK routine SNRM2, written 25 October 1982, modified on 14
-// October 1993 by Sven Hammarling of Nag Ltd.
-
-static void normvec(int n, double *Z, double *wi)
-{
- for (int j = 0; j < n; j++){ // Go through the columns of the array
- double scale = 0.0;
- double ssq = 1.0;
-
- for (int i = 0; i < n; i++){
- double absxi = fabs(Z[j*n+i]);
- if (absxi > PREC){
- double dummy = scale/absxi;
- if (scale < absxi){
- ssq = 1.0 + ssq*dummy*dummy;
- scale = absxi;
- }
- else
- ssq += 1.0/dummy/dummy;
- }
- }
-
- if (fabs(wi[j]) > PREC){
- // If complex eigenvalue, take into account imaginary part of
- // eigenvector
- for (int i = 0; i < n; i++){
- double absxi = fabs(Z[(j + 1)*n+i]);
- if (absxi > PREC){
- double dummy = scale/absxi;
- if (scale < absxi){
- ssq = 1.0 + ssq*dummy*dummy;
- scale = absxi;
- }
- else
- ssq += 1.0/dummy/dummy;
- }
- }
- }
-
- // this is norm of the (possibly complex) vector
- double norm = scale*sqrt(ssq);
-
- for (int i = 0; i < n; i++)
- Z[j*n+i] /= norm;
-
- if (fabs(wi[j]) > PREC){
- // If complex eigenvalue, also scale imaginary part of
- // eigenvector
- j++;
- for (int i = 0; i < n; i++)
- Z[j*n+i] /= norm;
- }
- }
-}
-
-int EigSolve(int nm,int n,double *A,double *wr,double *wi,
- double *v,int *work1,double *work2)
-{
- int is1,is2,ierr;
-
- balanc(nm,n,A,is1,is2,work2);
- elmhes(nm,n,is1,is2,A,work1);
- eltran(nm,n,is1,is2,A,work1,v);
- ierr = hqr(nm,nm,is1,is2,A,wr,wi,v);
- if (ierr) return 0;
- balbak(nm,n,is1,is2,work2,n,v);
- normvec(n,v,wi);
- return 1;
-}
-
-// Solve eigen problem for a real, symmetric matrix using the Jacobi
-// algorithm (based on a routine from Laurent Stainier)
-
-int EigSolveSym(int n,int nm,double *A,double *d,double *V,
- double *b,double *z)
-{
- static const int NSWMAX = 50;
-
- int i,j,k,ii,ij,ki=0,kj=0,kki,kkj,ival,jval;
- int nrot,nsweep;
- double c,g,h,s,t,tau,theta,tresh,sum;
-
- // initialize eigenvectors
- for (j=0,jval=0; j < n; j++,jval+=nm) {
- for (i=0; i < n; i++) V[jval+i] = 0.e0;
- V[jval+j] = 1.e0;
- }
-
- // initialize b and d to the diagonal of A
- for (i=0, ii=0; i < n; i++, ii+=(i+1)) {
- b[i] = d[i] = A[ii];
- z[i] = 0.e0;
- }
-
- // begin sweeping process
- nrot = 0;
- for (nsweep=0; nsweep < NSWMAX; nsweep++) {
-
- // Sum off-diagonal elements
- sum = 0.e0;
- for (i=0, ij=0; i < n; i++, ij++)
- for (j=0; j < i; j++, ij++)
- sum += fabs(A[ij]);
-
- if (sum < PREC) break;
-
- // compute a treshold on the first 3 sweeps
- if (nsweep < 4)
- tresh = 0.2*sum/(n*n);
- else
- tresh = 0.e0;
-
- // browse the matrix
- for (i=0, ij=0, ival=0; i < n; i++, ij++, ival+=i)
- for (j=0, jval=0; j < i; j++, ij++, jval+=j) {
-
- // test off-diagonal element
- g = 100.*fabs(A[ij]);
- if ((nsweep > 4)
- && (fabs(d[i]+g) == fabs(d[i]))
- && (fabs(d[j]+g) == fabs(d[j]))) {
-
- A[ij] = 0.e0;
- }
- else if (fabs(A[ij]) > tresh) {
-
- // compute the rotation
- h = d[j]-d[i];
- if ((fabs(h)+g) == fabs(h))
- t = A[ij]/h;
- else {
- theta = 0.5*h/A[ij];
- t = 1.e0/(fabs(theta)+sqrt(1.0+theta*theta));
- if (theta < 0.e0) t = -t;
- }
- c = 1.e0/sqrt(1.+t*t);
- s = t*c;
- tau = s/(1.e0+c);
-
- // rotate rows i and j
- h = t*A[ij];
- z[i] -= h;
- z[j] += h;
- d[i] -= h;
- d[j] += h;
- A[ij] = 0.e0;
-
- for (k=0,kki=i*nm,kkj=j*nm; k < n; k++, kki++, kkj++) {
-
- g = V[kki];
- h = V[kkj];
- V[kki] = g-s*(h+g*tau);
- V[kkj] = h+s*(g-h*tau);
-
- ki = (k<=i)?(ival+k):(ki+k);
- kj = (k<=j)?(jval+k):(kj+k);
- if (k == i || k == j) continue;
-
- g = A[ki];
- h = A[kj];
- A[ki] = g-s*(h+g*tau);
- A[kj] = h+s*(g-h*tau);
- }
-
- nrot++;
- }
- }
-
- // update d and reinitialize z
- for (i=0; i < n; i++) {
- b[i] += z[i];
- d[i] = b[i];
- z[i] = 0.e0;
- }
- }
-
- if(nsweep >= NSWMAX)
- return 0;
- else
- return nrot+1;
-}
-
-// Sort the eigenvalues/vectors in ascending order according to their
-// real part. Warning: this will screw up the ordering of complex
-// eigenvectors.
-
-void EigSort(int n, double *wr, double *wi, double *B)
-{
- for (int i = 0; i < n - 1; i++){
- int k = i;
- double ek = wr[i];
- // search for something to swap
- for (int j = i + 1; j < n; j++){
- const double ej = wr[j];
- if(ej < ek){
- k = j;
- ek = ej;
- }
- }
- if (k != i){
- swap(&wr[i], 1, &wr[k], 1, 1);
- swap(&wi[i], 1, &wi[k], 1, 1);
- swap(&B[n*i], 1, &B[n*k], 1, n);
- }
- }
-}
-
-// Small driver for 3x3 matrices, that does not modify the input
-// matrix
-
-int EigSolve3x3(const double A[9], double wr[3], double wi[3], double v[9])
-{
- int ierr;
- if(fabs(A[1]-A[3]) < PREC &&
- fabs(A[2]-A[6]) < PREC &&
- fabs(A[5]-A[7]) < PREC){
- double work1[3], work2[3], S[6] = { A[0],
- A[1], A[4],
- A[2], A[5], A[8]};
- ierr = EigSolveSym(3, 3, S, wr, v, work1, work2);
- wi[0] = wi[1] = wi[2] = 0.0;
- }
- else{
- int work1[3];
- double work2[3], M[9] = { A[0], A[1], A[2],
- A[3], A[4], A[5],
- A[6], A[7], A[8]};
- ierr = EigSolve(3, 3, M, wr, wi, v, work1, work2);
- }
- EigSort(3, wr, wi, v);
- return ierr;
-}
-
-# if 0
-int main ()
-{
- //double A[9] = {-0.00299043,-8.67362e-19,0,
- // -8.67362e-19,-0.00299043,-1.73472e-18,
- // 0,-1.73472e-18,0.01};
- double A[9] = {1, 2, 3, 2, 4, 5, 3, 5, 6};
- //double A[9] = {1, 2, 3, 1, 4, 5, 3, 5, 6};
- double wr[3], wi[3], B[9];
-
- for(int i = 0; i < 1; i++) // perf test
- if(!EigSolve3x3(A, wr, wi, B))
- printf("EigSolve did not converge!\n");
-
- if(wi[0] != 0.0 || wi[1] != 0.0 || wi[2] != 0.0)
- printf("imaginary eigenvalues\n");
-
- printf ("real= %g %g %g \n", wr[0], wr[1], wr[2]);
- printf ("imag= %.16g %.16g %.16g \n", wi[0], wi[1], wi[2]);
- printf ("eigvec1 = %g %g %g \n", B[0], B[1], B[2]);
- printf ("eigvec2 = %g %g %g \n", B[3], B[4], B[5]);
- printf ("eigvec3 = %g %g %g \n", B[6], B[7], B[8]);
-}
-#endif
diff --git a/Numeric/EigSolve.h b/Numeric/EigSolve.h
deleted file mode 100644
index d4eb1aa6fa..0000000000
--- a/Numeric/EigSolve.h
+++ /dev/null
@@ -1,17 +0,0 @@
-// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-
-#ifndef _EIGSOLVE_H_
-#define _EIGSOLVE_H_
-
-int EigSolve(int nm,int n,double *A,double *wr,double *wi,
- double *v,int *work1,double *work2);
-int EigSolveSym(int n,int nm,double *A,double *d,double *V,
- double *b,double *z);
-void EigSort(int n, double *wr, double *wi, double *B);
-
-int EigSolve3x3(const double A[9], double wr[3], double wi[3], double v[9]);
-
-#endif
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index a9007871b0..bcff7a120f 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -552,7 +552,7 @@ bool newton_fd(void (*func)(fullVector<double> &, fullVector<double> &, void *),
if (N == 1)
dx(0) = f(0) / J(0, 0);
else
- J.lu_solve(f, dx);
+ J.luSolve(f, dx);
for (int i = 0; i < N; i++)
x(i) -= relax * dx(i);
diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index 69928a6f35..36dde4c57c 100644
--- a/Numeric/fullMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -105,21 +105,17 @@ void fullMatrix<std::complex<double> >::mult(const fullVector<std::complex<doubl
#if defined(HAVE_LAPACK)
extern "C" {
- void F77NAME(dgesv)(int *N, int *nrhs, double *A, int *lda, int *ipiv,
+ void F77NAME(dgesv)(int *N, int *nrhs, double *A, int *lda, int *ipiv,
double *b, int *ldb, int *info);
void F77NAME(dgetrf)(int *M, int *N, double *A, int *lda, int *ipiv, int *info);
- void F77NAME(dgetri)(int *M, double *A, int *lda, int *ipiv, double *work, int *lwork,
- int *info);
+ void F77NAME(dgetri)(int *M, double *A, int *lda, int *ipiv, double *work,
+ int *lwork, int *info);
void F77NAME(dgesvd)(const char* jobu, const char *jobvt, int *M, int *N,
double *A, int *lda, double *S, double* U, int *ldu,
double *VT, int *ldvt, double *work, int *lwork, int *info);
- void F77NAME(dgeev)(const char *jobvl, const char *jobvr,
- int *n, double *a, int *lda,
- double *wr, double *wi,
- double *vl, int *ldvl,
- double *vr, int *ldvr,
- double *work, int *lwork,
- int *info);
+ void F77NAME(dgeev)(const char *jobvl, const char *jobvr, int *n, double *a,
+ int *lda, double *wr, double *wi, double *vl, int *ldvl,
+ double *vr, int *ldvr, double *work, int *lwork, int *info);
}
template<>
@@ -148,17 +144,15 @@ bool fullMatrix<double>::invertInPlace()
template<>
-bool fullMatrix<double>::eig(fullMatrix<double> &VL, fullVector<double> &DR,
- fullVector<double> &DI, fullMatrix<double> &VR,
+bool fullMatrix<double>::eig(fullVector<double> &DR, fullVector<double> &DI,
+ fullMatrix<double> &VL, fullMatrix<double> &VR,
bool sortRealPart)
{
int N = size1(), info;
- int LWORK = 10 * N;
- double *work = new double[LWORK];
-
+ int lwork = 10 * N;
+ double *work = new double[lwork];
F77NAME(dgeev)("V", "V", &N, _data, &N, DR._data, DI._data,
- VL._data, &N, VR._data, &N, work, &LWORK, &info);
-
+ VL._data, &N, VR._data, &N, work, &lwork, &info);
delete [] work;
if(info > 0)
@@ -166,25 +160,25 @@ bool fullMatrix<double>::eig(fullMatrix<double> &VL, fullVector<double> &DR,
else if(info < 0)
Msg::Error("Wrong %d-th argument in eig", -info);
- if (sortRealPart) {
+ if(sortRealPart) {
double tmp[8];
// do permutations
- for (int i = 0; i < size1() - 1; i++) {
+ for(int i = 0; i < size1() - 1; i++) {
int minR = i;
- for (int j = i + 1; j < size1(); j++)
- if (fabs(DR(j)) < fabs(DR(minR))) minR = j;
- if ( minR != i ){
+ for(int j = i + 1; j < size1(); j++)
+ if(fabs(DR(j)) < fabs(DR(minR))) minR = j;
+ if(minR != i){
tmp[0] = DR(i); tmp[1] = DI(i);
- tmp[2] = VL(0,i); tmp[3] = VL(1,i); tmp[4] = VL(2,i);
- tmp[5] = VR(0,i); tmp[6] = VR(1,i); tmp[7] = VR(2,i);
+ tmp[2] = VL(0, i); tmp[3] = VL(1, i); tmp[4] = VL(2, i);
+ tmp[5] = VR(0, i); tmp[6] = VR(1, i); tmp[7] = VR(2, i);
DR(i) = DR(minR); DI(i) = DI(minR);
- VL(0,i) = VL(0,minR); VL(1,i) = VL(1,minR); VL(2,i) = VL(2,minR);
- VR(0,i) = VR(0,minR); VR(1,i) = VR(1,minR); VR(2,i) = VR(2,minR);
+ VL(0,i) = VL(0, minR); VL(1, i) = VL(1, minR); VL(2, i) = VL(2, minR);
+ VR(0,i) = VR(0, minR); VR(1, i) = VR(1, minR); VR(2, i) = VR(2, minR);
DR(minR) = tmp[0]; DI(minR) = tmp[1];
- VL(0,minR) = tmp[2]; VL(1,minR) = tmp[3]; VL(2,minR) = tmp[4];
- VR(0,minR) = tmp[5]; VR(1,minR) = tmp[6]; VR(2,minR) = tmp[7];
+ VL(0, minR) = tmp[2]; VL(1, minR) = tmp[3]; VL(2, minR) = tmp[4];
+ VR(0, minR) = tmp[5]; VR(1, minR) = tmp[6]; VR(2, minR) = tmp[7];
}
}
}
@@ -193,7 +187,7 @@ bool fullMatrix<double>::eig(fullMatrix<double> &VL, fullVector<double> &DR,
}
template<>
-bool fullMatrix<double>::lu_solve(const fullVector<double> &rhs, fullVector<double> &result)
+bool fullMatrix<double>::luSolve(const fullVector<double> &rhs, fullVector<double> &result)
{
int N = size1(), nrhs = 1, lda = N, ldb = N, info;
int *ipiv = new int[N];
@@ -257,10 +251,10 @@ bool fullMatrix<double>::svd(fullMatrix<double> &V, fullVector<double> &S)
{
fullMatrix<double> 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));
- fullVector<double> WORK(LWORK);
+ int lwork = std::max(3 * std::min(M, N) + std::max(M, N), 5 * std::min(M, N));
+ fullVector<double> WORK(lwork);
F77NAME(dgesvd)("O", "A", &M, &N, _data, &LDA, S._data, _data, &LDA,
- VT._data, &LDVT, WORK._data, &LWORK, &info);
+ VT._data, &LDVT, WORK._data, &lwork, &info);
V = VT.transpose();
if(info == 0) return true;
if(info > 0)
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index a96d7326b7..4447b51283 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -146,7 +146,7 @@ class fullMatrix
}
#endif
;
- inline void set_all(const scalar &m)
+ inline void setAll(const scalar &m)
{
for(int i = 0; i < _r * _c; i++) _data[i] = m;
}
@@ -198,7 +198,7 @@ class fullMatrix
_data[j + _r * i] = t;
}
}
- bool lu_solve(const fullVector<scalar> &rhs, fullVector<scalar> &result)
+ bool luSolve(const fullVector<scalar> &rhs, fullVector<scalar> &result)
#if !defined(HAVE_LAPACK)
{
Msg::Error("LU factorization requires LAPACK");
@@ -214,11 +214,9 @@ class fullMatrix
}
#endif
;
- bool eig(fullMatrix<scalar> &VL, // left eigenvectors
- fullVector<double> &DR, // Real part of eigenvalues
- fullVector<double> &DI, // Im part of eigen
- fullMatrix<scalar> &VR,
- bool sortRealPart=false) // if true: sorted from min 'DR' to max 'DR'
+ bool eig(fullVector<double> &eigenValReal, fullVector<double> &eigenValImag,
+ fullMatrix<scalar> &leftEigenVect, fullMatrix<scalar> &rightEigenVect,
+ bool sortRealPart=false)
#if !defined(HAVE_LAPACK)
{
Msg::Error("Eigenvalue computations requires LAPACK");
diff --git a/Numeric/functionSpace.cpp b/Numeric/functionSpace.cpp
index d54c96c461..c535a49c53 100644
--- a/Numeric/functionSpace.cpp
+++ b/Numeric/functionSpace.cpp
@@ -191,7 +191,7 @@ static fullMatrix<double> generatePascalSerendipityTetrahedron(int order)
4 * std::max(0, (order - 2) * (order - 1) / 2);
fullMatrix<double> monomials(nbMonomials, 3);
- monomials.set_all(0);
+ monomials.setAll(0);
int index = 1;
for (int p = 1; p < order; p++) {
for (int i = 0; i < 3; i++) {
diff --git a/Plugin/Eigenvectors.cpp b/Plugin/Eigenvectors.cpp
index 14cda42fc6..338aac5612 100644
--- a/Plugin/Eigenvectors.cpp
+++ b/Plugin/Eigenvectors.cpp
@@ -5,7 +5,7 @@
#include "Eigenvectors.h"
#include "Numeric.h"
-#include "EigSolve.h"
+#include "fullMatrix.h"
#include "GmshDefines.h"
StringXNumber EigenvectorsOptions_Number[] = {
@@ -91,7 +91,8 @@ PView *GMSH_EigenvectorsPlugin::execute(PView *v)
PViewDataList *dmax = getDataList(max);
int nbcomplex = 0;
-
+ fullMatrix<double> mat(3, 3), vl(3, 3), vr(3, 3);
+ fullVector<double> dr(3), di(3);
for(int ent = 0; ent < data1->getNumEntities(0); ent++){
for(int ele = 0; ele < data1->getNumElements(0, ent); ele++){
if(data1->skipElement(0, ent, ele)) continue;
@@ -115,21 +116,21 @@ PView *GMSH_EigenvectorsPlugin::execute(PView *v)
}
for(int step = 0; step < data1->getNumTimeSteps(); step++){
for(int nod = 0; nod < numNodes; nod++){
- double val[9];
- for(int comp = 0; comp < numComp; comp++)
- data1->getValue(step, ent, ele, nod, comp, val[comp]);
- double wr[3], wi[3], B[9];
- if(!EigSolve3x3(val, wr, wi, B))
- Msg::Error("Eigensolver failed to converge");
- nbcomplex += nonzero(wi);
- if(!scale) wr[0] = wr[1] = wr[2] = 1.;
- for(int i = 0; i < 3; i++){
- double res;
- // wrong if there are complex eigenvals (B contains both
- // real and imag parts: cf. explanation in EigSolve.cpp)
- res = wr[0] * B[i]; outmin->push_back(res);
- res = wr[1] * B[3 + i]; outmid->push_back(res);
- res = wr[2] * B[6 + i]; outmax->push_back(res);
+ for(int i = 0; i < 3; i++)
+ for(int j = 0; j < 3; j++)
+ data1->getValue(step, ent, ele, nod, 3 * i + j, mat(i, j));
+ if(mat.eig(dr, di, vl, vr, true)){
+ if(!scale) dr(0) = dr(1) = dr(2) = 1.;
+ for(int i = 0; i < 3; i++){
+ double res;
+ res = dr(0) * vr(i, 0); outmin->push_back(res);
+ res = dr(1) * vr(i, 1); outmid->push_back(res);
+ res = dr(2) * vr(i, 2); outmax->push_back(res);
+ }
+ if(di(0) || di(1) || di(2)) nbcomplex++;
+ }
+ else{
+ Msg::Error("Could not compute eigenvalues/vectors");
}
}
}
@@ -137,8 +138,7 @@ PView *GMSH_EigenvectorsPlugin::execute(PView *v)
}
if(nbcomplex)
- Msg::Error("%d tensors have complex eigenvalues/eigenvectors",
- nbcomplex);
+ Msg::Error("%d tensors have complex eigenvalues", nbcomplex);
for(int i = 0; i < data1->getNumTimeSteps(); i++){
double time = data1->getTime(i);
diff --git a/Solver/convexCombinationTerm.h b/Solver/convexCombinationTerm.h
index dcaef8fb3b..6a2088756c 100644
--- a/Solver/convexCombinationTerm.h
+++ b/Solver/convexCombinationTerm.h
@@ -39,7 +39,7 @@ class convexCombinationTerm : public femTerm<double,double> {
{
MElement *e = se->getMeshElement();
- m.set_all(0.);
+ m.setAll(0.);
const int nbNodes = e->getNumVertices();
const double _diff = 1.0;
for (int j = 0; j < nbNodes; j++){
diff --git a/Solver/crossConfTerm.h b/Solver/crossConfTerm.h
index 2b0488321c..52521de117 100644
--- a/Solver/crossConfTerm.h
+++ b/Solver/crossConfTerm.h
@@ -55,7 +55,7 @@ class crossConfTerm : public femTerm<double, double> {
double grads[256][3];
e->getIntegrationPoints(integrationOrder, &npts, &GP);
- m.set_all(0.);
+ m.setAll(0.);
for (int i = 0; i < npts; i++){
const double u = GP[i].pt[0];
diff --git a/Solver/eigenSolve.cpp b/Solver/eigenSolve.cpp
new file mode 100644
index 0000000000..2a61b17162
--- /dev/null
+++ b/Solver/eigenSolve.cpp
@@ -0,0 +1,104 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include "eigenSolve.h"
+
+#if defined(HAVE_SLEPC)
+
+#include <slepceps.h>
+
+eigenSolve::eigenSolve(dofManager *manager, const char *A, const char *B)
+{
+ if(A){
+ _A = manager->getLinearSystem(A);
+ if(!_A) Msg::Error("Could not find system A");
+ }
+ if(B){
+ _B = manager->getLinearSystem(B);
+ if(!_B) Msg::Error("Could not find system B");
+ }
+}
+
+void eigenSolve::solve()
+{
+ if(!_A) return;
+ Mat A = _A->getMatrix();
+ Mat B = _B ? _b->getMatrix() : PETSC_NULL;
+
+ // treatment of a generalized eigenvalue problem A x - \lambda B x = 0
+ EPS eps;
+ _try(EPSCreate(PETSC_COMM_WORLD, &eps));
+ _try(EPSSetOperators(eps, A, B));
+ _try(EPSSetProblemType(eps, EPS_GNHEP));
+
+ // set solver parameters at runtime
+ _try(EPSSetFromOptions(eps));
+
+ Msg::Info("SLEPc solving...");
+ _try(EPSSolve(eps));
+ int its;
+ _try(EPSGetIterationNumber(eps, &its));
+ Msg::Info("Number of iterations of the method: %d", its);
+
+ // get some information from the solver and display it
+ EPSType type;
+ _try(EPSGetType(eps, &type));
+ Msg::Info("Solution method: %s", type);
+ int nev;
+ _try(EPSGetDimensions(eps, &nev, PETSC_NULL));
+ Msg::Info("Number of requested eigenvalues: %d", nev);
+ PetscReal tol;
+ int maxit;
+ _try(EPSGetTolerances(eps, &tol, &maxit));
+ Msg::Info("Stopping condition: tol=%.4g, maxit=%d", tol, maxit);
+
+ // get number of converged approximate eigenpairs
+ int nconv;
+ _try(EPSGetConverged(eps, &nconv));
+ Msg::Info("Number of converged eigenpairs: %d", nconv);
+
+ Msg::Info("Re[Eigenvalue] Im[Eigenvalue] Rel. error = ||Ax - eig*Bx||/||eig*x||");
+
+ // if number of converged solutions is greated than requested
+ // discarded surplus solutions
+ if(nconv > nev) nconv = nev;
+ // if number of converged solutions is less than requested display
+ // only converged solutions
+ if(nconv < nev) nev = nconv;
+
+ Vec xr, xi;
+ _try(MatGetVecs(A, PETSC_NULL, &xr));
+ _try(MatGetVecs(A, PETSC_NULL, &xi));
+ for (int i = nconv - 1; i > nconv - nev - 1; i--){
+ PetscScalar kr, ki;
+ _try(EPSGetEigenpair(eps, i, &kr, &ki, xr, xi));
+ PetscReal error;
+ _try(EPSComputeRelativeError(eps, i, &error));
+#if defined(PETSC_USE_COMPLEX)
+ PetscReal re = PetscRealPart(kr);
+ PetscReal im = PetscImaginaryPart(kr);
+#else
+ PetscReal re = kr;
+ PetscReal im = ki;
+#endif
+ // Output eigenvalues
+ Msg::Info("EIG %03d %s %.16e %s %.16e %3.6e",
+ (nconv-i), (re<0)? "" : " ", re, (im<0)? "" : " ", im, error);
+
+ // Output eigenvector
+ PetscScalar *real_tmp, *imag_tmp;
+ _try(VecGetArray(xr, &real_tmp));
+ _try(VecGetArray(xi, &imag_tmp));
+ }
+
+ // Free work space
+ _try(EPSDestroy(eps));
+ _try(VecDestroy(xr));
+ _try(VecDestroy(xi));
+ _try(SlepcFinalize());
+ Msg::Info("SLEPc done");
+}
+
+#endif
diff --git a/Solver/eigenSolve.h b/Solver/eigenSolve.h
new file mode 100644
index 0000000000..1e02ce35c7
--- /dev/null
+++ b/Solver/eigenSolve.h
@@ -0,0 +1,47 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _EIGEN_SOLVE_H_
+#define _EIGEN_SOLVE_H_
+
+#include <complex>
+#include "GmshConfig.h"
+#include "GmshMessage.h"
+#include "dofManager.h"
+
+#if defined(HAVE_SLEPC)
+
+#include <slepc.h>
+
+class eigenSolve{
+ private:
+ linearSolver *_A, *_B;
+ std::vector<std::complex<double> > _eigenValues;
+ std::vector<std::vector<std::complex<double> > > _eigenVectors;
+ void _try(int ierr) const { CHKERRABORT(PETSC_COMM_WORLD, ierr); }
+ public:
+ eigenSolve(dofManager<double, double> *manager, const char *A,
+ const char *B=0);
+ solve();
+ std::complex<double> getEigenValue(int num){ return _eigenValuesp[num]; }
+ std::vector<std::complex<double> > &getEigenVector(int num){ return _eigenVectors[num]; }
+};
+
+#else
+
+class eigenSolve{
+ private:
+ std::vector<std::complex<double> > _dummy;
+ public:
+ eigenSolve(dofManager<double, double> *manager, const char *A,
+ const char *B=0){}
+ void solve(){}{ Msg::Error("Eigen solver requires SLEPc"); }
+ std::complex<double> getEigenValue(int num){ return 0.; }
+ std::vector<complex<double> > &getEigenVector(int num){ return _dummy; }
+};
+
+#endif
+
+#endif
diff --git a/Solver/elasticityTerm.cpp b/Solver/elasticityTerm.cpp
index 62351800c0..7fd8bbec30 100644
--- a/Solver/elasticityTerm.cpp
+++ b/Solver/elasticityTerm.cpp
@@ -22,7 +22,7 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
double GradsT[100][3];
e->getIntegrationPoints(integrationOrder, &npts, &GP);
- m.set_all(0.);
+ m.setAll(0.);
double FACT = _E / (1 + _nu);
double C11 = FACT * (1 - _nu) / (1 - 2 * _nu);
@@ -53,8 +53,8 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
inv3x3(jac, invjac);
se->gradNodalShapeFunctions(u, v, w, invjac,Grads);
- B.set_all(0.);
- BT.set_all(0.);
+ B.setAll(0.);
+ BT.setAll(0.);
if (se->getShapeEnrichement() == se->getTestEnrichement()){
for (int j = 0; j < nbNodes; j++){
@@ -92,7 +92,7 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
}
}
- BTH.set_all(0.);
+ BTH.setAll(0.);
BTH.gemm(BT, H);
m.gemm(BTH, B, weight * detJ, 1.);
}
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index c120bbb2ad..09af7cf50f 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -63,7 +63,7 @@ class helmholtzTerm : public femTerm<scalar, scalar> {
double Grads[100][3], grads[100][3];
double sf[100];
// set the local matrix to 0
- m.set_all(0.);
+ m.setAll(0.);
// loop over integration points
for (int i = 0; i < npts; i++){
// compute stuff at this point
diff --git a/Solver/linearSystemFull.h b/Solver/linearSystemFull.h
index 327f730320..85a1f0ba93 100644
--- a/Solver/linearSystemFull.h
+++ b/Solver/linearSystemFull.h
@@ -64,7 +64,7 @@ class linearSystemFull : public linearSystem<scalar> {
}
virtual void zeroMatrix()
{
- _a->set_all(0.);
+ _a->setAll(0.);
}
virtual void zeroRightHandSide()
{
@@ -73,7 +73,7 @@ class linearSystemFull : public linearSystem<scalar> {
virtual int systemSolve()
{
if (_b->size())
- _a->lu_solve(*_b, *_x);
+ _a->luSolve(*_b, *_x);
return 1;
}
};
--
GitLab