help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Help-gsl] Newbie eigenvalue issues


From: Dimitrova, Maria
Subject: [Help-gsl] Newbie eigenvalue issues
Date: Fri, 23 Sep 2016 11:23:57 +0000

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


address@hidden




/************************************************************************************************
 *                                                                              
                *
 *   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;
}





reply via email to

[Prev in Thread] Current Thread [Next in Thread]