C++ VecGetSize函数代码示例

static PetscErrorCode TSInterpolate_Sundials(TS ts,PetscReal t,Vec X)
{
  TS_Sundials    *cvode = (TS_Sundials*)ts->data;
  N_Vector       y;
  PetscErrorCode ierr;
  PetscScalar    *x_data;
  PetscInt       glosize,locsize;

  PetscFunctionBegin;
  /* get the vector size */
  ierr = VecGetSize(X,&glosize);CHKERRQ(ierr);
  ierr = VecGetLocalSize(X,&locsize);CHKERRQ(ierr);

  /* allocate the memory for N_Vec y */
  y = N_VNew_Parallel(cvode->comm_sundials,locsize,glosize);
  if (!y) SETERRQ(PETSC_COMM_SELF,1,"Interpolated y is not allocated");

  ierr = VecGetArray(X,&x_data);CHKERRQ(ierr);
  N_VSetArrayPointer((realtype*)x_data,y);
  ierr = CVodeGetDky(cvode->mem,t,0,y);CHKERRQ(ierr);
  ierr = VecRestoreArray(X,&x_data);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode  SNESMonitorRange_Private(SNES snes,PetscInt it,PetscReal *per)
{
  PetscErrorCode ierr;
  Vec            resid;
  PetscReal      rmax,pwork;
  PetscInt       i,n,N;
  PetscScalar    *r;

  PetscFunctionBegin;
  ierr  = SNESGetFunction(snes,&resid,0,0);CHKERRQ(ierr);
  ierr  = VecNorm(resid,NORM_INFINITY,&rmax);CHKERRQ(ierr);
  ierr  = VecGetLocalSize(resid,&n);CHKERRQ(ierr);
  ierr  = VecGetSize(resid,&N);CHKERRQ(ierr);
  ierr  = VecGetArray(resid,&r);CHKERRQ(ierr);
  pwork = 0.0;
  for (i=0; i<n; i++) {
    pwork += (PetscAbsScalar(r[i]) > .20*rmax);
  }
  ierr = MPI_Allreduce(&pwork,per,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));CHKERRQ(ierr);
  ierr = VecRestoreArray(resid,&r);CHKERRQ(ierr);
  *per = *per/N;
  PetscFunctionReturn(0);
}

void ReplicatableVector::ReplicatePetscVector(Vec vec)
{
    // If the size has changed then we'll need to make a new context
    PetscInt isize;
    VecGetSize(vec, &isize);
    unsigned size = isize;

    if (this->GetSize() != size)
    {
        Resize(size);
    }
    if (mReplicated == NULL)
    {
        // This creates mToAll (the scatter context) and mReplicated (to store values)
        VecScatterCreateToAll(vec, &mToAll, &mReplicated);
    }

    // Replicate the data
//PETSc-3.x.x or PETSc-2.3.3
#if ( (PETSC_VERSION_MAJOR == 3) || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 3 && PETSC_VERSION_SUBMINOR == 3)) //2.3.3 or 3.x.x
    VecScatterBegin(mToAll, vec, mReplicated, INSERT_VALUES, SCATTER_FORWARD);
    VecScatterEnd  (mToAll, vec, mReplicated, INSERT_VALUES, SCATTER_FORWARD);
#else
//PETSc-2.3.2 or previous
    VecScatterBegin(vec, mReplicated, INSERT_VALUES, SCATTER_FORWARD, mToAll);
    VecScatterEnd  (vec, mReplicated, INSERT_VALUES, SCATTER_FORWARD, mToAll);
#endif

    // Information is now in mReplicated PETSc vector
    // Copy into mData
    double* p_replicated;
    VecGetArray(mReplicated, &p_replicated);
    for (unsigned i=0; i<size; i++)
    {
        mpData[i] = p_replicated[i];
    }
}

