

#include <stdio.h>
#include <mpl.h>
#include <math.h>

#define MAX_PROWS 16
#define MAX_PCOLS 8

#define MAXROWS 1024
#define MAXCOLS 1024

#define MAXBETAS 50

plural int Spins[MAX_PROWS][MAX_PCOLS];    /* spins */
plural int J[MAX_PROWS][MAX_PCOLS][2];     /* couplings */

double Beta_array[MAXBETAS];
double Delta_G[MAXBETAS];

int pRows, pCols;

extern double d_rand();

extern init_random();
extern plural double prand();

plural int Spin_posi ();
plural int Spin_posj ();
plural int Spin_negi ();
plural int Spin_negj ();

plural int J_negi ();
plural int J_negj ();


/*******************************************************************************
                                
                                MAIN PROGRAM

*******************************************************************************/

void main()
{
   int Sweep;
   int Lattice;  
   int NTherm;
   int Cur_Bet_Ptr;
   int Num_Betas;
   int Rows, Cols;     
   int Initial_Conf;     
   int Num_Spins;      
   int Freq;          
   int Num_Sweeps;     
   int Model;
   int Seed;
   int Beta_Accept[MAXBETAS];
   int Beta_Calcs[MAXBETAS];
     
   double Beta;          
   double Energy;       
   double Magnet;         
   double Accept; 
   double Total_Energy = 0.0;  
   double Total_Magnet = 0.0;   
   double Total_Accept = 0.0;   
   double Tot_Energy_Sq = 0.0;  
   double Tot_Magnet_Sq = 0.0;   
   double Tot_Accept_Sq = 0.0;   
   double Error_Energy;  
   double Error_Magnet;   
   double Error_Accept;   

   void Input_Parameters ();
   void Metropolis_Update ();
   void Beta_Update ();
   void Initialize_Lattice (); 
   void Calculate ();

   FILE *Outfile;


   Outfile = fopen ( "info", "w+");
   
   printf("Running on a %d x %d array of %d processors\n", 
	   nxproc, nyproc, nproc);
   fflush(stdout);

   Input_Parameters( &Num_Spins, &Rows, &Cols, &Initial_Conf,
		     &NTherm, &Freq, &Num_Sweeps, &Num_Betas, &Seed, &Model );

   /* Initialize random numbers */
   init_random( Seed );

   /* Check data decomposition */
   if ((Rows % nyproc) != 0) {
      printf("WARNING! Rows is not a multiple of nyproc\n");
   }
   if ((Cols % nxproc) != 0) {
      printf("WARNING! Cols is not a multiple of nxproc\n");
   }
   pCols = Cols / nxproc;
   pRows = Rows / nyproc;
   if (pRows < 2) {
      printf("WARNING! pRows is not > 2!!!\n");
   }
   if (pCols < 2) {
      printf("WARNING! pCols is not > 2!!!\n");
   }
   printf("Running a %d x %d subarray on each processors\n", pRows, pCols);
   printf("Initial_Conf = %d\n", Initial_Conf);
   fflush(stdout);


   Initialize_Lattice ( pRows, pCols, Initial_Conf, Model );

   Calculate ( pRows, pCols, &Energy, &Magnet );

   Energy /= (double)Num_Spins;
   Magnet /= (double)Num_Spins;

   printf("mag = %f    energy = %f\n", Magnet, Energy);
   fflush(stdout);

   fprintf ( Outfile, "THERMALIZATION SWEEPS :\n\n" );
   fprintf ( Outfile, "Acceptance \n" );
   fprintf ( Outfile, "------------\n\n" );
   fflush ( Outfile );

   for ( Sweep = 0; Sweep < NTherm; Sweep++ )
   {
         Beta = Beta_array[0];
         Metropolis_Update ( Beta, pRows, pCols, &Accept );

         fprintf ( Outfile, "%f\n", ( Accept/(double)Num_Spins ) ); 
   }

   fflush ( Outfile );
   fprintf ( Outfile, "\n\n" );
   fprintf ( Outfile, "BEGIN CALCULATIONS :\n\n" );
   fprintf ( Outfile, "LATTICE SIZE :  %d X %d = %d\n\n", 
		       Rows, Cols, Num_Spins );

   Calculate ( pRows, pCols, &Energy, &Magnet );

   Energy /= (double)Num_Spins;
   Magnet /= (double)Num_Spins;
   
   for ( Cur_Bet_Ptr = 0; Cur_Bet_Ptr < Num_Betas; Cur_Bet_Ptr++ )
   {
      Beta_Accept[Cur_Bet_Ptr] = 0.0;
      Beta_Calcs[Cur_Bet_Ptr] = 0.0;
   }
   Cur_Bet_Ptr = 0;

   for ( Sweep = 0; Sweep < Num_Sweeps; Sweep++ )
   {

      Energy *= (double)Num_Spins;
      Beta_Update ( &Cur_Bet_Ptr, Num_Betas, Energy, Beta_Accept );
      Beta_Calcs[Cur_Bet_Ptr] ++;
      Beta = Beta_array[Cur_Bet_Ptr];

      Total_Energy = 0.0;
      Total_Magnet = 0.0;
      Total_Accept = 0.0;

      Tot_Energy_Sq = 0.0;
      Tot_Magnet_Sq = 0.0;
      Tot_Accept_Sq = 0.0;

      for ( Lattice = 0; Lattice < Freq; Lattice++ )
      {
         Metropolis_Update ( Beta, pRows, pCols, &Accept );
      }
 
      Calculate ( pRows, pCols, &Energy, &Magnet );

      Energy /= (double)Num_Spins;
      Magnet /= (double)Num_Spins;
      Accept /= (double)Num_Spins;

      Total_Energy += Energy;
      Total_Magnet += Magnet;
      Total_Accept += Accept;

      Tot_Energy_Sq += Energy * Energy;
      Tot_Magnet_Sq += Magnet * Magnet;
      Tot_Accept_Sq += Accept * Accept;

      fprintf ( Outfile, "%f  %6d  %10.6lf\n", Beta, Sweep+1,  Energy ); 
      /*
      fprintf ( Outfile, "%6d \t %10.6lf \t %10.6lf \t %10.6lf\n", 
                         Sweep+1,  Energy, Magnet, Accept ); 
      */
      fflush ( Outfile );
   }

   Total_Energy /= Num_Sweeps;
   Total_Magnet /= Num_Sweeps;
   Total_Accept /= Num_Sweeps;

   Tot_Energy_Sq /= Num_Sweeps;
   Tot_Magnet_Sq /= Num_Sweeps;
   Tot_Accept_Sq /= Num_Sweeps;

   Error_Energy = Tot_Energy_Sq - (Total_Energy * Total_Energy);
   Error_Magnet = Tot_Magnet_Sq - (Total_Magnet * Total_Magnet);
   Error_Accept = Tot_Accept_Sq - (Total_Accept * Total_Accept);
   Error_Energy = sqrt ( Error_Energy / Num_Sweeps );
   Error_Magnet = sqrt ( Error_Magnet / Num_Sweeps );
   Error_Accept = sqrt ( Error_Accept / Num_Sweeps );
      
   fprintf ( Outfile, "\n\n" );
   fprintf ( Outfile, "Energy         = %10.8lf +/- %10.8lf\n", 
		      Total_Energy, Error_Energy ); 
   fprintf ( Outfile, "Magnetization  =  %10.8lf +/- %10.8lf\n", 
		      Total_Magnet, Error_Magnet ); 
   fprintf ( Outfile, "Acceptance     =  %10.8lf +/- %10.8lf\n\n", 
		      Total_Accept, Error_Accept ); 

   for ( Cur_Bet_Ptr = 0; Cur_Bet_Ptr < Num_Betas; Cur_Bet_Ptr++ )
      if ( Beta_Calcs[Cur_Bet_Ptr] > 0 ) 
      {
         fprintf ( Outfile, "%d  %f  %10.8lf\n", 
		            Cur_Bet_Ptr, Beta_array[Cur_Bet_Ptr],
			    (double)Beta_Accept[Cur_Bet_Ptr] / 
			    (double)Beta_Calcs[Cur_Bet_Ptr]); 
      }

   fclose ( Outfile );

}       
           


