!!!MM行列の逆行列を計算

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

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

なお、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 行でデータが格納されます（行数と列数は同じ）。