ryujimiyaの日記

C#を使って数値解析したい

clapackのラッパー(KrdLab.Lisys CLW.dll)を更新しました

KrdLab.Lisysのryujimiya追加分を更新しました。

HPlaneWGSimulator X DelFEMのversion1.0.0.2、HPlaneWGSimulatorのversion1.2.0.9で取り入れましたのでよろしければご利用ください。

  KrdLab Lisys

   ryujimiya追加分  ※本家は、こちら

変更したのは複素数連立方程式求解関数(zgesv)です。

いままでこの関数内部でメモリの再確保を行っていたのですがそれを解消しました。

どうしてメモリの再確保が必要だったかというと、C#のSystem.Numerics.Complexからclapackの複素数構造体doublecomplexへ渡すとき、両者のメモリ配置が全く異なるのでdoublecomplexを新たに確保するしかなかったからです。

いままではこんな感じ。

            zgesv(
                     ....
                array<System::Numerics::Complex>^  A,
                int  a_row,
                int  a_col,
                     ....
                )
            {
                
                // Aをnativeの複素数構造体に変換
                doublecomplex* a = nullptr;
                try
                {
                    a = new doublecomplex[a_row * a_col];
                    for (int i = 0; i < a_row * a_col; i++)
                    {
                        a[i].r = A[i].Real;
                        a[i].i = A[i].Imaginary;
                    }
                }
                catch (System::Exception^ /*exception*/)
                {
                    if (a != nullptr)
                    {
                        delete[] a;
                        a = nullptr;
                    }
                    throw;
                }
            }

今回、KrdLab.clapack.Complexというvalue classを追加しました。

        /// <summary>
        /// 複素数構造体
        ///   clapackのdoublecomplexに対応します
        ///   pin_ptrを介してdoublecomplex型でアクセスするのがこの構造体の主要用途です
        /// </summary>
        public value class Complex
        {
        private:
            /// <summary>
            /// 実数部
            /// </summary>
            double r;
            /// <summary>
            /// 虚数部
            /// </summary>
            double i;
        public:
            /// <summary>
            /// コンストラクタ
            /// </summary>
            /// <param name="r_">実数部</param>
            /// <param name="i_">虚数部</param>
            Complex(double r_, double i_)
            {
                r = r_;
                i = i_;
            }
            /// <summary>
            /// 実数部
            /// </summary>
            property double Real
            {
                double get() { return r; }
                void set(double value) { r = value; }
            }
            /// <summary>
            /// 虚数部
            /// </summary>
            property double Imaginary
            {
                double get() { return i; }
                void set(double value) { i = value; }
            }
            <以下省略>
        }

 一方、doublecomplexは

        typedef struct { doublereal r, i; } doublecomplex;

 なので、pin_ptrを使えばKrdLab.clapack.Complexからdoublecomplexへキャストできます。

            array<KrdLab.clapack.Complex>^ X;
            
            pin_ptr<void> a_ptr = &X[0];
            doublecomplex* a = (doublecomplex *)(void *)a_ptr;

 これでメモリの再確保をしなくて済んだので、アプリ側でメモリ確保に成功すればclapackはエラーにならず解いてくれるようになりました。