mylu.m

CG共轭梯度法迭代解矩阵方程

function [l,u,p]=mylu(a)
siz=size(a);
if siz(1)~=siz(2)
error('矩阵非为方阵')
return
end
if det(a)==0
 error('矩阵不能被三角分解')
end
si=siz(1);
u=a;p=eye(si);l=eye(si);
for i=1:si
for j=i:si
t(j)=u(j,i);
for k=1:i-1
t(j)=t(j)-u(j,k)*u(k,i);
end
end
 pos=i;ma=abs(t(i));
for j=i+1:si
if ma<abs(t(j))
ma=abs(t(j));
 pos=j;
end
end
if pos~=i
for j=1:si
temp=u(i,j);
u(i,j)=u(pos,j);
 u(pos,j)=temp;
end
for j=1:si
temp=p(i,j);
           p(i,j)=p(pos,j);
            p(pos,j)=temp;
         end
       temp=t(pos);
       t(pos)=t(i);
       t(i)=temp;
    end
    u(i,i)=t(i);
    for j=i+1:si
        u(j,i)=t(j)/t(i);
    end
    for j=i+1:si
        for k=1:i-1
            u(i,j)=u(i,j)-u(i,k)*u(k,j);
        end
    end
end
l=tril(u,-1)+eye(si);
u=triu(u,0);