/*
 * 最大最大激励化 第 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;
}

Result merge(Result* results) {
  Result r1 = results[0];
  Result r2 = results[1];

  if (r1.C + r1.n * r1.CM == r2.C) { // split n
    Result r = {r1.m, 2*r1.n, r1.CM, r1.C};
    return r;

  } else if (r1.C + r1.m == r2.C) { // split m
    Result r = {2*r1.m, r1.n, r1.CM, r1.C};
    return r;

  } else { // split k
    int x;
    for (x = 0; x < r1.n; x++) {
      cblas_daxpy(r2.m, 1, r2.C + r2.m * x, 1, r1.C + r1.CM * x, 1);
    }
    free(r2.C);
    Result r = {r1.m, r1.n, r1.CM, r1.C};
    return r;
  }
}

// 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);
        }
    }
}

int main()
{
    int n = 10;
    int in_x =1;
    int in_y =1;

    std::vector<double> x(n);
    std::vector<double> y(n);

    double alpha = 10;

    std::fill(x.begin(),x.end(),1.0);
    std::fill(y.begin(),y.end(),2.0);

    cblas_daxpy( n, alpha, &x[0], in_x, &y[0], in_y);

    //Print y 
    for(int j=0;j<n;j++)
        std::cout << y[j] << "\t";

    std::cout << std::endl;
}

int main () {
  typedef boost::multi_array<double, 1> vector;
  typedef vector::index vector_index;

  int N = 100000000;

  vector x(boost::extents[N]);
  vector y(boost::extents[N]);
  vector a(boost::extents[N]);

#pragma omp parallel for
  for (vector_index i = 0; i < N; i++) {
    x[i] = 1;
    y[i] = 1;
    a[i] = y[i];
  }

  cblas_daxpy(N, 1.0, &x[0], 1, &a[0], 1);

  std::cout << a[0] << '\n';

  return 0;
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_dmcsrm 
(JNIEnv * env, jobject calling_obj, jint M, jint N, jdoubleArray j_A, jint lda, 
 jdoubleArray j_B, jintArray j_ir, jintArray j_jc, jdoubleArray j_C, jint ldc){
	jdouble * A = (*env)->GetPrimitiveArrayCritical(env, j_A, JNI_FALSE);
	jdouble * B = (*env)->GetPrimitiveArrayCritical(env, j_B, JNI_FALSE);
	jint * ir = (*env)->GetPrimitiveArrayCritical(env, j_ir, JNI_FALSE);
	jint * jc = (*env)->GetPrimitiveArrayCritical(env, j_jc, JNI_FALSE);
	jdouble * C = (*env)->GetPrimitiveArrayCritical(env, j_C, JNI_FALSE);

    int ioff = jc[0];
    int i, j, k;
    for (i = 0; i < N; i++) {
      for (j = jc[i]-ioff; j < jc[i+1]-ioff; j++) {
        k = ir[j]-ioff;
        cblas_daxpy(M, B[j], A+(i*lda), 1, C+(k*ldc), 1);
      }
    }

	(*env)->ReleasePrimitiveArrayCritical(env, j_C, C, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_jc, jc, 0);	
    (*env)->ReleasePrimitiveArrayCritical(env, j_ir, ir, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_B, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_A, A, 0);
}

CPS_START_NAMESPACE
//--------------------------------------------------------------------
//  CVS keywords
//
//  $Source: /home/chulwoo/CPS/repo/CVS/cps_only/cps_pp/src/util/dirac_op/d_op_mobius/noarch/mobius_m.C,v $
//  $State: Exp $
//
//--------------------------------------------------------------------
//------------------------------------------------------------------
// mobius_m.C
//
// mobius_m is the fermion matrix.
// The in, out fields are defined on the checkerboard lattice
//
//------------------------------------------------------------------

CPS_END_NAMESPACE
#include<util/dwf.h>
#include<util/mobius.h>
#include<util/gjp.h>
#include<util/vector.h>
#include<util/verbose.h>
#include<util/error.h>
#include<util/dirac_op.h>
#include<util/time_cps.h>

#include "blas-subs.h"

CPS_START_NAMESPACE


//4d precond. mobius Dirac op:
// M_5 - kappa_b^2 M4eo M_5^-1 M4oe
void  mobius_m(Vector *out,
               Matrix *gauge_field,
               Vector *in,
               Float mass,
               Dwf *mobius_lib_arg)
{

    //------------------------------------------------------------------
    // Initializations
    //------------------------------------------------------------------
    const int f_size = 24 * mobius_lib_arg->vol_4d * mobius_lib_arg->ls / 2;
    const Float kappa_ratio = mobius_lib_arg->mobius_kappa_b/mobius_lib_arg->mobius_kappa_c;
    const Float minus_kappa_b_sq = -mobius_lib_arg->mobius_kappa_b * mobius_lib_arg->mobius_kappa_b;
    Vector  *frm_tmp2 = (Vector *) mobius_lib_arg->frm_tmp2;
    //Vector *temp = (Vector *) smalloc(f_size * sizeof(Float));
    Float norm;


    //  out = [ 1 + kappa_b/kappa_c 1/2 dslash_5  - kappa_b^2 Meo M5inv Moe] in
    // (dslash_5 uses (1+-g5), not P_R,L, i.e. no factor of 1/2 which is here out front)
    //    1. ftmp2 = Meo M5inv Moe in
    //    2. out <-  in
    //    3. out += -kappa_b^2 ftmp2
    //    4. out +=  -kappa_b/kappa_c /2 dslash_5 in
    //         (done by the dslash_5 with a5_inv = -kappa_b/kappa_c/2 *GJP.MobiusA5Inv() )


    //--------------------------------------------------------------
    //    1. ftmp2 = Meo M5inv Moe in
    //--------------------------------------------------------------
    // Apply Dslash O <- E
    //------------------------------------------------------------------
    time_elapse();
    mobius_dslash_4(out, gauge_field, in, 0, 0, mobius_lib_arg, mass);
    DEBUG_MOBIUS_DSLASH("mobius_dslash_4 %e\n", time_elapse());

    //------------------------------------------------------------------
    // Apply M_5^-1 (hopping in 5th dir + diagonal)
    //------------------------------------------------------------------
    mobius_m5inv(out, mass, 0, mobius_lib_arg);
    DEBUG_MOBIUS_DSLASH("mobius_m5inv %e\n", time_elapse());

    //------------------------------------------------------------------
    // Apply Dslash E <- O
    //------------------------------------------------------------------
    mobius_dslash_4(frm_tmp2, gauge_field, out, 1, 0, mobius_lib_arg, mass);
    DEBUG_MOBIUS_DSLASH("mobius_dslash_4 %e\n", time_elapse());

    //------------------------------------------------------------------
    //    2. out <-  in
    //------------------------------------------------------------------
#ifndef USE_BLAS
    moveFloat((IFloat*)out, (IFloat*)in, f_size);
#else
    cblas_dcopy(f_size, (IFloat*)in, 1, (IFloat*)out, 1);
#endif
    DEBUG_MOBIUS_DSLASH("out <- in %e\n", time_elapse());

    //------------------------------------------------------------------
    //    3. out += -kap2 ftmp2
    //------------------------------------------------------------------
#ifndef USE_BLAS
    fTimesV1PlusV2((IFloat*)out, minus_kappa_b_sq, (IFloat*)frm_tmp2,
                   (IFloat *)out, f_size);
#else

    cblas_daxpy(f_size, minus_kappa_b_sq, (IFloat*)frm_tmp2,1,
//.........这里部分代码省略.........

void plotMerit(double *z, double psi_k, double descentCondition)
{
  int incx = 1, incy = 1;
  double q_0, q_tk, qp_tk, merit_k;
  /* double tmin = 1e-12; */
  double tk = 1, aux;
  double m1 = 1e-4;
  double Nstep = 0;
  int i = 0;

  FILE *fp;

  (*sFphi)(sN, z, sphi_z, 0);
  aux = cblas_dnrm2(sN, sphi_z, 1);
  /* Computes merit function */
  aux = 0.5 * aux * aux;
  printf("plot psi_z %e\n", aux);


  if (!sPlotMerit)
    return;

  if (sPlotMerit)
  {
    /*    sPlotMerit=0;*/
    strcpy(fileName, "outputLS");


    (*sFphi)(sN, z, sphi_z, 0);
    q_0 =  cblas_dnrm2(sN, sphi_z , incx);
    q_0 = 0.5 * q_0 * q_0;

    fp = fopen(fileName, "w");

    /*    sPlotMerit=0;*/
    tk = 5e-7;
    aux = -tk;
    Nstep = 1e4;
    for (i = 0; i < 2 * Nstep; i++)
    {
      cblas_dcopy(sN, z, incx, sz2, incx);
      cblas_daxpy(sN , aux , sdir_descent , incx , sz2 , incy);
      (*sFphi)(sN, sz2, sphi_z, 0);
      q_tk =  cblas_dnrm2(sN, sphi_z , incx);
      q_tk = 0.5 * q_tk * q_tk;


      (*sFjacobianPhi)(sN, sz2, sjacobianPhi_z, 1);
      /* Computes the jacobian of the merit function, jacobian_psi = transpose(jacobianPhiMatrix).phiVector */
      cblas_dgemv(CblasColMajor,CblasTrans, sN, sN, 1.0, sjacobianPhi_z, sN, sphi_z, incx, 0.0, sgrad_psi_z, incx);
      qp_tk = cblas_ddot(sN, sgrad_psi_z, 1, sdir_descent, 1);

      merit_k = psi_k + m1 * aux * descentCondition;


      fprintf(fp, "%e %.16e %.16e %e\n", aux, q_tk, merit_k, qp_tk);
      if (i == Nstep - 1)
        aux = 0;
      else
        aux += tk / Nstep;
    }

    fclose(fp);
  }
}

