src/math/SquareMatrix.hpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
 
/**
 * @file SquareMatrix.hpp
 * @author Teng Lin
 * @date 10/11/2004
 * @version 1.0
 */
#ifndef MATH_SQUAREMATRIX_HPP 
#define MATH_SQUAREMATRIX_HPP 

#include "math/RectMatrix.hpp"

namespace oopse {

  /**
   * @class SquareMatrix SquareMatrix.hpp "math/SquareMatrix.hpp"
   * @brief A square matrix class
   * @template Real the element type
   * @template Dim the dimension of the square matrix
   */
  template<typename Real, int Dim>
  class SquareMatrix : public RectMatrix<Real, Dim, Dim> {
  public:
    typedef Real ElemType;
    typedef Real* ElemPoinerType;

    /** default constructor */
    SquareMatrix() {
      for (unsigned int i = 0; i < Dim; i++)
        for (unsigned int j = 0; j < Dim; j++)
          this->data_[i][j] = 0.0;
    }

    /** Constructs and initializes every element of this matrix to a scalar */ 
    SquareMatrix(Real s) : RectMatrix<Real, Dim, Dim>(s){
    }

    /** Constructs and initializes from an array */ 
    SquareMatrix(Real* array) : RectMatrix<Real, Dim, Dim>(array){
    }


    /** copy constructor */
    SquareMatrix(const RectMatrix<Real, Dim, Dim>& m) : RectMatrix<Real, Dim, Dim>(m) {
    }
            
    /** copy assignment operator */
    SquareMatrix<Real, Dim>& operator =(const RectMatrix<Real, Dim, Dim>& m) {
      RectMatrix<Real, Dim, Dim>::operator=(m);
      return *this;
    }
                                   
    /** Retunrs  an identity matrix*/

    static SquareMatrix<Real, Dim> identity() {
      SquareMatrix<Real, Dim> m;
                
      for (unsigned int i = 0; i < Dim; i++) 
        for (unsigned int j = 0; j < Dim; j++) 
          if (i == j)
            m(i, j) = 1.0;
          else
            m(i, j) = 0.0;

      return m;
    }

    /** 
     * Retunrs  the inversion of this matrix. 
     * @todo need implementation
     */
    SquareMatrix<Real, Dim>  inverse() {
      SquareMatrix<Real, Dim> result;

      return result;
    }        

    /**
     * Returns the determinant of this matrix.
     * @todo need implementation
     */
    Real determinant() const {
      Real det;
      return det;
    }

    /** Returns the trace of this matrix. */
    Real trace() const {
      Real tmp = 0;
               
      for (unsigned int i = 0; i < Dim ; i++)
        tmp += this->data_[i][i];

      return tmp;
    }

    /** Tests if this matrix is symmetrix. */            
    bool isSymmetric() const {
      for (unsigned int i = 0; i < Dim - 1; i++)
        for (unsigned int j = i; j < Dim; j++)
          if (fabs(this->data_[i][j] - this->data_[j][i]) > oopse::epsilon) 
            return false;
                        
      return true;
    }

    /** Tests if this matrix is orthogonal. */            
    bool isOrthogonal() {
      SquareMatrix<Real, Dim> tmp;

      tmp = *this * transpose();

      return tmp.isDiagonal();
    }

    /** Tests if this matrix is diagonal. */
    bool isDiagonal() const {
      for (unsigned int i = 0; i < Dim ; i++)
        for (unsigned int j = 0; j < Dim; j++)
          if (i !=j && fabs(this->data_[i][j]) > oopse::epsilon) 
            return false;
                        
      return true;
    }

    /** Tests if this matrix is the unit matrix. */
    bool isUnitMatrix() const {
      if (!isDiagonal())
        return false;
                
      for (unsigned int i = 0; i < Dim ; i++)
        if (fabs(this->data_[i][i] - 1) > oopse::epsilon)
          return false;
                    
      return true;
    }         

    /** Return the transpose of this matrix */
    SquareMatrix<Real,  Dim> transpose() const{
      SquareMatrix<Real,  Dim> result;
                
      for (unsigned int i = 0; i < Dim; i++)
        for (unsigned int j = 0; j < Dim; j++)              
          result(j, i) = this->data_[i][j];

      return result;
    }
            
    /** @todo need implementation */
    void diagonalize() {
      //jacobi(m, eigenValues, ortMat);
    }

    /**
     * Jacobi iteration routines for computing eigenvalues/eigenvectors of 
     * real symmetric matrix
     *
     * @return true if success, otherwise return false
     * @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. 
     */
           
