Hello, I am getting started with the GSL libraries and tested how finding the eigenvalues and eigenvectors of matrices works. I spotted some weird results, so I must have made a mistake. The output from my program looks like this:
Matrix S 2.000000 2.000000 3.000000 0.000000 View matrix S 2.000000 2.000000 3.000000 0.000000 Eigenvalues of S: -2.162278 4.162278 However, online calculators do not agree. The eigenvalues are given as 3.65 and -1.64. Here are the references to them http://www.wolframalpha.com/input/?i=eigenvalues+%7B%7B2,2%7D,%7B3,0%7D%7D www.mathportal.org/calculators/matrices-calculators/matrix-calculator.php?formId=1&val1=2%3A2%3Anull%3Anull%3Anull%3Anull%3B3%3A0%3Anull%3Anull%3Anull%3Anull%3Bnull%3Anull%3Anull%3Anull%3Anull%3Anull%3Bnull%3Anull%3Anull%3Anull%3Anull%3Anull%3Bnull%3Anull%3Anull%3Anull%3Anull%3Anull%3Bnull%3Anull%3Anull%3Anull%3Anull%3Anull&val2=2&val3=2&rb1=evec My code is below. Thank you for the help. Best regards, Maria [email protected] /************************************************************************************************ * * * COMPILED WITH THE LINE: * * gcc -Wall -Wextra eigenval-vect.c -lm -lgsl -lgslcblas && ./a.out * * * * * ************************************************************************************************/ #include <stdio.h> #include <stdlib.h> #include <math.h> #include <gsl/gsl_math.h> #include <gsl/gsl_blas.h> #include <gsl/gsl_vector.h> #include <gsl/gsl_matrix.h> #include <gsl/gsl_eigen.h> #define MatrixOrder 2 #define MSG(msg) printf( "\n" #msg "\n") #define DIV printf("\n=======================================================================================\n\n") #define DEBUG(msg, var, fmt) printf( #msg "\n" #var "=%" #fmt "\n", var) #define PRINTF(var, fmt) printf("\n**DEBUG: " #var "=%" #fmt "\n", var) void printArr(gsl_matrix *array, int nrows, int ncols) { int i,j; printf("\n"); for (i=0;i<nrows;i++) { printf("\n"); for (j=0;j<ncols;j++) { printf("% f ",gsl_matrix_get(array,i,j)); } } printf("\n\n"); } int main() { int i=0; double data[] = { 2.0, 2.0, 3.0, 0.0 }; gsl_matrix *S, *C; gsl_vector *E; gsl_eigen_symmv_workspace *wrkDiag; // for diagonalization S = gsl_matrix_calloc(MatrixOrder,MatrixOrder); C = gsl_matrix_calloc(MatrixOrder,MatrixOrder); E = gsl_vector_calloc(MatrixOrder); wrkDiag = gsl_eigen_symmv_alloc(MatrixOrder); // for NxN matrix wrk is O(4N) /*****************************************************************************************/ gsl_matrix_view m = gsl_matrix_view_array (data, 2, 2); MSG(Matrix S); printArr(&m.matrix,2,2); MSG(View matrix S); gsl_matrix_fprintf(stdout,&m.matrix,"%f"); /*****************************************************************************************/ DIV; gsl_eigen_symmv (&m.matrix,E,C,wrkDiag); gsl_eigen_symmv_sort(E,C,GSL_EIGEN_SORT_VAL_ASC); MSG(Eigenvalues of S:); for (i=0;i<2;i++) { printf("% f\t",gsl_vector_get(E,i)); } printf("\n\n"); MSG(Eigenvectors of S:); printArr(C,2,2); DIV; /*****************************************************************************************/ gsl_matrix_free(S); gsl_vector_free(E); gsl_matrix_free(C); gsl_eigen_symmv_free(wrkDiag); return 0; }
