C++ VecAYPX函数代码示例

本文整理汇总了C++中VecAYPX函数的典型用法代码示例。如果您正苦于以下问题：C++ VecAYPX函数的具体用法？C++ VecAYPX怎么用？C++ VecAYPX使用的例子？那么, 这里精选的函数代码示例或许可以为您提供帮助。

在下文中一共展示了VecAYPX函数的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。

示例1: KSPBuildSolutionDefault

/*
   KSPBuildResidualDefault - Default code to compute the residual.

   Input Parameters:
.  ksp - iterative context
.  t   - pointer to temporary vector
.  v   - pointer to user vector

   Output Parameter:
.  V - pointer to a vector containing the residual

   Level: advanced

   Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations

.keywords:  KSP, build, residual, default

.seealso: KSPBuildSolutionDefault()
*/
PetscErrorCode KSPBuildResidualDefault(KSP ksp,Vec t,Vec v,Vec *V)
{
  PetscErrorCode ierr;
  Mat            Amat,Pmat;

  PetscFunctionBegin;
  if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);}
  ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr);
  ierr = KSPBuildSolution(ksp,t,NULL);CHKERRQ(ierr);
  ierr = KSP_MatMult(ksp,Amat,t,v);CHKERRQ(ierr);
  ierr = VecAYPX(v,-1.0,ksp->vec_rhs);CHKERRQ(ierr);
  *V   = v;
  PetscFunctionReturn(0);
}

开发者ID:mchandra，项目名称:petsc，代码行数:33，代码来源:iterativ.c

示例2: IFunction_Hull1972A1

PetscErrorCode IFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
{
  PetscErrorCode ierr;
  PetscScalar    *y,*f;

  PetscFunctionBegin;
  ierr = VecGetArray(Y,&y);CHKERRQ(ierr);
  ierr = VecGetArray(F,&f);CHKERRQ(ierr);
  f[0] = -y[0];
  ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr);
  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
  /* Left hand side = ydot - f(y) */
  ierr = VecAYPX(F,-1.0,Ydot);
  PetscFunctionReturn(0);
}

开发者ID:00liujj，项目名称:petsc，代码行数:15，代码来源:ex31.c

示例3: MonitorError

static PetscErrorCode MonitorError(TS ts,PetscInt step,PetscReal t,Vec x,void *ctx)
{
  PetscErrorCode ierr;
  MonitorCtx     *mon = (MonitorCtx*)ctx;
  PetscReal      h,nrm_x,nrm_exact,nrm_diff;

  PetscFunctionBeginUser;
  if (!mon->problem->solution) PetscFunctionReturn(0);
  ierr = (*mon->problem->solution)(t,mon->x,mon->problem->data);CHKERRQ(ierr);
  ierr = VecNorm(x,NORM_2,&nrm_x);CHKERRQ(ierr);
  ierr = VecNorm(mon->x,NORM_2,&nrm_exact);CHKERRQ(ierr);
  ierr = VecAYPX(mon->x,-1,x);CHKERRQ(ierr);
  ierr = VecNorm(mon->x,NORM_2,&nrm_diff);CHKERRQ(ierr);
  ierr = TSGetTimeStep(ts,&h);CHKERRQ(ierr);
  ierr = PetscPrintf(mon->comm,"step %4D t=%12.8e h=% 8.2e  |x|=%9.2e  |x_e|=%9.2e  |x-x_e|=%9.2e\n",step,t,h,nrm_x,nrm_exact,nrm_diff);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

开发者ID:ZJLi2013，项目名称:petsc，代码行数:17，代码来源:ex8.c

示例4: PetscLogStagePush

// 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

开发者ID:barbagroup，项目名称:PetIBM，代码行数:18，代码来源:rigidkinematics.cpp

示例5: _GetResidual

/* these functions were associated with the PETScMatrixSolver class, before it was
 * depreciated */
Vec _GetResidual( MGSolver_PETScData* mgData ) {
	if( mgData->expiredResidual ) {
		VecDuplicate( mgData->curSolution, &mgData->residual );	
		VecSetFromOptions( mgData->residual );
#if( PETSC_VERSION_MAJOR <= 2 && PETSC_VERSION_MINOR >= 3 && PETSC_VERSION_SUBMINOR >= 3 )
		VecSetOption( mgData->residual, VEC_IGNORE_NEGATIVE_INDICES );
#elif( PETSC_VERSION_MAJOR >= 3 )
		VecSetOption( mgData->residual, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE );
#endif
		MatMult( mgData->matrix, mgData->curSolution, mgData->residual );
		VecAYPX( mgData->residual, -1.0, mgData->curRHS );

		mgData->expiredResidual = False;
	}

	return mgData->residual;
}

开发者ID:OlympusMonds，项目名称:EarthByte_Underworld，代码行数:19，代码来源:MultigridSolver.c

示例6: IFunction_Hull1972C34

PetscErrorCode IFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
{
  PetscErrorCode ierr;
  PetscScalar    *y,*f;
  PetscInt       N,i;

  PetscFunctionBegin;
  ierr = VecGetSize (Y,&N);CHKERRQ(ierr);
  ierr = VecGetArray(Y,&y);CHKERRQ(ierr);
  ierr = VecGetArray(F,&f);CHKERRQ(ierr);
  f[0] = -2.0*y[0] + y[1];
  for (i = 1; i < N-1; i++) {
    f[i] = y[i-1] - 2.0*y[i] + y[i+1];
  }
  f[N-1] = y[N-2] - 2.0*y[N-1];
  ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr);
  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
  /* Left hand side = ydot - f(y) */
  ierr = VecAYPX(F,-1.0,Ydot);
  PetscFunctionReturn(0);
}