    static int jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& d, 
                      SquareMatrix<Real, Dim>& v);
  };//end SquareMatrix


  /*=========================================================================

  Program:   Visualization Toolkit
  Module:    $RCSfile: SquareMatrix.hpp,v $

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

  This software is distributed WITHOUT ANY WARRANTY; without even
  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  PURPOSE.  See the above copyright notice for more information.

  =========================================================================*/

#define VTK_ROTATE(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 VTK_MAX_ROTATIONS 20

  // 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<typename Real, int Dim>
  int SquareMatrix<Real, Dim>::jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w, 
                                      SquareMatrix<Real, Dim>& 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;

    // only allocate memory if the matrix is large
    if (n > 4) {
      b = new Real[n];
      z = new Real[n]; 
    }

    // initialize
    for (ip=0; ip<n; ip++) {
      for (iq=0; iq<n; iq++) {
        v(ip, iq) = 0.0;
      }
      v(ip, ip) = 1.0;
    }
    for (ip=0; ip<n; ip++) {
      b[ip] = w[ip] = a(ip, ip);
      z[ip] = 0.0;
    }

    // begin rotation sequence
    for (i=0; i<VTK_MAX_ROTATIONS; i++) {
      sm = 0.0;
      for (ip=0; ip<n-1; ip++) {
        for (iq=ip+1; iq<n; iq++) {
          sm += fabs(a(ip, iq));
        }
      }
      if (sm == 0.0) {
        break;
      }

      if (i < 3) {                                // first 3 sweeps
        tresh = 0.2*sm/(n*n);
      } else {
        tresh = 0.0;
      }

      for (ip=0; ip<n-1; ip++) {
        for (iq=ip+1; iq<n; iq++) {
          g = 100.0*fabs(a(ip, iq));

          // after 4 sweeps
          if (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<n; j++) {
              VTK_ROTATE(a,ip,j,iq,j);
            }
            for (j=0; j<n; j++) {
              VTK_ROTATE(v,j,ip,j,iq);
            }
          }
        }
      }

      for (ip=0; ip<n; ip++) {
        b[ip] += z[ip];
        w[ip] = b[ip];
        z[ip] = 0.0;
      }
    }

    //// this is NEVER called
    if ( i >= 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<n-1; j++) {                  // boundary incorrect
      k = j;
      tmp = w[k];
      for (i=j+1; i<n; i++) {                // boundary incorrect, shifted already
        if (w[i] >= tmp) {                   // why exchage if same?
          k = i;
          tmp = w[k];
        }
      }
      if (k != j) {
        w[k] = w[j];
        w[j] = tmp;
        for (i=0; i<n; i++) {
          tmp = v(i, j);
          v(i, j) = v(i, k);
          v(i, k) = tmp;
        }
      }
    }
    // 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 ceil_half_n = (n >> 1) + (n & 1);
    for (j=0; j<n; j++) {
      for (numPos=0, i=0; i<n; i++) {
        if ( v(i, j) >= 0.0 ) {
          numPos++;
        }
      }
      //    if ( numPos < ceil(double(n)/double(2.0)) )
      if ( numPos < ceil_half_n) {
        for (i=0; i<n; i++) {
          v(i, j) *= -1.0;
        }
      }
    }

    if (n > 4) {
      delete [] b;
      delete [] z;
    }
    return 1;
  }


}
#endif //MATH_SQUAREMATRIX_HPP 

