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

             return result;
        }        

        /**
         * Returns the determinant of this matrix.
         * @todo need implementation
         */
        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;
        }         

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