C++ PetscSqrtScalar函数代码示例

OStack程序员社区-中国程序员成长平台 › 门户 › 编程› C++›C++教程

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

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

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

示例1: BSSCR_PCScGtKGUseStandardScaling

PetscErrorCode BSSCR_PCScGtKGUseStandardScaling( PC pc ) 
{
	PC_SC_GtKG ctx = (PC_SC_GtKG)pc->data;
	Mat K,G,D,C;
	Vec rG;
	PetscScalar rg2, rg, ra;  
	PetscInt N;
	Vec rA, rC;
	Vec L1,L2, R1,R2;
	
	BSSCR_BSSCR_pc_error_ScGtKG( pc, __func__ );
	
	
	L1 = ctx->X1;
	L2 = ctx->X2;
	
	R1 = ctx->Y1;
	R2 = ctx->Y2;
	
	
	rA = L1;
	rC = L2;
	
	K = ctx->F;
	G = ctx->Bt;
	D = ctx->B;
	C = ctx->C;
	
	VecDuplicate( rA, &rG );
	
	/* Get magnitude of K */  
	MatGetRowMax( K, rA, PETSC_NULL );
	
	VecSqrt( rA );  
	VecReciprocal( rA );
	
	VecDot( rA,rA, &ra );
	VecGetSize( rA, &N );
	ra = PetscSqrtScalar( ra/N );
	
	
	/* Get magnitude of G */
	MatGetRowMax( G, rG, PETSC_NULL );
	
	VecDot( rG, rG, &rg2 );
	VecGetSize( rG, &N );
	rg = PetscSqrtScalar(rg2/N);
	//      printf("rg = %f \n", rg );
	
	VecSet( rC, 1.0/(rg*ra) );
	
	Stg_VecDestroy(&rG );
	
	VecCopy( L1, R1 );
	VecCopy( L2, R2 );
	
	PetscFunctionReturn(0);
}

开发者ID:AngelValverdePerez，项目名称:underworld2，代码行数:58，代码来源:pc_ScaledGtKG.c

示例2: sol_true