LData::LData(const std::string& name,
             Vec vec,
             const std::vector<int>& nonlocal_petsc_indices,
             const bool manage_petsc_vec)
    : d_name(name), d_global_node_count(0), d_local_node_count(0), d_ghost_node_count(0), d_depth(0),
      d_nonlocal_petsc_indices(nonlocal_petsc_indices), d_global_vec(vec), d_managing_petsc_vec(manage_petsc_vec),
      d_array(NULL), d_boost_array(NULL), d_boost_local_array(NULL), d_boost_vec_array(NULL),
      d_boost_local_vec_array(NULL), d_ghosted_local_vec(NULL), d_ghosted_local_array(NULL),
      d_boost_ghosted_local_array(NULL), d_boost_vec_ghosted_local_array(NULL)
{
    int ierr;
    int depth;
    ierr = VecGetBlockSize(d_global_vec, &depth);
    IBTK_CHKERRQ(ierr);
#if !defined(NDEBUG)
    TBOX_ASSERT(depth >= 0);
#endif
    d_depth = depth;
    int global_node_count;
    ierr = VecGetSize(d_global_vec, &global_node_count);
    IBTK_CHKERRQ(ierr);
#if !defined(NDEBUG)
    TBOX_ASSERT(global_node_count >= 0);
#endif
    d_global_node_count = global_node_count;
    d_global_node_count /= d_depth;
    int local_node_count;
    ierr = VecGetLocalSize(d_global_vec, &local_node_count);
    IBTK_CHKERRQ(ierr);
#if !defined(NDEBUG)
    TBOX_ASSERT(local_node_count >= 0);
#endif
    d_local_node_count = local_node_count;
    d_local_node_count /= d_depth;
    d_ghost_node_count = static_cast<int>(d_nonlocal_petsc_indices.size());
    return;
} // LData

PetscErrorCode loadVector(char * type_v,Vec * b){
	char file[PETSC_MAX_PATH_LEN];
	char err[PETSC_MAX_PATH_LEN];
	PetscErrorCode ierr;
	PetscBool flag;
	PetscViewer fd;
	PetscInt size;

	// check if there is a vec file, vector is not mandatory
	ierr=PetscOptionsGetString(PETSC_NULL,type_v,file,PETSC_MAX_PATH_LEN-1,&flag);CHKERRQ(ierr);
	if (!flag) {		
		PetscPrintf(PETSC_COMM_WORLD,"Error : %s is not properly set\n",type_v);
		*b = NULL;
	}else{
		PetscPrintf(PETSC_COMM_WORLD,"Loading Vector : %s\n",file);
		ierr=PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr);
		ierr=VecLoad(*b,fd);CHKERRQ(ierr);
		ierr=PetscViewerDestroy(&fd);CHKERRQ(ierr);
		ierr=VecGetSize(*b,&size);CHKERRQ(ierr);
		PetscPrintf(PETSC_COMM_WORLD,"Loaded Vector of size : %d\n",size);
	}

	return 0;
}

PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
{
  Vec               g = (Vec)ctx;
  const PetscScalar *xx,*gg;
  PetscScalar       *ff,d;
  PetscErrorCode    ierr;
  PetscInt          i,n;

  /*
     Get pointers 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 = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
  ierr = VecGetArray(f,&ff);CHKERRQ(ierr);
  ierr = VecGetArrayRead(g,&gg);CHKERRQ(ierr);

  /*
     Compute function
  */
  ierr  = VecGetSize(x,&n);CHKERRQ(ierr);
  d     = (PetscReal)(n - 1); d = d*d;
  ff[0] = xx[0];
  for (i=1; i<n-1; i++) ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - gg[i];
  ff[n-1] = xx[n-1] - 1.0;

  /*
     Restore vectors
  */
  ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
  ierr = VecRestoreArray(f,&ff);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(g,&gg);CHKERRQ(ierr);
  return 0;
}

PetscErrorCode gauss_seidel(PC pc,Vec bb,Vec xx,Vec w,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt m,PetscBool guesszero,PetscInt *its,PCRichardsonConvergedReason *reason)
{
  PetscInt          i,n1;
  PetscErrorCode    ierr;
  PetscScalar       *x;
  const PetscScalar *b;

  PetscFunctionBegin;
  *its    = m;
  *reason = PCRICHARDSON_CONVERGED_ITS;
  ierr    = VecGetSize(bb,&n1);CHKERRQ(ierr); n1--;
  ierr    = VecGetArrayRead(bb,&b);CHKERRQ(ierr);
  ierr    = VecGetArray(xx,&x);CHKERRQ(ierr);
  while (m--) {
    x[0] =  .5*(x[1] + b[0]);
    for (i=1; i<n1; i++) x[i] = .5*(x[i+1] + x[i-1] + b[i]);
    x[n1] = .5*(x[n1-1] + b[n1]);
    for (i=n1-1; i>0; i--) x[i] = .5*(x[i+1] + x[i-1] + b[i]);
    x[0] =  .5*(x[1] + b[0]);
  }
  ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr);
  ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@C
  TaoDefaultComputeGradient - computes the gradient using finite differences.

  Collective on Tao

  Input Parameters:
