[Help-gsl] Real generalized symmetric-definite eigensystems

[Dimitrova, Maria]

Hello,

I have yet another question about matrix operations, this time regarding 
eigensystems. In a sample program I define two matrices and try to solve the 
eigenvalue problem. Both matrices are symmetric and all of their elements are 
positive. However, on calling the function gsl_eigen_gensymmv I get an error 
that they have to be positive definite. And well, they are. The function 
gsl_matrix_isnonneg returns 1.

Here is the output, followed by the source code. Thank you very much.

$ gcc -Wall -Wextra diagonalization.c -lm -lgsl -lgslcblas && ./a.out
=======================================================================================

Initialized matrices
Matrix H

 1.000000   2.000000   4.000000
 2.000000   3.000000   5.000000
 4.000000   5.000000   6.000000

Matrix S

 7.000000   8.000000   10.000000
 8.000000   9.000000   11.000000
 10.000000   11.000000   12.000000

=======================================================================================

Check for positive definiteness
"H: "
stat=1
"S: "
stat=1
=======================================================================================

Generalized eigenvalue problem
gsl: cholesky.c:157: ERROR: matrix must be positive definite
Default GSL error handler invoked.
Aborted (core dumped)


/************************************************************************************************
 *   COMPILED WITH THE LINE:                                                    
                *
 *   gcc -Wall -Wextra diagonalization.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 3
#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 elem=0; // used in the initialization of the arrays
    int stat=0;
    int i,j;
    gsl_matrix *C, *H, *S;
    gsl_vector *E;
    gsl_eigen_gensymmv_workspace *wrkEig; // for generalized eigenvalue problem
    S = gsl_matrix_calloc(MatrixOrder,MatrixOrder);
    H = gsl_matrix_calloc(MatrixOrder,MatrixOrder);
    C = gsl_matrix_calloc(MatrixOrder,MatrixOrder);
    E = gsl_vector_calloc(MatrixOrder);
    wrkEig = gsl_eigen_gensymmv_alloc(MatrixOrder); // for NxN matrix wrk is 
O(4N)
/*****************************************************************************************/
    // Initialization of the matrices
    for (i=0;i<MatrixOrder;i++)
    {
        for (j=0;j<=i;j++)
        {
            elem++;
            gsl_matrix_set(H,i,j,elem);
            gsl_matrix_set(H,j,i,elem);
        }
    } // H[i][j] initialized
    DIV;
    MSG(Initialized matrices);
    MSG(Matrix H);
    printArr(H,MatrixOrder,MatrixOrder);
    for (i=0;i<MatrixOrder;i++)
    {
        for (j=0;j<=i;j++)
        {
            elem++;
            gsl_matrix_set(S,i,j,elem);
            gsl_matrix_set(S,j,i,elem);
        }
    } // S[i][j] initialized
    MSG(Matrix S);
    printArr(S,MatrixOrder,MatrixOrder);
/*****************************************************************************************/

    DIV;
    MSG(Check for positive definiteness);
    stat=gsl_matrix_isnonneg(H); //return 1 if all the elements of the matrix m 
are non-negative, and 0 otherwise
    DEBUG("H: ",stat,d);
    DEBUG("S: ",stat,d);

    DIV;
    MSG(Generalized eigenvalue problem);
    gsl_eigen_gensymmv(H,S,E,C,wrkEig);

/*****************************************************************************************/
    gsl_matrix_free(H);
    gsl_matrix_free(S);
    gsl_vector_free(E);
    gsl_matrix_free(C);
    gsl_eigen_gensymmv_free(wrkEig);
    return 0;
}



Best regards,

Maria Dimitrova

address@hidden