static PetscErrorCode sol_true(PetscReal t,Vec U)
{
  PetscErrorCode ierr;
  PetscScalar    *u;

  PetscFunctionBegin;
  ierr = VecGetArray(U,&u);CHKERRQ(ierr);
  u[0] = PetscSqrtScalar(1.0+PetscCosScalar(t));
  u[1] = PetscSqrtScalar(2.0+PetscCosScalar(5.0*t));
  ierr = VecRestoreArray(U,&u);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

开发者ID:petsc，项目名称:petsc，代码行数:12，代码来源:ex3.c

示例3: RHSFunction_Hull1972B4

PetscErrorCode RHSFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, 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[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
  f[1] =  y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
  f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
  ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr);
  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

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

示例4: SphericalDistribution

// http://mathworld.wolfram.com/SpherePointPicking.html
PetscErrorCode SphericalDistribution(PetscRandom rnd, PetscReal a, Coor *d)
{
  PetscReal u,q;

  PetscRandomGetValue(rnd,&u); // [0,1]
  u = (1-a) * u + a; // [a,1]
  PetscRandomGetValue(rnd,&q); // [0,1]
  q = 2 * PETSC_PI * q; // [0, 2pi]

  d->x = PetscSqrtScalar( 1 - u*u ) * cos(q);
  d->y = PetscSqrtScalar( 1 - u*u ) * sin(q);
  d->z = u;

  return 0;
}

开发者ID:adrielb，项目名称:DCell，代码行数:16，代码来源:Fiber.c

示例5: KSPSolve_SpecEst

static PetscErrorCode  KSPSolve_SpecEst(KSP ksp)
{
  PetscErrorCode ierr;
  KSP_SpecEst    *spec = (KSP_SpecEst*)ksp->data;

  PetscFunctionBegin;
  if (spec->current) {
    ierr = KSPSolve(spec->kspcheap,ksp->vec_rhs,ksp->vec_sol);CHKERRQ(ierr);
    ierr = KSPSpecEstPropagateUp(ksp,spec->kspcheap);CHKERRQ(ierr);
  } else {
    PetscInt  i,its,neig;
    PetscReal *real,*imag,rad = 0;
    ierr = KSPSolve(spec->kspest,ksp->vec_rhs,ksp->vec_sol);CHKERRQ(ierr);
    ierr = KSPSpecEstPropagateUp(ksp,spec->kspest);CHKERRQ(ierr);
    ierr = KSPComputeExtremeSingularValues(spec->kspest,&spec->max,&spec->min);CHKERRQ(ierr);

    ierr = KSPGetIterationNumber(spec->kspest,&its);CHKERRQ(ierr);
    ierr = PetscMalloc2(its,PetscReal,&real,its,PetscReal,&imag);CHKERRQ(ierr);
    ierr = KSPComputeEigenvalues(spec->kspest,its,real,imag,&neig);CHKERRQ(ierr);
    for (i=0; i<neig; i++) {
      /* We would really like to compute w (nominally 1/radius) to minimize |1-wB|.  Empirically it
         is better to compute rad = |1-B| than rad = |B|.  There must be a cheap way to do better. */
      rad = PetscMax(rad,PetscRealPart(PetscSqrtScalar((PetscScalar)(PetscSqr(real[i]-1.) + PetscSqr(imag[i])))));
    }
    ierr = PetscFree2(real,imag);CHKERRQ(ierr);
    spec->radius = rad;

    ierr = KSPChebyshevSetEigenvalues(spec->kspcheap,spec->max*spec->maxfactor,spec->min*spec->minfactor);CHKERRQ(ierr);
    ierr = KSPRichardsonSetScale(spec->kspcheap,spec->richfactor/spec->radius);
    ierr = PetscInfo3(ksp,"Estimated singular value min=%G max=%G, spectral radius=%G",spec->min,spec->max,spec->radius);CHKERRQ(ierr);
    spec->current = PETSC_TRUE;
  }
  PetscFunctionReturn(0);
}

开发者ID:Kun-Qu，项目名称:petsc，代码行数:34，代码来源:specest.c

示例6: VecNorm_Nest

static PetscErrorCode VecNorm_Nest(Vec xin,NormType type,PetscReal *z)
{
  Vec_Nest       *bx = (Vec_Nest*)xin->data;
  PetscInt       i,nr;
  PetscReal      z_i;
  PetscReal      _z;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  nr = bx->nb;
  _z = 0.0;

  if (type == NORM_2) {
    PetscScalar dot;
    ierr = VecDot(xin,xin,&dot);CHKERRQ(ierr);
    _z = PetscAbsScalar(PetscSqrtScalar(dot));
  } else if (type == NORM_1) {
    for (i=0; i<nr; i++) {
      ierr = VecNorm(bx->v[i],type,&z_i);CHKERRQ(ierr);
      _z = _z + z_i;
    }
  } else if (type == NORM_INFINITY) {
    for (i=0; i<nr; i++) {
      ierr = VecNorm(bx->v[i],type,&z_i);CHKERRQ(ierr);
      if (z_i > _z) _z = z_i;
    }
  }

  *z = _z;
  PetscFunctionReturn(0);
}

开发者ID:plguhur，项目名称:petsc，代码行数:31，代码来源:vecnest.c

示例7: StokesCalcError

PetscErrorCode StokesCalcError(Stokes *s)
{
  PetscScalar    scale = PetscSqrtScalar(s->nx*s->ny);
  PetscReal      val;
  Vec            y0, y1;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  /* error y-x */
  ierr = VecAXPY(s->y, -1.0, s->x);CHKERRQ(ierr);
  /* ierr = VecView(s->y, (PetscViewer)PETSC_VIEWER_DEFAULT);CHKERRQ(ierr); */

  /* error in velocity */
  ierr = VecGetSubVector(s->y, s->isg[0], &y0);CHKERRQ(ierr);
  ierr = VecNorm(y0, NORM_2, &val);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD," discretization error u = %g\n",(double)(PetscRealPart(val/scale)));CHKERRQ(ierr);
  ierr = VecRestoreSubVector(s->y, s->isg[0], &y0);CHKERRQ(ierr);

  /* error in pressure */
  ierr = VecGetSubVector(s->y, s->isg[1], &y1);CHKERRQ(ierr);
  ierr = VecNorm(y1, NORM_2, &val);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD," discretization error p = %g\n",(double)(PetscRealPart(val/scale)));CHKERRQ(ierr);
  ierr = VecRestoreSubVector(s->y, s->isg[1], &y1);CHKERRQ(ierr);

  /* total error */
  ierr = VecNorm(s->y, NORM_2, &val);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD," discretization error [u,p] = %g\n", (double)PetscRealPart((val/scale)));CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

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

