/* iccbin -tpp2 -openmp -L/opt/intel/mkl/10.1.1.019/lib/64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -lpthread -lm */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <omp.h>

#define ITMAX 10000

double ddot_(int* , double* , int* , double* , int* );

main(int argc, char **argv)
{  if(argc < 2 || argc > 2)
   {  printf(" <N> \n");
      exit(0);
   }
   int N=atof(argv[1]); 
   if (2*(N/2) == N)
   {  printf("\n N must be odd!\n");
      exit(0);
   }
    

   int N2=N*N, Nb=N2/2, Nr=Nb+1, NAr=2*N*(N-1), NAb=NAr, NA=2*NAr;
   double s;
   double zb[Nb], zr[Nr], resb[Nb], resr[Nr];
/* Data structure for the matrix: simplified Morse data 
   structure since we know that the diagonal entries are 4
   and the off-diagonals are all -1*/
   double four=4, one=1;
   
/* jc: column location, j_num: number of nonzero elements */  
   int jcr[NAr], j_numr[Nr], jcb[NAb], j_numb[Nb];
/* j_ad: address of the beginning of each line */
   int j_adr[Nr+1], j_adb[Nb+1];
/* For the structure of the original matrix */
   int jc[NA], j_num[N2], p[N2];
   
   clock_t clock0, clock1;
  
   int i,j,k,l,m,it=0;
   double norm=N;
   int inc=1;

   /* Initializations: */
   for(i=0; i<Nr; i++)  {zr[i] = 0; j_adr[i]=0;}
   for(i=0; i<Nb; i++)  {zb[i] = 0; j_adb[i]=0;}
   j_adr[Nr]=0; j_adb[Nb]=0;
   
/* Morse structure of the off-diagonals blocks for the original matrix */
   i=0; k=0;
   for (j=0; j<N; j++)
   {  if(j==0)
	  {  j_num[N*i+j]=2;
	     jc[k++]=N*i+j+1; jc[k++]=N*(i+1)+j;
	  }
	  else if(j<N-1)
	  {  j_num[N*i+j]=3;
	     jc[k++]=N*i+j-1; jc[k++]=N*i+j+1; jc[k++]=N*(i+1)+j;
	  }	        
	  else
	  {  j_num[N*i+j]=2;
	     jc[k++]=N*i+j-1; jc[k++]=N*(i+1)+j;
	  }
   }   
   for (i=1; i<N-1; i++)
   {  for (j=0; j<N; j++)
	  {  if(j==0)
	     {  j_num[N*i+j]=3;
	        jc[k++]=N*(i-1)+j; jc[k++]=N*i+j+1; jc[k++]=N*(i+1)+j;
	     }
	     else if(j<N-1)
	     {  j_num[N*i+j]=4;
	        jc[k++]=N*(i-1)+j; jc[k++]=N*i+j-1; jc[k++]=N*i+j+1; jc[k++]=N*(i+1)+j;
	     }	        
	     else
	     {  j_num[N*i+j]=3;
	        jc[k++]=N*(i-1)+j; jc[k++]=N*i+j-1; jc[k++]=N*(i+1)+j;
	     }
      }   
   }   
   i=N-1;
   for (j=0; j<N; j++)
   {  if(j==0)
	  {  j_num[N*i+j]=2;
	     jc[k++]=N*(i-1)+j; jc[k++]=N*i+j+1;
	  }
	  else if(j<N-1)
	  {  j_num[N*i+j]=3;
	     jc[k++]=N*(i-1)+j; jc[k++]=N*i+j-1; jc[k++]=N*i+j+1;
	  }
	  else
	  {  j_num[N*i+j]=2;
	     jc[k++]=N*(i-1)+j; jc[k++]=N*i+j-1;
	  }	        
   }   

/* Permutation of original matrix => red-black matrix */
   for (i=0; i<Nb; i++)  
   {  p[2*i]=i; p[2*i+1]=Nr+i;}
   i=Nb; p[2*i]=i;
   
   k=0; l=0, m=0;
   for (i=0; i<Nb; i++)
   {  j_numr[i]=j_num[2*i]; j_numb[i]=j_num[2*i+1];
      for (j=0; j<j_numr[i]; j++) jcr[l++]=p[jc[k++]]-Nr;
      for (j=0; j<j_numb[i]; j++) jcb[m++]=p[jc[k++]];
   } 
   i=Nb;
   j_numr[i]=j_num[2*i]; 
   for (j=0; j<j_numr[i]; j++) jcr[l++]=p[jc[k++]]-Nr;

   int num_threads=omp_get_max_threads();
   
   /* Compute the "adresses" of the beginning of each line */
   for (i=1; i<=Nr; i++)
      for (j=0; j<i; j++) j_adr[i]+=j_numr[j];
   for (i=1; i<=Nb; i++)
      for (j=0; j<i; j++) j_adb[i]+=j_numb[j];
      
   
   clock0 = clock();
  
  /* Red-Black Gauss-Seidel */

   while (it<ITMAX && norm>1.e-6)
   {  
/* Update red nodes */ 
      for(i=0; i<Nr; i++)
      {  zr[i] = one;
	 for(j=j_adr[i]; j<j_adr[i+1]; j++) zr[i] += one*zb[jcr[j]];
	 zr[i] /= four;
      }
      
/* Update black nodes */ 
      for(i=0; i<Nb; i++)
      {  zb[i] = one;
	 for(j=j_adb[i]; j<j_adb[i+1]; j++) zb[i] += one*zr[jcb[j]]; 
	 zb[i] /= four;
      }
  
    /* residual */ 
      for(i=0; i<Nr; i++)
      {  resr[i] = one - four*zr[i];
	 for(j=j_adr[i]; j<j_adr[i+1]; j++) resr[i] += one*zb[jcr[j]]; 
      }
      norm = ddot_(&Nr,resr,&inc,resr,&inc);
      
      for(i=0; i<Nb; i++)
      {  resb[i] = one - four*zb[i];
	 for(j=j_adb[i]; j<j_adb[i+1]; j++) resb[i] += one*zr[jcb[j]];
      }
      norm += ddot_(&Nb,resb,&inc,resb,&inc);
      norm = sqrt(norm)/N;
            
      it++;
/*      printf("\n norm=%e    it=%d \n",norm,it);*/
   }

   clock1 = clock();
    
   printf("\n Norm of the residual: %e    Number of iterations: %d \n",norm,it);
   printf("\n zb[%d]=%e",Nb/2,zb[Nb/2]);
/*   for (i=0; i<Nb; i++) printf("\n zb[%d]=%e",i,zb[i]);*/
      
   printf("\n nproc:%d  CPU-time :  %e miliseconds\n",num_threads,
	 1.0e+3*(double)(clock1 - clock0)/CLOCKS_PER_SEC);
}