c++ - 使用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;
}

最佳答案

DSBTRD不计算特征值。它将基质还原为对角线形式；您正在拉出所得三对角矩阵的主要对角线，并假装这些是特征值，但事实并非如此。

您需要在生成的对角线形式上调用DSTERF（或其他一些例程之一）以获取特征值。

有关更多详细信息，请参见LAPACK User's Guide。

关于c++ - 使用LAPACK计算无限带状矩阵的特征值，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/13217177/