פרויקט תוכנה, שאלות ותשובות!

void dominant_eig(sparse_matrix_arr* matrix, float precision, float* eig_vec, float *eig_val,int nnz)
{
elem *X; //a random vector that is used in the power method
elem *Y; //a random vector that is used in the power method
elem temp_eig_val; //will save the current evaluation of the eigen value
int i;
int N=matrix->n; //the matrix size;
int m=0;
elem eig_num; //temp_eig_val = A(A^m*X)*Y / (A^m*X)*Y; eig_num
elem eig_denom; //and eig_denom are the upper and the bottom
//parts of the fraction corresponding (num=upper,denom=bottom
elem *temp_vec_1;
elem *temp_vec_2;
elem dif; // dif= || A(A^m*X)-lambda1*(A^m*X)||



X=(elem *)calloc(N,sizeof(elem),');
for (i=0;i {
X[i]=(rand() % 10,');
}




for (m=1;m>0;m++)
{

/*
Y=A^m*X
lambda1= A(A^m*X)Y / (A^m*X)Y //eig_val
V1=A^m*X //eig_vec

or in other words:
lambda1=(AY)*Y / Y*Y
V1=Y

condition:

|| A(A^m*X)-lambda1*(A^m*X)|| <=epsilon
or in other words:

||A*Y - lamba1*Y||<=epsilon



we need to compute: A*X, A^2*X, A^3*X,... and we compute it as A*X,
A(A*X), A(A(A*X)),...

*/


Y=mat_vector_mult(matrix,Y,'); //Y=A*Y = A(A(A(A...(A*X))))...) = A^m*X
temp_vec_1=mat_vector_mult(matrix,Y,'); //temp_vec1=AY = A^(m+1) * X
eig_num=vec_mult_vec(temp_vec_1,Y,N,'); //eig_num=(AY)*Y
eig_denom=vec_mult_vec(Y,Y,N,'); //eig_denom=Y*Y
temp_eig_val = eig_num / eig_denom;

//temp_vec_2=lambda1*Y
for (i=0; i {
temp_vec_2[i]=temp_eig_val*Y[i];
}
dif=dif_precision(temp_vec_1,temp_vec_2,N,'); // dif=||A*Y - lamba1*Y||

if (dif <= precision)
{
for (i=0;i {
eig_vec[i]=Y[i];
}
*eig_val=temp_eig_val;


free(temp_vec_1,');
free(temp_vec_2,');
free(X,');
free(Y,');
return;
}

normalize_vec(Y,N,');
}
}

0 0 1
0 0 1
1 1 0

N=3, nnz=4;
A->values[0]=1;
A->values[1]=1;
A->values[2]=1;
A->values[3]=1;

A->rowind[0]=2;
A->rowind[1]=2;
A->rowind[2]=0;
A->rowind[3]=1;

A->colptr[0]=0;
A->colptr[1]=1;
A->colptr[2]=2;
A->colptr[3]=4;






0 1 0
1 0 1
0 1 0

N=3,nnz=4;
A->values[0]=1;
A->values[1]=1;
A->values[2]=1;
A->values[3]=1;

A->rowind[0]=1;
A->rowind[1]=0;
A->rowind[2]=2;
A->rowind[3]=1;

A->colptr[0]=0;
A->colptr[1]=1;
A->colptr[2]=3;
A->colptr[3]=4;




0 1 1
1 0 0
1 0 0

N=3,nnz=4;
A->values[0]=1;
A->values[1]=1;
A->values[2]=1;
A->values[3]=1;

A->rowind[0]=1;
A->rowind[1]=2;
A->rowind[2]=0;
A->rowind[3]=0;

A->colptr[0]=0;
A->colptr[1]=2;
A->colptr[2]=3;
A->colptr[3]=4;