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

            /** 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:	146
Committed:	Fri Oct 22 23:09:57 2004 UTC (20 years, 6 months ago) by tim
Original Path:	*trunk/src/math/SquareMatrix.hpp*
File size:	12399 byte(s)
Log Message:	more work in Snapshot
#	User	Rev	Content
1	tim	70	/*
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	123	#ifndef MATH_SQUAREMATRIX_HPP
33	tim	70	#define MATH_SQUAREMATRIX_HPP
34
35	tim	74	#include "math/RectMatrix.hpp"
36	tim	70
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	74	class SquareMatrix : public RectMatrix<Real, Dim, Dim> {
47	tim	70	public:
48	tim	137	typedef Real ElemType;
49			typedef Real* ElemPoinerType;
50	tim	70
51	tim	137	/** 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	70
58	tim	137	/** copy constructor */
59			SquareMatrix(const RectMatrix<Real, Dim, Dim>& m) : RectMatrix<Real, Dim, Dim>(m) {
60			}
61	tim	70
62	tim	137	/** copy assignment operator */
63			SquareMatrix<Real, Dim>& operator =(const RectMatrix<Real, Dim, Dim>& m) {
64			RectMatrix<Real, Dim, Dim>::operator=(m);
65			return *this;
66			}
67
68			/** Retunrs an identity matrix*/
69	tim	74
70	tim	137	static SquareMatrix<Real, Dim> identity() {
71			SquareMatrix<Real, Dim> m;
72
73			for (unsigned int i = 0; i < Dim; i++)
74			for (unsigned int j = 0; j < Dim; j++)
75			if (i == j)
76			m(i, j) = 1.0;
77			else
78			m(i, j) = 0.0;
79	tim	70
80	tim	137	return m;
81			}
82	tim	74
83	tim	137	/**
84			* Retunrs the inversion of this matrix.
85			* @todo need implementation
86			*/
87			SquareMatrix<Real, Dim> inverse() {
88			SquareMatrix<Real, Dim> result;
89	tim	70
90	tim	137	return result;
91			}
92	tim	70
93	tim	137	/**
94			* Returns the determinant of this matrix.
95			* @todo need implementation
96			*/
97			Real determinant() const {
98			Real det;
99			return det;
100			}
101	tim	70
102	tim	137	/** Returns the trace of this matrix. */
103			Real trace() const {
104			Real tmp = 0;
105
106			for (unsigned int i = 0; i < Dim ; i++)
107			tmp += data_[i][i];
108	tim	70
109	tim	137	return tmp;
110			}
111	tim	70
112	tim	137	/** Tests if this matrix is symmetrix. */
113			bool isSymmetric() const {
114			for (unsigned int i = 0; i < Dim - 1; i++)
115			for (unsigned int j = i; j < Dim; j++)
116			if (fabs(data_[i][j] - data_[j][i]) > oopse::epsilon)
117			return false;
118
119			return true;
120			}
121	tim	70
122	tim	137	/** Tests if this matrix is orthogonal. */
123			bool isOrthogonal() {
124			SquareMatrix<Real, Dim> tmp;
125	tim	70
126	tim	137	tmp = this transpose();
127	tim	70
128	tim	137	return tmp.isDiagonal();
129			}
130	tim	70
131	tim	137	/** Tests if this matrix is diagonal. */
132			bool isDiagonal() const {
133			for (unsigned int i = 0; i < Dim ; i++)
134			for (unsigned int j = 0; j < Dim; j++)
135			if (i !=j && fabs(data_[i][j]) > oopse::epsilon)
136			return false;
137
138			return true;
139			}
140
141			/** Tests if this matrix is the unit matrix. */
142			bool isUnitMatrix() const {
143			if (!isDiagonal())
144	tim	70	return false;
145
146	tim	137	for (unsigned int i = 0; i < Dim ; i++)
147			if (fabs(data_[i][i] - 1) > oopse::epsilon)
148			return false;
149
150			return true;
151			}
152	tim	70
153	tim	137	/** @todo need implementation */
154			void diagonalize() {
155			//jacobi(m, eigenValues, ortMat);
156			}
157	tim	76
158	tim	137	/**
159			* Jacobi iteration routines for computing eigenvalues/eigenvectors of
160			* real symmetric matrix
161			*
162			* @return true if success, otherwise return false
163			* @param a symmetric matrix whose eigenvectors are to be computed. On return, the matrix is
164			* overwritten
165			* @param w will contain the eigenvalues of the matrix On return of this function
166			* @param v the columns of this matrix will contain the eigenvectors. The eigenvectors are
167			* normalized and mutually orthogonal.
168			*/
169
170			static int jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& d,
171			SquareMatrix<Real, Dim>& v);
172	tim	70	};//end SquareMatrix
173
174	tim	76
175	tim	123	/*=========================================================================
176	tim	76
177	tim	123	Program: Visualization Toolkit
178			Module: $RCSfile: SquareMatrix.hpp,v $
179	tim	76
180	tim	123	Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
181			All rights reserved.
182			See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
183
184			This software is distributed WITHOUT ANY WARRANTY; without even
185			the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
186			PURPOSE. See the above copyright notice for more information.
187
188			=========================================================================*/
189
190			#define VTK_ROTATE(a,i,j,k,l) g=a(i, j);h=a(k, l);a(i, j)=g-s(h+gtau);\
191			a(k, l)=h+s(g-htau)
192
193			#define VTK_MAX_ROTATIONS 20
194
195			// Jacobi iteration for the solution of eigenvectors/eigenvalues of a nxn
196			// real symmetric matrix. Square nxn matrix a; size of matrix in n;
197			// output eigenvalues in w; and output eigenvectors in v. Resulting
198			// eigenvalues/vectors are sorted in decreasing order; eigenvectors are
199			// normalized.
200			template<typename Real, int Dim>
201			int SquareMatrix<Real, Dim>::jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
202	tim	146	SquareMatrix<Real, Dim>& v) {
203			const int n = Dim;
204			int i, j, k, iq, ip, numPos;
205			Real tresh, theta, tau, t, sm, s, h, g, c, tmp;
206			Real bspace[4], zspace[4];
207			Real *b = bspace;
208			Real *z = zspace;
209	tim	123
210	tim	146	// only allocate memory if the matrix is large
211			if (n > 4) {
212			b = new Real[n];
213			z = new Real[n];
214	tim	123	}
215
216	tim	146	// initialize
217			for (ip=0; ip<n; ip++) {
218			for (iq=0; iq<n; iq++) {
219			v(ip, iq) = 0.0;
220			}
221			v(ip, ip) = 1.0;
222	tim	123	}
223	tim	146	for (ip=0; ip<n; ip++) {
224			b[ip] = w[ip] = a(ip, ip);
225			z[ip] = 0.0;
226	tim	123	}
227	tim	76
228	tim	146	// begin rotation sequence
229			for (i=0; i<VTK_MAX_ROTATIONS; i++) {
230			sm = 0.0;
231			for (ip=0; ip<n-1; ip++) {
232			for (iq=ip+1; iq<n; iq++) {
233			sm += fabs(a(ip, iq));
234			}
235	tim	123	}
236	tim	146	if (sm == 0.0) {
237			break;
238			}
239	tim	76
240	tim	146	if (i < 3) { // first 3 sweeps
241			tresh = 0.2sm/(nn);
242			} else {
243			tresh = 0.0;
244			}
245	tim	76
246	tim	146	for (ip=0; ip<n-1; ip++) {
247			for (iq=ip+1; iq<n; iq++) {
248			g = 100.0*fabs(a(ip, iq));
249	tim	76
250	tim	146	// after 4 sweeps
251			if (i > 3 && (fabs(w[ip])+g) == fabs(w[ip])
252			&& (fabs(w[iq])+g) == fabs(w[iq])) {
253			a(ip, iq) = 0.0;
254			} else if (fabs(a(ip, iq)) > tresh) {
255			h = w[iq] - w[ip];
256			if ( (fabs(h)+g) == fabs(h)) {
257			t = (a(ip, iq)) / h;
258			} else {
259			theta = 0.5*h / (a(ip, iq));
260			t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
261			if (theta < 0.0) {
262			t = -t;
263			}
264			}
265			c = 1.0 / sqrt(1+t*t);
266			s = t*c;
267			tau = s/(1.0+c);
268			h = t*a(ip, iq);
269			z[ip] -= h;
270			z[iq] += h;
271			w[ip] -= h;
272			w[iq] += h;
273			a(ip, iq)=0.0;
274	tim	76
275	tim	146	// ip already shifted left by 1 unit
276			for (j = 0;j <= ip-1;j++) {
277			VTK_ROTATE(a,j,ip,j,iq);
278			}
279			// ip and iq already shifted left by 1 unit
280			for (j = ip+1;j <= iq-1;j++) {
281			VTK_ROTATE(a,ip,j,j,iq);
282			}
283			// iq already shifted left by 1 unit
284			for (j=iq+1; j<n; j++) {
285			VTK_ROTATE(a,ip,j,iq,j);
286			}
287			for (j=0; j<n; j++) {
288			VTK_ROTATE(v,j,ip,j,iq);
289			}
290			}
291	tim	93	}
292			}
293
294	tim	146	for (ip=0; ip<n; ip++) {
295			b[ip] += z[ip];
296			w[ip] = b[ip];
297			z[ip] = 0.0;
298			}
299	tim	93	}
300
301	tim	146	//// this is NEVER called
302			if ( i >= VTK_MAX_ROTATIONS ) {
303			std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions" << std::endl;
304			return 0;
305	tim	123	}
306	tim	76
307	tim	146	// sort eigenfunctions these changes do not affect accuracy
308			for (j=0; j<n-1; j++) { // boundary incorrect
309			k = j;
310	tim	123	tmp = w[k];
311	tim	146	for (i=j+1; i<n; i++) { // boundary incorrect, shifted already
312			if (w[i] >= tmp) { // why exchage if same?
313			k = i;
314			tmp = w[k];
315			}
316	tim	93	}
317	tim	146	if (k != j) {
318			w[k] = w[j];
319			w[j] = tmp;
320			for (i=0; i<n; i++) {
321			tmp = v(i, j);
322			v(i, j) = v(i, k);
323			v(i, k) = tmp;
324			}
325	tim	123	}
326	tim	93	}
327	tim	146	// insure eigenvector consistency (i.e., Jacobi can compute vectors that
328			// are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
329			// reek havoc in hyperstreamline/other stuff. We will select the most
330			// positive eigenvector.
331			int ceil_half_n = (n >> 1) + (n & 1);
332			for (j=0; j<n; j++) {
333			for (numPos=0, i=0; i<n; i++) {
334			if ( v(i, j) >= 0.0 ) {
335			numPos++;
336			}
337	tim	93	}
338	tim	146	// if ( numPos < ceil(double(n)/double(2.0)) )
339			if ( numPos < ceil_half_n) {
340			for (i=0; i<n; i++) {
341			v(i, j) *= -1.0;
342			}
343	tim	123	}
344	tim	93	}
345	tim	76
346	tim	146	if (n > 4) {
347			delete [] b;
348			delete [] z;
349	tim	123	}
350	tim	146	return 1;
351	tim	76	}
352
353
354			}
355	tim	123	#endif //MATH_SQUAREMATRIX_HPP
356	tim	76