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:

        /** 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. */
         SquareMatrix<Real, Dim>  inverse() {
             SquareMatrix<Real, Dim> result;

             return result;
        }        

        /** Returns the determinant of this matrix. */
        double determinant() const {
            double det;
            return det;
        }

        /** Returns the trace of this matrix. */
        double trace() const {
           double 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;
        }         

        void diagonalize() {
            jacobi(m, eigenValues, ortMat);
        }

        /**
         * Finds the eigenvalues and eigenvectors of a symmetric matrix
         * @param eigenvals a reference to a vector3 where the
         * eigenvalues will be stored. The eigenvalues are ordered so
         * that eigenvals[0] <= eigenvals[1] <= eigenvals[2].
         * @return an orthogonal matrix whose ith column is an
         * eigenvector for the eigenvalue eigenvals[i]
         */
        SquareMatrix<Real, Dim>  findEigenvectors(Vector<Real, Dim>& eigenValues) {
            SquareMatrix<Real, Dim> ortMat;
            
            if ( !isSymmetric()){
                throw();
            }
            
            SquareMatrix<Real, Dim> m(*this);
            jacobi(m, eigenValues, ortMat);

            return ortMat;
        }
        /**
         * Jacobi iteration routines for computing eigenvalues/eigenvectors of 
         * real symmetric matrix
         *
         * @return true if success, otherwise return false
         * @param a source matrix
         * @param w output eigenvalues 
         * @param v output eigenvectors 
         */
        bool jacobi(const SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w, 
                              SquareMatrix<Real, Dim>& v);
    };//end SquareMatrix


#define ROT(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 MAX_ROTATIONS 60

