CODEWITHSUNDEEP
c++lapackeigenvalue
opes
opes

Reputation: 1828

Calculating eigenvalues of an infinite banded matrix using LAPACK

I'm new to C++ and am trying to figure out how to use LAPACK to find the eigenvalues of an infinite banded matrix (anharmonic oscillator problem). I know that I'm calculating the matrix correctly as I've checked the values and they all match up. However, I'm not sure if I'm passing the values to the subroutine correctly or if I've got something mixed up as the eigenvalues that are being returned are not what I am expecting. I'm using the dsbtrd subroutine to compute this. Here's the manual for that: http://www.netlib.org/lapack/explore-html/d0/d62/dsbtrd_8f.html

Any ideas on where I might be going wrong?

#include <iostream>
#include <algorithm>
#include <string>
#include <math.h>
using namespace std;

//SUBROUTINE DSBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO )

extern "C" {
  void dsbtrd_(const char *vect, const char *uplo, int *n, int *kd, double *ab, int *ldab, double d[], double e[], double *q, int *ldq, double work[], int *info);
}

#define MAX 14

int main(){
    // Values needed for dsbtrd
    const char *vect = "V";
    const char *uplo = "U";
    int n;
    int ldab = MAX;
    int ldq = MAX;
    int info;
    double ab[MAX][ldab];
    double d[MAX];
    double e[MAX];
    double q[MAX][ldq];
    double work[MAX];

    //other values needed
    int i,j,delta;
    double eps;
    double a[MAX][MAX];
    double g[MAX][MAX];

    //Read in value of eps and n
    cout <<"Enter epsilon: \n";
    cin >> eps;
    cout << "Epsilon = " << eps << "\n";

    cout <<"Enter n: \n";
    cin >> n;
    cout << "n = " << n << "\n";

    if(n >= MAX){
         cerr << "n is great than max \n";
     }

     //Build matrix g
     for(j = 0; j < n; j++){
         for(i = 0; i < n; i++){
             int m = min(i,j);

             if(i == j){
                 g[i][j] = 1.5*(pow(m,2) + m +.5);
             }else if( i == j + 2 || i == j - 2){
                 g[i][j] = (m + 1.5)*sqrt((m+1)*(m+2));
             }else if(i == j + 4 || i == j -4){
                 g[i][j] = .25*sqrt((m+1)*(m+2)*(m+3)*(m+4));
             }else{
                 g[i][j] = 0;
             }
         }
    }

    //Build the starting matrix a
    //row i, column j

    for(j = 0; j < n; j++){
        for(i = 0; i < n; i++){
            if(i == j){
                delta = 1;
            }else{
                delta = 0;
            }

            a[i][j] = (i + .5)*delta + eps*g[i][j];
        }
    }

    //Build the matrix ab
    // ab(kd+1+i-j,j) = a(i,j) for max(1,j-kd)<=i<=j
    int kd = n - 1;

    for(j = 1; j <= n; j++){
        for(i = max(1,j-kd); i <= j; i++){
            ab[j-1][kd+i-j] = a[j-1][i-1];
        }
    }

    //Solve for eigenvalues
    dsbtrd_(vect, uplo, &n, &kd, &ab[0][0], &ldab, d, e, &q[0][0], &ldq, work, &info);

    //Check for success
    if(info == 0)
    {
        //Write answer
        for(i = 0; i < n; i++){
            cout << "Eigenvalue " << i << ": " << d[i] << "\n";
        }
    }
    else
    {
        //Write error
        cerr << "dsbtrd returned error " << info << "\n";
    }
    return info;
}

Upvotes: 0

Views: 920

Answers (1)

Stephen Canon
Stephen Canon

Reputation: 106267

DSBTRD doesn't calculate eigenvalues. It reduces the matrix to tridiagonal form; you're pulling out the main diagonal of the resulting tridiagonal matrix and pretending that those are the eigenvalues, but they aren't.

You need to call DSTERF (or one of a few other routines) on the resulting tridiagonal form to get the eigenvalues.

For more details, consult the LAPACK User's Guide.

Upvotes: 2

Related Questions