C++ VecRestoreArray函数代码示例

//.........这里部分代码省略.........
        xt=user->top[i+1];
        xlt = user->top[i];
      }else {
        xt = x[row+mx];
      }

      if (i>0 && j+1<my){
        xlt = x[row-1+mx];
      }
      if (j>0 && i+1<mx){
        xrb = x[row+1-mx];
      }

      d1 = (xc-xl);
      d2 = (xc-xr);
      d3 = (xc-xt);
      d4 = (xc-xb);
      d5 = (xr-xrb);
      d6 = (xrb-xb);
      d7 = (xlt-xl);
      d8 = (xt-xlt);

      df1dxc = d1*hydhx;
      df2dxc = ( d1*hydhx + d4*hxdhy );
      df3dxc = d3*hxdhy;
      df4dxc = ( d2*hydhx + d3*hxdhy );
      df5dxc = d2*hydhx;
      df6dxc = d4*hxdhy;

      d1 *= rhx;
      d2 *= rhx;
      d3 *= rhy;
      d4 *= rhy;
      d5 *= rhy;
      d6 *= rhx;
      d7 *= rhy;
      d8 *= rhx;

      f1 = PetscSqrtScalar( 1.0 + d1*d1 + d7*d7);
      f2 = PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
      f3 = PetscSqrtScalar( 1.0 + d3*d3 + d8*d8);
      f4 = PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
      f5 = PetscSqrtScalar( 1.0 + d2*d2 + d5*d5);
      f6 = PetscSqrtScalar( 1.0 + d4*d4 + d6*d6);

      ft = ft + (f2 + f4);

      df1dxc /= f1;
      df2dxc /= f2;
      df3dxc /= f3;
      df4dxc /= f4;
      df5dxc /= f5;
      df6dxc /= f6;

      g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc )/2.0;
    }
  }

  for (j=0; j<my; j++){   /* left side */
    d3=(user->left[j+1] - user->left[j+2])*rhy;
    d2=(user->left[j+1] - x[j*mx])*rhx;
    ft = ft+PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
  }

  for (i=0; i<mx; i++){ /* bottom */
    d2=(user->bottom[i+1]-user->bottom[i+2])*rhx;
    d3=(user->bottom[i+1]-x[i])*rhy;
    ft = ft+PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
  }

  for (j=0; j< my; j++){ /* right side */
    d1=(x[(j+1)*mx-1]-user->right[j+1])*rhx;
    d4=(user->right[j]-user->right[j+1])*rhy;
    ft = ft+PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
  }

  for (i=0; i<mx; i++){ /* top side */
    d1=(x[(my-1)*mx + i] - user->top[i+1])*rhy;
    d4=(user->top[i+1] - user->top[i])*rhx;
    ft = ft+PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
  }

  /* Bottom left corner */
  d1=(user->left[0]-user->left[1])*rhy;
  d2=(user->bottom[0]-user->bottom[1])*rhx;
  ft +=PetscSqrtScalar( 1.0 + d1*d1 + d2*d2);

  /* Top right corner */
  d1=(user->right[my+1] - user->right[my])*rhy;
  d2=(user->top[mx+1] - user->top[mx])*rhx;
  ft +=PetscSqrtScalar( 1.0 + d1*d1 + d2*d2);

  (*fcn)=ft*area;

  /* Restore vectors */
  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
  ierr = VecRestoreArray(G,&g);CHKERRQ(ierr);
  ierr = PetscLogFlops(67*mx*my);CHKERRQ(ierr);
  return 0;
}

