• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    公众号

C++ PetscAbsScalar函数代码示例

原作者: [db:作者] 来自: [db:来源] 收藏 邀请

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



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

示例1: KSPSolve_PIPEFCG_cycle

static PetscErrorCode KSPSolve_PIPEFCG_cycle(KSP ksp)
{
  PetscErrorCode ierr;
  PetscInt       i,j,k,idx,kdx,mi;
  KSP_PIPEFCG    *pipefcg;
  PetscScalar    alpha=0.0,gamma,*betas,*dots;
  PetscReal      dp=0.0, delta,*eta,*etas;
  Vec            B,R,Z,X,Qcurr,W,ZETAcurr,M,N,Pcurr,Scurr,*redux;
  Mat            Amat,Pmat;

  PetscFunctionBegin;

  /* We have not checked these routines for use with complex numbers. The inner products
     are likely not defined correctly for that case */
#if (defined(PETSC_USE_COMPLEX) && !defined(PETSC_SKIP_COMPLEX))
  SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"PIPEFGMRES has not been implemented for use with complex scalars");
#endif

#define VecXDot(x,y,a)         (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecDot       (x,y,a)   : VecTDot       (x,y,a))
#define VecXDotBegin(x,y,a)    (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecDotBegin  (x,y,a)   : VecTDotBegin  (x,y,a))
#define VecXDotEnd(x,y,a)      (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecDotEnd    (x,y,a)   : VecTDotEnd    (x,y,a))
#define VecMXDot(x,n,y,a)      (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecMDot      (x,n,y,a) : VecMTDot      (x,n,y,a))
#define VecMXDotBegin(x,n,y,a) (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecMDotBegin (x,n,y,a) : VecMTDotBegin (x,n,y,a))
#define VecMXDotEnd(x,n,y,a)   (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecMDotEnd   (x,n,y,a) : VecMTDotEnd   (x,n,y,a))

  pipefcg       = (KSP_PIPEFCG*)ksp->data;
  X             = ksp->vec_sol;
  B             = ksp->vec_rhs;
  R             = ksp->work[0];
  Z             = ksp->work[1];
  W             = ksp->work[2];
  M             = ksp->work[3];
  N             = ksp->work[4];

  redux = pipefcg->redux;
  dots  = pipefcg->dots;
  etas  = pipefcg->etas;
  betas = dots;        /* dots takes the result of all dot products of which the betas are a subset */

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

  /* Compute cycle initial residual */
  ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr);
  ierr = VecAYPX(R,-1.0,B);CHKERRQ(ierr);                   /* r <- b - Ax */
  ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr);                /* z <- Br     */

  Pcurr = pipefcg->Pvecs[0];
  Scurr = pipefcg->Svecs[0];
  Qcurr = pipefcg->Qvecs[0];
  ZETAcurr = pipefcg->ZETAvecs[0];
  ierr  = VecCopy(Z,Pcurr);CHKERRQ(ierr);
  ierr  = KSP_MatMult(ksp,Amat,Pcurr,Scurr);CHKERRQ(ierr);  /* S = Ap     */
  ierr  = VecCopy(Scurr,W);CHKERRQ(ierr);                   /* w = s = Az */

  /* Initial state of pipelining intermediates */
  redux[0] = R;
  redux[1] = W;
  ierr     = VecMXDotBegin(Z,2,redux,dots);CHKERRQ(ierr);
  ierr     = PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)Z));CHKERRQ(ierr); /* perform asynchronous reduction */
  ierr     = KSP_PCApply(ksp,W,M);CHKERRQ(ierr);            /* m = B(w) */
  ierr     = KSP_MatMult(ksp,Amat,M,N);CHKERRQ(ierr);       /* n = Am   */
  ierr     = VecCopy(M,Qcurr);CHKERRQ(ierr);                /* q = m    */
  ierr     = VecCopy(N,ZETAcurr);CHKERRQ(ierr);             /* zeta = n */
  ierr     = VecMXDotEnd(Z,2,redux,dots);CHKERRQ(ierr);
  gamma    = dots[0];
  delta    = PetscRealPart(dots[1]);
  etas[0]  = delta;
  alpha    = gamma/delta;

  i = 0;
  do {
    ksp->its++;

    /* Update X, R, Z, W */
    ierr = VecAXPY(X,+alpha,Pcurr);CHKERRQ(ierr);           /* x <- x + alpha * pi    */
    ierr = VecAXPY(R,-alpha,Scurr);CHKERRQ(ierr);           /* r <- r - alpha * si    */
    ierr = VecAXPY(Z,-alpha,Qcurr);CHKERRQ(ierr);           /* z <- z - alpha * qi    */
    ierr = VecAXPY(W,-alpha,ZETAcurr);CHKERRQ(ierr);        /* w <- w - alpha * zetai */

    /* Compute norm for convergence check */
    switch (ksp->normtype) {
      case KSP_NORM_PRECONDITIONED:
        ierr = VecNorm(Z,NORM_2,&dp);CHKERRQ(ierr);         /* dp <- sqrt(z'*z) = sqrt(e'*A'*B'*B*A*e) */
        break;
      case KSP_NORM_UNPRECONDITIONED:
        ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr);         /* dp <- sqrt(r'*r) = sqrt(e'*A'*A*e)      */
        break;
      case KSP_NORM_NATURAL:
        dp = PetscSqrtReal(PetscAbsScalar(gamma));          /* dp <- sqrt(r'*z) = sqrt(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]);
    }

    /* Check for convergence */
    ksp->rnorm = dp;
    KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
    ierr = KSPMonitor(ksp,ksp->its,dp);CHKERRQ(ierr);
//.........这里部分代码省略.........
开发者ID:tom-klotz,项目名称:petsc,代码行数:101,代码来源:pipefcg.c


示例2: PCGAMGProlongator_Classical


//.........这里部分代码省略.........
          if (PetscRealPart(rval[k]) > 0) {
            a_pos += rval[k];
          } else {
            a_neg += rval[k];
          }
          ntotal++;
        } else diag = rval[k];
      }
      ierr = MatRestoreRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);

      /* ghosted all connections */
      if (gA) {
        ierr = MatGetRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
        for (k = 0; k < ncols; k++) {
          if (PetscRealPart(rval[k]) > 0.) {
            a_pos += PetscRealPart(rval[k]);
          } else {
            a_neg += PetscRealPart(rval[k]);
          }
          ntotal++;
        }
        ierr = MatRestoreRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
      }

      if (g_neg == 0.) {
        alpha = 0.;
      } else {
        alpha = -a_neg/g_neg;
      }

      if (g_pos == 0.) {
        diag += a_pos;
        beta = 0.;
      } else {
        beta = -a_pos/g_pos;
      }
      if (diag == 0.) {
        invdiag = 0.;
      } else invdiag = 1. / diag;
      /* on */
      ierr = MatGetRow(lG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
      idx = 0;
      for (j = 0;j < ncols;j++) {
        col = icols[j];
        if (lcid[col] >= 0 && vcols[j] != 0.) {
          row_f = i + fs;
          row_c = lcid[col];
          /* set the values for on-processor ones */
          if (PetscRealPart(vcols[j]) < 0.) {
            pij = vcols[j]*alpha*invdiag;
          } else {
            pij = vcols[j]*beta*invdiag;
          }
          if (PetscAbsScalar(pij) != 0.) {
            pvals[idx] = pij;
            pcols[idx] = row_c;
            idx++;
          }
        }
      }
      ierr = MatRestoreRow(lG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
      /* off */
      if (gG) {
        ierr = MatGetRow(gG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
        for (j = 0; j < ncols; j++) {
          col = icols[j];
          if (gcid[col] >= 0 && vcols[j] != 0.) {
            row_f = i + fs;
            row_c = gcid[col];
            /* set the values for on-processor ones */
            if (PetscRealPart(vcols[j]) < 0.) {
              pij = vcols[j]*alpha*invdiag;
            } else {
              pij = vcols[j]*beta*invdiag;
            }
            if (PetscAbsScalar(pij) != 0.) {
              pvals[idx] = pij;
              pcols[idx] = row_c;
              idx++;
            }
          }
        }
        ierr = MatRestoreRow(gG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
      }
      ierr = MatSetValues(*P,1,&row_f,idx,pcols,pvals,INSERT_VALUES);CHKERRQ(ierr);
    }
  }

  ierr = MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  ierr = PetscFree(lsparse);CHKERRQ(ierr);
  ierr = PetscFree(gsparse);CHKERRQ(ierr);
  ierr = PetscFree(pcols);CHKERRQ(ierr);
  ierr = PetscFree(pvals);CHKERRQ(ierr);
  ierr = PetscFree(lcid);CHKERRQ(ierr);
  if (gG) {ierr = PetscFree(gcid);CHKERRQ(ierr);}

  PetscFunctionReturn(0);
}
开发者ID:hsahasra,项目名称:petsc-magma-dense-mat,代码行数:101,代码来源:classical.c


示例3: 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;
  MatStructure   pflag;
  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,&pflag);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;
