/* This class implements the 1D Poisson Matrix operators, i.e., 
 * matrix-vector multiply and psolve.
 */

public class Poisson1D implements Sparse {
    
    public Poisson1D (int n) {
	this.n = n;
    } 

    public String toString () {
	String result = "";

	for (int i = 0; i < n; i++) {
	    for (int j = 0; j < n; j++) {
		if (i == j) {
		    result += "2\t";
		}
		else if ((j == i-1) || (j == i+1)) {
		    result += "-1\t";
		}
		else {
		    result += "0\t";
		}
	    }
	    result += "\n";
	}
	return result;
    }
   
    public double [2d] toDense () {
	int high = n-1;
	double [2d] m = new double [0:high,0:high];
	m[0,0] = 2;
	m[0,1] = -1;
	foreach (i in [1:high-1]) {
	    m[i[1],i[1]] = 2;
	    m[i[1],i[1]-1] = -1;
	    m[i[1],i[1]+1] = -1;
	}
	m[high,high-1] = -1;
	m[high,high] = 2;
	return m;
    }

    public void matvec (double [1d] single [1d] y, double [1d] single [1d] x) {
	/* for simplicity, may want to assume at least 2 elems per proc */
	/* not implemented -- good warmup exercise */
	y[0][0] = 2*x[0][0] - x[0][1];
	foreach (p in y.domain().shrink(1,-1)) {
	    y[p][0] = 2*x[p][0] - x[p][1] - x[p-1][x[p-1].domain().max()];
	}
	foreach (p in y.domain()) {
	    foreach (i in y.domain().shrink(1)) {
		y[i] = -x[i+[-1]] + 2*x[i] - x[i+[1]];
	}
	Point<1> high = x.domain().max();
	y[high] = 2*x[high] - x[high+[-1]];
    } 

    public int dim () {
	return n;
    }

    /* Internal State */
    protected int n;  // n should be used for error checking in matvec, etc.

    /* Test code */
    public static void tester (String [] args) {
	if (args.length < 1) {
	    printUsage();
	    return;
	}
	int sz = Integer.parseInt(args[0]);

	// Allocate source and destination vectors
	double [1d] vec1 = new double [0:sz-1];
	double [1d] vec2 = new double [0:sz-1];
	foreach (i in vec1.domain()) {
	    vec1[i] = i[1];
	}

	// Test: Poisson matrix times vector
	Poisson1D p = new Poisson1D(sz);
       	System.out.println("Poisson matrix: \n" + p);
	p.matvec(vec2, vec1);
	System.out.println("Source vector v1: \n" + Vector.toString(vec1) 
			   + "\n");
	System.out.println("Result vector v2: \n" + Vector.toString(vec2) 
			   + "\n"); 
    }

    private static void printUsage() {
	System.out.println("Poisson1D <n> \n for a test run on an nxn matrix");
    }
}