本文整理汇总了C++中cblas_daxpy函数的典型用法代码示例。如果您正苦于以下问题:C++ cblas_daxpy函数的具体用法?C++ cblas_daxpy怎么用?C++ cblas_daxpy使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了cblas_daxpy函数的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: squareNorm
/*
* 最大最大激励化 第 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;
}
开发者ID:chengjunwen,项目名称:my_deeplearning,代码行数:44,代码来源:LayerWiseRBMs.cpp
示例2: merge
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;
}
}
开发者ID:pbirsinger,项目名称:CARDIO,代码行数:22,代码来源:carma.c
示例3: cg_solver_calc_p
// 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);
}
}
}
开发者ID:UoB-HPC,项目名称:TeaLeaf,代码行数:23,代码来源:cg_solver.c
示例4: main
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;
}
开发者ID:JhonM,项目名称:playground,代码行数:22,代码来源:blas_test.cpp
示例5: main
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;
}
开发者ID:barrymoo,项目名称:multi-array-testing,代码行数:23,代码来源:ma-cublas.cpp
示例6: cblas_daxpy
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);
}
开发者ID:Dcep,项目名称:BIDMat,代码行数:24,代码来源:BIDMat_CBLAS.c
示例7: mobius_m
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,
//.........这里部分代码省略.........
开发者ID:DeanHowarth,项目名称:QUDA-CPS,代码行数:101,代码来源:mobius_m.C
示例8: plotMerit
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);
}
}
开发者ID:bremond,项目名称:siconos,代码行数:65,代码来源:NonSmoothNewtonNeighbour.c
示例9: lineSearch_Wolfe
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;
}
开发者ID:bremond,项目名称:siconos,代码行数:66,代码来源:NonSmoothNewtonNeighbour.c
示例10: cblas_daxpy
void caffe_axpy<double>(const int N, const double alpha, const double* X,
double* Y) {
cblas_daxpy(N, alpha, X, 1, Y, 1);
}
开发者ID:flair2005,项目名称:mmd-caffe,代码行数:4,代码来源:math_functions.cpp
示例11: eblas_daxpy_sub
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);
}
开发者ID:yalcinozhabes,项目名称:pythonJDFTx,代码行数:3,代码来源:BlasExtra.cpp
示例12: lcp_latin
//.........这里部分代码省略.........
*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;
开发者ID:radarsat1,项目名称:siconos,代码行数:67,代码来源:lcp_latin.c
示例13: get_top_delta
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);
}
开发者ID:roles,项目名称:deep_learning,代码行数:6,代码来源:dpblas_da.c
示例14: train_da
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);
}
开发者ID:roles,项目名称:deep_learning,代码行数:76,代码来源:dpblas_da.c
示例15: bi_conjugate_gradient_sparse
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);
//.........这里部分代码省略.........
开发者ID:KwnstantinaLazaridou,项目名称:circuit_simulation,代码行数:101,代码来源:sparse_matrix_solution.c
示例16: cgs
/* 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)
//.........这里部分代码省略.........
开发者ID:ryseto,项目名称:demsd,代码行数:101,代码来源:cgs.cpp
示例17: cblas_daxpy
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);
}
}
开发者ID:ZhitingHu,项目名称:NN,代码行数:6,代码来源:math_functions.cpp
示例18: FrictionContact2D_latin
//.........这里部分代码省略.........
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++)
{
开发者ID:bremond,项目名称:siconos,代码行数:67,代码来源:FrictionContact2D_latin.c
示例19: _globalLineSearchSparseGP
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;
}
开发者ID:siconos,项目名称:siconos,代码行数:98,代码来源:gfc3d_nonsmooth_Newton_AlartCurnier.c
示例20: fc3d_ProjectedGradientOnCylinder
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;
//.........这里部分代码省略.........
开发者ID:xhub,项目名称:siconos,代码行数:101,代码来源:fc3d_ProjectedGradientOnCylinder.c
注:本文中的cblas_daxpy函数示例由纯净天空整理自Github/MSDocs等源码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。 |
请发表评论