0

私は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;
}
4

1 に答える 1

2

DSBTRD は固有値を計算しません。行列を三重対角形に縮小します。結果の三重対角行列の主対角線を引き出し、それらが固有値であるふりをしていますが、そうではありません。

固有値を取得するには、結果の三重対角形式で DSTERF (またはいくつかの他のルーチンの 1 つ) を呼び出す必要があります。

詳細については、LAPACK ユーザー ガイドを参照してください。

于 2012-11-16T13:27:39.943 に答える