template<typename Real, int Dim>
bool SquareMatrix<Real, Dim>::jacobi(const SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w, 
                              SquareMatrix<Real, Dim>& v) {
    const int N = Dim;                                                                       
    int i, j, k, iq, ip;
    double tresh, theta, tau, t, sm, s, h, g, c;
    double tmp;
    Vector<Real, Dim> b, z;

    // 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<MAX_ROTATIONS; i++) {
        sm = 0.0;
        for (ip=0; ip<2; ip++) {
            for (iq=ip+1; iq<N; iq++)
                sm += fabs(a(ip, iq));
        }
        
        if (sm == 0.0)
            break;

        if (i < 4)
            tresh = 0.2*sm/(9);
        else
            tresh = 0.0;

        for (ip=0; ip<2; ip++) {
            for (iq=ip+1; iq<N; iq++) {
                g = 100.0*fabs(a(ip, iq));
                if (i > 4 && (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;
                    
                    for (j=0;j<ip-1;j++) 
                        ROT(a,j,ip,j,iq);

                    for (j=ip+1;j<iq-1;j++) 
                        ROT(a,ip,j,j,iq);

                    for (j=iq+1; j<N; j++) 
                        ROT(a,ip,j,iq,j);
                    for (j=0; j<N; j++) 
                        ROT(v,j,ip,j,iq);
                }
            }
        }//for (ip=0; ip<2; ip++) 

        for (ip=0; ip<N; ip++) {
            b(ip) += z(ip);
            w(ip) = b(ip);
            z(ip) = 0.0;
        }
        
    } // end for (i=0; i<MAX_ROTATIONS; i++) 

    if ( i >= MAX_ROTATIONS )
        return false;

    // sort eigenfunctions
    for (j=0; j<N; j++) {
        k = j;
        tmp = w(k);
        for (i=j; i<N; i++) {
            if (w(i) >= tmp) {
            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 numPos;
    for (j=0; j<N; j++) {
        for (numPos=0, i=0; i<N; i++) if ( v(i, j) >= 0.0 ) numPos++;
        if ( numPos < 2 ) for(i=0; i<N; i++) v(i, j) *= -1.0;
    }

    return true;
}

#undef ROT
#undef MAX_ROTATIONS

}

#endif //MATH_SQUAREMATRIX_HPP 
Revision:	93
Committed:	Sun Oct 17 01:19:11 2004 UTC (21 years ago) by tim
File size:	9746 byte(s)
Log Message:	math library in progress
#	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			#ifndef MATH_SQUAREMATRIX_HPP
33			#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
49			/** default constructor */
50			SquareMatrix() {
51			for (unsigned int i = 0; i < Dim; i++)
52			for (unsigned int j = 0; j < Dim; j++)
53			data_[i][j] = 0.0;
54			}
55
56			/** copy constructor */
57	tim	74	SquareMatrix(const RectMatrix<Real, Dim, Dim>& m) : RectMatrix<Real, Dim, Dim>(m) {
58	tim	70	}
59
60			/** copy assignment operator */
61	tim	74	SquareMatrix<Real, Dim>& operator =(const RectMatrix<Real, Dim, Dim>& m) {
62			RectMatrix<Real, Dim, Dim>::operator=(m);
63			return *this;
64	tim	70	}
65	tim	74
66			/** Retunrs an identity matrix*/
67	tim	70
68	tim	74	static SquareMatrix<Real, Dim> identity() {
69			SquareMatrix<Real, Dim> m;
70	tim	70
71			for (unsigned int i = 0; i < Dim; i++)
72	tim	74	for (unsigned int j = 0; j < Dim; j++)
73	tim	70	if (i == j)
74	tim	74	m(i, j) = 1.0;
75	tim	70	else
76	tim	74	m(i, j) = 0.0;
77
78			return m;
79	tim	70	}
80
81	tim	74	/** Retunrs the inversion of this matrix. */
82			SquareMatrix<Real, Dim> inverse() {
83			SquareMatrix<Real, Dim> result;
84
85			return result;
86	tim	76	}
87	tim	70
88			/** Returns the determinant of this matrix. */
89			double determinant() const {
90	tim	74	double det;
91			return det;
92	tim	70	}
93
94			/** Returns the trace of this matrix. */
95			double trace() const {
96			double tmp = 0;
97
98			for (unsigned int i = 0; i < Dim ; i++)
99			tmp += data_[i][i];
100
101			return tmp;
102			}
103
104			/** Tests if this matrix is symmetrix. */
105			bool isSymmetric() const {
106			for (unsigned int i = 0; i < Dim - 1; i++)
107			for (unsigned int j = i; j < Dim; j++)
108	tim	74	if (fabs(data_[i][j] - data_[j][i]) > oopse::epsilon)
109	tim	70	return false;
110
111			return true;
112			}
113
114	tim	76	/** Tests if this matrix is orthogonal. */
115	tim	74	bool isOrthogonal() {
116			SquareMatrix<Real, Dim> tmp;
117	tim	70
118	tim	74	tmp = this transpose();
119	tim	70
120	tim	76	return tmp.isDiagonal();
121	tim	70	}
122
123			/** Tests if this matrix is diagonal. */
124			bool isDiagonal() const {
125			for (unsigned int i = 0; i < Dim ; i++)
126			for (unsigned int j = 0; j < Dim; j++)
127	tim	74	if (i !=j && fabs(data_[i][j]) > oopse::epsilon)
128	tim	70	return false;
129
130			return true;
131			}
132
133			/** Tests if this matrix is the unit matrix. */
134			bool isUnitMatrix() const {
135			if (!isDiagonal())
136			return false;
137
138			for (unsigned int i = 0; i < Dim ; i++)
139	tim	74	if (fabs(data_[i][i] - 1) > oopse::epsilon)
140	tim	70	return false;
141
142			return true;
143	tim	74	}
144	tim	70
145	tim	76	void diagonalize() {
146			jacobi(m, eigenValues, ortMat);
147			}
148
149			/**
150			* Finds the eigenvalues and eigenvectors of a symmetric matrix
151			* @param eigenvals a reference to a vector3 where the
152			* eigenvalues will be stored. The eigenvalues are ordered so
153			* that eigenvals[0] <= eigenvals[1] <= eigenvals[2].
154			* @return an orthogonal matrix whose ith column is an
155			* eigenvector for the eigenvalue eigenvals[i]
156			*/
157			SquareMatrix<Real, Dim> findEigenvectors(Vector<Real, Dim>& eigenValues) {
158			SquareMatrix<Real, Dim> ortMat;
159
160			if ( !isSymmetric()){
161			throw();
162			}
163
164			SquareMatrix<Real, Dim> m(*this);
165			jacobi(m, eigenValues, ortMat);
166
167			return ortMat;
168			}
169			/**
170			* Jacobi iteration routines for computing eigenvalues/eigenvectors of
171			* real symmetric matrix
172			*
173			* @return true if success, otherwise return false
174			* @param a source matrix
175			* @param w output eigenvalues
176			* @param v output eigenvectors
177			*/
178	tim	83	bool jacobi(const SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
179	tim	76	SquareMatrix<Real, Dim>& v);
180	tim	70	};//end SquareMatrix
181
182	tim	76
183			#define ROT(a,i,j,k,l) g=a(i, j);h=a(k, l);a(i, j)=g-s(h+gtau);a(k, l)=h+s(g-htau)
184			#define MAX_ROTATIONS 60
185
186	tim	83	template<typename Real, int Dim>
187			bool SquareMatrix<Real, Dim>::jacobi(const SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
188			SquareMatrix<Real, Dim>& v) {
189	tim	76	const int N = Dim;
190			int i, j, k, iq, ip;
191			double tresh, theta, tau, t, sm, s, h, g, c;
192			double tmp;
193			Vector<Real, Dim> b, z;
194
195			// initialize
196	tim	93	for (ip=0; ip<N; ip++) {
197			for (iq=0; iq<N; iq++)
198			v(ip, iq) = 0.0;
199			v(ip, ip) = 1.0;
200	tim	76	}
201	tim	93
202			for (ip=0; ip<N; ip++) {
203			b(ip) = w(ip) = a(ip, ip);
204			z(ip) = 0.0;
205	tim	76	}
206
207			// begin rotation sequence
208	tim	93	for (i=0; i<MAX_ROTATIONS; i++) {
209			sm = 0.0;
210			for (ip=0; ip<2; ip++) {
211			for (iq=ip+1; iq<N; iq++)
212			sm += fabs(a(ip, iq));
213			}
214
215			if (sm == 0.0)
216			break;
217	tim	76
218	tim	93	if (i < 4)
219			tresh = 0.2*sm/(9);
220			else
221			tresh = 0.0;
222	tim	76
223	tim	93	for (ip=0; ip<2; ip++) {
224			for (iq=ip+1; iq<N; iq++) {
225			g = 100.0*fabs(a(ip, iq));
226			if (i > 4 && (fabs(w(ip))+g) == fabs(w(ip))
227			&& (fabs(w(iq))+g) == fabs(w(iq))) {
228			a(ip, iq) = 0.0;
229			} else if (fabs(a(ip, iq)) > tresh) {
230			h = w(iq) - w(ip);
231			if ( (fabs(h)+g) == fabs(h)) {
232			t = (a(ip, iq)) / h;
233			} else {
234			theta = 0.5*h / (a(ip, iq));
235			t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
236	tim	76
237	tim	93	if (theta < 0.0)
238			t = -t;
239			}
240	tim	76
241	tim	93	c = 1.0 / sqrt(1+t*t);
242			s = t*c;
243			tau = s/(1.0+c);
244			h = t*a(ip, iq);
245			z(ip) -= h;
246			z(iq) += h;
247			w(ip) -= h;
248			w(iq) += h;
249			a(ip, iq)=0.0;
250
251			for (j=0;j<ip-1;j++)
252			ROT(a,j,ip,j,iq);
253
254			for (j=ip+1;j<iq-1;j++)
255			ROT(a,ip,j,j,iq);
256
257			for (j=iq+1; j<N; j++)
258			ROT(a,ip,j,iq,j);
259			for (j=0; j<N; j++)
260			ROT(v,j,ip,j,iq);
261			}
262			}
263			}//for (ip=0; ip<2; ip++)
264
265			for (ip=0; ip<N; ip++) {
266			b(ip) += z(ip);
267			w(ip) = b(ip);
268			z(ip) = 0.0;
269			}
270
271			} // end for (i=0; i<MAX_ROTATIONS; i++)
272
273	tim	76	if ( i >= MAX_ROTATIONS )
274	tim	93	return false;
275	tim	76
276			// sort eigenfunctions
277	tim	93	for (j=0; j<N; j++) {
278			k = j;
279			tmp = w(k);
280			for (i=j; i<N; i++) {
281			if (w(i) >= tmp) {
282			k = i;
283			tmp = w(k);
284			}
285			}
286
287			if (k != j) {
288			w(k) = w(j);
289			w(j) = tmp;
290			for (i=0; i<N; i++) {
291			tmp = v(i, j);
292			v(i, j) = v(i, k);
293			v(i, k) = tmp;
294			}
295			}
296	tim	76	}
297
298			// insure eigenvector consistency (i.e., Jacobi can compute
299			// vectors that are negative of one another (.707,.707,0) and
300			// (-.707,-.707,0). This can reek havoc in
301			// hyperstreamline/other stuff. We will select the most
302			// positive eigenvector.
303			int numPos;
304	tim	93	for (j=0; j<N; j++) {
305			for (numPos=0, i=0; i<N; i++) if ( v(i, j) >= 0.0 ) numPos++;
306			if ( numPos < 2 ) for(i=0; i<N; i++) v(i, j) *= -1.0;
307	tim	76	}
308
309			return true;
310	tim	70	}
311	tim	76
312			#undef ROT
313			#undef MAX_ROTATIONS
314
315			}
316
317	tim	70	#endif //MATH_SQUAREMATRIX_HPP