本文整理汇总了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;未经允许,请勿转载。 |
请发表评论