Revision:	507
Committed:	Fri Apr 15 22:04:00 2005 UTC (20 years ago) by gezelter
Original Path:	*trunk/src/math/SquareMatrix.hpp*
File size:	11236 byte(s)
Log Message:	xemacs has been drafted to perform our indentation services
#	User	Rev	Content
1	gezelter	507	/*
2	gezelter	246	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	tim	70	*
4	gezelter	246	* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9			* 1. Acknowledgement of the program authors must be made in any
10			* publication of scientific results based in part on use of the
11			* program. An acceptable form of acknowledgement is citation of
12			* the article in which the program was described (Matthew
13			* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			* Parallel Simulation Engine for Molecular Dynamics,"
16			* J. Comput. Chem. 26, pp. 252-271 (2005))
17			*
18			* 2. Redistributions of source code must retain the above copyright
19			* notice, this list of conditions and the following disclaimer.
20			*
21			* 3. Redistributions in binary form must reproduce the above copyright
22			* notice, this list of conditions and the following disclaimer in the
23			* documentation and/or other materials provided with the
24			* distribution.
25			*
26			* This software is provided "AS IS," without a warranty of any
27			* kind. All express or implied conditions, representations and
28			* warranties, including any implied warranty of merchantability,
29			* fitness for a particular purpose or non-infringement, are hereby
30			* excluded. The University of Notre Dame and its licensors shall not
31			* be liable for any damages suffered by licensee as a result of
32			* using, modifying or distributing the software or its
33			* derivatives. In no event will the University of Notre Dame or its
34			* licensors be liable for any lost revenue, profit or data, or for
35			* direct, indirect, special, consequential, incidental or punitive
36			* damages, however caused and regardless of the theory of liability,
37			* arising out of the use of or inability to use software, even if the
38			* University of Notre Dame has been advised of the possibility of
39			* such damages.
40	tim	70	*/
41	gezelter	246
42	tim	70	/**
43			* @file SquareMatrix.hpp
44			* @author Teng Lin
45			* @date 10/11/2004
46			* @version 1.0
47			*/
48	gezelter	507	#ifndef MATH_SQUAREMATRIX_HPP
49	tim	70	#define MATH_SQUAREMATRIX_HPP
50
51	tim	74	#include "math/RectMatrix.hpp"
52	tim	70
53			namespace oopse {
54
55	gezelter	507	/**
56			* @class SquareMatrix SquareMatrix.hpp "math/SquareMatrix.hpp"
57			* @brief A square matrix class
58			* @template Real the element type
59			* @template Dim the dimension of the square matrix
60			*/
61			template<typename Real, int Dim>
62			class SquareMatrix : public RectMatrix<Real, Dim, Dim> {
63			public:
64			typedef Real ElemType;
65			typedef Real* ElemPoinerType;
66	tim	70
67	gezelter	507	/** default constructor */
68			SquareMatrix() {
69			for (unsigned int i = 0; i < Dim; i++)
70			for (unsigned int j = 0; j < Dim; j++)
71			this->data_[i][j] = 0.0;
72			}
73	tim	70
74	gezelter	507	/** Constructs and initializes every element of this matrix to a scalar */
75			SquareMatrix(Real s) : RectMatrix<Real, Dim, Dim>(s){
76			}
77	tim	151
78	gezelter	507	/** Constructs and initializes from an array */
79			SquareMatrix(Real* array) : RectMatrix<Real, Dim, Dim>(array){
80			}
81	tim	151
82
83	gezelter	507	/** copy constructor */
84			SquareMatrix(const RectMatrix<Real, Dim, Dim>& m) : RectMatrix<Real, Dim, Dim>(m) {
85			}
86	tim	70
87	gezelter	507	/** copy assignment operator */
88			SquareMatrix<Real, Dim>& operator =(const RectMatrix<Real, Dim, Dim>& m) {
89			RectMatrix<Real, Dim, Dim>::operator=(m);
90			return *this;
91			}
92	tim	137
93	gezelter	507	/** Retunrs an identity matrix*/
94	tim	74
95	gezelter	507	static SquareMatrix<Real, Dim> identity() {
96			SquareMatrix<Real, Dim> m;
97	tim	137
98	gezelter	507	for (unsigned int i = 0; i < Dim; i++)
99			for (unsigned int j = 0; j < Dim; j++)
100			if (i == j)
101			m(i, j) = 1.0;
102			else
103			m(i, j) = 0.0;
104	tim	70
105	gezelter	507	return m;
106			}
107	tim	74
108	gezelter	507	/**
109			* Retunrs the inversion of this matrix.
110			* @todo need implementation
111			*/
112			SquareMatrix<Real, Dim> inverse() {
113			SquareMatrix<Real, Dim> result;
114	tim	70
115	gezelter	507	return result;
116			}
117	tim	70
118	gezelter	507	/**
119			* Returns the determinant of this matrix.
120			* @todo need implementation
121			*/
122			Real determinant() const {
123			Real det;
124			return det;
125			}
126	tim	70
127	gezelter	507	/** Returns the trace of this matrix. */
128			Real trace() const {
129			Real tmp = 0;
130	tim	137
131	gezelter	507	for (unsigned int i = 0; i < Dim ; i++)
132			tmp += this->data_[i][i];
133	tim	70
134	gezelter	507	return tmp;
135			}
136	tim	70
137	gezelter	507	/** Tests if this matrix is symmetrix. */
138			bool isSymmetric() const {
139			for (unsigned int i = 0; i < Dim - 1; i++)
140			for (unsigned int j = i; j < Dim; j++)
141			if (fabs(this->data_[i][j] - this->data_[j][i]) > oopse::epsilon)
142			return false;
143	tim	137
144	gezelter	507	return true;
145			}
146	tim	70
147	gezelter	507	/** Tests if this matrix is orthogonal. */
148			bool isOrthogonal() {
149			SquareMatrix<Real, Dim> tmp;
150	tim	70
151	gezelter	507	tmp = this transpose();
152	tim	70
153	gezelter	507	return tmp.isDiagonal();
154			}
155	tim	70
156	gezelter	507	/** Tests if this matrix is diagonal. */
157			bool isDiagonal() const {
158			for (unsigned int i = 0; i < Dim ; i++)
159			for (unsigned int j = 0; j < Dim; j++)
160			if (i !=j && fabs(this->data_[i][j]) > oopse::epsilon)
161			return false;
162	tim	137
163	gezelter	507	return true;
164			}
165	tim	137
166	gezelter	507	/** Tests if this matrix is the unit matrix. */
167			bool isUnitMatrix() const {
168			if (!isDiagonal())
169			return false;
170	tim	70
171	gezelter	507	for (unsigned int i = 0; i < Dim ; i++)
172			if (fabs(this->data_[i][i] - 1) > oopse::epsilon)
173			return false;
174	tim	137
175	gezelter	507	return true;
176			}
177	tim	70
178	gezelter	507	/** Return the transpose of this matrix */
179			SquareMatrix<Real, Dim> transpose() const{
180			SquareMatrix<Real, Dim> result;
181	tim	273
182	gezelter	507	for (unsigned int i = 0; i < Dim; i++)
183			for (unsigned int j = 0; j < Dim; j++)
184			result(j, i) = this->data_[i][j];
185	tim	273
186	gezelter	507	return result;
187			}
188	tim	273
189	gezelter	507	/** @todo need implementation */
190			void diagonalize() {
191			//jacobi(m, eigenValues, ortMat);
192			}
193	tim	76
194	gezelter	507	/**
195			* Jacobi iteration routines for computing eigenvalues/eigenvectors of
196			* real symmetric matrix
197			*
198			* @return true if success, otherwise return false
199			* @param a symmetric matrix whose eigenvectors are to be computed. On return, the matrix is
200			* overwritten
201			* @param w will contain the eigenvalues of the matrix On return of this function
202			* @param v the columns of this matrix will contain the eigenvectors. The eigenvectors are
203			* normalized and mutually orthogonal.
204			*/
205	tim	137
206	gezelter	507	static int jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& d,
207			SquareMatrix<Real, Dim>& v);
208			};//end SquareMatrix
209	tim	70
210	tim	76
211	gezelter	507	/*=========================================================================
212	tim	76
213	tim	123	Program: Visualization Toolkit
214			Module: $RCSfile: SquareMatrix.hpp,v $
215	tim	76
216	tim	123	Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
217			All rights reserved.
218			See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
219
220	gezelter	507	This software is distributed WITHOUT ANY WARRANTY; without even
221			the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
222			PURPOSE. See the above copyright notice for more information.
223	tim	123
224	gezelter	507	=========================================================================*/
225	tim	123
226	gezelter	507	#define VTK_ROTATE(a,i,j,k,l) g=a(i, j);h=a(k, l);a(i, j)=g-s(h+gtau); \
227			a(k, l)=h+s(g-htau)
228	tim	123
229			#define VTK_MAX_ROTATIONS 20
230
231	gezelter	507	// Jacobi iteration for the solution of eigenvectors/eigenvalues of a nxn
232			// real symmetric matrix. Square nxn matrix a; size of matrix in n;
233			// output eigenvalues in w; and output eigenvectors in v. Resulting
234			// eigenvalues/vectors are sorted in decreasing order; eigenvectors are
235			// normalized.
236			template<typename Real, int Dim>
237			int SquareMatrix<Real, Dim>::jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
238			SquareMatrix<Real, Dim>& v) {
239			const int n = Dim;
240			int i, j, k, iq, ip, numPos;
241			Real tresh, theta, tau, t, sm, s, h, g, c, tmp;
242			Real bspace[4], zspace[4];
243			Real *b = bspace;
244			Real *z = zspace;
245	tim	123
246	gezelter	507	// only allocate memory if the matrix is large
247			if (n > 4) {
248			b = new Real[n];
249			z = new Real[n];
250			}
251	tim	123
252	gezelter	507	// initialize
253			for (ip=0; ip<n; ip++) {
254			for (iq=0; iq<n; iq++) {
255			v(ip, iq) = 0.0;
256			}
257			v(ip, ip) = 1.0;
258			}
259			for (ip=0; ip<n; ip++) {
260			b[ip] = w[ip] = a(ip, ip);
261			z[ip] = 0.0;
262			}
263	tim	76
264	gezelter	507	// begin rotation sequence
265			for (i=0; i<VTK_MAX_ROTATIONS; i++) {
266			sm = 0.0;
267			for (ip=0; ip<n-1; ip++) {
268			for (iq=ip+1; iq<n; iq++) {
269			sm += fabs(a(ip, iq));
270			}
271			}
272			if (sm == 0.0) {
273			break;
274			}
275	tim	76
276	gezelter	507	if (i < 3) { // first 3 sweeps
277			tresh = 0.2sm/(nn);
278			} else {
279			tresh = 0.0;
280			}
281	tim	76
282	gezelter	507	for (ip=0; ip<n-1; ip++) {
283			for (iq=ip+1; iq<n; iq++) {
284			g = 100.0*fabs(a(ip, iq));
285	tim	76
286	gezelter	507	// after 4 sweeps
287			if (i > 3 && (fabs(w[ip])+g) == fabs(w[ip])
288			&& (fabs(w[iq])+g) == fabs(w[iq])) {
289			a(ip, iq) = 0.0;
290			} else if (fabs(a(ip, iq)) > tresh) {
291			h = w[iq] - w[ip];
292			if ( (fabs(h)+g) == fabs(h)) {
293			t = (a(ip, iq)) / h;
294			} else {
295			theta = 0.5*h / (a(ip, iq));
296			t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
297			if (theta < 0.0) {
298			t = -t;
299			}
300			}
301			c = 1.0 / sqrt(1+t*t);
302			s = t*c;
303			tau = s/(1.0+c);
304			h = t*a(ip, iq);
305			z[ip] -= h;
306			z[iq] += h;
307			w[ip] -= h;
308			w[iq] += h;
309			a(ip, iq)=0.0;
310	tim	76
311	gezelter	507	// ip already shifted left by 1 unit
312			for (j = 0;j <= ip-1;j++) {
313			VTK_ROTATE(a,j,ip,j,iq);
314			}
315			// ip and iq already shifted left by 1 unit
316			for (j = ip+1;j <= iq-1;j++) {
317			VTK_ROTATE(a,ip,j,j,iq);
318			}
319			// iq already shifted left by 1 unit
320			for (j=iq+1; j<n; j++) {
321			VTK_ROTATE(a,ip,j,iq,j);
322			}
323			for (j=0; j<n; j++) {
324			VTK_ROTATE(v,j,ip,j,iq);
325			}
326			}
327			}
328			}
329	tim	93
330	gezelter	507	for (ip=0; ip<n; ip++) {
331			b[ip] += z[ip];
332			w[ip] = b[ip];
333			z[ip] = 0.0;
334			}
335			}
336	tim	93
337	gezelter	507	//// this is NEVER called
338			if ( i >= VTK_MAX_ROTATIONS ) {
339			std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions" << std::endl;
340			return 0;
341			}
342	tim	76
343	gezelter	507	// sort eigenfunctions these changes do not affect accuracy
344			for (j=0; j<n-1; j++) { // boundary incorrect
345			k = j;
346			tmp = w[k];
347			for (i=j+1; i<n; i++) { // boundary incorrect, shifted already
348			if (w[i] >= tmp) { // why exchage if same?
349			k = i;
350			tmp = w[k];
351			}
352			}
353			if (k != j) {
354			w[k] = w[j];
355			w[j] = tmp;
356			for (i=0; i<n; i++) {
357			tmp = v(i, j);
358			v(i, j) = v(i, k);
359			v(i, k) = tmp;
360			}
361			}
362			}
363			// insure eigenvector consistency (i.e., Jacobi can compute vectors that
364			// are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
365			// reek havoc in hyperstreamline/other stuff. We will select the most
366			// positive eigenvector.
367			int ceil_half_n = (n >> 1) + (n & 1);
368			for (j=0; j<n; j++) {
369			for (numPos=0, i=0; i<n; i++) {
370			if ( v(i, j) >= 0.0 ) {
371			numPos++;
372			}
373			}
374			// if ( numPos < ceil(double(n)/double(2.0)) )
375			if ( numPos < ceil_half_n) {
376			for (i=0; i<n; i++) {
377			v(i, j) *= -1.0;
378			}
379			}
380			}
381	tim	76
382	gezelter	507	if (n > 4) {
383			delete [] b;
384			delete [] z;
385	tim	76	}
386	gezelter	507	return 1;
387			}
388	tim	76
389
390			}
391	tim	123	#endif //MATH_SQUAREMATRIX_HPP
392	tim	76