+ tao - the Tao context
. X - compute gradient at this point
- dummy - not used

  Output Parameters:
. G - Gradient Vector

   Options Database Key:
+  -tao_fd_gradient - Activates TaoDefaultComputeGradient()
-  -tao_fd_delta <delta> - change in x used to calculate finite differences

   Level: advanced

   Note:
   This routine is slow and expensive, and is not currently optimized
   to take advantage of sparsity in the problem.  Although
   TaoAppDefaultComputeGradient is not recommended for general use
   in large-scale applications, It can be useful in checking the
   correctness of a user-provided gradient.  Use the tao method TAOTEST
   to get an indication of whether your gradient is correct.


   Note:
   This finite difference gradient evaluation can be set using the routine TaoSetGradientRoutine() or by using the command line option -tao_fd_gradient

.seealso: TaoSetGradientRoutine()

@*/
PetscErrorCode TaoDefaultComputeGradient(Tao tao,Vec X,Vec G,void *dummy)
{
  PetscScalar    *x,*g;
  PetscReal      f, f2;
  PetscErrorCode ierr;
  PetscInt       low,high,N,i;
  PetscBool      flg;
  PetscReal      h=PETSC_SQRT_MACHINE_EPSILON;

  PetscFunctionBegin;
  ierr = TaoComputeObjective(tao, X,&f);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(NULL,"-tao_fd_delta",&h,&flg);CHKERRQ(ierr);
  ierr = VecGetSize(X,&N);CHKERRQ(ierr);
  ierr = VecGetOwnershipRange(X,&low,&high);CHKERRQ(ierr);
  ierr = VecGetArray(G,&g);CHKERRQ(ierr);
  for (i=0;i<N;i++) {
    if (i>=low && i<high) {
      ierr = VecGetArray(X,&x);CHKERRQ(ierr);
      x[i-low] += h;
      ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
    }

    ierr = TaoComputeObjective(tao,X,&f2);CHKERRQ(ierr);

    if (i>=low && i<high) {
      ierr = VecGetArray(X,&x);CHKERRQ(ierr);
      x[i-low] -= h;
      ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
    }
    if (i>=low && i<high) {
      g[i-low]=(f2-f)/h;
    }
  }
  ierr = VecRestoreArray(G,&g);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *dummy)
{
  const PetscScalar *xx;
  PetscScalar       *ff,*FF,d,d2;  
  PetscErrorCode    ierr;
  PetscInt          i,n;

  ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
  ierr = VecGetArray(f,&ff);CHKERRQ(ierr);
  ierr = VecGetArray((Vec)dummy,&FF);CHKERRQ(ierr);
  ierr = VecGetSize(x,&n);CHKERRQ(ierr);
  d    = (PetscReal)(n - 1); d2 = d*d;

  if (second_order) ff[0] = d*(0.5*d*(-xx[2] + 4.*xx[1] - 3.*xx[0]) - X0DOT);
  else ff[0] = d*(d*(xx[1] - xx[0]) - X0DOT);

  for (i=1; i<n-1; i++) ff[i] = d2*(xx[i-1] - 2.*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];

  ff[n-1] = d*d*(xx[n-1] - X1);
  ierr    = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
  ierr    = VecRestoreArray(f,&ff);CHKERRQ(ierr);
  ierr    = VecRestoreArray((Vec)dummy,&FF);CHKERRQ(ierr);
  return 0;
}

