| 49 |
|
* length 3 vectors |
| 50 |
|
*/ |
| 51 |
|
|
| 52 |
< |
void identityMat3(double A[3][3]) { |
| 52 |
> |
void identityMat3(RealType A[3][3]) { |
| 53 |
|
int i; |
| 54 |
|
for (i = 0; i < 3; i++) { |
| 55 |
|
A[i][0] = A[i][1] = A[i][2] = 0.0; |
| 57 |
|
} |
| 58 |
|
} |
| 59 |
|
|
| 60 |
< |
void swapVectors3(double v1[3], double v2[3]) { |
| 60 |
> |
void swapVectors3(RealType v1[3], RealType v2[3]) { |
| 61 |
|
int i; |
| 62 |
|
for (i = 0; i < 3; i++) { |
| 63 |
< |
double tmp = v1[i]; |
| 63 |
> |
RealType tmp = v1[i]; |
| 64 |
|
v1[i] = v2[i]; |
| 65 |
|
v2[i] = tmp; |
| 66 |
|
} |
| 67 |
|
} |
| 68 |
|
|
| 69 |
< |
double normalize3(double x[3]) { |
| 70 |
< |
double den; |
| 69 |
> |
RealType normalize3(RealType x[3]) { |
| 70 |
> |
RealType den; |
| 71 |
|
int i; |
| 72 |
|
if ( (den = norm3(x)) != 0.0 ) { |
| 73 |
|
for (i=0; i < 3; i++) |
| 78 |
|
return den; |
| 79 |
|
} |
| 80 |
|
|
| 81 |
< |
void matMul3(double a[3][3], double b[3][3], double c[3][3]) { |
| 82 |
< |
double r00, r01, r02, r10, r11, r12, r20, r21, r22; |
| 81 |
> |
void matMul3(RealType a[3][3], RealType b[3][3], RealType c[3][3]) { |
| 82 |
> |
RealType r00, r01, r02, r10, r11, r12, r20, r21, r22; |
| 83 |
|
|
| 84 |
|
r00 = a[0][0]*b[0][0] + a[0][1]*b[1][0] + a[0][2]*b[2][0]; |
| 85 |
|
r01 = a[0][0]*b[0][1] + a[0][1]*b[1][1] + a[0][2]*b[2][1]; |
| 98 |
|
c[2][0] = r20; c[2][1] = r21; c[2][2] = r22; |
| 99 |
|
} |
| 100 |
|
|
| 101 |
< |
void matVecMul3(double m[3][3], double inVec[3], double outVec[3]) { |
| 102 |
< |
double a0, a1, a2; |
| 101 |
> |
void matVecMul3(RealType m[3][3], RealType inVec[3], RealType outVec[3]) { |
| 102 |
> |
RealType a0, a1, a2; |
| 103 |
|
|
| 104 |
|
a0 = inVec[0]; a1 = inVec[1]; a2 = inVec[2]; |
| 105 |
|
|
| 108 |
|
outVec[2] = m[2][0]*a0 + m[2][1]*a1 + m[2][2]*a2; |
| 109 |
|
} |
| 110 |
|
|
| 111 |
< |
double matDet3(double a[3][3]) { |
| 111 |
> |
RealType matDet3(RealType a[3][3]) { |
| 112 |
|
int i, j, k; |
| 113 |
< |
double determinant; |
| 113 |
> |
RealType determinant; |
| 114 |
|
|
| 115 |
|
determinant = 0.0; |
| 116 |
|
|
| 124 |
|
return determinant; |
| 125 |
|
} |
| 126 |
|
|
| 127 |
< |
void invertMat3(double a[3][3], double b[3][3]) { |
| 127 |
> |
void invertMat3(RealType a[3][3], RealType b[3][3]) { |
| 128 |
|
|
| 129 |
|
int i, j, k, l, m, n; |
| 130 |
< |
double determinant; |
| 130 |
> |
RealType determinant; |
| 131 |
|
|
| 132 |
|
determinant = matDet3( a ); |
| 133 |
|
|
| 150 |
|
} |
| 151 |
|
} |
| 152 |
|
|
| 153 |
< |
void transposeMat3(double in[3][3], double out[3][3]) { |
| 154 |
< |
double temp[3][3]; |
| 153 |
> |
void transposeMat3(RealType in[3][3], RealType out[3][3]) { |
| 154 |
> |
RealType temp[3][3]; |
| 155 |
|
int i, j; |
| 156 |
|
|
| 157 |
|
for (i = 0; i < 3; i++) { |
| 166 |
|
} |
| 167 |
|
} |
| 168 |
|
|
| 169 |
< |
void printMat3(double A[3][3] ){ |
| 169 |
> |
void printMat3(RealType A[3][3] ){ |
| 170 |
|
|
| 171 |
|
fprintf(stderr, "[ %g, %g, %g ]\n[ %g, %g, %g ]\n[ %g, %g, %g ]\n", |
| 172 |
|
A[0][0] , A[0][1] , A[0][2], |
| 174 |
|
A[2][0] , A[2][1] , A[2][2]) ; |
| 175 |
|
} |
| 176 |
|
|
| 177 |
< |
void printMat9(double A[9] ){ |
| 177 |
> |
void printMat9(RealType A[9] ){ |
| 178 |
|
|
| 179 |
|
fprintf(stderr, "[ %g, %g, %g ]\n[ %g, %g, %g ]\n[ %g, %g, %g ]\n", |
| 180 |
|
A[0], A[1], A[2], |
| 182 |
|
A[6], A[7], A[8]); |
| 183 |
|
} |
| 184 |
|
|
| 185 |
< |
double matTrace3(double m[3][3]){ |
| 186 |
< |
double trace; |
| 185 |
> |
RealType matTrace3(RealType m[3][3]){ |
| 186 |
> |
RealType trace; |
| 187 |
|
trace = m[0][0] + m[1][1] + m[2][2]; |
| 188 |
|
|
| 189 |
|
return trace; |
| 190 |
|
} |
| 191 |
|
|
| 192 |
< |
void crossProduct3(double a[3],double b[3], double out[3]){ |
| 192 |
> |
void crossProduct3(RealType a[3],RealType b[3], RealType out[3]){ |
| 193 |
|
|
| 194 |
|
out[0] = a[1] * b[2] - a[2] * b[1]; |
| 195 |
|
out[1] = a[2] * b[0] - a[0] * b[2] ; |
| 197 |
|
|
| 198 |
|
} |
| 199 |
|
|
| 200 |
< |
double dotProduct3(double a[3], double b[3]){ |
| 200 |
> |
RealType dotProduct3(RealType a[3], RealType b[3]){ |
| 201 |
|
return a[0]*b[0] + a[1]*b[1]+ a[2]*b[2]; |
| 202 |
|
} |
| 203 |
|
|
| 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]) { |
| 210 |
> |
void diagonalize3x3(const RealType A[3][3], RealType w[3], RealType V[3][3]) { |
| 211 |
|
int i,j,k,maxI; |
| 212 |
< |
double tmp, maxVal; |
| 212 |
> |
RealType tmp, maxVal; |
| 213 |
|
|
| 214 |
|
/* do the matrix[3][3] to **matrix conversion for Jacobi*/ |
| 215 |
< |
double C[3][3]; |
| 216 |
< |
double *ATemp[3],*VTemp[3]; |
| 215 |
> |
RealType C[3][3]; |
| 216 |
> |
RealType *ATemp[3],*VTemp[3]; |
| 217 |
|
for (i = 0; i < 3; i++) |
| 218 |
|
{ |
| 219 |
|
C[i][0] = A[i][0]; |
| 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) { |
| 354 |
> |
int JacobiN(RealType **a, int n, RealType *w, RealType **v) { |
| 355 |
|
|
| 356 |
|
int i, j, k, iq, ip, numPos; |
| 357 |
|
int ceil_half_n; |
| 358 |
< |
double tresh, theta, tau, t, sm, s, h, g, c, tmp; |
| 359 |
< |
double bspace[4], zspace[4]; |
| 360 |
< |
double *b = bspace; |
| 361 |
< |
double *z = zspace; |
| 358 |
> |
RealType tresh, theta, tau, t, sm, s, h, g, c, tmp; |
| 359 |
> |
RealType bspace[4], zspace[4]; |
| 360 |
> |
RealType *b = bspace; |
| 361 |
> |
RealType *z = zspace; |
| 362 |
|
|
| 363 |
|
|
| 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)); |
| 367 |
> |
b = (RealType *) calloc(n, sizeof(RealType)); |
| 368 |
> |
z = (RealType *) calloc(n, sizeof(RealType)); |
| 369 |
|
} |
| 370 |
|
|
| 371 |
|
/* initialize*/ |
| 526 |
|
numPos++; |
| 527 |
|
} |
| 528 |
|
} |
| 529 |
< |
/* if ( numPos < ceil(double(n)/double(2.0)) )*/ |
| 529 |
> |
/* if ( numPos < ceil(RealType(n)/RealType(2.0)) )*/ |
| 530 |
|
if ( numPos < ceil_half_n) |
| 531 |
|
{ |
| 532 |
|
for(i=0; i<n; i++) |