int main(int argc,char **argv)
{
  Vec                p;
  PetscScalar        *x_ptr;
  PetscErrorCode     ierr;
  PetscMPIInt        size;
  AppCtx             ctx;
  Vec                lowerb,upperb;
  Tao                tao;
  TaoConvergedReason reason;
  KSP                ksp;
  PC                 pc;

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Initialize program
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  PetscInitialize(&argc,&argv,NULL,help);
  PetscFunctionBeginUser;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    Set runtime options
    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Swing equation options","");CHKERRQ(ierr);
  {
    ctx.beta    = 2;
    ctx.c       = 10000.0;
    ctx.u_s     = 1.0;
    ctx.omega_s = 1.0;
    ctx.omega_b = 120.0*PETSC_PI;
    ctx.H       = 5.0;
    ierr        = PetscOptionsScalar("-Inertia","","",ctx.H,&ctx.H,NULL);CHKERRQ(ierr);
    ctx.D       = 5.0;
    ierr        = PetscOptionsScalar("-D","","",ctx.D,&ctx.D,NULL);CHKERRQ(ierr);
    ctx.E       = 1.1378;
    ctx.V       = 1.0;
    ctx.X       = 0.545;
    ctx.Pmax    = ctx.E*ctx.V/ctx.X;;
    ierr        = PetscOptionsScalar("-Pmax","","",ctx.Pmax,&ctx.Pmax,NULL);CHKERRQ(ierr);
    ctx.Pm     = 0.4;
    ierr        = PetscOptionsScalar("-Pm","","",ctx.Pm,&ctx.Pm,NULL);CHKERRQ(ierr);
    ctx.tf      = 0.1;
    ctx.tcl     = 0.2;
    ierr        = PetscOptionsReal("-tf","Time to start fault","",ctx.tf,&ctx.tf,NULL);CHKERRQ(ierr);
    ierr        = PetscOptionsReal("-tcl","Time to end fault","",ctx.tcl,&ctx.tcl,NULL);CHKERRQ(ierr);

  }
  ierr = PetscOptionsEnd();CHKERRQ(ierr);

  /* Create TAO solver and set desired solution method */
  ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr);
  ierr = TaoSetType(tao,TAOBLMVM);CHKERRQ(ierr);

  /*
     Optimization starts
  */
  /* Set initial solution guess */
  ierr = VecCreateSeq(PETSC_COMM_WORLD,1,&p);CHKERRQ(ierr);
  ierr = VecGetArray(p,&x_ptr);CHKERRQ(ierr);
  x_ptr[0]   = ctx.Pm;
  ierr = VecRestoreArray(p,&x_ptr);CHKERRQ(ierr);

  ierr = TaoSetInitialVector(tao,p);CHKERRQ(ierr);
  /* Set routine for function and gradient evaluation */
  ierr = TaoSetObjectiveRoutine(tao,FormFunction,(void *)&ctx);CHKERRQ(ierr);
  ierr = TaoSetGradientRoutine(tao,TaoDefaultComputeGradient,(void *)&ctx);CHKERRQ(ierr);

  /* Set bounds for the optimization */
  ierr = VecDuplicate(p,&lowerb);CHKERRQ(ierr);
  ierr = VecDuplicate(p,&upperb);CHKERRQ(ierr);
  ierr = VecGetArray(lowerb,&x_ptr);CHKERRQ(ierr);
  x_ptr[0] = 0.;
  ierr = VecRestoreArray(lowerb,&x_ptr);CHKERRQ(ierr);
  ierr = VecGetArray(upperb,&x_ptr);CHKERRQ(ierr);
  x_ptr[0] = 1.1;;
  ierr = VecRestoreArray(upperb,&x_ptr);CHKERRQ(ierr);
  ierr = TaoSetVariableBounds(tao,lowerb,upperb);

  /* Check for any TAO command line options */
  ierr = TaoSetFromOptions(tao);CHKERRQ(ierr);
  ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr);
  if (ksp) {
    ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
    ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
  }

  ierr = TaoSetTolerances(tao,1e-15,1e-15,1e-15,1e-15,1e-15);
  /* SOLVE THE APPLICATION */
  ierr = TaoSolve(tao); CHKERRQ(ierr);

  /* Get information on termination */
  ierr = TaoGetConvergedReason(tao,&reason);CHKERRQ(ierr);
  if (reason <= 0){
    ierr=PetscPrintf(MPI_COMM_WORLD, "Try another method! \n");CHKERRQ(ierr);
  }

  ierr = VecView(p,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  /* Free TAO data structures */
  ierr = TaoDestroy(&tao);CHKERRQ(ierr);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    ierr = VecDuplicate(user.xloc,&user.rloc);CHKERRQ(ierr);

    /* Create the scatter between the global x and local xloc */
    ierr = ISCreateStride(MPI_COMM_SELF,2,0,1,&islocal);CHKERRQ(ierr);
    ierr = ISCreateStride(MPI_COMM_SELF,2,0,1,&isglobal);CHKERRQ(ierr);
    ierr = VecScatterCreate(x,isglobal,user.xloc,islocal,&user.scatter);CHKERRQ(ierr);
    ierr = ISDestroy(&isglobal);CHKERRQ(ierr);
    ierr = ISDestroy(&islocal);CHKERRQ(ierr);
  }

  /*
     Create Jacobian matrix data structure
  */
  ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr);
  ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr);
  ierr = MatSetFromOptions(J);CHKERRQ(ierr);
  ierr = MatSetUp(J);CHKERRQ(ierr);

  ierr = PetscOptionsHasName(PETSC_NULL,"-hard",&flg);CHKERRQ(ierr);
  if (!flg) {
    /* 
     Set function evaluation routine and vector.
    */
    ierr = SNESSetFunction(snes,r,FormFunction1,&user);CHKERRQ(ierr);

    /* 
     Set Jacobian matrix data structure and Jacobian evaluation routine
    */
    ierr = SNESSetJacobian(snes,J,J,FormJacobian1,PETSC_NULL);CHKERRQ(ierr);
  } else {
    if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This case is a uniprocessor example only!");
    ierr = SNESSetFunction(snes,r,FormFunction2,PETSC_NULL);CHKERRQ(ierr);
    ierr = SNESSetJacobian(snes,J,J,FormJacobian2,PETSC_NULL);CHKERRQ(ierr);
  }

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Customize nonlinear solver; set runtime options
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /* 
     Set linear solver defaults for this problem. By extracting the
     KSP, KSP, and PC contexts from the SNES context, we can then
     directly call any KSP, KSP, and PC routines to set various options.
  */
  ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
  ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
  ierr = KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);CHKERRQ(ierr);

  /* 
     Set SNES/KSP/KSP/PC runtime options, e.g.,
         -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
     These options will override those specified above as long as
     SNESSetFromOptions() is called _after_ any other customization
     routines.
  */
  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Evaluate initial guess; then solve nonlinear system
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  if (!flg) {
    ierr = VecSet(x,pfive);CHKERRQ(ierr);
  } else {
    ierr = VecGetArray(x,&xx);CHKERRQ(ierr);
    xx[0] = 2.0; xx[1] = 3.0;
    ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr);
  }
  /*
     Note: The user should initialize the vector, x, with the initial guess
     for the nonlinear solver prior to calling SNESSolve().  In particular,
     to employ an initial guess of zero, the user should explicitly set
     this vector to zero by calling VecSet().
  */

  ierr = SNESSolve(snes,PETSC_NULL,x);CHKERRQ(ierr);
  ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
  if (flg) {
    Vec f;
    ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = SNESGetFunction(snes,&f,0,0);CHKERRQ(ierr);
    ierr = VecView(r,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  }

  ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SNES iterations = %D\n",its);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Free work space.  All PETSc objects should be destroyed when they
     are no longer needed.
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&r);CHKERRQ(ierr);
  ierr = MatDestroy(&J);CHKERRQ(ierr); ierr = SNESDestroy(&snes);CHKERRQ(ierr);
  if (size > 1){
    ierr = VecDestroy(&user.xloc);CHKERRQ(ierr); 
    ierr = VecDestroy(&user.rloc);CHKERRQ(ierr);
    ierr = VecScatterDestroy(&user.scatter);CHKERRQ(ierr);
  }
  ierr = PetscFinalize();
  return 0;
}

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

  // Now set the activation ...
  unsigned int numSteps = (unsigned int)(ceil(( ti.stop - ti.start)/ti.step));

 // Vec tauVec;

  // PetscScalar ***tauArray;

  unsigned int paramTimeSteps = (unsigned int)(ceil(( (double)(numSteps))/ ((double)(2*numParams)) ));

  /*
  for (int b=0; b<numParams; b++) {
    std::vector<Vec> tau;
    unsigned int tBegin = paramTimeSteps*b;
    unsigned int tEnd   = tBegin + numSteps/2; // paramTimeSteps*(b+2);

    // std::cout << "For param " << b << ": Time step range is " << tBegin << " -> " << tEnd << std::endl; 
    for (unsigned int t=0; t<numSteps+1; t++) {
      double newTime = (dt*(t-tBegin)*numSteps)/((double)(paramTimeSteps));
     // double fff = 0.0;
      CHKERRQ( DACreateGlobalVector(da, &tauVec) );
      CHKERRQ( VecSet( tauVec, 0.0));

      if ( (t>=tBegin) && (t<=tEnd)) {
        CHKERRQ(DAVecGetArray(da, tauVec, &tauArray));
        for (int k = z; k < z + p ; k++) {
          for (int j = y; j < y + n; j++) {
            for (int i = x; i < x + m; i++) {
              acx = (i)*hx; acy = (j)*hx; acz = (k)*hx;
              tauArray[k][j][i] = sin(M_PI*newTime)*cos(2*M_PI*acx)*cos(2*M_PI*acy)*cos(2*M_PI*acz);
            }
          }
        }
        CHKERRQ( DAVecRestoreArray ( da, tauVec, &tauArray ) );
      }
      tau.push_back(tauVec);
    }
    fBasis.push_back(tau);
  }
  */
 // std::cout << "Finished setting basis" << std::endl;

  /*
  // Set initial velocity ...
  CHKERRQ(DAVecGetArray(da, initialVelocity, &solArray));

  for (int k = z; k < z + p ; k++) {
  for (int j = y; j < y + n; j++) {
  for (int i = x; i < x + m; i++) {
  acx = (i)*hx; acy = (j)*hx; acz = (k)*hx;
  solArray[k][j][i] = M_PI*cos(2*M_PI*acx)*cos(2*M_PI*acy)*cos(2*M_PI*acz);
  }
  }
  }
  CHKERRQ( DAVecRestoreArray ( da, initialVelocity, &solArray ) );
  */

  std::vector<Vec> newF;

  Vec alpha;
  PetscScalar *avec;

  VecCreateSeq(PETSC_COMM_SELF, numParams, &alpha);
  /*
  VecCreate(PETSC_COMM_WORLD, &alpha);
  VecSetSizes(alpha, numParams, PETSC_DECIDE);

/*
   FormInitialGuess - Forms initial approximation.

   Input Parameters:
   user - user-defined application context
   X - vector

   Output Parameter:
   X - vector
 */
int FormInitialGuess(AppCtx *user,Vec X)
{
  int         i,j,row,mx,my,ierr;
  PetscReal   one = 1.0,lambda,temp1,temp,hx,hy,hxdhy,hydhx,sc;
  PetscScalar *x;

  /*
      Process 0 has to wait for all other processes to get here
   before proceeding to write in the shared vector
  */
  ierr = PetscBarrier((PetscObject)X);CHKERRQ(ierr);
  if (user->rank) {
    /*
       All the non-busy processors have to wait here for process 0 to finish
       evaluating the function; otherwise they will start using the vector values
       before they have been computed
    */
    ierr = PetscBarrier((PetscObject)X);CHKERRQ(ierr);
    return 0;
  }

  mx = user->mx;               my = user->my;        lambda = user->param;
  hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
  sc = hx*hy*lambda;        hxdhy = hx/hy;            hydhx = hy/hx;

  temp1 = lambda/(lambda + one);

  /*
     Get a pointer to vector data.
       - For default PETSc vectors, VecGetArray() returns a pointer to
         the data array.  Otherwise, the routine is implementation dependent.
       - You MUST call VecRestoreArray() when you no longer need access to
         the array.
  */
  ierr = VecGetArray(X,&x);CHKERRQ(ierr);

  /*
     Compute initial guess over the locally owned part of the grid
  */
#pragma arl(4)
#pragma distinct (*x,*f)
#pragma no side effects (sqrt)
  for (j=0; j<my; j++) {
    temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
    for (i=0; i<mx; i++) {
      row = i + j*mx;
      if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
        x[row] = 0.0;
        continue;
      }
      x[row] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
    }
  }

  /*
     Restore vector
  */
  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);

  ierr = PetscBarrier((PetscObject)X);CHKERRQ(ierr);
  return 0;
}

void Stokes_SLE_PenaltySolver_MakePenalty( Stokes_SLE_PenaltySolver* self, Stokes_SLE* sle, Vec* _penalty ) {
    Vec fVec = sle->fForceVec->vector, hVec = sle->hForceVec->vector, penalty, lambda;
    Mat kMat = sle->kStiffMat->matrix;
    FeMesh *mesh = sle->kStiffMat->rowVariable->feMesh;
    FeVariable *velField = sle->kStiffMat->rowVariable;
    SolutionVector* uVec = sle->uSolnVec;
    FeEquationNumber *eqNum = uVec->eqNum;
    IArray *inc;
    PetscScalar *lambdaVals, lambdaMin, *penaltyVals;
    int numDofs, numLocalElems, nodeCur, numLocalNodes, rank, eq;
    SolutionVector *solVec = sle->uSolnVec;
    double *velBackup;
    Vec vecBackup;
    int ii, jj, kk;

    MPI_Comm_rank( MPI_COMM_WORLD, &rank );

    numDofs = Mesh_GetDimSize( mesh );
    numLocalElems = FeMesh_GetElementLocalSize( mesh );
    numLocalNodes = FeMesh_GetNodeLocalSize( mesh );

    velBackup = (double*)malloc( numLocalNodes*numDofs*sizeof(double) );
    for( ii = 0; ii < numLocalNodes; ii++ )
        FeVariable_GetValueAtNode( velField, ii, velBackup + ii*numDofs );

    VecDuplicate( hVec, &penalty );
    VecGetArray( penalty, &penaltyVals );

    VecDuplicate( fVec, &lambda );
    MatGetDiagonal( kMat, lambda );
    {
        PetscInt idx;
        PetscReal min, max;

        VecMin( lambda, &idx, &min );
        VecMax( lambda, &idx, &max );
        if( rank == 0 ) {
           printf( "LAMBDA RANGE:\n" );
           printf( "  MIN: %e\n", min );
           printf( "  MAX: %e\n", max );
        }
    }

    vecBackup = solVec->vector;
    solVec->vector = lambda;
    SolutionVector_UpdateSolutionOntoNodes( solVec );

    inc = IArray_New();
    lambdaVals = (double*)malloc( numDofs*sizeof(double) );

    for( ii = 0; ii < numLocalElems; ii++ ) {

        lambdaMin = DBL_MAX;

        FeMesh_GetElementNodes( mesh, ii, inc );
        for( jj = 0; jj < inc->size; jj++ ) {

            nodeCur = inc->ptr[jj];
            FeVariable_GetValueAtNode( velField, nodeCur, lambdaVals );

            for( kk = 0; kk < numDofs; kk++ ) {

                eq = eqNum->mapNodeDof2Eq[nodeCur][kk];
                if( eq == -1 )
                    continue;

/*
                eq = *(int*)STreeMap_Map( eqNum->ownedMap, &eq );

                VecGetValues( lambda, 1, &eq, &lambdaVal );
*/

                if( lambdaVals[kk] < 0.0 )
                    printf( "%g\n",  lambdaVals[kk] );
                if( lambdaVals[kk] < lambdaMin )
                    lambdaMin = lambdaVals[kk];

            }
        }

        penaltyVals[ii] = lambdaMin;

    }

    if( lambdaVals ) free( lambdaVals );
    Stg_Class_Delete( inc );

    solVec->vector = vecBackup;

    for( ii = 0; ii < numLocalNodes; ii++ )
        FeVariable_SetValueAtNode( velField, ii, velBackup + ii*numDofs );
    if( velBackup ) free( velBackup );
    FeVariable_SyncShadowValues( velField );

    Stg_VecDestroy(&lambda );

    VecRestoreArray( penalty, &penaltyVals );
    VecAssemblyBegin( penalty );
    VecAssemblyEnd( penalty );

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

PetscErrorCode MatSolve_SeqSpooles(Mat A,Vec b,Vec x)
{
  Mat_Spooles      *lu = (Mat_Spooles*)A->spptr;
  PetscScalar      *array;
  DenseMtx         *mtxY, *mtxX ;
  PetscErrorCode   ierr;
  PetscInt         irow,neqns=A->cmap->n,nrow=A->rmap->n,*iv;
#if defined(PETSC_USE_COMPLEX)
  double           x_real,x_imag;
#else
  double           *entX;
#endif

  PetscFunctionBegin;
  mtxY = DenseMtx_new();
  DenseMtx_init(mtxY, lu->options.typeflag, 0, 0, nrow, 1, 1, nrow); /* column major */
  ierr = VecGetArray(b,&array);CHKERRQ(ierr);

  if (lu->options.useQR) {   /* copy b to mtxY */
    for ( irow = 0 ; irow < nrow; irow++ )  
#if !defined(PETSC_USE_COMPLEX)
      DenseMtx_setRealEntry(mtxY, irow, 0, *array++); 
#else
      DenseMtx_setComplexEntry(mtxY, irow, 0, PetscRealPart(array[irow]), PetscImaginaryPart(array[irow]));
#endif
  } else {                   /* copy permuted b to mtxY */
    iv = IV_entries(lu->oldToNewIV); 
    for ( irow = 0 ; irow < nrow; irow++ ) 
#if !defined(PETSC_USE_COMPLEX)
      DenseMtx_setRealEntry(mtxY, *iv++, 0, *array++); 
#else
      DenseMtx_setComplexEntry(mtxY,*iv++,0,PetscRealPart(array[irow]),PetscImaginaryPart(array[irow]));
#endif
  }
  ierr = VecRestoreArray(b,&array);CHKERRQ(ierr);

  mtxX = DenseMtx_new();
  DenseMtx_init(mtxX, lu->options.typeflag, 0, 0, neqns, 1, 1, neqns);
  if (lu->options.useQR) {
    FrontMtx_QR_solve(lu->frontmtx, lu->mtxA, mtxX, mtxY, lu->mtxmanager,
                  lu->cpus, lu->options.msglvl, lu->options.msgFile);
  } else {
    FrontMtx_solve(lu->frontmtx, mtxX, mtxY, lu->mtxmanager, 
                 lu->cpus, lu->options.msglvl, lu->options.msgFile);
  }
  if ( lu->options.msglvl > 2 ) {
    int err;
    ierr = PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n right hand side matrix after permutation");CHKERRQ(ierr);
    DenseMtx_writeForHumanEye(mtxY, lu->options.msgFile); 
    ierr = PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n solution matrix in new ordering");CHKERRQ(ierr);
    DenseMtx_writeForHumanEye(mtxX, lu->options.msgFile);
    err = fflush(lu->options.msgFile);
    if (err) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SYS,"fflush() failed on file");    
  }

  /* permute solution into original ordering, then copy to x */  
  DenseMtx_permuteRows(mtxX, lu->newToOldIV);
  ierr = VecGetArray(x,&array);CHKERRQ(ierr); 

#if !defined(PETSC_USE_COMPLEX)
  entX = DenseMtx_entries(mtxX);
  DVcopy(neqns, array, entX);
#else
  for (irow=0; irow<nrow; irow++){
    DenseMtx_complexEntry(mtxX,irow,0,&x_real,&x_imag);
    array[irow] = x_real+x_imag*PETSC_i;   
  }
#endif

  ierr = VecRestoreArray(x,&array);CHKERRQ(ierr);
  
  /* free memory */
  DenseMtx_free(mtxX);
  DenseMtx_free(mtxY);
  PetscFunctionReturn(0);
}

