1 | #include <math.h> |
---|
2 | #include <stdio.h> |
---|
3 | #include <stdlib.h> |
---|
4 | #define ROTATE(a,i,j,k,l) g=a[i*n + j];h=a[k*n + l];a[i*n + j]=g-s*(h+g*tau);\ |
---|
5 | a[k*n + l]=h+s*(g-h*tau); |
---|
6 | |
---|
7 | /* |
---|
8 | Computes all eigenvalues and eigenvectors of a real symmetric matrix a[1..n][1..n]. On |
---|
9 | output, elements of a above the diagonal are destroyed. d[1..n] returns the eigenvalues of a. |
---|
10 | v[1..n][1..n] is a matrix whose columns contain, on output, the normalized eigenvectors of |
---|
11 | a. nrot returns the number of Jacobi rotations that were required. |
---|
12 | */ |
---|
13 | void jacobi(float *a, int n, float d[], float *v, int *nrot) |
---|
14 | { |
---|
15 | int j,iq,ip,i; |
---|
16 | float tresh,theta,tau,t,sm,s,h,g,c,*b,*z; |
---|
17 | |
---|
18 | b = (float *) malloc(n * sizeof(float)); |
---|
19 | z = (float *) malloc(n * sizeof(float)); |
---|
20 | |
---|
21 | /* Fork a team of threads giving them their own copies of variables */ |
---|
22 | #pragma omp parallel private(nthreads, tid, i, j, ip, iq) shared(tresh, theta, tau, sm, s, h, g, c, b, z, a, n, d, v, nrot) |
---|
23 | { |
---|
24 | |
---|
25 | /* Obtain thread number */ |
---|
26 | tid = omp_get_thread_num(); |
---|
27 | printf("Hello World from thread = %d\n", tid); |
---|
28 | |
---|
29 | /* Only master thread does this */ |
---|
30 | if (tid == 0) |
---|
31 | { |
---|
32 | nthreads = omp_get_num_threads(); |
---|
33 | printf("Number of threads = %d\n", nthreads); |
---|
34 | } |
---|
35 | |
---|
36 | |
---|
37 | #pragma omp parallel for |
---|
38 | for (ip=0;ip<n;ip++) { // Initialize to the identity matrix. |
---|
39 | for (iq=0;iq<n;iq++) v[ip*n + iq]=0.0; |
---|
40 | v[ip*n + ip]=1.0; |
---|
41 | } |
---|
42 | #pragma omp parallel for |
---|
43 | for (ip=0;ip<n;ip++) { // Initialize b and d to the diagonal of a. |
---|
44 | b[ip]=d[ip]=a[ip*n + ip]; |
---|
45 | z[ip]=0.0; // This vector will accumulate terms of the form t*a[pq] as in equation (11.1.14). |
---|
46 | } |
---|
47 | } |
---|
48 | *nrot=0; |
---|
49 | for (i=1;i<=50;i++) { |
---|
50 | sm=0.0; |
---|
51 | #pragma omp parallel for reduction(+:sm) |
---|
52 | for (ip=0;ip<n-1;ip++) { // Sum off-diagonal elements. |
---|
53 | for (iq=ip+1;iq<n;iq++) |
---|
54 | sm += fabs(a[ip*n + iq]); |
---|
55 | } |
---|
56 | if (sm == 0.0) { // The normal return, which relies on quadratic convergence to machine underflow. |
---|
57 | free(z); |
---|
58 | free(b); |
---|
59 | return; |
---|
60 | } |
---|
61 | if (i < 4) |
---|
62 | tresh=0.2*sm/(n*n); // ...on the first three sweeps. |
---|
63 | else |
---|
64 | tresh=0.0; // ...thereafter. |
---|
65 | for (ip=0;ip<n-1;ip++) { |
---|
66 | for (iq=ip+1;iq<n;iq++) { |
---|
67 | g=100.0*fabs(a[ip*n + iq]); |
---|
68 | if (i > 4 && (float)(fabs(d[ip])+g) == (float)fabs(d[ip]) // After four sweeps, skip the rotation if the off-diagonal element is small. |
---|
69 | && (float)(fabs(d[iq])+g) == (float)fabs(d[iq])) |
---|
70 | a[ip*n + iq]=0.0; |
---|
71 | else if (fabs(a[ip*n + iq]) > tresh) { |
---|
72 | h=d[iq]-d[ip]; |
---|
73 | if ((float)(fabs(h)+g) == (float)fabs(h)) |
---|
74 | t=(a[ip*n + iq])/h; // t = 1/(2*theta) |
---|
75 | else { |
---|
76 | theta=0.5*h/(a[ip*n + iq]); // Equation (11.1.10). |
---|
77 | t=1.0/(fabs(theta)+sqrt(1.0+theta*theta)); |
---|
78 | if (theta < 0.0) t = -t; |
---|
79 | } |
---|
80 | |
---|
81 | c=1.0/sqrt(1+t*t); |
---|
82 | s=t*c; |
---|
83 | tau=s/(1.0+c); |
---|
84 | h=t*a[ip*n + iq]; |
---|
85 | z[ip] -= h; |
---|
86 | z[iq] += h; |
---|
87 | d[ip] -= h; |
---|
88 | d[iq] += h; |
---|
89 | a[ip*n + iq]=0.0; |
---|
90 | |
---|
91 | #pragma omp parallel for |
---|
92 | for (j=0;j<=ip-1;j++) { // Case of rotations 1 <= j < p. |
---|
93 | ROTATE(a,j,ip,j,iq) |
---|
94 | } |
---|
95 | #pragma omp parallel for |
---|
96 | for (j=ip+1;j<=iq-1;j++) { // Case of rotations p < j < q. |
---|
97 | ROTATE(a,ip,j,j,iq) |
---|
98 | } |
---|
99 | #pragma omp parallel for |
---|
100 | for (j=iq+1;j<n;j++) { // Case of rotations q < j <= n. |
---|
101 | ROTATE(a,ip,j,iq,j) |
---|
102 | } |
---|
103 | #pragma omp parallel for |
---|
104 | for (j=0;j<n;j++) { |
---|
105 | ROTATE(v,j,ip,j,iq) |
---|
106 | } |
---|
107 | |
---|
108 | ++(*nrot); |
---|
109 | } |
---|
110 | } |
---|
111 | } |
---|
112 | |
---|
113 | for (ip=0;ip<n;ip++) { |
---|
114 | b[ip] += z[ip]; |
---|
115 | d[ip]=b[ip]; // Updte d with the sum of t*a[pq], |
---|
116 | z[ip]=0.0; // and reinitialize z. |
---|
117 | } |
---|
118 | } |
---|
119 | |
---|
120 | printf("Too many iterations in routine jacobi\n"); |
---|
121 | } |
---|
122 | |
---|
123 | /* |
---|
124 | Given the eigenvalues d[1..n] and eigenvectors v[1..n][1..n] as output from jacobi |
---|
125 | (x11.1) or tqli (x11.3), this routine sorts the eigenvalues into descending order, and rearranges |
---|
126 | the columns of v correspondingly. The method is straight insertion. |
---|
127 | */ |
---|
128 | void eigsrt(float d[], float *v, int n) |
---|
129 | { |
---|
130 | int k,j,i; |
---|
131 | float p; |
---|
132 | for (i=0;i<n-1;i++) { |
---|
133 | p=d[k=i]; |
---|
134 | for (j=i+1;j<n;j++) |
---|
135 | if (d[j] >= p) p=d[k=j]; |
---|
136 | if (k != i) { |
---|
137 | d[k]=d[i]; |
---|
138 | d[i]=p; |
---|
139 | for (j=0;j<n;j++) { |
---|
140 | p=v[j*n + i]; |
---|
141 | v[j*n + i]=v[j*n + k]; |
---|
142 | v[j*n + k]=p; |
---|
143 | } |
---|
144 | } |
---|
145 | } |
---|
146 | } |
---|