一道大学里的题，请大家帮忙看一看

[guowuji]

109 方阵A的Crout分解
1 功能 对方阵A作LU分解,其中L为下三角阵, U为单位上三角阵.
2 数学方法简介 对奇异方阵A通常是消去主对角线以下的元素,使A变为下三角阵,记未变换的A为A(0) .
第一步, 若aij(0) 0, 令uij =aij(0) /a11(0) , j=2,3,…..,n .构造初等方阵
1 u12 u13 . . .u1n
U1= 1
1
.
1
则
1 -u12 -u13 . . . –u1n
U1(-1) = 1
1
1
用u1(-1) 右乘A,
a11(0) 0 0 0
AU1(-1) = a21(0) a22(1) a23(1) … a2n(1)
… …
an1(0) an2(1) an3(1) … a nn(1)
第二步,若a22(1) 0, 令u2j =a2j(1) /a22(1) , j=3,4,……,n. 构造初等方阵
1 0 0 … 0
u2(-1)= 1 -u23 -u2n
1
1
用u2(-1) 右乘A(1),则
a11(0) 0 0 … 0
A(1)U2(-1)=A(0) U1(-1)U2(-1)= a21(0) a22(0) 0 … 0
a31(0) a32(1) a33(2) a3n(2)
… … …
an1(0) an2(1) an3(2) ann(2)
若a33(2)0,可按上述步骤继续往下作,akk(k-1) 0,k=1,2,……,n,作完n-1步,A U1(-1) U2(-1)…Un-1(-1)=A(n-1) =L
则L为下三角阵,令U(-1) =U1(-1)U2(-1)…Un-1(-1) ,它为单位上三角阵,则
A=A(n-1) U=LU
这里L为下三角阵,U为单位上三角阵.
为了省去一些工作单元,用紧凑格式编写程序,即将L存放在A的左下方(含对角线),将单位上三角阵U存放在A的右上方(不含对角线元,对角线元全为1不必存放).
上述变换过程的计算方法步骤为,对I=1,2,…,n, 有
lji =aji –∑ljk uki , j=i,i+1,…,n
uij=(aij –∑ljk ukj)/ lii , j=i +1,i+2,…,n
每步先算i列的元素,再算i行的元素.
3 使用说明
输入参数
N 整变量,方阵A的阶数.
A(N,N) N*N个元素的二维实数组,按行存放矩阵A.
输出参数
W 标志,W=0分解完毕,W=1分解进行不下去.
L(N,N) 下三角阵L;
U(N,N) 单位上三角阵U.
10 ‘*************************************’
20 ‘* 109 方阵的A的CROUT分解 *’
30 ‘*************************'************’
40 INPUT "N,E="; N, E
50 DIM A(N, N)
60 PRINT TAB(3); "EXAMPLE"; TAB(8); "JUZHEN A": PRINT
70 FOR I = 1 TO N
80 FOR J = 1 TO N
90 READ A(I, J): PRINT USING; "####.#"; A(I, J);
100 NEXT J
110 PRINT
120 NEXT I
130 PRINT: GOSUB 500
140 PRINT TAB(5); "W="; W: PRINT
150 IF (W = 1) THEN GOTO 380
160 PRINT TAB(15); "JUZHEN L"
170 PRINT
180 FOR I = 1 TO N
190 FOR J = 1 TO I
200 PRINT USING; "##.#####"; A(I, J);: PRINT " ";
210 NEXT J
220 PRINT
230 NEXT I
240 FOR I = 1 TO N
250 A(I, I) = 1
260 NEXT I
270 PRINT
280 PRINT TAB(15); "JUZHEN U"
290 PRINT
300 FOR I = 1 TO N
310 FOR J = I TO N
320 PRINT TAB(10 * (J - 1)); USING; "##.#####"; A(I, J);
330 NEXT J
340 PRINT
350 NEXT I
360 DATA 1, 2, 3, 4, 1, 4, 9, 16, 1, 8, 27, 64, 1, 16, 81, 256
370 DATA
380 END
500 'SON'
510 IF ABS（A（1，1））〈=E THEN GOTO 750
520 FOR J=2 TO N
530 A（1，J）=A（1，J）/A（1，1）
540 NEXT J
550 P=0
560 FOR K=2 TO N
570 FOR I=K TO N
580 FOR 1TO K-1
590 P=P+A（I，R）*A（R，K）
600 NEXT R
610 A（I，K）= A（I，K）-P
620 P=0
630 NEXT I
640 P=0
650 FOR J=K +1TO N
660 FOR R=1TO K-1
670 P=P+A（K，R）*A（R，J）
680 NEXT R
690 IF ABS（A（K，K））〈=E THEN GOTO 750
700 A（K，J）=（A（K，J）-P）/A（K，K）
710 P=0
720 NEXT J
730 NEXT K
740 W=0：RETURN
750 W=1：RETURN
====================================================================
要求用C语言或者pascal做出来，请大家务必帮忙。

