// Ejemplo del Metodo de Diferencias Finitas para resolver la ecuacion diferencial parcial:
//    Uxx=-Pi*Pi*cos(Pix)
//        xi<=U<=xf
//        U(xi)=vi y U(xf)=vf
// con solucion analitica: 
//     cos(pix)
//
// Para conocer mas del metodo ver:
//    http://132.248.182.159/acl/MDF
//
// Autor: Antonio Carrillo Ledesma
//    http://mmc.geofisica.unam.mx/acl




#include <gmm/gmm.h>
#include <math.h>

const double  pi = 3.141592653589793;

// Lado derecho
double LD(double x)
{
   return ( -pi * pi * cos(pi * x));
}
 
// Solucion analitica
double SA(double x)
{
   return (cos(pi * x));
}



int main(void)
{
   int M=11;                 // Particion
   int N=M-2;                // Nodos interiores
   double xi=0.0;            // Inicio dominio
   double xf=1.0;            // Fin dominio
   double h=(xf-xi)/(M-1);   // Incremento en la malla
   double vi=1.0;            // Condicion inicial en el inicio del dominio
   double vf=-1.0;           // Condicion inicial en el fin del dominio


   // Matriz densa
   gmm::dense_matrix<double> AA(N, N);


   // Matriz dispersa
   gmm::row_matrix< gmm::rsvector<double> > A(N, N);

   // Vectores
   std::vector<double> x(N), b(N);


   int i;
   // Stencil
   double P = -2 / (h * h);
   double Q = 1 / (h * h);
   double R = 1 / (h * h);



   A(0, 0) = P; // Primer renglon de la matriz A y vector b
   A(0, 1) = Q;
   b[0] = LD(xi + h) - (vi / (h * h));

   for(i = 1; i < N - 1; i++) // Renglones intermedios de la matriz A y vector b
   {
      A(i, i - 1) = R;
      A(i, i) = P;
      A(i, i + 1) = Q;
      b[i] = LD(xi + (i + 1) * h);
   }

   A(N - 1, N - 2) = R; // Relglon final de la matriz A y vector b
   A(N - 1, N - 1) = P;
   b[N - 1] = LD(xi + (i + 1) * h) - (vf / (h * h));

   // Copia la matriz dispersa a la densa para usarla en LU
   gmm::copy(A,AA);

   // Visualiza la matriz y el vector
   std::cout << "Matriz A"<< AA << gmm::endl;
   std::cout << "Vector b"<< b << gmm::endl;


   // LU para matrices densa
   gmm::lu_solve(AA, x, b);
   std::cout << "LU"<< x << gmm::endl;



   gmm::identity_matrix PS;   // Optional scalar product for cg
   gmm::identity_matrix PR;   // Optional preconditioner
   gmm::iteration iter(10E-6);// Iteration object with the max residu
   size_t restart = 50;       // restart parameter for GMRES



   gmm::cg(A, x, b, PS, PR, iter); // Conjugate gradient
   std::cout << "CGM"<< x << std::endl;



   gmm::bicgstab(A, x, b, PR, iter); // BICGSTAB BiConjugate Gradient Stabilized
   std::cout << "BICGSTAB"<< x << std::endl;



   gmm::gmres(A, x, b, PR, restart, iter); // GMRES generalized minimum residual
   std::cout << "GMRES"<< x << std::endl;



   gmm::qmr(A, x, b, PR, iter); // Quasi-Minimal Residual method.
   std::cout << "Quasi-Minimal"<< x << std::endl;



   // Visualiza la solucion numerica
   std::cout << "Solucion Numerica"<< std::endl;
   std::cout << xi << "       " << vi << gmm::endl; // Valor de U en la frontera izquierda
   for(i = 0; i < N; i++)
   {
      std::cout << (i + 1)*h << "       " << x[i] << gmm::endl;
   }
   std::cout << xf << "       " << vf << gmm::endl; // Valor de U en la frontera derecha


   // Visualiza la solucion analitica
   std::cout << "Solucion Analitica"<< std::endl;
   std::cout << xi << "       " << SA(xi) << gmm::endl; // Valor de U en la frontera izquierda
   for(i = 0; i < N; i++)
   {
      std::cout << (i + 1)*h << "       " << SA((xi + (i + 1)*h)) << gmm::endl;
   }
   std::cout << xf << "       " << SA(xf) << gmm::endl; // Valor de U en la frontera derecha


   // Visualiza el error en valor absoluto en cada nodo
   std::cout << "Error en el calculo"<< std::endl;
   std::cout << xi << "       " << abs(vi - SA(xi))  << gmm::endl; // Valor de U en la frontera izquierda
   for(i = 0; i < N; i++)
   {
      std::cout << (i + 1)*h << "	" << fabs(x[i] - SA((xi + (i + 1)*h))) << gmm::endl;
   }
   std::cout << xf << "       " << abs(vf - SA(xf)) << gmm::endl; // Valor de U en la frontera derecha

   return 0;
}