int main(int argc,char **args)
{
  typedef enum {RANDOM, CONSTANT, TANH, NUM_FUNCS} FuncType;
  const char      *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
  PetscMPIInt     size;
  int             n = 10,N,Ny,ndim=4,i,dim[4],DIM;
  Vec             x,y,z;
  PetscScalar     s;
  PetscRandom     rdm;
  PetscReal       enorm;
  PetscInt        func     = RANDOM;
  FuncType        function = RANDOM;
  PetscBool       view     = PETSC_FALSE;
  PetscErrorCode  ierr;
  PetscScalar     *x_array,*y_array,*z_array;
  fftw_plan       fplan,bplan;

  ierr = PetscInitialize(&argc,&args,(char*)0,help);CHKERRQ(ierr);
#if defined(PETSC_USE_COMPLEX)
  SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires real numbers");
#endif

  ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);
  if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This is a uniprocessor example only!");
  ierr     = PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex142");CHKERRQ(ierr);
  ierr     = PetscOptionsEList("-function", "Function type", "ex142", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);CHKERRQ(ierr);
  ierr     = PetscOptionsBool("-vec_view draw", "View the functions", "ex142", view, &view, NULL);CHKERRQ(ierr);
  function = (FuncType) func;
  ierr     = PetscOptionsEnd();CHKERRQ(ierr);

  for (DIM = 0; DIM < ndim; DIM++) {
    dim[DIM] = n;  /* size of real space vector in DIM-dimension */
  }
  ierr = PetscRandomCreate(PETSC_COMM_SELF, &rdm);CHKERRQ(ierr);
  ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr);

  for (DIM = 1; DIM < 5; DIM++) {
    /* create vectors of length N=dim[0]*dim[1]* ...*dim[DIM-1] */
    /*----------------------------------------------------------*/
    N = Ny = 1;
    for (i = 0; i < DIM-1; i++) {
      N *= dim[i];
    }
    Ny = N; Ny *= 2*(dim[DIM-1]/2 + 1); /* add padding elements to output vector y */
    N *= dim[DIM-1];


    ierr = PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n",DIM,N);CHKERRQ(ierr);
    ierr = VecCreateSeq(PETSC_COMM_SELF,N,&x);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject) x, "Real space vector");CHKERRQ(ierr);

    ierr = VecCreateSeq(PETSC_COMM_SELF,Ny,&y);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject) y, "Frequency space vector");CHKERRQ(ierr);

    ierr = VecDuplicate(x,&z);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject) z, "Reconstructed vector");CHKERRQ(ierr);

    /* Set fftw plan                    */
    /*----------------------------------*/
    ierr = VecGetArray(x,&x_array);CHKERRQ(ierr);
    ierr = VecGetArray(y,&y_array);CHKERRQ(ierr);
    ierr = VecGetArray(z,&z_array);CHKERRQ(ierr);

    unsigned int flags = FFTW_ESTIMATE; /*or FFTW_MEASURE */
    /* The data in the in/out arrays is overwritten during FFTW_MEASURE planning, so such planning
     should be done before the input is initialized by the user. */
    ierr = PetscPrintf(PETSC_COMM_SELF,"DIM: %d, N %d, Ny %d\n",DIM,N,Ny);CHKERRQ(ierr);

    switch (DIM) {
    case 1:
      fplan = fftw_plan_dft_r2c_1d(dim[0], (double*)x_array, (fftw_complex*)y_array, flags);
      bplan = fftw_plan_dft_c2r_1d(dim[0], (fftw_complex*)y_array, (double*)z_array, flags);
      break;
    case 2:
      fplan = fftw_plan_dft_r2c_2d(dim[0],dim[1],(double*)x_array, (fftw_complex*)y_array,flags);
      bplan = fftw_plan_dft_c2r_2d(dim[0],dim[1],(fftw_complex*)y_array,(double*)z_array,flags);
      break;
    case 3:
      fplan = fftw_plan_dft_r2c_3d(dim[0],dim[1],dim[2],(double*)x_array, (fftw_complex*)y_array,flags);
      bplan = fftw_plan_dft_c2r_3d(dim[0],dim[1],dim[2],(fftw_complex*)y_array,(double*)z_array,flags);
      break;
    default:
      fplan = fftw_plan_dft_r2c(DIM,(int*)dim,(double*)x_array, (fftw_complex*)y_array,flags);
      bplan = fftw_plan_dft_c2r(DIM,(int*)dim,(fftw_complex*)y_array,(double*)z_array,flags);
      break;
    }

    ierr = VecRestoreArray(x,&x_array);CHKERRQ(ierr);
    ierr = VecRestoreArray(y,&y_array);CHKERRQ(ierr);
    ierr = VecRestoreArray(z,&z_array);CHKERRQ(ierr);

    /* Initialize Real space vector x:
       The data in the in/out arrays is overwritten during FFTW_MEASURE planning, so planning
       should be done before the input is initialized by the user.
    --------------------------------------------------------*/
    if (function == RANDOM) {
      ierr = VecSetRandom(x, rdm);CHKERRQ(ierr);
    } else if (function == CONSTANT) {
      ierr = VecSet(x, 1.0);CHKERRQ(ierr);
    } else if (function == TANH) {
//.........这里部分代码省略.........

PetscErrorCode SetInitialGuess(Vec X,AppCtx *user)
{
  PetscErrorCode    ierr;
  PetscInt          nele,nen,n,i;
  const PetscInt    *ele;
  Vec               coords, rand1, rand2;
  const PetscScalar *_coords;
  PetscScalar       x[3],y[3];
  PetscInt          idx[3];
  PetscScalar       *xx,*w1,*w2,*u1,*u2,*u3;
  PetscViewer       view_out;

  PetscFunctionBeginUser;
  /* Get ghosted coordinates */
  ierr = DMGetCoordinatesLocal(user->da,&coords);CHKERRQ(ierr);
  ierr = VecDuplicate(user->u1,&rand1);
  ierr = VecDuplicate(user->u1,&rand2);
  ierr = VecSetRandom(rand1,NULL);
  ierr = VecSetRandom(rand2,NULL);

  ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr);
  ierr = VecGetArrayRead(coords,&_coords);CHKERRQ(ierr);
  ierr = VecGetArray(X,&xx);CHKERRQ(ierr);
  ierr = VecGetArray(user->work1,&w1);
  ierr = VecGetArray(user->work2,&w2);
  ierr = VecGetArray(user->u1,&u1);
  ierr = VecGetArray(user->u2,&u2);
  ierr = VecGetArray(user->u3,&u3);

  /* Get local element info */
  ierr = DMDAGetElements(user->da,&nele,&nen,&ele);CHKERRQ(ierr);
  for (i=0; i < nele; i++) {
    idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
    x[0]   = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
    x[1]   = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
    x[2]   = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];

    PetscScalar vals1[3],vals2[3],valsrand[3];
    PetscInt    r;
    for (r=0; r<3; r++) {
      valsrand[r]=5*x[r]*(1-x[r])*y[r]*(1-y[r]);
      if (x[r]>=0.5 && y[r]>=0.5) {
        vals1[r]=0.75;
        vals2[r]=0.0;
      }
      if (x[r]>=0.5 && y[r]<0.5) {
        vals1[r]=0.0;
        vals2[r]=0.0;
      }
      if (x[r]<0.5 && y[r]>=0.5) {
        vals1[r]=0.0;
        vals2[r]=0.75;
      }
      if (x[r]<0.5 && y[r]<0.5) {
        vals1[r]=0.75;
        vals2[r]=0.0;
      }
    }

    ierr = VecSetValues(user->work1,3,idx,vals1,INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecSetValues(user->work2,3,idx,vals2,INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecSetValues(user->work3,3,idx,valsrand,INSERT_VALUES);CHKERRQ(ierr);
  }

  ierr = VecAssemblyBegin(user->work1);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(user->work1);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(user->work2);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(user->work2);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(user->work3);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(user->work3);CHKERRQ(ierr);

  ierr = VecAXPY(user->work1,1.0,user->work3);CHKERRQ(ierr);
  ierr = VecAXPY(user->work2,1.0,user->work3);CHKERRQ(ierr);

  for (i=0; i<n/4; i++) {
    xx[4*i] = w1[i];
    if (xx[4*i]>1) xx[4*i]=1;

    xx[4*i+1] = w2[i];
    if (xx[4*i+1]>1) xx[4*i+1]=1;

    if (xx[4*i]+xx[4*i+1]>1) xx[4*i+1] = 1.0 - xx[4*i];

    xx[4*i+2] = 1.0 - xx[4*i] - xx[4*i+1];
    xx[4*i+3] = 0.0;

    u1[i] = xx[4*i];
    u2[i] = xx[4*i+1];
    u3[i] = xx[4*i+2];
  }

  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view_out);CHKERRQ(ierr);
  ierr = VecView(user->u1,view_out);CHKERRQ(ierr);
  ierr = VecView(user->u2,view_out);CHKERRQ(ierr);
  ierr = VecView(user->u3,view_out);CHKERRQ(ierr);
  PetscViewerDestroy(&view_out);

  ierr = DMDARestoreElements(user->da,&nele,&nen,&ele);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(coords,&_coords);CHKERRQ(ierr);
  ierr = VecRestoreArray(X,&xx);CHKERRQ(ierr);
//.........这里部分代码省略.........

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  PetscMPIInt    rank,size;
  PetscInt       rstart,rend,i,k,N,numPoints=1000000;
  PetscScalar    dummy,result=0,h=1.0/numPoints,*xarray;
  Vec            x,xend;

  ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);

  /*
     Create a parallel vector.
       Here we set up our x vector which will be given values below.
       The xend vector is a dummy vector to find the value of the
         elements at the endpoints for use in the trapezoid rule.
  */
  ierr   = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
  ierr   = VecSetSizes(x,PETSC_DECIDE,numPoints);CHKERRQ(ierr);
  ierr   = VecSetFromOptions(x);CHKERRQ(ierr);
  ierr   = VecGetSize(x,&N);CHKERRQ(ierr);
  ierr   = VecSet(x,result);CHKERRQ(ierr);
  ierr   = VecDuplicate(x,&xend);CHKERRQ(ierr);
  result = 0.5;
  if (!rank) {
    i    = 0;
    ierr = VecSetValues(xend,1,&i,&result,INSERT_VALUES);CHKERRQ(ierr);
  }
  if (rank == size-1) {
    i    = N-1;
    ierr = VecSetValues(xend,1,&i,&result,INSERT_VALUES);CHKERRQ(ierr);
  }
  /*
     Assemble vector, using the 2-step process:
       VecAssemblyBegin(), VecAssemblyEnd()
     Computations can be done while messages are in transition
     by placing code between these two statements.
  */
  ierr = VecAssemblyBegin(xend);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(xend);CHKERRQ(ierr);

  /*
     Set the x vector elements.
      i*h will return 0 for i=0 and 1 for i=N-1.
      The function evaluated (2x/(1+x^2)) is defined above.
      Each evaluation is put into the local array of the vector without message passing.
  */
  ierr = VecGetOwnershipRange(x,&rstart,&rend);CHKERRQ(ierr);
  ierr = VecGetArray(x,&xarray);CHKERRQ(ierr);
  k    = 0;
  for (i=rstart; i<rend; i++) {
    xarray[k] = i*h;
    xarray[k] = func(xarray[k]);
    k++;
  }
  ierr = VecRestoreArray(x,&xarray);CHKERRQ(ierr);

  /*
     Evaluates the integral.  First the sum of all the points is taken.
     That result is multiplied by the step size for the trapezoid rule.
     Then half the value at each endpoint is subtracted,
     this is part of the composite trapezoid rule.
  */
  ierr   = VecSum(x,&result);CHKERRQ(ierr);
  result = result*h;
  ierr   = VecDot(x,xend,&dummy);CHKERRQ(ierr);
  result = result-h*dummy;

  /*
      Return the value of the integral.
  */
  ierr = PetscPrintf(PETSC_COMM_WORLD,"ln(2) is %g\n",(double)result);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&xend);CHKERRQ(ierr);

  ierr = PetscFinalize();
  return ierr;
}

