//		-*- 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:	10/12/2016
//
// Purpose:	Gauss-Seidel method for approximating solution to Laplace's
//		equation on a square. This program generates data for
//		Figure 14.10.
//
// To compile:	javac exampleLaplace.java
// To execute:	java exampleJLaplace
//
import java.io.*;

public class exampleLaplace
{
    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 nelts = 9*9	;	// dimension of linear system
	int maxiters = 200+1;	// hard stop at 200 iterations of method
	double [][] a = new double[nelts][nelts];	// matrix coefficients
	double [] b = {	// RHS of equation
	    2.09517, 1.2214, 1.34986, 1.49182, 1.64872, 1.82212, 2.01375,
	    2.22554, 4.91921, 0.96, 0., 0., 0., 0., 0., 0., 0., 2.22554, 0.91,
	    0., 0., 0., 0., 0., 0., 0., 2.01375, 0.84, 0., 0., 0., 0., 0., 0.,
	    0., 1.82212, 0.75, 0., 0., 0., 0., 0., 0., 0., 1.64872, 0.64, 0.,
	    0., 0., 0., 0., 0., 0., 1.49182, 0.51, 0., 0., 0., 0., 0., 0., 0.,
	    1.34986, 0.36, 0., 0., 0., 0., 0., 0., 0., 1.2214, 0.346434,
	    0.309017, 0.45399, 0.587785, 0.707107, 0.809017, 0.891007,
	    0.951057, 2.09286
	}	;
	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<nelts ; i++ ) {
	    a[i][i] = 4.0	;	// diagonal entries
	    if ( i>(9-1) ) a[i][i-9] = -1.0	;	// 10th sub-diagonal
	    if ( i<9*(9-1) ) a[i][i+9] = -1.0	;	// 10th sup-diagonal
	    // every 9th row starting with 1st has pattern 0 ... 0 -4 1 0 ...
	    if ( (i%9)==0 ) a[i][i+1] = -1.0	;
	    // every 9th row starting with 8th has pattern 0 ... 0 1 -4 0 ...
	    else if ( (i%9)==(9-1) ) a[i][i-1] = -1.0;
	    // all other rows have pattern 0 ... 0 1 -4 1 0 ...
	    else {
		a[i][i-1] = -1.0	;
		a[i][i+1] = -1.0	;
	    }	// end of else clause
	}	// end of for ( i ) loop
	// Initial approximation
	for ( i=0 ; i<nelts ; i++ ) soln[0][i] = 1.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+1 ; j<nelts ; j++ ) soln[k%2][i] -= a[i][j]*soln[(k-1)%2][j];
		soln[k%2][i] /= a[i][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
//