/*******************************************************************************

                           PROCEDURE INITIALIZE_LATTICE

*******************************************************************************/

void
Initialize_Lattice( Rows, Cols, Initial_Conf, Model )
int Rows, Cols, Initial_Conf, Model;
{
   int i, j, dim;

   switch (Initial_Conf)
   {

      case 1  :
           for ( i = 0; i < Rows; i++ )
              for ( j = 0; j < Cols; j++ )           
                 Spins[i][j] = 1;
	   printf("All spins 1\n");
           break;

      case 2  : 
           for ( i = 0; i < Rows; i++ )
              for ( j = 0; j < Cols; j++ )           
                 Spins[i][j] = -1;
	   printf("All spins -1\n");
           break;

      case 3  :
           for ( i = 0; i < Rows; i++ )
              for ( j = 0; j < Cols; j++ )           
                 if ( prand() < 0.5  )
                    Spins[i][j] = 1;
                 else 
                    Spins[i][j] = -1;
	   printf("Random spins\n");
           break;
 
      default  :
           printf("initialization of spins is screwed up!\n");
	   fflush(stdout);

   }

   switch (Model)
   {

      case 1  :  /* ferromagnet */
           for ( i = 0; i < Rows; i++ )
              for ( j = 0; j < Cols; j++ )           
                 for ( dim = 0; dim < 2; dim++ )           
                    J[i][j][dim] = 1;
	   printf("Ferromagnet\n");
           break;

      case 2  :  /* spin glass */
           for ( i = 0; i < Rows; i++ )
              for ( j = 0; j < Cols; j++ )           
                 for ( dim = 0; dim < 2; dim++ )           
                    if ( prand() < 0.5  )
                       J[i][j][dim] = 1;
                    else 
                       J[i][j][dim] = -1;
	   printf("Spin Glass\n");
           break;
 
      default  :
           printf("initialization of spins is screwed up!\n");
	   fflush(stdout);

   }

}


