MM行列の逆行列を計算_clapack

MM行列の逆行列を計算

行と列の要素数が同じ（ここではM個とする）場合に限り、CLAPACKの「dgetrf_」と「dgetri_」の２つを使うことで逆行列の計算が可能です。

行と列の要素数が異なる場合は、擬似逆行列?を使うようにします。

なお、CLAPACKのヘッダファイルである「f2c.h」「clapack.h」を忘れずにインクルードするようにしてください。また、インクルードの指定順は「f2c.h」「clapack.h」の順に指定しないとエラーになります。

#include "f2c.h"
#include "clapack.h" 

/**
 *  M x M 行列の逆行列を計算する
 *  @param[in]  pM     対象のM x N行列（列行の順である点に注意） 
 *  @param[in]  mCou   行列Mの行数
 *  @param[in]  nCou   行列Mの列数
 *  @param[out] pRetM  計算後の行列要素が返る配列（列行の順である点に注意）
 *  @return 処理に成功すればtrueが返る
 */
bool MatrixInverse(double *pM, const int mCou, const int nCou, double *pRetM)
{
   integer m, n, lda, info;
   integer piv[500];
   integer lwork;
   double *pMat;
   double work[500];

   if(mCou != nCou) return false;

   m     = mCou;
   n     = nCou;
   lda   = m;
   lwork = 500;

   pMat = (double *)alloca(sizeof(double) * mCou * nCou);

   // 行列要素を作業バッファにコピー
   memcpy(pMat, pM, sizeof(double) * mCou * nCou);

   // 逆行列の計算
   dgetrf_(&m, &n, pMat, &lda, piv, &info);
   if(info != 0) return false;
   dgetri_(&n, pMat, &lda, piv, work, &lwork, &info);
   if(info != 0) return false;

   // 行列をpRetMにコピー
   memcpy(pRetM, pMat, sizeof(double) * mCou * nCou);

   return true;
}

途中で作業用の領域piv、work とそのサイズlworkが出てきますが、適当にmCou * nCouより十分大きな値を入れてます。計算途中で容赦なく変数の値を書き換えるようですので、一度pMatというローカルメモリを確保してそこに計算データを格納するようにしています。

結果のpRetMatには列 x 行でデータが格納されます（行数と列数は同じ）。