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 
         */
        void 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<Real, int Dim>
void SquareMatrix<Real, int Dim>::jacobi(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<N; ip++) 
        {
            b(ip) += z(ip);
            w(ip) = b(ip);
            z(ip) = 0.0;
        }
    }

    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:	1569
Committed:	Thu Oct 14 23:28:09 2004 UTC (19 years, 8 months ago) by tim
File size:	8995 byte(s)
Log Message:	math library in progress
#	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			#ifndef MATH_SQUAREMATRIX_HPP
33			#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
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	1567	SquareMatrix(const RectMatrix<Real, Dim, Dim>& m) : RectMatrix<Real, Dim, Dim>(m) {
58	tim	1563	}
59
60			/** copy assignment operator */
61	tim	1567	SquareMatrix<Real, Dim>& operator =(const RectMatrix<Real, Dim, Dim>& m) {
62			RectMatrix<Real, Dim, Dim>::operator=(m);
63			return *this;
64	tim	1563	}
65	tim	1567
66			/** Retunrs an identity matrix*/
67	tim	1563
68	tim	1567	static SquareMatrix<Real, Dim> identity() {
69			SquareMatrix<Real, Dim> m;
70	tim	1563
71			for (unsigned int i = 0; i < Dim; i++)
72	tim	1567	for (unsigned int j = 0; j < Dim; j++)
73	tim	1563	if (i == j)
74	tim	1567	m(i, j) = 1.0;
75	tim	1563	else
76	tim	1567	m(i, j) = 0.0;
77
78			return m;
79	tim	1563	}
80
81	tim	1567	/** Retunrs the inversion of this matrix. */
82			SquareMatrix<Real, Dim> inverse() {
83			SquareMatrix<Real, Dim> result;
84
85			return result;
86	tim	1569	}
87	tim	1563
88			/** Returns the determinant of this matrix. */
89			double determinant() const {
90	tim	1567	double det;
91			return det;
92	tim	1563	}
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	1567	if (fabs(data_[i][j] - data_[j][i]) > oopse::epsilon)
109	tim	1563	return false;
110
111			return true;
112			}
113
114	tim	1569	/** Tests if this matrix is orthogonal. */
115	tim	1567	bool isOrthogonal() {
116			SquareMatrix<Real, Dim> tmp;
117	tim	1563
118	tim	1567	tmp = this transpose();
119	tim	1563
120	tim	1569	return tmp.isDiagonal();
121	tim	1563	}
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	1567	if (i !=j && fabs(data_[i][j]) > oopse::epsilon)
128	tim	1563	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	1567	if (fabs(data_[i][i] - 1) > oopse::epsilon)
140	tim	1563	return false;
141
142			return true;
143	tim	1567	}
144	tim	1563
145	tim	1569	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			void jacobi(const SquareMatrix<Real, Dim>& a,
179			Vector<Real, Dim>& w,
180			SquareMatrix<Real, Dim>& v);
181	tim	1563	};//end SquareMatrix
182
183	tim	1569
184			#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)
185			#define MAX_ROTATIONS 60
186
187			template<Real, int Dim>
188			void SquareMatrix<Real, int Dim>::jacobi(SquareMatrix<Real, Dim>& a,
189			Vector<Real, Dim>& w,
190			SquareMatrix<Real, Dim>& v) {
191			const int N = Dim;
192			int i, j, k, iq, ip;
193			double tresh, theta, tau, t, sm, s, h, g, c;
194			double tmp;
195			Vector<Real, Dim> b, z;
196
197			// initialize
198			for (ip=0; ip<N; ip++)
199			{
200			for (iq=0; iq<N; iq++) v(ip, iq) = 0.0;
201			v(ip, ip) = 1.0;
202			}
203			for (ip=0; ip<N; ip++)
204			{
205			b(ip) = w(ip) = a(ip, ip);
206			z(ip) = 0.0;
207			}
208
209			// begin rotation sequence
210			for (i=0; i<MAX_ROTATIONS; i++)
211			{
212			sm = 0.0;
213			for (ip=0; ip<2; ip++)
214			{
215			for (iq=ip+1; iq<N; iq++) sm += fabs(a(ip, iq));
216			}
217			if (sm == 0.0) break;
218
219			if (i < 4) tresh = 0.2*sm/(9);
220			else tresh = 0.0;
221
222			for (ip=0; ip<2; ip++)
223			{
224			for (iq=ip+1; iq<N; iq++)
225			{
226			g = 100.0*fabs(a(ip, iq));
227			if (i > 4 && (fabs(w(ip))+g) == fabs(w(ip))
228			&& (fabs(w(iq))+g) == fabs(w(iq)))
229			{
230			a(ip, iq) = 0.0;
231			}
232			else if (fabs(a(ip, iq)) > tresh)
233			{
234			h = w(iq) - w(ip);
235			if ( (fabs(h)+g) == fabs(h)) t = (a(ip, iq)) / h;
236			else
237			{
238			theta = 0.5*h / (a(ip, iq));
239			t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
240			if (theta < 0.0) t = -t;
241			}
242			c = 1.0 / sqrt(1+t*t);
243			s = t*c;
244			tau = s/(1.0+c);
245			h = t*a(ip, iq);
246			z(ip) -= h;
247			z(iq) += h;
248			w(ip) -= h;
249			w(iq) += h;
250			a(ip, iq)=0.0;
251			for (j=0;j<ip-1;j++)
252			{
253			ROT(a,j,ip,j,iq);
254			}
255			for (j=ip+1;j<iq-1;j++)
256			{
257			ROT(a,ip,j,j,iq);
258			}
259			for (j=iq+1; j<N; j++)
260			{
261			ROT(a,ip,j,iq,j);
262			}
263			for (j=0; j<N; j++)
264			{
265			ROT(v,j,ip,j,iq);
266			}
267			}
268			}
269			}
270
271			for (ip=0; ip<N; ip++)
272			{
273			b(ip) += z(ip);
274			w(ip) = b(ip);
275			z(ip) = 0.0;
276			}
277			}
278
279			if ( i >= MAX_ROTATIONS )
280			return false;
281
282			// sort eigenfunctions
283			for (j=0; j<N; j++)
284			{
285			k = j;
286			tmp = w(k);
287			for (i=j; i<N; i++)
288			{
289			if (w(i) >= tmp)
290			{
291			k = i;
292			tmp = w(k);
293			}
294			}
295			if (k != j)
296			{
297			w(k) = w(j);
298			w(j) = tmp;
299			for (i=0; i<N; i++)
300			{
301			tmp = v(i, j);
302			v(i, j) = v(i, k);
303			v(i, k) = tmp;
304			}
305			}
306			}
307
308			// insure eigenvector consistency (i.e., Jacobi can compute
309			// vectors that are negative of one another (.707,.707,0) and
310			// (-.707,-.707,0). This can reek havoc in
311			// hyperstreamline/other stuff. We will select the most
312			// positive eigenvector.
313			int numPos;
314			for (j=0; j<N; j++)
315			{
316			for (numPos=0, i=0; i<N; i++) if ( v(i, j) >= 0.0 ) numPos++;
317			if ( numPos < 2 ) for(i=0; i<N; i++) v(i, j) *= -1.0;
318			}
319
320			return true;
321	tim	1563	}
322	tim	1569
323			#undef ROT
324			#undef MAX_ROTATIONS
325
326			}
327
328
329			}
330	tim	1563	#endif //MATH_SQUAREMATRIX_HPP