PetscErrorCode AddMuAbsorption(double *muinv, Vec muinvpml, double Qabs, int add)
{
  //compute muinvpml/(1+i/Qabs)
  double Qinv = (add==0) ? 0.0: (1.0/Qabs);
  double d=1 + pow(Qinv,2);
  PetscErrorCode ierr;
  int N;
  ierr=VecGetSize(muinvpml,&N);CHKERRQ(ierr);

  double *ptmuinvpml;
  ierr=VecGetArray(muinvpml, &ptmuinvpml);CHKERRQ(ierr);

  int i;
  double a,b;
  for(i=0;i<N/2;i++)
    {
      a=ptmuinvpml[i];
      b=ptmuinvpml[i+N/2];      
      muinv[i]= (a+b*Qinv)/d;
      muinv[i+N/2]=(b-a*Qinv)/d;
    }
  ierr=VecRestoreArray(muinvpml,&ptmuinvpml);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscInt main(PetscInt argc,char **args)
{
  PetscErrorCode  ierr;
  PetscMPIInt     rank,size;
  PetscInt        N0=3,N1=3,N2=3,N=N0*N1*N2;
  PetscRandom     rdm;
  PetscScalar     a;
  PetscReal       enorm;
  Vec             x,y,z,input,output;
  PetscBool       view=PETSC_FALSE,use_interface=PETSC_TRUE;
  Mat             A;
  PetscInt        DIM, dim[3],vsize;
  PetscReal       fac;

  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);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);


  ierr = PetscRandomCreate(PETSC_COMM_WORLD, &rdm);CHKERRQ(ierr);
  ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr);

  ierr = VecCreate(PETSC_COMM_WORLD,&input);CHKERRQ(ierr);
  ierr = VecSetSizes(input,PETSC_DECIDE,N);CHKERRQ(ierr);
  ierr = VecSetFromOptions(input);CHKERRQ(ierr);
  ierr = VecSetRandom(input,rdm);CHKERRQ(ierr);
  ierr = VecDuplicate(input,&output);
