# | 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. Acknowledgement of the program authors must be made in any |
10 | + | * publication of scientific results based in part on use of the |
11 | + | * program. An acceptable form of acknowledgement is citation of |
12 | + | * the article in which the program was described (Matthew |
13 | + | * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher |
14 | + | * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented |
15 | + | * Parallel Simulation Engine for Molecular Dynamics," |
16 | + | * J. Comput. Chem. 26, pp. 252-271 (2005)) |
17 | + | * |
18 | + | * 2. Redistributions of source code must retain the above copyright |
19 | + | * notice, this list of conditions and the following disclaimer. |
20 | + | * |
21 | + | * 3. Redistributions in binary form must reproduce the above copyright |
22 | + | * notice, this list of conditions and the following disclaimer in the |
23 | + | * documentation and/or other materials provided with the |
24 | + | * distribution. |
25 | + | * |
26 | + | * This software is provided "AS IS," without a warranty of any |
27 | + | * kind. All express or implied conditions, representations and |
28 | + | * warranties, including any implied warranty of merchantability, |
29 | + | * fitness for a particular purpose or non-infringement, are hereby |
30 | + | * excluded. The University of Notre Dame and its licensors shall not |
31 | + | * be liable for any damages suffered by licensee as a result of |
32 | + | * using, modifying or distributing the software or its |
33 | + | * derivatives. In no event will the University of Notre Dame or its |
34 | + | * licensors be liable for any lost revenue, profit or data, or for |
35 | + | * direct, indirect, special, consequential, incidental or punitive |
36 | + | * damages, however caused and regardless of the theory of liability, |
37 | + | * arising out of the use of or inability to use software, even if the |
38 | + | * University of Notre Dame has been advised of the possibility of |
39 | + | * such damages. |
40 | */ | |
41 | < | |
41 | > | |
42 | /** | |
43 | * @file SquareMatrix.hpp | |
44 | * @author Teng Lin | |
# | Line 36 | Line 52 | namespace oopse { | |
52 | ||
53 | namespace oopse { | |
54 | ||
55 | < | /** |
56 | < | * @class SquareMatrix SquareMatrix.hpp "math/SquareMatrix.hpp" |
57 | < | * @brief A square matrix class |
58 | < | * @template Real the element type |
59 | < | * @template Dim the dimension of the square matrix |
60 | < | */ |
61 | < | template<typename Real, int Dim> |
62 | < | class SquareMatrix : public RectMatrix<Real, Dim, Dim> { |
63 | < | public: |
55 | > | /** |
56 | > | * @class SquareMatrix SquareMatrix.hpp "math/SquareMatrix.hpp" |
57 | > | * @brief A square matrix class |
58 | > | * @template Real the element type |
59 | > | * @template Dim the dimension of the square matrix |
60 | > | */ |
61 | > | template<typename Real, int Dim> |
62 | > | class SquareMatrix : public RectMatrix<Real, Dim, Dim> { |
63 | > | public: |
64 | > | typedef Real ElemType; |
65 | > | typedef Real* ElemPoinerType; |
66 | ||
67 | < | /** default constructor */ |
68 | < | SquareMatrix() { |
69 | < | for (unsigned int i = 0; i < Dim; i++) |
70 | < | for (unsigned int j = 0; j < Dim; j++) |
71 | < | data_[i][j] = 0.0; |
72 | < | } |
67 | > | /** default constructor */ |
68 | > | SquareMatrix() { |
69 | > | for (unsigned int i = 0; i < Dim; i++) |
70 | > | for (unsigned int j = 0; j < Dim; j++) |
71 | > | this->data_[i][j] = 0.0; |
72 | > | } |
73 | ||
74 | < | /** copy constructor */ |
75 | < | SquareMatrix(const RectMatrix<Real, Dim, Dim>& m) : RectMatrix<Real, Dim, Dim>(m) { |
76 | < | } |
59 | < | |
60 | < | /** copy assignment operator */ |
61 | < | SquareMatrix<Real, Dim>& operator =(const RectMatrix<Real, Dim, Dim>& m) { |
62 | < | RectMatrix<Real, Dim, Dim>::operator=(m); |
63 | < | return *this; |
64 | < | } |
65 | < | |
66 | < | /** Retunrs an identity matrix*/ |
74 | > | /** Constructs and initializes every element of this matrix to a scalar */ |
75 | > | SquareMatrix(Real s) : RectMatrix<Real, Dim, Dim>(s){ |
76 | > | } |
77 | ||
78 | < | static SquareMatrix<Real, Dim> identity() { |
79 | < | SquareMatrix<Real, Dim> m; |
80 | < | |
71 | < | for (unsigned int i = 0; i < Dim; i++) |
72 | < | for (unsigned int j = 0; j < Dim; j++) |
73 | < | if (i == j) |
74 | < | m(i, j) = 1.0; |
75 | < | else |
76 | < | m(i, j) = 0.0; |
78 | > | /** Constructs and initializes from an array */ |
79 | > | SquareMatrix(Real* array) : RectMatrix<Real, Dim, Dim>(array){ |
80 | > | } |
81 | ||
78 | – | return m; |
79 | – | } |
82 | ||
83 | < | /** Retunrs the inversion of this matrix. */ |
84 | < | SquareMatrix<Real, Dim> inverse() { |
85 | < | SquareMatrix<Real, Dim> result; |
83 | > | /** copy constructor */ |
84 | > | SquareMatrix(const RectMatrix<Real, Dim, Dim>& m) : RectMatrix<Real, Dim, Dim>(m) { |
85 | > | } |
86 | > | |
87 | > | /** copy assignment operator */ |
88 | > | SquareMatrix<Real, Dim>& operator =(const RectMatrix<Real, Dim, Dim>& m) { |
89 | > | RectMatrix<Real, Dim, Dim>::operator=(m); |
90 | > | return *this; |
91 | > | } |
92 | > | |
93 | > | /** Retunrs an identity matrix*/ |
94 | ||
95 | < | return result; |
96 | < | } |
95 | > | static SquareMatrix<Real, Dim> identity() { |
96 | > | SquareMatrix<Real, Dim> m; |
97 | > | |
98 | > | for (unsigned int i = 0; i < Dim; i++) |
99 | > | for (unsigned int j = 0; j < Dim; j++) |
100 | > | if (i == j) |
101 | > | m(i, j) = 1.0; |
102 | > | else |
103 | > | m(i, j) = 0.0; |
104 | ||
105 | < | |
105 | > | return m; |
106 | > | } |
107 | ||
108 | < | /** Returns the determinant of this matrix. */ |
109 | < | double determinant() const { |
110 | < | double det; |
111 | < | return det; |
112 | < | } |
108 | > | /** |
109 | > | * Retunrs the inversion of this matrix. |
110 | > | * @todo need implementation |
111 | > | */ |
112 | > | SquareMatrix<Real, Dim> inverse() { |
113 | > | SquareMatrix<Real, Dim> result; |
114 | ||
115 | < | /** Returns the trace of this matrix. */ |
116 | < | double trace() const { |
98 | < | double tmp = 0; |
99 | < | |
100 | < | for (unsigned int i = 0; i < Dim ; i++) |
101 | < | tmp += data_[i][i]; |
115 | > | return result; |
116 | > | } |
117 | ||
118 | < | return tmp; |
119 | < | } |
118 | > | /** |
119 | > | * Returns the determinant of this matrix. |
120 | > | * @todo need implementation |
121 | > | */ |
122 | > | Real determinant() const { |
123 | > | Real det; |
124 | > | return det; |
125 | > | } |
126 | ||
127 | < | /** Tests if this matrix is symmetrix. */ |
128 | < | bool isSymmetric() const { |
129 | < | for (unsigned int i = 0; i < Dim - 1; i++) |
130 | < | for (unsigned int j = i; j < Dim; j++) |
131 | < | if (fabs(data_[i][j] - data_[j][i]) > oopse::epsilon) |
132 | < | return false; |
112 | < | |
113 | < | return true; |
114 | < | } |
127 | > | /** Returns the trace of this matrix. */ |
128 | > | Real trace() const { |
129 | > | Real tmp = 0; |
130 | > | |
131 | > | for (unsigned int i = 0; i < Dim ; i++) |
132 | > | tmp += this->data_[i][i]; |
133 | ||
134 | < | /** Tests if this matrix is orthogona. */ |
135 | < | bool isOrthogonal() { |
118 | < | SquareMatrix<Real, Dim> tmp; |
134 | > | return tmp; |
135 | > | } |
136 | ||
137 | < | tmp = *this * transpose(); |
137 | > | /** Tests if this matrix is symmetrix. */ |
138 | > | bool isSymmetric() const { |
139 | > | for (unsigned int i = 0; i < Dim - 1; i++) |
140 | > | for (unsigned int j = i; j < Dim; j++) |
141 | > | if (fabs(this->data_[i][j] - this->data_[j][i]) > oopse::epsilon) |
142 | > | return false; |
143 | > | |
144 | > | return true; |
145 | > | } |
146 | ||
147 | < | return tmp.isUnitMatrix(); |
148 | < | } |
147 | > | /** Tests if this matrix is orthogonal. */ |
148 | > | bool isOrthogonal() { |
149 | > | SquareMatrix<Real, Dim> tmp; |
150 | ||
151 | < | /** Tests if this matrix is diagonal. */ |
126 | < | bool isDiagonal() const { |
127 | < | for (unsigned int i = 0; i < Dim ; i++) |
128 | < | for (unsigned int j = 0; j < Dim; j++) |
129 | < | if (i !=j && fabs(data_[i][j]) > oopse::epsilon) |
130 | < | return false; |
131 | < | |
132 | < | return true; |
133 | < | } |
151 | > | tmp = *this * transpose(); |
152 | ||
153 | < | /** Tests if this matrix is the unit matrix. */ |
154 | < | bool isUnitMatrix() const { |
137 | < | if (!isDiagonal()) |
138 | < | return false; |
139 | < | |
140 | < | for (unsigned int i = 0; i < Dim ; i++) |
141 | < | if (fabs(data_[i][i] - 1) > oopse::epsilon) |
142 | < | return false; |
143 | < | |
144 | < | return true; |
145 | < | } |
153 | > | return tmp.isDiagonal(); |
154 | > | } |
155 | ||
156 | < | };//end SquareMatrix |
156 | > | /** Tests if this matrix is diagonal. */ |
157 | > | bool isDiagonal() const { |
158 | > | for (unsigned int i = 0; i < Dim ; i++) |
159 | > | for (unsigned int j = 0; j < Dim; j++) |
160 | > | if (i !=j && fabs(this->data_[i][j]) > oopse::epsilon) |
161 | > | return false; |
162 | > | |
163 | > | return true; |
164 | > | } |
165 | ||
166 | + | /** Tests if this matrix is the unit matrix. */ |
167 | + | bool isUnitMatrix() const { |
168 | + | if (!isDiagonal()) |
169 | + | return false; |
170 | + | |
171 | + | for (unsigned int i = 0; i < Dim ; i++) |
172 | + | if (fabs(this->data_[i][i] - 1) > oopse::epsilon) |
173 | + | return false; |
174 | + | |
175 | + | return true; |
176 | + | } |
177 | + | |
178 | + | /** Return the transpose of this matrix */ |
179 | + | SquareMatrix<Real, Dim> transpose() const{ |
180 | + | SquareMatrix<Real, Dim> result; |
181 | + | |
182 | + | for (unsigned int i = 0; i < Dim; i++) |
183 | + | for (unsigned int j = 0; j < Dim; j++) |
184 | + | result(j, i) = this->data_[i][j]; |
185 | + | |
186 | + | return result; |
187 | + | } |
188 | + | |
189 | + | /** @todo need implementation */ |
190 | + | void diagonalize() { |
191 | + | //jacobi(m, eigenValues, ortMat); |
192 | + | } |
193 | + | |
194 | + | /** |
195 | + | * Jacobi iteration routines for computing eigenvalues/eigenvectors of |
196 | + | * real symmetric matrix |
197 | + | * |
198 | + | * @return true if success, otherwise return false |
199 | + | * @param a symmetric matrix whose eigenvectors are to be computed. On return, the matrix is |
200 | + | * overwritten |
201 | + | * @param w will contain the eigenvalues of the matrix On return of this function |
202 | + | * @param v the columns of this matrix will contain the eigenvectors. The eigenvectors are |
203 | + | * normalized and mutually orthogonal. |
204 | + | */ |
205 | + | |
206 | + | static int jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& d, |
207 | + | SquareMatrix<Real, Dim>& v); |
208 | + | };//end SquareMatrix |
209 | + | |
210 | + | |
211 | + | /*========================================================================= |
212 | + | |
213 | + | Program: Visualization Toolkit |
214 | + | Module: $RCSfile: SquareMatrix.hpp,v $ |
215 | + | |
216 | + | Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen |
217 | + | All rights reserved. |
218 | + | See Copyright.txt or http://www.kitware.com/Copyright.htm for details. |
219 | + | |
220 | + | This software is distributed WITHOUT ANY WARRANTY; without even |
221 | + | the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR |
222 | + | PURPOSE. See the above copyright notice for more information. |
223 | + | |
224 | + | =========================================================================*/ |
225 | + | |
226 | + | #define VTK_ROTATE(a,i,j,k,l) g=a(i, j);h=a(k, l);a(i, j)=g-s*(h+g*tau); \ |
227 | + | a(k, l)=h+s*(g-h*tau) |
228 | + | |
229 | + | #define VTK_MAX_ROTATIONS 20 |
230 | + | |
231 | + | // Jacobi iteration for the solution of eigenvectors/eigenvalues of a nxn |
232 | + | // real symmetric matrix. Square nxn matrix a; size of matrix in n; |
233 | + | // output eigenvalues in w; and output eigenvectors in v. Resulting |
234 | + | // eigenvalues/vectors are sorted in decreasing order; eigenvectors are |
235 | + | // normalized. |
236 | + | template<typename Real, int Dim> |
237 | + | int SquareMatrix<Real, Dim>::jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w, |
238 | + | SquareMatrix<Real, Dim>& v) { |
239 | + | const int n = Dim; |
240 | + | int i, j, k, iq, ip, numPos; |
241 | + | Real tresh, theta, tau, t, sm, s, h, g, c, tmp; |
242 | + | Real bspace[4], zspace[4]; |
243 | + | Real *b = bspace; |
244 | + | Real *z = zspace; |
245 | + | |
246 | + | // only allocate memory if the matrix is large |
247 | + | if (n > 4) { |
248 | + | b = new Real[n]; |
249 | + | z = new Real[n]; |
250 | + | } |
251 | + | |
252 | + | // initialize |
253 | + | for (ip=0; ip<n; ip++) { |
254 | + | for (iq=0; iq<n; iq++) { |
255 | + | v(ip, iq) = 0.0; |
256 | + | } |
257 | + | v(ip, ip) = 1.0; |
258 | + | } |
259 | + | for (ip=0; ip<n; ip++) { |
260 | + | b[ip] = w[ip] = a(ip, ip); |
261 | + | z[ip] = 0.0; |
262 | + | } |
263 | + | |
264 | + | // begin rotation sequence |
265 | + | for (i=0; i<VTK_MAX_ROTATIONS; i++) { |
266 | + | sm = 0.0; |
267 | + | for (ip=0; ip<n-1; ip++) { |
268 | + | for (iq=ip+1; iq<n; iq++) { |
269 | + | sm += fabs(a(ip, iq)); |
270 | + | } |
271 | + | } |
272 | + | if (sm == 0.0) { |
273 | + | break; |
274 | + | } |
275 | + | |
276 | + | if (i < 3) { // first 3 sweeps |
277 | + | tresh = 0.2*sm/(n*n); |
278 | + | } else { |
279 | + | tresh = 0.0; |
280 | + | } |
281 | + | |
282 | + | for (ip=0; ip<n-1; ip++) { |
283 | + | for (iq=ip+1; iq<n; iq++) { |
284 | + | g = 100.0*fabs(a(ip, iq)); |
285 | + | |
286 | + | // after 4 sweeps |
287 | + | if (i > 3 && (fabs(w[ip])+g) == fabs(w[ip]) |
288 | + | && (fabs(w[iq])+g) == fabs(w[iq])) { |
289 | + | a(ip, iq) = 0.0; |
290 | + | } else if (fabs(a(ip, iq)) > tresh) { |
291 | + | h = w[iq] - w[ip]; |
292 | + | if ( (fabs(h)+g) == fabs(h)) { |
293 | + | t = (a(ip, iq)) / h; |
294 | + | } else { |
295 | + | theta = 0.5*h / (a(ip, iq)); |
296 | + | t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta)); |
297 | + | if (theta < 0.0) { |
298 | + | t = -t; |
299 | + | } |
300 | + | } |
301 | + | c = 1.0 / sqrt(1+t*t); |
302 | + | s = t*c; |
303 | + | tau = s/(1.0+c); |
304 | + | h = t*a(ip, iq); |
305 | + | z[ip] -= h; |
306 | + | z[iq] += h; |
307 | + | w[ip] -= h; |
308 | + | w[iq] += h; |
309 | + | a(ip, iq)=0.0; |
310 | + | |
311 | + | // ip already shifted left by 1 unit |
312 | + | for (j = 0;j <= ip-1;j++) { |
313 | + | VTK_ROTATE(a,j,ip,j,iq); |
314 | + | } |
315 | + | // ip and iq already shifted left by 1 unit |
316 | + | for (j = ip+1;j <= iq-1;j++) { |
317 | + | VTK_ROTATE(a,ip,j,j,iq); |
318 | + | } |
319 | + | // iq already shifted left by 1 unit |
320 | + | for (j=iq+1; j<n; j++) { |
321 | + | VTK_ROTATE(a,ip,j,iq,j); |
322 | + | } |
323 | + | for (j=0; j<n; j++) { |
324 | + | VTK_ROTATE(v,j,ip,j,iq); |
325 | + | } |
326 | + | } |
327 | + | } |
328 | + | } |
329 | + | |
330 | + | for (ip=0; ip<n; ip++) { |
331 | + | b[ip] += z[ip]; |
332 | + | w[ip] = b[ip]; |
333 | + | z[ip] = 0.0; |
334 | + | } |
335 | + | } |
336 | + | |
337 | + | //// this is NEVER called |
338 | + | if ( i >= VTK_MAX_ROTATIONS ) { |
339 | + | std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions" << std::endl; |
340 | + | return 0; |
341 | + | } |
342 | + | |
343 | + | // sort eigenfunctions these changes do not affect accuracy |
344 | + | for (j=0; j<n-1; j++) { // boundary incorrect |
345 | + | k = j; |
346 | + | tmp = w[k]; |
347 | + | for (i=j+1; i<n; i++) { // boundary incorrect, shifted already |
348 | + | if (w[i] >= tmp) { // why exchage if same? |
349 | + | k = i; |
350 | + | tmp = w[k]; |
351 | + | } |
352 | + | } |
353 | + | if (k != j) { |
354 | + | w[k] = w[j]; |
355 | + | w[j] = tmp; |
356 | + | for (i=0; i<n; i++) { |
357 | + | tmp = v(i, j); |
358 | + | v(i, j) = v(i, k); |
359 | + | v(i, k) = tmp; |
360 | + | } |
361 | + | } |
362 | + | } |
363 | + | // insure eigenvector consistency (i.e., Jacobi can compute vectors that |
364 | + | // are negative of one another (.707,.707,0) and (-.707,-.707,0). This can |
365 | + | // reek havoc in hyperstreamline/other stuff. We will select the most |
366 | + | // positive eigenvector. |
367 | + | int ceil_half_n = (n >> 1) + (n & 1); |
368 | + | for (j=0; j<n; j++) { |
369 | + | for (numPos=0, i=0; i<n; i++) { |
370 | + | if ( v(i, j) >= 0.0 ) { |
371 | + | numPos++; |
372 | + | } |
373 | + | } |
374 | + | // if ( numPos < ceil(double(n)/double(2.0)) ) |
375 | + | if ( numPos < ceil_half_n) { |
376 | + | for (i=0; i<n; i++) { |
377 | + | v(i, j) *= -1.0; |
378 | + | } |
379 | + | } |
380 | + | } |
381 | + | |
382 | + | if (n > 4) { |
383 | + | delete [] b; |
384 | + | delete [] z; |
385 | + | } |
386 | + | return 1; |
387 | + | } |
388 | + | |
389 | + | |
390 | + | typedef SquareMatrix<double, 6> Mat6x6d; |
391 | } | |
392 | #endif //MATH_SQUAREMATRIX_HPP | |
393 | + |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |