/**
 * @brief Normalizes all FSR scalar fluxes and Track boundary angular
 *        fluxes to the total fission source (times \f$ \nu \f$).
 */
void VectorizedSolver::normalizeFluxes() {

  FP_PRECISION* nu_sigma_f;
  FP_PRECISION volume;
  FP_PRECISION tot_fission_source;
  FP_PRECISION norm_factor;

  /* Compute total fission source for each FSR, energy group */
  #pragma omp parallel for private(volume, nu_sigma_f)  \
    reduction(+:tot_fission_source) schedule(guided)
  for (int r=0; r < _num_FSRs; r++) {

    /* Get pointers to important data structures */
    nu_sigma_f = _FSR_materials[r]->getNuSigmaF();
    volume = _FSR_volumes[r];

    /* Loop over energy group vector lengths */
    for (int v=0; v < _num_vector_lengths; v++) {

      /* Loop over each energy group within this vector */
      #pragma simd vectorlength(VEC_LENGTH)
      for (int e=v*VEC_LENGTH; e < (v+1)*VEC_LENGTH; e++) {
        _fission_sources(r,e) = nu_sigma_f[e] * _scalar_flux(r,e);
        _fission_sources(r,e) *= volume;
      }
    }
  }

  /* Compute the total fission source */
  int size = _num_FSRs * _num_groups;
  #ifdef SINGLE
  tot_fission_source = cblas_sasum(size, _fission_sources, 1);
  #else
  tot_fission_source = cblas_dasum(size, _fission_sources, 1);
  #endif

  /* Compute the normalization factor */
  norm_factor = 1.0 / tot_fission_source;

  log_printf(DEBUG, "Tot. Fiss. Src. = %f, Normalization factor = %f",
             tot_fission_source, norm_factor);

  /* Normalize the FSR scalar fluxes */
  #ifdef SINGLE
  cblas_sscal(size, norm_factor, _scalar_flux, 1);
  #else
  cblas_dscal(size, norm_factor, _scalar_flux, 1);
  #endif

  /* Normalize the Track angular boundary fluxes */
  size = 2 * _tot_num_tracks * _num_polar * _num_groups;

  #ifdef SINGLE
  cblas_sscal(size, norm_factor, _boundary_flux, 1);
  #else
  cblas_dscal(size, norm_factor, _boundary_flux, 1);
  #endif

  return;
}

double maxeig(double *xmat, mwSignedIndex n) 
{
	// xmat is symmetric n x n matrix
	mwSignedIndex incx=1,indmax,maxloop=10000,k=0;
	double alpha, beta, dmax, dmax_temp,dmax_tol;
	double *bufveca=(double *) calloc(n,sizeof(double));
	double *bufvecb=(double *) calloc(n,sizeof(double));

	dmax_tol=.001;
	// do power series to get approximation to max eigenvalue of A+X
	alpha=0.0;cblas_dscal(n,alpha,bufveca,incx);bufveca[0]=1.0; // x_0 = [1,0,0,...] do something better later
	beta=0.0;dmax=1.0;dmax_temp=0.0;
	while ((dabsf(dmax-dmax_temp)>dmax_tol)&&(k<=maxloop)){
		dmax_temp=dmax;
		alpha=1.0;
		cblas_dgemv(CblasColMajor,CblasNoTrans,n,n,alpha,xmat,n,bufveca,incx,beta,bufvecb,incx);
		indmax=idxmax(bufvecb,n);dmax=bufvecb[indmax];
		alpha=1.0/dmax;cblas_dscal(n,alpha,bufvecb,incx);
		cblas_dcopy(n,bufvecb,incx,bufveca,incx);
		k++;
	}
	alpha=1.0;	
	// compute Rayleigh Quotient to approximate max eigenvalue of A+X
	cblas_dgemv(CblasColMajor,CblasNoTrans,n,n,alpha,xmat,n,bufvecb,incx,beta,bufveca,incx);
	dmax=doubdot(bufvecb,bufveca,n);alpha=doubnorm2(bufvecb,n);dmax=dmax/alpha/alpha;
	free(bufveca);free(bufvecb);
	return dmax;
}

