DelFEMで行列の固有値を求める（３） [一般化固有値問題の場合] 【2012-10-12修正】
前回の記事では[A]{x}=λ{x}の固有値を求めましたが、実はこの形式の固有値問題をそのまま数値解析に適用できたことはあまりありません。たいてい[A]{x}=λ[B]{x}の形式の一般化固有値問題を解くことになります。
単純な方法としては[B]の逆行列を求めてそれを両辺に掛けて計算すればいいのですが、下記に分かりやすい記事がありましたのでこれを元に[A]{x}=λ[B]{x}を計算してみました。
http://kentiku-kouzou.jp/taishin-kouuchi.html
一般固有値問題、(1)式の場合には、(30)式の代わりに、
を用いる。詳しく述べると、
つまり、べき乗数法の更新ベクトル残差ベクトルに質量行列を掛けてやるだけのようです。
逆反復法でも同じようです。以下一般化固有値問題用に修正したソースコード。
/* 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="massMat">質量行列バッファ(行列[i,j]を[i + j * RowSize]にマップしたバッファ)</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, double[[]] massMat, 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, massMat); { // 更新を正規化 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, massMat); // 更新ベクトルを解ベクトルにコピー 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, massMat); // レリー積の計算 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, double[[]] massMat) { DelFEM4NetMatVec.CVector_Blk_Ptr resPtr = ls.GetResidualPtr(0); uint nblk = resPtr.NBlk(); double[] saveResVec = new double[nblk]; for (uint iblk = 0; iblk < nblk; iblk++) { saveResVec[iblk] = resPtr.GetValue(iblk, 0); } System.Diagnostics.Debug.Assert(nblk * nblk == massMat.Length); for (int ino = 0; ino < nblk; ino++) { double sum = 0; for (int jno = 0; jno < nblk; jno++) { // M[ino, jno] sum += massMat[ino + jno * nblk] * saveResVec[jno]; } resPtr.SetValue((uint)ino, 0, sum); } } } }
ソース中のproductMassMatToUpdateVec productResidualMatToUpdateVecというのが質量行列を更新ベクトル残差ベクトルに掛けている処理です。
合っているのか確かめてみます。例題が下記にありました。
http://sake03.shinshu-u.ac.jp/dynamics/note/Dyna02-03exp.html
[K] = { 610e+5, -260e+5}
{-260e+5, 260e+5}
[M] = {0.20e+5, 0.0}
{ 0.0, 0.15e+5}
のとき、([K] - ω^2[M]){x} = {0}の固有値を求める問題です。
ω^2 = 753と4030ということなので、753を取得できればOKです。
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 = 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;
}
}
// 作成したインデックスをリニアシステムの対角行列？に登録
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");
for (uint iblk = 0; iblk < nblk; iblk++)
{
DoubleArrayIndexer ptr = femMat.GetPtrValDia(iblk);
System.Diagnostics.Debug.Assert(ptr.Count == 1);
ptr[0] = stiffnessMat[iblk, iblk];
Console.WriteLine(" ( " + iblk + " ) = " + ptr[0]);
uint npsup = 0;
ConstUIntArrayIndexer ptrInd = femMat.GetPtrIndPSuP(iblk, out npsup);
DoubleArrayIndexer ptrVal = femMat.GetPtrValPSuP(iblk, out npsup);
System.Diagnostics.Debug.Assert(ptrInd.Count == ptrVal.Count);
for (int i = 0; i < ptrVal.Count; i++)
{
ptrVal[i] = stiffnessMat[iblk, ptrInd[i]];
Console.WriteLine(" ( " + iblk + ", " + ptrInd[i] + " ) = " + ptrVal[i]);
}
}
// 更新ベクトル(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);
// 行列の固有値を取得する
// 最大値が取得できるようです。→嘘かもしれない。最小値の求まる行列もある。
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>(); //何をセットすればいいのか?
/*
uint ntmp = ls.GetTmpVectorArySize();
for (uint itmp = 0; itmp < ntmp; itmp++)
{
aIdVec.Add(itmp);
}
*/
Console.WriteLine("massMat");
for (int ino = 0; ino < massMat.GetLength(0); ino++)
{
for (int jno = 0; jno < massMat.GetLength(1); jno++)
{
Console.WriteLine(" M({0},{1}) = {2}", ino, jno, massMat[ino, jno]);
}
}
double[] massMatBuffer = matrix_ToBuffer(massMat);
double eigenVal = CEigenLanczos_CSharp.MinimumEigenValueVector_InvPower_General(
ls,
prec,
massMatBuffer,
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));
}
ここでソースコード中にあるmatrix_ToBufferというのは２次元配列を一次元に直す関数です。
private static double[] matrix_ToBuffer(double[,] mat)
{
// rowから先に埋めていく
for (int j = 0; j < mat.GetLength(1); j++) //col
{
for (int i = 0; i < mat.GetLength(0); i++) // row
{
mat_[j * mat.GetLength(0) + i] = mat[i, j];
}
}
return mat_;
}
実行結果：
＜削除＞
あれ、誤差かな？と思ったのですが、上記引用元は四捨五入しているようです。
うまくいったようです。
更新ベクトルに掛けたとき固有値は正しいのですが、固有ベクトルが変でした。
下記が修正後の実行結果です。KrdLabさんのLisysでも計算して比較してみたので
今度はうまくいっていると思います（思いたい）。