//.........这里部分代码省略.........
开发者ID:hsahasra,项目名称:petsc-magma-dense-mat,代码行数:101,代码来源:cg.c


示例4: KSPLGMRESUpdateHessenberg

static PetscErrorCode KSPLGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool  hapend,PetscReal *res)
{
    PetscScalar   *hh,*cc,*ss,tt;
    PetscInt      j;
    KSP_LGMRES    *lgmres = (KSP_LGMRES *)(ksp->data);

    PetscFunctionBegin;
    hh  = HH(0,it);  /* pointer to beginning of column to update - so
                      incrementing hh "steps down" the (it+1)th col of HH*/
    cc  = CC(0);     /* beginning of cosine rotations */
    ss  = SS(0);     /* beginning of sine rotations */

    /* Apply all the previously computed plane rotations to the new column
       of the Hessenberg matrix */
    /* Note: this uses the rotation [conj(c)  s ; -s   c], c= cos(theta), s= sin(theta) */

    for (j=1; j<=it; j++) {
        tt  = *hh;
        *hh = PetscConj(*cc) * tt + *ss * *(hh+1);
        hh++;
        *hh = *cc++ * *hh - (*ss++ * tt);
        /* hh, cc, and ss have all been incremented one by end of loop */
    }

    /*
      compute the new plane rotation, and apply it to:
       1) the right-hand-side of the Hessenberg system (GRS)
          note: it affects GRS(it) and GRS(it+1)
       2) the new column of the Hessenberg matrix
          note: it affects HH(it,it) which is currently pointed to
          by hh and HH(it+1, it) (*(hh+1))
      thus obtaining the updated value of the residual...
    */

    /* compute new plane rotation */

    if (!hapend) {
        tt        = PetscSqrtScalar(PetscConj(*hh) * *hh + PetscConj(*(hh+1)) * *(hh+1));
        if (tt == 0.0) {
            ksp->reason = KSP_DIVERGED_NULL;
            PetscFunctionReturn(0);
        }
        *cc       = *hh / tt;   /* new cosine value */
        *ss       = *(hh+1) / tt;  /* new sine value */

        /* apply to 1) and 2) */
        *GRS(it+1) = - (*ss * *GRS(it));
        *GRS(it)   = PetscConj(*cc) * *GRS(it);
        *hh        = PetscConj(*cc) * *hh + *ss * *(hh+1);

        /* residual is the last element (it+1) of right-hand side! */
        *res      = PetscAbsScalar(*GRS(it+1));

    } else {
        /* happy breakdown: HH(it+1, it) = 0, therfore we don't need to apply
                  another rotation matrix (so RH doesn't change).  The new residual is
                  always the new sine term times the residual from last time (GRS(it)),
                  but now the new sine rotation would be zero...so the residual should
                  be zero...so we will multiply "zero" by the last residual.  This might
                  not be exactly what we want to do here -could just return "zero". */

        *res = 0.0;
    }
    PetscFunctionReturn(0);
}
开发者ID:erdc-cm,项目名称:petsc-dev,代码行数:65,代码来源:lgmres.c


示例5: FormFunctionLocal

PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
{
  AppCtx         *user = (AppCtx*)ptr;
  PetscErrorCode ierr;
  PetscInt       xints,xinte,yints,yinte,i,j;
  PetscReal      hx,hy,dhx,dhy,hxdhy,hydhx;
  PetscReal      grashof,prandtl,lid;
  PetscScalar    u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

  PetscFunctionBeginUser;
  grashof = user->grashof;
  prandtl = user->prandtl;
  lid     = user->lidvelocity;

  /*
     Define mesh intervals ratios for uniform grid.

     Note: FD formulae below are normalized by multiplying through by
     local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.


  */
  dhx   = (PetscReal)(info->mx-1);  dhy = (PetscReal)(info->my-1);
  hx    = 1.0/dhx;                   hy = 1.0/dhy;
  hxdhy = hx*dhy;                 hydhx = hy*dhx;

  xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;

  /* Test whether we are on the bottom edge of the global array */
  if (yints == 0) {
    j     = 0;
    yints = yints + 1;
    /* bottom edge */
    for (i=info->xs; i<info->xs+info->xm; i++) {
      f[j][i].u     = x[j][i].u;
      f[j][i].v     = x[j][i].v;
      f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
      f[j][i].temp  = x[j][i].temp-x[j+1][i].temp;
    }
  }

  /* Test whether we are on the top edge of the global array */
  if (yinte == info->my) {
    j     = info->my - 1;
    yinte = yinte - 1;
    /* top edge */
    for (i=info->xs; i<info->xs+info->xm; i++) {
      f[j][i].u     = x[j][i].u - lid;
      f[j][i].v     = x[j][i].v;
      f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
      f[j][i].temp  = x[j][i].temp-x[j-1][i].temp;
    }
  }

  /* Test whether we are on the left edge of the global array */
  if (xints == 0) {
    i     = 0;
    xints = xints + 1;
    /* left edge */
    for (j=info->ys; j<info->ys+info->ym; j++) {
      f[j][i].u     = x[j][i].u;
      f[j][i].v     = x[j][i].v;
      f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
      f[j][i].temp  = x[j][i].temp;
    }
  }

  /* Test whether we are on the right edge of the global array */
  if (xinte == info->mx) {
    i     = info->mx - 1;
    xinte = xinte - 1;
    /* right edge */
    for (j=info->ys; j<info->ys+info->ym; j++) {
      f[j][i].u     = x[j][i].u;
      f[j][i].v     = x[j][i].v;
      f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
      f[j][i].temp  = x[j][i].temp - (PetscReal)(grashof>0);
    }
  }

  /* Compute over the interior points */
  for (j=yints; j<yinte; j++) {
    for (i=xints; i<xinte; i++) {

      /*
       convective coefficients for upwinding
      */
      vx  = x[j][i].u; avx = PetscAbsScalar(vx);
      vxp = .5*(vx+avx); vxm = .5*(vx-avx);
      vy  = x[j][i].v; avy = PetscAbsScalar(vy);
      vyp = .5*(vy+avy); vym = .5*(vy-avy);

      /* U velocity */
      u         = x[j][i].u;
      uxx       = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
      uyy       = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
      f[j][i].u = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

      /* V velocity */
      u         = x[j][i].v;
//.........这里部分代码省略.........
开发者ID:wgapl,项目名称:petsc,代码行数:101,代码来源:ex19.c


示例6: SNESSetFunction

/*@C
   SNESComputeJacobianDefault - Computes the Jacobian using finite differences.

   Collective on SNES

   Input Parameters:
+  x1 - compute Jacobian at this point
-  ctx - application's function context, as set with SNESSetFunction()

   Output Parameters:
+  J - Jacobian matrix (not altered in this routine)
-  B - newly computed Jacobian matrix to use with preconditioner (generally the same as J)

   Options Database Key:
+  -snes_fd - Activates SNESComputeJacobianDefault()
.  -snes_test_err - Square root of function error tolerance, default square root of machine
                    epsilon (1.e-8 in double, 3.e-4 in single)
-  -mat_fd_type - Either wp or ds (see MATMFFD_WP or MATMFFD_DS)

   Notes:
   This routine is slow and expensive, and is not currently optimized
   to take advantage of sparsity in the problem.  Although
   SNESComputeJacobianDefault() is not recommended for general use
   in large-scale applications, It can be useful in checking the
   correctness of a user-provided Jacobian.

   An alternative routine that uses coloring to exploit matrix sparsity is
   SNESComputeJacobianDefaultColor().

   Level: intermediate

.keywords: SNES, finite differences, Jacobian

.seealso: SNESSetJacobian(), SNESComputeJacobianDefaultColor(), MatCreateSNESMF()
@*/
PetscErrorCode  SNESComputeJacobianDefault(SNES snes,Vec x1,Mat J,Mat B,void *ctx)
{
  Vec               j1a,j2a,x2;
  PetscErrorCode    ierr;
  PetscInt          i,N,start,end,j,value,root;
  PetscScalar       dx,*y,wscale;
  const PetscScalar *xx;
  PetscReal         amax,epsilon = PETSC_SQRT_MACHINE_EPSILON;
  PetscReal         dx_min = 1.e-16,dx_par = 1.e-1,unorm;
  MPI_Comm          comm;
  PetscErrorCode    (*eval_fct)(SNES,Vec,Vec)=0;
  PetscBool         assembled,use_wp = PETSC_TRUE,flg;
  const char        *list[2] = {"ds","wp"};
  PetscMPIInt       size;
  const PetscInt    *ranges;

  PetscFunctionBegin;
  ierr     = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_test_err",&epsilon,0);CHKERRQ(ierr);
  eval_fct = SNESComputeFunction;

  ierr = PetscObjectGetComm((PetscObject)x1,&comm);CHKERRQ(ierr);
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
  ierr = MatAssembled(B,&assembled);CHKERRQ(ierr);
  if (assembled) {
    ierr = MatZeroEntries(B);CHKERRQ(ierr);
  }
  if (!snes->nvwork) {
    snes->nvwork = 3;

    ierr = VecDuplicateVecs(x1,snes->nvwork,&snes->vwork);CHKERRQ(ierr);
    ierr = PetscLogObjectParents(snes,snes->nvwork,snes->vwork);CHKERRQ(ierr);
  }
  j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2];

  ierr = VecGetSize(x1,&N);CHKERRQ(ierr);
  ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr);
  ierr = (*eval_fct)(snes,x1,j1a);CHKERRQ(ierr);

  ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing options","SNES");CHKERRQ(ierr);
  ierr = PetscOptionsEList("-mat_fd_type","Algorithm to compute difference parameter","SNESComputeJacobianDefault",list,2,"wp",&value,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsEnd();CHKERRQ(ierr);
  if (flg && !value) use_wp = PETSC_FALSE;

  if (use_wp) {
    ierr = VecNorm(x1,NORM_2,&unorm);CHKERRQ(ierr);
  }
  /* Compute Jacobian approximation, 1 column at a time.
      x1 = current iterate, j1a = F(x1)
      x2 = perturbed iterate, j2a = F(x2)
   */
  for (i=0; i<N; i++) {
    ierr = VecCopy(x1,x2);CHKERRQ(ierr);
    if (i>= start && i<end) {
      ierr = VecGetArrayRead(x1,&xx);CHKERRQ(ierr);
      if (use_wp) dx = PetscSqrtReal(1.0 + unorm);
      else        dx = xx[i-start];
      ierr = VecRestoreArrayRead(x1,&xx);CHKERRQ(ierr);
      if (PetscAbsScalar(dx) < dx_min) dx = (PetscRealPart(dx) < 0. ? -1. : 1.) * dx_par;
      dx    *= epsilon;
      wscale = 1.0/dx;
      ierr   = VecSetValues(x2,1,&i,&dx,ADD_VALUES);CHKERRQ(ierr);
    } else {
      wscale = 0.0;
    }
    ierr = VecAssemblyBegin(x2);CHKERRQ(ierr);
//.........这里部分代码省略.........
开发者ID:haubentaucher,项目名称:petsc,代码行数:101,代码来源:snesj.c


示例7: KSPSolve_FCG

PetscErrorCode KSPSolve_FCG(KSP ksp)
{
  PetscErrorCode ierr;
  PetscInt       i,k,idx,mi;
  KSP_FCG        *fcg = (KSP_FCG*)ksp->data;
  PetscScalar    alpha=0.0,beta = 0.0,dpi,s;
  PetscReal      dp=0.0;
  Vec            B,R,Z,X,Pcurr,Ccurr;
  Mat            Amat,Pmat;
  PetscInt       eigs = ksp->calc_sings; /* Variables for eigen estimation - START*/
  PetscInt       stored_max_it = ksp->max_it;
  PetscScalar    alphaold = 0,betaold = 1.0,*e = 0,*d = 0;/* Variables for eigen estimation  - FINISH */

  PetscFunctionBegin;

#define VecXDot(x,y,a) (((fcg->type) == (KSP_CG_HERMITIAN)) ? VecDot(x,y,a) : VecTDot(x,y,a))
#define VecXMDot(a,b,c,d) (((fcg->type) == (KSP_CG_HERMITIAN)) ? VecMDot(a,b,c,d) : VecMTDot(a,b,c,d))

  X             = ksp->vec_sol;
  B             = ksp->vec_rhs;
  R             = ksp->work[0];
  Z             = ksp->work[1];

  ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr);
  if (eigs) {e = fcg->e; d = fcg->d; e[0] = 0.0; }
  /* Compute initial residual needed for convergence check*/
  ksp->its = 0;
  if (!ksp->guess_zero) {
    ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr);
    ierr = VecAYPX(R,-1.0,B);CHKERRQ(ierr);                    /*   r <- b - Ax     */
  } 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 <- dqrt(z'*z) = sqrt(e'*A'*B'*B*A*e)     */
      break;
    case KSP_NORM_UNPRECONDITIONED:
      ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr);              /*   dp <- sqrt(r'*r) = sqrt(e'*A'*A*e)     */
      break;
    case KSP_NORM_NATURAL:
      ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr);               /*   z <- Br         */
      ierr = VecXDot(R,Z,&s);CHKERRQ(ierr);
      dp = PetscSqrtReal(PetscAbsScalar(s));                   /*   dp <- sqrt(r'*z) = sqrt(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]);
  }

  /* Initial Convergence Check */
  ierr       = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
  ierr       = KSPMonitor(ksp,0,dp);CHKERRQ(ierr);
  ksp->rnorm = dp;
  if (ksp->normtype == KSP_NORM_NONE) {
    ierr = KSPConvergedSkip(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  } else {
    ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  }
  if (ksp->reason) PetscFunctionReturn(0);

  /* Apply PC if not already done for convergence check */
  if(ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->normtype == KSP_NORM_NONE){
    ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr);               /*   z <- Br         */
  }

  i = 0;
  do {
    ksp->its = i+1;

    /*  If needbe, allocate a new chunk of vectors in P and C */
    ierr = KSPAllocateVectors_FCG(ksp,i+1,fcg->vecb);CHKERRQ(ierr);

    /* Note that we wrap around and start clobbering old vectors */
    idx = i % (fcg->mmax+1);
    Pcurr = fcg->Pvecs[idx];
    Ccurr = fcg->Cvecs[idx];

    /* Compute a new column of P (Currently does not support modified G-S or iterative refinement)*/
    switch(fcg->truncstrat){
      case KSP_FCG_TRUNC_TYPE_NOTAY :
        mi = PetscMax(1,i%(fcg->mmax+1));
        break;
      case KSP_FCG_TRUNC_TYPE_STANDARD :
        mi = fcg->mmax;
        break;
      default:
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unrecognized FCG Truncation Strategy");CHKERRQ(ierr);
    }
    ierr = VecCopy(Z,Pcurr);CHKERRQ(ierr);

    {
      PetscInt l,ndots;

      l = PetscMax(0,i-mi);
      ndots = i-l;
      if (ndots){
        PetscInt    j;
//.........这里部分代码省略.........
开发者ID:pombredanne,项目名称:petsc,代码行数:101,代码来源:fcg.c


示例8: MatGetOrdering_AWBM

PETSC_EXTERN PetscErrorCode MatGetOrdering_AWBM(Mat A, MatOrderingType type, IS *permR, IS *permC)
{
  Vec *scalR, *scalC, scalRVec, scalCVec;
  scalR = &scalRVec; scalC = &scalCVec;

  /* EVERYTHING IS WRITTEN AS IF THE MATRIX WERE COLUMN-MAJOR */
  Mat_SeqAIJ      *aij = (Mat_SeqAIJ *) A->data;
  PetscInt         n   = A->rmap->n; /* Number of local columns */
  PetscInt         m   = A->cmap->n; /* Number of local rows */
  PetscInt        *match;            /* The row matched to each column, and inverse column permutation */
  PetscInt        *matchR;           /* The column matched to each row */
  PetscInt        *p;                /* The column permutation */
  const PetscInt  *ia  = aij->i;
  const PetscInt  *ja  = aij->j;
  const MatScalar *a   = aij->a;
  Vec              colMax;
  PetscScalar     *a_j, *sr, *sc;
  PetscReal       *weights /* c_ij */, *u /* u_i */, *v /* v_j */, eps = PETSC_SQRT_MACHINE_EPSILON;
  PetscInt         debug = 0, r, c, r1, c1;
  PetscErrorCode   ierr;

  PetscFunctionBegin;
  ierr = PetscOptionsGetInt(NULL, "-debug", &debug, NULL);CHKERRQ(ierr);
  ierr = MatGetVecs(A, NULL, &colMax);CHKERRQ(ierr);
  ierr = MatGetRowMaxAbs(A, colMax, NULL);CHKERRQ(ierr);
  ierr = PetscMalloc2(n, &match, m, &matchR);CHKERRQ(ierr);
  ierr = PetscMalloc1(n, &p);CHKERRQ(ierr);
  ierr = PetscCalloc3(m, &u, n, &v, ia[n], &weights);CHKERRQ(ierr);
  for (c = 0; c < n; ++c) match[c] = -1;
  /* Compute weights */
  ierr = VecGetArray(colMax, &a_j);CHKERRQ(ierr);
  for (c = 0; c < n; ++c) {
    for (r = ia[c]; r < ia[c+1]; ++r) {
      PetscReal ar = PetscAbsScalar(a[r]);

      if (ar == 0.0) weights[r] = PETSC_MAX_REAL;
      else           weights[r] = log(a_j[c]/ar);
    }
  }
  /* Compute local row weights */
  for (r = 0; r < m; ++r) u[r] = PETSC_MAX_REAL;
  for (c = 0; c < n; ++c) {
    for (r = ia[c]; r < ia[c+1]; ++r) {
      u[ja[r]] = PetscMin(u[ja[r]], weights[r]);
    }
  }
  /* Compute local column weights */
  for (c = 0; c < n; ++c) {
    v[c] = PETSC_MAX_REAL;
    for (r = ia[c]; r < ia[c+1]; ++r) {
      v[c] = PetscMin(v[c], weights[r] - u[ja[r]]);
    }
  }
  for (r = 0; r < m; ++r) matchR[r] = -1;
  /* Match columns */
  ierr = CheckUnmatched(n, match, matchR);CHKERRQ(ierr);
  for (c = 0; c < n; ++c) {
    /* if (match[c] >= 0) continue; */
    if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "Row %d\n  Weights:", c);CHKERRQ(ierr);}
    for (r = ia[c]; r < ia[c+1]; ++r) {
      PetscReal weight = weights[r] - u[ja[r]] - v[c];
      if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " %g", weight);CHKERRQ(ierr);}
      if ((weight <= eps) && (matchR[ja[r]] < 0)) {
        if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "Matched %d -- %d\n", c, ja[r]);CHKERRQ(ierr);}
        match[c]      = ja[r];
        matchR[ja[r]] = c;
        break;
      }
    }
    if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);}
  }
  /* Deal with unmatched columns */
  ierr = CheckUnmatched(n, match, matchR);CHKERRQ(ierr);
  for (c = 0; c < n; ++c) {
    if (match[c] >= 0) continue;
    for (r = ia[c]; r < ia[c+1]; ++r) {
      PetscReal weight = weights[r] - u[ja[r]] - v[c];
      if (weight > eps) continue;
      /* \bar c_ij = 0 and (r, j1) \in M */
      c1 = matchR[ja[r]];
      for (r1 = ia[c1]; r1 < ia[c1+1]; ++r1) {
        PetscReal weight1 = weights[r1] - u[ja[r1]] - v[c1];
        if ((matchR[ja[r1]] < 0) && (weight1 <= eps)) {
          /* (r, c1) in M is replaced by (r, c) and (r1, c1) */
          if (debug) {
            ierr = PetscPrintf(PETSC_COMM_SELF, "Replaced match %d -- %d\n", c1, ja[r]);CHKERRQ(ierr);
            ierr = PetscPrintf(PETSC_COMM_SELF, "  Added  match %d -- %d\n", c,  ja[r]);CHKERRQ(ierr);
            ierr = PetscPrintf(PETSC_COMM_SELF, "  Added  match %d -- %d\n", c1, ja[r1]);CHKERRQ(ierr);
          }
          match[c]       = ja[r];
          matchR[ja[r]]  = c;
          match[c1]      = ja[r1];
          matchR[ja[r1]] = c1;
          break;
        }
      }
      if (match[c] >= 0) break;
    }
  }
  /* Allow matching with non-optimal rows */
