/*====================================================================================
! Purpose(s): 
!
! Contain(s): 
!
! Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
!             Gabriel Barajas   (barajasg@unican.es)
!====================================================================================*/

#include <malloc.h>
#include <stdlib.h>
#include <memory.h>

extern "C"  void csolve(double *AC, double *AL, double *AR, double *AT, double *AB, double *AX, double *AZ, 
                        int AWidth, int AHeight)
/*!=======================================================================
! Purpose(s): Solves the equation A*X=Z, where A is defined by the coefficients at the center (AC), left (AL),
!             Right (AR), Top (AT) and bottom (AB).
!
! Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
!             Gabriel Barajas   (barajasg@unican.es)
!=======================================================================*/
{
  static double *ME=NULL, *MXn=NULL;
  double *C, *L, *R, *T, *B, *Z, *X, *E, *Xn;
  double prov;
  int Tot=AWidth*AHeight;
  
  // The first time reserves memory space for working arrays.
  if (ME==NULL)  ME=(double *)malloc(Tot*sizeof(double));
  if (MXn==NULL) MXn=(double *)malloc(Tot*sizeof(double));
  
  for (int n=1;n<30000;n++)
  {
  
    // Copies X over Xn.
    memcpy(MXn, AX, Tot*sizeof(double));
  
    C=AC;
    L=AL;
    R=AR;
    T=AT;
    B=AB;
    Z=AZ;
    X=AX;
    E=ME;
    Xn=MXn;
    for (int i=1;i<=Tot;i++)
    {
  	  prov=*Z;
      if (*L!=0.0) prov=prov-*L*(*(Xn-1));
      if (*R!=0.0) prov=prov-*R*(*(Xn+1));
      if (*T!=0.0) prov=prov-*T*(*(Xn+AWidth));
      if (*B!=0.0) prov=prov-*B*(*(Xn-AWidth));
    
      *X=prov/ *C;

      C++;
      L++;
      R++;
      T++;
      B++;
      Z++;
      X++;
      E++;
      Xn++;

    } //for
  } //for n
}


//----------------------------------------------------------------------------------------


extern "C"  void cmatmul9(double *AE, double *AD, double *AC, double *AB, double *AA, double *ABU,
                          double *ACU, double *ADU, double *AEU, int AWidth, int AHeight,  double *AX, double *AY)
{
/*!=======================================================================
! Purpose(s): The first and last AWidth+1 rows doesn´t have complete diagonals, so we multiply them in a different
!             fashion.
!
! Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
!             Gabriel Barajas   (barajasg@unican.es)
!=======================================================================*/  
  double register *x, *y;
  double *a, *b, *c, *d, *e, *bu, *cu, *du, *eu;
} 
  
//----------------------------------------------------------------------------------------