void get_class(crbm *m, double *h, double *py, int batch_size){
    int i;
    double sum;

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                batch_size, m->ncat, m->nhidden,
                1.0, h, m->nhidden, m->u, m->ncat,
                0, py, m->ncat);

    cblas_dger(CblasRowMajor, batch_size, m->ncat, 1,
               I, 1, m->by, 1, py, m->ncat);

    for(i = 0; i < batch_size * m->ncat; i++){
        py[i] = exp(py[i]);
    }

    //sum
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                batch_size, 1, m->ncat,
                1, py, m->ncat, I, 1,
                0, a, 1);

    for(i = 0; i < batch_size; i++){
        cblas_dscal(m->ncat, 1.0 / a[i], py + i * m->ncat, 1);
        //printf("sum:%.2lf\n", cblas_dasum(m->ncat, py + i * m->ncat, 1));
    }
}

void THBlas_(scal)(long n, real a, real *x, long incx)
{
  if(n == 1)
    incx = 1;

#if defined(USE_BLAS) && (defined(TH_REAL_IS_DOUBLE) || defined(TH_REAL_IS_FLOAT))
  if( (n <= INT_MAX) && (incx <= INT_MAX) )
  {
    int i_n = (int)n;
    int i_incx = (int)incx;

#if defined(TH_REAL_IS_DOUBLE)
    cblas_dscal(i_n, a, x, i_incx);
#else
    cblas_sscal(i_n, a, x, i_incx);
#endif
    return;
  }
#endif
  {
    long i;
    for(i = 0; i < n; i++)
      x[i*incx] *= a;
  }
}

JNIEXPORT void JNICALL Java_uncomplicate_neanderthal_CBLAS_dscal
(JNIEnv *env, jclass clazz,
 jint N, jdouble alpha,
 jobject X, jint offsetX, jint incX) {

    double *cX = (double *) (*env)->GetDirectBufferAddress(env, X);
    cblas_dscal(N, alpha, cX + offsetX, incX);
};

 /*
  * Class:     com_intel_analytics_bigdl_mkl_MKL
  * Method:    vdscal
  * Signature: (II[DII)V
  */
JNIEXPORT void JNICALL Java_com_intel_analytics_bigdl_mkl_MKL_vdscal
   (JNIEnv * env, jclass cls, jint n, jdouble a, jdoubleArray x, jint xOffset, jint incx) {

   jdouble * jni_x = (*env)->GetPrimitiveArrayCritical(env, x, JNI_FALSE);

   cblas_dscal(n, a, jni_x + xOffset, incx);

   (*env)->ReleasePrimitiveArrayCritical(env, x, jni_x, 0);
}

void eblas_dscal(int N, double a, double* x, int incx)
{
	#ifdef MKL_PROVIDES_BLAS
	cblas_dscal(N, a, x, incx);
	#else
	threadLaunch((N<100000) ? 1 : 0,
		eblas_dscal_sub, N, a, x, incx);
	#endif
}

/* Computes: y <- alpha A^T*x + beta y */
int mfiles_dgemv1(double alpha, const mxArray *A, const mxArray *x,
                  double beta, mxArray *y) {
    size_t rA = mxGetM(A);
    size_t cA = mxGetN(A);
    size_t rx = mxGetM(x);
    size_t cx = mxGetN(x);
    size_t ry = mxGetM(y);
    size_t cy = mxGetN(y);

    if (mxIsSparse(x) || mxIsSparse(y)) {
        mexErrMsgIdAndTxt("mfiles:BadType",
                          "Sparse vectors are not supported.");
    }

    if (mxIsComplex(A) || mxIsComplex(x) || mxIsComplex(y)) {
        mexErrMsgIdAndTxt("mfiles:BadType",
                          "Complex data is not supported.");
    }

    if ((rA != rx) || (cA != ry) || (cx != 1) || (cy != 1)) {
        mexErrMsgIdAndTxt("mfiles:BadDim",
                          "Dimensions of matrices do not match.");
    }

    if (mxIsSparse(A)) {
        double *px = mxGetPr(x);
        double *py = mxGetPr(y);
        double *pz = mxCalloc(ry, sizeof (double));

        cs *cs_A = cs_calloc(1, sizeof (cs));
        mfiles_mx2cs(A, cs_A);

        /* Transpose A */
        cs *cs_AT = cs_transpose(cs_A, 1);
        /* Compute z <- A^T*x */
        cs_gaxpy(cs_AT, px, pz);

        /* Compute y <- beta y */
        cblas_dscal(ry, beta, py, 1);
        /* Compute y <- alpha*z+y */
        cblas_daxpy(ry, alpha, pz, 1, py, 1);

        cs_free(cs_A); /* Check this cs_free and cs_spfree ? */
        cs_spfree(cs_AT);
        mxFree(pz);
    } else {
        double *pA = mxGetPr(A);
        double *px = mxGetPr(x);
        double *py = mxGetPr(y);
        cblas_dgemv(CblasRowMajor, CblasTrans,
                    rA, cA, alpha, pA, rA, px, 1, beta, py, 1);

    }

    return EXIT_SUCCESS;
}

 void Scope::process(const double* inputs)
 {
     double max = 1.;
     m_decoder->process(inputs, m_matrix);
     max = fabs(m_matrix[cblas_idamax(m_number_of_points, m_matrix, 1)]);
     if(max > 1.)
     {
         cblas_dscal(m_number_of_points, (1. / max), m_matrix, 1.);
     }
 }

