--- trunk/src/math/SquareMatrix.hpp 2004/10/14 23:28:09 76 +++ trunk/src/math/SquareMatrix.hpp 2004/10/20 18:07:08 123 @@ -29,7 +29,7 @@ * @date 10/11/2004 * @version 1.0 */ -#ifndef MATH_SQUAREMATRIX_HPP + #ifndef MATH_SQUAREMATRIX_HPP #define MATH_SQUAREMATRIX_HPP #include "math/RectMatrix.hpp" @@ -78,22 +78,28 @@ namespace oopse { return m; } - /** Retunrs the inversion of this matrix. */ + /** + * Retunrs the inversion of this matrix. + * @todo need implementation + */ SquareMatrix inverse() { SquareMatrix result; return result; } - /** Returns the determinant of this matrix. */ - double determinant() const { - double det; + /** + * Returns the determinant of this matrix. + * @todo need implementation + */ + Real determinant() const { + Real det; return det; } /** Returns the trace of this matrix. */ - double trace() const { - double tmp = 0; + Real trace() const { + Real tmp = 0; for (unsigned int i = 0; i < Dim ; i++) tmp += data_[i][i]; @@ -142,189 +148,244 @@ namespace oopse { return true; } + /** @todo need implementation */ void diagonalize() { - jacobi(m, eigenValues, ortMat); + //jacobi(m, eigenValues, ortMat); } /** - * Finds the eigenvalues and eigenvectors of a symmetric matrix - * @param eigenvals a reference to a vector3 where the - * eigenvalues will be stored. The eigenvalues are ordered so - * that eigenvals[0] <= eigenvals[1] <= eigenvals[2]. - * @return an orthogonal matrix whose ith column is an - * eigenvector for the eigenvalue eigenvals[i] - */ - SquareMatrix findEigenvectors(Vector& eigenValues) { - SquareMatrix ortMat; - - if ( !isSymmetric()){ - throw(); - } - - SquareMatrix m(*this); - jacobi(m, eigenValues, ortMat); - - return ortMat; - } - /** * Jacobi iteration routines for computing eigenvalues/eigenvectors of * real symmetric matrix * * @return true if success, otherwise return false - * @param a source matrix - * @param w output eigenvalues - * @param v output eigenvectors + * @param a symmetric matrix whose eigenvectors are to be computed. On return, the matrix is + * overwritten + * @param w will contain the eigenvalues of the matrix On return of this function + * @param v the columns of this matrix will contain the eigenvectors. The eigenvectors are + * normalized and mutually orthogonal. */ - void jacobi(const SquareMatrix& a, - Vector& w, + + static int jacobi(SquareMatrix& a, Vector& d, SquareMatrix& v); };//end SquareMatrix -#define ROT(a,i,j,k,l) g=a(i, j);h=a(k, l);a(i, j)=g-s*(h+g*tau);a(k, l)=h+s*(g-h*tau) -#define MAX_ROTATIONS 60 +/*========================================================================= -template -void SquareMatrix::jacobi(SquareMatrix& a, - Vector& w, - SquareMatrix& v) { - const int N = Dim; - int i, j, k, iq, ip; - double tresh, theta, tau, t, sm, s, h, g, c; - double tmp; - Vector b, z; + Program: Visualization Toolkit + Module: $RCSfile: SquareMatrix.hpp,v $ - // initialize - for (ip=0; ip 4 && (fabs(w(ip))+g) == fabs(w(ip)) - && (fabs(w(iq))+g) == fabs(w(iq))) - { - a(ip, iq) = 0.0; - } - else if (fabs(a(ip, iq)) > tresh) - { - h = w(iq) - w(ip); - if ( (fabs(h)+g) == fabs(h)) t = (a(ip, iq)) / h; - else - { - theta = 0.5*h / (a(ip, iq)); - t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta)); - if (theta < 0.0) t = -t; - } - c = 1.0 / sqrt(1+t*t); - s = t*c; - tau = s/(1.0+c); - h = t*a(ip, iq); - z(ip) -= h; - z(iq) += h; - w(ip) -= h; - w(iq) += h; - a(ip, iq)=0.0; - for (j=0;j= MAX_ROTATIONS ) - return false; + // Jacobi iteration for the solution of eigenvectors/eigenvalues of a nxn + // real symmetric matrix. Square nxn matrix a; size of matrix in n; + // output eigenvalues in w; and output eigenvectors in v. Resulting + // eigenvalues/vectors are sorted in decreasing order; eigenvectors are + // normalized. + template + int SquareMatrix::jacobi(SquareMatrix& a, Vector& w, + SquareMatrix& v) { + const int n = Dim; + int i, j, k, iq, ip, numPos; + Real tresh, theta, tau, t, sm, s, h, g, c, tmp; + Real bspace[4], zspace[4]; + Real *b = bspace; + Real *z = zspace; - // sort eigenfunctions - for (j=0; j= tmp) - { - k = i; - tmp = w(k); - } - } - if (k != j) - { - w(k) = w(j); - w(j) = tmp; - for (i=0; i 4) + { + b = new Real[n]; + z = new Real[n]; + } - // insure eigenvector consistency (i.e., Jacobi can compute - // vectors that are negative of one another (.707,.707,0) and - // (-.707,-.707,0). This can reek havoc in - // hyperstreamline/other stuff. We will select the most - // positive eigenvector. - int numPos; - for (j=0; j= 0.0 ) numPos++; - if ( numPos < 2 ) for(i=0; i 3 && (fabs(w[ip])+g) == fabs(w[ip]) + && (fabs(w[iq])+g) == fabs(w[iq])) + { + a(ip, iq) = 0.0; + } + else if (fabs(a(ip, iq)) > tresh) + { + h = w[iq] - w[ip]; + if ( (fabs(h)+g) == fabs(h)) + { + t = (a(ip, iq)) / h; + } + else + { + theta = 0.5*h / (a(ip, iq)); + t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta)); + if (theta < 0.0) + { + t = -t; + } + } + c = 1.0 / sqrt(1+t*t); + s = t*c; + tau = s/(1.0+c); + h = t*a(ip, iq); + z[ip] -= h; + z[iq] += h; + w[ip] -= h; + w[iq] += h; + a(ip, iq)=0.0; + // ip already shifted left by 1 unit + for (j = 0;j <= ip-1;j++) + { + VTK_ROTATE(a,j,ip,j,iq); + } + // ip and iq already shifted left by 1 unit + for (j = ip+1;j <= iq-1;j++) + { + VTK_ROTATE(a,ip,j,j,iq); + } + // iq already shifted left by 1 unit + for (j=iq+1; j= VTK_MAX_ROTATIONS ) + { + std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions" << std::endl; + return 0; + } + + // sort eigenfunctions these changes do not affect accuracy + for (j=0; j= tmp) // why exchage if same? + { + k = i; + tmp = w[k]; + } + } + if (k != j) + { + w[k] = w[j]; + w[j] = tmp; + for (i=0; i> 1) + (n & 1); + for (j=0; j= 0.0 ) + { + numPos++; + } + } + // if ( numPos < ceil(double(n)/double(2.0)) ) + if ( numPos < ceil_half_n) + { + for(i=0; i 4) + { + delete [] b; + delete [] z; + } + return 1; + } + + } #endif //MATH_SQUAREMATRIX_HPP +