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++)
                        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 += 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(data_[i][j] - 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(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(data_[i][i] - 1) > oopse::epsilon)
                        return false;
                    
                return true;
            }         

            /** @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:	1930
Committed:	Wed Jan 12 22:41:40 2005 UTC (19 years, 5 months ago) by gezelter
File size:	13650 byte(s)
Log Message:	merging new_design branch into OOPSE-2.0
#	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			data_[i][j] = 0.0;
72			}
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			tmp += 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			if (fabs(data_[i][j] - data_[j][i]) > oopse::epsilon)
142			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			if (i !=j && fabs(data_[i][j]) > oopse::epsilon)
161			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			if (fabs(data_[i][i] - 1) > oopse::epsilon)
173			return false;
174
175			return true;
176			}
177	tim	1563
178	tim	1630	/** @todo need implementation */
179			void diagonalize() {
180			//jacobi(m, eigenValues, ortMat);
181			}
182	tim	1569
183	tim	1630	/**
184			* Jacobi iteration routines for computing eigenvalues/eigenvectors of
185			* real symmetric matrix
186			*
187			* @return true if success, otherwise return false
188			* @param a symmetric matrix whose eigenvectors are to be computed. On return, the matrix is
189			* overwritten
190			* @param w will contain the eigenvalues of the matrix On return of this function
191			* @param v the columns of this matrix will contain the eigenvectors. The eigenvectors are
192			* normalized and mutually orthogonal.
193			*/
194
195			static int jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& d,
196			SquareMatrix<Real, Dim>& v);
197	tim	1563	};//end SquareMatrix
198
199	tim	1569
200	tim	1616	/*=========================================================================
201	tim	1569
202	tim	1616	Program: Visualization Toolkit
203			Module: $RCSfile: SquareMatrix.hpp,v $
204	tim	1569
205	tim	1616	Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
206			All rights reserved.
207			See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
208
209			This software is distributed WITHOUT ANY WARRANTY; without even
210			the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
211			PURPOSE. See the above copyright notice for more information.
212
213			=========================================================================*/
214
215			#define VTK_ROTATE(a,i,j,k,l) g=a(i, j);h=a(k, l);a(i, j)=g-s(h+gtau);\
216			a(k, l)=h+s(g-htau)
217
218			#define VTK_MAX_ROTATIONS 20
219
220			// Jacobi iteration for the solution of eigenvectors/eigenvalues of a nxn
221			// real symmetric matrix. Square nxn matrix a; size of matrix in n;
222			// output eigenvalues in w; and output eigenvectors in v. Resulting
223			// eigenvalues/vectors are sorted in decreasing order; eigenvectors are
224			// normalized.
225			template<typename Real, int Dim>
226			int SquareMatrix<Real, Dim>::jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
227	tim	1639	SquareMatrix<Real, Dim>& v) {
228			const int n = Dim;
229			int i, j, k, iq, ip, numPos;
230			Real tresh, theta, tau, t, sm, s, h, g, c, tmp;
231			Real bspace[4], zspace[4];
232			Real *b = bspace;
233			Real *z = zspace;
234	tim	1616
235	tim	1639	// only allocate memory if the matrix is large
236			if (n > 4) {
237			b = new Real[n];
238			z = new Real[n];
239	tim	1616	}
240
241	tim	1639	// initialize
242			for (ip=0; ip<n; ip++) {
243			for (iq=0; iq<n; iq++) {
244			v(ip, iq) = 0.0;
245			}
246			v(ip, ip) = 1.0;
247	tim	1616	}
248	tim	1639	for (ip=0; ip<n; ip++) {
249			b[ip] = w[ip] = a(ip, ip);
250			z[ip] = 0.0;
251	tim	1616	}
252	tim	1569
253	tim	1639	// begin rotation sequence
254			for (i=0; i<VTK_MAX_ROTATIONS; i++) {
255			sm = 0.0;
256			for (ip=0; ip<n-1; ip++) {
257			for (iq=ip+1; iq<n; iq++) {
258			sm += fabs(a(ip, iq));
259			}
260	tim	1616	}
261	tim	1639	if (sm == 0.0) {
262			break;
263			}
264	tim	1569
265	tim	1639	if (i < 3) { // first 3 sweeps
266			tresh = 0.2sm/(nn);
267			} else {
268			tresh = 0.0;
269			}
270	tim	1569
271	tim	1639	for (ip=0; ip<n-1; ip++) {
272			for (iq=ip+1; iq<n; iq++) {
273			g = 100.0*fabs(a(ip, iq));
274	tim	1569
275	tim	1639	// after 4 sweeps
276			if (i > 3 && (fabs(w[ip])+g) == fabs(w[ip])
277			&& (fabs(w[iq])+g) == fabs(w[iq])) {
278			a(ip, iq) = 0.0;
279			} else if (fabs(a(ip, iq)) > tresh) {
280			h = w[iq] - w[ip];
281			if ( (fabs(h)+g) == fabs(h)) {
282			t = (a(ip, iq)) / h;
283			} else {
284			theta = 0.5*h / (a(ip, iq));
285			t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
286			if (theta < 0.0) {
287			t = -t;
288			}
289			}
290			c = 1.0 / sqrt(1+t*t);
291			s = t*c;
292			tau = s/(1.0+c);
293			h = t*a(ip, iq);
294			z[ip] -= h;
295			z[iq] += h;
296			w[ip] -= h;
297			w[iq] += h;
298			a(ip, iq)=0.0;
299	tim	1569
300	tim	1639	// ip already shifted left by 1 unit
301			for (j = 0;j <= ip-1;j++) {
302			VTK_ROTATE(a,j,ip,j,iq);
303			}
304			// ip and iq already shifted left by 1 unit
305			for (j = ip+1;j <= iq-1;j++) {
306			VTK_ROTATE(a,ip,j,j,iq);
307			}
308			// iq already shifted left by 1 unit
309			for (j=iq+1; j<n; j++) {
310			VTK_ROTATE(a,ip,j,iq,j);
311			}
312			for (j=0; j<n; j++) {
313			VTK_ROTATE(v,j,ip,j,iq);
314			}
315			}
316	tim	1586	}
317			}
318
319	tim	1639	for (ip=0; ip<n; ip++) {
320			b[ip] += z[ip];
321			w[ip] = b[ip];
322			z[ip] = 0.0;
323			}
324	tim	1586	}
325
326	tim	1639	//// this is NEVER called
327			if ( i >= VTK_MAX_ROTATIONS ) {
328			std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions" << std::endl;
329			return 0;
330	tim	1616	}
331	tim	1569
332	tim	1639	// sort eigenfunctions these changes do not affect accuracy
333			for (j=0; j<n-1; j++) { // boundary incorrect
334			k = j;
335	tim	1616	tmp = w[k];
336	tim	1639	for (i=j+1; i<n; i++) { // boundary incorrect, shifted already
337			if (w[i] >= tmp) { // why exchage if same?
338			k = i;
339			tmp = w[k];
340			}
341	tim	1586	}
342	tim	1639	if (k != j) {
343			w[k] = w[j];
344			w[j] = tmp;
345			for (i=0; i<n; i++) {
346			tmp = v(i, j);
347			v(i, j) = v(i, k);
348			v(i, k) = tmp;
349			}
350	tim	1616	}
351	tim	1586	}
352	tim	1639	// insure eigenvector consistency (i.e., Jacobi can compute vectors that
353			// are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
354			// reek havoc in hyperstreamline/other stuff. We will select the most
355			// positive eigenvector.
356			int ceil_half_n = (n >> 1) + (n & 1);
357			for (j=0; j<n; j++) {
358			for (numPos=0, i=0; i<n; i++) {
359			if ( v(i, j) >= 0.0 ) {
360			numPos++;
361			}
362	tim	1586	}
363	tim	1639	// if ( numPos < ceil(double(n)/double(2.0)) )
364			if ( numPos < ceil_half_n) {
365			for (i=0; i<n; i++) {
366			v(i, j) *= -1.0;
367			}
368	tim	1616	}
369	tim	1586	}
370	tim	1569
371	tim	1639	if (n > 4) {
372			delete [] b;
373			delete [] z;
374	tim	1616	}
375	tim	1639	return 1;
376	tim	1569	}
377
378
379			}
380	tim	1616	#endif //MATH_SQUAREMATRIX_HPP
381	tim	1569