高斯-约旦法 求矩阵的逆（急）

[lhj0532]

我在网上找了一个高斯-约旦法求矩阵的逆，结果不对，请问他的原理是什么？一些程序错在什么地方？
void inverse(double n[4][4],double m[4][4]){
int is[4]={0,0,0,0};
int js[4]={0,0,0,0};
double fDet = 1.0;
double f =1.0;
for(int x=0;x<4;x++)
for(int y=0;y<4;y++)
m[x][y]=n[x][y];
for (int k = 0; k < 4; k ++){
double fMax = 0.0;
//戞堦庡尦
for (int i = k; i < 4; i ++){
for (int j = k; j < 4; j ++){
double f = abs(m[i][j]);
if (f > fMax){
fMax = f;
is[k] = i;//嵟戝擵強嵼揑峴
js[k] = j;//嵟戝擵強嵼揑楍
}
}
}
cout<<endl<<"is["<<k<<"]="<<is[k];
cout<<endl<<"js["<<k<<"]="<<js[k]<<endl;
if (is[k] != k){//晄嵼峴
f = -f;
//峴??
swap(m[k][0],m[is[k]][0]);
swap(m[k][1],m[is[k]][1]);
swap(m[k][2],m[is[k]][2]);
swap(m[k][3],m[is[k]][3]);
}
if (js[k] != k){ //晄嵼?楍
f = -f;
//楍??
swap(m[0][k], m[0][js[k]]);
swap(m[1][k], m[1][js[k]]);
swap(m[2][k], m[2][js[k]]);
swap(m[3][k], m[3][js[k]]);
}
fDet *= m[k][k];
m[k][k] = 1.0/m[k][k];
for (int j= 0;j<4;j++){
if (j!= k)
m[k][j] *= m[k][k];
}
for (i = 0; i < 4; i ++)
if (i != k)
for (int j = 0; j < 4; j ++)
if (j != k)
m[i][j] = m[i][j] - m[i][k] * m[k][j];
for (i = 0; i < 4; i ++){
if (i != k)
m[i][k] *= -m[k][k];
}
}
for (k = 3; k >= 0; k --){
if (js[k] != k){
swap(m[k][0], m[js[k]][0]);
swap(m[k][1], m[js[k]][1]);
swap(m[k][2], m[js[k]][2]);
swap(m[k][3], m[js[k]][3]);
}
if (is[k] != k){
swap(m[0][k], m[0][is[k]]);
swap(m[1][k], m[1][is[k]]);
swap(m[2][k], m[2][is[k]]);
swap(m[3][k], m[3][is[k]]);
}
}
}

[layman2008]

高斯-约旦法（全选主元）求逆的步骤如下：

首先，对于 k从0到n - 1 作如下几步：

从第 k 行、第 k 列开始的右下角子阵中选取绝对值最大的元素，并记住次元素所在的行号和列号，在通过行交换和列交换将它交换到主元素位置上。这一步称为全选主元。
m(k, k) = 1 / m(k, k)
m(k, j) = m(k, j) * m(k, k)，j = 0, 1, ..., n-1；j != k
m(i, j) = m(i, j) - m(i, k) * m(k, j)，i, j = 0, 1, ..., n-1；i, j != k
m(i, k) = -m(i, k) * m(k, k)，i = 0, 1, ..., n-1；i != k
最后，根据在全选主元过程中所记录的行、列交换的信息进行恢复，恢复的原则如下：在全选主元过程中，先交换的行（列）后进行恢复；原来的行（列）交换用列（行）交换来恢复。

[layman2008]

float Inverse(CLAYMATRIX& mOut, const CLAYMATRIX& rhs)
{
using std::swap;

CLAYMATRIX m = rhs;

int niSign[4];
int njSign[4];
float fDet = 1;
int nSgn = 1;

for (int k = 0; k < 4; ++ k)
{
// 第一步，全选主元
float fMax = 0;
for (int i = k; i < 4; ++ i)
{
for (int j = k; j < 4; ++ j)
{
const float fTmp = fabs(m(i, j));
if (fTmp > fMax)
{
fMax = fTmp;
is[k] = i;
js[k] = j;
}
}
}
if (fabs(fMax) < 0.0001f)
return 0;

if (niSign[k] != k)
{
nSgn = -nSgn;
swap(m(k, 0), m(niSign[k], 0));
swap(m(k, 1), m(niSign[k], 1));
swap(m(k, 2), m(niSign[k], 2));
swap(m(k, 3), m(niSign[k], 3));
}
if (njSign[k] != k)
{
nSgn = -nSgn;
swap(m(0, k), m(0, njSign[k]));
swap(m(1, k), m(1, njSign[k]));
swap(m(2, k), m(2, njSign[k]));
swap(m(3, k), m(3, njSign[k]));
}

// 计算行列值
fDet *= m(k, k);

// 计算逆矩阵

// 第二步
m(k, k) = 1.0f / m(k, k);
// 第三步
for (int j = 0; j < 4; ++ j)
{
if (j != k)
m(k, j) *= m(k, k);
}
// 第四步
for (int i = 0; i < 4; ++ i)
{
if (i != k)
{
for (int j = 0; j < 4; ++ j)
{
if (j != k)
m(i, j) = m(i, j) - m(i, k) * m(k, j);
}
}
}
// 第五步
for (int i = 0; i < 4; ++ i)
{
if (i != k)
m(i, k) *= -m(k, k);
}
}

for (int k = 3; k >= 0; -- k)
{
if (njSign[k] != k)
{
swap(m(k, 0), m(njSign[k], 0));
swap(m(k, 1), m(njSign[k], 1));
swap(m(k, 2), m(njSign[k], 2));
swap(m(k, 3), m(njSign[k], 3));
}
if (niSign[k] != k)
{
swap(m(0, k), m(0, niSign[k]));
swap(m(1, k), m(1, niSign[k]));
swap(m(2, k), m(2, niSign[k]));
swap(m(3, k), m(3, niSign[k]));
}
}

mOut = m;
return fDet * f;
}

[lhj0532]

多谢大家！！！
我找到了毛病所在，初始化时有错误：
int is[4]={0,0,0,0};
int js[4]={0,0,0,0};
应该是
int is[4]={0,1,2,3};
int js[4]={0,1,2,3};
还是感谢大家！！！