C++ ddot函数代码示例

void Update_lbfgs
   (int   n,               /* I  num unknowns              */
    int   t,               /* I  num vectors to store      */
    real  *s_vec,          /* I  x_new - x_old             */
    real  *y_vec,          /* I  g_new - g_old             */
    real  *Bs_vec,         /* IO scratch space             */
    int   *NUPDT,          /* IO num updates made to B/H   */
    real  *CMPS,           /* IO storage for t s_vecs      */
    real  *CMPY)           /* IO storage for t y_vecs      */
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*   This routine updates the limited memory BFGS approximations for
*   the Hessian B and its inverse H.  The update uses s_vec and
*   y_vec, whose product should be negative to preserve the positive
*   definiteness of the BFGS approximations.  If s'y is not a descent
*   direction, then y is "damped", an idea due to Powell.
*   Then y_vec becomes:
*        y_damped = alpha*y + (1-alpha)*Bs
*
*                                      0.8*s'Bs
*                       where alpha = ----------
*                                     s'Bs - s'y
*
*   If 0 <= s'y <= 1.0E-8 y'y then the update is skipped to prevent
*   division by a small number.
*
*   The L-BFGS structures are currently FORTRAN subroutines that
*   require the variables NUPDT, CMPS, and CMPY.  These are used as
*   work space storage and should not be altered between FORTRAN calls.
*********************************************************************/
{
  real  s_dot_y, sBs, y_dot_y, alpha;


/*-- Damp the estimate if necessary. */

  s_dot_y = ddot (n, s_vec, 1, y_vec, 1);
  multbv_ (&n, &t, s_vec, Bs_vec, NUPDT, CMPS, CMPY);
  sBs = ddot (n, s_vec, 1, Bs_vec, 1);

  if (s_dot_y < (0.2 * sBs)) {
    fprintf (bfgs_fp, "--- damping L-BFGS update\n");
    alpha = 0.8 * sBs / (sBs - s_dot_y);
    dscal (n, alpha, y_vec, 1);
    daxpy (n, (1.0 - alpha), Bs_vec, 1, y_vec, 1);
    s_dot_y = ddot (n, s_vec, 1, y_vec, 1);
  }

/*-- Decide whether to skip the update. */

  y_dot_y = ddot (n, y_vec, 1, y_vec, 1);
  if ((s_dot_y >= 0.0) && (s_dot_y <= (sqrt(MCHEPS) * y_dot_y))) {
    fprintf (bfgs_fp, "--- skipping L-BFGS update\n");
    return;
  }

/*-- Make the updates. */
  updtbh_ (&n, &t, s_vec, y_vec, NUPDT, CMPS, CMPY);

  return;
}

/*
 * Computes residuals.
 *
 * hrx = -A'*y - G'*z;  rx = hrx - c.*tau;  hresx = norm(rx,2);
 * hry = A*x;           ry = hry - b.*tau;  hresy = norm(ry,2);
 * hrz = s + G*x;       rz = hrz - h.*tau;  hresz = norm(rz,2);
 * rt = kappa + c'*x + b'*y + h'*z;
 */ 
void computeResiduals(pwork *w)
{
	/* rx = -A'*y - G'*z - c.*tau */
	if( w->p > 0 ) {
        sparseMtVm(w->A, w->y, w->rx, 1, 0);
        sparseMtVm(w->G, w->z, w->rx, 0, 0);
    } else {
        sparseMtVm(w->G, w->z, w->rx, 1, 0);
    }
	w->hresx = norm2(w->rx, w->n);
	vsubscale(w->n, w->tau, w->c, w->rx);
		
	/* ry = A*x - b.*tau */
	if( w->p > 0 ){
        sparseMV(w->A, w->x, w->ry, 1, 1);
        w->hresy = norm2(w->ry, w->p);
        vsubscale(w->p, w->tau, w->b, w->ry);
    } else {
        w->hresy = 0;
        w->ry = NULL;
	}
    
	/* rz = s + G*x - h.*tau */
	sparseMV(w->G, w->x, w->rz, 1, 1);
	vadd(w->m, w->s, w->rz);
	w->hresz = norm2(w->rz, w->m);
	vsubscale(w->m, w->tau, w->h, w->rz);

	/* rt = kappa + c'*x + b'*y + h'*z; */
	w->cx = ddot(w->n, w->c, w->x);
	w->by = w->p > 0 ? ddot(w->p, w->b, w->y) : 0.0;
	w->hz = ddot(w->m, w->h, w->z);
	w->rt = w->kap + w->cx + w->by + w->hz;    
}

