部分选主元matlab,部分选主元的Doolittle分解 | 学步园

步骤：

假设用紧凑格式的Doolittle法已经完成了第 r-1 (1<=r<=n) 步分解，第 r 步分解，首先在数组 A 的第 r 列主对角元以下(含主对角元)

选主元，具体步骤：

1、计算中间量 Si ,并存入 A( i , r ) ( r = 1 : n-1 ).

A( i , r )

(i = r , ... , n ).

2、挑选绝对值最大的 Si，即确定行号 ir , 满足

| Sir | = max ( | Si | , r<=i<=n )

3、换行。如果 ir != r ,则交换数组 A 中的第 r 行与 ir 行的对应元素。交换位置的元素以它在 A 中所处的新位置记。此时

urr = A( r, r ) .

4、分解计算

A( i , r )

A( r , r)

按上述方法完成 n 步分解，便实现了系数矩阵分解 PA=LU.

C语言代码(未考虑无解情况)：

#include

#include

#include

#define M 128

double mat[M][M];

double ans[M],y[M];

int n;

bool Read()

{

int i,j;

if(1!=scanf("%d",&n)) return false;

for(i=0;i

for(j=0;j<=n;j++)

scanf("%lf",&mat[i][j]);

}

return true;

}

void Print()

{

int i,j;

for(i=0;i

for(j=0;j<=n;j++)

printf("%.2f/t",mat[i][j]);

puts("");

}

puts("");

}

void Work()

{

int i,k,r,t;

double d,tmp[M];

for(r=0;r

//

for(i=t=r;i

for(k=0,d=0.0;k

d+=mat[i][k]*mat[k][r];

mat[i][r]-=d;

if(fabs(mat[i][r])>fabs(mat[t][r]))

t=i;

}

Print();

memmove(tmp,mat[t],sizeof(tmp));

memmove(mat[t],mat[r],sizeof(tmp));

memmove(mat[r],tmp,sizeof(tmp));

Print();

for(i=r+1;i

mat[i][r]/=mat[r][r];

for(i=r+1;i<=n;i++){

for(d=0.0,k=0;k

d+=mat[r][k]*mat[k][i];

mat[r][i]-=d;

}

Print();

}

}

void Result()

{

int r,i;

double tmp;

for(r=n-1;r>=0;r--){

for(i=r+1,tmp=0.0;i

tmp+=mat[r][i]*ans[i];

ans[r]=(mat[r][n]-tmp)/mat[r][r];

}

for(i=0;i

printf("%lf/t",ans[i]);

puts("");

}

int main()

{

while(Read()){

Work();

Print();

Result();

}

return 0;

}

/*

3

1 -1 3 1

2 -4 6 4

4 -9 2 1

*/