src/math/SquareMatrix.hpp

/*
 * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
 * 
 * Contact: oopse@oopse.org
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */

/**
 * @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:	1644
Committed:	Mon Oct 25 22:46:19 2004 UTC (19 years, 8 months ago) by tim
File size:	12720 byte(s)
Log Message:	add getArray function to RectMatrix and Vector classes
#	User	Rev	Content
1	tim	1563	/*
2			* Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project
3			*
4			* Contact: oopse@oopse.org
5			*
6			* This program is free software; you can redistribute it and/or
7			* modify it under the terms of the GNU Lesser General Public License
8			* as published by the Free Software Foundation; either version 2.1
9			* of the License, or (at your option) any later version.
10			* All we ask is that proper credit is given for our work, which includes
11			* - but is not limited to - adding the above copyright notice to the beginning
12			* of your source code files, and to any copyright notice that you may distribute
13			* with programs based on this work.
14			*
15			* This program is distributed in the hope that it will be useful,
16			* but WITHOUT ANY WARRANTY; without even the implied warranty of
17			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18			* GNU Lesser General Public License for more details.
19			*
20			* You should have received a copy of the GNU Lesser General Public License
21			* along with this program; if not, write to the Free Software
22			* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
23			*
24			*/
25
26			/**
27			* @file SquareMatrix.hpp
28			* @author Teng Lin
29			* @date 10/11/2004
30			* @version 1.0
31			*/
32	tim	1616	#ifndef MATH_SQUAREMATRIX_HPP
33	tim	1563	#define MATH_SQUAREMATRIX_HPP
34
35	tim	1567	#include "math/RectMatrix.hpp"
36	tim	1563
37			namespace oopse {
38
39			/**
40			* @class SquareMatrix SquareMatrix.hpp "math/SquareMatrix.hpp"
41			* @brief A square matrix class
42			* @template Real the element type
43			* @template Dim the dimension of the square matrix
44			*/
45			template<typename Real, int Dim>
46	tim	1567	class SquareMatrix : public RectMatrix<Real, Dim, Dim> {
47	tim	1563	public:
48	tim	1630	typedef Real ElemType;
49			typedef Real* ElemPoinerType;
50	tim	1563
51	tim	1630	/** default constructor */
52			SquareMatrix() {
53			for (unsigned int i = 0; i < Dim; i++)
54			for (unsigned int j = 0; j < Dim; j++)
55			data_[i][j] = 0.0;
56			}
57	tim	1563
58	tim	1644	/** Constructs and initializes every element of this matrix to a scalar */
59			SquareMatrix(Real s) : RectMatrix<Real, Dim, Dim>(s){
60			}
61
62			/** Constructs and initializes from an array */
63			SquareMatrix(Real* array) : RectMatrix<Real, Dim, Dim>(array){
64			}
65
66
67	tim	1630	/** copy constructor */
68			SquareMatrix(const RectMatrix<Real, Dim, Dim>& m) : RectMatrix<Real, Dim, Dim>(m) {
69			}
70	tim	1563
71	tim	1630	/** copy assignment operator */
72			SquareMatrix<Real, Dim>& operator =(const RectMatrix<Real, Dim, Dim>& m) {
73			RectMatrix<Real, Dim, Dim>::operator=(m);
74			return *this;
75			}
76
77			/** Retunrs an identity matrix*/
78	tim	1567
79	tim	1630	static SquareMatrix<Real, Dim> identity() {
80			SquareMatrix<Real, Dim> m;
81
82			for (unsigned int i = 0; i < Dim; i++)
83			for (unsigned int j = 0; j < Dim; j++)
84			if (i == j)
85			m(i, j) = 1.0;
86			else
87			m(i, j) = 0.0;
88	tim	1563
89	tim	1630	return m;
90			}
91	tim	1567
92	tim	1630	/**
93			* Retunrs the inversion of this matrix.
94			* @todo need implementation
95			*/
96			SquareMatrix<Real, Dim> inverse() {
97			SquareMatrix<Real, Dim> result;
98	tim	1563
99	tim	1630	return result;
100			}
101	tim	1563
102	tim	1630	/**
103			* Returns the determinant of this matrix.
104			* @todo need implementation
105			*/
106			Real determinant() const {
107			Real det;
108			return det;
109			}
110	tim	1563
111	tim	1630	/** Returns the trace of this matrix. */
112			Real trace() const {
113			Real tmp = 0;
114
115			for (unsigned int i = 0; i < Dim ; i++)
116			tmp += data_[i][i];
117	tim	1563
118	tim	1630	return tmp;
119			}
120	tim	1563
121	tim	1630	/** Tests if this matrix is symmetrix. */
122			bool isSymmetric() const {
123			for (unsigned int i = 0; i < Dim - 1; i++)
124			for (unsigned int j = i; j < Dim; j++)
125			if (fabs(data_[i][j] - data_[j][i]) > oopse::epsilon)
126			return false;
127
128			return true;
129			}
130	tim	1563
131	tim	1630	/** Tests if this matrix is orthogonal. */
132			bool isOrthogonal() {
133			SquareMatrix<Real, Dim> tmp;
134	tim	1563
135	tim	1630	tmp = this transpose();
136	tim	1563
137	tim	1630	return tmp.isDiagonal();
138			}
139	tim	1563
140	tim	1630	/** Tests if this matrix is diagonal. */
141			bool isDiagonal() const {
142			for (unsigned int i = 0; i < Dim ; i++)
143			for (unsigned int j = 0; j < Dim; j++)
144			if (i !=j && fabs(data_[i][j]) > oopse::epsilon)
145			return false;
146
147			return true;
148			}
149
150			/** Tests if this matrix is the unit matrix. */
151			bool isUnitMatrix() const {
152			if (!isDiagonal())
153	tim	1563	return false;
154
155	tim	1630	for (unsigned int i = 0; i < Dim ; i++)
156			if (fabs(data_[i][i] - 1) > oopse::epsilon)
157			return false;
158
159			return true;
160			}
161	tim	1563
162	tim	1630	/** @todo need implementation */
163			void diagonalize() {
164			//jacobi(m, eigenValues, ortMat);
165			}
166	tim	1569
167	tim	1630	/**
168			* Jacobi iteration routines for computing eigenvalues/eigenvectors of
169			* real symmetric matrix
170			*
171			* @return true if success, otherwise return false
172			* @param a symmetric matrix whose eigenvectors are to be computed. On return, the matrix is
173			* overwritten
174			* @param w will contain the eigenvalues of the matrix On return of this function
175			* @param v the columns of this matrix will contain the eigenvectors. The eigenvectors are
176			* normalized and mutually orthogonal.
177			*/
178
179			static int jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& d,
180			SquareMatrix<Real, Dim>& v);
181	tim	1563	};//end SquareMatrix
182
183	tim	1569
184	tim	1616	/*=========================================================================
185	tim	1569
186	tim	1616	Program: Visualization Toolkit
187			Module: $RCSfile: SquareMatrix.hpp,v $
188	tim	1569
189	tim	1616	Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
190			All rights reserved.
191			See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
192
193			This software is distributed WITHOUT ANY WARRANTY; without even
194			the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
195			PURPOSE. See the above copyright notice for more information.
196
197			=========================================================================*/
198
199			#define VTK_ROTATE(a,i,j,k,l) g=a(i, j);h=a(k, l);a(i, j)=g-s(h+gtau);\
200			a(k, l)=h+s(g-htau)
201
202			#define VTK_MAX_ROTATIONS 20
203
204			// Jacobi iteration for the solution of eigenvectors/eigenvalues of a nxn
205			// real symmetric matrix. Square nxn matrix a; size of matrix in n;
206			// output eigenvalues in w; and output eigenvectors in v. Resulting
207			// eigenvalues/vectors are sorted in decreasing order; eigenvectors are
208			// normalized.
209			template<typename Real, int Dim>
210			int SquareMatrix<Real, Dim>::jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
211	tim	1639	SquareMatrix<Real, Dim>& v) {
212			const int n = Dim;
213			int i, j, k, iq, ip, numPos;
214			Real tresh, theta, tau, t, sm, s, h, g, c, tmp;
215			Real bspace[4], zspace[4];
216			Real *b = bspace;
217			Real *z = zspace;
218	tim	1616
219	tim	1639	// only allocate memory if the matrix is large
220			if (n > 4) {
221			b = new Real[n];
222			z = new Real[n];
223	tim	1616	}
224
225	tim	1639	// initialize
226			for (ip=0; ip<n; ip++) {
227			for (iq=0; iq<n; iq++) {
228			v(ip, iq) = 0.0;
229			}
230			v(ip, ip) = 1.0;
231	tim	1616	}
232	tim	1639	for (ip=0; ip<n; ip++) {
233			b[ip] = w[ip] = a(ip, ip);
234			z[ip] = 0.0;
235	tim	1616	}
236	tim	1569
237	tim	1639	// begin rotation sequence
238			for (i=0; i<VTK_MAX_ROTATIONS; i++) {
239			sm = 0.0;
240			for (ip=0; ip<n-1; ip++) {
241			for (iq=ip+1; iq<n; iq++) {
242			sm += fabs(a(ip, iq));
243			}
244	tim	1616	}
245	tim	1639	if (sm == 0.0) {
246			break;
247			}
248	tim	1569
249	tim	1639	if (i < 3) { // first 3 sweeps
250			tresh = 0.2sm/(nn);
251			} else {
252			tresh = 0.0;
253			}
254	tim	1569
255	tim	1639	for (ip=0; ip<n-1; ip++) {
256			for (iq=ip+1; iq<n; iq++) {
257			g = 100.0*fabs(a(ip, iq));
258	tim	1569
259	tim	1639	// after 4 sweeps
260			if (i > 3 && (fabs(w[ip])+g) == fabs(w[ip])
261			&& (fabs(w[iq])+g) == fabs(w[iq])) {
262			a(ip, iq) = 0.0;
263			} else if (fabs(a(ip, iq)) > tresh) {
264			h = w[iq] - w[ip];
265			if ( (fabs(h)+g) == fabs(h)) {
266			t = (a(ip, iq)) / h;
267			} else {
268			theta = 0.5*h / (a(ip, iq));
269			t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
270			if (theta < 0.0) {
271			t = -t;
272			}
273			}
274			c = 1.0 / sqrt(1+t*t);
275			s = t*c;
276			tau = s/(1.0+c);
277			h = t*a(ip, iq);
278			z[ip] -= h;
279			z[iq] += h;
280			w[ip] -= h;
281			w[iq] += h;
282			a(ip, iq)=0.0;
283	tim	1569
284	tim	1639	// ip already shifted left by 1 unit
285			for (j = 0;j <= ip-1;j++) {
286			VTK_ROTATE(a,j,ip,j,iq);
287			}
288			// ip and iq already shifted left by 1 unit
289			for (j = ip+1;j <= iq-1;j++) {
290			VTK_ROTATE(a,ip,j,j,iq);
291			}
292			// iq already shifted left by 1 unit
293			for (j=iq+1; j<n; j++) {
294			VTK_ROTATE(a,ip,j,iq,j);
295			}
296			for (j=0; j<n; j++) {
297			VTK_ROTATE(v,j,ip,j,iq);
298			}
299			}
300	tim	1586	}
301			}
302
303	tim	1639	for (ip=0; ip<n; ip++) {
304			b[ip] += z[ip];
305			w[ip] = b[ip];
306			z[ip] = 0.0;
307			}
308	tim	1586	}
309
310	tim	1639	//// this is NEVER called
311			if ( i >= VTK_MAX_ROTATIONS ) {
312			std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions" << std::endl;
313			return 0;
314	tim	1616	}
315	tim	1569
316	tim	1639	// sort eigenfunctions these changes do not affect accuracy
317			for (j=0; j<n-1; j++) { // boundary incorrect
318			k = j;
319	tim	1616	tmp = w[k];
320	tim	1639	for (i=j+1; i<n; i++) { // boundary incorrect, shifted already
321			if (w[i] >= tmp) { // why exchage if same?
322			k = i;
323			tmp = w[k];
324			}
325	tim	1586	}
326	tim	1639	if (k != j) {
327			w[k] = w[j];
328			w[j] = tmp;
329			for (i=0; i<n; i++) {
330			tmp = v(i, j);
331			v(i, j) = v(i, k);
332			v(i, k) = tmp;
333			}
334	tim	1616	}
335	tim	1586	}
336	tim	1639	// insure eigenvector consistency (i.e., Jacobi can compute vectors that
337			// are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
338			// reek havoc in hyperstreamline/other stuff. We will select the most
339			// positive eigenvector.
340			int ceil_half_n = (n >> 1) + (n & 1);
341			for (j=0; j<n; j++) {
342			for (numPos=0, i=0; i<n; i++) {
343			if ( v(i, j) >= 0.0 ) {
344			numPos++;
345			}
346	tim	1586	}
347	tim	1639	// if ( numPos < ceil(double(n)/double(2.0)) )
348			if ( numPos < ceil_half_n) {
349			for (i=0; i<n; i++) {
350			v(i, j) *= -1.0;
351			}
352	tim	1616	}
353	tim	1586	}
354	tim	1569
355	tim	1639	if (n > 4) {
356			delete [] b;
357			delete [] z;
358	tim	1616	}
359	tim	1639	return 1;
360	tim	1569	}
361
362
363			}
364	tim	1616	#endif //MATH_SQUAREMATRIX_HPP
365	tim	1569