double *cgsolve(int k)
{
	int i, first_i, last_i;
	int n = k * k;
	int maxiters = 1000 > 5*k ? 1000 : k;

	// partition data
	if (n % size) {
		first_i = (n / size + 1) * rank;
		last_i = (rank != size-1 ? first_i+n/size+1 : n);
	} else {
		first_i = n / size * rank;
		last_i = n / size * (rank + 1);
	}

	double *b_vec = (double *)malloc(n * sizeof(double));
	double *r_vec = (double *)malloc(n * sizeof(double));
	double *d_vec = (double *)malloc(n * sizeof(double));
	double *A_vec = (double *)malloc(n * sizeof(double));
	double *x_vec = (double *)malloc(n * sizeof(double));

	for (i=0; i<n; i++) {
		double tmp = cs240_getB(i, n);
		b_vec[i] = tmp;
		r_vec[i] = tmp;
		d_vec[i] = tmp;
		x_vec[i] = 0;
	}


	double normb = sqrt(ddot(b_vec+first_i, b_vec+first_i, last_i-first_i));
	double rtr = ddot(r_vec+first_i, r_vec+first_i, last_i-first_i);
	double relres = 1;

	i = 0;
	while (relres > 1e-6 && i++ < maxiters) {
	/*while (i++ < 1) {*/

		matvec(A_vec, d_vec, k);
		double alpha = rtr / ddot(d_vec+first_i, A_vec+first_i, last_i-first_i);
		daxpy(x_vec, d_vec, 1, alpha, n);
		daxpy(r_vec, A_vec, 1, -1*alpha, n);
		double rtrold = rtr;
		rtr = ddot(r_vec+first_i, r_vec+first_i, last_i-first_i);
		double beta = rtr / rtrold;
		daxpy(d_vec, r_vec, beta, 1, n);
		relres = sqrt(rtr) / normb;
	}
	return x_vec;
}

void loglike (int N, int W, int D, int T, double alpha, double beta, int *w, int *d, int **Nwt, int **Ndt, int *Nt, int *Nd) //
{
	int    i, j, t;
	double llike;
	static int init = 0;
	static double **prob_w_given_t;
	static double **prob_t_given_d;
	static double *Nd_;
	double Nt_;

	if (init==0) {
		init = 1;
		prob_w_given_t = dmat(W,T);
		prob_t_given_d = dmat(D,T);
		Nd_ = dvec(D);
		for (j = 0; j < D; j++) Nd_[j] = Nd[j] + T*alpha;
	}

	for (t = 0; t < T; t++) {
		Nt_ = Nt[t] + W*beta;
		for (i = 0; i < W; i++) prob_w_given_t[i][t] = (Nwt[i][t]+beta) / Nt_;
		for (j = 0; j < D; j++) prob_t_given_d[j][t] = (Ndt[j][t]+alpha)/ Nd_[j];
	}

	llike = 0;
	for (i = 0; i < N; i++)
		llike += log(ddot(T, prob_w_given_t[w[i]], prob_t_given_d[d[i]]));

	printf(">>> llike = %.6e    ", llike);
	printf("pplex = %.4f\n", exp(-llike/N));
}

static void task_ddot_xy_work_blocking( void * arg , TPI_ThreadPool pool )
{
    int p_size , p_rank ;

    if ( ! TPI_Rank( pool , & p_rank , & p_size ) ) {

        struct TaskXY * const t = (struct TaskXY *) arg ;

        const unsigned block_size   = t->block ;
        const unsigned block_start  = block_size * p_rank ;
        const unsigned block_stride = block_size * p_size ;

        const double * const x_end = t->x_beg + t->number ;
        const double * x = t->x_beg + block_start ;
        const double * y = t->x_beg + block_start ;

        double local = 0.0 ;

        for ( ; x < x_end ; x += block_stride , y += block_stride ) {
            const unsigned n = x_end - x ;
            local += ddot( ( block_size < n ? block_size : n ) , x , y );
        }

        t->xy_sum[ p_rank ] = local ;
    }
}

static int DvechmatDot(void* AA, double x[], int nn, int n, double *v){
  dvechmat* A=(dvechmat*)AA;
  ffinteger ione=1,nnn=nn;
  double dd,*val=A->AA->val;
  dd=ddot(&nnn,val,&ione,x,&ione);
  *v=2*dd*A->alpha;
  return 0;
}

double  KNITRO_EXPORT  KTR_ddot  (const int             n,
                                  const double * const  x,
                                  const int             incx,
                                  const double * const  y,
                                  const int             incy)
{
    return( ddot (n, x, incx, y, incy) );
}

void dgemv (double * A, double * x, double * y, int nrows, int ncols) {
// Calcula o produto de uma matriz por um vetor
	while (nrows -- > 0) {
		* y ++ = ddot(A, x, ncols);
		A += ncols;
		}
	return ;
	}

double dangle( LWDVector a, LWDVector b )
{
   LWDVector na, nb;

   dcopyv( na, a );
   dnormalize( na );
   dcopyv( nb, b );
   dnormalize( nb );
   return acos( ddot( na, nb ));
}

double inline
cblas_ddot(
  const int     n,
  const double* x,
  const int     incx,
  const double* y,
  const int     incy)
{
  return ddot(n, x, incx, y, incy);
}

/* ************************************************************
   TIME-CRITICAL PROCEDURE -- r=realssqr(x,n)
   Computes r=sum(x_i^2) using BLAS.
   ************************************************************ */
double realssqr(const double *x, const mwIndex n)
{
    mwIndex one=1;
    #ifdef PC
    return ddot(&n,x,&one,x,&one);
    #endif

    #ifdef UNIX
    return ddot_(&n,x,&one,x,&one);
    #endif    
}

double innerproduct(const Vector x, const Vector y)
{
  double res=ddot(&x->len, x->data, &x->stride, y->data, &y->stride);
#ifdef HAVE_MPI
  if (x->comm_size > 1) {
    double r2=res;
    MPI_Allreduce(&r2, &res, 1, MPI_DOUBLE, MPI_SUM, *x->comm);
  }
#endif

  return res;
}

double dotproduct(Vector u, Vector v)
{
  double locres;
  locres = ddot(&u->len, u->data, &u->stride, v->data, &v->stride);
  return locres;

#ifdef HAVE_MPI
  if (u->comm_size > 1) {
    MPI_Allreduce(&locres, &res, 1, MPI_DOUBLE, MPI_SUM, *u->comm);
    return res;
  }
#endif
}

/* compute a part of X*Y */
void mexFunction
(
    int nargout,
    mxArray *pargout [ ],
    int nargin,
    const mxArray *pargin  [ ]
)
{
    double *Xt, *Y, *Z, *I, *J, *v, LL;
    ptrdiff_t m, n, r, L, p, ir, jr, k;
    ptrdiff_t inc = 1;

    if (nargin != 5 || nargout > 1)
        mexErrMsgTxt ("Usage: v = partXY (Xt, Y, I, J, L)") ;

    /* ---------------------------------------------------------------- */
    /* inputs */
    /* ---------------------------------------------------------------- */
    
    Xt = mxGetPr( pargin [0] );
    Y  = mxGetPr( pargin [1] );
    I  = mxGetPr( pargin [2] );
    J  = mxGetPr( pargin [3] );
    LL = mxGetScalar( pargin [4] ); L = (ptrdiff_t) LL;
    m  = mxGetN( pargin [0] );
    n  = mxGetN( pargin [1] );
    r  = mxGetM( pargin [0] ); 
    if ( r != mxGetM( pargin [1] ))
        mexErrMsgTxt ("rows of Xt must be equal to rows of Y") ;
    if ( r > m || r > n )
        mexErrMsgTxt ("rank must be r <= min(m,n)") ;
    
    /* ---------------------------------------------------------------- */
    /* output */
    /* ---------------------------------------------------------------- */

    pargout [0] = mxCreateDoubleMatrix(1, L, mxREAL);
    v = mxGetPr( pargout [0] );
    
    /* C array indices start from 0 */
    for (p = 0; p < L; p++) {
        ir = ( I[p] - 1 ) * r;
        jr = ( J[p] - 1 ) * r;
        /* v[p] = 0.0;
        for (k = 0; k < r; k++)
            v[p] += Xt[ ir + k ] * Y[ jr + k ];*/
        v[p] = ddot(&r, Xt+ir, &inc, Y+jr, &inc);
    }
    
    return;
}

/* compute a part of X*Y */
void mexFunction
(
    int nargout,
    mxArray *pargout [ ],
    int nargin,
    const mxArray *pargin  [ ]
)
{
    if (nargin != 5 || nargout > 1)
        mexErrMsgTxt ("Usage: v = partXY (Xt, Y, I, J, L)") ;

    /* ---------------------------------------------------------------- */
    /* inputs */
    /* ---------------------------------------------------------------- */
    
    double *Xt = mxGetPr( pargin [0] );
    double *Y  = mxGetPr( pargin [1] );
    double *I  = mxGetPr( pargin [2] );
    double *J  = mxGetPr( pargin [3] );
    double  LL = mxGetScalar( pargin [4] ); 
    
    ptrdiff_t L = (ptrdiff_t) LL;
    ptrdiff_t m = mxGetN( pargin [0] );
    ptrdiff_t n = mxGetN( pargin [1] );
    ptrdiff_t r = mxGetM( pargin [0] ); 
    
    if ( r != mxGetM( pargin [1] ))
        mexErrMsgTxt ("rows of Xt must be equal to rows of Y") ;
    if ( r > m || r > n )
        mexErrMsgTxt ("rank must be r <= min(m,n)") ;
    
    /* ---------------------------------------------------------------- */
    /* output */
    /* ---------------------------------------------------------------- */

    pargout [0] = mxCreateDoubleMatrix(1, L, mxREAL);
    double *v = mxGetPr( pargout [0] );
    ptrdiff_t inc = 1;
    
    /* C array indices start from 0 */
    for (ptrdiff_t p = 0; p < L; p++) 
    {
        ptrdiff_t ir = ( I[p] - 1 ) * r;
        ptrdiff_t jr = ( J[p] - 1 ) * r;
        v[p] = ddot(&r, Xt+ir, &inc, Y+jr, &inc);
    }
    
    return;
}