矩阵算法描述

矩阵乘法：

设C=A*B;

A的列数(A.n)必须等于B的行数(B.m);

C的列数=B的列数,C的行数=A的行数;

C[i][j]为A的第i行与B的第j列乘积累加;

具体实现:

for(i=0;i

for(j=0;j

for(k=0;k

C[i][j] += A[i][k] * B[k][j];

}}}

代码:

int multMatrix(matrix dat,matrix *res)//矩阵左乘

{

matrix mul;

int i,j,k;

input_matrix(&mul);

if(dat.n!=mul.m){

printf("所输入矩阵列数必须等于第一个矩阵的行数\n");

return 0;

}

res->m = dat.m;

res->n = mul.n;

for(i=0; i

for(j=0; j

res->array[i][j] = 0;//因为指针指向了一个奇怪的数值

for(k=0; k

res->array[i][j] += dat.array[i][k] * mul.array[k][j];

}

}

}

return 1;

}

矩阵求逆:

必须为方阵;

对于k从0到n-1作如下几步,n为方阵维数;

1.从第k行,第k列,开始的右下角子阵中选取绝对值最大的元素a[r][c],并记住此元素所在的行号和列号,分别以一维数组is[n],与js[n]保存, 设r为绝对值最大的行号,c为绝对值最大的列号,is[k]=r , js[k]=c, 如果a[k][k]的右下角矩阵全为0,则该矩阵无法求逆;

再通过行交换与列交换将它交换到主元素位置上.这一步称全选主元;

a[k][j]<->a[r][j]; j =0,1,2…n-1; r != k;

a[i][k]<->a[i][c]; i =0,1,2…n-1; c != k;

2.a[k][k] = 1/a[k][k]

3.a[k][j] = a[k][j]*a[k][k], j=0,1,...n-1; i,j!=k

4.a[i][j]=a[i][j]-a[i][k]*a[k][j], i,j=0,1,...,n-1; i,j!=k

5.a[i][k]=-a[i][k]*a[k][k], i=0,1,...,n-1; i!=k;

对于k从n-1到0作恢复,原则如下

在全选主元过程中,先交换的行,列后进行恢复;原来的行(列)交换用列(行)交换来恢复;

a[k][j]<->a[c][j];a[i][k]<->a[i][r]; i , j =0,1,2…n-1; c , r != k;

代码:

int inv(matrix dat,matrix *res)

{

int is[MAX];//记录绝对值最大的行号

int js[MAX];//记录绝对值最大的列号

int i,j,k;

double d,p;

if(dat.m!=dat.n){

printf("行数与列数不相等");

return 0;

}else{

//*res = dat;

for(k=0 ; k

{

d = 0.0;

for(i = k; i < dat.m ;i++)

{

for(j = k; j < dat.n ;j++)

{

p = fabs(dat.array[i][j]);//找出a[k][k]右下角绝对值最大的元素

if(p>d)

{

d = p;is[k]=i;js[k]=j;//记录绝对值最大的元素行号和列号,设为a[r][c]

}

}

}

if(d+1.0 == 1.0)//失败

{

printf("error");

return 0;

}

//通过行交换与列交换将它交换到主元素位置上

if(is[k]!=k)//如果绝对值最大的元素正好在k行上，就不用交换了

{

for(j=0; ja[r][j]

{

p = dat.array[k][j];

dat.array[k][j] = dat.array[is[k]][j];

dat.array[is[k]][j] = p;

}

}

if(js[k]!=k)//如果绝对值最大的元素正好在k列上，就不用交换了

{

for(i=0; ia[i][c]

{

p = dat.array[i][k];

dat.array[i][k] = dat.array[i][js[k]];

dat.array[i][js[k]] = p;

}

}

dat.array[k][k] = 1/dat.array[k][k];//a[k][k] = 1/a[k][k]

for(j=0; j

{

if(j!=k)//a[k][k]不要改

{

dat.array[k][j] = dat.array[k][j]*dat.array[k][k];

}

}

for(i=0; i

{

if(i!=k)

for(j=0 ; j

if(j!=k)

{

dat.array[i][j] = dat.array[i][j]-dat.array[i][k]*dat.array[k][j];

}

}

for(i=0 ; i

{

if(i!=k)

{

dat.array[i][k] = -dat.array[i][k]*dat.array[k][k];

}

}

}

for(k=dat.n-1;k>=0;k--)

{

if(js[k]!=k)

{

for(j=0;j<=dat.n-1;j++)//a[k][j]<->a[c][j]

{

p = dat.array[k][j];

dat.array[k][j] = dat.array[js[k]][j];

dat.array[js[k]][j] = p;

}

}

if(is[k]!=k)

{

for(i=0;i<=dat.m-1;i++)//a[i][k]<->a[i][r]

{

p = dat.array[i][k];

dat.array[i][k] = dat.array[i][is[k]];

dat.array[i][is[k]] = p;

}

}

}

*res = dat;

return 1;

}

}

求行列式:

必须为方阵;

对于k从0到n-1作如下几步,n为方阵维数

1. 从第k行,第k列,开始的右下角子阵中,找出绝对值最大的元素a[is][js],把其下标保存在is,js中,如果a[k][k]的右下角矩阵全为0,则行列式的值为0;

2. 再通过行交换与列交换将它交换到主元素位置上.这一步称全选主元;

a[k][j]<->a[is][j]; j =0,1,2…n-1; is != k;

a[i][k]<->a[i][js]; i =0,1,2…n-1; js != k;

交换时记录交换次数,交换次数为奇数时,要把所求得的行列式*(-1);

3. det = det * a[k][k];

4. d = a[i][k] / a[k][k] ; i = k+1,k+2,…n-1; d为第i行第k列元素a[i][k]与a[k][k]的比值;

5. a[i][j] = a[i][j] – d * a[k][j]; i,j = k+1,k+2…n-1; i行减去k行的d倍数(第j列);

结束k的循环

det = det * a[n-1][n-1];//若交换次数为奇数,则*-1;

det为返回的行列式值;

代码:

double det(matrix dat)

{

int i,j,k,is = 0,js = 0 ;

double f,det,d,q;

f=1.0;det=1.0;

if(dat.m!=dat.n)

{

printf("行数与列数不相等");

return 0;

}

for(k=0; k<=dat.m-2; k++)

{

q=0.0;//从a[k][k]的右下角开始，找出绝对值最大的元素,把其下标保存在is,js中for(i=k; i

{

for(j=k; j

{

d = fabs(dat.array[i][j]);

if(d>q)

{

q=d;is=i;js=j;

}

}

}

//如果右下角矩阵全为0,则行列式值为0

if(q==0.0)

{

det =0.0;

return det;

}

//不是当前行,则置换第k行和第is行

if(is != k)

{

f = -f;//f记录交换次数,为奇数时,返回结果为负

for(j=k; j

{

d = dat.array[k][j];

dat.array[k][j] = dat.array[is][j];

dat.array[is][j] = d;

}

}

//不是当前列,则置换第k行和第js行

if(js != k)

{

f = -f;

for(i=k; i

{

d = dat.array[i][js];

dat.array[i][js] = dat.array[i][k];

dat.array[i][k] = d;

}

}

det = det*dat.array[k][k];

for(i=k+1;i

{

d=dat.array[i][k]/dat.array[k][k];//d为第i行第k列元素a[i][k]与a[k][k]的比值

for(j=k+1; j

{

dat.array[i][j]=dat.array[i][j]-d*dat.array[k][j];//i行减去k行的d倍数(第j 列)

}

}

}

det = f*det*dat.array[dat.m-1][dat.n-1];

return det;

}

矩阵的轶:

设m为行,n 为列

for(ci=0; ci

for(i=ri; i < m; i++){

当a[i][ci]!=0,且i!=ri,执行下面几步:

1. 从第ci列开始交换ri行与i行(fz为1)

2. a[i][ci]=a[i][ci]*a[ri][ci]; ci=ci,ci+1,…n-1;

a[i][ci] += a[ri][ci] * -a[i][ci]; ci = ci,ci+1,…,n-1;

(fz为0)

3. fz=0; (无论i是否等于ri)

}ri++; (fz==0)

}返回结果为轶ri

代码:

//交换两行

void exchang_row(double *a,double *b,int n)

{

int i;

double t;

for(i=0;i

{ t=a[i]; a[i]=b[i]; b[i]=t;}

}

//将某一行乘以k

void mul_row(double *a,double k,int n)

{

int i;

for(i=0;i

a[i]*=k;

}

//将a2行的k倍加到a1行上

void add_row(double *a1,double *a2,int k,int n)

{

int i;

for(i=0;i

a1[i]+=a2[i]*k;

}

int rank_matrix(matrix dat,matrix *res) //求秩(返回值为秩，第2个参数为变换得到的阶梯阵) {

int m,n,i;;

double t;

int ri,ci; //行标记与列标记

int f_z; //某行是否全为0的标志，为1表示全为0

m=dat.m;n=dat.n;

for(ri=ci=0;ci

{

f_z=1;

for(i=ri;i

{

if(dat.array[i][ci]!=0)

{

if(i!=ri)

if(f_z)

exchang_row(&(dat.array[ri][ci]),&(dat.array[i][ci]),n-ci);

else

{

t=dat.array[i][ci];

mul_row(&(dat.array[i][ci]),dat.array[ri][ci],n-ci);

add_row(&(dat.array[i][ci]),&(dat.array[ri][ci]),-t,n-ci);

}

f_z=0;

}

}

if(!f_z) ri++;

}

*res=dat;

return ri; }