static PetscErrorCode PCSetUp_SVD(PC pc)
{
#if defined(PETSC_MISSING_LAPACK_GESVD)
  SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"GESVD - Lapack routine is unavailable\nNot able to provide singular value estimates.");
#else
  PC_SVD         *jac = (PC_SVD*)pc->data;
  PetscErrorCode ierr;
  PetscScalar    *a,*u,*v,*d,*work;
  PetscBLASInt   nb,lwork;
  PetscInt       i,n;
  PetscMPIInt    size;

  PetscFunctionBegin;
  ierr = MatDestroy(&jac->A);CHKERRQ(ierr);
  ierr = MPI_Comm_size(((PetscObject)pc->pmat)->comm,&size);CHKERRQ(ierr);
  if (size > 1) {
    Mat      redmat;
    PetscInt M;
    ierr = MatGetSize(pc->pmat,&M,NULL);CHKERRQ(ierr);
    ierr = MatGetRedundantMatrix(pc->pmat,size,PETSC_COMM_SELF,M,MAT_INITIAL_MATRIX,&redmat);CHKERRQ(ierr);
    ierr = MatConvert(redmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);CHKERRQ(ierr);
    ierr = MatDestroy(&redmat);CHKERRQ(ierr);
  } else {
    ierr = MatConvert(pc->pmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);CHKERRQ(ierr);
  }
  if (!jac->diag) {    /* assume square matrices */
    ierr = MatGetVecs(jac->A,&jac->diag,&jac->work);CHKERRQ(ierr);
  }
  if (!jac->U) {
    ierr = MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->U);CHKERRQ(ierr);
    ierr = MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->Vt);CHKERRQ(ierr);
  }
  ierr  = MatGetSize(pc->pmat,&n,NULL);CHKERRQ(ierr);
  ierr  = PetscBLASIntCast(n,&nb);CHKERRQ(ierr);
  lwork = 5*nb;
  ierr  = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr);
  ierr  = MatDenseGetArray(jac->A,&a);CHKERRQ(ierr);
  ierr  = MatDenseGetArray(jac->U,&u);CHKERRQ(ierr);
  ierr  = MatDenseGetArray(jac->Vt,&v);CHKERRQ(ierr);
  ierr  = VecGetArray(jac->diag,&d);CHKERRQ(ierr);
#if !defined(PETSC_USE_COMPLEX)
  {
    PetscBLASInt lierr;
    ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
    PetscStackCall("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,d,u,&nb,v,&nb,work,&lwork,&lierr));
    if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr);
    ierr = PetscFPTrapPop();CHKERRQ(ierr);
  }
#else
  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not coded for complex");
