C++ cblas_daxpy函数代码示例

/*
 * 最大最大激励化 第 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);
//.........这里部分代码省略.........