int lineSearch_Wolfe(double *z, double qp_0)
{
  int incx = 1, incy = 1;
  double q_0, q_tk, qp_tk;
  double tmin = 1e-12;
  int maxiter = 100;
  int niter = 0;
  double tk = 1;
  double tg, td;
  double m1 = 0.1;
  double m2 = 0.9;


  (*sFphi)(sN, z, sphi_z, 0);
  q_0 =  cblas_dnrm2(sN, sphi_z , incx);
  q_0 = 0.5 * q_0 * q_0;

  tg = 0;
  td = 10e5;

  tk = (tg + td) / 2.0;

  while (niter < maxiter || (td - tg) < tmin)
  {
    niter++;
    /*q_tk = 0.5*|| phi(z+tk*d) ||*/
    cblas_dcopy(sN, z, incx, sz2, incx);
    cblas_daxpy(sN , tk , sdir_descent , incx , sz2 , incy);
    (*sFphi)(sN, sz2, sphi_z, 0);
    q_tk =  cblas_dnrm2(sN, sphi_z , incx);
    q_tk = 0.5 * q_tk * q_tk;

    (*sFjacobianPhi)(sN, sz2, sjacobianPhi_z, 1);
    /* Computes the jacobian of the merit function, jacobian_psi = transpose(jacobianPhiMatrix).phiVector */
    cblas_dgemv(CblasColMajor,CblasTrans, sN, sN, 1.0, sjacobianPhi_z, sN, sphi_z, incx, 0.0, sgrad_psi_z, incx);
    qp_tk = cblas_ddot(sN, sgrad_psi_z, 1, sdir_descent, 1);
    if (qp_tk <  m2 * qp_0 && q_tk < q_0 + m1 * tk * qp_0)
    {
      /*too small*/
      if (niter == 1)
        break;
      tg = tk;
      tk = (tg + td) / 2.0;
      continue;
    }
    else if (q_tk > q_0 + m1 * tk * qp_0)
    {
      /*too big*/
      td = tk;
      tk = (tg + td) / 2.0;
      continue;
    }
    else
      break;
  }

  cblas_dcopy(sN, sz2, incx, z, incx);

  if ((td - tg) <= tmin)
  {
    printf("NonSmoothNewton2::lineSearchWolfe warning, resulting tk < tmin, linesearch stopped.\n");
    return 0;
  }
  return 1;

}