#endif
  ierr = MatDenseRestoreArray(jac->A,&a);CHKERRQ(ierr);
  ierr = MatDenseRestoreArray(jac->U,&u);CHKERRQ(ierr);
  ierr = MatDenseRestoreArray(jac->Vt,&v);CHKERRQ(ierr);
  for (i=n-1; i>=0; i--) if (PetscRealPart(d[i]) > jac->zerosing) break;
  jac->nzero = n-1-i;
  if (jac->monitor) {
    ierr = PetscViewerASCIIAddTab(jac->monitor,((PetscObject)pc)->tablevel);CHKERRQ(ierr);
    ierr = PetscViewerASCIIPrintf(jac->monitor,"    SVD: condition number %14.12e, %D of %D singular values are (nearly) zero\n",(double)PetscRealPart(d[0]/d[n-1]),jac->nzero,n);CHKERRQ(ierr);
    if (n >= 10) {              /* print 5 smallest and 5 largest */
      ierr = PetscViewerASCIIPrintf(jac->monitor,"    SVD: smallest singular values: %14.12e %14.12e %14.12e %14.12e %14.12e\n",(double)PetscRealPart(d[n-1]),(double)PetscRealPart(d[n-2]),(double)PetscRealPart(d[n-3]),(double)PetscRealPart(d[n-4]),(double)PetscRealPart(d[n-5]));CHKERRQ(ierr);
      ierr = PetscViewerASCIIPrintf(jac->monitor,"    SVD: largest singular values : %14.12e %14.12e %14.12e %14.12e %14.12e\n",(double)PetscRealPart(d[4]),(double)PetscRealPart(d[3]),(double)PetscRealPart(d[2]),(double)PetscRealPart(d[1]),(double)PetscRealPart(d[0]));CHKERRQ(ierr);
    } else {                    /* print all singular values */
      char     buf[256],*p;
      size_t   left = sizeof(buf),used;
      PetscInt thisline;
      for (p=buf,i=n-1,thisline=1; i>=0; i--,thisline++) {
        ierr  = PetscSNPrintfCount(p,left," %14.12e",&used,(double)PetscRealPart(d[i]));CHKERRQ(ierr);
        left -= used;
        p    += used;
        if (thisline > 4 || i==0) {
          ierr     = PetscViewerASCIIPrintf(jac->monitor,"    SVD: singular values:%s\n",buf);CHKERRQ(ierr);
          p        = buf;
          thisline = 0;
        }
      }
    }
    ierr = PetscViewerASCIISubtractTab(jac->monitor,((PetscObject)pc)->tablevel);CHKERRQ(ierr);
  }
  ierr = PetscInfo2(pc,"Largest and smallest singular values %14.12e %14.12e\n",(double)PetscRealPart(d[0]),(double)PetscRealPart(d[n-1]));CHKERRQ(ierr);
  for (i=0; i<n-jac->nzero; i++) d[i] = 1.0/d[i];
  for (; i<n; i++) d[i] = 0.0;
  if (jac->essrank > 0) for (i=0; i<n-jac->nzero-jac->essrank; i++) d[i] = 0.0; /* Skip all but essrank eigenvalues */
  ierr = PetscInfo1(pc,"Number of zero or nearly singular values %D\n",jac->nzero);CHKERRQ(ierr);
  ierr = VecRestoreArray(jac->diag,&d);CHKERRQ(ierr);
#if defined(foo)
  {
    PetscViewer viewer;
    ierr = PetscViewerBinaryOpen(PETSC_COMM_SELF,"joe",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
    ierr = MatView(jac->A,viewer);CHKERRQ(ierr);
    ierr = MatView(jac->U,viewer);CHKERRQ(ierr);
    ierr = MatView(jac->Vt,viewer);CHKERRQ(ierr);
    ierr = VecView(jac->diag,viewer);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
  }
#endif
  ierr = PetscFree(work);CHKERRQ(ierr);
  PetscFunctionReturn(0);
#endif
}

