Вычисление собственных значений бесконечной полосчатой ​​матрицы с использованием LAPACK

Я новичок в C ++ и пытаюсь выяснить, как использовать LAPACK, чтобы найти собственные значения бесконечной полосчатой ​​матрицы (проблема ангармонического осциллятора). Я знаю, что правильно вычисляю матрицу, так как проверил значения, и все они совпадают. Однако я не уверен, правильно ли я передаю значения подпрограмме или я что-то перепутал, поскольку возвращаемые собственные значения не те, которые я ожидаю. Я использую подпрограмму dsbtrd для вычисления этого. Вот руководство для этого: http://www.netlib.org/lapack/explore-html/d0/d62/dsbtrd_8f.html

Любые идеи о том, где я могу пойти не так?

#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;
}