开发者ID:00liujj，项目名称:petsc，代码行数:21，代码来源:ex31.c

示例7: MatMultAdd_SchurComplement

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

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

开发者ID:erdc-cm，项目名称:petsc-dev，代码行数:22，代码来源:schurm.c

示例8: Metos3DTimeStepFunction

PetscErrorCode
Metos3DTimeStepFunction(SNES snes, Vec ynBD, Vec fnBD, void *ctx)
{
    Metos3D     *metos3d    = (Metos3D*)ctx;
    // geometry
    PetscInt    nvec        = metos3d->vectorLength;
    // bgc
    PetscInt    ntracer     = metos3d->tracerCount;
    // load
    PetscInt    nvecloc     = metos3d->vectorLengthLocal;
    // parameter
    PetscInt    nparam      = metos3d->parameterCount;
    PetscReal   *u0         = metos3d->u0;
    // work vars
    PetscInt    itracer;
    Vec         *yin, *yinold, *yout;
    PetscFunctionBegin;
    // create work vectors
    Metos3DUtilVecCreateAndSetValue(metos3d, ntracer, nvec, nvecloc, &yin, 0.0);
    Metos3DUtilVecCreateAndSetValue(metos3d, ntracer, nvec, nvecloc, &yinold, 0.0);
    Metos3DUtilVecCreateAndSetValue(metos3d, ntracer, nvec, nvecloc, &yout, 0.0);
    // yin      = ynBD
    // yinold   = ynBD
    // yout     = Phi(yin)
    // yout     = - yout + yinold;
    // fnBD     = yout
    Metos3DUtilVecCopyDiagonalToSeparate(metos3d, ntracer, nvecloc, &ynBD, yin);
    Metos3DUtilVecCopyDiagonalToSeparate(metos3d, ntracer, nvecloc, &ynBD, yinold);
    Metos3DTimeStepPhi(metos3d, yin, yout, nparam, u0);
    for (itracer = 0; itracer < ntracer; itracer++) VecAYPX(yout[itracer], -1.0, yinold[itracer]);
    Metos3DUtilVecCopySeparateToDiagonal(metos3d, ntracer, nvecloc, yout, &fnBD);
    // free work vectors
    VecDestroyVecs(ntracer, &yin);
    VecDestroyVecs(ntracer, &yinold);
    VecDestroyVecs(ntracer, &yout);
    // debug
    Metos3DDebug(metos3d, kDebugLevel, "Metos3DTimeStepFunction\n");
    PetscFunctionReturn(0);
}

开发者ID:neeljp，项目名称:simpack，代码行数:39，代码来源:metos3d_timestep.c

示例9: VecAYPX_MultiVec

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

开发者ID:mesnardo，项目名称:IBAMR，代码行数:23，代码来源:PETScMultiVec.cpp

示例10: SNESSolve

/*@
   SNESApplyNPC - Calls SNESSolve() on preconditioner for the SNES

   Collective on SNES

   Input Parameters:
+  snes - the SNES context
.  x - input vector
-  f - optional; the function evaluation on x

   Output Parameter:
.  y - function vector, as set by SNESSetFunction()

   Notes:
   SNESComputeFunction() should be called on x before SNESApplyNPC() is called, as it is
   with SNESComuteJacobian().

   Level: developer

.keywords: SNES, nonlinear, compute, function

.seealso: SNESGetNPC(),SNESSetNPC(),SNESComputeFunction()
@*/
PetscErrorCode  SNESApplyNPC(SNES snes,Vec x,Vec f,Vec y)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
  PetscValidHeaderSpecific(x,VEC_CLASSID,2);
  PetscValidHeaderSpecific(y,VEC_CLASSID,3);
  PetscCheckSameComm(snes,1,x,2);
  PetscCheckSameComm(snes,1,y,3);
  ierr = VecValidValues(x,2,PETSC_TRUE);CHKERRQ(ierr);
  if (snes->pc) {
    if (f) {
      ierr = SNESSetInitialFunction(snes->pc,f);CHKERRQ(ierr);
    }
    ierr = VecCopy(x,y);CHKERRQ(ierr);
    ierr = PetscLogEventBegin(SNES_NPCSolve,snes->pc,x,y,0);CHKERRQ(ierr);
    ierr = SNESSolve(snes->pc,snes->vec_rhs,y);CHKERRQ(ierr);
    ierr = PetscLogEventEnd(SNES_NPCSolve,snes->pc,x,y,0);CHKERRQ(ierr);
    ierr = VecAYPX(y,-1.0,x);CHKERRQ(ierr);
    PetscFunctionReturn(0);
  }
  PetscFunctionReturn(0);
}

开发者ID:lw4992，项目名称:petsc，代码行数:47，代码来源:snespc.c

示例11: KSPSolve_GCR

PetscErrorCode KSPSolve_GCR( KSP ksp )
{
  KSP_GCR        *ctx = (KSP_GCR*)ksp->data;
  PetscErrorCode ierr;
  Mat            A, B;
  Vec            r,b,x;
  PetscReal      norm_r;

  PetscFunctionBegin;
  ierr = KSPGetOperators( ksp, &A, &B, PETSC_NULL );CHKERRQ(ierr);
  x = ksp->vec_sol;
  b = ksp->vec_rhs;
  r = ctx->R;

  /* compute initial residual */
  ierr = MatMult( A, x, r );CHKERRQ(ierr);
  ierr = VecAYPX( r, -1.0, b );CHKERRQ(ierr); /* r = b - A x  */
  ierr = VecNorm( r, NORM_2, &norm_r );CHKERRQ(ierr);

  ksp->its = 0;
  ksp->rnorm0 = norm_r;

  KSPLogResidualHistory(ksp,ksp->rnorm0);
  ierr = KSPMonitor(ksp,ksp->its,ksp->rnorm0);CHKERRQ(ierr);
  ierr = (*ksp->converged)(ksp,ksp->its,ksp->rnorm0,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  if (ksp->reason) PetscFunctionReturn(0);

  do {
    ierr = KSPSolve_GCR_cycle( ksp );CHKERRQ(ierr);
    if (ksp->reason) break; /* catch case when convergence occurs inside the cycle */
  } while( ksp->its < ksp->max_it );CHKERRQ(ierr);
  if ( ksp->its >= ksp->max_it) {
    ksp->reason = KSP_DIVERGED_ITS;
  }
  PetscFunctionReturn(0);
}

开发者ID:erdc-cm，项目名称:petsc-dev，代码行数:36，代码来源:gcr.c

示例12: KSPSolve_CG

PetscErrorCode  KSPSolve_CG(KSP ksp)
{
  PetscErrorCode ierr;
  PetscInt       i,stored_max_it,eigs;
  PetscScalar    dpi = 0.0,a = 1.0,beta,betaold = 1.0,b = 0,*e = 0,*d = 0,delta,dpiold;
  PetscReal      dp  = 0.0;
  Vec            X,B,Z,R,P,S,W;
  KSP_CG         *cg;
  Mat            Amat,Pmat;
  PetscBool      diagonalscale;

  PetscFunctionBegin;
  ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr);
  if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);

  cg            = (KSP_CG*)ksp->data;
  eigs          = ksp->calc_sings;
  stored_max_it = ksp->max_it;
  X             = ksp->vec_sol;
  B             = ksp->vec_rhs;
  R             = ksp->work[0];
  Z             = ksp->work[1];
  P             = ksp->work[2];
  if (cg->singlereduction) {
    S = ksp->work[3];
    W = ksp->work[4];
  } else {
    S = 0;                      /* unused */
    W = Z;
  }

