#include #include #include #include using namespace std; float** makematrix(int n, int m); void freematrix(int n, float **a); float powermethod(float** A, int n, float tolerance, long* idum, int k) { //our current and previous guess for the eigenvalue float guess1 = 1, guess2 = 10; //our guess for the eigenvector float *eigenvec1 = (float*)calloc(n, sizeof(float)); //vectors to work with float *eigenvec2 = (float*)calloc(n, sizeof(float)); float *vec3 = (float*)calloc(n, sizeof(float)); //initialize our first guess for(int i = 0; i < n; i++) eigenvec1[i] = (float)rand()/RAND_MAX; while(guess1-guess2 <= -tolerance || guess1-guess2 >= tolerance) { float c = 0, d = 0, normconst = 0,test = 0; // set the value of c and zero out vec3 for(i = 0; i < n; i++) { c += eigenvec1[i]; vec3[i] = 0; } c = c/n; //set the value of d for(i = 0; i < n/2; i++) d += eigenvec1[i]; for(i = n/2; i < n; i++) d -= eigenvec1[i]; d = d/n; //set the value of eigenvec2 for(i = 0; i < n/2; i++) eigenvec2[i] = eigenvec1[i]-c-d; for(i = n/2; i < n; i++) eigenvec2[i] = eigenvec1[i]-c+d; //normalize eigenvec2 with L2 norm for(i = 0; i < n; i++) normconst += eigenvec2[i]*eigenvec2[i]; normconst = sqrt(normconst); for(i = 0; i < n; i++) eigenvec2[i] = eigenvec2[i]/normconst; //compute A(A(eigenvec2)), place results in eigenvec1. for(i = 0; i < n; i++) for(int j = 0; j < k; j++) vec3[i] += eigenvec2[(int)A[i][j]]; for(i = 0; i < n; i++) eigenvec1[i] = 0; for(i = 0; i < n; i++) for(int j = 0; j < k; j++) eigenvec1[i] += vec3[(int)A[i][j]]; guess2 = guess1; guess1 = 0; //compute eigenvalue guess that corresponds to our //eigenvector guess for(i = 0; i < n; i++) guess1 += eigenvec2[i]*eigenvec1[i]; guess1 = sqrt(guess1); } free(eigenvec1); free(eigenvec2); free(vec3); return guess1; }