【GSLで数値計算7】
エルミート行列の固有値の固有ベクトル
文責：遠藤 理平 （2018年7月 5日）
カテゴリ：計算物理学(165)
本稿では、GSL（GNU科学技術計算ライブラリ）を用いてエルミート行列の固有値と固有ベクトルの計算を行うサンプルプログラムを示します。→ GNU科学技術計算ライブラリのインストール（Windows, VisualStudio2017）はこちら
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 | //////////////////////////////////////////////////////////////// エルミート行列の固有値の固有ベクトル//////////////////////////////////////////////////////////////#include <iostream>#include <iomanip>#include <string>#include <gsl/gsl_math.h>#include <gsl/gsl_blas.h>#include <gsl/gsl_vector.h>#include <gsl/gsl_matrix.h>#include <gsl/gsl_eigen.h>#include <gsl/gsl_sort_vector.h>#include <gsl/gsl_complex.h>#include <gsl/gsl_complex_math.h>//行列の次元const int DIM = 2;int main(void) { std::cout << std::setprecision(5); //GSL複素数ベクトルのメモリ確保（固有ベクトル用） gsl_vector_complex *eig_vec = gsl_vector_complex_alloc(DIM); //GSL複素数ベクトルのメモリ確保（固有値ベクトル用） gsl_vector *eig = gsl_vector_alloc(DIM); //GSL複素数行列のメモリ確保（エルミート行列用） gsl_matrix_complex *matrix = gsl_matrix_complex_alloc(DIM, DIM); //GSL複素数行列のメモリ確保（固有ベクトルのまとめ行列用） gsl_matrix_complex *eig_vecs = gsl_matrix_complex_alloc(DIM, DIM); //エルミート行列の設定 gsl_matrix_complex *frank_a = gsl_matrix_complex_alloc(DIM, DIM); gsl_complex A[DIM][DIM]; //GSL_SET_COMPLEX(&aa, 1.0, 0.0); // A[0][0] = 1 + 0i A[0][0] = { 1.0, 0.0 }; gsl_matrix_complex_set(frank_a, 0, 0, A[0][0]); // A[0][1] = 1 - 1i A[0][1] = { 1.0, -1.0 }; gsl_matrix_complex_set(frank_a, 0, 1, A[0][1]); // A[1][0] = 1 + 1i A[1][0] = { 1.0, 1.0 }; gsl_matrix_complex_set(frank_a, 1, 0, A[1][0]); // A[1][1] = 2 + 0i A[1][1] = { 2.0, 0.0 }; gsl_matrix_complex_set(frank_a, 1, 1, A[1][1]); // 行列の表示 std::cout << "元の行列を表示" << std::endl; for (int i = 0; i < DIM; i++) { for (int j = 0; j < DIM; j++) { gsl_complex aa = gsl_matrix_complex_get(frank_a, i, j); std::cout << GSL_REAL(aa) << " + " << GSL_IMAG(aa) << "i "; } std::cout << std::endl; } std::cout << std::endl; ///////////////////////////////////////// // 固有値のみ計算 ///////////////////////////////////////// // 作業領域の確保 gsl_eigen_herm_workspace *workspace = gsl_eigen_herm_alloc(DIM); // 行列のコピー gsl_matrix_complex_memcpy(matrix, frank_a); // 固有値の計算 gsl_eigen_herm(matrix, eig, workspace); // 固有値の並び替え gsl_sort_vector(eig); std::cout << "-----固有値のみを計算-----" << std::endl; for (int i = 0; i < DIM; i++) { std::cout << i << "番目の固有値:" << gsl_vector_get(eig, i) << std::endl; } std::cout << std::endl; //領域の開放 gsl_eigen_herm_free(workspace); ///////////////////////////////////////// // 固有値と固有ベクトルの計算 ///////////////////////////////////////// // 作業領域の確保 gsl_eigen_hermv_workspace *workspace_v = gsl_eigen_hermv_alloc(DIM); // 行列のコピー gsl_matrix_complex_memcpy(matrix, frank_a); // 固有値・固有ベクトルの計算 gsl_eigen_hermv(matrix, eig, eig_vecs, workspace_v); // 置換インデックスの保存 gsl_permutation *perm = gsl_permutation_alloc(DIM); gsl_sort_vector_index(perm, eig); // 固有値の書き出し std::cout << "-----固有値と固有ベクトルを計算-----" << std::endl; for (int i = 0; i < DIM; i++) { std::cout << i << "番目の固有値: " << gsl_vector_get(eig, gsl_permutation_get(perm, i)); std::cout << " 固有ベクトル: "; for (int j = 0; j < DIM; j++) { gsl_complex z_temp = gsl_matrix_complex_get(eig_vecs, i, j); std::cout << GSL_REAL(z_temp); if (GSL_IMAG(z_temp) >= 0) std::cout << "+"; std::cout << GSL_IMAG(z_temp) << "i "; } std::cout << std::endl; } std::cout << std::endl; //領域の開放 gsl_permutation_free(perm); ///////////////////////////////////////// // 計算チェック ///////////////////////////////////////// // AX - EX == 0 ? gsl_matrix_complex *tmp_matrix = gsl_matrix_complex_alloc(DIM, DIM); // 単位行列に設定 gsl_matrix_complex_set_identity(tmp_matrix); //固有値による対角ベクトルを準備（matrix） for (int i = 0; i < DIM; i++) { //実部と虚部 gsl_complex aa = { gsl_vector_get(eig, i) , 0 }; for (int j = 0; j < DIM; j++) { gsl_matrix_complex_set(matrix, j, i, aa); } } //EXを計算 gsl_matrix_complex_mul_elements(matrix, eig_vecs); gsl_complex al = { 1.0, 0.0 }; gsl_complex beta = { -1.0, 0.0 }; //AX-EXを計算 gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, al, frank_a, eig_vecs, beta, matrix); std::cout << "検算 AX-XE = " << std::endl; for (int i = 0; i < DIM; i++) { for (int j = 0; j < DIM; j++) { gsl_complex z_temp = gsl_matrix_complex_get(matrix, i, j); std::cout << GSL_REAL(z_temp); if (GSL_IMAG(z_temp) >= 0) std::cout << "+"; std::cout << GSL_IMAG(z_temp) << "i "; } std::cout << std::endl; } std::cout << std::endl; //領域の開放 gsl_matrix_complex_free(eig_vecs); gsl_matrix_complex_free(tmp_matrix); gsl_eigen_hermv_free(workspace_v); gsl_matrix_complex_free(frank_a); gsl_matrix_complex_free(matrix); gsl_vector_complex_free(eig_vec); gsl_vector_free(eig); std::string s; std::cin >> s; return 0;} |
実行結果
いよいよ次はシュレディンガー方程式の固有値問題に取り組みます。
