MKLのLAPACKでLU分解

MKLで逆行列を求めるには、行列を一度LU分解しておく必要がある。
ダウンロードしてきたMKLのC#のラッパーには、LAPACKのLU分解パッケージ"dgetrf"がなかったので、記載されていたものを参考に作成してみた。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Security;
using System.Runtime.InteropServices;

namespace mkl
{
    public sealed class LAPACK
    {
        private LAPACK() { }

        /* Constants for LAPACK_JOB */
        public sealed class JOB
        {
            private JOB() { }
            public static int JobN = 201;      /* job ='N' */
            public static int JobV = 202;      /* job ='V' */
            public static int JobA = 203;      /* job ='A' */
            public static int JobS = 204;      /* job ='S' */
            public static int JobO = 205;      /* job ='O' */
        };


        // LAPACK dgetrf wrapper
        public static int dgetrf(int m, int n, double[] a, int lda,int[] ipiv)
        {
            int info = 0;
            //int[] ipiv = new int[Math.Min(m,n)];

            LAPACKNative.dgetrf(ref m, ref n, a, ref lda, ipiv, ref info);

            return info;
        }
    }

    /** LAPACK native declarations */
    [SuppressUnmanagedCodeSecurity]
    internal sealed class LAPACKNative
    {
        private LAPACKNative() { }

 
        // LAPACK native dgetrf declaration
        [DllImport("mkl_rt", CallingConvention = CallingConvention.Cdecl,
            ExactSpelling = true, SetLastError = false)]
        internal static extern void dgetrf(
            ref int m,                  //行列の行数(m)
            ref int n,                  //行列の列数(n)
            [In, Out] double[] a,       //行列へのポインタ
            ref int lda,                //行列の行数
            [In, Out] int[] ipiv,       //ピボット選択の添字を格納する長さ(min(m,n))の配列
            ref int info                //実行後，例外情報が格納
            );
    }
}

でもって、これを呼ぶプログラムを作ってみた。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using mkl;

namespace dgetri_test
{
    class test_dgetri
    {
        static void Main(string[] args)
        {
            Console.WriteLine("--- dgetrf test---");
            int info;
            double[] A = new double[] {3,1,4,3,2,4,6,6,6,};
            int m = 3, n = 3;
            int lda = m;
            int[] ipiv = new int[Math.Min(m, n)];
            
            //print the parameters
            printMatrix("Matrix A:",A, m, n);

            //compute
            info = LAPACK.dgetrf(m, n, A, lda, ipiv);

            //print result
            printMatrix("Matrix A on exit:",A,m,n);
            printVector("Vector P on exit:",ipiv,m);
            Console.WriteLine("info on exit:" + info);
            Console.WriteLine("TEST PASSED");
            Console.WriteLine();
        }

        /** Print the matrix X assuming row-major order of elements. */
        private static void printMatrix(String prompt, double[] X, int I, int J)
        {
            Console.WriteLine(prompt);
            for (int i = 0; i < I; i++)
            {
                for (int j = 0; j < J; j++)
                    Console.Write("\t" + X[i * J + j]);
                Console.WriteLine();
            }
        }
        /** Print the vector X */
        private static void printVector(string prompt, int[] X, int N)
        {
            Console.WriteLine(prompt);
            for (int n = 0; n < N; n++)
                Console.Write("\t" + X[n]);
            Console.WriteLine();
        }
    }
}

結果を見てみると、「下三角行列L」と、「上三角行列Uの対角成分以外」が格納されているようだ。しかもこの行列について、行列の積LUを計算すると適当な置換がかかっている。

dgetrfの呼出し方がおかしいのか、それともこういうもんなのか謎。