#define VecXDot(x,y,a) (((cg->type) == (KSP_CG_HERMITIAN)) ? VecDot(x,y,a) : VecTDot(x,y,a))

  if (eigs) {e = cg->e; d = cg->d; e[0] = 0.0; }
  ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr);

  ksp->its = 0;
  if (!ksp->guess_zero) {
    ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr);            /*     r <- b - Ax     */
    ierr = VecAYPX(R,-1.0,B);CHKERRQ(ierr);
  } else {
    ierr = VecCopy(B,R);CHKERRQ(ierr);                         /*     r <- b (x is 0) */
  }

  switch (ksp->normtype) {
  case KSP_NORM_PRECONDITIONED:
    ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr);                   /*     z <- Br         */
    ierr = VecNorm(Z,NORM_2,&dp);CHKERRQ(ierr);                /*    dp <- z'*z = e'*A'*B'*B*A'*e'     */
    break;
  case KSP_NORM_UNPRECONDITIONED:
    ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr);                /*    dp <- r'*r = e'*A'*A*e            */
    break;
  case KSP_NORM_NATURAL:
    ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr);                   /*     z <- Br         */
    if (cg->singlereduction) {
      ierr = KSP_MatMult(ksp,Amat,Z,S);CHKERRQ(ierr);
      ierr = VecXDot(Z,S,&delta);CHKERRQ(ierr);
    }
    ierr = VecXDot(Z,R,&beta);CHKERRQ(ierr);                     /*  beta <- z'*r       */
    if (PetscIsInfOrNanScalar(beta)) {
      if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf inner product");
      else {
        ksp->reason = KSP_DIVERGED_NANORINF;
        PetscFunctionReturn(0);
      }
    }
    dp = PetscSqrtReal(PetscAbsScalar(beta));                           /*    dp <- r'*z = r'*B*r = e'*A'*B*A*e */
    break;
  case KSP_NORM_NONE:
    dp = 0.0;
    break;
  default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
  }
  ierr       = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
  ierr       = KSPMonitor(ksp,0,dp);CHKERRQ(ierr);
  ksp->rnorm = dp;

  ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);      /* test for convergence */
  if (ksp->reason) PetscFunctionReturn(0);

  if (ksp->normtype != KSP_NORM_PRECONDITIONED && (ksp->normtype != KSP_NORM_NATURAL)) {
    ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr);                   /*     z <- Br         */
  }
  if (ksp->normtype != KSP_NORM_NATURAL) {
    if (cg->singlereduction) {
      ierr = KSP_MatMult(ksp,Amat,Z,S);CHKERRQ(ierr);
      ierr = VecXDot(Z,S,&delta);CHKERRQ(ierr);
    }
    ierr = VecXDot(Z,R,&beta);CHKERRQ(ierr);         /*  beta <- z'*r       */
    if (PetscIsInfOrNanScalar(beta)) {
      if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf inner product");
      else {
        ksp->reason = KSP_DIVERGED_NANORINF;
        PetscFunctionReturn(0);
      }
    }
  }

  i = 0;
  do {
//.........这里部分代码省略.........

开发者ID:PeiLiu90，项目名称:petsc，代码行数:101，代码来源:cg.c

示例13: KSPSolve_PIPECG

PetscErrorCode  KSPSolve_PIPECG(KSP ksp)
{
  PetscErrorCode ierr;
  PetscInt       i;
  PetscScalar    alpha = 0.0,beta = 0.0,gamma = 0.0,gammaold = 0.0,delta = 0.0;
  PetscReal      dp    = 0.0;
  Vec            X,B,Z,P,W,Q,U,M,N,R,S;
  Mat            Amat,Pmat;
  PetscBool      diagonalscale;

  PetscFunctionBegin;
  ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr);
  if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);

  X = ksp->vec_sol;
  B = ksp->vec_rhs;
  M = ksp->work[0];
  Z = ksp->work[1];
  P = ksp->work[2];
  N = ksp->work[3];
  W = ksp->work[4];
  Q = ksp->work[5];
  U = ksp->work[6];
  R = ksp->work[7];
  S = ksp->work[8];

  ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr);

  ksp->its = 0;
  if (!ksp->guess_zero) {
    ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr);            /*     r <- b - Ax     */
    ierr = VecAYPX(R,-1.0,B);CHKERRQ(ierr);
  } else {
    ierr = VecCopy(B,R);CHKERRQ(ierr);                         /*     r <- b (x is 0) */
  }

  ierr = KSP_PCApply(ksp,R,U);CHKERRQ(ierr);                   /*     u <- Br   */

  switch (ksp->normtype) {
  case KSP_NORM_PRECONDITIONED:
    ierr = VecNormBegin(U,NORM_2,&dp);CHKERRQ(ierr);                /*     dp <- u'*u = e'*A'*B'*B*A'*e'     */
    ierr = PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)U));CHKERRQ(ierr);
    ierr = KSP_MatMult(ksp,Amat,U,W);CHKERRQ(ierr);              /*     w <- Au   */
    ierr = VecNormEnd(U,NORM_2,&dp);CHKERRQ(ierr);
    break;
  case KSP_NORM_UNPRECONDITIONED:
    ierr = VecNormBegin(R,NORM_2,&dp);CHKERRQ(ierr);                /*     dp <- r'*r = e'*A'*A*e            */
    ierr = PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)R));CHKERRQ(ierr);
    ierr = KSP_MatMult(ksp,Amat,U,W);CHKERRQ(ierr);              /*     w <- Au   */
    ierr = VecNormEnd(R,NORM_2,&dp);CHKERRQ(ierr);
    break;
  case KSP_NORM_NATURAL:
    ierr = VecDotBegin(R,U,&gamma);CHKERRQ(ierr);                  /*     gamma <- u'*r       */
    ierr = PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)R));CHKERRQ(ierr);
    ierr = KSP_MatMult(ksp,Amat,U,W);CHKERRQ(ierr);              /*     w <- Au   */
    ierr = VecDotEnd(R,U,&gamma);CHKERRQ(ierr);
    KSPCheckDot(ksp,gamma);
    dp = PetscSqrtReal(PetscAbsScalar(gamma));                  /*     dp <- r'*u = r'*B*r = e'*A'*B*A*e */
    break;
  case KSP_NORM_NONE:
    ierr = KSP_MatMult(ksp,Amat,U,W);CHKERRQ(ierr);
    dp   = 0.0;
    break;
  default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
  }
  ierr       = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
  ierr       = KSPMonitor(ksp,0,dp);CHKERRQ(ierr);
  ksp->rnorm = dp;
  ierr       = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); /* test for convergence */
  if (ksp->reason) PetscFunctionReturn(0);

  i = 0;
  do {
    if (i > 0 && ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
      ierr = VecNormBegin(R,NORM_2,&dp);CHKERRQ(ierr);
    } else if (i > 0 && ksp->normtype == KSP_NORM_PRECONDITIONED) {
      ierr = VecNormBegin(U,NORM_2,&dp);CHKERRQ(ierr);
    }
    if (!(i == 0 && ksp->normtype == KSP_NORM_NATURAL)) {
      ierr = VecDotBegin(R,U,&gamma);CHKERRQ(ierr);
    }
    ierr = VecDotBegin(W,U,&delta);CHKERRQ(ierr);
    ierr = PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)R));CHKERRQ(ierr);

    ierr = KSP_PCApply(ksp,W,M);CHKERRQ(ierr);           /*   m <- Bw       */
    ierr = KSP_MatMult(ksp,Amat,M,N);CHKERRQ(ierr);      /*   n <- Am       */

    if (i > 0 && ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
      ierr = VecNormEnd(R,NORM_2,&dp);CHKERRQ(ierr);
    } else if (i > 0 && ksp->normtype == KSP_NORM_PRECONDITIONED) {
      ierr = VecNormEnd(U,NORM_2,&dp);CHKERRQ(ierr);
    }
    if (!(i == 0 && ksp->normtype == KSP_NORM_NATURAL)) {
      ierr = VecDotEnd(R,U,&gamma);CHKERRQ(ierr);
    }
    ierr = VecDotEnd(W,U,&delta);CHKERRQ(ierr);

    if (i > 0) {
      if (ksp->normtype == KSP_NORM_NATURAL) dp = PetscSqrtReal(PetscAbsScalar(gamma));
      else if (ksp->normtype == KSP_NORM_NONE) dp = 0.0;
//.........这里部分代码省略.........

开发者ID:pombredanne，项目名称:petsc，代码行数:101，代码来源:pipecg.c

示例14: SNESSolve_NCG


//.........这里部分代码省略.........
    }
    if (snes->nfuncs >= snes->max_funcs) {
      snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
      PetscFunctionReturn(0);
    }
    if (snes->domainerror) {
      snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
      PetscFunctionReturn(0);
    }
    ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr);
    /* Monitor convergence */
    ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
    snes->iter = i;
    snes->norm = fnorm;
    ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
    ierr       = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr);
    ierr       = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr);

    /* Test for convergence */
    ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
    if (snes->reason) PetscFunctionReturn(0);

    /* Call general purpose update function */
    if (snes->ops->update) {
      ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr);
    }
    if (snes->pc) {
      if (snes->functype == SNES_FUNCTION_PRECONDITIONED) {
        ierr = SNESApplyNPC(snes,X,NULL,dX);CHKERRQ(ierr);
        ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr);
        if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
          snes->reason = SNES_DIVERGED_INNER;
          PetscFunctionReturn(0);
        }
        ierr = VecCopy(dX,F);CHKERRQ(ierr);
      } else {
        ierr = SNESApplyNPC(snes,X,F,dX);CHKERRQ(ierr);
        ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr);
        if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
          snes->reason = SNES_DIVERGED_INNER;
          PetscFunctionReturn(0);
        }
      }
    } else {
      ierr = VecCopy(F,dX);CHKERRQ(ierr);
    }

    /* compute the conjugate direction lX = dX + beta*lX with beta = ((dX, dX) / (dX_old, dX_old) (Fletcher-Reeves update)*/
    switch (ncg->type) {
    case SNES_NCG_FR: /* Fletcher-Reeves */
      dXolddotdXold= dXdotdX;
      ierr         = VecDot(dX, dX, &dXdotdX);CHKERRQ(ierr);
      beta         = PetscRealPart(dXdotdX / dXolddotdXold);
      break;
    case SNES_NCG_PRP: /* Polak-Ribiere-Poylak */
      dXolddotdXold= dXdotdX;
      ierr         = VecDotBegin(dX, dX, &dXdotdXold);CHKERRQ(ierr);
      ierr         = VecDotBegin(dXold, dX, &dXdotdXold);CHKERRQ(ierr);
      ierr         = VecDotEnd(dX, dX, &dXdotdX);CHKERRQ(ierr);
      ierr         = VecDotEnd(dXold, dX, &dXdotdXold);CHKERRQ(ierr);
      beta         = PetscRealPart(((dXdotdX - dXdotdXold) / dXolddotdXold));
      if (beta < 0.0) beta = 0.0; /* restart */
      break;
    case SNES_NCG_HS: /* Hestenes-Stiefel */
      ierr = VecDotBegin(dX, dX, &dXdotdX);CHKERRQ(ierr);
      ierr = VecDotBegin(dX, dXold, &dXdotdXold);CHKERRQ(ierr);
      ierr = VecDotBegin(lX, dX, &lXdotdX);CHKERRQ(ierr);
      ierr = VecDotBegin(lX, dXold, &lXdotdXold);CHKERRQ(ierr);
      ierr = VecDotEnd(dX, dX, &dXdotdX);CHKERRQ(ierr);
      ierr = VecDotEnd(dX, dXold, &dXdotdXold);CHKERRQ(ierr);
      ierr = VecDotEnd(lX, dX, &lXdotdX);CHKERRQ(ierr);
      ierr = VecDotEnd(lX, dXold, &lXdotdXold);CHKERRQ(ierr);
      beta = PetscRealPart((dXdotdX - dXdotdXold) / (lXdotdX - lXdotdXold));
      break;
    case SNES_NCG_DY: /* Dai-Yuan */
      ierr = VecDotBegin(dX, dX, &dXdotdX);CHKERRQ(ierr);
      ierr = VecDotBegin(lX, dX, &lXdotdX);CHKERRQ(ierr);
      ierr = VecDotBegin(lX, dXold, &lXdotdXold);CHKERRQ(ierr);
      ierr = VecDotEnd(dX, dX, &dXdotdX);CHKERRQ(ierr);
      ierr = VecDotEnd(lX, dX, &lXdotdX);CHKERRQ(ierr);
      ierr = VecDotEnd(lX, dXold, &lXdotdXold);CHKERRQ(ierr);
      beta = PetscRealPart(dXdotdX / (lXdotdXold - lXdotdX));CHKERRQ(ierr);
      break;
    case SNES_NCG_CD: /* Conjugate Descent */
      ierr = VecDotBegin(dX, dX, &dXdotdX);CHKERRQ(ierr);
      ierr = VecDotBegin(lX, dXold, &lXdotdXold);CHKERRQ(ierr);
      ierr = VecDotEnd(dX, dX, &dXdotdX);CHKERRQ(ierr);
      ierr = VecDotEnd(lX, dXold, &lXdotdXold);CHKERRQ(ierr);
      beta = PetscRealPart(dXdotdX / lXdotdXold);CHKERRQ(ierr);
      break;
    }
    if (ncg->monitor) {
      ierr = PetscViewerASCIIPrintf(ncg->monitor, "beta = %e\n", (double)beta);CHKERRQ(ierr);
    }
    ierr = VecAYPX(lX, beta, dX);CHKERRQ(ierr);
  }
  ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);CHKERRQ(ierr);
  if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
  PetscFunctionReturn(0);
}

