From 14e965b3f90ecf885f078b04019c5c73e2a1d02d Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Tue, 14 Dec 2004 20:05:44 +0000
Subject: [PATCH] - Replaced the general eigenvalue solver with a version from
Laurent's MatLib code (a bit more stable thanks to floating point checks with
adjustable tolerance)
- Added a solver for symmetric matrices (again from Laurent. Thanks, dude!)
- Added a small driver for 3x3 matrices, that selects the appropriate
solver (sym/nonsym) automatically
---
Numeric/EigSolve.cpp | 1606 +++++++++++++++++++--------------------
Numeric/Makefile | 6 +-
Numeric/Numeric.h | 3 -
Plugin/Eigenvectors.cpp | 35 +-
Plugin/Makefile | 5 +-
doc/CREDITS | 10 +-
6 files changed, 820 insertions(+), 845 deletions(-)
diff --git a/Numeric/EigSolve.cpp b/Numeric/EigSolve.cpp
index 4b73c31c5f..76ecc0ddf3 100644
--- a/Numeric/EigSolve.cpp
+++ b/Numeric/EigSolve.cpp
@@ -1,4 +1,4 @@
-// $Id: EigSolve.cpp,v 1.2 2004-12-08 16:44:33 geuzaine Exp $
+// $Id: EigSolve.cpp,v 1.3 2004-12-14 20:05:43 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -18,565 +18,809 @@
// USA.
//
// Please report all bugs and problems to <gmsh@geuz.org>.
-
-// This file contains a routine that numerically finds the eigenvalues
-// and eigenvectors of a N X N Real Matrix. It is a C++ translation of
-// a Javascript translation of Fortran routines from EISPACK. The
-// Javascript implementation comes from http://www.akiti.ca/Eig3Solv.html
-//
-// References:
//
-// Smith, B.T.; J.M. Boyle; J.J. Dongarra; B.S. Garbow; Y. Ikebe;
-// V.C. Klema; and C.B. Moler.
-// "Matrix Eigensystem Routines--(EISPACK) Guide"
-// Springer-Verlag, Berlin.
-// 1976
-//
-// Garbow, B.S.; J.M. Boyle; J.J. Dongarra; and C.B. Moler.
-// "Matrix Eigensystem Routines--(EISPACK) Guide Extension"
-// Springer-Verlag, Berlin.
-// 1977
+// 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:
+// IMPORTANT! How to use the output of EigSolve():
//
-// If the i-th eigenvalue is real, the i-th COLUMN of the eigenvector
-// Matrix contains the corresponding eigenvector.
+// 1) If the i-th eigenvalue is real, the i-th COLUMN of the
+// eigenvector Matrix contains the corresponding eigenvector.
//
-// If the i-th eigenvalue is complex with positive imaginary part,
+// 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 Error Code. If it does not equal -1, some eigenvalues and
+// Note the return value. If it does not equal 1, some eigenvalues and
// all eigenvectors are meaningless.
-#include "Gmsh.h"
-#include "Numeric.h"
+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 cdivA(int N, double ar, double ai, double br, double bi,
- double **A, int in1, int in2, int in3)
+static void cdiv(double ar,double ai,double br,double bi,double& cr,double& ci)
{
- // Division routine for dividing one complex number into another:
- // This routine does (ar + ai)/(br + bi) and returns the results in
- // the specified elements of the A matrix.
-
- double s, ars, ais, brs, bis;
-
- s = fabs(br) + fabs(bi);
- ars = ar/s;
- ais = ai/s;
- brs = br/s;
- bis = bi/s;
- s = brs*brs + bis*bis;
- A[in1][in2] = (ars*brs + ais*bis)/s;
- A[in1][in3] = (-(ars*bis) + ais*brs)/s;
- return;
-} // End cdivA
-
-static void hqr2(int N, double **A, double **B, int low, int igh,
- double *wi, double *wr, int ierr)
+ 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)
{
- // Computes the eigenvalues and eigenvectors of a real
- // upper-Hessenberg Matrix using the QR method.
-
- double norm = 0.0, p = 0.0, q = 0.0, ra, s = 0.0, sa, t = 0.0;
- double tst1, tst2, vi, vr, w, x, y, zz = 0.0, r = 0.0;
- int k = 0, l = 0, m, mp2, en = igh, incrFlag = 1, its = 0;
- int itn = 30*N, enm2 = 0, na = 0, notlas;
+ static const double RADIX = 16.e0;
+
+ double b2 = RADIX*RADIX;
+ int i,j,ij,ji,k=0,l=n-1;
+ bool last = false;
- for (int i = 0; i < N; i++){ // Store eigenvalues already isolated
- // and compute matrix norm.
- for (int j = k; j < N; j++)
- norm += fabs(A[i][j]);
- k = i;
- if ((i < low) || (i > igh)){
- wi[i] = 0.0;
- wr[i] = A[i][i];
- } //End if (i < low or i > igh)
- } // End for i
-
- // Search next eigenvalues
- while (en >= low){
-
- if (incrFlag) { // Skip this part if incrFlag is set to 0 at very
- // end of while loop
- its = 0;
- na = en - 1;
- enm2 = na - 1;
- } //End if (incrFlag)
+ // 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
- incrFlag = 1;
-
- // Look for single small sub-diagonal element for l = en step -1
- // until low
-
- for (int i = low; i <= en; i++){
- l = en + low - i;
- if (l == low)
- break;
- s = fabs(A[l - 1][l - 1]) + fabs(A[l][l]);
- if (s == 0.0)
- s = norm;
- tst1 = s;
- tst2 = tst1 + fabs(A[l][l - 1]);
- if (tst2 == tst1)
- break;
- } //End for i
-
- x = A[en][en];
-
- if (l == en){ //One root found
- wr[en] = A[en][en] = x + t;
- wi[en] = 0.0;
- en--;
- continue;
- } //End if (l == en)
-
- y = A[na][na];
- w = A[en][na]*A[na][en];
-
- if (l == na){ //Two roots found
- p = (-x + y)/2;
- q = p*p + w;
- zz = sqrt(fabs(q));
- x = A[en][en] = x + t;
- A[na][na] = y + t;
- if (q >= 0.0){//Real Pair
- zz = ((p < 0.0) ? (-zz + p) : (p + zz));
- wr[en] = wr[na] = x + zz;
- if (zz != 0.0)
- wr[en] = -(w/zz) + x;
- wi[en] = wi[na] = 0.0;
- x = A[en][na];
- s = fabs(x) + fabs(zz);
- p = x/s;
- q = zz/s;
- r = sqrt(p*p + q*q);
- p /= r;
- q /= r;
- for (int j = na; j < N; j++){ //Row modification
- zz = A[na][j];
- A[na][j] = q*zz + p*A[en][j];
- A[en][j] = -(p*zz) + q*A[en][j];
- }//End for j
- for (int j = 0; j <= en; j++){ // Column modification
- zz = A[j][na];
- A[j][na] = q*zz + p*A[j][en];
- A[j][en] = -(p*zz) + q*A[j][en] ;
- }//End for j
- for (int j = low; j <= igh; j++){//Accumulate transformations
- zz = B[j][na];
- B[j][na] = q*zz + p*B[j][en];
- B[j][en] = -(p*zz) + q*B[j][en];
- }//End for j
- } //End if (q >= 0.0)
- else {//else q < 0.0
- wr[en] = wr[na] = x + p;
- wi[na] = zz;
- wi[en] = -zz;
- } //End else
- en--;
- en--;
- continue;
- }//End if (l == na)
-
- if (itn == 0){ //Set error; all eigenvalues have not converged
- //after 30 iterations.
- ierr = en + 1;
- return;
- }//End if (itn == 0)
-
- if ((its == 10) || (its == 20)){ //Form exceptional shift
- t += x;
- for (int i = low; i <= en; i++)
- A[i][i] += -x;
- s = fabs(A[en][na]) + fabs(A[na][enm2]);
- y = x = 0.75*s;
- w = -(0.4375*s*s);
- } //End if (its equals 10 or 20)
-
- its++;
- itn--;
-
- // Look for two consecutive small sub-diagonal elements. Do m = en
- // - 2 to l in increments of -1
-
- for (m = enm2; m >= l; m--){
- zz = A[m][m];
- r = -zz + x;
- s = -zz + y;
- p = (-w + r*s)/A[m + 1][m] + A[m][m + 1];
- q = -(zz + r + s) + A[m + 1][m + 1];
- r = A[m + 2][m + 1];
- s = fabs(p) + fabs(q) + fabs(r);
- p /= s;
- q /= s;
- r /= s;
- if (m == l)
- break;
- tst1 = fabs(p) * (fabs(A[m - 1][m - 1]) + fabs(zz) + fabs(A[m + 1][m + 1]));
- tst2 = tst1 + fabs(A[m][m - 1]) * (fabs(q) + fabs(r));
- if (tst1 == tst2)
- break;
- }//End for m
-
- mp2 = m + 2;
-
- for (int i = mp2; i <= en; i++){
- A[i][i - 2] = 0.0;
- if (i == mp2)
- continue;
- A[i][i - 3] = 0.0;
- }//End for i
-
- // Double qr step involving rows l to en and columns m to en.
-
- for (int i = m; i <= na; i++){
- notlas = ((i != na) ? 1 : 0);
- if (i != m){
- p = A[i][i - 1];
- q = A[i + 1][i - 1];
- r = ((notlas) ? A[i + 2][i - 1] : 0.0);
- x = fabs(p) + fabs(q) + fabs(r);
- if (x == 0.0) //Drop through rest of for i loop
- continue;
- p /= x;
- q /= x;
- r /= x;
- } //End if (i != m)
-
- s = sqrt(p*p + q*q + r*r);
- if (p < 0.0)
- s = -s;
-
- if (i != m)
- A[i][i - 1] = -(s*x);
- else {
- if (l != m)
- A[i][i - 1] = -A[i][i - 1];
+ 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--;
}
-
- p += s;
- x = p/s;
- y = q/s;
- zz = r/s;
- q /= p;
- r /= p;
- k = ((i + 3 < en) ? i + 3 : en);
-
- if (notlas){ //Do row modification
- for (int j = i; j < N; j++) {
- p = A[i][j] + q*A[i + 1][j] + r*A[i + 2][j];
- A[i][j] += -(p*x);
- A[i + 1][j] += -(p*y);
- A[i + 2][j] += -(p*zz);
- }//End for j
-
- for (int j = 0; j <= k; j++) {//Do column modification
- p = x*A[j][i] + y*A[j][i + 1] + zz*A[j][i + 2];
- A[j][i] += -p;
- A[j][i + 1] += -(p*q);
- A[j][i + 2] += -(p*r);
- }//End for j
-
- for (int j = low; j <= igh; j++) {//Accumulate transformations
- p = x*B[j][i] + y*B[j][i + 1] + zz*B[j][i + 2];
- B[j][i] += -p;
- B[j][i + 1] += -(p*q);
- B[j][i + 2] += -(p*r);
- } // End for j
- }//End if notlas
-
else {
- for (int j = i; j < N; j++) {//Row modification
- p = A[i][j] + q*A[i + 1][j];
- A[i][j] += -(p*x);
- A[i + 1][j] += -(p*y);
- }//End for j
-
- for (int j = 0; j <= k; j++){//Column modification
- p = x*A[j][i] + y*A[j][i +1];
- A[j][i] += -p;
- A[j][i + 1] += -(p*q);
- }//End for j
-
- for (int j = low; j <= igh; j++){//Accumulate transformations
- p = x*B[j][i] + y*B[j][i +1];
- B[j][i] += -p;
- B[j][i + 1] += -(p*q);
- }//End for j
-
- } //End else if notlas
- }//End for i
- incrFlag = 0;
- }//End while (en >= low)
-
- if (norm == 0.0)
- return;
-
- //Step from (N - 1) to 0 in steps of -1
-
- for (en = (N - 1); en >= 0; en--){
- p = wr[en];
- q = wi[en];
- na = en - 1;
-
- if (q > 0.0)
- continue;
-
- if (q == 0.0){//Real vector
- m = en;
- A[en][en] = 1.0;
-
- for (int j = na; j >= 0; j--){
- w = -p + A[j][j];
- r = 0.0;
- for (int ii = m; ii <= en; ii++)
- r += A[j][ii]*A[ii][en];
-
- if (wi[j] < 0.0){
- zz = w;
- s = r;
- }//End wi[j] < 0.0
- else {//wi[j] >= 0.0
- m = j;
- if (wi[j] == 0.0){
- t = w;
- if (t == 0.0){
- t = tst1 = norm;
- do {
- t *= 0.01;
- tst2 = norm + t;
- } while (tst2 > tst1);
- } //End if t == 0.0
- A[j][en] = -(r/t);
- }//End if wi[j] == 0.0
- else { //wi[j] > 0.0; Solve real equations
- x = A[j][j + 1];
- y = A[j + 1][j];
- q = (-p + wr[j])*(-p + wr[j]) + wi[j]*wi[j];
- A[j][en] = t = (-(zz*r) + x*s)/q;
- A[j + 1][en] = ((fabs(x) > fabs(zz)) ? -(r + w*t)/x : -(s + y*t)/zz);
- }//End else wi[j] > 0.0
-
- // Overflow control
- t = fabs(A[j][en]);
- if (t == 0.0)
- continue; //go up to top of for j loop
- tst1 = t;
- tst2 = tst1 + 1.0/tst1;
- if (tst2 > tst1)
- continue; //go up to top of for j loop
- for (int ii = j; ii <= en; ii++)
- A[ii][en] /= t;
-
- }//End else wi[j] >= 0.0
-
- }//End for j
-
- } //End q == 0.0
-
- else {//else q < 0.0, complex vector
- //Last vector component chosen imaginary so that eigenvector
- //matrix is triangular
- m = na;
-
- if (fabs(A[en][na]) <= fabs(A[na][en]))
- cdivA(N, 0.0, -A[na][en], -p + A[na][na], q, A, na, na, en);
+ 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 {
- A[na][na] = q/A[en][na];
- A[na][en] = -(-p + A[en][en])/A[en][na];
- } //End else (fabs(A[en][na] > fabs(A[na][en])
-
- A[en][na] = 0.0;
- A[en][en] = 1.0;
-
- for (int j = (na - 1); j >= 0; j--) {
- w = -p + A[j][j];
- sa = ra = 0.0;
-
- for (int ii = m; ii <= en; ii++) {
- ra += A[j][ii]*A[ii][na];
- sa += A[j][ii]*A[ii][en];
- } //End for ii
-
- if (wi[j] < 0.0){
- zz = w;
- r = ra;
- s = sa;
- continue;
- } //End if (wi[j] < 0.0)
-
- //else wi[j] >= 0.0
- m = j;
- if (wi[j] == 0.0)
- cdivA(N, -ra, -sa, w, q, A, j, na, en);
- else {//wi[j] > 0.0; solve complex equations
- x = A[j][j + 1];
- y = A[j + 1][j];
- vr = -(q*q) + (-p + wr[j])*(-p + wr[j]) + wi[j]*wi[j];
- vi = (-p + wr[j])*2.0*q;
- if ((vr == 0.0) && (vi == 0.0)){
- tst1 = norm*(fabs(w) + fabs(q) + fabs(x) + fabs(y) + fabs(zz));
- vr = tst1;
- do {
- vr *= 0.01;
- tst2 = tst1 + vr;
- } while (tst2 > tst1);
- } //End if vr and vi == 0.0
- cdivA(N, -(zz*ra) + x*r + q*sa, -(zz*sa + q*ra) + x*s, vr, vi, A, j, na, en);
-
- if (fabs(x) > (fabs(zz) + fabs(q))){
- A[j + 1][na] = (-(ra + w*A[j][na]) + q*A[j][en])/x;
- A[j + 1][en] = -(sa + w*A[j][en] + q*A[j][na])/x;
- }//End if
- else
- cdivA(N, -(r + y*A[j][na]), -(s + y*A[j][en]), zz, q, A, j + 1, na, en);
-
- }//End else wi[j] > 0.0
-
- t = ((fabs(A[j][na]) >= fabs(A[j][en])) ? fabs(A[j][na]) : fabs(A[j][en]));
-
- if (t == 0.0) continue; // go to top of for j loop
-
- tst1 = t;
- tst2 = tst1 + 1.0/tst1;
- if (tst2 > tst1) continue; //go to top of for j loop
-
- for (int ii = j; ii <= en; ii++){
- A[ii][na] /= t;
- A[ii][en] /= t;
- } //End for ii loop
-
- } // End for j
-
- }//End else q < 0.0
-
- }//End for en
+ H[nn1*nm+nn1] = q/H[nn*nm+nn1];
+ H[nn1*nm+nn] = -(H[nn*nm+nn]-p)/H[nn*nm+nn1];
+ }
- //End back substitution. Vectors of isolated roots.
+ for (i=nn-2; i >= 0; i--) {
+ in = i*nm+nn;
+ in1 = i*nm+nn1;
+ w = H[i*nm+i]-p;
- for (int i = 0; i < N; i++){
- if ((i < low) || (i > igh)) {
- for (int j = i; j < N; j++)
- B[i][j] = A[i][j];
- }//End if i
- }//End for i
+ 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];
+ }
- // Multiply by transformation matrix to give vectors of original
- // full matrix.
-
- //Step from (N - 1) to low in steps of -1.
+ if (wi[i] < 0.e0) {
+ z = w;
+ r = ra;
+ s = sa;
+ continue;
+ }
- for (int i = (N - 1); i >= low; i--){
-
- m = ((i < igh) ? i : igh);
-
- for (int ii = low; ii <= igh; ii++){
- zz = 0.0;
- for (int jj = low; jj <= m; jj++)
- zz += B[ii][jj]*A[jj][i];
- B[ii][i] = zz;
- }//End for ii
- }//End of for i loop
-
- return;
-} //End of function hqr2
-
-static void norVecC(int N, double **Z, double *wi)
+ 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)
{
- // 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. I have further modified it for use in this Javascript
- // routine, for use with a column of an array rather than a column
- // vector.
- //
- // Z - eigenvector Matrix
- // wi - eigenvalue vector
-
- double scale, ssq, absxi, dummy, norm;
-
- for (int j = 0; j < N; j++){ //Go through the columns of the vector
- //array
- scale = 0.0;
- ssq = 1.0;
-
- for (int i = 0; i < N; i++){
- if (Z[i][j] != 0){
- absxi = fabs(Z[i][j]);
- dummy = scale/absxi;
+ 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;
- }//End if (scale < absxi)
+ }
else
ssq += 1.0/dummy/dummy;
- }//End if (Z[i][j] != 0)
- } //End for i
+ }
+ }
- if (wi[j] != 0){// If complex eigenvalue, take into account
- // imaginary part of eigenvector
- for (int i = 0; i < N; i++){
- if (Z[i][j + 1] != 0){
- absxi = fabs(Z[i][j + 1]);
- dummy = scale/absxi;
+ 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;
- }//End if (scale < absxi)
+ }
else
ssq += 1.0/dummy/dummy;
- }//End if (Z[i][j + 1] != 0)
- } //End for i
- }//End if (wi[j] != 0)
+ }
+ }
+ }
- norm = scale*sqrt(ssq); //This is the norm of the (possibly
- //complex) vector
+ // this is norm of the (possibly complex) vector
+ double norm = scale*sqrt(ssq);
- for (int i = 0; i < N; i++)
- Z[i][j] /= norm;
+ for (int i = 0; i < n; i++)
+ Z[j*n+i] /= norm;
- if (wi[j] != 0){// If complex eigenvalue, also scale imaginary
- // part of eigenvector
+ if (fabs(wi[j]) > PREC){
+ // If complex eigenvalue, also scale imaginary part of
+ // eigenvector
j++;
- for (int i = 0; i < N; i++)
- Z[i][j] /= norm;
- }//End if (wi[j] != 0)
-
- }// End for j
-
- return;
-} // End norVecC
+ for (int i = 0; i < n; i++)
+ Z[j*n+i] /= norm;
+ }
+ }
+}
-static void swap_elements(int N, double *v, int i, int k)
+int EigSolve(int nm,int n,double *A,double *wr,double *wi,
+ double *v,int *work1,double *work2)
{
- double tmp = v[k];
- v[k] = v[i];
- v[i] = tmp;
+ 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;
}
-static void swap_columns(int N, double **v, int i, int k)
+// 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)
{
- for(int n = 0; n < N; n++){
- double tmp = v[n][k];
- v[n][k] = v[n][i];
- v[n][i] = tmp;
+ static const int NSWMAX = 50;
+
+ int i,j,k,ii,ij,ik,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;
}
-void EigSort(int N, double *wr, double *wi, double **B)
-{
- // Sorts the eigenvalues/vectors in ascending order according to
- // their real part. Warning: this might screw up the ordering of
- // complex eigenvectors --- VERIFY THIS
+// Sort the eigenvalues/vectors in ascending order according to their
+// real part. Warning: this will screw up the ordering of complex
+// eigenvectors.
- for (int i = 0; i < N - 1; i++){
+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++){
+ for (int j = i + 1; j < n; j++){
const double ej = wr[j];
if(ej < ek){
k = j;
@@ -584,310 +828,60 @@ void EigSort(int N, double *wr, double *wi, double **B)
}
}
if (k != i){
- swap_elements(N, wr, i, k);
- swap_elements(N, wi, i, k);
- swap_columns(N, B, i, k);
+ 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);
}
}
}
-int EigSolve(int N, double **A, double *wr, double *wi, double **B,
- int *itmp, double *dtmp)
-{
- // Computes the N eigenvalues of the square, real matrix A. All
- // vectors have dimension N and all matrices have dimension
- // NxN. Warning: the matrix A gets modified during the
- // computation.
-
- // Balance the matrix to improve accuracy of
- // eigenvalues. (Introduces no rounding errors, since it scales A by
- // powers of the radix.)
-
- double *scale = dtmp; //Contains information about transformations.
- int *trace = itmp; //Records row and column interchanges
-
- double radix = 2; //Base of machine floating point representation.
- double c, f, g, r, s, b2 = radix*radix;
- int ierr = -1, igh, low, k = 0, l = N - 1, noconv;
-
- //Search through rows, isolating eigenvalues and pushing them down.
-
- noconv = l;
-
- while (noconv >= 0){
- r = 0;
-
- for (int j = 0; j <= l; j++) {
- if (j == noconv) continue;
- if (A[noconv][j] != 0.0){
- r = 1;
- break;
- }
- } //End for j
-
- if (r == 0){
- scale[l] = noconv;
-
- if (noconv != l){
- for (int i = 0; i <= l; i++){
- f = A[i][noconv];
- A[i][noconv] = A[i][l];
- A[i][l] = f;
- }//End for i
- for (int i = 0; i < N; i++){
- f = A[noconv][i];
- A[noconv][i] = A[l][i];
- A[l][i] = f;
- }//End for i
- }//End if (noconv != l)
-
- if (l == 0)
- break; //break out of while loop
-
- noconv = --l;
-
- }//End if (r == 0)
-
- else //else (r != 0)
- noconv--;
-
- }//End while noconv
-
- if (l > 0) { //Search through columns, isolating eigenvalues and
- //pushing them left.
-
- noconv = 0;
-
- while (noconv <= l){
- c = 0;
-
- for (int i = k; i <= l; i++){
- if (i == noconv) continue;
- if (A[i][noconv] != 0.0){
- c = 1;
- break;
- }
- }//End for i
-
- if (c == 0){
- scale[k] = noconv;
-
- if (noconv != k){
- for (int i = 0; i <= l; i++){
- f = A[i][noconv];
- A[i][noconv] = A[i][k];
- A[i][k] = f;
- }//End for i
- for (int i = k; i < N; i++){
- f = A[noconv][i];
- A[noconv][i] = A[k][i];
- A[k][i] = f;
- }//End for i
-
- }//End if (noconv != k)
-
- noconv = ++k;
- }//End if (c == 0)
-
- else //else (c != 0)
- noconv++;
-
- }//End while noconv
-
- //Balance the submatrix in rows k through l.
-
- for (int i = k; i <= l; i++)
- scale[i] = 1.0;
-
- //Iterative loop for norm reduction
- do {
- noconv = 0;
- for (int i = k; i <= l; i++) {
- c = r = 0.0;
- for (int j = k; j <= l; j++){
- if (j == i) continue;
- c += fabs(A[j][i]);
- r += fabs(A[i][j]);
- } // End for j
- //guard against zero c or r due to underflow:
- if ((c == 0.0) || (r == 0.0)) continue;
- g = r/radix;
- f = 1.0;
- s = c + r;
- while (c < g) {
- f *= radix;
- c *= b2;
- } // End while (c < g)
- g = r*radix;
- while (c >= g) {
- f /= radix;
- c /= b2;
- } // End while (c >= g)
-
- //Now balance
- if ((c + r)/f < 0.95*s) {
- g = 1.0/f;
- scale[i] *= f;
- noconv = 1;
- for (int j = k; j < N; j++)
- A[i][j] *= g;
- for (int j = 0; j <= l; j++)
- A[j][i] *= f;
- } //End if ((c + r)/f < 0.95*s)
- } // End for i
- } while (noconv); // End of do-while loop.
-
- } //End if (l > 0)
-
- low = k;
- igh = l;
-
- //End of balanc
-
- // Now reduce the real general Matrix to upper-Hessenberg form using
- // stabilized elementary similarity transformations.
-
- for (int i = (low + 1); i < igh; i++){
- k = i;
- c = 0.0;
-
- for (int j = i; j <= igh; j++){
- if (fabs(A[j][i - 1]) > fabs(c)){
- c = A[j][i - 1];
- k = j;
- }//End if
- }//End for j
-
- trace[i] = k;
-
- if (k != i){
- for (int j = (i - 1); j < N; j++){
- r = A[k][j];
- A[k][j] = A[i][j];
- A[i][j] = r;
- }//End for j
-
- for (int j = 0; j <= igh; j++){
- r = A[j][k];
- A[j][k] = A[j][i];
- A[j][i] = r;
- }//End for j
- }//End if (k != i)
-
- if (c != 0.0){
- for (int m = (i + 1); m <= igh; m++){
- r = A[m][i - 1];
-
- if (r != 0.0){
- r = A[m][i - 1] = r/c;
- for (int j = i; j < N; j++)
- A[m][j] += -(r*A[i][j]);
- for (int j = 0; j <= igh; j++)
- A[j][i] += r*A[j][m];
- }//End if (r != 0.0)
- }//End for m
- }//End if (c != 0)
- } //End for i.
-
- // Accumulate the stabilized elementary similarity transformations
- // used in the reduction of A to upper-Hessenberg form. Introduces
- // no rounding errors since it only transfers the multipliers used
- // in the reduction process into the eigenvector matrix.
-
- for (int i = 0; i < N; i++){ //Initialize B to the identity Matrix.
- for (int j = 0; j < N; j++)
- B[i][j] = 0.0;
- B[i][i] = 1.0;
- } //End for i
-
- for (int i = (igh - 1); i >= (low + 1); i--){
- k = trace[i];
- for (int j = (i + 1); j <= igh; j++)
- B[j][i] = A[j][i - 1];
-
- if (i == k)
- continue;
-
- for (int j = i; j <= igh; j++){
- B[i][j] = B[k][j];
- B[k][j] = 0.0;
- } //End for j
-
- B[k][i] = 1.0;
-
- } // End for i
-
- hqr2(N, A, B, low, igh, wi, wr, ierr);
-
- if (ierr == -1){
-
- if (low != igh){
- for (int i = low; i <= igh; i++){
- s = scale[i];
- for (int j = 0; j < N; j++)
- B[i][j] *= s;
- }//End for i
- }//End if (low != igh)
-
- for (int i = (low - 1); i >= 0; i--){
- k = (int)scale[i];//FIXME: I added the cast
- if (k != i){
- for (int j = 0; j < N; j++){
- s = B[i][j];
- B[i][j] = B[k][j];
- B[k][j] = s;
- }//End for j
- }//End if k != i
- }//End for i
-
- for (int i = (igh + 1); i < N; i++){
- k = (int)scale[i];//FIXME: Iadded the casr
- if (k != i){
- for (int j = 0; j < N; j++){
- s = B[i][j];
- B[i][j] = B[k][j];
- B[k][j] = s;
- }//End for j
- }//End if k != i
- }//End for i
-
- norVecC(N, B, wi); //Normalize the eigenvectors
-
- }//End if ierr = -1
+// 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;
-} //End of Eig3Solve
+}
# if 0
-int main()
+int main ()
{
- const int N = 3;
- double *wr = new double[N];
- double *wi = new double[N];
- double **A = new double*[N];
- double **B = new double*[N];
- for(int i = 0; i < N; i++){
- A[i] = new double[N];
- B[i] = new double[N];
- }
- double *dtmp = new double[N];
- int *itmp = new int[N];
-
- for(int i = 0; i < 1000000; i++){ // perf test
- A[0][0] = 1.; A[0][1] = 2.; A[0][2] = 3.;
- A[1][0] = 2.; A[1][1] = 4.; A[1][2] = 5.;
- A[2][0] = 3.; A[2][1] = 5.; A[2][2] = 6.;
- EigSolve(3, A, wr, wi, B, itmp, dtmp);
- }
- /*
- printf("%g %g %g\n", A[0][0], A[0][1], A[0][2]);
- printf("%g %g %g\n", A[1][0], A[1][1], A[1][2]);
- printf("%g %g %g\n", A[2][0], A[2][1], A[2][2]);
- */
- printf("real= %g %g %g \n", wr[0], wr[1], wr[2]);
- printf("imag= %g %g %g \n", wi[0], wi[1], wi[2]);
- printf("eigvec1 = %g %g %g \n", B[0][0], B[1][0], B[2][0]);
- printf("eigvec2 = %g %g %g \n", B[0][1], B[1][1], B[2][1]);
- printf("eigvec3 = %g %g %g \n", B[0][2], B[1][2], B[2][2]);
+ //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/Makefile b/Numeric/Makefile
index 66720267cd..aad54994c8 100644
--- a/Numeric/Makefile
+++ b/Numeric/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.16 2004-12-08 05:38:56 geuzaine Exp $
+# $Id: Makefile,v 1.17 2004-12-14 20:05:43 geuzaine Exp $
#
# Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
#
@@ -56,9 +56,7 @@ depend:
Numeric.o: Numeric.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h Numeric.h
-EigSolve.o: EigSolve.cpp ../Common/Gmsh.h ../Common/Message.h \
- ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
- ../DataStr/avl.h ../DataStr/Tools.h Numeric.h
+EigSolve.o: EigSolve.cpp
gsl_newt.o: gsl_newt.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h Numeric.h
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 4e59c3bb82..d065d004ba 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -73,9 +73,6 @@ double angle_02pi(double A3);
void eigenvalue(double mat[3][3], double re[3]);
void FindCubicRoots(const double coeff[4], double re[3], double im[3]);
void eigsort(double d[3]);
-int EigSolve(int N, double **A, double *wr, double *wi, double **B,
- int *itmp, double *dtmp);
-void EigSort(int N, double *wr, double *wi, double **B);
double InterpolateIso(double *X, double *Y, double *Z,
double *Val, double V, int I1, int I2,
diff --git a/Plugin/Eigenvectors.cpp b/Plugin/Eigenvectors.cpp
index 69bbd5a04c..9c50ae69b7 100644
--- a/Plugin/Eigenvectors.cpp
+++ b/Plugin/Eigenvectors.cpp
@@ -1,4 +1,4 @@
-// $Id: Eigenvectors.cpp,v 1.2 2004-12-10 03:42:29 geuzaine Exp $
+// $Id: Eigenvectors.cpp,v 1.3 2004-12-14 20:05:43 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -26,6 +26,7 @@
#include "Context.h"
#include "Malloc.h"
#include "Numeric.h"
+#include "EigSolve.h"
extern Context_T CTX;
@@ -89,7 +90,7 @@ void GMSH_EigenvectorsPlugin::catchErrorMessage(char *errorMessage) const
static int nonzero(double v[3])
{
for(int i = 0; i < 3; i++)
- if(fabs(v[i]) > 1.e-15) return 1;
+ if(fabs(v[i]) > 1.e-16) return 1;
return 0;
}
@@ -101,14 +102,7 @@ static void eigenvectors(List_T *inList, int inNb,
{
if(!inNb) return;
- int itmp[3], nbcomplex = 0;
- double wr[3], wi[3], dtmp[3];
- double **A = new double*[3];
- double **B = new double*[3];
- for(int i = 0; i < 3; i++){
- A[i] = new double[3];
- B[i] = new double[3];
- }
+ int nbcomplex = 0;
int nb = List_Nbr(inList) / inNb;
for(int i = 0; i < List_Nbr(inList); i += nb) {
@@ -121,11 +115,9 @@ static void eigenvectors(List_T *inList, int inNb,
for(int k = 0; k < nbNod; k++){
double *val = (double *)List_Pointer_Fast(inList, i + 3 * nbNod +
nbNod * 9 * j + 9 * k);
- A[0][0] = val[0]; A[0][1] = val[1]; A[0][2] = val[2];
- A[1][0] = val[3]; A[1][1] = val[4]; A[1][2] = val[5];
- A[2][0] = val[6]; A[2][1] = val[7]; A[2][2] = val[8];
- EigSolve(3, A, wr, wi, B, itmp, dtmp);
- EigSort(3, wr, wi, B);
+ double wr[3], wi[3], B[9];
+ if(!EigSolve3x3(val, wr, wi, B))
+ Msg(GERROR, "Eigensolver failed to converge");
nbcomplex += nonzero(wi);
if(!scale)
wr[0] = wr[1] = wr[2] = 1.;
@@ -133,9 +125,9 @@ static void eigenvectors(List_T *inList, int inNb,
double res;
// wrong if there are complex eigenvals (B contains both
// real and imag parts: cf. explanation in EigSolve.cpp)
- res = wr[0] * B[l][0]; List_Add(minList, &res);
- res = wr[1] * B[l][1]; List_Add(midList, &res);
- res = wr[2] * B[l][2]; List_Add(maxList, &res);
+ res = wr[0] * B[l]; List_Add(minList, &res);
+ res = wr[1] * B[3+l]; List_Add(midList, &res);
+ res = wr[2] * B[6+l]; List_Add(maxList, &res);
}
}
}
@@ -144,13 +136,6 @@ static void eigenvectors(List_T *inList, int inNb,
(*maxNb)++;
}
- for(int i = 0; i < 3; i++){
- delete [] A[i];
- delete [] B[i];
- }
- delete [] A;
- delete [] B;
-
if(nbcomplex)
Msg(GERROR, "%d tensors have complex eigenvalues/eigenvectors", nbcomplex);
}
diff --git a/Plugin/Makefile b/Plugin/Makefile
index 06ff391c00..b8c9c2afc1 100644
--- a/Plugin/Makefile
+++ b/Plugin/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.64 2004-12-10 03:35:29 geuzaine Exp $
+# $Id: Makefile,v 1.65 2004-12-14 20:05:44 geuzaine Exp $
#
# Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
#
@@ -133,7 +133,8 @@ Eigenvectors.o: Eigenvectors.cpp Plugin.h ../Common/Options.h \
../Common/Message.h ../Common/Views.h ../Common/ColorTable.h \
../DataStr/List.h ../Common/VertexArray.h ../Common/SmoothNormals.h \
../Common/GmshMatrix.h ../Common/AdaptiveViews.h Eigenvectors.h \
- ../Common/Context.h ../DataStr/Malloc.h ../Numeric/Numeric.h
+ ../Common/Context.h ../DataStr/Malloc.h ../Numeric/Numeric.h \
+ ../Numeric/EigSolve.h
Octree.o: Octree.cpp Octree.h OctreeInternals.h
OctreeInternals.o: OctreeInternals.cpp ../Common/Message.h \
OctreeInternals.h
diff --git a/doc/CREDITS b/doc/CREDITS
index 0cc576a75e..b46045990c 100644
--- a/doc/CREDITS
+++ b/doc/CREDITS
@@ -1,4 +1,4 @@
-$Id: CREDITS,v 1.21 2004-12-06 07:13:47 geuzaine Exp $
+$Id: CREDITS,v 1.22 2004-12-14 20:05:44 geuzaine Exp $
Gmsh is copyright (C) 1997-2004
@@ -17,10 +17,10 @@ to elementary geometry" interface, Netgen integration).
Other code contributors include: David Colignon <David.Colignon at
univ.u-3mrs.fr> for new colormaps; Patrick Dular <patrick.dular at
ulg.ac.be> for transfinite mesh bug fixes; Laurent Stainier
-<l.stainier at ulg.ac.be> for help with the MacOS port and the tensor
-display code; Pierre Badel <badel at freesurf.fr> for help with the
-GSL integration; and Marc Ume <Marc.Ume at digitalgraphics.be> for the
-original list code.
+<l.stainier at ulg.ac.be> for the eigenvalue solvers and for help with
+the MacOS port and the tensor display code; Pierre Badel <badel at
+freesurf.fr> for help with the GSL integration; and Marc Ume <Marc.Ume
+at digitalgraphics.be> for the original list code.
The AVL tree code (DataStr/avl.*) and the YUV image code
(Graphics/gl2yuv.*) are copyright (C) 1988-1993, 1995 The Regents of
--
GitLab