ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/trunk/src/math/SquareMatrix3.hpp

Comparing trunk/src/math/SquareMatrix3.hpp (file contents):
Revision 93 by tim, Sun Oct 17 01:19:11 2004 UTC vs.
Revision 2000 by gezelter, Sat May 31 22:35:05 2014 UTC

#	Line 1 | Line 1
1		/*
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.
2	>	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3		*
4	+	* The University of Notre Dame grants you ("Licensee") a
5	+	* non-exclusive, royalty free, license to use, modify and
6	+	* redistribute this software in source and binary code form, provided
7	+	* that the following conditions are met:
8	+	*
9	+	* 1. Redistributions of source code must retain the above copyright
10	+	*    notice, this list of conditions and the following disclaimer.
11	+	*
12	+	* 2. Redistributions in binary form must reproduce the above copyright
13	+	*    notice, this list of conditions and the following disclaimer in the
14	+	*    documentation and/or other materials provided with the
15	+	*    distribution.
16	+	*
17	+	* This software is provided "AS IS," without a warranty of any
18	+	* kind. All express or implied conditions, representations and
19	+	* warranties, including any implied warranty of merchantability,
20	+	* fitness for a particular purpose or non-infringement, are hereby
21	+	* excluded.  The University of Notre Dame and its licensors shall not
22	+	* be liable for any damages suffered by licensee as a result of
23	+	* using, modifying or distributing the software or its
24	+	* derivatives. In no event will the University of Notre Dame or its
25	+	* licensors be liable for any lost revenue, profit or data, or for
26	+	* direct, indirect, special, consequential, incidental or punitive
27	+	* damages, however caused and regardless of the theory of liability,
28	+	* arising out of the use of or inability to use software, even if the
29	+	* University of Notre Dame has been advised of the possibility of
30	+	* such damages.
31	+	*
32	+	* SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33	+	* research, please cite the appropriate papers when you publish your
34	+	* work.  Good starting points are:
35	+	*
36	+	* [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37	+	* [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	+	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39	+	* [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40	+	* [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41		*/
42	<
42	>
43		/**
44		* @file SquareMatrix3.hpp
45		* @author Teng Lin
46		* @date 10/11/2004
47		* @version 1.0
48		*/
49	<	#ifndef MATH_SQUAREMATRIX_HPP
50	<	#define  MATH_SQUAREMATRIX_HPP
51	<
49	>	#ifndef MATH_SQUAREMATRIX3_HPP
50	>	#define  MATH_SQUAREMATRIX3_HPP
51	>	#include "config.h"
52	>	#include <cmath>
53	>	#include <vector>
54		#include "Quaternion.hpp"
55		#include "SquareMatrix.hpp"
56		#include "Vector3.hpp"
57	+	#include "utils/NumericConstant.hpp"
58	+	namespace OpenMD {
59
60	<	namespace oopse {
60	>	template<typename Real>
61	>	class SquareMatrix3 : public SquareMatrix<Real, 3> {
62	>	public:
63
64	<	template<typename Real>
65	<	class SquareMatrix3 : public SquareMatrix<Real, 3> {
43	<	public:
64	>	typedef Real ElemType;
65	>	typedef Real* ElemPoinerType;
66
67	<	/** default constructor */
68	<	SquareMatrix3() : SquareMatrix<Real, 3>() {
69	<	}
67	>	/** default constructor */
68	>	SquareMatrix3() : SquareMatrix<Real, 3>() {
69	>	}
70
71	<	/** copy  constructor */
72	<	SquareMatrix3(const SquareMatrix<Real, 3>& m)  : SquareMatrix<Real, 3>(m) {
73	<	}
71	>	/** Constructs and initializes every element of this matrix to a scalar */
72	>	SquareMatrix3(Real s) : SquareMatrix<Real,3>(s){
73	>	}
74
75	<	SquareMatrix3( const Vector3<Real>& eulerAngles) {
76	<	setupRotMat(eulerAngles);
77	<	}
56	<
57	<	SquareMatrix3(Real phi, Real theta, Real psi) {
58	<	setupRotMat(phi, theta, psi);
59	<	}
75	>	/** Constructs and initializes from an array */
76	>	SquareMatrix3(Real* array) : SquareMatrix<Real,3>(array){
77	>	}
78
61	–	SquareMatrix3(const Quaternion<Real>& q) {
62	–	*this = q.toRotationMatrix3();
63	–	}
79
80	<	SquareMatrix3(Real w, Real x, Real y, Real z) {
81	<	Quaternion<Real> q(w, x, y, z);
82	<	*this = q.toRotationMatrix3();
68	<	}
80	>	/** copy  constructor */
81	>	SquareMatrix3(const SquareMatrix<Real, 3>& m)  : SquareMatrix<Real, 3>(m) {
82	>	}
83
84	<	/** copy assignment operator */
85	<	SquareMatrix3<Real>& operator =(const SquareMatrix<Real, 3>& m) {
86	<	if (this == &m)
87	<	return *this;
88	<	SquareMatrix<Real, 3>::operator=(m);
89	<	}
84	>	SquareMatrix3( const Vector3<Real>& eulerAngles) {
85	>	setupRotMat(eulerAngles);
86	>	}
87	>
88	>	SquareMatrix3(Real phi, Real theta, Real psi) {
89	>	setupRotMat(phi, theta, psi);
90	>	}
91
92	<	/**
93	<	* Sets this matrix to a rotation matrix by three euler angles
79	<	* @ param euler
80	<	*/
81	<	void setupRotMat(const Vector3<Real>& eulerAngles) {
82	<	setupRotMat(eulerAngles[0], eulerAngles[1], eulerAngles[2]);
83	<	}
92	>	SquareMatrix3(const Quaternion<Real>& q) {
93	>	setupRotMat(q);
94
95	<	/**
86	<	* Sets this matrix to a rotation matrix by three euler angles
87	<	* @param phi
88	<	* @param theta
89	<	* @psi theta
90	<	*/
91	<	void setupRotMat(Real phi, Real theta, Real psi) {
92	<	Real sphi, stheta, spsi;
93	<	Real cphi, ctheta, cpsi;
95	>	}
96
97	<	sphi = sin(phi);
98	<	stheta = sin(theta);
99	<	spsi = sin(psi);
100	<	cphi = cos(phi);
101	<	ctheta = cos(theta);
102	<	cpsi = cos(psi);
97	>	SquareMatrix3(Real w, Real x, Real y, Real z) {
98	>	setupRotMat(w, x, y, z);
99	>	}
100	>
101	>	/** copy assignment operator */
102	>	SquareMatrix3<Real>& operator =(const SquareMatrix<Real, 3>& m) {
103	>	if (this == &m)
104	>	return *this;
105	>	SquareMatrix<Real, 3>::operator=(m);
106	>	return *this;
107	>	}
108
102	–	data_[0][0] = cpsi * cphi - ctheta * sphi * spsi;
103	–	data_[0][1] = cpsi * sphi + ctheta * cphi * spsi;
104	–	data_[0][2] = spsi * stheta;
105	–
106	–	data_[1][0] = -spsi * ctheta - ctheta * sphi * cpsi;
107	–	data_[1][1] = -spsi * stheta + ctheta * cphi * cpsi;
108	–	data_[1][2] = cpsi * stheta;
109
110	<	data_[2][0] = stheta * sphi;
111	<	data_[2][1] = -stheta * cphi;
112	<	data_[2][2] = ctheta;
113	<	}
110	>	SquareMatrix3<Real>& operator =(const Quaternion<Real>& q) {
111	>	this->setupRotMat(q);
112	>	return *this;
113	>	}
114
115	+
116	+	/**
117	+	* Sets this matrix to a rotation matrix by three euler angles
118	+	* @ param euler
119	+	*/
120	+	void setupRotMat(const Vector3<Real>& eulerAngles) {
121	+	setupRotMat(eulerAngles[0], eulerAngles[1], eulerAngles[2]);
122	+	}
123
124	<	/**
125	<	* Sets this matrix to a rotation matrix by quaternion
126	<	* @param quat
127	<	*/
128	<	void setupRotMat(const Quaternion<Real>& quat) {
129	<	*this = quat.toRotationMatrix3();
130	<	}
124	>	/**
125	>	* Sets this matrix to a rotation matrix by three euler angles
126	>	* @param phi
127	>	* @param theta
128	>	* @param psi
129	>	*/
130	>	void setupRotMat(Real phi, Real theta, Real psi) {
131	>	Real sphi, stheta, spsi;
132	>	Real cphi, ctheta, cpsi;
133
134	<	/**
135	<	* Sets this matrix to a rotation matrix by quaternion
136	<	* @param w the first element
137	<	* @param x the second element
138	<	* @param y the third element
139	<	* @parma z the fourth element
130	<	*/
131	<	void setupRotMat(Real w, Real x, Real y, Real z) {
132	<	Quaternion<Real> q(w, x, y, z);
133	<	*this = q.toRotationMatrix3();
134	<	}
134	>	sphi = sin(phi);
135	>	stheta = sin(theta);
136	>	spsi = sin(psi);
137	>	cphi = cos(phi);
138	>	ctheta = cos(theta);
139	>	cpsi = cos(psi);
140
141	<	/**
142	<	* Returns the quaternion from this rotation matrix
143	<	* @return the quaternion from this rotation matrix
144	<	* @exception invalid rotation matrix
145	<	*/
146	<	Quaternion<Real> toQuaternion() {
147	<	Quaternion<Real> q;
143	<	Real t, s;
144	<	Real ad1, ad2, ad3;
145	<	t = data_[0][0] + data_[1][1] + data_[2][2] + 1.0;
141	>	this->data_[0][0] = cpsi * cphi - ctheta * sphi * spsi;
142	>	this->data_[0][1] = cpsi * sphi + ctheta * cphi * spsi;
143	>	this->data_[0][2] = spsi * stheta;
144	>
145	>	this->data_[1][0] = -spsi * ctheta - ctheta * sphi * cpsi;
146	>	this->data_[1][1] = -spsi * stheta + ctheta * cphi * cpsi;
147	>	this->data_[1][2] = cpsi * stheta;
148
149	<	if( t > 0.0 ){
149	>	this->data_[2][0] = stheta * sphi;
150	>	this->data_[2][1] = -stheta * cphi;
151	>	this->data_[2][2] = ctheta;
152	>	}
153
149	–	s = 0.5 / sqrt( t );
150	–	q[0] = 0.25 / s;
151	–	q[1] = (data_[1][2] - data_[2][1]) * s;
152	–	q[2] = (data_[2][0] - data_[0][2]) * s;
153	–	q[3] = (data_[0][1] - data_[1][0]) * s;
154	–	} else {
154
155	<	ad1 = fabs( data_[0][0] );
156	<	ad2 = fabs( data_[1][1] );
157	<	ad3 = fabs( data_[2][2] );
155	>	/**
156	>	* Sets this matrix to a rotation matrix by quaternion
157	>	* @param quat
158	>	*/
159	>	void setupRotMat(const Quaternion<Real>& quat) {
160	>	setupRotMat(quat.w(), quat.x(), quat.y(), quat.z());
161	>	}
162
163	<	if( ad1 >= ad2 && ad1 >= ad3 ){
163	>	/**
164	>	* Sets this matrix to a rotation matrix by quaternion
165	>	* @param w the first element
166	>	* @param x the second element
167	>	* @param y the third element
168	>	* @param z the fourth element
169	>	*/
170	>	void setupRotMat(Real w, Real x, Real y, Real z) {
171	>	Quaternion<Real> q(w, x, y, z);
172	>	*this = q.toRotationMatrix3();
173	>	}
174
175	<	s = 2.0 * sqrt( 1.0 + data_[0][0] - data_[1][1] - data_[2][2] );
176	<	q[0] = (data_[1][2] + data_[2][1]) / s;
177	<	q[1] = 0.5 / s;
165	<	q[2] = (data_[0][1] + data_[1][0]) / s;
166	<	q[3] = (data_[0][2] + data_[2][0]) / s;
167	<	} else if ( ad2 >= ad1 && ad2 >= ad3 ) {
168	<	s = sqrt( 1.0 + data_[1][1] - data_[0][0] - data_[2][2] ) * 2.0;
169	<	q[0] = (data_[0][2] + data_[2][0]) / s;
170	<	q[1] = (data_[0][1] + data_[1][0]) / s;
171	<	q[2] = 0.5 / s;
172	<	q[3] = (data_[1][2] + data_[2][1]) / s;
173	<	} else {
175	>	void setupSkewMat(Vector3<Real> v) {
176	>	setupSkewMat(v[0], v[1], v[2]);
177	>	}
178
179	<	s = sqrt( 1.0 + data_[2][2] - data_[0][0] - data_[1][1] ) * 2.0;
180	<	q[0] = (data_[0][1] + data_[1][0]) / s;
181	<	q[1] = (data_[0][2] + data_[2][0]) / s;
182	<	q[2] = (data_[1][2] + data_[2][1]) / s;
183	<	q[3] = 0.5 / s;
184	<	}
185	<	}
179	>	void setupSkewMat(Real v1, Real v2, Real v3) {
180	>	this->data_[0][0] = 0;
181	>	this->data_[0][1] = -v3;
182	>	this->data_[0][2] = v2;
183	>	this->data_[1][0] = v3;
184	>	this->data_[1][1] = 0;
185	>	this->data_[1][2] = -v1;
186	>	this->data_[2][0] = -v2;
187	>	this->data_[2][1] = v1;
188	>	this->data_[2][2] = 0;
189	>
190	>
191	>	}
192
183	–	return q;
184	–
185	–	}
193
194	<	/**
195	<	* Returns the euler angles from this rotation matrix
196	<	* @return the euler angles in a vector
197	<	* @exception invalid rotation matrix
198	<	* We use so-called "x-convention", which is the most common definition.
199	<	* In this convention, the rotation given by Euler angles (phi, theta, psi), where the first
200	<	* rotation is by an angle phi about the z-axis, the second is by an angle
201	<	* theta (0 <= theta <= 180)about the x-axis, and thethird is by an angle psi about the
202	<	* z-axis (again).
203	<	*/
197	<	Vector3<Real> toEulerAngles() {
198	<	Vector<Real> myEuler;
199	<	Real phi,theta,psi,eps;
200	<	Real ctheta,stheta;
201	<
202	<	// set the tolerance for Euler angles and rotation elements
194	>	/**
195	>	* Returns the quaternion from this rotation matrix
196	>	* @return the quaternion from this rotation matrix
197	>	* @exception invalid rotation matrix
198	>	*/
199	>	Quaternion<Real> toQuaternion() {
200	>	Quaternion<Real> q;
201	>	Real t, s;
202	>	Real ad1, ad2, ad3;
203	>	t = this->data_[0][0] + this->data_[1][1] + this->data_[2][2] + 1.0;
204
205	<	theta = acos(min(1.0,max(-1.0,data_[2][2])));
205	<	ctheta = data_[2][2];
206	<	stheta = sqrt(1.0 - ctheta * ctheta);
205	>	if( t > NumericConstant::epsilon ){
206
207	<	// when sin(theta) is close to 0, we need to consider singularity
208	<	// In this case, we can assign an arbitary value to phi (or psi), and then determine
209	<	// the psi (or phi) or vice-versa. We'll assume that phi always gets the rotation, and psi is 0
210	<	// in cases of singularity.
211	<	// we use atan2 instead of atan, since atan2 will give us -Pi to Pi.
212	<	// Since 0 <= theta <= 180, sin(theta) will be always non-negative. Therefore, it never
214	<	// change the sign of both of the parameters passed to atan2.
207	>	s = 0.5 / sqrt( t );
208	>	q[0] = 0.25 / s;
209	>	q[1] = (this->data_[1][2] - this->data_[2][1]) * s;
210	>	q[2] = (this->data_[2][0] - this->data_[0][2]) * s;
211	>	q[3] = (this->data_[0][1] - this->data_[1][0]) * s;
212	>	} else {
213
214	<	if (fabs(stheta) <= oopse::epsilon){
215	<	psi = 0.0;
216	<	phi = atan2(-data_[1][0], data_[0][0]);
219	<	}
220	<	// we only have one unique solution
221	<	else{
222	<	phi = atan2(data_[2][0], -data_[2][1]);
223	<	psi = atan2(data_[0][2], data_[1][2]);
224	<	}
214	>	ad1 = this->data_[0][0];
215	>	ad2 = this->data_[1][1];
216	>	ad3 = this->data_[2][2];
217
218	<	//wrap phi and psi, make sure they are in the range from 0 to 2*Pi
227	<	if (phi < 0)
228	<	phi += M_PI;
218	>	if( ad1 >= ad2 && ad1 >= ad3 ){
219
220	<	if (psi < 0)
221	<	psi += M_PI;
222	<
223	<	myEuler[0] = phi;
224	<	myEuler[1] = theta;
225	<	myEuler[2] = psi;
220	>	s = 0.5 / sqrt( 1.0 + this->data_[0][0] - this->data_[1][1] - this->data_[2][2] );
221	>	q[0] = (this->data_[1][2] - this->data_[2][1]) * s;
222	>	q[1] = 0.25 / s;
223	>	q[2] = (this->data_[0][1] + this->data_[1][0]) * s;
224	>	q[3] = (this->data_[0][2] + this->data_[2][0]) * s;
225	>	} else if ( ad2 >= ad1 && ad2 >= ad3 ) {
226	>	s = 0.5 / sqrt( 1.0 + this->data_[1][1] - this->data_[0][0] - this->data_[2][2] );
227	>	q[0] = (this->data_[2][0] - this->data_[0][2] ) * s;
228	>	q[1] = (this->data_[0][1] + this->data_[1][0]) * s;
229	>	q[2] = 0.25 / s;
230	>	q[3] = (this->data_[1][2] + this->data_[2][1]) * s;
231	>	} else {
232
233	<	return myEuler;
233	>	s = 0.5 / sqrt( 1.0 + this->data_[2][2] - this->data_[0][0] - this->data_[1][1] );
234	>	q[0] = (this->data_[0][1] - this->data_[1][0]) * s;
235	>	q[1] = (this->data_[0][2] + this->data_[2][0]) * s;
236	>	q[2] = (this->data_[1][2] + this->data_[2][1]) * s;
237	>	q[3] = 0.25 / s;
238	>	}
239	>	}
240	>
241	>	return q;
242	>
243	>	}
244	>
245	>	/**
246	>	* Returns the euler angles from this rotation matrix
247	>	* @return the euler angles in a vector
248	>	* @exception invalid rotation matrix
249	>	* We use so-called "x-convention", which is the most common definition.
250	>	* In this convention, the rotation given by Euler angles (phi, theta,
251	>	* psi), where the first rotation is by an angle phi about the z-axis,
252	>	* the second is by an angle theta (0 <= theta <= 180) about the x-axis,
253	>	* and the third is by an angle psi about the z-axis (again).
254	>	*/
255	>	Vector3<Real> toEulerAngles() {
256	>	Vector3<Real> myEuler;
257	>	Real phi;
258	>	Real theta;
259	>	Real psi;
260	>	Real ctheta;
261	>	Real stheta;
262	>
263	>	// set the tolerance for Euler angles and rotation elements
264	>
265	>	theta = acos(std::min((RealType)1.0, std::max((RealType)-1.0,this->data_[2][2])));
266	>	ctheta = this->data_[2][2];
267	>	stheta = sqrt(1.0 - ctheta * ctheta);
268	>
269	>	// when sin(theta) is close to 0, we need to consider
270	>	// singularity In this case, we can assign an arbitary value to
271	>	// phi (or psi), and then determine the psi (or phi) or
272	>	// vice-versa. We'll assume that phi always gets the rotation,
273	>	// and psi is 0 in cases of singularity.
274	>	// we use atan2 instead of atan, since atan2 will give us -Pi to Pi.
275	>	// Since 0 <= theta <= 180, sin(theta) will be always
276	>	// non-negative. Therefore, it will never change the sign of both of
277	>	// the parameters passed to atan2.
278	>
279	>	if (fabs(stheta) < 1e-6){
280	>	psi = 0.0;
281	>	phi = atan2(-this->data_[1][0], this->data_[0][0]);
282	>	}
283	>	// we only have one unique solution
284	>	else{
285	>	phi = atan2(this->data_[2][0], -this->data_[2][1]);
286	>	psi = atan2(this->data_[0][2], this->data_[1][2]);
287	>	}
288	>
289	>	//wrap phi and psi, make sure they are in the range from 0 to 2*Pi
290	>	if (phi < 0)
291	>	phi += 2.0 * M_PI;
292	>
293	>	if (psi < 0)
294	>	psi += 2.0 * M_PI;
295	>
296	>	myEuler[0] = phi;
297	>	myEuler[1] = theta;
298	>	myEuler[2] = psi;
299	>
300	>	return myEuler;
301	>	}
302	>
303	>	/** Returns the determinant of this matrix. */
304	>	Real determinant() const {
305	>	Real x,y,z;
306	>
307	>	x = this->data_[0][0] * (this->data_[1][1] * this->data_[2][2] - this->data_[1][2] * this->data_[2][1]);
308	>	y = this->data_[0][1] * (this->data_[1][2] * this->data_[2][0] - this->data_[1][0] * this->data_[2][2]);
309	>	z = this->data_[0][2] * (this->data_[1][0] * this->data_[2][1] - this->data_[1][1] * this->data_[2][0]);
310	>
311	>	return(x + y + z);
312	>	}
313	>
314	>	/** Returns the trace of this matrix. */
315	>	Real trace() const {
316	>	return this->data_[0][0] + this->data_[1][1] + this->data_[2][2];
317	>	}
318	>
319	>	/**
320	>	* Sets the value of this matrix to  the inversion of itself.
321	>	* @note since simple algorithm can be applied to inverse the 3 by 3 matrix, we hide the
322	>	* implementation of inverse in SquareMatrix class
323	>	*/
324	>	SquareMatrix3<Real>  inverse() const {
325	>	SquareMatrix3<Real> m;
326	>	RealType det = determinant();
327	>	if (fabs(det) <= OpenMD::epsilon) {
328	>	//"The method was called on a matrix with |determinant| <= 1e-6.",
329	>	//"This is a runtime or a programming error in your application.");
330	>	std::vector<int> zeroDiagElementIndex;
331	>	for (int i =0; i < 3; ++i) {
332	>	if (fabs(this->data_[i][i]) <= OpenMD::epsilon) {
333	>	zeroDiagElementIndex.push_back(i);
334		}
335	+	}
336	+
337	+	if (zeroDiagElementIndex.size() == 2) {
338	+	int index = zeroDiagElementIndex[0];
339	+	m(index, index) = 1.0 / this->data_[index][index];
340	+	}else if (zeroDiagElementIndex.size() == 1) {
341	+
342	+	int a = (zeroDiagElementIndex[0] + 1) % 3;
343	+	int b = (zeroDiagElementIndex[0] + 2) %3;
344	+	RealType denom = this->data_[a][a] * this->data_[b][b] - this->data_[b][a]*this->data_[a][b];
345	+	m(a, a) = this->data_[b][b] /denom;
346	+	m(b, a) = -this->data_[b][a]/denom;
347	+
348	+	m(a,b) = -this->data_[a][b]/denom;
349	+	m(b, b) = this->data_[a][a]/denom;
350	+
351	+	}
352	+
353	+	/*
354	+	for(std::vector<int>::iterator iter = zeroDiagElementIndex.begin(); iter != zeroDiagElementIndex.end() ++iter) {
355	+	if (this->data_[*iter][0] > OpenMD::epsilon || this->data_[*iter][1] ||this->data_[*iter][2] ||
356	+	this->data_[0][*iter] > OpenMD::epsilon || this->data_[1][*iter] ||this->data_[2][*iter] ) {
357	+	std::cout << "can not inverse matrix" << std::endl;
358	+	}
359	+	}
360	+	*/
361	+	} else {
362	+
363	+	m(0, 0) = this->data_[1][1]*this->data_[2][2] - this->data_[1][2]*this->data_[2][1];
364	+	m(1, 0) = this->data_[1][2]*this->data_[2][0] - this->data_[1][0]*this->data_[2][2];
365	+	m(2, 0) = this->data_[1][0]*this->data_[2][1] - this->data_[1][1]*this->data_[2][0];
366	+	m(0, 1) = this->data_[2][1]*this->data_[0][2] - this->data_[2][2]*this->data_[0][1];
367	+	m(1, 1) = this->data_[2][2]*this->data_[0][0] - this->data_[2][0]*this->data_[0][2];
368	+	m(2, 1) = this->data_[2][0]*this->data_[0][1] - this->data_[2][1]*this->data_[0][0];
369	+	m(0, 2) = this->data_[0][1]*this->data_[1][2] - this->data_[0][2]*this->data_[1][1];
370	+	m(1, 2) = this->data_[0][2]*this->data_[1][0] - this->data_[0][0]*this->data_[1][2];
371	+	m(2, 2) = this->data_[0][0]*this->data_[1][1] - this->data_[0][1]*this->data_[1][0];
372	+
373	+	m /= det;
374	+	}
375	+	return m;
376	+	}
377	+
378	+	SquareMatrix3<Real> transpose() const{
379	+	SquareMatrix3<Real> result;
380	+
381	+	for (unsigned int i = 0; i < 3; i++)
382	+	for (unsigned int j = 0; j < 3; j++)
383	+	result(j, i) = this->data_[i][j];
384	+
385	+	return result;
386	+	}
387	+	/**
388	+	* Extract the eigenvalues and eigenvectors from a 3x3 matrix.
389	+	* The eigenvectors (the columns of V) will be normalized.
390	+	* The eigenvectors are aligned optimally with the x, y, and z
391	+	* axes respectively.
392	+	* @param a symmetric matrix whose eigenvectors are to be computed. On return, the matrix is
393	+	*     overwritten
394	+	* @param w will contain the eigenvalues of the matrix On return of this function
395	+	* @param v the columns of this matrix will contain the eigenvectors. The eigenvectors are
396	+	*    normalized and mutually orthogonal.
397	+	* @warning a will be overwritten
398	+	*/
399	+	static void diagonalize(SquareMatrix3<Real>& a, Vector3<Real>& w, SquareMatrix3<Real>& v);
400	+	};
401	+	/*=========================================================================
402	+
403	+	Program:   Visualization Toolkit
404	+	Module:    $RCSfile: SquareMatrix3.hpp,v $
405	+
406	+	Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
407	+	All rights reserved.
408	+	See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
409	+
410	+	This software is distributed WITHOUT ANY WARRANTY; without even
411	+	the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
412	+	PURPOSE.  See the above copyright notice for more information.
413	+
414	+	=========================================================================*/
415	+	template<typename Real>
416	+	void SquareMatrix3<Real>::diagonalize(SquareMatrix3<Real>& a, Vector3<Real>& w,
417	+	SquareMatrix3<Real>& v) {
418	+	int i,j,k,maxI;
419	+	Real tmp, maxVal;
420	+	Vector3<Real> v_maxI, v_k, v_j;
421	+
422	+	// diagonalize using Jacobi
423	+	SquareMatrix3<Real>::jacobi(a, w, v);
424	+	// if all the eigenvalues are the same, return identity matrix
425	+	if (w[0] == w[1] && w[0] == w[2] ) {
426	+	v = SquareMatrix3<Real>::identity();
427	+	return;
428	+	}
429	+
430	+	// transpose temporarily, it makes it easier to sort the eigenvectors
431	+	v = v.transpose();
432	+
433	+	// if two eigenvalues are the same, re-orthogonalize to optimally line
434	+	// up the eigenvectors with the x, y, and z axes
435	+	for (i = 0; i < 3; i++) {
436	+	if (w((i+1)%3) == w((i+2)%3)) {// two eigenvalues are the same
437	+	// find maximum element of the independant eigenvector
438	+	maxVal = fabs(v(i, 0));
439	+	maxI = 0;
440	+	for (j = 1; j < 3; j++) {
441	+	if (maxVal < (tmp = fabs(v(i, j)))){
442	+	maxVal = tmp;
443	+	maxI = j;
444	+	}
445	+	}
446
447	<	/**
448	<	* Sets the value of this matrix to  the inversion of itself.
449	<	* @note since simple algorithm can be applied to inverse the 3 by 3 matrix, we hide the
450	<	* implementation of inverse in SquareMatrix class
451	<	*/
245	<	void  inverse();
447	>	// swap the eigenvector into its proper position
448	>	if (maxI != i) {
449	>	tmp = w(maxI);
450	>	w(maxI) = w(i);
451	>	w(i) = tmp;
452
453	<	void diagonalize();
453	>	v.swapRow(i, maxI);
454	>	}
455	>	// maximum element of eigenvector should be positive
456	>	if (v(maxI, maxI) < 0) {
457	>	v(maxI, 0) = -v(maxI, 0);
458	>	v(maxI, 1) = -v(maxI, 1);
459	>	v(maxI, 2) = -v(maxI, 2);
460	>	}
461
462	<	};
462	>	// re-orthogonalize the other two eigenvectors
463	>	j = (maxI+1)%3;
464	>	k = (maxI+2)%3;
465
466	<	typedef template SquareMatrix3<double> Mat3x3d
467	<	typedef template SquareMatrix3<double> RotMat3x3d;
466	>	v(j, 0) = 0.0;
467	>	v(j, 1) = 0.0;
468	>	v(j, 2) = 0.0;
469	>	v(j, j) = 1.0;
470
471	<	} //namespace oopse
471	>	/** @todo */
472	>	v_maxI = v.getRow(maxI);
473	>	v_j = v.getRow(j);
474	>	v_k = cross(v_maxI, v_j);
475	>	v_k.normalize();
476	>	v_j = cross(v_k, v_maxI);
477	>	v.setRow(j, v_j);
478	>	v.setRow(k, v_k);
479	>
480	>
481	>	// transpose vectors back to columns
482	>	v = v.transpose();
483	>	return;
484	>	}
485	>	}
486	>
487	>	// the three eigenvalues are different, just sort the eigenvectors
488	>	// to align them with the x, y, and z axes
489	>
490	>	// find the vector with the largest x element, make that vector
491	>	// the first vector
492	>	maxVal = fabs(v(0, 0));
493	>	maxI = 0;
494	>	for (i = 1; i < 3; i++) {
495	>	if (maxVal < (tmp = fabs(v(i, 0)))) {
496	>	maxVal = tmp;
497	>	maxI = i;
498	>	}
499	>	}
500	>
501	>	// swap eigenvalue and eigenvector
502	>	if (maxI != 0) {
503	>	tmp = w(maxI);
504	>	w(maxI) = w(0);
505	>	w(0) = tmp;
506	>	v.swapRow(maxI, 0);
507	>	}
508	>	// do the same for the y element
509	>	if (fabs(v(1, 1)) < fabs(v(2, 1))) {
510	>	tmp = w(2);
511	>	w(2) = w(1);
512	>	w(1) = tmp;
513	>	v.swapRow(2, 1);
514	>	}
515	>
516	>	// ensure that the sign of the eigenvectors is correct
517	>	for (i = 0; i < 2; i++) {
518	>	if (v(i, i) < 0) {
519	>	v(i, 0) = -v(i, 0);
520	>	v(i, 1) = -v(i, 1);
521	>	v(i, 2) = -v(i, 2);
522	>	}
523	>	}
524	>
525	>	// set sign of final eigenvector to ensure that determinant is positive
526	>	if (v.determinant() < 0) {
527	>	v(2, 0) = -v(2, 0);
528	>	v(2, 1) = -v(2, 1);
529	>	v(2, 2) = -v(2, 2);
530	>	}
531	>
532	>	// transpose the eigenvectors back again
533	>	v = v.transpose();
534	>	return ;
535	>	}
536	>
537	>	/**
538	>	* Return the multiplication of two matrixes  (m1 * m2).
539	>	* @return the multiplication of two matrixes
540	>	* @param m1 the first matrix
541	>	* @param m2 the second matrix
542	>	*/
543	>	template<typename Real>
544	>	inline SquareMatrix3<Real> operator *(const SquareMatrix3<Real>& m1, const SquareMatrix3<Real>& m2) {
545	>	SquareMatrix3<Real> result;
546	>
547	>	for (unsigned int i = 0; i < 3; i++)
548	>	for (unsigned int j = 0; j < 3; j++)
549	>	for (unsigned int k = 0; k < 3; k++)
550	>	result(i, j)  += m1(i, k) * m2(k, j);
551	>
552	>	return result;
553	>	}
554	>
555	>	template<typename Real>
556	>	inline SquareMatrix3<Real> outProduct(const Vector3<Real>& v1, const Vector3<Real>& v2) {
557	>	SquareMatrix3<Real> result;
558	>
559	>	for (unsigned int i = 0; i < 3; i++) {
560	>	for (unsigned int j = 0; j < 3; j++) {
561	>	result(i, j)  = v1[i] * v2[j];
562	>	}
563	>	}
564	>
565	>	return result;
566	>	}
567	>
568	>
569	>	typedef SquareMatrix3<RealType> Mat3x3d;
570	>	typedef SquareMatrix3<RealType> RotMat3x3d;
571	>
572	>	const Mat3x3d M3Zero(0.0);
573	>
574	>
575	>	} //namespace OpenMD
576		#endif // MATH_SQUAREMATRIX_HPP
577	+

Comparing trunk/src/math/SquareMatrix3.hpp (property svn:keywords):
Revision 93 by tim, Sun Oct 17 01:19:11 2004 UTC vs.
Revision 2000 by gezelter, Sat May 31 22:35:05 2014 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines