DelFEMで行列の固有値を求める（２）

DelFEMの固有値解析ライブラリMinimumEigenValueVector_InvPowerをC#に移植してみました。

前回のテストプログラムでは、 DelFEM4NetLsSol.CEigenLanczos. MinimumEigenValueVector_InvPower というメソッドを使用していますが、これはDelFEMの関数を呼び出すだけのラッパーで、具体的な処理は下記のように何もありません。

 DelFEM4Net/src/ls/eigen_lanczos.cpp

#include "DelFEM4Net/stub/clr_stub.h"
#include "DelFEM4Net/ls/eigen_lanczos.h"
#include "DelFEM4Net/ls/preconditioner.h"
#include "DelFEM4Net/ls/linearsystem.h"
#include "DelFEM4Net/ls/solver_ls_iter.h"

using namespace DelFEM4NetLsSol;

double CEigenLanczos::MinimumEigenValueVector_InvPower(
        DelFEM4NetLsSol::CLinearSystem^ ls,
        DelFEM4NetLsSol::CPreconditioner^ pls,じ
        IList^ aIdVec, 
        unsigned int itr_invp,    // 逆べき乗法の最大反復回数
        unsigned int itr_lssol,   // ICCG法の最大反復回数
        double conv_res_lssol,    // ICCG法の収束基準の相対残差
        [Out] int% iflag_conv )   // 0:正常に終了  1:ICCGが収束しなかった
{
    std::vector aIdVec_;
    DelFEM4NetCom::ClrStub::ListToVector(aIdVec, aIdVec_);
    int iflag_conv_ = 0;
    
    double ret = LsSol::MinimumEigenValueVector_InvPower(
        *(ls->Self),
        *(pls->PrecondSelf),
        aIdVec_,
        itr_invp,
        itr_lssol,
        conv_res_lssol,
        iflag_conv_);
        
    iflag_conv = iflag_conv_;
    
    return ret;
}

LsSol:: MinimumEigenValueVector_InvPower というのがDelFEM本体の関数です。

今回、この関数の中身をC#に移植してみました。

逆べき乗法による固有値解析ソースコード

(DelFEM\/src/ls/eigen_lanczos.cppのライブラリを移植)

/*
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
    {
        ///
        /// 逆べき乗法による最小固有値の取得
        ///
        ///リニアシステムオブジェクト
        ///プリコンディショナ―オブジェクト
        ///無視ベクトル
        ///逆べき乗法の最大反復回数
        ///ICCG法の最大反復回数
        ///ICCG法の収束基準の相対残差
        ///0:正常に終了  1:ICCGが収束しなかった
        /// 
        public static double MinimumEigenValueVector_InvPower(
            DelFEM4NetLsSol.CLinearSystem ls,
            DelFEM4NetLsSol.CPreconditioner pls,
            IList 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);
            }
            {  // 更新を正規化
                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);
                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);
                }
                // レリー積の計算
                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;
        }
    }        
}

使用方法は、前回の固有値解析ソースのDelFEM4NetLsSol.CEigenLanczos. MinimumEigenValueVector_InvPower をDelFEM4NetLsSol.CEigenLanczos_CSharp. MinimumEigenValueVector_InvPowerに書き換えるだけです。同じ結果が得られます。

DelFEMの処理をほぼそのままコピーして作成したのですが、DelFEM本体の低レベルI/Fの使い方の参考になりそうです。