int main(int argc,char **argv)
{
  PetscMPIInt    rank,size;
  PetscInt       nlocal = 6,nghost = 2,ifrom[2],i,rstart,rend;
  PetscErrorCode ierr;
  PetscBool      flg,flg2;
  PetscScalar    value,*array,*tarray=0;
  Vec            lx,gx,gxs;

  PetscInitialize(&argc,&argv,(char*)0,help);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  if (size != 2) SETERRQ(PETSC_COMM_SELF,1,"Must run example with two processors\n");

  /*
     Construct a two dimensional graph connecting nlocal degrees of
     freedom per processor. From this we will generate the global
     indices of needed ghost values

     For simplicity we generate the entire graph on each processor:
     in real application the graph would stored in parallel, but this
     example is only to demonstrate the management of ghost padding
     with VecCreateGhost().

     In this example we consider the vector as representing
     degrees of freedom in a one dimensional grid with periodic
     boundary conditions.

        ----Processor  1---------  ----Processor 2 --------
         0    1   2   3   4    5    6    7   8   9   10   11
                               |----|
         |-------------------------------------------------|

  */

  if (!rank) {
    ifrom[0] = 11; ifrom[1] = 6;
  } else {
    ifrom[0] = 0;  ifrom[1] = 5;
  }

  /*
     Create the vector with two slots for ghost points. Note that both
     the local vector (lx) and the global vector (gx) share the same
     array for storing vector values.
  */
  ierr = PetscOptionsHasName(NULL,"-allocate",&flg);CHKERRQ(ierr);
  ierr = PetscOptionsHasName(NULL,"-vecmpisetghost",&flg2);CHKERRQ(ierr);
  if (flg) {
    ierr = PetscMalloc1((nlocal+nghost),&tarray);CHKERRQ(ierr);
    ierr = VecCreateGhostWithArray(PETSC_COMM_WORLD,nlocal,PETSC_DECIDE,nghost,ifrom,tarray,&gxs);CHKERRQ(ierr);
  } else if (flg2) {
    ierr = VecCreate(PETSC_COMM_WORLD,&gxs);CHKERRQ(ierr);
    ierr = VecSetType(gxs,VECMPI);CHKERRQ(ierr);
    ierr = VecSetSizes(gxs,nlocal,PETSC_DECIDE);CHKERRQ(ierr);
    ierr = VecMPISetGhost(gxs,nghost,ifrom);CHKERRQ(ierr);
  } else {
    ierr = VecCreateGhost(PETSC_COMM_WORLD,nlocal,PETSC_DECIDE,nghost,ifrom,&gxs);CHKERRQ(ierr);
  }

  /*
      Test VecDuplicate()
  */
  ierr = VecDuplicate(gxs,&gx);CHKERRQ(ierr);
  ierr = VecDestroy(&gxs);CHKERRQ(ierr);

  /*
     Access the local representation
  */
  ierr = VecGhostGetLocalForm(gx,&lx);CHKERRQ(ierr);

  /*
     Set the values from 0 to 12 into the "global" vector
  */
  ierr = VecGetOwnershipRange(gx,&rstart,&rend);CHKERRQ(ierr);
  for (i=rstart; i<rend; i++) {
    value = (PetscScalar) i;
    ierr  = VecSetValues(gx,1,&i,&value,INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = VecAssemblyBegin(gx);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(gx);CHKERRQ(ierr);

  ierr = VecGhostUpdateBegin(gx,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecGhostUpdateEnd(gx,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);

  /*
     Print out each vector, including the ghost padding region.
  */
  ierr = VecGetArray(lx,&array);CHKERRQ(ierr);
  for (i=0; i<nlocal+nghost; i++) {
    ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"%D %g\n",i,(double)PetscRealPart(array[i]));CHKERRQ(ierr);
  }
  ierr = VecRestoreArray(lx,&array);CHKERRQ(ierr);
  ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);

  ierr = VecGhostRestoreLocalForm(gx,&lx);CHKERRQ(ierr);
  ierr = VecDestroy(&gx);CHKERRQ(ierr);
  if (flg) {ierr = PetscFree(tarray);CHKERRQ(ierr);}
  ierr = PetscFinalize();
  return 0;
//.........这里部分代码省略.........

static PetscErrorCode PCApply_Kaczmarz(PC pc,Vec x,Vec y)
{
  PC_Kaczmarz       *jac = (PC_Kaczmarz*)pc->data;
  PetscInt          xs,xe,ys,ye,ncols,i,j;
  const PetscInt    *cols;
  const PetscScalar *vals;
  PetscErrorCode    ierr;
  PetscScalar       r;
  PetscReal         anrm;
  PetscScalar       *xarray,*yarray;
  PetscReal         lambda=jac->lambda;

  PetscFunctionBegin;
  ierr = MatGetOwnershipRange(pc->pmat,&xs,&xe);CHKERRQ(ierr);
  ierr = MatGetOwnershipRangeColumn(pc->pmat,&ys,&ye);CHKERRQ(ierr);
  ierr = VecSet(y,0.);CHKERRQ(ierr);
  ierr = VecGetArray(x,&xarray);CHKERRQ(ierr);
  ierr = VecGetArray(y,&yarray);CHKERRQ(ierr);
  for (i=xs;i<xe;i++) {
    /* get the maximum row width and row norms */
    ierr = MatGetRow(pc->pmat,i,&ncols,&cols,&vals);CHKERRQ(ierr);
    r = xarray[i-xs];
    anrm = 0.;
    for (j=0;j<ncols;j++) {
      if (cols[j] >= ys && cols[j] < ye) {
        r -= yarray[cols[j]-ys]*vals[j];
      }
      anrm += PetscRealPart(PetscSqr(vals[j]));
    }
    if (anrm > 0.) {
      for (j=0;j<ncols;j++) {
        if (cols[j] >= ys && cols[j] < ye) {
          yarray[cols[j]-ys] += vals[j]*lambda*r/anrm;
        }
      }
    }
    ierr = MatRestoreRow(pc->pmat,i,&ncols,&cols,&vals);CHKERRQ(ierr);
  }
  if (jac->symmetric) {
    for (i=xe-1;i>=xs;i--) {
      ierr = MatGetRow(pc->pmat,i,&ncols,&cols,&vals);CHKERRQ(ierr);
      r = xarray[i-xs];
      anrm = 0.;
      for (j=0;j<ncols;j++) {
        if (cols[j] >= ys && cols[j] < ye) {
          r -= yarray[cols[j]-ys]*vals[j];
        }
        anrm += PetscRealPart(PetscSqr(vals[j]));
      }
      if (anrm > 0.) {
        for (j=0;j<ncols;j++) {
          if (cols[j] >= ys && cols[j] < ye) {
            yarray[cols[j]-ys] += vals[j]*lambda*r/anrm;
          }
        }
      }
      ierr = MatRestoreRow(pc->pmat,i,&ncols,&cols,&vals);CHKERRQ(ierr);
    }
  }
  ierr = VecRestoreArray(y,&yarray);CHKERRQ(ierr);
  ierr = VecRestoreArray(x,&xarray);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

//.........这里部分代码省略.........
        xb = x[row-mx];
      }

      if (i+1 == mx){
        xr=user->right[j+1];
        xrb = user->right[j];
      } else {
        xr = x[row+1];
      }

      if (j+1==my){
        xt=user->top[i+1];
        xlt = user->top[i];
      }else {
        xt = x[row+mx];
      }

      if (i>0 && j+1<my){
        xlt = x[row-1+mx];
      }
      if (j>0 && i+1<mx){
        xrb = x[row+1-mx];
      }


      d1 = (xc-xl)*rhx;
      d2 = (xc-xr)*rhx;
      d3 = (xc-xt)*rhy;
      d4 = (xc-xb)*rhy;
      d5 = (xrb-xr)*rhy;
      d6 = (xrb-xb)*rhx;
      d7 = (xlt-xl)*rhy;
      d8 = (xlt-xt)*rhx;

      f1 = PetscSqrtScalar( 1.0 + d1*d1 + d7*d7);
      f2 = PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
      f3 = PetscSqrtScalar( 1.0 + d3*d3 + d8*d8);
      f4 = PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
      f5 = PetscSqrtScalar( 1.0 + d2*d2 + d5*d5);
      f6 = PetscSqrtScalar( 1.0 + d4*d4 + d6*d6);


      hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+(-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
      hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+(-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
      ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+(-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
      hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+(-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

      hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
      htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);

      hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) + hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
           (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +  (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

      hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5;  hc*=0.5;

      k=0;
      if (j>0){
        v[k]=hb; col[k]=row - mx; k++;
      }

      if (j>0 && i < mx -1){
        v[k]=hbr; col[k]=row - mx+1; k++;
      }

      if (i>0){
        v[k]= hl; col[k]=row - 1; k++;
      }

      v[k]= hc; col[k]=row; k++;

      if (i < mx-1 ){
        v[k]= hr; col[k]=row+1; k++;
      }

      if (i>0 && j < my-1 ){
        v[k]= htl; col[k] = row+mx-1; k++;
      }

      if (j < my-1 ){
        v[k]= ht; col[k] = row+mx; k++;
      }

      /*
         Set matrix values using local numbering, which was defined
         earlier, in the main routine.
      */
      ierr = MatSetValues(Hessian,1,&row,k,col,v,INSERT_VALUES);CHKERRQ(ierr);
    }
  }

  /* Restore vectors */
  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);

  /* Assemble the matrix */
  ierr = MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  ierr = PetscLogFlops(199*mx*my);CHKERRQ(ierr);
  return 0;
}

void UniqueVec::restoreVec(PetscScalar *local_array) {
  VecRestoreArray(*data_.get(), &local_array);
}