手計算でやっていると変数が増えたときに苦労しますが，行列を使うと計算機ですぐに解を求められます．
このような連立一次方程式があったときに，次のような行列を考えます．
ここで行列の行基本変形を使って変形します．行基本変形とは，
です．これらを使って次のように変形します．この変形を前進消去といいます．
これは式で現すと，
と同じ事なので，後は割算と代入を下から順に行えば解を求めることができます．これを後退代入といいます．
前進消去は上の例でいうと一行目を巧くα倍して二行目，三行目から引くことで一列目の二行目以降をゼロにし，二行目を巧くα倍して三行目から引くことで二列目の三行目をゼロにしています．
この事から，計算途中で対角成分がゼロになってしまうと解が求まらない事がわかります．そのため，巧く行を入れ換え対角成分がゼロにならないようにしがら計算しています．それでも対角成分がゼロになってしまうときは解が存在しないもしくは不定ということになります．
ソースコードを見てみます．計算にあたってベクトルや行列を扱うために構造体を宣言し，領域確保や初期化するための関数と開放するための関数を作成しています．
typedef struct{
int size;
double *v;
}Vector;
typedef struct{
int height;
int width;
double *a;
}Matrix;
Vector *CreateVector(int size);
void FreeVector(Vector* vector);
Matrix *CreateMatrix(int width, int height);
void FreeMatrix(Matrix* matrix);
それでは前進消去です．
Matrix *temp; temp = CreateMatrix(matrix->width, matrix->height); for(i=0; i<temp->height*temp->width; i++) temp->a[i] = matrix->a[i];
アルゴリズムには関係ないのですが，解が存在しないときにreturnしてしまうため行列の値が壊れないようにコピーして計算しています．
if(temp->height > temp->width) size = temp->width; else size = temp->height;
連立一次方程式を解くときに縦長の行列になることは無いと思うのですが，念の為に対処しています．
d = i;
max = fabs(temp->a[i*temp->width + i]);
for(j=i+1; j<temp->height; j++){
if(max < fabs(temp->a[j*temp->width + i])){
d = j;
max = fabs(temp->a[j*temp->width + i]);
}
}
for(j=i; j<temp->width; j++){
copy = temp->a[i*temp->width + j];
temp->a[i*temp->width + j] = temp->a[d*temp->width + j];
temp->a[d*temp->width + j] = copy;
}
解が存在するときの連立一次方程式である場合に，対角成分がゼロにならないように対処しています．例えば二列目であればニ列目の二行目以降の中で絶対値が最大である成分を持つ行番号をdに確保し，その行と2行目を交換します．
後退代入に関しては特に難しい部分は無いと思います．