//  ierr = VecGetSize(input,&vsize);CHKERRQ(ierr);
//  printf("Size of the input Vector is %d\n",vsize);

  DIM = 3;
  dim[0] = N0; dim[1] = N1; dim[2] = N2;

  ierr = MatCreateFFT(PETSC_COMM_WORLD,DIM,dim,MATFFTW,&A);CHKERRQ(ierr);
  ierr = MatGetVecs(A,&x,&y);CHKERRQ(ierr);
  ierr = MatGetVecs(A,&z,PETSC_NULL);CHKERRQ(ierr);
  ierr = VecGetSize(y,&vsize);CHKERRQ(ierr);
  printf("The vector size from the main routine is %d\n",vsize);

  ierr = InputTransformFFT(A,input,x);CHKERRQ(ierr);
  ierr = MatMult(A,x,y);CHKERRQ(ierr);
  ierr = MatMultTranspose(A,y,z);CHKERRQ(ierr);
  ierr = OutputTransformFFT(A,z,output);CHKERRQ(ierr);

  fac = 1.0/(PetscReal)N;
  ierr = VecScale(output,fac);CHKERRQ(ierr);

  ierr = VecAssemblyBegin(input);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(input);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(output);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(output);CHKERRQ(ierr);

  ierr = VecView(input,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = VecView(output,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  ierr = VecAXPY(output,-1.0,input);CHKERRQ(ierr);
  ierr = VecNorm(output,NORM_1,&enorm);CHKERRQ(ierr);
//  if (enorm > 1.e-14){
    if (!rank)
      ierr = PetscPrintf(PETSC_COMM_SELF,"  Error norm of |x - z| %e\n",enorm);CHKERRQ(ierr);
//      }




// ierr = MatGetVecs(A,&z,PETSC_NULL);CHKERRQ(ierr);
//  printf("Vector size from ex148 %d\n",vsize);
//  ierr = PetscObjectSetName((PetscObject) x, "Real space vector");CHKERRQ(ierr);
//      ierr = PetscObjectSetName((PetscObject) y, "Frequency space vector");CHKERRQ(ierr);
//      ierr = PetscObjectSetName((PetscObject) z, "Reconstructed vector");CHKERRQ(ierr);

  ierr = PetscFinalize();
  return 0;

}

//.........这里部分代码省略.........
        PetscViewerASCIIPrintf( v, "\n");
    }
    if (h) {
        BSSCR_VecInfoLog(v,h,"stokes_b2");
        PetscViewerASCIIPrintf( v, "\n");
    }

    PetscViewerASCIIPopTab( v );


    /*--------------------------------------------------------------------------------------------*/
    PetscViewerASCIIPrintf( v, "--------------------------------------------------\n");
    PetscViewerASCIIPrintf( v, "Solution summary:\n");
    PetscViewerASCIIPushTab( v );

    if (u) {
        BSSCR_VecInfoLog(v,u,"x1");
        PetscViewerASCIIPrintf( v, "\n");
    }
    if (p) {
        BSSCR_VecInfoLog(v,p,"x2");
        PetscViewerASCIIPrintf( v, "\n");
    }

    {
        PetscScalar s,sum;
        PetscReal r,max,min;
        PetscInt N, loc;
        double r1,r2;
        Vec K_d;
        PetscInt loc_max, loc_min;


        VecGetSize( u, &N );
        VecMax( u, &loc, &r );
        PetscViewerASCIIPrintf( v, "u_max: %1.12e [%d] \n", r, loc );
        VecMin( u, &loc, &r );
        PetscViewerASCIIPrintf( v, "u_min: %1.12e [%d] \n", r, loc );

        VecDot( u,u, &s );
        PetscViewerASCIIPrintf( v, "u_rms: %1.12e \n", sqrt( PetscRealPart(s) )/N );

        VecDuplicate( u, &K_d );
        MatGetDiagonal( K, K_d );
        VecMax( K_d, &loc_max, &max );
        VecMin( K_d, &loc_min, &min );
        PetscViewerASCIIPrintf( v,"Algebraic contrast: max(K_d)=%.3e [%d] , min(K_d)=%.3e [%d] , max(K_d)/min(K_d) = %.8e \n", max,loc_max, min,loc_min, max/min );

        MatGetRowMax( K, K_d, PETSC_NULL );
        VecMax( K_d, &loc_max, &max );
        MatGetRowMin( K, K_d, PETSC_NULL );
        VecAbs( K_d );
        VecMin( K_d, &loc_min, &min );
        PetscViewerASCIIPrintf( v,"Algebraic contrast:   max(K)=%.3e [%d] , |min(K)|=%.3e [%d]  ,   max(K)/|min(K)| = %.8e \n", max,loc_max, min,loc_min, max/min );
        Stg_VecDestroy(&K_d );

        PetscViewerASCIIPrintf( v, "\n");

        VecGetSize( p, &N );
        VecMax( p, &loc, &r );
        PetscViewerASCIIPrintf( v, "p_max:  %1.12e [%d] \n", r, loc );
        VecMin( p, &loc, &r );
        PetscViewerASCIIPrintf( v, "p_min:  %1.12e [%d] \n", r, loc );

        VecDot( p,p, &s );
        PetscViewerASCIIPrintf( v, "p_rms:  %1.12e \n", sqrt( PetscRealPart(s) )/N );

int main(int argc, char **argv)
{
  MPI_Comm          comm;
  PetscMPIInt       rank;
  PetscErrorCode    ierr;
  User              user;
  PetscLogDouble       v1, v2;
  PetscInt          nplot = 0;
  char              filename1[2048], fileName[2048];
  PetscBool         set = PETSC_FALSE;
  PetscInt          steps_output;

  ierr = PetscInitialize(&argc, &argv, (char*) 0, help);CHKERRQ(ierr);
  comm = PETSC_COMM_WORLD;
  ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
  ierr = PetscNew(&user);CHKERRQ(ierr);
  ierr = PetscNew(&user->algebra);CHKERRQ(ierr);
  ierr = PetscNew(&user->model);CHKERRQ(ierr);
  ierr = PetscNew(&user->model->physics);CHKERRQ(ierr);

  Algebra   algebra = user->algebra;

  ierr = LoadOptions(comm, user);CHKERRQ(ierr);
  ierr = PetscTime(&v1);CHKERRQ(ierr);
  ierr = CreateMesh(comm, user);CHKERRQ(ierr);
  ierr = PetscTime(&v2);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
		       "Read and Distribute mesh takes %f sec \n", v2 - v1);CHKERRQ(ierr);
  ierr = SetUpLocalSpace(user);CHKERRQ(ierr); //Set up the dofs of each element
  ierr = ConstructGeometryFVM(&user->facegeom, &user->cellgeom, user);CHKERRQ(ierr);

  ierr = LimiterSetup(user);CHKERRQ(ierr);

  if(user->output_solution){
  // the output file options
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,0,"Options for output solution",0);CHKERRQ(ierr);
    ierr = PetscOptionsString("-solutionfile", "solution file", "AeroSim.c", filename1,filename1, 2048, &set);CHKERRQ(ierr);
    if(!set){SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,"please use option -solutionfile to specify solution file name \n");}
    ierr = PetscOptionsInt("-steps_output", "the number of time steps between two outputs", "", steps_output, &steps_output, &set);CHKERRQ(ierr);
    if(!set){ steps_output = 1;}
    ierr = PetscOptionsEnd();CHKERRQ(ierr);
  }

  if (user->TimeIntegralMethod == EXPLICITMETHOD) {
    if(user->myownexplicitmethod){
      ierr = PetscPrintf(PETSC_COMM_WORLD,"Using the fully explicit method based on my own routing\n");CHKERRQ(ierr);
      user->current_time = user->initial_time;
      user->current_step = 1;
      ierr = DMCreateGlobalVector(user->dm, &algebra->solution);CHKERRQ(ierr);
      ierr = PetscObjectSetName((PetscObject) algebra->solution, "solution");CHKERRQ(ierr);
      ierr = SetInitialCondition(user->dm, algebra->solution, user);CHKERRQ(ierr);
      ierr = VecDuplicate(algebra->solution, &algebra->fn);CHKERRQ(ierr);
      ierr = VecDuplicate(algebra->solution, &algebra->oldsolution);CHKERRQ(ierr);
      if(user->Explicit_RK2){
        ierr = PetscPrintf(PETSC_COMM_WORLD,"Use the second order Runge Kutta method \n");CHKERRQ(ierr);
      }else{
        ierr = PetscPrintf(PETSC_COMM_WORLD,"Use the first order forward Euler method \n");CHKERRQ(ierr);
      }
      nplot = 0; //the plot step
      while(user->current_time < (user->final_time - 0.05 * user->dt)){
        user->current_time = user->current_time + user->dt;

        ierr = FormTimeStepFunction(user, algebra, algebra->solution, algebra->fn);CHKERRQ(ierr);
        PetscReal fnnorm;
        ierr = VecNorm(algebra->fn,NORM_INFINITY,&fnnorm);CHKERRQ(ierr);
        if(0){
          PetscViewer    viewer;
          ierr = OutputVTK(user->dm, "function.vtk", &viewer);CHKERRQ(ierr);
          ierr = VecView(algebra->fn, viewer);CHKERRQ(ierr);
          ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
          ierr = PetscPrintf(PETSC_COMM_WORLD,"Step %D at time %g with founction norm = %g \n",
                                user->current_step, user->current_time, fnnorm);CHKERRQ(ierr);
          //break;
        }
        if(user->Explicit_RK2){
          ierr = VecCopy(algebra->solution, algebra->oldsolution);CHKERRQ(ierr);//U^n
          ierr = VecAXPY(algebra->solution, user->dt, algebra->fn);CHKERRQ(ierr);//U^{(1)}
          ierr = FormTimeStepFunction(user, algebra, algebra->solution, algebra->fn);CHKERRQ(ierr);//f(U^{(1)})
          ierr = VecAXPY(algebra->solution, 1.0, algebra->oldsolution);CHKERRQ(ierr);//U^n + U^{(1)}
          ierr = VecAXPY(algebra->solution, user->dt, algebra->fn);CHKERRQ(ierr);// + dt*f(U^{(1)})
          ierr = VecScale(algebra->solution, 0.5);CHKERRQ(ierr);
        }else{
          ierr = VecCopy(algebra->solution, algebra->oldsolution);CHKERRQ(ierr);
          ierr = VecAXPY(algebra->solution, user->dt, algebra->fn);CHKERRQ(ierr);
        }

        {// Monitor the solution and function norms
          PetscReal         norm;
          PetscLogDouble    space =0;
          PetscInt          size;

          ierr = VecNorm(algebra->solution,NORM_INFINITY,&norm);CHKERRQ(ierr);
          ierr = VecGetSize(algebra->solution, &size);CHKERRQ(ierr);
          norm = norm/size;
          if (norm>1.e5) {
            SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_LIB,
            "The norm of the solution is: %f (current time: %f). The explicit method is going to DIVERGE!!!", norm, user->current_time);
          }
          if (user->current_step%10==0) {
            ierr = PetscPrintf(PETSC_COMM_WORLD,"Step %D at time %g with solution norm = %g and founction norm = %g \n",
//.........这里部分代码省略.........

static PetscErrorCode KSPSolve_LSQR(KSP ksp)
{
  PetscErrorCode ierr;
  PetscInt       i,size1,size2;
  PetscScalar    rho,rhobar,phi,phibar,theta,c,s,tmp,tau;
  PetscReal      beta,alpha,rnorm;
  Vec            X,B,V,V1,U,U1,TMP,W,W2,SE,Z = NULL;
  Mat            Amat,Pmat;
  MatStructure   pflag;
  KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;
  PetscBool      diagonalscale,nopreconditioner;

  PetscFunctionBegin;
  ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr);
  if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);

  ierr = PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr);
  ierr = PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);CHKERRQ(ierr);

  /*  nopreconditioner =PETSC_FALSE; */
  /* Calculate norm of right hand side */
  ierr = VecNorm(ksp->vec_rhs,NORM_2,&lsqr->rhs_norm);CHKERRQ(ierr);

  /* mark norm of matrix with negative number to indicate it has not yet been computed */
  lsqr->anorm = -1.0;

  /* vectors of length m, where system size is mxn */
  B  = ksp->vec_rhs;
  U  = lsqr->vwork_m[0];
  U1 = lsqr->vwork_m[1];

  /* vectors of length n */
  X  = ksp->vec_sol;
  W  = lsqr->vwork_n[0];
  V  = lsqr->vwork_n[1];
  V1 = lsqr->vwork_n[2];
  W2 = lsqr->vwork_n[3];
  if (!nopreconditioner) Z = lsqr->vwork_n[4];

  /* standard error vector */
  SE = lsqr->se;
  if (SE) {
    ierr = VecGetSize(SE,&size1);CHKERRQ(ierr);
    ierr = VecGetSize(X,&size2);CHKERRQ(ierr);
    if (size1 != size2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Standard error vector (size %d) does not match solution vector (size %d)",size1,size2);
    ierr = VecSet(SE,0.0);CHKERRQ(ierr);
  }

  /* Compute initial residual, temporarily use work vector u */
  if (!ksp->guess_zero) {
    ierr = KSP_MatMult(ksp,Amat,X,U);CHKERRQ(ierr);       /*   u <- b - Ax     */
    ierr = VecAYPX(U,-1.0,B);CHKERRQ(ierr);
  } else {
    ierr = VecCopy(B,U);CHKERRQ(ierr);            /*   u <- b (x is 0) */
  }

  /* Test for nothing to do */
  ierr       = VecNorm(U,NORM_2,&rnorm);CHKERRQ(ierr);
  ierr       = PetscObjectAMSTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
  ksp->its   = 0;
  ksp->rnorm = rnorm;
  ierr       = PetscObjectAMSGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
  ierr = KSPLogResidualHistory(ksp,rnorm);CHKERRQ(ierr);
  ierr = KSPMonitor(ksp,0,rnorm);CHKERRQ(ierr);
  ierr = (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  if (ksp->reason) PetscFunctionReturn(0);

  beta = rnorm;
  ierr = VecScale(U,1.0/beta);CHKERRQ(ierr);
  ierr = KSP_MatMultTranspose(ksp,Amat,U,V);CHKERRQ(ierr);
  if (nopreconditioner) {
    ierr = VecNorm(V,NORM_2,&alpha);CHKERRQ(ierr);
  } else {
    ierr = PCApply(ksp->pc,V,Z);CHKERRQ(ierr);
    ierr = VecDotRealPart(V,Z,&alpha);CHKERRQ(ierr);
    if (alpha <= 0.0) {
      ksp->reason = KSP_DIVERGED_BREAKDOWN;
      PetscFunctionReturn(0);
    }
    alpha = PetscSqrtReal(alpha);
    ierr  = VecScale(Z,1.0/alpha);CHKERRQ(ierr);
  }
  ierr = VecScale(V,1.0/alpha);CHKERRQ(ierr);

  if (nopreconditioner) {
    ierr = VecCopy(V,W);CHKERRQ(ierr);
  } else {
    ierr = VecCopy(Z,W);CHKERRQ(ierr);
  }

  lsqr->arnorm = alpha * beta;
  phibar       = beta;
  rhobar       = alpha;
  i            = 0;
  do {
    if (nopreconditioner) {
      ierr = KSP_MatMult(ksp,Amat,V,U1);CHKERRQ(ierr);
    } else {
      ierr = KSP_MatMult(ksp,Amat,Z,U1);CHKERRQ(ierr);
    }
//.........这里部分代码省略.........

PetscErrorCode VecLoad_Binary(Vec vec, PetscViewer viewer)
{
  PetscMPIInt    size,rank,tag;
  int            fd;
  PetscInt       i,rows = 0,n,*range,N,bs;
  PetscErrorCode ierr;
  PetscBool      flag;
  PetscScalar    *avec,*avecwork;
  MPI_Comm       comm;
  MPI_Request    request;
  MPI_Status     status;
#if defined(PETSC_HAVE_MPIIO)
  PetscBool      useMPIIO;
#endif

  PetscFunctionBegin;
  ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);

  ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr);
  ierr = PetscViewerBinaryReadVecHeader_Private(viewer,&rows);CHKERRQ(ierr);
  /* Set Vec sizes,blocksize,and type if not already set. Block size first so that local sizes will be compatible. */
  ierr = PetscOptionsGetInt(((PetscObject)vec)->prefix, "-vecload_block_size", &bs, &flag);CHKERRQ(ierr);
  if (flag) {
    ierr = VecSetBlockSize(vec, bs);CHKERRQ(ierr);
  }
  if (vec->map->n < 0 && vec->map->N < 0) {
    ierr = VecSetSizes(vec,PETSC_DECIDE,rows);CHKERRQ(ierr);
  }

  /* If sizes and type already set,check if the vector global size is correct */
  ierr = VecGetSize(vec, &N);CHKERRQ(ierr);
  if (N != rows) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%d) then input vector (%d)", rows, N);

#if defined(PETSC_HAVE_MPIIO)
  ierr = PetscViewerBinaryGetMPIIO(viewer,&useMPIIO);CHKERRQ(ierr);
  if (useMPIIO) {
    ierr = VecLoad_Binary_MPIIO(vec, viewer);CHKERRQ(ierr);
    PetscFunctionReturn(0);
  }
#endif

  ierr = VecGetLocalSize(vec,&n);CHKERRQ(ierr);
  ierr = PetscObjectGetNewTag((PetscObject)viewer,&tag);CHKERRQ(ierr);
  ierr = VecGetArray(vec,&avec);CHKERRQ(ierr);
  if (!rank) {
    ierr = PetscBinaryRead(fd,avec,n,PETSC_SCALAR);CHKERRQ(ierr);

    if (size > 1) {
      /* read in other chuncks and send to other processors */
      /* determine maximum chunck owned by other */
      range = vec->map->range;
      n = 1;
      for (i=1; i<size; i++) n = PetscMax(n,range[i+1] - range[i]);

      ierr = PetscMalloc(n*sizeof(PetscScalar),&avecwork);CHKERRQ(ierr);
      for (i=1; i<size; i++) {
        n    = range[i+1] - range[i];
        ierr = PetscBinaryRead(fd,avecwork,n,PETSC_SCALAR);CHKERRQ(ierr);
        ierr = MPI_Isend(avecwork,n,MPIU_SCALAR,i,tag,comm,&request);CHKERRQ(ierr);
        ierr = MPI_Wait(&request,&status);CHKERRQ(ierr);
      }
      ierr = PetscFree(avecwork);CHKERRQ(ierr);
    }
  } else {
    ierr = MPI_Recv(avec,n,MPIU_SCALAR,0,tag,comm,&status);CHKERRQ(ierr);
  }

  ierr = VecRestoreArray(vec,&avec);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(vec);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(vec);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}