/*******************************************************************************

                            PROCEDURE CALCULATE

*******************************************************************************/

void
Calculate ( Rows, Cols, Energy, Magnet )
int Rows, Cols;
double *Energy, *Magnet;

{
   plural int Sum_Spins;
   int i, j;

   *Magnet = 0.0;
   *Energy = 0.0;

   for ( i = 0; i < Rows; i++ )
      for ( j = 0; j < Cols; j++ ) 
      {
         Sum_Spins = Spins[i][j] * ( J[i][j][0] * Spin_posi(i,j) 
				   + J[i][j][1] * Spin_posj(i,j) );

         *Magnet += (double) reduceAdd32(Spins[i][j]);
         *Energy += -1.0 * (double) reduceAdd32(Sum_Spins);
      }

}



/*******************************************************************************

                            PROCEDURE METROPOLIS_UPDATE

*******************************************************************************/

void
Metropolis_Update ( Beta, Rows, Cols, Accept )
double Beta;
int Rows, Cols;
double *Accept;
{
   plural int Old_Spin, New_Spin, Spin_Sum;
   plural double Old_Energy, New_Energy, Energy_Diff;
   plural int Acceptance;
   int i, j;

   *Accept = 0.0;
   
   for ( i = 0; i < Rows; i++ )
      for ( j = 0; j < Cols; j++ ) 
      {
         Acceptance = 0;

         Old_Spin = Spins[i][j];
         New_Spin = -Old_Spin;

         Spin_Sum = J[i][j][0]  * Spin_posi(i,j) + 
		    J[i][j][1]  * Spin_posj(i,j) + 
		    J_negi(i,j) * Spin_negi(i,j) + 
		    J_negj(i,j) * Spin_negj(i,j);

         Old_Energy = -Beta * Old_Spin * Spin_Sum;
         New_Energy = -Beta * New_Spin * Spin_Sum;
         Energy_Diff = New_Energy - Old_Energy;

         if ( (Energy_Diff <= 0.0) || ( p_exp(-Energy_Diff) > prand() ) )
         { 
            Spins[i][j] = New_Spin;
            Acceptance = Acceptance + 1;
         }

         *Accept += (double)reduceAdd32(Acceptance);
      }


}



/*******************************************************************************

			   PROCEDURE BETA_UPDATE

*******************************************************************************/

