//		-*- Mode: Java -*-
// Author(s):	J. Robert Buchanan and Zhoude Shao
// Address:	Department of Mathematics
//		Millersville University
//      	P.O. Box 1002
//		Millersville, PA 17551-0302
// Phone:	717-871-7305
// FAX:		717-871-7948
// Email:	Robert.Buchanan@millersville.edu, Zhoude.Shao@millersville.edu
//
// Date:	11/03/2016
//
// Purpose:	Gauss-Seidel method for approximating solution to Laplace's
//		equation on a rectangle. This program accompanies
//		Exercise 14.6.10.
//
// To compile:	javac exercise20.java
// To execute:	java exercise20
//
import java.io.*;

public class exercise20
{
    public static boolean stopcond(
				   int k	,
				   int nelts	,
				   double [][] soln
				   )
    {
	double [] norm = new double[2]	;
	int i	;

	norm[0] = 0.0	;
	norm[1] = 0.0	;
	for ( i=0 ; i<nelts ; i++ ) {
	    norm[0] += (soln[0][i]-soln[1][i])*(soln[0][i]-soln[1][i])	;
	    norm[1] += soln[k%2][i]*soln[k%2][i]	;
	}	// end of for ( i ) loop
	norm[0] = Math.sqrt(norm[0])	;
	norm[1] = Math.sqrt(norm[1])	;

	return( norm[0]/norm[1] < 1.0E-06 )	;
    }	// end of function stopcond( )
    
    public static void main(String args[])
    {
	int i, j, k;	// counter(s)
	int m = 10, n = 10	;
	int nelts = m*n;	// dimension of linear system
	int maxiters = 2000+1;	// hard stop at 2000 iterations of method
	double [][] a = new double[nelts][nelts];	// matrix coefficients
	double [] b = {	// RHS of equation
	    -4.91921, -2.22554, -2.01375, -1.82212, -1.64872, -1.49182,
	    -1.34986, -1.2214, -1.10517, -1.0, -2.22554, 0.0, 0.0, 0.0,
	    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -2.01375, 0.0, 0.0, 0.0, 0.0,
	    0.0, 0.0, 0.0, 0.0, 0.0, -1.82212, 0.0, 0.0, 0.0, 0.0, 0.0,
	    0.0, 0.0, 0.0, 0.0, -1.64872, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
	    0.0, 0.0, 0.0, -1.49182, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
	    0.0, 0.0, -1.34986, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
	    0.0, -1.2214, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
	    -1.10517, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0,
	    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0	}	;
	double [][] soln = new double[2][nelts];// 2 recent approximations

	// Populate matrix a, since a is sparse start by zero-ing out
	for ( i=0 ; i<nelts ; i++ ) for ( j=0 ; j<nelts ; j++ ) a[i][j] = 0.0;
	// Set the nonzero entries
	for ( i=0 ; i<n ; i++ ) {	// first block row
	    a[i][i] = -4.0	;	// diagonal entries
	    a[i][i+n] = 1.0	;	// nth super diagonal entries
	    if ( i==0 ) a[i][i+1] = 1.0;
	    else if ( i==(n-1) ) a[i][i-1] = 2.0;
	    else {
		a[i][i-1] = 1.0;
		a[i][i+1] = 1.0;
	    }
	}	// end for ( i ) loop
	for ( i=n ; i<nelts-n ; i++ ) {	// bulk m-2 blocks
	    a[i][i] = -4.0	;	// diagonal entries
	    a[i][i+n] = 1.0	;
	    a[i][i-n] = 1.0	;
	    if ( 0==(i%n) ) a[i][i+1] = 1.0	;
	    else if ( (n-1)==(i%n) ) a[i][i-1] = 2.0	;
	    else {
		a[i][i-1] = 1.0	;
		a[i][i+1] = 1.0	;
	    }
	}	// end of for ( i ) loop
	for ( i=nelts-n ; i<nelts ; i++ ) {	// last block row
	    a[i][i] = -4.0	;	// diagonal entries
	    a[i][i-n] = 2.0	;	// nth sub diagonal entries
	    if ( i==(nelts-n) ) a[i][i+1] = 1.0	;
	    else if ( i==nelts-1 ) a[i][i-1] = 2.0	;
	    else {
		a[i][i-1] = 1.0;
		a[i][i+1] = 1.0;
	    }
	}	// end for ( i ) loop
	// Initial approximation
	for ( i=0 ; i<nelts ; i++ ) soln[0][i] = 0.0	;
	for ( k=1 ; k<maxiters ; k++ ) {
	    for ( i=0 ; i<nelts ; i++ ) {
		soln[k%2][i] = b[i]	;
		for ( j=0 ; j<i ; j++ ) soln[k%2][i] -= a[i][j]*soln[k%2][j];
		for ( j=i ; j<nelts ; j++ ) soln[k%2][i] -= a[i][j]*soln[(k-1)%2][j];
		soln[k%2][i] /= a[i][i]	;
		soln[k%2][i] += soln[(k-1)%2][i];
	    }	// end of for ( i ) loop
	    if ( stopcond(k,nelts,soln) ) break	;
	}	// end of for ( k ) loop
	System.out.printf("Gauss-Seidel method halted on iteration %d.\n",k);
	for ( i=0 ; i<nelts ; i++ ) System.out.printf("%8.4f ", soln[k%2][i]);
	System.out.printf("\n")	;
    }	// end of function main( )
}

//
//	EOF
//