[guowuji]

已经有人作出一部分了，请大家把后面的继续完成即可，谢谢了
program JuZhen;

{$APPTYPE CONSOLE}

uses
SysUtils;

const
//360 DATA 1, 2, 3, 4, 1, 4, 9, 16, 1, 8, 27, 64, 1, 16, 81, 256
//370 DATA
OriData:array[0..15] of Integer=(1, 2, 3, 4, 1, 4, 9, 16, 1, 8, 27, 64, 1, 16, 81, 256);

function Tab(N:Integer):String;
var
i:Integer;
begin
Result:='';
for i:=1 to N do
Result:=Result+' ';
end;

var
N,E,W:Integer;
A:array of array of Integer;
i,j:Integer;
begin
Writeln('N,E='); //40 INPUT "N,E="; N, E
Readln(N);
Readln(E);
SetLength(A,N); //50 DIM A(N, N)
for i:=0 to N-1 do
SetLength(A[i],N);
//60 PRINT TAB(3); "EXAMPLE"; TAB(8); "JUZHEN A": PRINT
Writeln(Tab(3)+'EXAMPLE'+Tab(8)+'JUZHEN A');
for i:=0 to N-1 do //70 FOR I = 1 TO N
begin
for j:=0 to N-1 do //80 FOR J = 1 TO N
begin //90 READ A(I, J): PRINT USING; "####.#"; A(I, J);
A[i,j]:=OriData[i*N+j];
Write(Format('%4.1f',[A[i,j]*1.0]));
end; //100 NEXT J
Writeln; //110 PRINT
end; //120 NEXT I

//130 PRINT: GOSUB 500
//GOSUB 500 ?? ******
//140 PRINT TAB(5); "W="; W: PRINT
Writeln(Tab(5)+'W='+IntToStr(W));
if W=1 then //150 IF (W = 1) THEN GOTO 380
exit;
//160 PRINT TAB(15); "JUZHEN L"
//170 PRINT
Writeln(Tab(15)+'JUZHEN L');
for i:=0 to N-1 do //180 FOR I = 1 TO N
begin
for j:=0 to i do //190 FOR J = 1 TO I
//200 PRINT USING; "##.#####"; A(I, J);: PRINT " ";
Write(Format('%2.5f ',[A[i,j]*1.0]));
//210 NEXT J
Writeln; //220 PRINT
end; //230 NEXT I
for i:=0 to N-1 do //240 FOR I = 1 TO N
A[i,i]:=1; //250 A(I, I) = 1
//260 NEXT I
Writeln; //270 PRINT
//280 PRINT TAB(15); "JUZHEN U"
//290 PRINT
Writeln(' JUZHEN U');
for i:=0 to N-1 do //300 FOR I = 1 TO N
begin
for j:=i to N-1 do //310 FOR J = I TO N
//320 PRINT TAB(10 * (J - 1)); USING; "##.#####"; A(I, J);
Writeln(Tab(10*j)+Format('%2.5f',[A[i,j]*1.0]));
//330 NEXT J
Writeln; //340 PRINT
end; //350 NEXT I

