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:	2069
Committed:	Tue Mar 1 20:10:14 2005 UTC (19 years, 4 months ago) by tim
File size:	14108 byte(s)
Log Message:	fix compilation problem for g++ 3.4
#	User	Rev	Content
1	gezelter	1930	/*
2			* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	tim	1563	*
4	gezelter	1930	* 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	1563	*/
41	gezelter	1930
42	tim	1563	/**
43			* @file SquareMatrix.hpp
44			* @author Teng Lin
45			* @date 10/11/2004
46			* @version 1.0
47			*/
48	tim	1616	#ifndef MATH_SQUAREMATRIX_HPP
49	tim	1563	#define MATH_SQUAREMATRIX_HPP
50
51	tim	1567	#include "math/RectMatrix.hpp"
52	tim	1563
53			namespace oopse {
54
55			/**
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	tim	1567	class SquareMatrix : public RectMatrix<Real, Dim, Dim> {
63	tim	1563	public:
64	tim	1630	typedef Real ElemType;
65			typedef Real* ElemPoinerType;
66	tim	1563
67	tim	1630	/** default constructor */
68			SquareMatrix() {
69			for (unsigned int i = 0; i < Dim; i++)
70			for (unsigned int j = 0; j < Dim; j++)
71	tim	2069	this->data_[i][j] = 0.0;
72	tim	1630	}
73	tim	1563
74	tim	1644	/** Constructs and initializes every element of this matrix to a scalar */
75			SquareMatrix(Real s) : RectMatrix<Real, Dim, Dim>(s){
76			}
77
78			/** Constructs and initializes from an array */
79			SquareMatrix(Real* array) : RectMatrix<Real, Dim, Dim>(array){
80			}
81
82
83	tim	1630	/** copy constructor */
84			SquareMatrix(const RectMatrix<Real, Dim, Dim>& m) : RectMatrix<Real, Dim, Dim>(m) {
85			}
86	tim	1563
87	tim	1630	/** 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
93			/** Retunrs an identity matrix*/
94	tim	1567
95	tim	1630	static SquareMatrix<Real, Dim> identity() {
96			SquareMatrix<Real, Dim> m;
97
98			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	1563
105	tim	1630	return m;
106			}
107	tim	1567
108	tim	1630	/**
109			* Retunrs the inversion of this matrix.
110			* @todo need implementation
111			*/
112			SquareMatrix<Real, Dim> inverse() {
113			SquareMatrix<Real, Dim> result;
114	tim	1563
115	tim	1630	return result;
116			}
117	tim	1563
118	tim	1630	/**
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	1563
127	tim	1630	/** Returns the trace of this matrix. */
128			Real trace() const {
129			Real tmp = 0;
130
131			for (unsigned int i = 0; i < Dim ; i++)
132	tim	2069	tmp += this->data_[i][i];
133	tim	1563
134	tim	1630	return tmp;
135			}
136	tim	1563
137	tim	1630	/** 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	tim	2069	if (fabs(this->data_[i][j] - this->data_[j][i]) > oopse::epsilon)
142	tim	1630	return false;
143
144			return true;
145			}
146	tim	1563
147	tim	1630	/** Tests if this matrix is orthogonal. */
148			bool isOrthogonal() {
149			SquareMatrix<Real, Dim> tmp;
150	tim	1563
151	tim	1630	tmp = this transpose();
152	tim	1563
153	tim	1630	return tmp.isDiagonal();
154			}
155	tim	1563
156	tim	1630	/** 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	tim	2069	if (i !=j && fabs(this->data_[i][j]) > oopse::epsilon)
161	tim	1630	return false;
162
163			return true;
164			}
165
166			/** Tests if this matrix is the unit matrix. */
167			bool isUnitMatrix() const {
168			if (!isDiagonal())
169	tim	1563	return false;
170
171	tim	1630	for (unsigned int i = 0; i < Dim ; i++)
172	tim	2069	if (fabs(this->data_[i][i] - 1) > oopse::epsilon)
173	tim	1630	return false;
174
175			return true;
176			}
177	tim	1563
178	tim	1957	/** Return the transpose of this matrix */
179			SquareMatrix<Real, Dim> transpose() const{
180			SquareMatrix<Real, Dim> result;
181
182			for (unsigned int i = 0; i < Dim; i++)
183			for (unsigned int j = 0; j < Dim; j++)
184	tim	2069	result(j, i) = this->data_[i][j];
185	tim	1957
186			return result;
187			}
188
189	tim	1630	/** @todo need implementation */
190			void diagonalize() {
191			//jacobi(m, eigenValues, ortMat);
192			}
193	tim	1569
194	tim	1630	/**
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
206			static int jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& d,
207			SquareMatrix<Real, Dim>& v);
208	tim	1563	};//end SquareMatrix
209
210	tim	1569
211	tim	1616	/*=========================================================================
212	tim	1569
213	tim	1616	Program: Visualization Toolkit
214			Module: $RCSfile: SquareMatrix.hpp,v $
215	tim	1569
216	tim	1616	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			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
224			=========================================================================*/
225
226			#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
229			#define VTK_MAX_ROTATIONS 20
230
231			// 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	tim	1639	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	1616
246	tim	1639	// only allocate memory if the matrix is large
247			if (n > 4) {
248			b = new Real[n];
249			z = new Real[n];
250	tim	1616	}
251
252	tim	1639	// 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	tim	1616	}
259	tim	1639	for (ip=0; ip<n; ip++) {
260			b[ip] = w[ip] = a(ip, ip);
261			z[ip] = 0.0;
262	tim	1616	}
263	tim	1569
264	tim	1639	// 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	tim	1616	}
272	tim	1639	if (sm == 0.0) {
273			break;
274			}
275	tim	1569
276	tim	1639	if (i < 3) { // first 3 sweeps
277			tresh = 0.2sm/(nn);
278			} else {
279			tresh = 0.0;
280			}
281	tim	1569
282	tim	1639	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	1569
286	tim	1639	// 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	1569
311	tim	1639	// 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	tim	1586	}
328			}
329
330	tim	1639	for (ip=0; ip<n; ip++) {
331			b[ip] += z[ip];
332			w[ip] = b[ip];
333			z[ip] = 0.0;
334			}
335	tim	1586	}
336
337	tim	1639	//// this is NEVER called
338			if ( i >= VTK_MAX_ROTATIONS ) {
339			std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions" << std::endl;
340			return 0;
341	tim	1616	}
342	tim	1569
343	tim	1639	// sort eigenfunctions these changes do not affect accuracy
344			for (j=0; j<n-1; j++) { // boundary incorrect
345			k = j;
346	tim	1616	tmp = w[k];
347	tim	1639	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	tim	1586	}
353	tim	1639	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	tim	1616	}
362	tim	1586	}
363	tim	1639	// 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	tim	1586	}
374	tim	1639	// 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	tim	1616	}
380	tim	1586	}
381	tim	1569
382	tim	1639	if (n > 4) {
383			delete [] b;
384			delete [] z;
385	tim	1616	}
386	tim	1639	return 1;
387	tim	1569	}
388
389
390			}
391	tim	1616	#endif //MATH_SQUAREMATRIX_HPP
392	tim	1569