//	library  matvec.cc (<-- printvec.cc)

#include <stdlib.h> 		// exit
#include <iostream.h>
#include <fstream.h>		// needed for ifstream and ofstream
#include <math.h>		// sqrt

//-------------------------------------------------------------------
//		Print elements of a vector (length n)
//
void PrintVec(const int n, const double x[])
{
 int i;
 
 cout << endl;
 for (i=0; i<n; i++)
  {
   cout << "  " << x[i];
  }
 cout << endl << endl; 
 return;
}
//-------------------------------------------------------------------
//		Read dimension from file, allocate vector, 
//				read vector components from file
void ConstructVector(const char filename[], int & n, double * & vec)
{
 int i;
 ifstream infile(filename);
 
 infile >> n;		// read dimension
 
 vec = new double [n];	// allocate memory for vector
 
 for (i=0; i<n; i++)
  {
    infile >> vec[i]; 
  }
 
 return;
}
//-------------------------------------------------------------------
//		Print elements of a matrix 
//			(nrow rows and ncol columns)
//
void PrintMat(const int nrow, const int ncol, const double a[])
{
 int i,j;
 
 cout << endl;
 for (i=0; i<nrow; i++)
  {
    cout << "Row " << i << ":";
    for (j=0; j<ncol; j++)
     {
       cout << "  " << a[i*ncol+j] ;
     }
    cout << endl; 
  }
 cout << endl << endl; 
 return;
}


//-------------------------------------------------------------------
//		Read dimensions from file, allocate matrix, 
//				read matrix components from file
void ConstructMatrix(const char filename[], 
			int & nrow, int & ncol, double * & mat)
{
 int i,j;
 ifstream infile(filename);
 
 infile >> nrow;			// read dimensions
 infile >> ncol;
 
 mat = new double [nrow*ncol];	// allocate memory for matrix
 
 for (i=0; i<nrow; i++)
  {
    for (j=0; j<ncol; j++)
     {
       infile >> mat[i*ncol+j]; 
     }
  }
 
 return;
}
//-------------------------------------------------------------------
//		matrix times vector : erg := A * u
//			with a flexible number of columns in matrix
void MatMul(const int nrow, const int ncol, 
  	      const double a[], const double u[], double erg[])
{
 int i,j;
 for (i=0; i<nrow; i++)
  {
    erg[i] = 0.0;			// initialize erg[i]
    for (j=0; j<ncol; j++)
     {
       erg[i] = erg[i] + a[i*ncol+j] * u[j];
     }
  } 
 return;
} 

//######################################################################
//	Get vector of diagonal entries from a matrix
//		(the vector is already allocated)
void GetMatrixDiag(const int nrow, const int ncol, const double a[], 
			double diag[])
{
 int i;
 
 if (nrow!=ncol)
  {
   cout << "GetMatrixDiag is implemented only for square matrices." << endl;
   exit(-10);
  }
 else
  {
   for (i=0; i<nrow; i++)
    {
     diag[i] = a[i*ncol+i];
    }
  }
 return;
}

//-------------------------------------------------------------------
//	Calculate defect (defect vector is already allocated)
//		Time consuming version
void GetDefect_aux(const int nrow, const int ncol, const double a[], 
		const double f[], const double u[], double defect[])
{
 int i;
 double *aux;			// aux. vector
 aux = new double [nrow];	// allocate defect vector (time consuming!!)
 
 MatMul(nrow, ncol, a, u, aux);	//  aux := A*u
 
 for (i=0; i<nrow; i++)
 {
  defect[i] = f[i] - aux[i];  	// d := f -aux
 }
 
 delete [] aux;

 return;
}
//-------------------------------------------------------------------
//	Calculate defect (defect vector is already allocated)
//		Sophisticated version
void GetDefect(const int nrow, const int ncol, const double a[], 
		const double f[], const double u[], double defect[])
{
 int i;
					// use defect vector as aux.vector
 MatMul(nrow, ncol, a, u, defect);	//  aux := A*u
 
 for (i=0; i<nrow; i++)
 {
  defect[i] = f[i] - defect[i]; 	// d := f -aux
 }

 return;
}
//-------------------------------------------------------------------
//	Calculate L2-norm of a vector
double L2norm(const int n, const double vec[])
{
 double norm = 0.0;
    int i;
 
 for (i = 0; i < n; i++)
   {
     norm += vec[i]*vec[i];
   }
 return sqrt(norm);
}

//-------------------------------------------------------------------
//		Print elements of a vector (length n) to a file
//
void PrintVec2File(const int n, const double x[], const char name[])
{
 ofstream outfile(name); 
 int i;
 
 outfile << n << endl;
 outfile << endl;
 for (i=0; i<n; i++)
  {
   outfile << x[i] << endl;
  }
 outfile << endl << endl; 
 return;
}