iccg.m

来自「iccg不完全乔尔斯基分解共轭梯度求解大规模稀疏方程组」· M 代码 · 共 119 行

function bb=iccg(aa,ss,t,nd,nc)
nnd=6;
nnc=10;
lll=1;
ll=1;
for ii=1:nnd
    bb(ii)=t(ii);
end
sina(nnd,aa,bb,nd,nc)
ep0=0;
for ii=1:nnd
    c=abs(t(ii));
    ep0=ep0+c;
end
ep0
eps=ep0*1.0e-5;
%10
v=smi(ss,bb,nd,nc)
    for ii=1:nnd
        r(ii)=t(ii)-v(ii);
        p(ii)=r(ii);
    end
    sina(nnd,aa,p,nd,nc)
    tf=0;
    for ii=1:nnd
        tf=tf+r(ii)*p(ii);
    end
    for jj=1:nnd
        tr=0;
        ep=0;
        tn=0;
        v=smi(ss,p,nd,nc);
        for ii=1:nnd
            tr=tr+p(ii)*v(ii);
        end
        arf=tf/tr;
        for ii=1:nnd
            bb(ii)=bb(ii)+arf*p(ii);
            r(ii)=r(ii)-arf*v(ii);
            v(ii)=r(ii);
            c=abs(v(ii));
            ep=ep+c;
        end
        if(lll==1|mod(lll,1)==0)
               lll
               ep
        end
        if(ep<eps)
            break
            lll
            ep
        end
        sina(nnd,aa,v,nd,nc);
        for ii=1:nnd
            tn=tn+r(ii)*v(ii);
        end
        bat=tn/tf;
        tf=tn;
        for ii=1:nnd
            p(ii)=bat*p(ii)+v(ii);
        end
        lll=lll+1;
    end
    ll=ll+1;
    while(ep>eps)
        smi(ss,v,vv,nd,nc)
        for ii=1:nnd
            r(ii)=t(ii)-v(ii);
            p(ii)=r(ii);
        end
        sina(nnd,aa,p,nd,nc)
        tf=0;
        for ii=1:nnd
            tf=tf+r(ii)*p(ii);
        end
        for jj=1:nnd
            tr=0;
            ep=0;
            tn=0;
            smi(ss,v,p,nd,nc)
            for ii=1:nnd
                tr=tr+p(ii)*v(ii);
            end
            arf=tf/tr;
            for ii=1:nnd
                bb(ii)=bb(ii)+arf*p(ii);
                r(ii)=r(ii)-arf*v(ii);
                v(ii)=r(ii);
                c=abs(v(ii));
                ep=ep+c;
            end
            if(lll==1|mod(lll,1)==0)
                lll
                ep
            end
            if(ep<eps)
                break;
            end
            sina(nnd,aa,v,nd,nc);
            for ii=1:nnd
                tn=tn+r(ii)*v(ii);
            end
            bat=tn/tf;
            tf=tn;
            for ii=1:nnd
                p(ii)=bat*p(ii)+v(ii);
            end
            lll=lll+1;
        end
        ll=ll+1;
    end
   lll
    ep