C++ PetscSqrtScalar函数代码示例

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);
}

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);
}

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);
}

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

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);
}

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);
}

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);
}

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);
}

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);
}

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);
}

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);
}

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);
}

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;
}

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);
}

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);
}