// 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



#define PARTICION 5000  // Tama~no de la particion
#define VISUALIZA 0    // (0) No visualiza la salida, otro valor la visualiza



// Indica si se carga lo referente a OPENMP
#ifdef _OPENMP
#include <omp.h>
  int  threads=omp_get_num_threads();
#else
  int  threads=0;
#endif
#include <math.h>
#include <stdio.h>



// Lado derecho de la ecuacion diferencial parcial
double LadoDerecho(double x)
{
   double pi = 3.1415926535897932384626433832;
   return -pi * pi * cos(pi * x);
}

// Resuelve Ax=b usando el metodo Jacobi
void Jacobi(double **A, double *x, double *b, int n, int iter)
{
   int i, j, m;
   double sum;
   double *xt = new double[n];
   #pragma omp parallel for
   for (i = 0; i < n; i++) xt[i] = 0.0;


   for (m = 0; m < iter; m ++)
   {
      sum =  A[0][2] * x[1];
      xt[0] = (1.0 / A[0][1]) * (b[0] - sum);
      #pragma omp simd reduction(+:sum)
      for (i = 1; i < n-1; i++) {
         sum = A[i][0] * x[i-1] +  A[i][2] * x[i+1];
         xt[i] = (1.0 / A[i][1]) * (b[i] - sum);
      }
      sum = A[i][0] * x[i-1];
      xt[i] = (1.0 / A[i][1]) * (b[i] - sum);
      #pragma omp parallel for
      for (i = 0; i < n; i++) x[i] = xt[i];
   }
   if (VISUALIZA)
   {
      printf("\nMatriz\n");
      for (i = 0; i < n; i++)
      {
         for (j = 0; j < 3; j++) printf("%f ", A[i][j]);
         printf("\n");
      }
      printf("\nb\n");
      for (i = 0; i < n; i++) printf("%f ", b[i]);
      printf("\nx\n");
      for (i = 0; i < n; i++) printf("%f ", x[i]);
      printf("\n");
   }
   delete []xt;
}


// Resuelve Ax=b usando el metodo Gauss-Seidel
void Gauss_Seidel(double **A, double *x, double *b, int n, int iter)
{
   int i, j, m;
   double sum;

   for (m = 0; m < iter; m ++)
   {
      sum =  A[0][2] * x[1];
      x[0] = (1.0 / A[0][1]) * (b[0] - sum);
      #pragma omp simd reduction(+:sum)
      for (i = 1; i < n-1; i++) {
         sum = A[i][0] * x[i-1] +  A[i][2] * x[i+1];
         x[i] = (1.0 / A[i][1]) * (b[i] - sum);
      }
      sum = A[i][0] * x[i-1];
      x[i] = (1.0 / A[i][1]) * (b[i] - sum);
   }
   if (VISUALIZA)
   {
      printf("\nMatriz\n");
      for (i = 0; i < n; i++)
      {
         for (j = 0; j < 3; j++) printf("%f ", A[i][j]);
         printf("\n");
      }
      printf("\nb\n");
      for (i = 0; i < n; i++) printf("%f ", b[i]);
      printf("\nx\n");
      for (i = 0; i < n; i++) printf("%f ", x[i]);
      printf("\n");
   }
}



int main()
{
   double xi = 0.0;               // Inicio del diminio
   double xf = 1.0;                // Fin del dominio
   double vi = 1.0;               // Valor en la frontera xi
   double vf = -1.0;                // Valor en la frontera xf
   int n = PARTICION;              // Particion
   int N = n-2;                    // Incognitas
   double h = (xf - xi) / (n - 1); // Incremento en la malla
   int i;

   double *b = new double [N];     // Vector b
   double *x = new double [N];     // Vector x
   double **A = new double *[N];   // Matriz A
   for(i = 0; i < N; i++) A[i] = new double[3];

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

   // Primer renglon de la matriz A y vector b
   A[0][1] = P;
   A[0][2] = Q;
   b[0] = LadoDerecho(xi) - vi * R;
   // Renglones intermedios de la matriz A y vector b
   #pragma omp parallel for
   for (i = 1; i < N - 1; i++)
   {
      A[i][0] = R;
      A[i][1] = P;
      A[i][2] = Q;
      b[i] = LadoDerecho(xi + h * (i+1));
   }
   // Renglon final de la matriz A y vector b
   A[N - 1][0] = R;
   A[N - 1][1] = P;
   b[N - 1] = LadoDerecho(xi + h * N - 2) - vf * Q;


   // Resuleve el sistema lineal Ax=b por Gauss-Seidel
   #pragma omp parallel for
   for(i = 0; i < N; i++) x[i] = 0.0;
   Gauss_Seidel(A, x, b, N, N*7);
   printf("%f %f\n", xi, vi);
   for(i = 1; i < N+1; i++) printf("%f %f\n", xi + h*i, x[i-1]);
   printf("%f %f\n\n\n\n", xf, vf);


   // Resuleve el sistema lineal Ax=b por Jacobi
   #pragma omp parallel for
   for(i = 0; i < N; i++) x[i] = 0.0;
   Jacobi(A, x, b, N, N*14);
   printf("%f %f\n", xi, vi);
   for(i = 1; i < N+1; i++) printf("%f %f\n", xi + h*i, x[i-1]);
   printf("%f %f\n", xf, vf);


   // Liberacion de memoria de matrices y vectores
   delete []b;
   delete []x;
   for(i = 0; i < N; i++) delete []A[i];
   delete A;

   return 0;
}