void
Beta_Update ( Cur_Bet_Ptr, Num_Betas, Energy, Accept_Beta )
int *Cur_Bet_Ptr;
int Num_Betas;
double Energy;
int Accept_Beta[MAXBETAS];
{
   int New_Bet_Ptr, Old_Bet_Ptr;
   double Boltz_Diff, Beta_Diff, G_Diff;

   Old_Bet_Ptr = *Cur_Bet_Ptr;
   
   /* Choose new neighboring beta value */
   if ( d_rand() < 0.5 ) {
      New_Bet_Ptr = Old_Bet_Ptr + 1;
      if ( New_Bet_Ptr >= Num_Betas) {
         New_Bet_Ptr = Num_Betas - 1;
      }
      G_Diff = Delta_G[Old_Bet_Ptr];
   }
   else {
      New_Bet_Ptr = Old_Bet_Ptr - 1;
      if ( New_Bet_Ptr < 0 ) {
         New_Bet_Ptr = 0;
      }
      G_Diff = - Delta_G[New_Bet_Ptr];
   }

   Beta_Diff = Beta_array[New_Bet_Ptr] - Beta_array[Old_Bet_Ptr];
   Boltz_Diff = (Energy * Beta_Diff) - G_Diff;

printf("B = %f  E = %f  E_dB = %f   G_Diff = %f\n",
	Beta_array[Old_Bet_Ptr], Energy, Energy * Beta_Diff, G_Diff );
printf("Boltz_Diff = %f\n", Boltz_Diff);
fflush(stdout);

   if ( ( Boltz_Diff <= 0.0 ) || ( exp ( -Boltz_Diff) > d_rand() ) ) {
            *Cur_Bet_Ptr = New_Bet_Ptr;
	    Accept_Beta[Old_Bet_Ptr] ++;
   }

}

 


/*******************************************************************************

                            PROCEDURE INPUT_PARAMETERS

*******************************************************************************/

void 
Input_Parameters( Num_Spins, Rows, Cols, Initial_Conf, 
		  NTherm, Freq, Num_Sweeps, Num_Betas, Seed, Model )
int *Num_Spins, *Rows, *Cols, *Initial_Conf;
int *NTherm, *Freq, *Num_Sweeps;
int *Num_Betas;
int *Seed, *Model;

