DelFEMで行列の固有値を求める(4)[一般化固有値問題の場合(2)]
前回の記事に誤りがありましたのでお詫びいたします。ごめんなさい。
今回は前回の記事で作成したソースコードを少しDelFEM風に修正したので貼り付けたいと思います。
前回のメソッド：
public static double MinimumEigenValueVector_InvPower_General(
DelFEM4NetLsSol.CLinearSystem ls,
DelFEM4NetLsSol.CPreconditioner pls,
double[[]] massMat,
IList aIdVec,
uint itr_invp, // 逆べき乗法の最大反復回数
uint itr_lssol, // ICCG法の最大反復回数
double conv_ratio_lssol, // ICCG法の収束基準の相対残差
out int iflag_conv) // 0:正常に終了 1:ICCGが収束しなかった
剛性行列はCLinearSystemに隠されているのに対して、質量行列は生のdouble[[]]で渡していて違和感があります。今回はこれを"matvec"モジュールのCMatDia_BlkCrsに変えてみました。
/* DelFEM (Finite Element Analysis) Copyright (C) 2009 Nobuyuki Umetani n.umetani@gmail.com This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ // // eigen_lanczos.cppをC#に移植 // created by rujimiya (original codes created by Nobuyuki Umetani // using System; using System.Runtime.InteropServices; using System.Collections.Generic; namespace DelFEM4NetLsSol { public class CEigenLanczos_CSharp { /// <summary> /// 逆べき乗法による最小固有値の取得(一般化固有値問題) /// </summary> /// <param name="ls">リニアシステムオブジェクト</param> /// <param name="pls">プリコンディショナ―オブジェクト</param> /// <param name="diaMassMat">質量行列</param> /// <param name="aIdVec">無視ベクトル</param> /// <param name="itr_invp">逆べき乗法の最大反復回数</param> /// <param name="itr_lssol">ICCG法の最大反復回数</param> /// <param name="conv_ratio_lssol">ICCG法の収束基準の相対残差</param> /// <param name="iflag_conv">0:正常に終了 1:ICCGが収束しなかった</param> /// <returns></returns> public static double MinimumEigenValueVector_InvPower_General( DelFEM4NetLsSol.CLinearSystem ls, DelFEM4NetLsSol.CPreconditioner pls, DelFEM4NetMatVec.CMatDia_BlkCrs diaMassMat, IList<uint> aIdVec, uint itr_invp, // 逆べき乗法の最大反復回数 uint itr_lssol, // ICCG法の最大反復回数 double conv_ratio_lssol, // ICCG法の収束基準の相対残差 out int iflag_conv) // 0:正常に終了 1:ICCGが収束しなかった { uint max_itr_invp = itr_invp; if (ls.GetTmpVectorArySize() < 3) { ls.ReSizeTmpVecSolver(3); } int ires = -1; int iupd = -2; int itmp = 2; // 無視ベクトルのグラムシュミットの直交化 for (int iivec = 0; iivec < aIdVec.Count; iivec++) { int ivec0 = (int)aIdVec[iivec]; for (int jivec = 0; jivec < iivec; jivec++) { int jvec0 = (int)aIdVec[jivec]; double dot = ls.DOT(ivec0, jvec0); ls.AXPY(-dot, jvec0, ivec0); } double sq_normx = ls.DOT(ivec0, ivec0); double normx = Math.Sqrt(sq_normx); ls.SCAL(1.0 / normx, ivec0); } // 更新から無視ベクトル成分を引く for (int iivec = 0; iivec < aIdVec.Count; iivec++) { int ivec0 = (int)aIdVec[iivec]; double dot0 = ls.DOT(ivec0, iupd); ls.AXPY(-dot0, ivec0, iupd); } //BUGFIX 下記は間違いなので削除 //// 更新ベクトルに質量行列を左から掛ける //productMassMatToUpdateVec(ls, diaMassMat); { // 更新を正規化 double sq_normx = ls.DOT(iupd, iupd); double normx = Math.Sqrt(sq_normx); ls.SCAL(1.0 / normx, iupd); } double min_eigen = double.MaxValue; for (uint work_itr_invp = 0; work_itr_invp < max_itr_invp; work_itr_invp++) { // 更新ベクトルを残差ベクトルにコピー ls.COPY(iupd, ires); // 残差ベクトルに質量ベクトルを左から掛ける productMassMatToResidualVec(ls, diaMassMat); // 更新ベクトルを解ベクトルにコピー ls.COPY(iupd, itmp); { // 行列の逆を求める double conv_ratio = conv_ratio_lssol; uint num_iter = itr_lssol; DelFEM4NetLsSol.CLinearSystemPreconditioner lsp = new DelFEM4NetLsSol.CLinearSystemPreconditioner(ls, pls); bool res = DelFEM4NetLsSol.CSolverLsIter.Solve_PCG(ref conv_ratio, ref num_iter, lsp); //Console.WriteLine("PCG iter : " + num_iter + " " + conv_ratio); if (!res) { iflag_conv = 1; // ICCGが収束しなかったフラグ return 0; } } // 無視ベクトルの成分を引く for (int iivec = 0; iivec < aIdVec.Count; iivec++) { int ivec0 = (int)aIdVec[iivec]; double dot0 = ls.DOT(ivec0, iupd); ls.AXPY(-dot0, ivec0, iupd); } // BUGFIX 下記は間違いなので削除 //// 更新ベクトルに質量行列を左から掛ける //productMassMatToUpdateVec(ls, diaMassMat); // レリー積の計算 min_eigen = ls.DOT(itmp, iupd); // 更新ベクトルを正規化 double norm0 = Math.Sqrt(ls.DOT(iupd, iupd)); ls.SCAL(1.0 / norm0, iupd); } iflag_conv = 0; // 正常終了フラグ return 1.0 / min_eigen; } private static void productMassMatToResidualVec(CLinearSystem ls, DelFEM4NetMatVec.CMatDia_BlkCrs diaMassMat) { // 残差ベクトルのポインタを取得 DelFEM4NetMatVec.CVector_Blk_Ptr resPtr = ls.GetResidualPtr(0); uint nblk = resPtr.NBlk(); System.Diagnostics.Debug.Assert(nblk == diaMassMat.NBlkMatCol()); System.Diagnostics.Debug.Assert(nblk == diaMassMat.NBlkMatRow()); // 書換え用に残差ベクトルポインタを残差ベクトル（のリファレンス）にキャストする DelFEM4NetMatVec.CVector_Blk resVec = resPtr; // コピーを作成 DelFEM4NetMatVec.CVector_Blk saveResVec = new DelFEM4NetMatVec.CVector_Blk(resVec.NBlk(), (uint)resVec.Len()); saveResVec.op_Assign(resVec); // 質量行列と更新ベクトルの演算 // {y} = alpha * [A]{x} + beta * {y} diaMassMat.MatVec(1.0, saveResVec, 0.0, ref resVec); } } }
質量行列を残差ベクトルに掛ける処理 productMassMatToResidualVecがすっきりしました。
CMatDia_BlkCrsなどのDelFEM4NetMatVec名前空間にあるモジュールは、本家DelFEMのMatVec名前空間の線形方程式を計算するときの行列やベクトルの操作モジュール群のラッパーです。 DelFEM4Net version1.0.0.5でバグを修正したので使えると思います(というか、この実装をしていて不具合に気付きました、ごめんなさい）。
さて、作成したメソッドを実行してみます。呼び出し側はこんな感じ。
{
double[,] stiffnessMat = new double[2, 2]
{
{ 610.0e5, -260.0e5},
{-260.0e5, 260.0e5},
};
double[,] massMat = new double[2, 2]
{
{0.20e5, 0.00e5},
{0.00e5, 0.15e5},
};
Console.WriteLine("CEigenLanczos.MinimumEigenValueVector_InvPower_General");
// DelFEMの固有値解析ルーチンを使用
// リニア方程式ソルバ
CLinearSystem ls = new CLinearSystem();
// プリコンディショナ―
CPreconditioner_ILU prec = new CPreconditioner_ILU();
// 行列の０でない要素のある場所を指定し、インデックスを格納する
uint nblk = (uint)stiffnessMat.GetLength(0);
ls.AddLinSysSeg(nblk, 1);
// 行列の要素のインデックスを格納するインスタンスを生成
{
CIndexedArray crs = makePattern(stiffnessMat);
// 作成したインデックスをリニアシステムの対角行列？に登録
ls.AddMat_Dia(0, crs);
}
//prec.SetFillInLevel(1);
// ILU(0)のパターン初期化
prec.SetLinearSystem(ls);
// 行列の要素を格納する場所を取得
CMatDia_BlkCrs_Ptr femMat = ls.GetMatrixPtr(0);
// 行列へのマージ開始
ls.InitializeMarge();
//uint nblk = femMat.NBlkMatCol();
// 行列を要素を追加します。
System.Diagnostics.Debug.Assert(nblk == femMat.NBlkMatCol());
Console.WriteLine("srcMat merge");
CMatDia_BlkCrs femMatRef = femMat;
doMatMerge(ref femMatRef, stiffnessMat, true);
// 更新ベクトル(1, 1,..., 1)
CVector_Blk_Ptr updatePtr = ls.GetUpdatePtr(0);
System.Diagnostics.Debug.Assert(nblk == updatePtr.NBlk());
System.Diagnostics.Debug.Assert(1 == updatePtr.Len());
for (uint iblk = 0; iblk < updatePtr.NBlk(); iblk++)
{
double val = 1.0;
updatePtr.SetValue(iblk, 0, val);
//System.Console.WriteLine("update( " + iblk + " ) = " + updatePtr.GetValue(iblk, 0));
}
// マージ終了
double res = ls.FinalizeMarge();
//Console.WriteLine("Residual : " + res );
// プリコンディショナ―にリニアシステムの(行列の？)値を設定してILU分解を行う
bool ret = prec.SetValue(ls);
//Console.WriteLine("prec.SetValue(ls) ret = " + ret);
CMatDia_BlkCrs diaMassMat = new CMatDia_BlkCrs(nblk, 1);
{
CIndexedArray crs = makePattern(massMat);
diaMassMat.DeletePattern();
diaMassMat.AddPattern(crs);
//すべての要素を非０要素として埋めるならFillPatternが使える
//diaMassMat.FillPattern();
}
// 要素の値を０にする
// ここで値を格納する領域が確保されるので必須
diaMassMat.SetZero();
Console.WriteLine("massMat merge");
doMatMerge(ref diaMassMat, massMat);
// 行列の固有値を取得する
// 最大値が取得できるようです。→嘘かもしれない。最小値の求まる行列もある。
double conv_res_lssol = 1.0e-6;// ICCG法の収束基準の相対残差
uint itr_invp = 2000;// 逆べき乗法の最大反復回数
uint itr_lssol = 2000;// ICCG法の最大反復回数
int iflag_conv = 0; // 0:正常に終了 1:ICCGが収束しなかった
IList<uint> aIdVec = new List<uint>(); //何をセットすればいいのか?
double eigenVal = CEigenLanczos_CSharp.MinimumEigenValueVector_InvPower_General(
ls,
prec,
diaMassMat,
aIdVec,
itr_invp,
itr_lssol,
conv_res_lssol,
out iflag_conv);
// iflag_conv (0:正常に終了 1:ICCGが収束しなかった)
Console.WriteLine("////// Result of CEigenLanczos.MinimumEigenValueVector_InvPower //////");
Console.WriteLine("eigenVal = " + eigenVal);
Console.WriteLine("sqrt(eigenVal) = " + Math.Sqrt(eigenVal));
Console.WriteLine("iflag_conv = " + iflag_conv);
for (uint iblk = 0; iblk < updatePtr.NBlk(); iblk++)
{
System.Console.WriteLine("eigenVec( " + iblk + " ) = " + updatePtr.GetValue(iblk, 0));
}
Console.WriteLine("press any key");
Console.In.Read();
}
ここでmakePatternとdoMatMergeがでてきますが、それらは下記のとおりです。
/////////////////////////////////////////////////
/// <summary>
/// マトリクスのパターンを作成する
/// </summary>
/// <param name="mat"></param>
/// <returns></returns>
private static CIndexedArray makePattern(double[,] doubleMat)
{
uint nblk = (uint)doubleMat.GetLength(0);
CIndexedArray crs = new CIndexedArray();
crs.InitializeSize(nblk);
using (UIntVectorIndexer index = crs.index)
using (UIntVectorIndexer ary = crs.array)
{
index[0] = 0;
// 行を走査
for (int col = 0; col < nblk; col++)
{
// 列を走査
for (int row = 0; row < nblk; row++)
{
// ここではn x n の行列の要素すべてを格納できるようにします。
// 本当は非０の場合のみ追加するものだと思われます。
// 対角要素は除きます。
if (col != row)
{
ary.Add((uint)row);
}
}
index[col + 1] = (uint)ary.Count;
}
}
return crs;
}
/// <summary>
/// マトリクスをマージする
/// </summary>
/// <param name="mat"></param>
/// <param name="doubleMat"></param>
private static void doMatMerge(ref CMatDia_BlkCrs mat, double[,] doubleMat, bool printFlg= false)
{
uint nblk = mat.NBlkMatCol();
for (uint iblk = 0; iblk < nblk; iblk++)
{
DoubleArrayIndexer ptr = mat.GetPtrValDia(iblk);
System.Diagnostics.Debug.Assert(ptr.Count == 1);
ptr[0] = doubleMat[iblk, iblk];
if (printFlg)
{
Console.WriteLine(" ( " + iblk + " ) = " + ptr[0]);
}
uint npsup = 0;
ConstUIntArrayIndexer ptrInd = mat.GetPtrIndPSuP(iblk, out npsup);
DoubleArrayIndexer ptrVal = mat.GetPtrValPSuP(iblk, out npsup);
System.Diagnostics.Debug.Assert(ptrInd.Count == ptrVal.Count);
for (int i = 0; i < ptrVal.Count; i++)
{
ptrVal[i] = doubleMat[iblk, ptrInd[i]];
if (printFlg)
{
Console.WriteLine(" ( " + iblk + ", " + ptrInd[i] + " ) = " + ptrVal[i]);
}
}
}
}
DelFEMのmatvec、LinearSystemモジュールでは、行列の要素パターン作成→要素の値をマージというのがパターンのようなので、それぞれ関数にしました。この関数はまだ修正の余地がありそうです。それと、SetZeroというメソッドは必ず呼び出した方がいいと思います。これを呼ばないと”NULL"の配列に値をマージしようとして落ちます。
実行結果：