void caffe_axpy<double>(const int N, const double alpha, const double* X,
                        double* Y) {
    cblas_daxpy(N, alpha, X, 1, Y, 1);
}

void eblas_daxpy_sub(size_t iStart, size_t iStop, double a, const double* x, int incx, double* y, int incy)
{	cblas_daxpy(iStop-iStart, a, x+incx*iStart, incx, y+incy*iStart, incy);
}

//.........这里部分代码省略.........
    *info = 2;
    return;
  }

  /*            End of Cholesky          */




  /*            Iteration loops  */

  iter1 = 0;
  err1 = 1.;


  while ((iter1 < itt) && (err1 > errmax))
  {

    /*       Linear stage (zc,wc) -> (z,w)*/


    alpha = 1.;
    beta  = 1.;
    cblas_dgemv(CblasColMajor,CblasTrans, n, n, alpha, k, n, zc, incx, beta, wc, incy);


    cblas_dcopy(n, problem->q, incx, znum1, incy);


    alpha = -1.;
    cblas_dscal(n , alpha , znum1 , incx);

    alpha = 1.;
    cblas_daxpy(n, alpha, wc, incx, znum1, incy);
    nrhs = 1;
    DTRTRS(LA_UP, LA_TRANS, LA_NONUNIT, n, nrhs, DPO, n, znum1, n, &info2);
    DTRTRS(LA_UP, LA_NOTRANS, LA_NONUNIT, n, nrhs, DPO, n, znum1, n, &info2);

    cblas_dcopy(n, znum1, incx, z, incy);



    alpha = -1.;
    beta = 1.;
    cblas_dgemv(CblasColMajor,CblasTrans, n, n, alpha, k, n, z, incx, beta, wc, incy);

    cblas_dcopy(n, wc, incx, w, incy);


    /*         Local Stage                  */

    cblas_dcopy(n, w, incx, wt, incy);
    alpha = -1.;
    beta = 1.;
    cblas_dgemv(CblasColMajor,CblasTrans, n, n, alpha, k, n, z, incx, beta, wt, incy);

    for (i = 0; i < n; i++)
    {
      if (wt[i] > 0.0)
      {
        wc[i] = wt[i];
        zc[i] = 0.0;
      }
      else
      {
        wc[i] = 0.0;

void get_top_delta(const da *m, const double *y, 
                   const double *x, double *d, const int batch_size){
    cblas_dcopy(batch_size * m->n_in, y, 1, d, 1);
    cblas_daxpy(batch_size * m->n_in, -1,
                x, 1, d, 1);
}

void train_da(da *m, dataset_blas *train_set, dataset_blas *expected_set, 
              int mini_batch, int n_epcho, char* weight_filename){
    int i, j, k, p, q;
    int epcho;
    double cost, total_cost;
    time_t start_time, end_time;
    FILE *weight_file;

    //weight_file = fopen(weight_filename, "w");

    for(epcho = 0; epcho < n_epcho; epcho++){

        total_cost = 0.0;
        start_time = time(NULL);
        for(k = 0; k < train_set->N / mini_batch; k++){

            if((k+1) % 500 == 0){
                printf("epcho %d batch %d\n", epcho + 1, k + 1);
            }
            get_hidden_values(m, train_set->input + k * mini_batch * m->n_in, h_out, mini_batch);
            get_reconstruct_input(m, h_out, z_out, mini_batch);
            
            get_top_delta(m, z_out, expected_set->input + k * mini_batch * m->n_in, d_high, mini_batch);
            get_second_delta(m, h_out, d_high, d_low, mini_batch);

            /* modify weight matrix W */
            cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
                        m->n_out, m->n_in, mini_batch,
                        1, d_low, m->n_out,
                        train_set->input + k * mini_batch * m->n_in, m->n_in,
                        0, tr1, m->n_in);

            cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
                        m->n_out, m->n_in, mini_batch,
                        1, h_out, m->n_out,
                        d_high, m->n_in, 0, tr2, m->n_in);

            cblas_daxpy(m->n_out * m->n_in, 1, tr1, 1, tr2, 1);
            
            cblas_daxpy(m->n_out * m->n_in, -eta / mini_batch, tr2, 1, m->W, 1);

            /* modify bias vector */
            cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
                        m->n_out, 1, mini_batch,
                        1, d_low, m->n_out,
                        Ivec, 1, 0, tr1, 1);

            cblas_daxpy(m->n_out, -eta / mini_batch, tr1, 1, m->b, 1);
            
            cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
                        m->n_in, 1, mini_batch,
                        1, d_high, m->n_in,
                        Ivec, 1, 0, tr1, 1);

            cblas_daxpy(m->n_in, -eta / mini_batch, tr1, 1, m->c, 1);

            for(i = 0; i < mini_batch * m->n_in; i++){
                tr1[i] = log(z_out[i]);
            }
            total_cost -= cblas_ddot(mini_batch * m->n_in, expected_set->input + k * mini_batch * m->n_in, 1,
                                     tr1, 1) / mini_batch;
            for(i = 0; i < mini_batch * m->n_in; i++){
                tr1[i] = log(1.0 - z_out[i]);
            }
            cblas_dcopy(mini_batch * m->n_in, Ivec, 1, tr2, 1);
            cblas_daxpy(mini_batch * m->n_in, -1, expected_set->input + k * mini_batch * m->n_in,
                        1, tr2, 1);
            total_cost -= cblas_ddot(mini_batch * m->n_in, tr1, 1, tr2, 1) / mini_batch;

        }
        end_time = time(NULL);
        printf("epcho %d cost: %.5lf\ttime: %ds\n", epcho + 1, total_cost / train_set->N * mini_batch, (int)(end_time - start_time));
    }

    //fclose(weight_file);
}

