加载中…矩阵求逆——高斯约旦法（代码示例）
(2011-06-13 21:06:14)
标签：
矩阵求逆示例it |
分类： 编程经验--积跬步至千里 |
float inv(float *a, int row, int col, float *b){
int i,j,k;
int sign=1;
float fDet=1.0; //保存行列式的值
int *maxRow=(int
*)malloc(sizeof(int)*row);//存储下面全选主元过程中最大元素的坐标
int *maxCol=(int
*)malloc(sizeof(int)*row);//矩阵a应为方阵。所以row=col
float temp, fMax=0;
//copy matrix a to matrix b
for(i=0;i<row;++i)
for(j=0;j<col;++j)
b[i*col+j]=a[i*col+j];
for(i = 0; i< col; ++i){
maxRow[i] = i;
maxCol[i] = i;
}
for(k=0;k<row;++k){
//全选主元，即选取以b[k][k]为左上角的矩阵中最大的元素
fMax=0.0;
for(i=k;i<row;++i)
for(j=k;j<col;++j){
temp=(float)fabs(b[i*col+j]);
if(temp>fMax){
maxRow[k]=i;
maxCol[k]=j;
fMax=temp;
}
}
if(fMax<1e-30) return 0.0;
//行交换、列交换
if(maxRow[k]!=k){ swapRows(b, row, col, maxRow[k], k);
sign=-sign;}
if(maxCol[k]!=k){ swapCols(b, row, col, maxCol[k], k);
sign=-sign;}
//计算行列式的值
fDet *= b[k*col+k];
b[k*col+k]=1.0f/b[k*col+k];
for(j=0;j<col;++j)
if(j!=k)b[k*col+j] *= b[k*col+k];
for(i=0;i<row;++i)
for(j=0;j<col;++j)
if(i!=k&&j!=k)
b[i*col+j] -= (b[i*col+k]*b[k*col+j]);
for(i=0;i<row;++i)
if(i!=k)
b[i*col+k]*=-(b[k*col+k]);
}
for(k=col-1;k>=0;--k){
if(maxCol[k]!=k) swapRows(b, row, col, maxCol[k], k);
if(maxRow[k]!=k) swapCols(b, row, col, maxRow[k], k);
//先考虑列交换，如果发生过，要用相应的行交换回去。而且这里的k要从大到小。
}
free(maxRow);
free(maxCol);
return fDet * sign;
}
喜欢
0
赠金笔