end. //380 END

[guowuji]

就是把那段Basic翻译成C语言或者pascal就行。

[fei19790920]

没上过大学

[qiluping]

up

[shifan]

//LU分解
//A为n x n,index为n+1个元素的整形数组
//结束时A为复合的LU分解阵，即有上三角阵的全部和下三角阵的除去主对角线部分，下三
//角的主对角线为1，index返回行交换的信息，d为正负1，表示行交换为奇数还是偶数，
//用于求行列式时确定符号
//部分选主元，最后要用index把行交换回来
procedure LU_Decompose(var A:TMatrix;var index:IVector;var d:Extended);
var
i,imax,j,k,n:Integer;
big,dum,sum,temp:Extended;
vv:array of Extended;
begin
imax:=0;
if A.nRow<>A.nCol then
raise Exception.Create('LU_Decompose:Invalid Matrix A');
n:=A.nRow;
SetLength(vv,n+1);
SetLength(index,n+1);
//for i:=1 to n do index[i]:=i;
d:=1;
for i:=1 to n do
begin
big:=0;
for j:=1 to n do
begin
temp:=abs(A[i,j]);
if temp>big then
big:=temp;
end;
if big=0 then
raise Exception.Create('LU_Decompose: Singular Matrix');
vv[i]:=1.0/big;
end;
for j:=1 to n do
begin
for i:=1 to j-1 do
begin
sum:=A[i,j];
for k:=1 to i-1 do
sum:=sum-A[i,k]*A[k,j];
a[i,j]:=sum;
end;
big:=0;
for i:=j to n do
begin
sum:=a[i,j];
for k:=1 to j-1 do
sum:=sum-A[i,k]*A[k,j];
a[i,j]:=sum;
dum:=vv[i]*abs(sum);
if dum>=big then
begin
big:=dum;
imax:=i;
end;
end;
if j<>imax then
begin
for k:=1 to n do
begin
dum:=a[imax,k];
A[imax,k]:=A[j,k];
A[j,k]:=dum;
end;
d:=-d;
vv[imax]:=vv[j];
end;
index[j]:=imax;
if A[j,j]=0 then
A[j,j]:=TINY;
if j<>n then
begin
dum:=1/A[j,j];
for i:=j+1 to n do
A[i,j]:=A[i,j]*dum;
end;
end;
end;

[shifan]

补充一下，TMatrix是一个1基准的矩阵类，换成2维数组也可以
nRow，nCol返回TMatrix的行数和列数
type
IVector=array of Integer;
const
TINY=1e-20;

[shifan]

不选主元的话。。
//A为输入矩阵，A_L返回下三角，A_U返回上三角
//因为不选主元，所以有的矩阵会失败（抛出DivideByZero的异常）
procedure LU_Decompose(const A:TMatrix;var A_L,A_U:TMatrix);
var
n,i,j,k,l:Integer;
sum:Extended;
begin
sum:=0;
if A.nRow<>A.nCol then
raise Exception.Create('LU_Decompose:Invalid Matrix A');
n:=A.nRow;
A_L.nRow:=n;
A_L.nCol:=n;
A_U.nRow:=n;
A_U.nCol:=n;
for i:=1 to n do
begin
A_L[i,i]:=1;
for j:=i+1 to n do
A_L[i,j]:=0;
for j:=1 to i-1 do
A_U[i,j]:=0;
end;
for j:=1 to n do
begin
for i:=1 to j do
begin
sum:=0;
if i<>1 then
for k:=1 to i-1 do
sum:=sum+A_L[i,k]*A_U[k,j];
A_U[i,j]:=A[i,j]-sum;
end;
for i:=j+1 to n do
begin
sum:=0;
for k:=1 to i-1 do
sum:=sum+A_L[i,k]*A_U[k,j];
A_L[i,j]:=(A[i,j]-sum)/A_U[j,j];
end;
end;
end;

[WuChenCan]

好乱啊！