{
   int i; 
   char Sentinel;

   do   
   {
      Sentinel = 0;
      printf ( "Enter the vertical dimensions of the lattice (0 < Y <= %d): ",
	       MAXROWS );
      scanf ( "%d", Rows );

      if ( ( *Rows > MAXROWS ) || ( *Rows <= 0 ) )
      { 
         printf ( "\n\nInput out of range.\n\n" );
         Sentinel = 1;
      }
   } 
   while (Sentinel);

   do
   {
      Sentinel = 0;

      printf("\n\n");
      printf ( "Enter the horizontal dimensions of the lattice (0 < X <= %d): ",
	       MAXCOLS );
      scanf ( "%d", Cols );

      if ( ( *Cols > MAXCOLS ) || ( *Cols <= 0 ) )
      {
            printf ( "\n\nInput out of range.\n\n" );
            Sentinel = 1;
      }
   }
   while (Sentinel);

   *Num_Spins  = (*Rows) * (*Cols);
   printf ( "\n\nLATTICE SIZE :    %d X %d = %d\n", *Rows, *Cols, *Num_Spins );

   do 
   {
      Sentinel = 'n';

/*      printf ( "Enter the value of Beta: " );
      scanf ( "%lf", Beta );
*/
      printf ( "\n\nEnter the number of thermalization sweeps desired: " );
      scanf ( "%d", NTherm );
      if ( *NTherm < 0 ) *NTherm = 0;

      do
      {
         Sentinel = 'n';        

         printf ( "\n\nEnter the total number of calculation updates \n" );
         printf ( "desired for each Beta: : " );
         scanf ( "%d", Num_Sweeps );
         if ( *Num_Sweeps <= 0 )   
            {
               printf ( "\n\nInvalid input\n\n" );
               Sentinel = 'y';
            }
      } while (Sentinel == 'y');
     
      do
      { 
         Sentinel = 'n';

         printf ( "\n\nEnter the number of sweeps per calculation update: ");
         scanf ( "%d", Freq );
         if ( *Freq > *Num_Sweeps ) *Freq = *Num_Sweeps;
         if ( *Freq <= 0 ) 
         { 
            printf ( "\n\nInvalid input\n\n" );
            Sentinel = 'y';
         }   
      } while ( Sentinel == 'y' ); 
     
      printf ( "\n\nCurrent number of thermalization sweeps  = %d\n", *NTherm );
      printf ( "Current total number of sweeps           = %d\n", *Num_Sweeps );
      printf ( "Current number of sweeps/calculation     = %d\n\n", *Freq );
      printf ( "Are all values correct (y/n)?" );
      scanf ( "%1s",&Sentinel );
      printf ( "\n\n" );
   } while ( ( Sentinel != 'y' ) && ( Sentinel != 'Y' ) );
  
      printf ( "Enter the Beta values in order.\n" );
      printf ( "Input -1 to end input. ( # values <= %d )\n\n", MAXBETAS );

      *Num_Betas = 0;
      i = 0;
      do
      {
         scanf ( "%lf", Beta_array + i );
         (*Num_Betas)++;
         i++;
      } 
      while ( *( Beta_array + i - 1 ) > 0.0 );

      (*Num_Betas)--;
     
      printf ( "CURRENT BETA VALUES\n" );
      printf ( "-------------------\n" );
 
      for ( i = 0; i < *Num_Betas; i++ ) 
      {
         printf ( "%14.8lf\n", *( Beta_array + i ) );
      }
      fflush(stdout);
     
      printf ( "\n\nCURRENT G_DIFF VALUES\n" );
      printf ( "-------------  ------\n" );
 
      for ( i = 0; i < *Num_Betas; i++ ) 
      {
         scanf ( "%lf", Delta_G + i );
         printf ( "%14.8lf\n", *( Delta_G + i ) );
      }
 
   do 
   {
      Sentinel = 0;
      printf ( "\n\nEnter the number next to your choice of spin\n" );
      printf ( "arrangements for the initial lattice.\n\n" );
      printf ( "1 -- All spins up.\n" );
      printf ( "2 -- All spins down.\n" );
      printf ( "3 -- Random configuration.\n\n\n" );
      printf ( "Choice: " );
      scanf ( "%d", Initial_Conf );

      if ( (*Initial_Conf != 1) && (*Initial_Conf != 2) && (*Initial_Conf != 3) )
      {
         printf ( "\n\nInvalid input\n" );
         Sentinel = 1;
      }
   }
   while (Sentinel);
   
   do 
   {
      Sentinel = 0;
      printf ( "\n\nWhich model do you want to simulate?\n\n" );
      printf ( "1 -- Ferromagnet.\n" );
      printf ( "2 -- Spin Glass.\n" );
      printf ( "Choice: " );
      scanf ( "%d", Model );

      if ( (*Model != 1) && (*Model != 2) )
      {
         printf ( "\n\nInvalid input\n" );
         Sentinel = 1;
      }
   }
   while (Sentinel);
   fflush(stdout);

   printf ( "\n\nInput a seed for the random number generator:\n\n" );
   scanf ( "%d", Seed );

}


plural int Spin_posi ( i, j ) 
int i, j;
{
   if (i < pRows-1) return Spins[i+1][j];
   else             return xnetS[1].Spins[0][j];
}

plural int Spin_posj ( i, j ) 
int i, j;
{
   if (j < pCols-1) return Spins[i][j+1];
   else             return xnetE[1].Spins[i][0];
}


plural int Spin_negi ( i, j ) 
int i, j;
{
   if (i > 0) return Spins[i-1][j];
   else       return xnetN[1].Spins[pRows-1][j];
}

plural int Spin_negj ( i, j ) 
int i, j;
{
   if (j > 0) return Spins[i][j-1];
   else       return xnetW[1].Spins[i][pCols-1];
}


plural int J_negi ( i, j ) 
int i, j;
{
   if (i > 0) return J[i-1][j][0];
   else       return xnetN[1].J[pRows-1][j][0];
}

plural int J_negj ( i, j ) 
int i, j;
{
   if (j > 0) return J[i][j-1][1];
   else       return xnetW[1].J[i][pCols-1][1];
}

