/*============================================================================

 Function      vfiCPU

 Usage         vfiCPU(V, G)

 Arguments     V: pointer to array of doubles storing the value function.
                   
               G: pointer to array of doubles storing the policy function.
	              
 Description   This function performs value function iteration on the CPU,
               finding the maximum of the Bellman objective function for each
	       node in the state space and iterating until convergence.

 Dependencies  Global Variables: nk, nz, tol, howard (globalvas.h).

               Functions:        pow, fabs (math.h);
	                         cblas_daxpy, cblas_idamax (cblas.h);
	                         ar1, kGrid, vfInit, vfStep (auxfuncs.h).

 Return value  Returns 0 upon successful completion, 1 otherwise.

 =============================================================================

 Author:       Eric M. Aldrich

 Contact:      ealdrich@gmail.com

 Date:         20 April 2010

 ============================================================================*/

#include "globalvars.h"
#include "auxfuncs.h"
#include <math.h>
#include <ctime>
extern "C" {
  #include "cblas.h"
}

int vfiCPU(double* V, double* G)
{

  // admin
  int imax;
  double diff = 1.0;
  int i, j, l;
  clock_t start;

  // allocate variables in host memory
  double* K = new double[nk];
  double* Z = new double[nz];
  double* P = new double[(int)pow(nz,2)];
  double* V0 = new double[nk*nz];
  double* Vtemp;

  // compute TFP grid, capital grid and initial VF
  double lambda = 3;
  ar1(lambda, Z, P);
  kGrid(Z, K);
  vfInit(Z, V0);

  // various arguments to vfStep
  double* Exp = NULL;
  int klo, khi, nksub;
  double ydepK;
  double w, wmax; // if grid
  int windmax; // if grid
  int kslo, kshi, ksmid1, ksmid2; // if binary
  double w1, w2, w3; // if binary

  // iterate
  int count = 0;
  bool how = false;
  while(fabs(diff) > tol){
  //while(count < 8){
    if(count < 3 | count % howard == 0) how = false; else how = true;
    vfStep(klo, khi, nksub, kslo, kshi, ksmid1, ksmid2, w, wmax, windmax, w1,
	   w2, w3, i, j, l, ydepK, how, K, Z, P, Exp, V0, V, G);
    cblas_daxpy(nk*nz, -1.0, V, 1, V0, 1);
    imax = cblas_idamax(nk*nz, V0, 1);
    diff = *(V0+imax);
    Vtemp = V0;
    V0 = V;
    V = Vtemp;
    ++count;
  }
  //V = V0; // this assignment doesn't work sometimes (it works in GPU code)

  // i resort to a full copy in lieu of pointer reassignemnt
  for(i = 0 ; i < nk ; ++i){
    for(j = 0 ; j < nz ; ++j){
      V[i*nz+j] = V0[i*nz+j];
    }
  }

  return 0;

}
