DelFEMで複素エルミート行列の固有値を求める
このブログの最初の記事を改めてみると、DelFEMで実対称固有値を求めたときのソースコードの公開でした。
初心に戻って、今回は複素エルミート行列の固有値を求めてみました。[A]{x} = λ{x}のシンプルな固有値問題です。
内積の定義を複素数用に変える以外はほぼ実対称行列の固有値と同じでいけそうです。そこで「DelFEMで行列の固有値を求める（２）」で作成したルーチンを下記のように改造してみました。
逆反復法による複素エルミート行列の固有値解法ルーチン
/* 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> /// 逆べき乗法によるエルミート行列の最小固有値の取得 /// 固有値に負の固有値が含まれる場合は正しく計算されません。残差[A]{x} - λ{x}をチェックしてください。 /// </summary> /// <param name="ls">リニアシステムオブジェクト</param> /// <param name="pls">プリコンディショナ―オブジェクト</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 DelFEM4NetCom.Complex MinimumComplexEigenValueVector_InvPower( DelFEM4NetFem.Ls.CZLinearSystem ls, DelFEM4NetFem.Ls.CZPreconditioner_ILU prec, 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]; DelFEM4NetCom.Complex dot = ls.INPROCT(ivec0, jvec0); ls.AXPY(-dot, jvec0, ivec0); } double normx = Math.Sqrt(DelFEM4NetCom.Complex.Norm(ls.INPROCT(ivec0, ivec0))); ls.SCAL(1.0 / normx, ivec0); } // 更新から無視ベクトル成分を引く for (int iivec = 0; iivec < aIdVec.Count; iivec++) { int ivec0 = (int)aIdVec[iivec]; DelFEM4NetCom.Complex dot0 = ls.INPROCT(ivec0, iupd); ls.AXPY(-dot0, ivec0, iupd); } { // 更新を正規化 //DelFEM4NetCom.Complex normx = complex_Sqrt(ls.INPROCT(iupd, iupd)); double normx = Math.Sqrt(DelFEM4NetCom.Complex.Norm(ls.INPROCT(iupd, iupd))); ls.SCAL(1.0 / normx, iupd); } DelFEM4NetCom.Complex min_eigen = new DelFEM4NetCom.Complex(double.MaxValue, 0); for (uint work_itr_invp = 0; work_itr_invp < max_itr_invp; work_itr_invp++) { ls.COPY(iupd, ires); ls.COPY(iupd, itmp); { // 行列の逆を求める double conv_ratio = conv_ratio_lssol; uint num_iter = itr_lssol; //bool res = DelFEM4NetFem.Ls.CZSolverLsIter.Solve_PCG(ref conv_ratio, ref num_iter, ls, prec); //bool res = DelFEM4NetFem.Ls.CZSolverLsIter.Solve_CG(ref conv_ratio, ref num_iter, ls); bool res = DelFEM4NetFem.Ls.CZSolverLsIter.Solve_PCOCG(ref conv_ratio, ref num_iter, ls, prec); //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]; DelFEM4NetCom.Complex dot0 = ls.INPROCT(ivec0, iupd); ls.AXPY(-dot0, ivec0, iupd); } // レリー積の計算 //min_eigen = ls.INPROCT(itmp, iupd); min_eigen = ls.INPROCT(iupd, itmp); //DelFEM4NetCom.Complex normtmp = complex_Sqrt(ls.INPROCT(itmp, itmp)); double normtmp = Math.Sqrt(DelFEM4NetCom.Complex.Norm(ls.INPROCT(itmp, itmp))); min_eigen /= normtmp; // 更新ベクトルを正規化 //DelFEM4NetCom.Complex norm0 = complex_Sqrt(ls.INPROCT(iupd, iupd)); double norm0 = Math.Sqrt(DelFEM4NetCom.Complex.Norm(ls.INPROCT(iupd, iupd))); ls.SCAL(1.0 / norm0, iupd); } iflag_conv = 0; // 正常終了フラグ return 1.0 / min_eigen; } } }
また、途中で現れる連立方程式の解法は、PCGからPCOCGに変更しています。
最後、固有値を計算するときレイリ商を求めていますが、実対称行列では解のノルムが１に規格化されていましたが、PCOCGでは規格化されていないようなので、規格化の処理を追加しています。
では、実際に簡単な例で動作を確認してみます。
[A] = [ 3.0, -j1.0 ]
[ j1.0, 2.0]
このエルミート行列の 固有値は正の実数なので解けるでしょう。
ルーチンの呼び出し側処理はこちら。
double m1 = 3.0;
double m2 = 2.0;
double k = 1.0;
DelFEM4NetCom.Complex imagOne = new Complex(0.0, 1.0);
DelFEM4NetCom.Complex[,] srcMat = new DelFEM4NetCom.Complex[2, 2]
{
{ m1, -imagOne * k},
{imagOne * k, m2}
};
// DelFEMの固有値解析ルーチンを使用
// リニア方程式ソルバ
CZLinearSystem ls = new CZLinearSystem();
// プリコンディショナ―
CZPreconditioner_ILU prec = new CZPreconditioner_ILU();
// 行列の０でない要素のある場所を指定し、インデックスを格納する
uint nblk = (uint)srcMat.GetLength(0);
CFieldWorld world = null;
uint fieldValId = 0;
int[[]] tmpBuffer = null;
bool[,] matPattern = new bool[nblk, nblk];
for (int i = 0; i < matPattern.GetLength(0); i++)
{
for (int j = 0; j < matPattern.GetLength(1); j++)
{
matPattern[i, j] = true;
}
}
HPlaneWGSimulatorXDelFEM.DelFEMLsUtil.SetupLinearSystem(
matPattern,
out world,
out fieldValId,
out ls,
out prec,
out tmpBuffer);
// 行列の要素を格納する場所を取得
CZMatDia_BlkCrs_Ptr femMat = ls.GetMatrixPtr(fieldValId, ELSEG_TYPE.CORNER, world);
// 行列へのマージ開始
ls.InitializeMarge();
//uint nblk = femMat.NBlkMatCol();
// 行列を要素を追加します。
System.Diagnostics.Debug.Assert(nblk == femMat.NBlkMatCol());
Console.WriteLine("srcMat merge");
for (uint iblk = 0; iblk < nblk; iblk++)
{
ComplexArrayIndexer ptr = femMat.GetPtrValDia(iblk);
System.Diagnostics.Debug.Assert(ptr.Count == 1);
ptr[0] = srcMat[iblk, iblk];
Console.WriteLine(" ( " + iblk + " ) = " + ptr[0].Real + " + " + ptr[0].Imag + " i");
uint npsup = 0;
ConstUIntArrayIndexer ptrInd = femMat.GetPtrIndPSuP(iblk, out npsup);
ComplexArrayIndexer ptrVal = femMat.GetPtrValPSuP(iblk, out npsup);
System.Diagnostics.Debug.Assert(ptrInd.Count == ptrVal.Count);
for (int i = 0; i < ptrVal.Count; i++)
{
ptrVal[i] = srcMat[iblk, ptrInd[i]];
Console.WriteLine(" ( " + iblk + ", " + ptrInd[i] + " ) = " + ptrVal[i].Real + " + " + ptrVal[i].Imag + " i");
}
}
// 更新ベクトル(1, 1,..., 1)
CZVector_Blk_Ptr updatePtr = ls.GetUpdatePtr(fieldValId, ELSEG_TYPE.CORNER, world);
System.Diagnostics.Debug.Assert(nblk == updatePtr.BlkVecLen());
System.Diagnostics.Debug.Assert(1 == updatePtr.BlkLen());
for (uint iblk = 0; iblk < updatePtr.BlkVecLen(); iblk++)
{
DelFEM4NetCom.Complex val = new DelFEM4NetCom.Complex(1.0, 0.0);
updatePtr.SetValue(iblk, 0, val);
//System.Console.WriteLine("update( " + iblk + " ) = " + updatePtr.GetValue(iblk, 0));
}
// マージ終了
double res = ls.FinalizeMarge();
//Console.WriteLine("Residual : " + res );
// プリコンディショナ―にリニアシステムの(行列の？)値を設定してILU分解を行う
prec.SetFillInLevel(int.MaxValue);
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が収束しなかった
DelFEM4NetCom.Complex eigenVal = CEigenLanczos_CSharp.MinimumComplexEigenValueVector_InvPower(
ls,
prec,
aIdVec,
itr_invp,
itr_lssol,
conv_res_lssol,
out iflag_conv);
// iflag_conv (0:正常に終了 1:ICCGが収束しなかった)
Console.WriteLine("////// Result of CEigenLanczos.MinimumComplexEigenValueVector_InvPower //////");
Console.WriteLine("eigenVal = " + eigenVal.Real + " + " + eigenVal.Imag + " i");
Console.WriteLine("iflag_conv = " + iflag_conv);
for (uint iblk = 0; iblk < updatePtr.BlkVecLen(); iblk++)
{
DelFEM4NetCom.Complex cval = updatePtr.GetValue(iblk, 0);
System.Console.WriteLine("eigenVec( " + iblk + " ) = " + cval.Real + " + " + cval.Imag + " i");
}
// 固有ベクトルを先頭の値を基準に位相をずらす
System.Console.WriteLine("phaseShift");
DelFEM4NetCom.Complex cval0 = updatePtr.GetValue(0, 0);
DelFEM4NetCom.Complex phaseShift = cval0 / DelFEM4NetCom.Complex.Norm(cval0);
for (uint iblk = 0; iblk < updatePtr.BlkVecLen(); iblk++)
{
DelFEM4NetCom.Complex cval = updatePtr.GetValue(iblk, 0);
cval /= phaseShift;
System.Console.WriteLine("phaseShift eigenVec( " + iblk + " ) = " + cval.Real + " + " + cval.Imag + " i");
}
// 残差チェック
DelFEM4NetCom.Complex[] evec = new DelFEM4NetCom.Complex[updatePtr.BlkVecLen()];
for (uint iblk = 0; iblk < updatePtr.BlkVecLen(); iblk++)
{
evec[iblk] = updatePtr.GetValue(iblk, 0);
}
//DelFEM4NetCom.Complex[] residualVec = new DelFEM4NetCom.Complex[updatePtr.BlkVecLen()];
DelFEM4NetCom.Complex[] leftVec = product(srcMat, evec);
double residual = 0;
for (int i = 0; i < updatePtr.BlkVecLen(); i++)
{
double squaredNorm = DelFEM4NetCom.Complex.SquaredNorm(leftVec[i] - eigenVal * evec[i]);
residual += squaredNorm;
}
residual = Math.Sqrt(residual);
Console.WriteLine("residual: {0}", residual);
計算結果
LisysのKrdLab.clapack.FunctionExt.zgeevの結果も併せて表示していますが、両者は一致していると思います。
