Gmm++ の最初のステップ

どのように逆行列を計算しますか？

Gmm++ では全ての種類の行列の逆行列を計算することはできません．今のところ，逆行列を計算する唯一の方法は，密なLU分解(つまり，密行列の場合のみ)を使用することです．例えば

gmm::dense_matrix<double> M(3, 3), M2(3,3), M3(3,3);
gmm::copy(gmm::identity_matrix(), M);  // M = Id.
gmm::scale(M, 2.0);                    // M = 2 * Id.
M(1,2) = 1.0;

gmm::copy(M, M2);

gmm::lu_inverse(M);

gmm::mult(M, M2, M3);

std::cout << M << " times " << M2 << " is equal to " << M3 << endl;

詳細については，密なLU分解に対応するセクションを参照してください．型 gmm::dense_matrix<double> は， gmm::row_matrix< std::vector<double> > または gmm::col_matrix< std::vector<double> > で置き換えることができます．

線形システムをどのように解きますか？

線形システムを解く方法は複数あります．密行列の場合は，LU分解を使用するのが最善です．例えば

gmm::dense_matrix<double> M(3, 3);
gmm::clear(M);                  // M = 0.
M(0,0) = M(1,1) = M(2,2) = 2.0; // M = 2 * Id.
M(1,2) = 1.0;

std::vector<double> X(3), B(3), Bagain(3);
B[0] = 1.0; B[1] = 2.0; B[2] = 3.0;  // B = [1 2 3]

gmm::lu_solve(M, X, B);

gmm::mult(M, X, Bagain);

std::cout << M << " times " << gmm::vref(X) << " is equal to " << gmm::vref(Bagain) << endl;

今，偏微分方程式の離散化から来る疎なシステムを持っているなら，前処理行列の有無にかかわらず，さまざまな反復ソルバーがあります．これは，前処理GMRESによる例です．

int nbdof = 1000; // number of degrees of freedom.
gmm::row_matrix< gmm::rsvector<double> > M(nbdof, nbdof); // a sparse matrix
std::vector<double> X(nbdof), B(nbdof); // Unknown and left hand side.

... here the assembly of the pde discretization stiffness matrix ...
... and left hand side ...


// computation of a preconditioner (ILUT)
gmm::ilut_precond< gmm::row_matrix< gmm::rsvector<double> > > P(M, 10, 1e-4);

gmm::iteration iter(1E-8);  // defines an iteration object, with a max residu of 1E-8

gmm::gmres(M, X, B, P, 50, iter);  // execute the GMRES algorithm

std::cout << "The result " << gmm::vref(X) << endl;

どのようにベクトルを行列に変換し，それを変形することができますか？

Gmm++ では，ベクトルは行列とはみなされません．ベクトルを（1×n）行行列または（n×1）列行列として計算で使用する必要がある場合は，次を使用する必要があります．

gmm::row_vector(V) // gives a reference on V considered as
                   // a (1 by n) row matrix
gmm::col_vector(V) // gives a reference on V considered as
                   // a (n by 1) col matrix

たとえば，次のようにしてベクトルを密行列に変換できます．

std::vector<double> V(50);

// ... computation of V

gmm::dense_matrix<double> M(1, gmm::vect_size(V));
gmm::copy(gmm::row_vector(V), M);

次のようにして，行列 M の形状を変更することもできます．行列

gmm::reshape(M, 10, 5);

行列のサイズを変更するより良い方法は何ですか？

行列の大きさは，参照でない場合は次を使用して変更することができます．

gmm::resize(M, m, n);

この関数は元の行列とサイズ変更された行列の交点を尊重し，新しい要素はゼロに設定されます．重要なことは， std::vector のサイズ変更メソッドに基づいていることです．したがって，新しい行列のサイズが元のサイズよりも小さい場合，メモリフリーは行われません．

コンポーネントの古い値を保持する必要がない場合，または余分なメモリを解放したい場合は，次のように std::swap を使って行列のサイズを変更できます．

MATRIX_TYPE M(m1, n1);

... your code

{ MATRIX_TYPE(m2, n2) M2; std::swap(M, M2); } // resize matrix M.

もちろん，これはベクトルに対しても機能します．