void Space_trans::for_space()
{
  fftw_execute(pf);
  cblas_dscal(n3*md*2,1./md,realdata[0],1);

  /* for(int i = 0; i<n3*md; i++) */
  /* { */
  /*   realdata[i][0]/=md; */
  /*   realdata[i][1]/=md; */
  /* } */
  return;
}

/*
 * 最大最大激励化 第 unitdx 个单元
 *
 */
double LayerWiseRBMs::maximizeUnit(int layerIdx, int unitIdx, double * unitSample, double argvNorm, int epoch){

    int AMnumIn = layers[0]->numVis;                                            

    // unitsample 归一化
    double curNorm = squareNorm(unitSample, AMnumIn, 1);
    cblas_dscal(AMnumIn, argvNorm / curNorm, unitSample, 1);
	
	double maxValue =0;

	for(int k=0; k<epoch; k++){
	// forward
		for(int i=0; i<=layerIdx; i++){
			if(i==0)
				layers[i]->setInput(unitSample);
			else
				layers[i]->setInput(layers[i-1]->getOutput());
			layers[i]->setBatchSize(1);
			layers[i]->runBatch();
		}
		maxValue = layers[layerIdx]->getOutput()[unitIdx];
	//back propagate
		for(int i=layerIdx; i>=0; i--){
			if(i==layerIdx)
				layers[i]->getAMDelta(unitIdx, NULL)	;
			else
				layers[i]->getAMDelta(-1, layers[i+1]->AMDelta);
		}
        double lr = 0.01 * cblas_dasum(AMnumIn, unitSample, 1) /                
                    cblas_dasum(AMnumIn, layers[0]->AMDelta, 1);
		
	// update unitSample
		cblas_daxpy(AMnumIn, lr, layers[0]->AMDelta, 1, unitSample, 1);
	//归一化 unitSample
		curNorm = squareNorm(unitSample, AMnumIn, 1);
	    cblas_dscal(AMnumIn, argvNorm / curNorm, unitSample, 1);
	
	}
	return maxValue;
}

/* compute psi function */
void ACPsi(
  GlobalFrictionContactProblem* problem,
  AlartCurnierFun3x3Ptr computeACFun3x3,
  double *globalVelocity,
  double *reaction,
  double *velocity,
  double *rho,
  double *psi)
{

  assert(problem->H->size1 == problem->dimension * problem->numberOfContacts);
  unsigned int m = problem->H->size1;
  unsigned int n = problem->M->size0;
  unsigned int problem_size = n +  2*m ;

  cblas_dscal(problem_size, 0., psi, 1);

  /* -problem->M * globalVelocity + problem->H * reaction + problem->q ==> psi  */
  cblas_dscal(problem_size, 0., psi, 1);

  cblas_dcopy(n, problem->q, 1, psi, 1);
  NM_gemv(1., problem->H, reaction, 1, psi);
  NM_gemv(-1., problem->M, globalVelocity, 1, psi);

  /* -velocity + trans(problem->H) * globalVelocity + problem->b ==> psi + n */
  cblas_daxpy(m, -1., velocity, 1, psi + n, 1);
  cblas_daxpy(m, 1, problem->b, 1, psi + n, 1);
  NM_tgemv(1., problem->H, globalVelocity, 1, psi + n);

  /* compute AC function */
  fc3d_AlartCurnierFunction(m,
                            computeACFun3x3,
                            reaction,
                            velocity, problem->mu, rho,
                            psi+n+m,
                            NULL, NULL);

}

void Orthogonalize(OrthoContext* c, double* p, int numBases, double* orthonormalBases)
{
	memcpy(c->Pv->Data, p, c->Pv->Count * sizeof(double));
	memcpy(c->Bases->Data, orthonormalBases, numBases * c->Pv->Count * sizeof(double));
	c->Bases->RowCount = numBases;
	c->Dp->Count = numBases;

	int basisLen = c->Pv->Count;
	GEMV(1, c->Bases, c->Pv, 0, c->Dp);
	for (int i = 0, offset = 0; i < numBases; i++, offset += basisLen)
		AXPY2(-1 * c->Dp->Data[i], c->Bases->Data + offset, basisLen, c->Pv->Data);
	double mag = cblas_dnrm2(basisLen, c->Pv->Data, 1);
	cblas_dscal(basisLen, 1.0 / mag, c->Pv->Data, 1);
	memcpy(p, c->Pv->Data, basisLen * sizeof(double));
}

/* compute psi function */
void ACPsi(
  GlobalFrictionContactProblem* problem,
  AlartCurnierFun3x3Ptr computeACFun3x3,
  double *globalVelocity,
  double *reaction,
  double *velocity,
  double *rho,
  double *psi)
{

  assert(problem->H->size1 == problem->dimension * problem->numberOfContacts);

  unsigned int localProblemSize = problem->H->size1;

  unsigned int ACProblemSize = sizeOfPsi(NM_triplet(problem->M),
                                         NM_triplet(problem->H));

  unsigned int globalProblemSize = problem->M->size0;


  /* psi <-
       compute -problem->M * globalVelocity + problem->H * reaction + problem->q
     ... */
  cblas_dscal(ACProblemSize, 0., psi, 1);
  cblas_dcopy(globalProblemSize, problem->q, 1, psi, 1);
  NM_gemv(1., problem->H, reaction, 1, psi);
  NM_gemv(-1., problem->M, globalVelocity, 1, psi);


  /* psi + globalProblemSize <-
     compute -velocity + trans(problem->H) * globalVelocity + problem->b
   ... */
  cblas_daxpy(localProblemSize, -1., velocity, 1, psi + globalProblemSize, 1);
  cblas_daxpy(localProblemSize, 1, problem->b, 1, psi + globalProblemSize, 1);
  NM_tgemv(1., problem->H, globalVelocity, 1, psi + globalProblemSize);



  /* compute AC function */
  fc3d_AlartCurnierFunction(localProblemSize,
                                         computeACFun3x3,
                                         reaction,
                                         velocity, problem->mu, rho,
                                         psi+globalProblemSize+
                                         problem->H->size1,
                                         NULL, NULL);

}

// rescales the rows of A by the given weights
// weights only needs to be defined on the root process
void rescaleRows(double *localRowChunk, double *weights, distMatrixInfo *matInfo) {
    int mpi_rank = matInfo->mpi_rank;
    int numcols = matInfo->numcols;
    int localrows = matInfo->localrows;
    int * rowcounts = matInfo->rowcounts;
    int * rowoffsets = matInfo->rowoffsets;
    MPI_Comm *comm = matInfo->comm;

    double *localweights = (double *) malloc(localrows * sizeof(double));

    if(mpi_rank != 0) {
        MPI_Scatterv(NULL, rowcounts, rowoffsets, MPI_DOUBLE, localweights, localrows, MPI_DOUBLE, 0, *comm);
    } else {
        MPI_Scatterv(weights, rowcounts, rowoffsets, MPI_DOUBLE, localweights, localrows, MPI_DOUBLE, 0, *comm);
    }
    for(int rowIdx = 0; rowIdx < localrows; rowIdx = rowIdx + 1)
        cblas_dscal(numcols, localweights[rowIdx], localRowChunk + (rowIdx*numcols), 1);
    free(localweights);
}

/** Create the Householder symmetric matrix v*v^T
 *
 * This matrix is used to process the rows and columns of the input matrix.
 */
static void create_house_matrix_packed(size_t order, double shift, double *source, size_t incs, double *hhp) {
	double h[order];
	int i;

	/* zero out the destination hhp (it's packed aka triangular, so the size is
	 * non-square) */
	for (i = 0; i < order * (order + 1) / 2; i++) {
		hhp[i] = 0;
	}

	/* create and normalize householder vector h */
	cblas_dcopy(order, source, incs, h, 1);
	h[0] += MYSIGN(h[0]) * cblas_dnrm2(order, h, 1);
	h[0] -= shift;
	cblas_dscal(order, 1.0 / cblas_dnrm2(order, h, 1), h, 1);

	/* hhp = h h^T */
	cblas_dspr(CblasRowMajor, CblasUpper, order, 1.0, h, 1, hhp);
}

