LU分解を利用して逆行列を求めたい

お世話になっております。

C++でLU分解と連立一次方程式の解を求める関数SLE_bu_LUを利用して、逆行列の求める関数InvMatを作成中です。

「行列Bのj列目の列ベクトルをbj とすれば，行列Cのj列目の列ベクトル
xj は連立１次方程式 Axj = bj の解として求められる．これが各列について成立するので，
求めたx1, x2, x3，････, xmを各列としてC=[x1 x2 x3 ････ xm] を計算することがで
きる。その際，係数行列Aが同一の連立１次方程式をm 組解いて行列C が得られること
がわかる。従って，初めに１回だけAを LU 分解しておいて，前進代入と後退代入のプロセ
スをm回行うことによって行列C が求められる」
と議題には書かれており、Cを求める逆行列、Bを単位行列としてプログラミングを書いていたのですが、
関数InvMatではCとBのm列成分を関数SLE by LUに引数として渡すという方針であっていますでしょうか？
（そうなるとm回関数SLE_by_LUを呼び出し、m回LU分解してしまうので間違っていると思いました...）

関数SLE_bu_LUは正しく動き、追加した関数は一番下のInvMatのみでここが怪しいと思います。

よろしくお願いします。

#include <iostream> 
#include <iomanip> 
using namespace std; 
const int    N = 4;                    /*  連立方程式の元数      */ 
void SLE_by_LU(double *,double [][N],double *); // LU 分解により連立１次方程式 
                                                 // を解く関数のプロトタイプ宣言 
void InvMat(double [][N],double [][N]);         //逆行列を求める関数のプロトタイプ宣言
int main()  
{ 
    double   x[N]; 
   /*    入力データ （演習問題 2.1 の場合）  */ 
   double a[N][N]={{1.0,2.0,-12.0,8.0},{5.0,4.0,7.0,-2.0}, 
                   {-3.0,7.0,9.0,5.0},{6.0,-12.0,-9.0,3.0}}; 
   double ai[N][N]={}; //逆行列用
    double b[N]={27.0, 4.0,11.0,49.0}; 



SLE_by_LU(x,a,b);  // LU 分解により連立１次方程式を解く関数の呼び出し 


   /*   データ出力     */ 
   for(int i=0; i<N; i++)  {cout << 'x' << i+1 << '=' << setw(12) << 
                             fixed << setprecision(9) << x[i] << endl;
   }
   InvMat(a,ai);//逆行列の計算関数呼び出し

 return 0; 
}
/******  LU 分解により連立１次方程式を解く関数プログラム  **********/ 
/*      入力 ： a[N][N]=係数行列，b[N]=右辺ベクトル，N=元数     */ 
/*      出力 ： x[N]=方程式の解                                 */ 
/******************************************************************/ 


void SLE_by_LU(double x[],double a[][N],double b[]) 
{    int  i, j, k; 
     double  l[N][N], u[N][N], y[N];                    // LU 分解用配列 
     double  sigma;                                     // ∑演算用作業変数 
     /*      LU 分解     */ 
     for(j=0; j<N; j++) u[0][j]=a[0][j];                // ステップ 1.1 
     for(j=1; j<N; j++) l[j][0]=a[j][0]/u[0][0];        // ステップ 1.2 
     for(k=1; k<N; k++)                                 // ステップ 2 の反復 
     { 
        for(j=k; j<N; j++)                              // ステップ 2.1 
        { 
           u[k][j]=a[k][j]; 
           for(i=0; i<k; i++)  u[k][j]-=l[k][i]*u[i][j]; 
       } 
       for(j=k+1; j<N; j++)                            // ステップ 2.2 
       { 
          sigma=a[j][k]; 
          for(i=0; i<k; i++)  sigma-=l[j][i]*u[i][k]; 
          l[j][k]=sigma/u[k][k]; 
       } 
    } 
     
    /*     前進代入     */ 
    for(i=0; i<N; i++)                                 // i=1,2,･･･,n 
    { 
       y[i]=b[i];                                      // yi ← bi 
       for(j=0; j<i; j++)  y[i]-=l[i][j]*y[j];         // yi=bi-∑lijyj 
    } 
 
    /*     後退代入     */ 
    for(k=N-1; k>=0; k--) 
    { 
       sigma=y[k]; 
       for(j=k+1; j<N; j++) sigma-=u[k][j]*x[j]; 
       x[k]=sigma/u[k][k]; 
    } 
} 
void InvMat(double a[][N], double ai[][N]) //A　配列a[][N]に格納　逆行列AI　配列ai[][N]に格納
{//AC=BでBが単位行列Eの時、Cは逆行列AIとなる

 double array[N][N];//単位行列格納
 double array_m[N];//単位行列のm列目格納
 double ai_m[N];//AIのm列目格納
 for(int i = 0; i <N; i++){//単位行列を作る
        for(int j = 0; j <N; j++){
            if(i == j )
                {array[i][j] = 1.0;}
            else
                {array[i][j] = 0.0;}
        } 
    }

    //以下A.14より

for(int i=0;i<N;i++){//m列目処理
  for(int j=0;j<N;j++){
  array_m[i]=array[j][i];//m番目の単位行列の列成分
    ai_m[i]=ai[j][i];//m番目のAI行列の列成分
 SLE_by_LU(ai_m,a,array_m);
 cout<<ai[i][j];
  }
}
}