ガウスの掃き出し法によるC++プログラム

大学で、ガウスーヨルダンの掃き出し法により連立方程式を解き、係数マトリクスの逆行列と解を表示するプログラムを作れ、という課題が出ました。

4s+t+3u+2v=23
s+4t+3u+3v=30
5s+5t+10u+5v=65
4s+4t+2u+6v=42

という問題です。
まったく素人の状態から4時間ほどやったくらいでこの問題が出たので、解き方が全くわからなかったのですが、いろいろなサイトを巡ってこのようなプログラムを作りました。

#include <stdio.h>
#include <float.h>

#define N 4 /* 行列の行数および列数 */

double A[N][N] = { { 4.0,1.0,3.0,2.0}, /* 係数行列 A の定義 */
{1.0,4.0,3.0,3.0},
{5.0,5.0,10.0,5.0},
{4.0,4.0,2.0,6.0}};
double b[N] = {23.0,30.0,65.0,42.0}; /* 定数ベクトル b の定義 */

void Gauss_J( int, double*, double* );

void main( void )
{
int i; /* カウンタ */

printf( "%d元連立一次方程式\n", N );
for( i = 0; i < N ; i++ )
{
printf( "%g s + %g t + %g u + %g v = %g \n",
A[i][0], A[i][1], A[i][2],A[i][3], b[i] );
}
printf("の解は，\n" );
Gauss_J( N, (double *)A, (double *)b ); /* ガウス・ジョルダン法で解く */
printf( "s = %g \n", b[0] );
printf( "t = %g \n", b[1] );
printf( "u = %g \n", b[2] );
printf( "v = %g \n", b[3] );
}

void Gauss_J( int n, double *a, double *b )
{
int p, i, j, l ; /* カウンタ */
double pivot, c ; /* ピボット値 */

for ( p = 0 ; p < n ; p++ ) /* 1行目から n行目まで繰り返す */
{
pivot = a[ p*n + p ]; /* ピボットを取得する */
for ( i = p ; i < n ; i++ ) /* p行目の p列目から n列目まで */
{
a[ p*n + i ] /= pivot; /* 係数行列の p行を pivotで割る */
}
b[ p ] /= pivot; /* 定数ベクトルの p行を pivotで割る */

for ( l = 0 ; l < n ; l++ ) /* 1行目から n行目まで */
{
if ( l != p ) /* p行を除いて */
{
c = a[ l*n + p ]; /* 掃き出す */
for ( j = p ; j < n ; j++ )
{
a[ l*n + j ] -= c * a[ p*n + j ];
}
b[ l ] -= c * b[ p ];
}
}
}
return ;
}

これで行列の解は出るようになったのですが、逆行列が表示されてません。
どうすれば表示されるようになるのでしょうか？

No.2 ベストアンサー 回答者： foomufoomu 回答日時： 2011/06/10 15:06

ガウスの掃き出し方は、与えられた連立方程式を行列を、

[A][x]=[B]
の形にしておいて、A行列からピボットを取り出して割ったり引いたりして、最終的に[a]行列を単位行列（対角要素が１、それ以外が０の行列）に変形していくのですが、同時に同じ操作を[B]ベクトルにも施すと、最後には[B]が答えになっているというものです。

逆行列を求めるときは、[B]ベクトルの代わりに[A]と同じ大きさの単位行列[E]を置いて、[A]を単位行列に変形したのと同じ操作を[E]にも施すと、最後には[E]が逆行列になっています。

No.1 回答者： Tacosan 回答日時： 2011/06/10 12:28

逆行列を計算し表示するようにすれば表示されるようになります.