インテルMKLライブラリの一般化固有値問題が解けません。

#インテルの数値演算ライブラリを用いて、物理の固有振動を求めようとしています。
固有振動は一般化固有値問題に相当するため、大規模行列の固有値と固有ベクトルと同義です。
求めようとしている、固有方程式は　Ax = λBx でxがλが固有値、xが固有ベクトルです。
具体的には、インテルMKLライブラリのLAPACKルーチンの一つである、dsygvx()を使用して一般化固有値問題において指定した次数までの固有値と固有ベクトルを計算しようとしています。今回は全体の次数が100ですがそのうち30個を計算したいです。
現段階では、全ての固有値と固有ベクトルを求めるdspgv()を実装完了しています。

　コードは以下の通りです。

C++
1コード
2    //****入力パラメータ****
3	int itype = 1; //汎用対象固有値問題
4	char jobz = 'V'; //固有値と固有ベクトルを計算
5	char uplo = 'U';　　//下三角行列
6	int n = 100;　　//行列A,Bの次元数
7	int lda = n, ldb = n; //1次元
8	int lwork = 8 * n; 
9	double *a = new double[lda*n];
10	double *b = new double[ldb*n];
11	double *ap = new double[n*(n + 1) / 2]; //n*(n+1)/2
12	double *bp = new double[n*(n + 1) / 2];
13	double *work = new double[lwork];
14	char range = 'I';   //求めたい固有値の範囲
15	int il = 1, iu = 30; //固有値の個数( il < iu <= n) 
16	int ldz = n;
17	double abstol;  //求めた固有値の絶対許容誤差
18	int iwork = 5 * n;
19	
20	int count = 0;
21
22	for (int i = 0; i < n; i++) {　　//ap,bpに関しても同様に代入する(次元が違うだけ)が今回は省略。
23		for (int j = 0; j < i + 1; j++) {
24			a[count] = A.get(j, i);  //A行列を転置して、下三角形行列として代入
25			b[count] = B.get(j, i);  //B行列を転置して、下三角形行列として代入
26			count++;
27		}
28	}
29
30	//****出力パラメータ****
31	int info;
32	int m;
33	double *w = new double[n];
34	double *z = new double[ldz*n];
35	int *ifail = new int[n];
36
37	//全ての固有値の導出
38	//dspgv(&itype, &jobz, &uplo, &n, ap, bp, w, z, &ldz, work, &info);  //これは実行可能
39
40	//固有値の数指定
41	dsygvx(&itype, &jobz, &range, &uplo, &n, a, &lda, b, &ldb, &vl, &vu, &il, &iu, &abstol, &m, w, z, &ldz, work, &lwork, &iwork, ifail, &info);  //実行できない、、
42
43	if (info != 0) {　//固有値計算が終了していない
44		cout << "Error" << endl;
45		cout << "info= " << info << endl;
46		getchar();
47	}
48	else {
49		cout << "info= " << info << endl;
50		fp_MECH = fopen("data/6-29/eigenvalue_out_30.csv", "w");
51		for (int i = 0; i < n; i++) {
52			fprintf(fp_MECH, "%lf\n", w[i]);
53		}
54		fclose(fp_MECH);
55

以上のプログラムにおいて、30次まで指定した固有値を求めようとした際に、固有値計算が収束せずに計算が出来ません。
どなたかこのプログラムを正常に実行できるよう教えて欲しいです。

jbpb0

2021/07/01 23:25 編集

abstolを指定しないといけないのでは？参考 http://www.rcs.arch.t.u-tokyo.ac.jp/kusuhara/tips/memorandum/memo4.html#sec5 https://www.ktech.biz/doc/latest/ja/c/D4b1.html#a18ec2b64bd2db484b9c206c55fb59a95

Kam-GitHub

2021/07/02 01:14

jbpb0様　ご回答ありがとうございます。先ほど、abstolをいくつか設定して見たのですが、info = 101となり、固有値計算が収束しませんでした。 info = 101ということは、「info-n= 101-100=1次小行列が正定値でない. Bの分解が完了できず, 固有値・固有ベクトルは計算されなかった」とのことだと思うのですが、何か解決策はありますでしょうか。

jbpb0

2021/07/02 06:19

https://www.ktech.biz/doc/latest/ja/c/D4b1.html に記載の関数はどれも、そこに書かれているように、Bが正定値行列であることが使える条件にありますが、それは大丈夫でしょうか？参考 https://www.iwanttobeacat.com/entry/2019/12/30/122655 https://www.or.mist.i.u-tokyo.ac.jp/kanno/lecture/math_design/psd170420.pdf https://jp.mathworks.com/help/matlab/math/determine-whether-matrix-is-positive-definite.html

jbpb0

2021/07/03 11:58 編集

https://jp.mathworks.com/help/matlab/math/determine-whether-matrix-is-positive-definite.html の「方法 1: コレスキー分解の試行」に、「分解が失敗した場合、行列は対称正定値ではありません。」と書かれてます > Bの分解が完了できずもしかしたら、Bは正定値行列ではないのでは？もしそうなら、 https://www.ktech.biz/doc/latest/ja/c/D4b1.html#a0da976196662cd459f6ae53e6eb1b8bc に書かれてるように、「dspgv」も使えないはずなので、エラーが出てなくて計算結果が得られているとしても、その計算結果は間違っているのかも

Kam-GitHub

2021/07/03 13:39

ご返信頂きありがとうございます。遅くなってしまいすみません。教えていただいたURLを参考に、正定行列か判定するために試行錯誤していました。今まで、dspgv()で問題なく固有値計算が出来ているという認識のもとでしたが、そもそも間違っている可能性があるので、再度計算し直してみたいと思います。試しに、n=4として4×4の簡単な正定行列(下三角行列)として、dspgv()とdsygvx()どちらも計算してみたいと思います。

jbpb0

2021/07/03 15:50

> 正定行列か判定オープンソースのGNU Octaveでも、MATLABと同様にchol()で確認ができるようです参考 https://octave.org/doc/v4.0.0/Matrix-Factorizations.html