[90] | 1 | #include <stdio.h> |
---|
[89] | 2 | #include <math.h> |
---|
[90] | 3 | #include <malloc.h> |
---|
[91] | 4 | #include <omp.h> |
---|
| 5 | |
---|
[89] | 6 | #define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\ |
---|
| 7 | a[k][l]=h+s*(g-h*tau); |
---|
| 8 | |
---|
| 9 | /* |
---|
| 10 | Computes all eigenvalues and eigenvectors of a real symmetric matrix a[1..n][1..n]. On |
---|
| 11 | output, elements of a above the diagonal are destroyed. d[1..n] returns the eigenvalues of a. |
---|
| 12 | v[1..n][1..n] is a matrix whose columns contain, on output, the normalized eigenvectors of |
---|
| 13 | a. nrot returns the number of Jacobi rotations that were required. |
---|
| 14 | */ |
---|
| 15 | void jacobi(float **a, int n, float d[], float **v, int *nrot) |
---|
| 16 | { |
---|
| 17 | int j,iq,ip,i; |
---|
[90] | 18 | float tresh,theta,tau,sm,s,h,g,c,*b,*z; |
---|
| 19 | float t; |
---|
| 20 | b = (float *) malloc(n * sizeof(float)); |
---|
| 21 | z = (float *) malloc(n * sizeof(float)); |
---|
[89] | 22 | |
---|
[91] | 23 | int nthreads, tid; |
---|
| 24 | |
---|
| 25 | /* Fork a team of threads giving them their own copies of variables */ |
---|
[99] | 26 | #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) |
---|
[91] | 27 | { |
---|
| 28 | |
---|
| 29 | /* Obtain thread number */ |
---|
| 30 | tid = omp_get_thread_num(); |
---|
| 31 | printf("Hello World from thread = %d\n", tid); |
---|
| 32 | |
---|
| 33 | /* Only master thread does this */ |
---|
| 34 | if (tid == 0) |
---|
| 35 | { |
---|
| 36 | nthreads = omp_get_num_threads(); |
---|
| 37 | printf("Number of threads = %d\n", nthreads); |
---|
| 38 | } |
---|
| 39 | |
---|
[99] | 40 | |
---|
| 41 | #pragma omp parallel for |
---|
[90] | 42 | for (ip=0;ip<n;ip++) { // Initialize to the identity matrix. |
---|
| 43 | for (iq=0;iq<n;iq++) |
---|
| 44 | v[ip][iq]=0.0; |
---|
| 45 | v[ip][ip]=1.0; |
---|
[89] | 46 | } |
---|
[99] | 47 | |
---|
| 48 | #pragma omp parallel for |
---|
[90] | 49 | for (ip=0;ip<n;ip++) { // Initialize b and d to the diagonal of a. |
---|
| 50 | b[ip]=d[ip]=a[ip][ip]; |
---|
| 51 | z[ip]=0.0; // This vector will accumulate terms of the form t*a[pq] as in equation (11.1.14). |
---|
[89] | 52 | } |
---|
| 53 | |
---|
[99] | 54 | } /* All threads join master thread and disband */ |
---|
| 55 | |
---|
[89] | 56 | *nrot=0; |
---|
| 57 | for (i=1;i<=50;i++) { |
---|
[99] | 58 | |
---|
| 59 | sm = 0.0; |
---|
| 60 | |
---|
[90] | 61 | for (ip=0;ip<n-1;ip++) { // Sum off-diagonal elements. |
---|
[99] | 62 | #pragma omp parallel for reduction(+:sm) |
---|
| 63 | for (iq=ip+1;iq<n;iq++) { |
---|
| 64 | sm += fabs(a[ip][iq]); |
---|
| 65 | printf("Hello World from thread = %d\n", omp_get_thread_num()); |
---|
| 66 | } |
---|
[89] | 67 | } |
---|
[90] | 68 | if (sm == 0.0) { // The normal return, which relies on quadratic convergence to machine underflow. |
---|
| 69 | free(z); |
---|
| 70 | free(b); |
---|
[89] | 71 | return; |
---|
| 72 | } |
---|
| 73 | if (i < 4) |
---|
[90] | 74 | tresh=0.2*sm/(n*n); // ...on the first three sweeps. |
---|
[89] | 75 | else |
---|
[90] | 76 | tresh=0.0; // ...thereafter. |
---|
[99] | 77 | |
---|
[89] | 78 | for (ip=0;ip<n-1;ip++) { |
---|
| 79 | for (iq=ip+1;iq<n;iq++) { |
---|
| 80 | g=100.0*fabs(a[ip][iq]); |
---|
[90] | 81 | 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. |
---|
[89] | 82 | && (float)(fabs(d[iq])+g) == (float)fabs(d[iq])) |
---|
| 83 | a[ip][iq]=0.0; |
---|
| 84 | else if (fabs(a[ip][iq]) > tresh) { |
---|
| 85 | h=d[iq]-d[ip]; |
---|
| 86 | if ((float)(fabs(h)+g) == (float)fabs(h)) |
---|
[90] | 87 | t=(a[ip][iq])/h; // t = 1/(2*theta) |
---|
[89] | 88 | else { |
---|
[90] | 89 | theta=0.5*h/(a[ip][iq]); |
---|
[89] | 90 | t=1.0/(fabs(theta)+sqrt(1.0+theta*theta)); |
---|
| 91 | if (theta < 0.0) t = -t; |
---|
| 92 | } |
---|
| 93 | |
---|
| 94 | c=1.0/sqrt(1+t*t); |
---|
| 95 | s=t*c; |
---|
| 96 | tau=s/(1.0+c); |
---|
| 97 | h=t*a[ip][iq]; |
---|
| 98 | z[ip] -= h; |
---|
| 99 | z[iq] += h; |
---|
| 100 | d[ip] -= h; |
---|
| 101 | d[iq] += h; |
---|
| 102 | a[ip][iq]=0.0; |
---|
[99] | 103 | |
---|
| 104 | #pragma omp parallel for |
---|
[92] | 105 | for (j=0;j<=ip-1;j++) { // Case of rotations 1 <= j < p. |
---|
[89] | 106 | ROTATE(a,j,ip,j,iq) |
---|
[90] | 107 | } |
---|
[99] | 108 | #pragma omp parallel for |
---|
[92] | 109 | for (j=ip+1;j<=iq-1;j++) { // Case of rotations p < j < q. |
---|
[89] | 110 | ROTATE(a,ip,j,j,iq) |
---|
| 111 | } |
---|
[99] | 112 | #pragma omp parallel for |
---|
[90] | 113 | for (j=iq+1;j<n;j++) { // Case of rotations q < j <= n. |
---|
[89] | 114 | ROTATE(a,ip,j,iq,j) |
---|
| 115 | } |
---|
[99] | 116 | #pragma omp parallel for |
---|
[89] | 117 | for (j=0;j<n;j++) { |
---|
| 118 | ROTATE(v,j,ip,j,iq) |
---|
| 119 | } |
---|
| 120 | |
---|
| 121 | ++(*nrot); |
---|
| 122 | } |
---|
| 123 | } |
---|
| 124 | } |
---|
| 125 | |
---|
[99] | 126 | #pragma omp parallel for |
---|
[89] | 127 | for (ip=0;ip<n;ip++) { |
---|
| 128 | b[ip] += z[ip]; |
---|
[90] | 129 | d[ip]=b[ip]; // Update d with the sum of t*a[pq], |
---|
| 130 | z[ip]=0.0; // and reinitialize z. |
---|
[89] | 131 | } |
---|
| 132 | } |
---|
| 133 | |
---|
| 134 | printf("Too many iterations in routine jacobi\n"); |
---|
[90] | 135 | } |
---|
| 136 | |
---|
[91] | 137 | int main(int argc, char * argv[]) |
---|
[90] | 138 | { |
---|
| 139 | int n, i, j; |
---|
| 140 | float *d; |
---|
| 141 | float **a; |
---|
| 142 | float **v; |
---|
| 143 | int nrot = 0; |
---|
[91] | 144 | |
---|
| 145 | char * filename = argv[1]; |
---|
| 146 | FILE * f = fopen(filename, "r"); |
---|
[90] | 147 | fscanf(f,"%d", &n); |
---|
| 148 | |
---|
| 149 | a = (float **)malloc(n * sizeof(float*)); |
---|
| 150 | for(i = 0; i < n; i++) |
---|
| 151 | a[i] = (float *)malloc(n * sizeof(float)); |
---|
| 152 | |
---|
| 153 | v = (float **)malloc(n * sizeof(float*)); |
---|
| 154 | for(i = 0; i < n; i++) |
---|
| 155 | v[i] = (float *)malloc(n * sizeof(float)); |
---|
| 156 | |
---|
| 157 | d = (float *)malloc(n * sizeof(float)); |
---|
| 158 | |
---|
| 159 | for(i = 0; i < n; i++) |
---|
| 160 | for(j = 0; j < n; j++) |
---|
| 161 | fscanf(f,"%f", &a[i][j]); |
---|
| 162 | for(i = 0; i < n; i++) |
---|
| 163 | { |
---|
| 164 | for(j = 0; j < n; j++) |
---|
| 165 | printf("%f ", a[i][j]); |
---|
| 166 | printf("\n"); |
---|
| 167 | } |
---|
| 168 | |
---|
| 169 | jacobi(a, n, d, v, &nrot); |
---|
| 170 | |
---|
| 171 | printf("v:\n"); |
---|
| 172 | for(i = 0; i < n; i++) |
---|
| 173 | { |
---|
| 174 | for(j = 0; j < n; j++) |
---|
| 175 | printf("%f ", v[i][j]); |
---|
| 176 | printf("\n"); |
---|
| 177 | } |
---|
| 178 | |
---|
| 179 | printf("d:\n"); |
---|
| 180 | for(i = 0; i < n; i++) |
---|
| 181 | { |
---|
| 182 | printf("%f ", d[i]); |
---|
| 183 | } |
---|
[91] | 184 | printf("\n"); |
---|
[90] | 185 | |
---|
| 186 | return 0; |
---|
| 187 | } |
---|
| 188 | |
---|