void bi_conjugate_gradient_sparse(cs *A, double *b, double* x, int n, double itol){
   
    int i,j,iter;
     
    double rho,rho1,alpha,beta,omega;
     
    double r[n], r_t[n];
    double z[n], z_t[n];
    double q[n], q_t[n], temp_q[n];
    double p[n], p_t[n], temp_p[n];
    double res[n];                  //NA VGEI!
    double precond[n];
     
    //Initializations      
    memset(precond, 0, n*sizeof(double));
    memset(r, 0, n*sizeof(double));
    memset(r_t, 0, n*sizeof(double));
    memset(z, 0, n*sizeof(double));
    memset(z_t, 0, n*sizeof(double));
    memset(q, 0, n*sizeof(double));
    memset(q_t, 0, n*sizeof(double));
    memset(temp_q, 0, n*sizeof(double));
    memset(p, 0, n*sizeof(double));
    memset(p_t, 0, n*sizeof(double));
    memset(temp_p, 0, n*sizeof(double));
    memset(res, 0, n*sizeof(double));
     
    /* Preconditioner */
    double max;
    int pp;
    for(j = 0; j < n; ++j){
        for(pp = A->p[j], max = fabs(A->x[pp]); pp < A->p[j+1]; pp++)
            if(fabs(A->x[pp]) > max)                  //vriskei to diagonio stoixeio
                max = fabs(A->x[pp]);
        precond[j] = 1/max;    
    }  
    cs *AT = cs_transpose (A, 1) ;
 
    cblas_dcopy (n, x, 1, res, 1);
 
    //r=b-Ax
    cblas_dcopy (n, b, 1, r, 1);
    memset(p, 0, n*sizeof(double));
    cs_gaxpy (A, x, p);
    for(i=0;i<n;i++){
        r[i]=r[i]-p[i];
     
    }
     
    cblas_dcopy (n, r, 1, r_t, 1);
     
    double r_norm = cblas_dnrm2 (n, r, 1);
    double b_norm = cblas_dnrm2 (n, b, 1);
    if(!b_norm)
        b_norm = 1;
 
    iter = 0;  
   
    while( r_norm/b_norm > itol && iter < n ){
       
        iter++;
 
        cblas_dcopy (n, r, 1, z, 1);            //gia na min allaksei o r
        cblas_dcopy (n, r_t, 1, z_t, 1);        //gia na min allaksei o r_t
        for(i=0;i<n;i++){
            z[i]=precond[i]*z[i];
            z_t[i]=precond[i]*z_t[i];
        }
     
        rho = cblas_ddot (n, z, 1, r_t, 1);    
        if (fpclassify(fabs(rho)) == FP_ZERO){
            printf("RHO aborting Bi-CG due to EPS...\n");
            exit(42);
        }
         
        if (iter == 1){
            cblas_dcopy (n, z, 1, p, 1);
            cblas_dcopy (n, z_t, 1, p_t, 1);
        }
        else{      
            //p = z + beta*p;
            beta = rho/rho1;           
 
            cblas_dscal (n, beta, p, 1);        //rescale p by beta
            cblas_dscal (n, beta, p_t, 1);      //rescale p_t by beta
         
            cblas_daxpy (n, 1, z, 1, p, 1);     //p = 1*z + p
            cblas_daxpy (n, 1, z_t, 1, p_t, 1); //p_t = 1*z_t + p_t
        }
         
        rho1 = rho;
         
        //q = Ap
        //q_t = trans(A)*p_t
        memset(q, 0, n*sizeof(double));
        cs_gaxpy (A, p, q);
        memset(q_t, 0, n*sizeof(double));
        cs_gaxpy(AT, p_t, q_t);        
         
        omega = cblas_ddot (n, p_t, 1, q, 1);
//.........这里部分代码省略.........

/* Ref: Weiss, Algorithm 11 CGS
 * INPUT
 *   n : dimension of the problem
 *   b [n] : r-h-s vector
 *   atimes (int n, static double *x, double *b, void *param) :
 *        calc matrix-vector product A.x = b.
 *   atimes_param : parameters for atimes().
 *   it : struct iter. following entries are used
 *        it->max = kend : max of iteration
 *        it->eps = eps  : criteria for |r^2|/|b^2|
 * OUTPUT
 *   returned value : 0 == success, otherwise (-1) == failed
 *   x [n] : solution
 *   it->niter : # of iteration
 *   it->res2  : |r^2| / |b^2|
 */
int
cgs (int n, const double *b, double *x,
     void (*atimes) (int, const double *, double *, void *),
     void *atimes_param,
     struct iter *it)
{
#ifndef HAVE_CBLAS_H
# ifdef HAVE_BLAS_H
  /* use Fortran BLAS routines */

  int i_1 = 1;
  double d_m1 = -1.0;
  double d_2 = 2.0;

# endif // !HAVE_BLAS_H
#endif // !HAVE_CBLAS_H

  int ret = -1;
  double eps2 = it->eps * it->eps;
  int itmax = it->max;

  double *r  = (double *)malloc (sizeof (double) * n);
  double *r0 = (double *)malloc (sizeof (double) * n);
  double *p  = (double *)malloc (sizeof (double) * n);
  double *u  = (double *)malloc (sizeof (double) * n);
  double *ap = (double *)malloc (sizeof (double) * n);
  double *q  = (double *)malloc (sizeof (double) * n);
  double *t  = (double *)malloc (sizeof (double) * n);
  CHECK_MALLOC (r,  "cgs");
  CHECK_MALLOC (r0, "cgs");
  CHECK_MALLOC (p,  "cgs");
  CHECK_MALLOC (u,  "cgs");
  CHECK_MALLOC (ap, "cgs");
  CHECK_MALLOC (q,  "cgs");
  CHECK_MALLOC (t,  "cgs");


  double r0ap;
  double rho, rho1;
  double delta;
  double beta;

  double res2 = 0.0;

#ifdef HAVE_CBLAS_H
  /**
   * ATLAS version
   */

  double b2 = cblas_ddot (n, b, 1, b, 1); // (b,b)
  eps2 *= b2;

  // initial residue
  atimes (n, x, r, atimes_param); // r = A.x
  cblas_daxpy (n, -1.0, b, 1, r, 1); // r = A.x - b

  cblas_dcopy (n, r, 1, r0, 1); // r0* = r
  cblas_dcopy (n, r, 1, p, 1); // p = r
  cblas_dcopy (n, r, 1, u, 1); // u = r

  rho = cblas_ddot (n, r0, 1, r, 1); // rho = (r0*, r)

  int i;
  for (i = 0; i < itmax; i ++)
    {
      atimes (n, p, ap, atimes_param); // ap = A.p
      r0ap = cblas_ddot (n, r0, 1, ap, 1); // r0ap = (r0*, A.p)
      delta = - rho / r0ap;

      cblas_dcopy (n, u, 1, q, 1); // q = u
      cblas_dscal (n, 2.0, q, 1); // q = 2 u
      cblas_daxpy (n, delta, ap, 1, q, 1); // q = 2 u + delta A.p

      atimes (n, q, t, atimes_param); // t = A.q

      cblas_daxpy (n, delta, t, 1, r, 1); // r = r + delta t
      cblas_daxpy (n, delta, q, 1, x, 1); // x = x + delta q

      res2 = cblas_ddot (n, r, 1, r, 1);
      if (it->debug == 2)
	{
	  fprintf (it->out, "libiter-cgs %d %e\n", i, res2 / b2);
	}
      if (res2 <= eps2)
//.........这里部分代码省略.........

void caffe_cpu_xpasv<double>(const int M, const int N, const double alpha,
    double* X, const double* a, const double* b) {
  for (int i = 0; i < M; ++i) {
    cblas_daxpy(N, alpha * a[i], b, 1, X + i * N, 1);
  }
}

//.........这里部分代码省略.........
    free(maxwt);
    free(zt);
    free(vectnt);
    free(ddln);
    free(ddlt);

    *info = 2;
    return;
  }

  /*                Iteration loops                  */


  iter1 = 0;
  err1  = 1.;



  while ((iter1 < itt) && (err1 > errmax))
  {

    /*               Linear stage (zc,wc) -> (z,w)         */

    alpha  = 1.;
    beta   = 1.;
    cblas_dgemv(CblasColMajor,CblasNoTrans, n, n, alpha, kfinv, n, zc, incx, beta, wc, incy);

    cblas_dcopy(n, qq, incx, znum1, incy);

    alpha = -1.;
    cblas_dscal(n , alpha , znum1 , incx);

    alpha = 1.;
    cblas_daxpy(n, alpha, wc, incx, znum1, incy);

    nrhs = 1;
    DTRTRS(LA_UP, LA_TRANS, LA_NONUNIT, n, nrhs, DPO, n, znum1, n, &info77);

    DTRTRS(LA_UP, LA_NOTRANS, LA_NONUNIT, n, nrhs, DPO, n, znum1, n, &info77);

    cblas_dcopy(n, znum1, incx, reaction, incy);

    alpha = -1.;
    beta = 1.;
    cblas_dgemv(CblasColMajor,CblasNoTrans, n, n, alpha, kfinv, n, reaction, incx, beta, wc, incy);

    cblas_dcopy(n, wc, incx, velocity, incy);



    /*               Local stage (z,w)->(zc,wc)          */


    for (i = 0; i < n; i++)
    {
      zc[i] = 0.;
      wc[i] = 0.0;
    }


    /*          Normal party                           */



    for (i = 0; i < nc; i++)
    {

int _globalLineSearchSparseGP(
  GlobalFrictionContactProblem *problem,
  AlartCurnierFun3x3Ptr computeACFun3x3,
  double *solution,
  double *direction,
  double *mu,
  double *rho,
  double *F,
  double *psi,
  CSparseMatrix *J,
  double *tmp,
  double alpha[1],
  unsigned int maxiter_ls)
{
  double inf = 1e10;
  double alphamin = 1e-16;
  double alphamax = inf;

  double m1 = 0.01, m2 = 0.99;

  unsigned int n = (unsigned)NM_triplet(problem->M)->m;

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

  unsigned int problem_size = n+2*m;

  // Computation of q(t) and q'(t) for t =0

  double q0 = 0.5 * cblas_ddot(problem_size, psi, 1, psi, 1);

  //  tmp <- J * direction
  cblas_dscal(problem_size, 0., tmp, 1);
  cs_gaxpy(J, direction, tmp);

  double dqdt0 = cblas_ddot(problem_size, psi, 1, tmp, 1);
  DEBUG_PRINTF("dqdt0=%e\n",dqdt0);
  DEBUG_PRINTF("q0=%e\n",q0);

  for(unsigned int iter = 0; iter < maxiter_ls; ++iter)
  {

    // tmp <- alpha*direction+solution
    cblas_dcopy(problem_size, solution, 1, tmp, 1);
    cblas_daxpy(problem_size, alpha[0], direction, 1, tmp, 1);

    ACPsi(
      problem,
      computeACFun3x3,
      tmp,  /* v */
      tmp+problem->M->size0+problem->H->size1, /* P */
      tmp+problem->M->size0, /* U */
      rho, psi);

    double q  = 0.5 * cblas_ddot(problem_size, psi, 1, psi, 1);

    assert(q >= 0);

    double slope = (q - q0) / alpha[0];

    int C1 = (slope >= m2 * dqdt0);
    int C2 = (slope <= m1 * dqdt0);

    DEBUG_PRINTF("C1=%i\t C2=%i\n",C1,C2);
    if(C1 && C2)
    {
      numerics_printf_verbose(1, "---- GFC3D - NSN_AC - global line search success. Number of ls iteration = %i  alpha = %.10e, q = %.10e",
                              iter,
                              alpha[0], q);
      
      return 0;

    }
    else if(!C1)
    {
      alphamin = alpha[0];
    }
    else
    {
      // not(C2)
      alphamax = alpha[0];
    }

    if(alpha[0] < inf)
    {
      alpha[0] = 0.5 * (alphamin + alphamax);
    }
    else
    {
      alpha[0] = alphamin;
    }

  }
  numerics_printf_verbose(1,"---- GFC3D - NSN_AC - global line search unsuccessful. Max number of ls iteration reached  = %i  with alpha = %.10e",
                  maxiter_ls, alpha[0]);
  

  return -1;
}

void fc3d_ProjectedGradientOnCylinder(FrictionContactProblem* problem, double *reaction, double *velocity, int* info, SolverOptions* options)
{
  /* int and double parameters */
  int* iparam = options->iparam;
  double* dparam = options->dparam;
  /* Number of contacts */
  int nc = problem->numberOfContacts;
  double* q = problem->q;
  NumericsMatrix* M = problem->M;
  /* Dimension of the problem */
  int n = 3 * nc;
  /* Maximum number of iterations */
  int itermax = iparam[0];
  /* Tolerance */
  double tolerance = dparam[0];




  /*****  Projected Gradient iterations *****/
  int j, iter = 0; /* Current iteration number */
  double error = 1.; /* Current error */
  int hasNotConverged = 1;
  int contact; /* Number of the current row of blocks in M */
  int nLocal = 3;
  dparam[0] = dparam[2]; // set the tolerance for the local solver
  double * velocitytmp = (double *)malloc(n * sizeof(double));

  double rho = 0.0;
  int isVariable = 0;
  double rhoinit, rhomin;
  if (dparam[3] > 0.0)
  {
    rho = dparam[3];
  }
  else
  {
    /* Variable step in fixed*/
    isVariable = 1;
    printf("Variable step (line search) in Projected Gradient iterations\n");
    rhoinit = dparam[3];
    rhomin = dparam[4];
  }

  double * reactionold;
  double * direction;
  if (isVariable)
  {
    reactionold = (double *)malloc(n * sizeof(double));
    direction = (double *)malloc(n * sizeof(double));
  }
  double alpha = 1.0;
  double beta = 1.0;

  /*   double minusrho  = -1.0*rho; */

  if (!isVariable)
  {
    while ((iter < itermax) && (hasNotConverged > 0))
    {
      ++iter;
      cblas_dcopy(n , q , 1 , velocitytmp, 1);
      prodNumericsMatrix(n, n, alpha, M, reaction, beta, velocitytmp);
      // projection for each contact
      cblas_daxpy(n, -1.0, velocitytmp, 1, reaction , 1);
      for (contact = 0 ; contact < nc ; ++contact)
        projectionOnCylinder(&reaction[ contact * nLocal],
                             options->dWork[contact]);

#ifdef VERBOSE_DEBUG

      printf("reaction before LS\n");
      for (contact = 0 ; contact < nc ; ++contact)
      {
        for (j = 0; j < 3; j++)
          printf("reaction[%i] = %le\t", 3 * contact + j, reaction[3 * contact + j]);
        printf("\n");
      }
      printf("velocitytmp before LS\n");
      for (contact = 0 ; contact < nc ; ++contact)
      {
        for (j = 0; j < 3; j++)
          printf("velocitytmp[%i] = %le\t", 3 * contact + j, velocitytmp[3 * contact + j]);
        printf("\n");
      }
#endif
      /* **** Criterium convergence **** */
      fc3d_Tresca_compute_error(problem, reaction , velocity, tolerance, options, &error);

      if (options->callback)
      {
        options->callback->collectStatsIteration(options->callback->env, nc * 3, 
                                        reaction, velocity, 
                                        error, NULL);
      }

      if (verbose > 0)
        printf("----------------------------------- FC3D - Projected Gradient On Cylinder (PGoC) - Iteration %i rho = %14.7e \tError = %14.7e\n", iter, rho, error);

      if (error < tolerance) hasNotConverged = 0;
//.........这里部分代码省略.........