//.........这里部分代码省略.........
开发者ID:spikegpu,项目名称:spike-petsc,代码行数:101,代码来源:petsc_mat_awbm.c


示例9: KSPPGMRESUpdateHessenberg

/*
.  it - column of the Hessenberg that is complete, PGMRES is actually computing two columns ahead of this
 */
static PetscErrorCode KSPPGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool *hapend,PetscReal *res)
{
  PetscScalar    *hh,*cc,*ss,*rs;
  PetscInt       j;
  PetscReal      hapbnd;
  KSP_PGMRES     *pgmres = (KSP_PGMRES*)(ksp->data);
  PetscErrorCode ierr;

  PetscFunctionBegin;
  hh = HH(0,it);   /* pointer to beginning of column to update */
  cc = CC(0);      /* beginning of cosine rotations */
  ss = SS(0);      /* beginning of sine rotations */
  rs = RS(0);      /* right hand side of least squares system */

  /* The Hessenberg matrix is now correct through column it, save that form for possible spectral analysis */
  for (j=0; j<=it+1; j++) *HES(j,it) = hh[j];

  /* check for the happy breakdown */
  hapbnd = PetscMin(PetscAbsScalar(hh[it+1] / rs[it]),pgmres->haptol);
  if (PetscAbsScalar(hh[it+1]) < hapbnd) {
    ierr    = PetscInfo4(ksp,"Detected happy breakdown, current hapbnd = %14.12e H(%D,%D) = %14.12e\n",(double)hapbnd,it+1,it,(double)PetscAbsScalar(*HH(it+1,it)));CHKERRQ(ierr);
    *hapend = PETSC_TRUE;
  }

  /* Apply all the previously computed plane rotations to the new column
     of the Hessenberg matrix */
  /* Note: this uses the rotation [conj(c)  s ; -s   c], c= cos(theta), s= sin(theta),
     and some refs have [c   s ; -conj(s)  c] (don't be confused!) */

  for (j=0; j<it; j++) {
    PetscScalar hhj = hh[j];
    hh[j]   = PetscConj(cc[j])*hhj + ss[j]*hh[j+1];
    hh[j+1] =          -ss[j] *hhj + cc[j]*hh[j+1];
  }

  /*
    compute the new plane rotation, and apply it to:
     1) the right-hand-side of the Hessenberg system (RS)
        note: it affects RS(it) and RS(it+1)
     2) the new column of the Hessenberg matrix
        note: it affects HH(it,it) which is currently pointed to
        by hh and HH(it+1, it) (*(hh+1))
    thus obtaining the updated value of the residual...
  */

  /* compute new plane rotation */

  if (!*hapend) {
    PetscReal delta = PetscSqrtReal(PetscSqr(PetscAbsScalar(hh[it])) + PetscSqr(PetscAbsScalar(hh[it+1])));
    if (delta == 0.0) {
      ksp->reason = KSP_DIVERGED_NULL;
      PetscFunctionReturn(0);
    }

    cc[it] = hh[it] / delta;    /* new cosine value */
    ss[it] = hh[it+1] / delta;  /* new sine value */

    hh[it]   = PetscConj(cc[it])*hh[it] + ss[it]*hh[it+1];
    rs[it+1] = -ss[it]*rs[it];
    rs[it]   = PetscConj(cc[it])*rs[it];
    *res     = PetscAbsScalar(rs[it+1]);
  } else { /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply
            another rotation matrix (so RH doesn't change).  The new residual is
            always the new sine term times the residual from last time (RS(it)),
            but now the new sine rotation would be zero...so the residual should
            be zero...so we will multiply "zero" by the last residual.  This might
            not be exactly what we want to do here -could just return "zero". */

    *res = 0.0;
  }
  PetscFunctionReturn(0);
}
开发者ID:firedrakeproject,项目名称:petsc,代码行数:75,代码来源:pgmres.c


示例10: KSPSolve_TCQMR

static PetscErrorCode KSPSolve_TCQMR(KSP ksp)
{
  PetscReal      rnorm0,rnorm,dp1,Gamma;
  PetscScalar    theta,ep,cl1,sl1,cl,sl,sprod,tau_n1,f;
  PetscScalar    deltmp,rho,beta,eptmp,ta,s,c,tau_n,delta;
  PetscScalar    dp11,dp2,rhom1,alpha,tmp;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ksp->its = 0;

  ierr = KSPInitialResidual(ksp,x,u,v,r,b);CHKERRQ(ierr);
  ierr = VecNorm(r,NORM_2,&rnorm0);CHKERRQ(ierr);          /*  rnorm0 = ||r|| */

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

  ierr  = VecSet(um1,0.0);CHKERRQ(ierr);
  ierr  = VecCopy(r,u);CHKERRQ(ierr);
  rnorm = rnorm0;
  tmp   = 1.0/rnorm; ierr = VecScale(u,tmp);CHKERRQ(ierr);
  ierr  = VecSet(vm1,0.0);CHKERRQ(ierr);
  ierr  = VecCopy(u,v);CHKERRQ(ierr);
  ierr  = VecCopy(u,v0);CHKERRQ(ierr);
  ierr  = VecSet(pvec1,0.0);CHKERRQ(ierr);
  ierr  = VecSet(pvec2,0.0);CHKERRQ(ierr);
  ierr  = VecSet(p,0.0);CHKERRQ(ierr);
  theta = 0.0;
  ep    = 0.0;
  cl1   = 0.0;
  sl1   = 0.0;
  cl    = 0.0;
  sl    = 0.0;
  sprod = 1.0;
  tau_n1= rnorm0;
  f     = 1.0;
  Gamma = 1.0;
  rhom1 = 1.0;

  /*
   CALCULATE SQUARED LANCZOS  vectors
   */
  ierr = (*ksp->converged)(ksp,ksp->its,rnorm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  while (!ksp->reason) {
    ierr = KSPMonitor(ksp,ksp->its,rnorm);CHKERRQ(ierr);
    ksp->its++;

    ierr   = KSP_PCApplyBAorAB(ksp,u,y,vtmp);CHKERRQ(ierr); /* y = A*u */
    ierr   = VecDot(y,v0,&dp11);CHKERRQ(ierr);
    ierr   = VecDot(u,v0,&dp2);CHKERRQ(ierr);
    alpha  = dp11 / dp2;                          /* alpha = v0'*y/v0'*u */
    deltmp = alpha;
    ierr   = VecCopy(y,z);CHKERRQ(ierr);
    ierr   = VecAXPY(z,-alpha,u);CHKERRQ(ierr); /* z = y - alpha u */
    ierr   = VecDot(u,v0,&rho);CHKERRQ(ierr);
    beta   = rho / (f*rhom1);
    rhom1  = rho;
    ierr   = VecCopy(z,utmp);CHKERRQ(ierr);    /* up1 = (A-alpha*I)*
                                                (z-2*beta*p) + f*beta*
                                                beta*um1 */
    ierr   = VecAXPY(utmp,-2.0*beta,p);CHKERRQ(ierr);
    ierr   = KSP_PCApplyBAorAB(ksp,utmp,up1,vtmp);CHKERRQ(ierr);
    ierr   = VecAXPY(up1,-alpha,utmp);CHKERRQ(ierr);
    ierr   = VecAXPY(up1,f*beta*beta,um1);CHKERRQ(ierr);
    ierr   = VecNorm(up1,NORM_2,&dp1);CHKERRQ(ierr);
    f      = 1.0 / dp1;
    ierr   = VecScale(up1,f);CHKERRQ(ierr);
    ierr   = VecAYPX(p,-beta,z);CHKERRQ(ierr);   /* p = f*(z-beta*p) */
    ierr   = VecScale(p,f);CHKERRQ(ierr);
    ierr   = VecCopy(u,um1);CHKERRQ(ierr);
    ierr   = VecCopy(up1,u);CHKERRQ(ierr);
    beta   = beta/Gamma;
    eptmp  = beta;
    ierr   = KSP_PCApplyBAorAB(ksp,v,vp1,vtmp);CHKERRQ(ierr);
    ierr   = VecAXPY(vp1,-alpha,v);CHKERRQ(ierr);
    ierr   = VecAXPY(vp1,-beta,vm1);CHKERRQ(ierr);
    ierr   = VecNorm(vp1,NORM_2,&Gamma);CHKERRQ(ierr);
    ierr   = VecScale(vp1,1.0/Gamma);CHKERRQ(ierr);
    ierr   = VecCopy(v,vm1);CHKERRQ(ierr);
    ierr   = VecCopy(vp1,v);CHKERRQ(ierr);

    /*
       SOLVE  Ax = b
     */
    /* Apply last two Given's (Gl-1 and Gl) rotations to (beta,alpha,Gamma) */
    if (ksp->its > 2) {
      theta =  sl1*beta;
      eptmp = -cl1*beta;
    }
    if (ksp->its > 1) {
      ep     = -cl*eptmp + sl*alpha;
      deltmp = -sl*eptmp - cl*alpha;
    }
    if (PetscAbsReal(Gamma) > PetscAbsScalar(deltmp)) {
      ta = -deltmp / Gamma;
      s  = 1.0 / PetscSqrtScalar(1.0 + ta*ta);
      c  = s*ta;
    } else {
      ta = -Gamma/deltmp;
      c  = 1.0 / PetscSqrtScalar(1.0 + ta*ta);
//.........这里部分代码省略.........
开发者ID:OpenCMISS-Dependencies,项目名称:petsc,代码行数:101,代码来源:tcqmr.c


示例11: KSPLGMRESCycle


//.........这里部分代码省略.........

        /*LGMRES_MOD: decide whether this is an arnoldi step or an aug step */
        if (loc_it < it_arnoldi) { /* Arnoldi */
            ierr = KSP_PCApplyBAorAB(ksp,VEC_VV(loc_it),VEC_VV(1+loc_it),VEC_TEMP_MATOP);
            CHKERRQ(ierr);
        } else { /*aug step */
            order = loc_it - it_arnoldi + 1; /* which aug step */
            for (ii=0; ii<aug_dim; ii++) {
                if (lgmres->aug_order[ii] == order) {
                    spot = ii;
                    break; /* must have this because there will be duplicates before aug_ct = aug_dim */
                }
            }

            ierr = VecCopy(A_AUGVEC(spot), VEC_VV(1+loc_it));
            CHKERRQ(ierr);
            /*note: an alternate implementation choice would be to only save the AUGVECS and
              not A_AUGVEC and then apply the PC here to the augvec */
        }

        /* update hessenberg matrix and do Gram-Schmidt - new direction is in
           VEC_VV(1+loc_it)*/
        ierr = (*lgmres->orthog)(ksp,loc_it);
        CHKERRQ(ierr);

        /* new entry in hessenburg is the 2-norm of our new direction */
        ierr = VecNorm(VEC_VV(loc_it+1),NORM_2,&tt);
        CHKERRQ(ierr);
        *HH(loc_it+1,loc_it)   = tt;
        *HES(loc_it+1,loc_it)  = tt;


        /* check for the happy breakdown */
        hapbnd  = PetscAbsScalar(tt / *GRS(loc_it));/* GRS(loc_it) contains the res_norm from the last iteration  */
        if (hapbnd > lgmres->haptol) hapbnd = lgmres->haptol;
        if (tt > hapbnd) {
            tmp = 1.0/tt;
            ierr = VecScale(VEC_VV(loc_it+1),tmp);
            CHKERRQ(ierr); /* scale new direction by its norm */
        } else {
            ierr = PetscInfo2(ksp,"Detected happy breakdown, current hapbnd = %G tt = %G\n",hapbnd,tt);
            CHKERRQ(ierr);
            hapend = PETSC_TRUE;
        }

        /* Now apply rotations to new col of hessenberg (and right side of system),
           calculate new rotation, and get new residual norm at the same time*/
        ierr = KSPLGMRESUpdateHessenberg(ksp,loc_it,hapend,&res);
        CHKERRQ(ierr);
        if (ksp->reason) break;

        loc_it++;
        lgmres->it  = (loc_it-1);  /* Add this here in case it has converged */

        ierr = PetscObjectTakeAccess(ksp);
        CHKERRQ(ierr);
        ksp->its++;
        ksp->rnorm = res;
        ierr = PetscObjectGrantAccess(ksp);
        CHKERRQ(ierr);

        ierr = (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);
        CHKERRQ(ierr);

        /* Catch error in happy breakdown and signal convergence and break from loop */
        if (hapend) {
开发者ID:erdc-cm,项目名称:petsc-dev,代码行数:67,代码来源:lgmres.c


示例12: PCISApplyInvSchur

PetscErrorCode  PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
{
  PetscErrorCode ierr;
  PC_IS          *pcis = (PC_IS*)(pc->data);

  PetscFunctionBegin;
  /*
    Neumann solvers.
    Applying the inverse of the local Schur complement, i.e, solving a Neumann
    Problem with zero at the interior nodes of the RHS and extracting the interface
    part of the solution. inverse Schur complement is applied to b and the result
    is stored in x.
  */
  /* Setting the RHS vec1_N */
  ierr = VecSet(vec1_N,0.0);CHKERRQ(ierr);
  ierr = VecScatterBegin(pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
  ierr = VecScatterEnd  (pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
  /* Checking for consistency of the RHS */
  {
    PetscBool flg = PETSC_FALSE;
    ierr = PetscOptionsGetBool(NULL,"-pc_is_check_consistency",&flg,NULL);CHKERRQ(ierr);
    if (flg) {
      PetscScalar average;
      PetscViewer viewer;
      ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);CHKERRQ(ierr);

      ierr    = VecSum(vec1_N,&average);CHKERRQ(ierr);
      average = average / ((PetscReal)pcis->n);
      ierr    = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr);
      if (pcis->pure_neumann) {
        ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is floating. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));CHKERRQ(ierr);
      } else {
        ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is fixed.    Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));CHKERRQ(ierr);
      }
      ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
      ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr);
    }
  }
  /* Solving the system for vec2_N */
  ierr = KSPSolve(pcis->ksp_N,vec1_N,vec2_N);CHKERRQ(ierr);
  /* Extracting the local interface vector out of the solution */
  ierr = VecScatterBegin(pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterEnd  (pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}
开发者ID:hsahasra,项目名称:petsc-magma-dense-mat,代码行数:45,代码来源:pcis.c


示例13: KSPLGMRESBuildSoln

static PetscErrorCode KSPLGMRESBuildSoln(PetscScalar* nrs,Vec vguess,Vec vdest,KSP ksp,PetscInt it)
{
    PetscScalar    tt;
    PetscErrorCode ierr;
    PetscInt       ii,k,j;
    KSP_LGMRES     *lgmres = (KSP_LGMRES *)(ksp->data);
    /*LGMRES_MOD */
    PetscInt       it_arnoldi, it_aug;
    PetscInt       jj, spot = 0;

    PetscFunctionBegin;
    /* Solve for solution vector that minimizes the residual */

    /* If it is < 0, no lgmres steps have been performed */
    if (it < 0) {
        ierr = VecCopy(vguess,vdest);
        CHKERRQ(ierr); /* VecCopy() is smart, exists immediately if vguess == vdest */
        PetscFunctionReturn(0);
    }

    /* so (it+1) lgmres steps HAVE been performed */

    /* LGMRES_MOD - determine if we need to use augvecs for the soln  - do not assume that
       this is called after the total its allowed for an approx space */
    if (lgmres->approx_constant) {
        it_arnoldi = lgmres->max_k - lgmres->aug_ct;
    } else {
        it_arnoldi = lgmres->max_k - lgmres->aug_dim;
    }
    if (it_arnoldi >= it +1) {
        it_aug = 0;
        it_arnoldi = it+1;
    } else {
        it_aug = (it + 1) - it_arnoldi;
    }

    /* now it_arnoldi indicates the number of matvecs that took place */
    lgmres->matvecs += it_arnoldi;


    /* solve the upper triangular system - GRS is the right side and HH is
       the upper triangular matrix  - put soln in nrs */
    if (*HH(it,it) == 0.0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"HH(it,it) is identically zero; it = %D GRS(it) = %G",it,PetscAbsScalar(*GRS(it)));
    if (*HH(it,it) != 0.0) {
        nrs[it] = *GRS(it) / *HH(it,it);
    } else {
        nrs[it] = 0.0;
    }

    for (ii=1; ii<=it; ii++) {
        k   = it - ii;
        tt  = *GRS(k);
        for (j=k+1; j<=it; j++) tt  = tt - *HH(k,j) * nrs[j];
        nrs[k]   = tt / *HH(k,k);
    }

    /* Accumulate the correction to the soln of the preconditioned prob. in VEC_TEMP */
    ierr = VecSet(VEC_TEMP,0.0);
    CHKERRQ(ierr); /* set VEC_TEMP components to 0 */

    /*LGMRES_MOD - if augmenting has happened we need to form the solution
      using the augvecs */
    if (!it_aug) { /* all its are from arnoldi */
        ierr = VecMAXPY(VEC_TEMP,it+1,nrs,&VEC_VV(0));
        CHKERRQ(ierr);
    } else { /*use aug vecs */
        /*first do regular krylov directions */
        ierr = VecMAXPY(VEC_TEMP,it_arnoldi,nrs,&VEC_VV(0));
        CHKERRQ(ierr);
        /*now add augmented portions - add contribution of aug vectors one at a time*/


        for (ii=0; ii<it_aug; ii++) {
            for (jj=0; jj<lgmres->aug_dim; jj++) {
                if (lgmres->aug_order[jj] == (ii+1)) {
                    spot = jj;
                    break; /* must have this because there will be duplicates before aug_ct = aug_dim */
                }
            }
            ierr = VecAXPY(VEC_TEMP,nrs[it_arnoldi+ii],AUGVEC(spot));
            CHKERRQ(ierr);
        }
    }
    /* now VEC_TEMP is what we want to keep for augmenting purposes - grab before the
       preconditioner is "unwound" from right-precondtioning*/
    ierr = VecCopy(VEC_TEMP, AUG_TEMP);
    CHKERRQ(ierr);

    ierr = KSPUnwindPreconditioner(ksp,VEC_TEMP,VEC_TEMP_MATOP);
    CHKERRQ(ierr);

    /* add solution to previous solution */
    /* put updated solution into vdest.*/
    ierr = VecCopy(vguess,vdest);
    CHKERRQ(ierr);
    ierr = VecAXPY(vdest,1.0,VEC_TEMP);
    CHKERRQ(ierr);

    PetscFunctionReturn(0);
}
开发者ID:erdc-cm,项目名称:petsc-dev,代码行数:100,代码来源:lgmres.c


示例14: KSPSolve_MINRES

该文章已有0人参与评论

请发表评论

全部评论

专题导读
上一篇:
C++ PetscCalloc1函数代码示例发布时间:2022-05-30
下一篇:
C++ PetSaveMode函数代码示例发布时间:2022-05-30
热门推荐
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap