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;
            }         

            /** 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) = 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:	273
Committed:	Tue Jan 25 17:45:23 2005 UTC (20 years, 9 months ago) by tim
File size:	14066 byte(s)
Log Message:	(1) complete section parser's error message (2) add GhostTorsion (3) accumulate inertial tensor from the directional atoms before calculate rigidbody's inertial tensor
#	Content
1	/*
2	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	*
4	* 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	*/
41
42	/**
43	* @file SquareMatrix.hpp
44	* @author Teng Lin
45	* @date 10/11/2004
46	* @version 1.0
47	*/
48	#ifndef MATH_SQUAREMATRIX_HPP
49	#define MATH_SQUAREMATRIX_HPP
50
51	#include "math/RectMatrix.hpp"
52
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	class SquareMatrix : public RectMatrix<Real, Dim, Dim> {
63	public:
64	typedef Real ElemType;
65	typedef Real* ElemPoinerType;
66
67	/** 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
74	/** 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	/** copy constructor */
84	SquareMatrix(const RectMatrix<Real, Dim, Dim>& m) : RectMatrix<Real, Dim, Dim>(m) {
85	}
86
87	/** 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
95	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
105	return m;
106	}
107
108	/**
109	* Retunrs the inversion of this matrix.
110	* @todo need implementation
111	*/
112	SquareMatrix<Real, Dim> inverse() {
113	SquareMatrix<Real, Dim> result;
114
115	return result;
116	}
117
118	/**
119	* Returns the determinant of this matrix.
120	* @todo need implementation
121	*/
122	Real determinant() const {
123	Real det;
124	return det;
125	}
126
127	/** 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
134	return tmp;
135	}
136
137	/** 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
147	/** Tests if this matrix is orthogonal. */
148	bool isOrthogonal() {
149	SquareMatrix<Real, Dim> tmp;
150
151	tmp = this transpose();
152
153	return tmp.isDiagonal();
154	}
155
156	/** 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	return false;
170
171	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
178	/** 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	result(j, i) = data_[i][j];
185
186	return result;
187	}
188
189	/** @todo need implementation */
190	void diagonalize() {
191	//jacobi(m, eigenValues, ortMat);
192	}
193
194	/**
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	};//end SquareMatrix
209
210
211	/*=========================================================================
212
213	Program: Visualization Toolkit
214	Module: $RCSfile: SquareMatrix.hpp,v $
215
216	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	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
246	// only allocate memory if the matrix is large
247	if (n > 4) {
248	b = new Real[n];
249	z = new Real[n];
250	}
251
252	// 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
264	// 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
276	if (i < 3) { // first 3 sweeps
277	tresh = 0.2sm/(nn);
278	} else {
279	tresh = 0.0;
280	}
281
282	for (ip=0; ip<n-1; ip++) {
283	for (iq=ip+1; iq<n; iq++) {
284	g = 100.0*fabs(a(ip, iq));
285
286	// 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
311	// 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
330	for (ip=0; ip<n; ip++) {
331	b[ip] += z[ip];
332	w[ip] = b[ip];
333	z[ip] = 0.0;
334	}
335	}
336
337	//// 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
343	// 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
382	if (n > 4) {
383	delete [] b;
384	delete [] z;
385	}
386	return 1;
387	}
388
389
390	}
391	#endif //MATH_SQUAREMATRIX_HPP
392