linbcg.java

PEPA模型性能分析工具

package matrix;

import gui.PepaConsole;

public class Linbcg {

    private static double EPS = 1.0E-15;
	private int[] ija;
	private double[] sa;
	private int size;
	private int vol;
	
	public Linbcg(int[] i,double[] d,int m)
	{
		ija=i;
		sa=d;
		vol=m+1;
		size=ija.length;
	}
	
	public void solve(int n, double[] b, double[] x, int itol, double tol, int itmax,
			int iter, double err, boolean gmres)
	{	double[] p,pp,r,rr,z,zz;
	   	double ak,akden,bk,bkden,bknum,bnrm,dxnrm,xnrm,zm1nrm,znrm;
	   	bkden=1.0;
	   	
	   	p = new double[vol];
        pp = new double[vol];
        r = new double[vol];
        rr = new double[vol];
        z = new double[vol];
        zz = new double[vol];
        
        for(int i= 1;i<=n;i++)
            x[i] = 0.0;  
       
        atimes(n, x, r, 0);
        for (int j = 1; j <= n; j ++) 
        {
            r[j] = b[j] - r[j];
            rr[j] = r[j];
        }        
        znrm = 1.0; //????
        // We are using the GMRES
        // generalised minimum residual method
        if(gmres)
        {
            atimes(n,r,rr,0);
        }
        
        if (itol == 1)
        {
            bnrm = snrm(n, b, itol);
            asolve(n, r, z, 0);
        }
        else if (itol == 2) 
        {
            asolve(n, b, z, 0);
            bnrm = snrm(n, z, itol);
            asolve(n, r, z, 0);
        }
        else if (itol == 3 || itol == 4) 
        {
            asolve(n, b, z, 0);
            bnrm = snrm(n, z, itol);
            asolve(n, r, z, 0);
            znrm = snrm(n, z, itol);
        }
        else throw new Error("Illegal itol Value");
        
        // MAIN WHILE lOOP
        while (iter <= itmax) 
        {            
            ++ iter;
            zm1nrm = znrm;
            asolve(n, rr, zz, 1);
            // for(bknum=0.0,j=1;j<=n;j++) bknum += z[j]*rr[j];
            bknum = 0.0;
            for (int j = 1; j <= n; j ++)
            {
                bknum += z[j] * rr[j];
            }
            // CALCULATE COEFFICIENT bk AND DIRECTION VECTORS p AND pp
            if (iter == 1) 
            {
                for (int j = 1; j <= n; j ++)
                {
                    p[j] = z[j];
                    pp[j] = zz[j];
                }
            }
            else
            {                
                bk = bknum / bkden;
                for (int j = 1; j <= n; j ++)
                {
                    p[j] = bk * p[j] + z[j];
                    pp[j] = bk * pp[j] + zz[j];
                }
            }            
            bkden = bknum;
            
            // CALCULATE COEFFICIENT ak, NEW ITERATE x and NEW RESIDUALS
            atimes(n, p, z, 0);
            
            // for (akden = 0.0,j=1;j<=n;j++) akden+= z[j]*pp[j];
            akden = 0.0;
            for (int j = 1; j <= n; j ++) akden += z[j] * pp[j];
            ak = bknum / akden;
            
            atimes(n, pp, zz, 1);
            
            for (int j = 1; j <= n; j ++) 
            {
                x[j] += ak * p[j];
                r[j] -= ak * z[j];
                rr[j] -= ak * zz[j];
            }
            asolve(n, r, z, 0);
            // Solve A.Z = r and check stopping criterion.
            if (itol == 1) 
            {
                znrm = 1.0;
                err =  snrm(n, r, itol)/ bnrm;
            }
            else if (itol == 2)
            {
                err = snrm(n, z, itol) / bnrm;
            }
            else if (itol == 3 || itol == 4) 
            {
                zm1nrm = znrm;
                znrm = snrm(n, z, itol);
                if (Math.abs(zm1nrm - znrm) > EPS * znrm) 
                {
                    dxnrm = Math.abs(ak) * snrm(n, p, itol);
                    err = znrm / Math.abs(zm1nrm - znrm) * dxnrm;
                }
                else 
                {
                    err = znrm / bnrm;
                    continue;
                }
                xnrm = snrm(n, x, itol);
                
                if (err <= 0.5 * xnrm)
                    err /= xnrm;
                else
                {
                    err = znrm / bnrm;
                    continue;
                }
            }       
            if (err <= tol) break; 
        }
        PepaConsole.addArgument("--迭代次数："+iter+"\n");
        if(iter>=itmax)
        	PepaConsole.addArgument("--已经迭代到最大次数，得到的结果将不够精确，请增加迭代次数!\n");
        PepaConsole.addArgument("--误差："+err+"\n");
        p = null;
        pp = null;
        z = null;
        zz = null;
        r = null;
        rr = null; 
        
        System.gc();
	}
	  private void atimes(int n, double[] x, double[] r, int itrnsp)
	  {
	        if (itrnsp == 1) 
	        	dsprstx(x, r, n); 
	        else dsprsax(x, r, n);
	    }
	    
	    /**
	     * solves ~Ax=b or ~A^T x=b for some preconditioner matrix ~A (possibly
	     * the trivial diagonal part of A)
	     */
	    private void asolve(int n, double[] b, double[] x, int useless)
	    {
	        for (int i = 1; i <= n; i ++) 
	        	x[i] = sa[i] != 0.0 ? b[i] / sa[i] : b[i];
	    }
	    
	    private void dsprsax(double[] x, double[] b, int n) {
	        int i;
	        int k;
	        if (ija[1] != n + 2) 
	        {
	            System.exit(- 1);
	        }
	        
	        for (i = 1; i <= n; i ++) 
	        {
	            b[i] = sa[i] * x[i];
	            for (k = ija[i]; k <= ija[i + 1] - 1; k ++)
	            {
	            	int ttt=ija[k];
	                b[i] += sa[k] * x[ija[k]];
	            }
	        }
	    }
	    
	    private void dsprstx(double[] x, double[] b, int n) 
	    {
	        int i,j,k;
	        if (ija[1] != n + 2) 
	        {
	            System.exit(- 1);
	        }
	        for (i = 1; i <= n; i ++) 
	        	b[i] = sa[i] * x[i];
	        for (i = 1; i <= n; i ++) 
	        {
	            for (k = ija[i]; k <= ija[i + 1] - 1; k ++)
	            {
	                j = ija[k];
	                b[j] += sa[k] * x[i];
	            }
	        }
	    }
	    
	    
	    /**
	     * Compute one of two norms for a vector sx[1..n], as signaled by itol
	     */
	    private double snrm(int n, double[] sx, int itol) 
	    {       
	        int i,isamax;
	        double ans;
	        
	        if (itol <= 3)
	        {
	            ans = 0.0;
	            for (i = 1; i <= n; i ++)
	            {                
	                ans += sx[i] * sx[i];
	            }
	            return Math.sqrt(ans);
	        }
	        else 
	        {
	            isamax = 1;
	            for (i = 1; i <= n; i ++) 
	            {
	                if (Math.abs(sx[i]) > Math.abs(sx[isamax]))
	                	isamax = i;
	            }
	            return Math.abs(sx[isamax]);
	        }
	    }
}