开发者ID:00liujj，项目名称:petsc，代码行数:101，代码来源:snesncg.c

示例15: KSPSolve_QCG


//.........这里部分代码省略.........
    ierr = VecDotRealPart(P,ASP,&ptasp);CHKERRQ(ierr);
    if (ptasp <= dzero) {

      /* Scaled negative curvature direction:  Compute a step so that
        ||w + step*p|| = delta and QS(w + step*p) is least */

      if (!i) {
        ierr = VecCopy(P,X);CHKERRQ(ierr);
        ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr);
        scal = pcgP->delta / xnorm;
        ierr = VecScale(X,scal);CHKERRQ(ierr);
      } else {
        /* Compute roots of quadratic */
        ierr = KSPQCGQuadraticRoots(W,P,pcgP->delta,&step1,&step2);CHKERRQ(ierr);
        ierr = VecDotRealPart(W,ASP,&wtasp);CHKERRQ(ierr);
        ierr = VecDotRealPart(BS,P,&bstp);CHKERRQ(ierr);
        ierr = VecCopy(W,X);CHKERRQ(ierr);
        q1   = step1*(bstp + wtasp + .5*step1*ptasp);
        q2   = step2*(bstp + wtasp + .5*step2*ptasp);
        if (q1 <= q2) {
          ierr = VecAXPY(X,step1,P);CHKERRQ(ierr);
        } else {
          ierr = VecAXPY(X,step2,P);CHKERRQ(ierr);
        }
      }
      pcgP->ltsnrm = pcgP->delta;                       /* convergence in direction of */
      ksp->reason  = KSP_CONVERGED_CG_NEG_CURVE;  /* negative curvature */
      if (!i) {
        ierr = PetscInfo1(ksp,"negative curvature: delta=%g\n",(double)pcgP->delta);CHKERRQ(ierr);
      } else {
        ierr = PetscInfo3(ksp,"negative curvature: step1=%g, step2=%g, delta=%g\n",(double)step1,(double)step2,(double)pcgP->delta);CHKERRQ(ierr);
      }

    } else {
      /* Compute step along p */
      step = rtr/ptasp;
      ierr = VecCopy(W,X);CHKERRQ(ierr);        /*  x = w  */
      ierr = VecAXPY(X,step,P);CHKERRQ(ierr);   /*  x <- step*p + x  */
      ierr = VecNorm(X,NORM_2,&pcgP->ltsnrm);CHKERRQ(ierr);

      if (pcgP->ltsnrm > pcgP->delta) {
        /* Since the trial iterate is outside the trust region,
            evaluate a constrained step along p so that
                    ||w + step*p|| = delta
          The positive step is always better in this case. */
        if (!i) {
          scal = pcgP->delta / pcgP->ltsnrm;
          ierr = VecScale(X,scal);CHKERRQ(ierr);
        } else {
          /* Compute roots of quadratic */
          ierr = KSPQCGQuadraticRoots(W,P,pcgP->delta,&step1,&step2);CHKERRQ(ierr);
          ierr = VecCopy(W,X);CHKERRQ(ierr);
          ierr = VecAXPY(X,step1,P);CHKERRQ(ierr);  /*  x <- step1*p + x  */
        }
        pcgP->ltsnrm = pcgP->delta;
        ksp->reason  = KSP_CONVERGED_CG_CONSTRAINED; /* convergence along constrained step */
        if (!i) {
          ierr = PetscInfo1(ksp,"constrained step: delta=%g\n",(double)pcgP->delta);CHKERRQ(ierr);
        } else {
          ierr = PetscInfo3(ksp,"constrained step: step1=%g, step2=%g, delta=%g\n",(double)step1,(double)step2,(double)pcgP->delta);CHKERRQ(ierr);
        }

      } else {
        /* Evaluate the current step */
        ierr = VecCopy(X,W);CHKERRQ(ierr);  /* update interior iterate */
        ierr = VecAXPY(R,-step,ASP);CHKERRQ(ierr); /* r <- -step*asp + r */
        ierr = VecNorm(R,NORM_2,&rnrm);CHKERRQ(ierr);

        ierr       = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
        ksp->rnorm = rnrm;
        ierr       = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
        ierr = KSPLogResidualHistory(ksp,rnrm);CHKERRQ(ierr);
        ierr = KSPMonitor(ksp,i+1,rnrm);CHKERRQ(ierr);
        ierr = (*ksp->converged)(ksp,i+1,rnrm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
        if (ksp->reason) {                 /* convergence for */
          ierr = PetscInfo3(ksp,"truncated step: step=%g, rnrm=%g, delta=%g\n",(double)PetscRealPart(step),(double)rnrm,(double)pcgP->delta);CHKERRQ(ierr);
        }
      }
    }
    if (ksp->reason) break;  /* Convergence has been attained */
    else {                   /* Compute a new AS-orthogonal direction */
      ierr = VecDot(R,R,&rntrn);CHKERRQ(ierr);
      beta = rntrn/rtr;
      ierr = VecAYPX(P,beta,R);CHKERRQ(ierr);  /*  p <- r + beta*p  */
      rtr  = PetscRealPart(rntrn);
    }
  }
  if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;

  /* Unscale x */
  ierr = VecCopy(X,WA2);CHKERRQ(ierr);
  ierr = PCApplySymmetricRight(pc,WA2,X);CHKERRQ(ierr);

  ierr = KSP_MatMult(ksp,Amat,X,WA);CHKERRQ(ierr);
  ierr = VecDotRealPart(B,X,&btx);CHKERRQ(ierr);
  ierr = VecDotRealPart(X,WA,&xtax);CHKERRQ(ierr);

  pcgP->quadratic = btx + .5*xtax;
  PetscFunctionReturn(0);
}

开发者ID:pombredanne，项目名称:petsc，代码行数:101，代码来源:qcg.c

注：本文中的VecAYPX函数示例由纯净天空整理自Github/MSDocs等开源代码及文档管理平台，相关代码片段筛选自各路编程大神贡献的开源项目，源码版权归原作者所有，传播和使用请参考对应项目的License；未经允许，请勿转载。