Having a few issues with Lapacke DSYGV

[ctcsback]

I am trying to use dsygv to solve a system A*x = (lambda)*B*x, where A and B are real symmetric matrices.

However, unlike the original lapack documentation, the lapacke documentation is quite lacking. I wrote some test code which attempts to solve for only eigenvectors. I used some example code which can be found here: http://www.google.com/url?sa=t&source=w ... qQ&cad=rja

I modified that and used the Lapacke interface to run their example and I'm double checking the numbers using Matlab, but I haven't had any luck. I've posted my code below. Maybe I'm missing something, but I can't find it at the moment.

Compiling is just: >>gcc main.c lapacke.a -llapack -lblas -I../include

and I'm not having any issues with that. Any help appreciated! :)

Code: Select all


#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <lapacke.h> // <--- header file

void dsygv(double *A, double *B, int n, double *E)
{
  char jobz, uplo;
  int itype, lda, ldb, info, i;
  double *work;

  itype = 1; /* This sets the type of generalized eigenvalue problem
      that we are solving.  We have the possible values
          1: Ax = lBx
          2: ABx = lx
          3: BAx = lx */

  jobz = 'N'; /* V/N indicates that eigenvectors should/should not
       be calculated. */

  uplo = 'U'; /* U/L indicated that the upper/lower triangle of the
       symmetric matrix is stored. */

  lda = n; // The leading dimension of the matrix A
  ldb = n; // The leading dimension of the matrix B

  work = malloc(sizeof(double)*(3*n-1)); // Allocate work vector

  info = LAPACKE_dsygv(LAPACK_ROW_MAJOR, itype, jobz, uplo, n, A, lda, B, ldb, work);

  assert(info==0);
  free(work);
}


void main(int nargs, char *args[])
{
  double *A, *B, *Evals;
  int i, j, k;

  int n=4;

  FILE *fp = fopen("debug.txt", "w");

  // Generate the matrices

  A = calloc(n*n,sizeof(double));
  for (i=0; i<n; i++)
  {
    for (j=i; j<n; j++)
    {
      A[i*n+j] = pow(i+1, j+1) + pow(j+1, i+1);
      A[j*n+i] = pow(i+1, j+1) + pow(j+1, i+1);
    }
  }

  B = calloc(n*n,sizeof(double));
  for (i=0; i<n; i++)
  {
    for (j=i; j<n; j++)
    {
      B[i*n+j] = pow(i+1, 2*(j+1)) + pow(j+1, 2*(i+1));
      B[j*n+i] = pow(i+1, 2*(j+1)) + pow(j+1, 2*(i+1));
    }
  }

  debug_print("Amat", A, n, n, 0);
  debug_print("Bmat", B, n, n, 0);

  Evals = calloc(n, sizeof(double));

  dsygv(A, B, n, Evals);

  // Output the results

  for (i=0; i<n; i++) {
    printf("Eigenvalue  at %d is %f\n", i, A[i]);
  }

  debug_print("Amat", A, n, n, 0);
  debug_print("Bmat", B, n, n, 0);
  free(Evals);
}



[admin]

indeed, you got a little mix up.

First thing to remember is that LAPACKE works the same way than LAPACK, as it is just a wrapper.
Second LAPACKE takes care of the workspace. Here w was the eigenvalue you were looking for!

Julie
Here is your program working:

Code: Select all

#include <stdio.h>
#include <math.h>
#include <assert.h>
#include "lapacke.h" // <--- header file

void debug_print(char *name, double *A, int m, int n, int k)
{
   int i,j;
   printf("====== %s =====\n",name);
   for (i=0; i< m; i++)
      {
      for (j=0; j< n; j++)
        printf ("%e ",A[i*n+j]);
      printf("\n");
      }
}

void dsygv(double *A, double *B, int n, double *E)
{
  char jobz, uplo;
  int itype, lda, ldb, info, i;

  itype = 1; /* This sets the type of generalized eigenvalue problem
      that we are solving.  We have the possible values
          1: Ax = lBx
          2: ABx = lx
          3: BAx = lx */

  jobz = 'N'; /* V/N indicates that eigenvectors should/should not
       be calculated. */

  uplo = 'U'; /* U/L indicated that the upper/lower triangle of the
       symmetric matrix is stored. */

  lda = n; // The leading dimension of the matrix A
  ldb = n; // The leading dimension of the matrix B


  info = LAPACKE_dsygv(LAPACK_ROW_MAJOR, itype, jobz, uplo, n, A, lda, B, ldb, E);

  assert(info==0);
}


void main(int nargs, char *args[])
{
  double *A, *B, *Evals;
  int i, j, k;

  int n=4;

  FILE *fp = fopen("debug.txt", "w");

  // Generate the matrices

  A = calloc(n*n,sizeof(double));
  for (i=0; i<n; i++)
  {
    for (j=i; j<n; j++)
    {
      A[i*n+j] = pow(i+1, j+1) + pow(j+1, i+1);
      A[j*n+i] = pow(i+1, j+1) + pow(j+1, i+1);
    }
  }

  B = calloc(n*n,sizeof(double));
  for (i=0; i<n; i++)
  {
    for (j=i; j<n; j++)
    {
      B[i*n+j] = pow(i+1, 2*(j+1)) + pow(j+1, 2*(i+1));
      B[j*n+i] = pow(i+1, 2*(j+1)) + pow(j+1, 2*(i+1));
    }
  }

  debug_print("Amat", A, n, n, 0);
  debug_print("Bmat", B, n, n, 0);

  Evals = calloc(n, sizeof(double));

  dsygv(A, B, n, Evals);

  // Output the results

  printf("=== Eigenvalue ====\n");
  for (i=0; i<n; i++) {
    printf("%f\n", Evals[i]);
  }

//  debug_print("Amat", A, n, n, 0);
//  debug_print("Bmat", B, n, n, 0);
  free(Evals);
}


[ctcsback]

Seems like a simple mistake on my part. Thanks for the help.