/* mB = mC   => C <- alpha*A + beta*C
   otherwise    C <- alpha*A + beta*B */
void geam(double alpha, Mat mA, double beta, Mat mB, Mat mC) {
  const int n = MatN(mA);
  const int n2 = MatN2(mA);
  const void* const a = MatElems(mA);
  const void* const b = MatElems(mB);
  void* const c = MatElems(mC);
  const bool dev = MatDev(mA);

  switch (MatElemSize(mA)) {
  case 4:
    if (dev) {
      float alpha32 = alpha, beta32 = beta;
      cublasSgeam(g_cublasHandle, CUBLAS_OP_N, CUBLAS_OP_N, n, n,
                  &alpha32, a, n, &beta32, b, n, c, n);
    } else {
      if (b == c) {
        cblas_sscal(n2, beta, c, 1);
      } else {
        memset(c, 0, MatSize(mC));
        cblas_saxpy(n2, beta, b, 1, c, 1);
      }
      cblas_saxpy(n2, alpha, a, 1, c, 1);
    }
    break;

  case 8:
    if (dev) {
      cublasDgeam(g_cublasHandle, CUBLAS_OP_N, CUBLAS_OP_N, n, n,
                  &alpha, a, n, &beta, b, n, c, n);
    } else {
      if (b == c) {
        cblas_dscal(n2, beta, c, 1);
      } else {
        memset(c, 0, MatSize(mC));
        cblas_daxpy(n2, beta, b, 1, c, 1);
      }
      cblas_daxpy(n2, alpha, a, 1, c, 1);
    }
    break;
  }
}

// The main chebyshev iteration
void cheby_solver_iterate(
        const int x,
        const int y,
        const int z,
        const int halo_depth,
        double alpha,
        double beta,
        double* vec_u,
        double* vec_u0,
        double* vec_p,
        double* vec_r,
        double* vec_w,
        double* vec_kx,
        double* vec_ky,
        double* vec_kz,
        int* a_row_index,
        int* a_col_index,
        double* a_non_zeros)
{
    int m = x*y*z;

    mkl_cspblas_dcsrgemv(
            "n", &m, a_non_zeros, a_row_index, a_col_index, vec_u, vec_w);

    int x_inner = x - 2*halo_depth;
    
#pragma omp parallel for
    for(int ii = halo_depth; ii < z-halo_depth; ++ii)
    {
        for(int jj = halo_depth; jj < y-halo_depth; ++jj)
        {
            const int offset = ii*x*y + jj*x + halo_depth;
            cblas_dcopy(x_inner, vec_u0 + offset, 1, vec_r + offset, 1);
            cblas_daxpy(x_inner, -1.0, vec_w + offset, 1, vec_r + offset, 1);
            cblas_dscal(x_inner, alpha, vec_p + offset, 1);
            cblas_daxpy(x_inner, beta, vec_r + offset, 1, vec_p + offset, 1);
        }
    }

    cheby_calc_u(x, y, z, halo_depth, vec_u, vec_p);
}

void THBlas_scale(long size, real alpha, real *y, long yStride)
{
  if(size == 1)
    yStride = 1;

#if USE_CBLAS
  if( (size < INT_MAX) && (yStride < INT_MAX) )
  {
#ifdef USE_DOUBLE
    cblas_dscal(size, alpha, y, yStride);
#else
    cblas_sscal(size, alpha, y, yStride);
#endif
    return;
  }
#endif
  {
    long i;
    for(i = 0; i < size; i++)
      y[i*yStride] *= alpha;
  }
}

// Calculates p
void cg_solver_calc_p(
        const int x,
        const int y,
        const int z,
        const int halo_depth,
        const double beta,
        double* vec_p,
        double* vec_r)
{
    int x_inner = x - 2*halo_depth;

#pragma omp parallel for
    for(int ii = halo_depth; ii < z-halo_depth; ++ii)
    {
        for(int jj = halo_depth; jj < y-halo_depth; ++jj)
        {
            const int offset = ii*x*y + jj*x + halo_depth;
            cblas_dscal(x_inner, beta, vec_p + offset, 1);
            cblas_daxpy(x_inner, 1.0, vec_r + offset, 1, vec_p + offset, 1);
        }
    }
}