//		-*- 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:	SOR method for approximating solution to Poisson's
//		equation on a rectangle. This program generates data for
//		Example 14.9.
//
// To compile:	javac SORPoisson.java
// To execute:	java SORPoisson
//
import java.io.*;

public class SORPoisson
{
    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 = 30, n = 40	;
	int nelts = (m-1)*(n-1);// dimension of linear system
	int maxiters = 200+1;	// hard stop at 200 iterations of method
	double w = 1.81621;	// SOR parameter
	double [][] a = new double[nelts][nelts];	// matrix coefficients
	double [] b = {	// RHS of equation
	    -0.0601, -0.0664208, -0.0734063, -0.0811265, -0.0896587, -0.0990881,
	    -0.109509, -0.121027, -0.133755, -0.147822, -0.163369, -0.18055,
	    -0.199539, -0.220525, -0.243718, -0.26935, -0.297677, -0.328984,
	    -0.363584, -0.401822, -0.444082, -0.490787, -0.542403, -0.599448,
	    -0.662493, -0.732168, -0.809171, -0.894272, -0.988323, -1.09227,
	    -1.20714, -1.3341, -1.47441, -1.62947, -1.80084, -1.99024, -2.19955,
	    -2.43088, -0.286541, -0.0543807, -0.0601, -0.0664208, -0.0734063,
	    -0.0811265, -0.0896587, -0.0990881, -0.109509, -0.121027, -0.133755,
	    -0.147822, -0.163369, -0.18055, -0.199539, -0.220525, -0.243718,
	    -0.26935, -0.297677, -0.328984, -0.363584, -0.401822, -0.444082,
	    -0.490787, -0.542403, -0.599448, -0.662493, -0.732168, -0.809171,
	    -0.894272, -0.988323, -1.09227, -1.20714, -1.3341, -1.47441,
	    -1.62947, -1.80084, -1.99024, -2.19955, 2.36912, -0.0492057,
	    -0.0543807, -0.0601, -0.0664208, -0.0734063, -0.0811265, -0.0896587,
	    -0.0990881, -0.109509, -0.121027, -0.133755, -0.147822, -0.163369,
	    -0.18055, -0.199539, -0.220525, -0.243718, -0.26935, -0.297677,
	    -0.328984, -0.363584, -0.401822, -0.444082, -0.490787, -0.542403,
	    -0.599448, -0.662493, -0.732168, -0.809171, -0.894272, -0.988323,
	    -1.09227, -1.20714, -1.3341, -1.47441, -1.62947, -1.80084, -1.99024,
	    5.00045, -0.0445232, -0.0492057, -0.0543807, -0.0601, -0.0664208,
	    -0.0734063, -0.0811265, -0.0896587, -0.0990881, -0.109509,-0.121027,
	    -0.133755, -0.147822, -0.163369, -0.18055, -0.199539, -0.220525,
	    -0.243718, -0.26935, -0.297677, -0.328984, -0.363584, -0.401822,
	    -0.444082, -0.490787, -0.542403, -0.599448, -0.662493, -0.732168,
	    -0.809171, -0.894272, -0.988323, -1.09227, -1.20714, -1.3341,
	    -1.47441, -1.62947, -1.80084, 7.60976, -0.0402862, -0.0445232,
	    -0.0492057, -0.0543807, -0.0601, -0.0664208, -0.0734063, -0.0811265,
	    -0.0896587, -0.0990881, -0.109509, -0.121027, -0.133755, -0.147822,
	    -0.163369, -0.18055, -0.199539, -0.220525, -0.243718, -0.26935,
	    -0.297677, -0.328984, -0.363584, -0.401822, -0.444082, -0.490787,
	    -0.542403, -0.599448, -0.662493, -0.732168, -0.809171, -0.894272,
	    -0.988323, -1.09227, -1.20714, -1.3341, -1.47441, -1.62947, 10.1992,
	    -0.0364525, -0.0402862, -0.0445232, -0.0492057, -0.0543807, -0.0601,
	    -0.0664208, -0.0734063, -0.0811265, -0.0896587, -0.0990881,
	    -0.109509, -0.121027, -0.133755, -0.147822, -0.163369, -0.18055,
	    -0.199539, -0.220525, -0.243718, -0.26935, -0.297677, -0.328984,
	    -0.363584, -0.401822, -0.444082, -0.490787, -0.542403, -0.599448,
	    -0.662493, -0.732168, -0.809171, -0.894272, -0.988323, -1.09227,
	    -1.20714, -1.3341, -1.47441, 12.7705, -0.0329836, -0.0364525,
	    -0.0402862, -0.0445232, -0.0492057, -0.0543807, -0.0601, -0.0664208,
	    -0.0734063, -0.0811265, -0.0896587, -0.0990881, -0.109509,-0.121027,
	    -0.133755, -0.147822, -0.163369, -0.18055, -0.199539, -0.220525,
	    -0.243718, -0.26935, -0.297677, -0.328984, -0.363584, -0.401822,
	    -0.444082, -0.490787, -0.542403, -0.599448, -0.662493, -0.732168,
	    -0.809171, -0.894272, -0.988323, -1.09227, -1.20714, -1.3341,
	    15.3256, -0.0298448, -0.0329836, -0.0364525, -0.0402862, -0.0445232,
	    -0.0492057, -0.0543807, -0.0601, -0.0664208, -0.0734063, -0.0811265,
	    -0.0896587, -0.0990881, -0.109509, -0.121027, -0.133755, -0.147822,
	    -0.163369, -0.18055, -0.199539, -0.220525, -0.243718, -0.26935,
	    -0.297677, -0.328984, -0.363584, -0.401822, -0.444082, -0.490787,
	    -0.542403, -0.599448, -0.662493, -0.732168, -0.809171, -0.894272,
	    -0.988323, -1.09227, -1.20714, 17.8659, -0.0270047, -0.0298448,
	    -0.0329836, -0.0364525, -0.0402862, -0.0445232, -0.0492057,
	    -0.0543807, -0.0601, -0.0664208, -0.0734063, -0.0811265, -0.0896587,
	    -0.0990881, -0.109509, -0.121027, -0.133755, -0.147822, -0.163369,
	    -0.18055, -0.199539, -0.220525, -0.243718, -0.26935, -0.297677,
	    -0.328984, -0.363584, -0.401822, -0.444082, -0.490787, -0.542403,
	    -0.599448, -0.662493, -0.732168, -0.809171, -0.894272, -0.988323,
	    -1.09227, 20.3929, -0.0244348, -0.0270047, -0.0298448, -0.0329836,
	    -0.0364525, -0.0402862, -0.0445232, -0.0492057, -0.0543807, -0.0601,
	    -0.0664208, -0.0734063, -0.0811265, -0.0896587, -0.0990881,
	    -0.109509, -0.121027, -0.133755, -0.147822, -0.163369, -0.18055,
	    -0.199539, -0.220525, -0.243718, -0.26935, -0.297677, -0.328984,
	    -0.363584, -0.401822, -0.444082, -0.490787, -0.542403, -0.599448,
	    -0.662493, -0.732168, -0.809171, -0.894272, -0.988323, 22.9077,
	    -0.0221096, -0.0244348, -0.0270047, -0.0298448, -0.0329836,
	    -0.0364525, -0.0402862, -0.0445232, -0.0492057, -0.0543807, -0.0601,
	    -0.0664208, -0.0734063, -0.0811265, -0.0896587, -0.0990881,
	    -0.109509, -0.121027, -0.133755, -0.147822, -0.163369, -0.18055,
	    -0.199539, -0.220525, -0.243718, -0.26935, -0.297677, -0.328984,
	    -0.363584, -0.401822, -0.444082, -0.490787, -0.542403, -0.599448,
	    -0.662493, -0.732168, -0.809171, -0.894272, 25.4117, -0.0200056,
	    -0.0221096, -0.0244348, -0.0270047, -0.0298448, -0.0329836,
	    -0.0364525, -0.0402862, -0.0445232, -0.0492057, -0.0543807, -0.0601,
	    -0.0664208, -0.0734063, -0.0811265, -0.0896587, -0.0990881,
	    -0.109509, -0.121027, -0.133755, -0.147822, -0.163369, -0.18055,
	    -0.199539, -0.220525, -0.243718, -0.26935, -0.297677, -0.328984,
	    -0.363584, -0.401822, -0.444082, -0.490787, -0.542403, -0.599448,
	    -0.662493, -0.732168, -0.809171, 27.9057, -0.0181018, -0.0200056,
	    -0.0221096, -0.0244348, -0.0270047, -0.0298448, -0.0329836,
	    -0.0364525, -0.0402862, -0.0445232, -0.0492057, -0.0543807, -0.0601,
	    -0.0664208, -0.0734063, -0.0811265, -0.0896587, -0.0990881,
	    -0.109509, -0.121027, -0.133755, -0.147822, -0.163369, -0.18055,
	    -0.199539, -0.220525, -0.243718, -0.26935, -0.297677, -0.328984,
	    -0.363584, -0.401822, -0.444082, -0.490787, -0.542403, -0.599448,
	    -0.662493, -0.732168, 30.3908, -0.0163792, -0.0181018, -0.0200056,
	    -0.0221096, -0.0244348, -0.0270047, -0.0298448, -0.0329836,
	    -0.0364525, -0.0402862, -0.0445232, -0.0492057, -0.0543807, -0.0601,
	    -0.0664208, -0.0734063, -0.0811265, -0.0896587, -0.0990881,
	    -0.109509, -0.121027, -0.133755, -0.147822, -0.163369, -0.18055,
	    -0.199539, -0.220525, -0.243718, -0.26935, -0.297677, -0.328984,
	    -0.363584, -0.401822, -0.444082, -0.490787, -0.542403, -0.599448,
	    -0.662493, 32.8678, -0.0148205, -0.0163792, -0.0181018, -0.0200056,
	    -0.0221096, -0.0244348, -0.0270047, -0.0298448, -0.0329836,
	    -0.0364525, -0.0402862, -0.0445232, -0.0492057, -0.0543807, -0.0601,
	    -0.0664208, -0.0734063, -0.0811265, -0.0896587, -0.0990881,
	    -0.109509, -0.121027, -0.133755, -0.147822, -0.163369, -0.18055,
	    -0.199539, -0.220525, -0.243718, -0.26935, -0.297677, -0.328984,
	    -0.363584, -0.401822, -0.444082, -0.490787, -0.542403, -0.599448,
	    35.3375, -0.0134101, -0.0148205, -0.0163792, -0.0181018, -0.0200056,
	    -0.0221096, -0.0244348, -0.0270047, -0.0298448, -0.0329836,
	    -0.0364525, -0.0402862, -0.0445232, -0.0492057, -0.0543807, -0.0601,
	    -0.0664208, -0.0734063, -0.0811265, -0.0896587, -0.0990881,
	    -0.109509, -0.121027, -0.133755, -0.147822, -0.163369, -0.18055,
	    -0.199539, -0.220525, -0.243718, -0.26935, -0.297677, -0.328984,
	    -0.363584, -0.401822, -0.444082, -0.490787, -0.542403, 37.8006,
	    -0.012134, -0.0134101, -0.0148205, -0.0163792, -0.0181018,
	    -0.0200056, -0.0221096, -0.0244348, -0.0270047, -0.0298448,
	    -0.0329836, -0.0364525, -0.0402862, -0.0445232, -0.0492057,
	    -0.0543807, -0.0601, -0.0664208, -0.0734063, -0.0811265, -0.0896587,
	    -0.0990881, -0.109509, -0.121027, -0.133755, -0.147822, -0.163369,
	    -0.18055, -0.199539, -0.220525, -0.243718, -0.26935, -0.297677,
	    -0.328984, -0.363584, -0.401822, -0.444082, -0.490787, 40.2576,
	    -0.0109793, -0.012134, -0.0134101, -0.0148205, -0.0163792,
	    -0.0181018, -0.0200056, -0.0221096, -0.0244348, -0.0270047,
	    -0.0298448, -0.0329836, -0.0364525, -0.0402862, -0.0445232,
	    -0.0492057, -0.0543807, -0.0601, -0.0664208, -0.0734063, -0.0811265,
	    -0.0896587, -0.0990881, -0.109509, -0.121027, -0.133755, -0.147822,
	    -0.163369, -0.18055, -0.199539, -0.220525, -0.243718, -0.26935,
	    -0.297677, -0.328984, -0.363584, -0.401822, -0.444082, 42.7092,
	    -0.00993446, -0.0109793, -0.012134, -0.0134101, -0.0148205,
	    -0.0163792, -0.0181018, -0.0200056, -0.0221096, -0.0244348,
	    -0.0270047, -0.0298448, -0.0329836, -0.0364525, -0.0402862,
	    -0.0445232, -0.0492057, -0.0543807, -0.0601, -0.0664208, -0.0734063,
	    -0.0811265, -0.0896587, -0.0990881, -0.109509, -0.121027, -0.133755,
	    -0.147822, -0.163369, -0.18055, -0.199539, -0.220525, -0.243718,
	    -0.26935, -0.297677, -0.328984, -0.363584, -0.401822, 45.1559,
	    -0.00898907, -0.00993446, -0.0109793, -0.012134, -0.0134101,
	    -0.0148205, -0.0163792, -0.0181018, -0.0200056, -0.0221096,
	    -0.0244348, -0.0270047, -0.0298448, -0.0329836, -0.0364525,
	    -0.0402862, -0.0445232, -0.0492057, -0.0543807, -0.0601, -0.0664208,
	    -0.0734063, -0.0811265, -0.0896587, -0.0990881, -0.109509,-0.121027,
	    -0.133755, -0.147822, -0.163369, -0.18055, -0.199539, -0.220525,
	    -0.243718, -0.26935, -0.297677, -0.328984, -0.363584, 47.5982,
	    -0.00813365, -0.00898907, -0.00993446, -0.0109793, -0.012134,
	    -0.0134101, -0.0148205, -0.0163792, -0.0181018, -0.0200056,
	    -0.0221096, -0.0244348, -0.0270047, -0.0298448, -0.0329836,
	    -0.0364525, -0.0402862, -0.0445232, -0.0492057, -0.0543807, -0.0601,
	    -0.0664208, -0.0734063, -0.0811265, -0.0896587, -0.0990881,
	    -0.109509, -0.121027, -0.133755, -0.147822, -0.163369, -0.18055,
	    -0.199539, -0.220525, -0.243718, -0.26935, -0.297677, -0.328984,
	    50.0364, -0.00735963, -0.00813365, -0.00898907, -0.00993446,
	    -0.0109793, -0.012134, -0.0134101, -0.0148205, -0.0163792,
	    -0.0181018, -0.0200056, -0.0221096, -0.0244348, -0.0270047,
	    -0.0298448, -0.0329836, -0.0364525, -0.0402862, -0.0445232,
	    -0.0492057, -0.0543807, -0.0601, -0.0664208, -0.0734063, -0.0811265,
	    -0.0896587, -0.0990881, -0.109509, -0.121027, -0.133755, -0.147822,
	    -0.163369, -0.18055, -0.199539, -0.220525, -0.243718, -0.26935,
	    -0.297677, 52.471, -0.00665927, -0.00735963, -0.00813365,
	    -0.00898907, -0.00993446, -0.0109793, -0.012134, -0.0134101,
	    -0.0148205, -0.0163792, -0.0181018, -0.0200056, -0.0221096,
	    -0.0244348, -0.0270047, -0.0298448, -0.0329836, -0.0364525,
	    -0.0402862, -0.0445232, -0.0492057, -0.0543807, -0.0601, -0.0664208,
	    -0.0734063, -0.0811265, -0.0896587, -0.0990881, -0.109509,-0.121027,
	    -0.133755, -0.147822, -0.163369, -0.18055, -0.199539, -0.220525,
	    -0.243718, -0.26935, 54.9023, -0.00602556, -0.00665927, -0.00735963,
	    -0.00813365, -0.00898907, -0.00993446, -0.0109793, -0.012134,
	    -0.0134101, -0.0148205, -0.0163792, -0.0181018, -0.0200056,
	    -0.0221096, -0.0244348, -0.0270047, -0.0298448, -0.0329836,
	    -0.0364525, -0.0402862, -0.0445232, -0.0492057, -0.0543807, -0.0601,
	    -0.0664208, -0.0734063, -0.0811265, -0.0896587, -0.0990881,
	    -0.109509, -0.121027, -0.133755, -0.147822, -0.163369, -0.18055,
	    -0.199539, -0.220525, -0.243718, 57.3307, -0.00545215, -0.00602556,
	    -0.00665927, -0.00735963, -0.00813365, -0.00898907, -0.00993446,
	    -0.0109793, -0.012134, -0.0134101, -0.0148205, -0.0163792,
	    -0.0181018, -0.0200056, -0.0221096, -0.0244348, -0.0270047,
	    -0.0298448, -0.0329836, -0.0364525, -0.0402862, -0.0445232,
	    -0.0492057, -0.0543807, -0.0601, -0.0664208, -0.0734063, -0.0811265,
	    -0.0896587, -0.0990881, -0.109509, -0.121027, -0.133755, -0.147822,
	    -0.163369, -0.18055, -0.199539, -0.220525, 59.7563, -0.00493331,
	    -0.00545215, -0.00602556, -0.00665927, -0.00735963, -0.00813365,
	    -0.00898907, -0.00993446, -0.0109793, -0.012134, -0.0134101,
	    -0.0148205, -0.0163792, -0.0181018, -0.0200056, -0.0221096,
	    -0.0244348, -0.0270047, -0.0298448, -0.0329836, -0.0364525,
	    -0.0402862, -0.0445232, -0.0492057, -0.0543807, -0.0601, -0.0664208,
	    -0.0734063, -0.0811265, -0.0896587, -0.0990881, -0.109509,-0.121027,
	    -0.133755, -0.147822, -0.163369, -0.18055, -0.199539, 62.1795,
	    -0.00446384, -0.00493331, -0.00545215, -0.00602556, -0.00665927,
	    -0.00735963, -0.00813365, -0.00898907, -0.00993446, -0.0109793,
	    -0.012134, -0.0134101, -0.0148205, -0.0163792, -0.0181018,
	    -0.0200056, -0.0221096, -0.0244348, -0.0270047, -0.0298448,
	    -0.0329836, -0.0364525, -0.0402862, -0.0445232, -0.0492057,
	    -0.0543807, -0.0601, -0.0664208, -0.0734063, -0.0811265, -0.0896587,
	    -0.0990881, -0.109509, -0.121027, -0.133755, -0.147822, -0.163369,
	    -0.18055, 64.6005, -0.00403905, -0.00446384, -0.00493331,
	    -0.00545215, -0.00602556, -0.00665927, -0.00735963, -0.00813365,
	    -0.00898907, -0.00993446, -0.0109793, -0.012134, -0.0134101,
	    -0.0148205, -0.0163792, -0.0181018, -0.0200056, -0.0221096,
	    -0.0244348, -0.0270047, -0.0298448, -0.0329836, -0.0364525,
	    -0.0402862, -0.0445232, -0.0492057, -0.0543807, -0.0601, -0.0664208,
	    -0.0734063, -0.0811265, -0.0896587, -0.0990881, -0.109509,-0.121027,
	    -0.133755, -0.147822, -0.163369, 67.0194, 1.79635, 3.59596, 5.39554,
	    7.19507, 8.99455, 10.794, 12.5933, 14.3926, 16.1919, 17.991,19.7901,
	    21.589, 23.3879, 25.1866, 26.9852, 28.7836, 30.5819, 32.38, 34.1779,
	    35.9756, 37.773, 39.5702, 41.367, 43.1635, 44.9597, 46.7555,48.5508,
	    50.3456, 52.1399, 53.9336, 55.7266, 57.5189, 59.3103, 61.1009,
	    62.8905, 64.679, 66.4662, 68.2522, 127.637
	}	;
	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-1 ; i++ ) {	// first block row
	    a[i][i] = 20.0	;	// diagonal entries
	    a[i][i+n-1] = -4.0	;	// super diagonal entries
	    if ( i==0 ) {	
		a[i][i+1] = -4.0;
		a[i][n] = -1.0	;
	    }
	    else if ( i==(n-2) ) {
		a[i][i-1] = -4.0;
		a[i][i+n-2]=-1.0;
	    }
	    else {
		a[i][i-1] = -4.0;
		a[i][i+1] = -4.0;
		a[i][i+n-2]=-1.0;
		a[i][i+n] = -1.0;
	    }
	}	// end for ( i ) loop
	for ( i=n-1 ; i<nelts-n+1 ; i++ ) {	// bulk m-3 blocks
	    a[i][i] = 20.0	;	// diagonal entries
	    a[i][i+n-1] = -4.0	;
	    a[i][i-n+1] = -4.0	;
	    if ( 0==(i%(n-1)) ) {
		a[i][i-n+2] = -1.0	;
		a[i][i+1] = -4.0	;
		a[i][i+n] = -1.0	;
	    }
	    else if ( (n-2)==(i%(n-1)) ) {
		a[i][i-n] = -1.0	;
		a[i][i-1] = -4.0	;
		a[i][i+n-2] = -1.0	;
	    }
	    else {
		a[i][i-n] = -1.0	;
		a[i][i-n+1] = -4.0	;
		a[i][i-n+2] = -1.0	;
		a[i][i-1] = -4.0	;
		a[i][i+1] = -4.0	;
		a[i][i+n-2] = -1.0	;
		a[i][i+n-1] = -4.0	;
		a[i][i+n] = -1.0	;
	    }
	}	// end of for ( i ) loop
	for ( i=nelts-n+1 ; i<nelts ; i++ ) {	// last block row
	    a[i][i] = 20.0	;
	    a[i][i-n+1] = -4.0	;
	    if ( i==(nelts-n+1) ) {
		a[i][i+1] = -4.0;
		a[i][i-n+2]=-1.0;
	    }
	    else if ( i==(nelts-1) ) {
		a[i][i-1] = -4.0;
		a[i][i-n] = -1.0;
	    }
	    else {
		a[i][i-1] = -4.0;
		a[i][i+1] = -4.0;
		a[i][i-n+2]=-1.0;
		a[i][i-n] = -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] *= w/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("SOR 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
//