示例8: InitialCondition

void InitialCondition(TS ts, Vec X)
{
    PetscScalar *x;

    VecGetArray(X, &x);

    for (PetscInt j=0; j<N2; j++) {
        for (PetscInt i=0; i<N1; i++) {
            for (PetscInt var=0; var<DOF; var++) {
            
                PetscScalar X1Coord = X1_MIN + DX1/2. + i*DX1;
                PetscScalar X2Coord = X2_MIN + DX2/2. + j*DX2;

                PetscScalar X1Center = (X1_MIN + X1_MAX)/2.;
                PetscScalar X2Center = (X2_MIN + X2_MAX)/2.;

                PetscScalar r = PetscSqrtScalar(
                                    PetscPowScalar(X1Coord-X1Center, 2.0) + 
                                    PetscPowScalar(X2Coord-X2Center, 2.0));

                x[INDEX_GLOBAL(i,j,var)] = exp(-r*r/.01);
            }
        }
    }

    VecRestoreArray(X, &x);
}

开发者ID:mchandra，项目名称:PetscOpenCL，代码行数:27，代码来源:petsc_opencl.cpp

示例9: IFunction_Hull1972B4

PetscErrorCode IFunction_Hull1972B4(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[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
  f[1] =  y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
  f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[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，代码行数:17，代码来源:ex31.c

示例10: ini_bou

PetscErrorCode ini_bou(Vec X,AppCtx* user)
{
  PetscErrorCode ierr;
  DM             cda;
  DMDACoor2d     **coors;
  PetscScalar    **p;
  Vec            gc;
  PetscInt       i,j;
  PetscInt       xs,ys,xm,ym,M,N;
  PetscScalar    xi,yi;
  PetscScalar    sigmax=user->sigmax,sigmay=user->sigmay;
  PetscScalar    rho   =user->rho;
  PetscScalar    mux   =user->mux,muy=user->muy;
  PetscMPIInt    rank;
  PetscScalar    sum;

  PetscFunctionBeginUser;
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
  ierr = DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
  user->dx = (user->xmax - user->xmin)/(M-1); user->dy = (user->ymax - user->ymin)/(N-1);
  ierr = DMGetCoordinateDM(user->da,&cda);CHKERRQ(ierr);
  ierr = DMGetCoordinates(user->da,&gc);CHKERRQ(ierr);
  ierr = DMDAVecGetArray(cda,gc,&coors);CHKERRQ(ierr);
  ierr = DMDAVecGetArray(user->da,X,&p);CHKERRQ(ierr);
  ierr = DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);

  /* mux and muy need to be grid points in the x and y-direction otherwise the solution goes unstable
     muy is set by choosing the y domain, no. of grid points along y-direction so that muy is a grid point
     in the y-direction. We only modify mux here
  */
  mux = user->mux = coors[0][M/2+10].x; /* For -pi < x < pi, this should be some angle between 0 and pi/2 */
  if (user->nonoiseinitial) {
    for (i=xs; i < xs+xm; i++) {
      for (j=ys; j < ys+ym; j++) {
        xi = coors[j][i].x; yi = coors[j][i].y;
        if ((xi == mux) && (yi == muy)) {
          p[j][i] = 1.0;
        }
      }
    }
  } else {
    /* Change PM_min accordingly */
    user->PM_min = user->Pmax*sin(mux);
    for (i=xs; i < xs+xm; i++) {
      for (j=ys; j < ys+ym; j++) {
        xi = coors[j][i].x; yi = coors[j][i].y;
        p[j][i] = (0.5/(PETSC_PI*sigmax*sigmay*PetscSqrtScalar(1.0-rho*rho)))*PetscExpScalar(-0.5/(1-rho*rho)*(PetscPowScalar((xi-mux)/sigmax,2) + PetscPowScalar((yi-muy)/sigmay,2) - 2*rho*(xi-mux)*(yi-muy)/(sigmax*sigmay)));
      }
    }
  }
  ierr = DMDAVecRestoreArray(cda,gc,&coors);CHKERRQ(ierr);
  ierr = DMDAVecRestoreArray(user->da,X,&p);CHKERRQ(ierr);
  ierr = VecSum(X,&sum);CHKERRQ(ierr);
  ierr = VecScale(X,1.0/sum);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

开发者ID:feelpp，项目名称:debian-petsc，代码行数:56，代码来源:ex7.c

示例11: TSAlphaAdaptDefault

PetscErrorCode  TSAlphaAdaptDefault(TS ts,PetscReal t,Vec X,Vec Xdot, PetscReal *nextdt,PetscBool *ok,void *ctx)
{
  TS_Alpha            *th;
  SNESConvergedReason snesreason;
  PetscReal           dt,normX,normE,Emax,scale;
  PetscErrorCode      ierr;
  PetscFunctionBegin;

  PetscValidHeaderSpecific(ts,TS_CLASSID,1);
#if PETSC_USE_DEBUG
  {
    PetscBool match;
    ierr = PetscObjectTypeCompare((PetscObject)ts,TSALPHA,&match);CHKERRQ(ierr);
    if (!match) SETERRQ(((PetscObject)ts)->comm,1,"Only for TSALPHA");
  }
#endif
  th = (TS_Alpha*)ts->data;

  ierr = SNESGetConvergedReason(ts->snes,&snesreason);CHKERRQ(ierr);
  if (snesreason < 0) {
    *ok = PETSC_FALSE;
    *nextdt *= th->scale_min;
    goto finally;
  }

  /* first-order aproximation to the local error */
  /* E = (X0 + dt*Xdot) - X */
  ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
  if (!th->E) {ierr = VecDuplicate(th->X0,&th->E);CHKERRQ(ierr);}
  ierr = VecWAXPY(th->E,dt,Xdot,th->X0);CHKERRQ(ierr);
  ierr = VecAXPY(th->E,-1,X);CHKERRQ(ierr);
  ierr = VecNorm(th->E,NORM_2,&normE);CHKERRQ(ierr);
  /* compute maximum allowable error */
  ierr = VecNorm(X,NORM_2,&normX);CHKERRQ(ierr);
  if (normX == 0) {ierr = VecNorm(th->X0,NORM_2,&normX);CHKERRQ(ierr);}
  Emax =  th->rtol * normX + th->atol;
  /* compute next time step */
  if (normE > 0) {
    scale = th->rho * PetscRealPart(PetscSqrtScalar((PetscScalar)(Emax/normE)));
    scale = PetscMax(scale,th->scale_min);
    scale = PetscMin(scale,th->scale_max);
    if (!(*ok))
      scale = PetscMin(1.0,scale);
    *nextdt *= scale;
  }
  /* accept or reject step */
  if (normE <= Emax)
    *ok = PETSC_TRUE;
  else
    *ok = PETSC_FALSE;

  finally:
  *nextdt = PetscMax(*nextdt,th->dt_min);
  *nextdt = PetscMin(*nextdt,th->dt_max);
  PetscFunctionReturn(0);
}

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

示例12: FormFunctionLocal

PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **t,PetscScalar **f,void *ptr)
{
  PetscInt     i,j;
  PetscScalar  hx,hy;
  PetscScalar  gradup,graddown,gradleft,gradright,gradx,grady;
  PetscScalar  coeffup,coeffdown,coeffleft,coeffright;

  PetscFunctionBeginUser;
  hx    = 1.0/(PetscReal)(info->mx-1);  hy    = 1.0/(PetscReal)(info->my-1);

  /* Evaluate function */
  for (j=info->ys; j<info->ys+info->ym; j++) {
    for (i=info->xs; i<info->xs+info->xm; i++) {

      if (i == 0 || i == info->mx-1 || j == 0 || j == info->my-1) {

        f[j][i] = t[j][i] - (1.0 - (2.0*hx*(PetscReal)i - 1.0)*(2.0*hx*(PetscReal)i - 1.0));

      } else {

        gradup     = (t[j+1][i] - t[j][i])/hy;
        graddown   = (t[j][i] - t[j-1][i])/hy;
        gradright  = (t[j][i+1] - t[j][i])/hx;
        gradleft   = (t[j][i] - t[j][i-1])/hx;

        gradx      = .5*(t[j][i+1] - t[j][i-1])/hx;
        grady      = .5*(t[j+1][i] - t[j-1][i])/hy;

        coeffup    = 1.0/PetscSqrtScalar(1.0 + gradup*gradup + gradx*gradx);
        coeffdown  = 1.0/PetscSqrtScalar(1.0 + graddown*graddown + gradx*gradx);

        coeffleft  = 1.0/PetscSqrtScalar(1.0 + gradleft*gradleft + grady*grady);
        coeffright = 1.0/PetscSqrtScalar(1.0 + gradright*gradright + grady*grady);

        f[j][i] = (coeffup*gradup - coeffdown*graddown)*hx + (coeffright*gradright - coeffleft*gradleft)*hy;

      }

    }
  }
  PetscFunctionReturn(0);
}

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

示例13: Rotation

PetscErrorCode Rotation(Coor x0, Coor x1, Coor *n, Coor *r, Coor *s)
{
  PetscReal mag;
  n->x = x1.x - x0.x;
  n->y = x1.y - x0.y;
  n->z = x1.z - x0.z;
  mag = PetscSqrtScalar(n->x*n->x + n->y*n->y + n->z*n->z);
  n->x /= mag;
  n->y /= mag;
  n->z /= mag;

  PetscReal nx2 = n->x*n->x,
            ny2 = n->y*n->y,
            nz2 = n->z*n->z;
  if( nx2+ ny2 > nx2 + nz2 ) {
    mag = PetscSqrtScalar(nx2 + ny2);
    r->x =  n->y / mag;
    r->y = -n->x / mag;
    r->z = 0;

    s->x = n->x * n->z;
    s->y = n->y * n->z;
    s->z = -nx2-ny2;
  } else {
    mag = PetscSqrtScalar(nx2 + nz2);
    r->x =  n->z / mag;
    r->y = 0;
    r->z = -n->x / mag;

    s->x = -n->x * n->y;
    s->y = nx2 + nz2;
    s->z = -n->y * n->z;
  }
  mag = PetscSqrtScalar(s->x*s->x + s->y*s->y + s->z*s->z);
  s->x /= mag;
  s->y /= mag;
  s->z /= mag;

  return 0;
}

开发者ID:adrielb，项目名称:DCell，代码行数:40，代码来源:Fiber.c

示例14: IPMComputeKKT

static PetscErrorCode IPMComputeKKT(Tao tao)
{
  TAO_IPM        *ipmP = (TAO_IPM *)tao->data;
  PetscScalar    norm;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = VecCopy(tao->gradient,ipmP->rd);CHKERRQ(ierr);

  if (ipmP->me > 0) {
    /* rd = gradient + Ae'*lamdae */
    ierr = MatMultTranspose(tao->jacobian_equality,ipmP->lamdae,ipmP->work);CHKERRQ(ierr);
    ierr = VecAXPY(ipmP->rd, 1.0, ipmP->work);CHKERRQ(ierr);

    /* rpe = ce(x) */
    ierr = VecCopy(tao->constraints_equality,ipmP->rpe);CHKERRQ(ierr);
  }
  if (ipmP->nb > 0) {
    /* rd = rd - Ai'*lamdai */
    ierr = MatMultTranspose(ipmP->Ai,ipmP->lamdai,ipmP->work);CHKERRQ(ierr);
    ierr = VecAXPY(ipmP->rd, -1.0, ipmP->work);CHKERRQ(ierr);

    /* rpi = cin - s */
    ierr = VecCopy(ipmP->ci,ipmP->rpi);CHKERRQ(ierr);
    ierr = VecAXPY(ipmP->rpi, -1.0, ipmP->s);CHKERRQ(ierr);

    /* com = s .* lami */
    ierr = VecPointwiseMult(ipmP->complementarity, ipmP->s,ipmP->lamdai);CHKERRQ(ierr);
  }
  /* phi = ||rd; rpe; rpi; com|| */
  ierr = VecDot(ipmP->rd,ipmP->rd,&norm);CHKERRQ(ierr);
  ipmP->phi = norm;
  if (ipmP->me > 0 ) {
    ierr = VecDot(ipmP->rpe,ipmP->rpe,&norm);CHKERRQ(ierr);
    ipmP->phi += norm;
  }
  if (ipmP->nb > 0) {
    ierr = VecDot(ipmP->rpi,ipmP->rpi,&norm);CHKERRQ(ierr);
    ipmP->phi += norm;
    ierr = VecDot(ipmP->complementarity,ipmP->complementarity,&norm);CHKERRQ(ierr);
    ipmP->phi += norm;
    /* mu = s'*lami/nb */
    ierr = VecDot(ipmP->s,ipmP->lamdai,&ipmP->mu);CHKERRQ(ierr);
    ipmP->mu /= ipmP->nb;
  } else {
    ipmP->mu = 1.0;
  }

  ipmP->phi = PetscSqrtScalar(ipmP->phi);
  PetscFunctionReturn(0);
}

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

示例15: KSPGMRESUpdateHessenberg

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

  PetscFunctionBegin;
  hh = HH(0,it);
  cc = CC(0);
  ss = SS(0);

  /* Apply all the previously computed plane rotations to the new column
     of the Hessenberg matrix */
  for (j=1; j<=it; j++) {
    tt  = *hh;
    *hh = PetscConj(*cc) * tt + *ss * *(hh+1);
    hh++;
    *hh = *cc++ * *hh - (*ss++ * tt);
  }

  /*
    compute the new plane rotation, and apply it to:
     1) the right-hand-side of the Hessenberg system
     2) the new column of the Hessenberg matrix
    thus obtaining the updated value of the residual
  */
  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;
    *ss        = *(hh+1) / tt;
    *GRS(it+1) = -(*ss * *GRS(it));
    *GRS(it)   = PetscConj(*cc) * *GRS(it);
    *hh        = PetscConj(*cc) * *hh + *ss * *(hh+1);
    *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:feelpp，项目名称:debian-petsc，代码行数:50，代码来源:gmres.c

示例16: KSPFischerGuessUpdate_Method2

PetscErrorCode  KSPFischerGuessUpdate_Method2(KSPFischerGuess_Method2 *itg,Vec x)
{
  PetscScalar    norm;
  PetscErrorCode ierr;
  int            curl = itg->curl,i;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(x,VEC_CLASSID,2);
  PetscValidPointer(itg,3);
  if (curl == itg->maxl) {
    ierr      = KSP_MatMult(itg->ksp,itg->mat,x,itg->Ax);CHKERRQ(ierr); /* norm = sqrt(x'Ax) */
    ierr      = VecDot(x,itg->Ax,&norm);CHKERRQ(ierr);
    ierr      = VecCopy(x,itg->xtilde[0]);CHKERRQ(ierr);
    ierr      = VecScale(itg->xtilde[0],1.0/PetscSqrtScalar(norm));CHKERRQ(ierr);
    itg->curl = 1;
  } else {
    if (!curl) {
      ierr = VecCopy(x,itg->xtilde[curl]);CHKERRQ(ierr);
    } else {
      ierr = VecWAXPY(itg->xtilde[curl],-1.0,itg->guess,x);CHKERRQ(ierr);
    }
    ierr = KSP_MatMult(itg->ksp,itg->mat,itg->xtilde[curl],itg->Ax);CHKERRQ(ierr);
    ierr = VecMDot(itg->Ax,curl,itg->xtilde,itg->alpha);CHKERRQ(ierr);
    for (i=0; i<curl; i++) itg->alpha[i] = -itg->alpha[i];
    ierr = VecMAXPY(itg->xtilde[curl],curl,itg->alpha,itg->xtilde);CHKERRQ(ierr);

    ierr = KSP_MatMult(itg->ksp,itg->mat,itg->xtilde[curl],itg->Ax);CHKERRQ(ierr); /* norm = sqrt(xtilde[curl]'Axtilde[curl]) */
    ierr = VecDot(itg->xtilde[curl],itg->Ax,&norm);CHKERRQ(ierr);
    if (PetscAbsScalar(norm) != 0.0) {
      ierr = VecScale(itg->xtilde[curl],1.0/PetscSqrtScalar(norm));CHKERRQ(ierr);
      itg->curl++;
    } else {
      ierr = PetscInfo(itg->ksp,"Not increasing dimension of Fischer space because new direction is identical to previous\n");CHKERRQ(ierr);
    }
  }
  PetscFunctionReturn(0);
}

开发者ID:fengyuqi，项目名称:petsc，代码行数:37，代码来源:iguess.c

示例17: ComputeSensiP

PetscErrorCode ComputeSensiP(Vec lambda,Vec mu,AppCtx *ctx)
{ 
  PetscErrorCode    ierr;
  PetscScalar       sensip;
  const PetscScalar *x,*y;
  
  PetscFunctionBegin;
  ierr = VecGetArrayRead(lambda,&x);CHKERRQ(ierr);
  ierr = VecGetArrayRead(mu,&y);CHKERRQ(ierr);
  sensip = 1./PetscSqrtScalar(1.-(ctx->Pm/ctx->Pmax)*(ctx->Pm/ctx->Pmax))/ctx->Pmax*x[0]+y[0];  
  ierr = PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt parameter pm: %.7f \n",(double)sensip);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(lambda,&x);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(mu,&y);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

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

示例18: SpeedOfSound_PG

/*
 For the ideal gas, the relationship of the speed of the sound "c",
 density "\rho" and the pressure "p" is:
 p = c^2 \rho
*/
PetscErrorCode SpeedOfSound_PG(User user,const Node *x,PetscScalar *c)
{
  PetscErrorCode ierr;
  PetscScalar p;

  PetscFunctionBeginUser;
  if (user->includeenergy){
    ierr = Pressure_Full(user,x,&p);CHKERRQ(ierr);
  }else{
    ierr = Pressure_Partial(user,x,&p);CHKERRQ(ierr);
  }

  (*c)=PetscSqrtScalar(user->adiabatic*PetscAbsScalar(p)/x->r);

  PetscFunctionReturn(0);
}

开发者ID:rlchen2008，项目名称:FVM-Rlchen，代码行数:21，代码来源:SetupFunctions.c

示例19: ini_bou

PetscErrorCode ini_bou(Vec X,AppCtx* user)
{
    PetscErrorCode ierr;
    DM             cda;
    DMDACoor2d     **coors;
    PetscScalar    **p;
    Vec            gc;
    PetscInt       i,j;
    PetscInt       xs,ys,xm,ym,M,N;
    PetscScalar    xi,yi;
    PetscScalar    sigmax=user->sigmax,sigmay=user->sigmay;
    PetscScalar    rho   =user->rho;
    PetscScalar    mux   =user->mux,muy=user->muy;
    PetscMPIInt    rank;

    PetscFunctionBeginUser;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    CHKERRQ(ierr);
    ierr = DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
    user->dx = (user->xmax - user->xmin)/(M-1);
    user->dy = (user->ymax - user->ymin)/(N-1);
    ierr = DMGetCoordinateDM(user->da,&cda);
    CHKERRQ(ierr);
    ierr = DMGetCoordinates(user->da,&gc);
    CHKERRQ(ierr);
    ierr = DMDAVecGetArray(cda,gc,&coors);
    CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->da,X,&p);
    CHKERRQ(ierr);
    ierr = DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);
    CHKERRQ(ierr);
    for (i=xs; i < xs+xm; i++) {
        for (j=ys; j < ys+ym; j++) {
            xi = coors[j][i].x;
            yi = coors[j][i].y;
            if (i == 0 || j == 0 || i == M-1 || j == N-1) p[j][i] = 0.0;
            else p[j][i] = (0.5/(PETSC_PI*sigmax*sigmay*PetscSqrtScalar(1.0-rho*rho)))*PetscExpScalar(-0.5/(1-rho*rho)*(PetscPowScalar((xi-mux)/sigmax,2) + PetscPowScalar((yi-muy)/sigmay,2) - 2*rho*(xi-mux)*(yi-muy)/(sigmax*sigmay)));
        }
    }
    /*  p[N/2+N%2][M/2+M%2] = 1/(user->dx*user->dy); */

    ierr = DMDAVecRestoreArray(cda,gc,&coors);
    CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->da,X,&p);
    CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

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

示例20: FormInitialSolution

PetscErrorCode FormInitialSolution(Vec U,void* ptr)
{
  AppCtx         *user=(AppCtx*)ptr;
  DM             da=user->da;
  PetscReal      c=user->c;
  PetscErrorCode ierr;
  PetscInt       i,j,xs,ys,xm,ym,Mx,My;
  PetscScalar    **u;
  PetscReal      hx,hy,x,y,r;

  PetscFunctionBeginUser;
  ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
                   PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

  hx     = 1.0/(PetscReal)(Mx-1);
  hy     = 1.0/(PetscReal)(My-1);

  /* Get pointers to vector data */
  ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);

  /* Get local grid boundaries */
  ierr = DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(ierr);

  /* Compute function over the locally owned part of the grid */
  for (j=ys; j<ys+ym; j++) {
    y = j*hy;
    for (i=xs; i<xs+xm; i++) {
      x = i*hx;
      r = PetscSqrtScalar((x-.5)*(x-.5) + (y-.5)*(y-.5));
      if (r < .125) {
        u[j][i] = PetscExpScalar(c*r*r*r);
      } else {
        u[j][i] = 0.0;
      }
    }
  }

  /* Restore vectors */
  ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

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

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