C++ VecScale函数代码示例

PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
{
  PetscErrorCode ierr;
  PetscInt       mx;
  PetscScalar    h;
  Vec            x;
  DM             da;

  PetscFunctionBeginUser;
  ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
  ierr = DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
  ierr = DMGetApplicationContext(da,&x);CHKERRQ(ierr);
  h    = 2.0*PETSC_PI/((mx));
  ierr = VecCopy(x,b);CHKERRQ(ierr);
  ierr = VecScale(b,h);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

// assemble the right-hand side of the system for the Lagrangian forces
PetscErrorCode RigidKinematicsSolver::assembleRHSForces()
{
    PetscErrorCode ierr;

    PetscFunctionBeginUser;

    ierr = PetscLogStagePush(stageRHSForces); CHKERRQ(ierr);

    // rhsf = UB - E u^{**}
    ierr = MatMult(E, solution->UGlobal, rhsf); CHKERRQ(ierr);
    ierr = VecScale(rhsf, -1.0); CHKERRQ(ierr);
    ierr = VecAYPX(rhsf, 1.0, UB); CHKERRQ(ierr);

    ierr = PetscLogStagePop(); CHKERRQ(ierr);

    PetscFunctionReturn(0);
}  // assembleRHSForces

PetscErrorCode PCApply_Noise(PC pc,Vec xin,Vec xout)
{
  PetscErrorCode ierr;
  PCNoise_Ctx    *ctx;
  PetscReal      nrmin, nrmnoise;

  PetscFunctionBeginUser;
  ierr = PCShellGetContext(pc,(void**)&ctx);CHKERRQ(ierr);

  /* xout is ||xin|| * ctx->eta*  f, where f is a pseudorandom unit vector
    (Note that this should always be combined additively with another PC) */
  ierr = VecSetRandom(xout,ctx->random);CHKERRQ(ierr);
  ierr = VecNorm(xin,NORM_2,&nrmin);CHKERRQ(ierr);
  ierr = VecNorm(xout,NORM_2,&nrmnoise);CHKERRQ(ierr);
  ierr = VecScale(xout,ctx->eta*(nrmin/nrmnoise));CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode MatMult_SchurComplement(Mat N,Vec x,Vec y)
{
  Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
  PetscErrorCode      ierr;

  PetscFunctionBegin;
  if (!Na->work1) {ierr = MatCreateVecs(Na->A,&Na->work1,NULL);CHKERRQ(ierr);}
  if (!Na->work2) {ierr = MatCreateVecs(Na->A,&Na->work2,NULL);CHKERRQ(ierr);}
  ierr = MatMult(Na->B,x,Na->work1);CHKERRQ(ierr);
  ierr = KSPSolve(Na->ksp,Na->work1,Na->work2);CHKERRQ(ierr);
  ierr = MatMult(Na->C,Na->work2,y);CHKERRQ(ierr);
  ierr = VecScale(y,-1.0);CHKERRQ(ierr);
  if (Na->D) {
    ierr = MatMultAdd(Na->D,x,y,y);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

static PetscErrorCode TronGradientProjections(Tao tao,TAO_TRON *tron)
{
  PetscErrorCode                 ierr;
  PetscInt                       i;
  TaoLineSearchConvergedReason ls_reason;
  PetscReal                      actred=-1.0,actred_max=0.0;
  PetscReal                      f_new;
  /*
     The gradient and function value passed into and out of this
     routine should be current and correct.

     The free, active, and binding variables should be already identified
  */
  PetscFunctionBegin;
  if (tron->Free_Local) {
    ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr);
  }
  ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);CHKERRQ(ierr);

  for (i=0;i<tron->maxgpits;i++){

    if ( -actred <= (tron->pg_ftol)*actred_max) break;

    tron->gp_iterates++; tron->total_gp_its++;
    f_new=tron->f;

    ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr);
    ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
    ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize);CHKERRQ(ierr);
    ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection,
                              &tron->pgstepsize, &ls_reason);CHKERRQ(ierr);
    ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);


    /* Update the iterate */
    actred = f_new - tron->f;
    actred_max = PetscMax(actred_max,-(f_new - tron->f));
    tron->f = f_new;
    if (tron->Free_Local) {
      ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr);
    }
    ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);CHKERRQ(ierr);
  }

  PetscFunctionReturn(0);
}

void FractalProjector::
new_vertex(Vertex* vrt, Edge* parent)
{
//	set the vertex to the center by calling the parents method
	RefinementCallbackLinear<APosition>::new_vertex(vrt, parent);

//	calculate the normal of the edge
	vector3 n;
	CalculateNormal(n, *m_pGrid, parent, m_aaPos);
//	first scale it to the same length as the edge, then scale it with the
//	given scaleFac
	number len = VecDistance(m_aaPos[parent->vertex(0)], m_aaPos[parent->vertex(1)]);
	VecScale (n, n, m_scaleFac * len);

//	offset the vertex by the calculated normal
	VecAdd(m_aaPos[vrt], m_aaPos[vrt], n);
}

PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
{
  Mat_Shell      *shell = (Mat_Shell*)A->data;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = (*shell->getdiagonal)(A,v);CHKERRQ(ierr);
  ierr = VecScale(v,shell->vscale);CHKERRQ(ierr);
  if (shell->dshift) {
    ierr = VecPointwiseMult(v,v,shell->dshift);CHKERRQ(ierr);
  } else {
    ierr = VecShift(v,shell->vshift);CHKERRQ(ierr);
  }
  if (shell->left)  {ierr = VecPointwiseMult(v,v,shell->left);CHKERRQ(ierr);}
  if (shell->right) {ierr = VecPointwiseMult(v,v,shell->right);CHKERRQ(ierr);}
  PetscFunctionReturn(0);
}

Mapping *
LinearMappingCreate(
        Vector          *center,
        Vector          *vaxis,
        Vector          *uaxis)
{
	Mapping *res;
	RSMatrix m;
	Vector n;

	res = (Mapping *)Malloc(sizeof(Mapping));
	res->flags = OBJSPACE;
	res->method= LinearMapping;

	if (center)
		res->center = *center;
	else
		res->center.x = res->center.y = res->center.z = 0.;

	if (uaxis && vaxis) {
		VecCross(uaxis, vaxis, &n);
		/* this is wrong, since uaxis and vaxis
		 * give U and V in world space, and we
		 * need the inverse.
		 */
		ArbitraryMatrix(
			uaxis->x, uaxis->y, uaxis->z,
			vaxis->x, vaxis->y, vaxis->z,
			n.x, n.y, n.z,
			res->center.x, res->center.y, res->center.z,
			&m);
		MatrixInvert(&m, &res->m);
		res->uaxis = *uaxis;
		res->vaxis = *vaxis;
		VecNormalize(&res->uaxis);
		VecNormalize(&res->vaxis);
	} else {
		VecScale(-1., res->center, &n);
		TranslationMatrix(n.x, n.y, n.z, &res->m);
		res->uaxis.x = res->vaxis.y = 1.;
		res->uaxis.y = res->uaxis.z = res->vaxis.x =
			res->vaxis.z = 0.;
	}
	return res;
}

PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
{
  Mat_Shell *shell = (Mat_Shell*)Y->data;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  if (shell->left || shell->right || shell->dshift) {
    if (!shell->dshift) {
      if (!shell->dshift_owned) {ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift_owned);CHKERRQ(ierr);}
      shell->dshift = shell->dshift_owned;
      ierr = VecSet(shell->dshift,shell->vshift+a);CHKERRQ(ierr);
    } else {ierr = VecScale(shell->dshift,a);CHKERRQ(ierr);}
    if (shell->left)  {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);CHKERRQ(ierr);}
    if (shell->right) {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);CHKERRQ(ierr);}
  } else shell->vshift += a;
  ierr = MatShellUseScaledMethods(Y);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode MatMultTranspose_SubMatrix(Mat N,Vec x,Vec y)
{
  Mat_SubMatrix  *Na = (Mat_SubMatrix*)N->data;
  Vec            xx  = 0;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PreScaleLeft(N,x,&xx);CHKERRQ(ierr);
  ierr = VecZeroEntries(Na->lwork);CHKERRQ(ierr);
  ierr = VecScatterBegin(Na->lrestrict,xx,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
  ierr = VecScatterEnd  (Na->lrestrict,xx,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
  ierr = MatMultTranspose(Na->A,Na->lwork,Na->rwork);CHKERRQ(ierr);
  ierr = VecScatterBegin(Na->rprolong,Na->rwork,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
  ierr = VecScatterEnd  (Na->rprolong,Na->rwork,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
  ierr = PostScaleRight(N,y);CHKERRQ(ierr);
  ierr = VecScale(y,Na->scale);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode GPCGGradProjections(Tao tao)
{
  PetscErrorCode                 ierr;
  TAO_GPCG                       *gpcg = (TAO_GPCG *)tao->data;
  PetscInt                       i;
  PetscReal                      actred=-1.0,actred_max=0.0, gAg,gtg=gpcg->gnorm,alpha;
  PetscReal                      f_new,gdx,stepsize;
  Vec                            DX=tao->stepdirection,XL=tao->XL,XU=tao->XU,Work=gpcg->Work;
  Vec                            X=tao->solution,G=tao->gradient;
  TaoLineSearchConvergedReason lsflag=TAOLINESEARCH_CONTINUE_ITERATING;

  /*
     The free, active, and binding variables should be already identified
  */
  PetscFunctionBegin;
  for (i=0;i<gpcg->maxgpits;i++){
    if ( -actred <= (gpcg->pg_ftol)*actred_max) break;
    ierr = VecBoundGradientProjection(G,X,XL,XU,DX);CHKERRQ(ierr);
    ierr = VecScale(DX,-1.0);CHKERRQ(ierr);
    ierr = VecDot(DX,G,&gdx);CHKERRQ(ierr);

    ierr = MatMult(tao->hessian,DX,Work);CHKERRQ(ierr);
    ierr = VecDot(DX,Work,&gAg);CHKERRQ(ierr);

    gpcg->gp_iterates++;
    gpcg->total_gp_its++;

    gtg=-gdx;
    alpha = PetscAbsReal(gtg/gAg);
    ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,alpha);CHKERRQ(ierr);
    f_new=gpcg->f;
    ierr = TaoLineSearchApply(tao->linesearch,X,&f_new,G,DX,&stepsize,&lsflag);CHKERRQ(ierr);

    /* Update the iterate */
    actred = f_new - gpcg->f;
    actred_max = PetscMax(actred_max,-(f_new - gpcg->f));
    gpcg->f = f_new;
    ierr = ISDestroy(&gpcg->Free_Local);CHKERRQ(ierr);
    ierr = VecWhichBetween(XL,X,XU,&gpcg->Free_Local);CHKERRQ(ierr);
  }

  gpcg->gnorm=gtg;
  PetscFunctionReturn(0);
} /* End gradient projections */

PetscErrorCode Update_q(Vec q,Vec u,Mat M_0,AppCtx *user)
{
  PetscErrorCode ierr;
  PetscScalar    *q_arr,*w_arr;
  PetscInt       i,n;

  PetscFunctionBeginUser;
  ierr = PetscLogEventBegin(event_update_q,0,0,0,0);CHKERRQ(ierr);
  ierr = MatMult(M_0,u,user->work1);CHKERRQ(ierr);
  ierr = VecScale(user->work1,-1.0);CHKERRQ(ierr);
  ierr = VecGetLocalSize(u,&n);CHKERRQ(ierr);
  ierr = VecGetArray(q,&q_arr);CHKERRQ(ierr);
  ierr = VecGetArray(user->work1,&w_arr);CHKERRQ(ierr);
  for (i=0; i<n; i++) q_arr[2*i]=q_arr[2*i+1] = w_arr[i];
  ierr = VecRestoreArray(q,&q_arr);CHKERRQ(ierr);
  ierr = VecRestoreArray(user->work1,&w_arr);CHKERRQ(ierr);
  ierr = PetscLogEventEnd(event_update_q,0,0,0,0);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode VecScale_MultiVec(Vec x, PetscScalar alpha)
{
#if !defined(NDEBUG)
    TBOX_ASSERT(x);
#endif
    Vec_MultiVec* mx = static_cast<Vec_MultiVec*>(x->data);
#if !defined(NDEBUG)
    TBOX_ASSERT(mx);
#endif
    PetscErrorCode ierr;
    for (PetscInt k = 0; k < mx->n; ++k)
    {
        ierr = VecScale(mx->array[k], alpha);
        CHKERRQ(ierr);
    }
    ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(x));
    CHKERRQ(ierr);
    PetscFunctionReturn(0);
} // VecScale_MultiVec

static PetscErrorCode TaoBQNLSComputeStep(Tao tao, PetscBool shift, KSPConvergedReason *ksp_reason, PetscInt *step_type)
{
  TAO_BNK        *bnk = (TAO_BNK *)tao->data;
  TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
  PetscErrorCode ierr;
  PetscInt       nupdates;

  PetscFunctionBegin;
  ierr = MatSolve(bqnk->B, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
  ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
  ierr = TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection);CHKERRQ(ierr);
  *ksp_reason = KSP_CONVERGED_ATOL;
  ierr = MatLMVMGetUpdateCount(bqnk->B, &nupdates);CHKERRQ(ierr);
  if (nupdates == 0) {
    *step_type = BNK_SCALED_GRADIENT;
  } else {
    *step_type = BNK_BFGS;
  }
  PetscFunctionReturn(0);
}

void ConvectionDiffusionFE<TDomain>::
lin_def_vector_source(const LocalVector& u,
                           std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
                           const size_t nip)
{
//	request geometry
	const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);

//	loop integration points
	for(size_t ip = 0; ip < geo.num_ip(); ++ip)
	{
	//	loop test spaces
		for(size_t i = 0; i < geo.num_sh(); ++i)
		{
		//	add to local defect
			VecScale(vvvLinDef[ip][_C_][i], geo.global_grad(ip, i),
			         	 	 	 	 	 	 geo.weight(ip));
		}
	}
}