//	Jacobi.cc
//		Solve a system of equations via Richardson iteration

#include <stdlib.h> 		// exit
#include <iostream.h>
#include <assert.h>
#include "matvec.hh"

void JacobiSolve(const int n, const double a[], const double f[], double v[]);

int main()
{
// const char matrixfilename[] = "matrix1.in";
// const char    rhsfilename[] = "rhs1.in";
 const char matrixfilename[] = "matrix2.in";
 const char    rhsfilename[] = "rhs2.in";

 double *matrix, *u, *f;
    int nrow, ncol, n;

 ConstructMatrix(matrixfilename, nrow, ncol, matrix); 
 
 assert(nrow == ncol);	// Square matrix ??

 ConstructVector(rhsfilename, n, f);
 
 assert(n == ncol);	// Compatible dimensions ??
 
 u = new double [n];	// Allocate solution vector
 JacobiSolve(n, matrix, f, u); 
 PrintVec2File(n,u,"sol.txt");

 delete [] u;
 delete [] f;
 delete [] matrix;
 return 0;
}

//----------------------------------------------------------------
//	Solve A*u = f  via Jacobi iteration
//		A is square matrix
 
void JacobiSolve(const int n, const double a[], const double f[], double v[])
{
 const double EPS=1e-8;
       double *d, *diag;
       double err0, err;
          int iter, i;

 diag = new double [n];		// Allocate memory for aux.vectors
 d    = new double [n];

 GetMatrixDiag(n,n,a,diag);	// Get Diagonal

 for (i=0; i<n; i++)		// Initialize solution u
     v[i] = 0.0;

 for (i=0; i<n; i++)		// Initialize defect d
     d[i] = f[i]; 

 err = err0 = L2norm(n,d);	// Initialize error variable 
 iter = 0;
 
 while (err > EPS*err0) 
 {
  iter++;
  for (i=0; i<n; i++)		// Update solution
   {
    v[i] = v[i] +d[i]/diag[i];
   }
  GetDefect(n, n, a, f, v, d); 	// Update defect
  err = L2norm(n,d);		// Update error variable
 }
 
 cout << iter << " iterations for epsilon = " << EPS << endl;

 delete [] d;
 delete [] diag; 

 return;
}