1 |
+ |
/* |
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 |
+ |
|
42 |
|
#include <stdio.h> |
43 |
|
#include <math.h> |
44 |
|
#include <stdlib.h> |
45 |
< |
#include "MatVec3.h" |
45 |
> |
#include "math/MatVec3.h" |
46 |
|
|
47 |
|
/* |
48 |
|
* Contains various utilities for dealing with 3x3 matrices and |
201 |
|
return a[0]*b[0] + a[1]*b[1]+ a[2]*b[2]; |
202 |
|
} |
203 |
|
|
204 |
< |
//---------------------------------------------------------------------------- |
205 |
< |
// Extract the eigenvalues and eigenvectors from a 3x3 matrix. |
206 |
< |
// The eigenvectors (the columns of V) will be normalized. |
207 |
< |
// The eigenvectors are aligned optimally with the x, y, and z |
208 |
< |
// axes respectively. |
204 |
> |
/*----------------------------------------------------------------------------*/ |
205 |
> |
/* Extract the eigenvalues and eigenvectors from a 3x3 matrix.*/ |
206 |
> |
/* The eigenvectors (the columns of V) will be normalized. */ |
207 |
> |
/* The eigenvectors are aligned optimally with the x, y, and z*/ |
208 |
> |
/* axes respectively.*/ |
209 |
|
|
210 |
|
void diagonalize3x3(const double A[3][3], double w[3], double V[3][3]) { |
211 |
|
int i,j,k,maxI; |
212 |
|
double tmp, maxVal; |
213 |
|
|
214 |
< |
// do the matrix[3][3] to **matrix conversion for Jacobi |
214 |
> |
/* do the matrix[3][3] to **matrix conversion for Jacobi*/ |
215 |
|
double C[3][3]; |
216 |
|
double *ATemp[3],*VTemp[3]; |
217 |
|
for (i = 0; i < 3; i++) |
223 |
|
VTemp[i] = V[i]; |
224 |
|
} |
225 |
|
|
226 |
< |
// diagonalize using Jacobi |
226 |
> |
/* diagonalize using Jacobi*/ |
227 |
|
JacobiN(ATemp,3,w,VTemp); |
228 |
|
|
229 |
< |
// if all the eigenvalues are the same, return identity matrix |
229 |
> |
/* if all the eigenvalues are the same, return identity matrix*/ |
230 |
|
if (w[0] == w[1] && w[0] == w[2]) |
231 |
|
{ |
232 |
|
identityMat3(V); |
233 |
|
return; |
234 |
|
} |
235 |
|
|
236 |
< |
// transpose temporarily, it makes it easier to sort the eigenvectors |
236 |
> |
/* transpose temporarily, it makes it easier to sort the eigenvectors*/ |
237 |
|
transposeMat3(V,V); |
238 |
|
|
239 |
< |
// if two eigenvalues are the same, re-orthogonalize to optimally line |
240 |
< |
// up the eigenvectors with the x, y, and z axes |
239 |
> |
/* if two eigenvalues are the same, re-orthogonalize to optimally line*/ |
240 |
> |
/* up the eigenvectors with the x, y, and z axes*/ |
241 |
|
for (i = 0; i < 3; i++) |
242 |
|
{ |
243 |
< |
if (w[(i+1)%3] == w[(i+2)%3]) // two eigenvalues are the same |
243 |
> |
if (w[(i+1)%3] == w[(i+2)%3]) /* two eigenvalues are the same*/ |
244 |
|
{ |
245 |
< |
// find maximum element of the independant eigenvector |
245 |
> |
/* find maximum element of the independant eigenvector*/ |
246 |
|
maxVal = fabs(V[i][0]); |
247 |
|
maxI = 0; |
248 |
|
for (j = 1; j < 3; j++) |
253 |
|
maxI = j; |
254 |
|
} |
255 |
|
} |
256 |
< |
// swap the eigenvector into its proper position |
256 |
> |
/* swap the eigenvector into its proper position*/ |
257 |
|
if (maxI != i) |
258 |
|
{ |
259 |
|
tmp = w[maxI]; |
261 |
|
w[i] = tmp; |
262 |
|
swapVectors3(V[i],V[maxI]); |
263 |
|
} |
264 |
< |
// maximum element of eigenvector should be positive |
264 |
> |
/* maximum element of eigenvector should be positive*/ |
265 |
|
if (V[maxI][maxI] < 0) |
266 |
|
{ |
267 |
|
V[maxI][0] = -V[maxI][0]; |
269 |
|
V[maxI][2] = -V[maxI][2]; |
270 |
|
} |
271 |
|
|
272 |
< |
// re-orthogonalize the other two eigenvectors |
272 |
> |
/* re-orthogonalize the other two eigenvectors*/ |
273 |
|
j = (maxI+1)%3; |
274 |
|
k = (maxI+2)%3; |
275 |
|
|
281 |
|
normalize3(V[k]); |
282 |
|
crossProduct3(V[k],V[maxI],V[j]); |
283 |
|
|
284 |
< |
// transpose vectors back to columns |
284 |
> |
/* transpose vectors back to columns*/ |
285 |
|
transposeMat3(V,V); |
286 |
|
return; |
287 |
|
} |
288 |
|
} |
289 |
|
|
290 |
< |
// the three eigenvalues are different, just sort the eigenvectors |
291 |
< |
// to align them with the x, y, and z axes |
290 |
> |
/* the three eigenvalues are different, just sort the eigenvectors*/ |
291 |
> |
/* to align them with the x, y, and z axes*/ |
292 |
|
|
293 |
< |
// find the vector with the largest x element, make that vector |
294 |
< |
// the first vector |
293 |
> |
/* find the vector with the largest x element, make that vector*/ |
294 |
> |
/* the first vector*/ |
295 |
|
maxVal = fabs(V[0][0]); |
296 |
|
maxI = 0; |
297 |
|
for (i = 1; i < 3; i++) |
302 |
|
maxI = i; |
303 |
|
} |
304 |
|
} |
305 |
< |
// swap eigenvalue and eigenvector |
305 |
> |
/* swap eigenvalue and eigenvector*/ |
306 |
|
if (maxI != 0) |
307 |
|
{ |
308 |
|
tmp = w[maxI]; |
310 |
|
w[0] = tmp; |
311 |
|
swapVectors3(V[maxI],V[0]); |
312 |
|
} |
313 |
< |
// do the same for the y element |
313 |
> |
/* do the same for the y element*/ |
314 |
|
if (fabs(V[1][1]) < fabs(V[2][1])) |
315 |
|
{ |
316 |
|
tmp = w[2]; |
319 |
|
swapVectors3(V[2],V[1]); |
320 |
|
} |
321 |
|
|
322 |
< |
// ensure that the sign of the eigenvectors is correct |
322 |
> |
/* ensure that the sign of the eigenvectors is correct*/ |
323 |
|
for (i = 0; i < 2; i++) |
324 |
|
{ |
325 |
|
if (V[i][i] < 0) |
329 |
|
V[i][2] = -V[i][2]; |
330 |
|
} |
331 |
|
} |
332 |
< |
// set sign of final eigenvector to ensure that determinant is positive |
332 |
> |
/* set sign of final eigenvector to ensure that determinant is positive*/ |
333 |
|
if (matDet3(V) < 0) |
334 |
|
{ |
335 |
|
V[2][0] = -V[2][0]; |
337 |
|
V[2][2] = -V[2][2]; |
338 |
|
} |
339 |
|
|
340 |
< |
// transpose the eigenvectors back again |
340 |
> |
/* transpose the eigenvectors back again*/ |
341 |
|
transposeMat3(V,V); |
342 |
|
} |
343 |
|
|
346 |
|
|
347 |
|
#define MAX_ROTATIONS 20 |
348 |
|
|
349 |
< |
// Jacobi iteration for the solution of eigenvectors/eigenvalues of a nxn |
350 |
< |
// real symmetric matrix. Square nxn matrix a; size of matrix in n; |
351 |
< |
// output eigenvalues in w; and output eigenvectors in v. Resulting |
352 |
< |
// eigenvalues/vectors are sorted in decreasing order; eigenvectors are |
353 |
< |
// normalized. |
349 |
> |
/* Jacobi iteration for the solution of eigenvectors/eigenvalues of a nxn*/ |
350 |
> |
/* real symmetric matrix. Square nxn matrix a; size of matrix in n;*/ |
351 |
> |
/* output eigenvalues in w; and output eigenvectors in v. Resulting*/ |
352 |
> |
/* eigenvalues/vectors are sorted in decreasing order; eigenvectors are*/ |
353 |
> |
/* normalized.*/ |
354 |
|
int JacobiN(double **a, int n, double *w, double **v) { |
355 |
|
|
356 |
|
int i, j, k, iq, ip, numPos; |
361 |
|
double *z = zspace; |
362 |
|
|
363 |
|
|
364 |
< |
// only allocate memory if the matrix is large |
364 |
> |
/* only allocate memory if the matrix is large*/ |
365 |
|
if (n > 4) |
366 |
|
{ |
367 |
|
b = (double *) calloc(n, sizeof(double)); |
368 |
|
z = (double *) calloc(n, sizeof(double)); |
369 |
|
} |
370 |
|
|
371 |
< |
// initialize |
371 |
> |
/* initialize*/ |
372 |
|
for (ip=0; ip<n; ip++) |
373 |
|
{ |
374 |
|
for (iq=0; iq<n; iq++) |
383 |
|
z[ip] = 0.0; |
384 |
|
} |
385 |
|
|
386 |
< |
// begin rotation sequence |
386 |
> |
/* begin rotation sequence*/ |
387 |
|
for (i=0; i<MAX_ROTATIONS; i++) |
388 |
|
{ |
389 |
|
sm = 0.0; |
399 |
|
break; |
400 |
|
} |
401 |
|
|
402 |
< |
if (i < 3) // first 3 sweeps |
402 |
> |
if (i < 3) /* first 3 sweeps*/ |
403 |
|
{ |
404 |
|
tresh = 0.2*sm/(n*n); |
405 |
|
} |
414 |
|
{ |
415 |
|
g = 100.0*fabs(a[ip][iq]); |
416 |
|
|
417 |
< |
// after 4 sweeps |
417 |
> |
/* after 4 sweeps*/ |
418 |
|
if (i > 3 && (fabs(w[ip])+g) == fabs(w[ip]) |
419 |
|
&& (fabs(w[iq])+g) == fabs(w[iq])) |
420 |
|
{ |
446 |
|
w[iq] += h; |
447 |
|
a[ip][iq]=0.0; |
448 |
|
|
449 |
< |
// ip already shifted left by 1 unit |
449 |
> |
/* ip already shifted left by 1 unit*/ |
450 |
|
for (j = 0;j <= ip-1;j++) |
451 |
|
{ |
452 |
|
MAT_ROTATE(a,j,ip,j,iq) |
453 |
|
} |
454 |
< |
// ip and iq already shifted left by 1 unit |
454 |
> |
/* ip and iq already shifted left by 1 unit*/ |
455 |
|
for (j = ip+1;j <= iq-1;j++) |
456 |
|
{ |
457 |
|
MAT_ROTATE(a,ip,j,j,iq) |
458 |
|
} |
459 |
< |
// iq already shifted left by 1 unit |
459 |
> |
/* iq already shifted left by 1 unit*/ |
460 |
|
for (j=iq+1; j<n; j++) |
461 |
|
{ |
462 |
|
MAT_ROTATE(a,ip,j,iq,j) |
477 |
|
} |
478 |
|
} |
479 |
|
|
480 |
< |
//// this is NEVER called |
480 |
> |
/*// this is NEVER called*/ |
481 |
|
if ( i >= MAX_ROTATIONS ) |
482 |
|
{ |
483 |
|
sprintf( painCave.errMsg, |
487 |
|
return 0; |
488 |
|
} |
489 |
|
|
490 |
< |
// sort eigenfunctions these changes do not affect accuracy |
491 |
< |
for (j=0; j<n-1; j++) // boundary incorrect |
490 |
> |
/* sort eigenfunctions these changes do not affect accuracy */ |
491 |
> |
for (j=0; j<n-1; j++) /* boundary incorrect*/ |
492 |
|
{ |
493 |
|
k = j; |
494 |
|
tmp = w[k]; |
495 |
< |
for (i=j+1; i<n; i++) // boundary incorrect, shifted already |
495 |
> |
for (i=j+1; i<n; i++) /* boundary incorrect, shifted already*/ |
496 |
|
{ |
497 |
< |
if (w[i] >= tmp) // why exchage if same? |
497 |
> |
if (w[i] >= tmp) /* why exchage if same?*/ |
498 |
|
{ |
499 |
|
k = i; |
500 |
|
tmp = w[k]; |
512 |
|
} |
513 |
|
} |
514 |
|
} |
515 |
< |
// insure eigenvector consistency (i.e., Jacobi can compute vectors that |
516 |
< |
// are negative of one another (.707,.707,0) and (-.707,-.707,0). This can |
517 |
< |
// reek havoc in hyperstreamline/other stuff. We will select the most |
518 |
< |
// positive eigenvector. |
515 |
> |
/* insure eigenvector consistency (i.e., Jacobi can compute vectors that*/ |
516 |
> |
/* are negative of one another (.707,.707,0) and (-.707,-.707,0). This can*/ |
517 |
> |
/* reek havoc in hyperstreamline/other stuff. We will select the most*/ |
518 |
> |
/* positive eigenvector.*/ |
519 |
|
ceil_half_n = (n >> 1) + (n & 1); |
520 |
|
for (j=0; j<n; j++) |
521 |
|
{ |
526 |
|
numPos++; |
527 |
|
} |
528 |
|
} |
529 |
< |
// if ( numPos < ceil(double(n)/double(2.0)) ) |
529 |
> |
/* if ( numPos < ceil(double(n)/double(2.0)) )*/ |
530 |
|
if ( numPos < ceil_half_n) |
531 |
|
{ |
532 |
|
for(i=0; i<n; i++) |