C++ KSPSetOperators函数代码示例

PetscErrorCode StokesSetupPC(Stokes *s, KSP ksp)
{
  KSP            *subksp;
  PC             pc;
  PetscInt       n = 1;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
  ierr = PCFieldSplitSetIS(pc, "0", s->isg[0]);CHKERRQ(ierr);
  ierr = PCFieldSplitSetIS(pc, "1", s->isg[1]);CHKERRQ(ierr);
  if (s->userPC) {
    ierr = PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, s->myS);CHKERRQ(ierr);
  }
  if (s->userKSP) {
    ierr = PCSetUp(pc);CHKERRQ(ierr);
    ierr = PCFieldSplitGetSubKSP(pc, &n, &subksp);CHKERRQ(ierr);
    ierr = KSPSetOperators(subksp[1], s->myS, s->myS);CHKERRQ(ierr);
    ierr = PetscFree(subksp);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

static PetscErrorCode ComputeKSPFETIDP(DomainData dd, KSP ksp_bddc, KSP *ksp_fetidp)
{
  PetscErrorCode ierr;
  KSP            temp_ksp;
  PC             pc,D;
  Mat            F;

  PetscFunctionBeginUser;
  ierr        = KSPGetPC(ksp_bddc,&pc);CHKERRQ(ierr);
  ierr        = PCBDDCCreateFETIDPOperators(pc,&F,&D);CHKERRQ(ierr);
  ierr        = KSPCreate(PetscObjectComm((PetscObject)F),&temp_ksp);CHKERRQ(ierr);
  ierr        = KSPSetOperators(temp_ksp,F,F);CHKERRQ(ierr);
  ierr        = KSPSetType(temp_ksp,KSPCG);CHKERRQ(ierr);
  ierr        = KSPSetPC(temp_ksp,D);CHKERRQ(ierr);
  ierr        = KSPSetComputeSingularValues(temp_ksp,PETSC_TRUE);CHKERRQ(ierr);
  ierr        = KSPSetFromOptions(temp_ksp);CHKERRQ(ierr);
  ierr        = KSPSetUp(temp_ksp);CHKERRQ(ierr);
  *ksp_fetidp = temp_ksp;
  ierr        = MatDestroy(&F);CHKERRQ(ierr);
  ierr        = PCDestroy(&D);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

HeatSolverBTCS::HeatSolverBTCS(int ny_, double dy_, int nx_, double dx_):
	Solver(ny_, dy_, nx_, dx_)
{
	// Create a sparse matrix A
	two_d_heat_BTCS(A, ny, dy, nx, dx);
	dt = 1.0;
	// set up linear solver context (ksp) and preconditioner (pc)
	KSPCreate(PETSC_COMM_WORLD,&ksp);
	KSPSetOperators(ksp,A,A);
	KSPGetPC(ksp,&pc);
	PCSetType(pc,PCLU);
	KSPSetFromOptions(ksp);

	// create rhs
	VecCreateSeq(PETSC_COMM_SELF, nx*ny, &rhs);
	VecSetFromOptions(rhs);
	VecSet(rhs,0.0);

	// prepare temp
	VecCreateSeq(PETSC_COMM_SELF, nx*ny, &temp);
	VecSetFromOptions(temp);
}

int main(int argc, char **argv){
	PetscErrorCode ierr;
	Vec x,b;
	Mat A;
	KSP ksp;

	ierr=PetscInitialize(&argc,&argv,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
	PetscPrintf(PETSC_COMM_WORLD,"]> Initializing PETSc/SLEPc\n");
	
	/*Load data*/
	ierr=loadInputs(&A,&b,&x);CHKERRQ(ierr);
	PetscPrintf(PETSC_COMM_WORLD,"]> Data loaded\n");

	/*Create the KSP context and setup*/
	ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);	
	ierr = KSPSetType(ksp,KSPFGMRES);CHKERRQ(ierr);	
	ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);	
	ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);	
	ierr = KSPSetUp(ksp);CHKERRQ(ierr);	
	PetscPrintf(PETSC_COMM_WORLD,"]> Krylov Solver settings done\n");

	/*Solve the system*/
	PetscPrintf(PETSC_COMM_WORLD,"]> Krylov Solver Launching solving process\n");
	ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
	PetscPrintf(PETSC_COMM_WORLD,"]> Krylov Solver System solved\n");

	/*Clean*/
	ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
	ierr = VecDestroy(&b);CHKERRQ(ierr);
	ierr = VecDestroy(&x);CHKERRQ(ierr);
	ierr = MatDestroy(&A);CHKERRQ(ierr);
	PetscPrintf(PETSC_COMM_WORLD,"]> Cleaned structures, finalizing\n");

	/*Finalize PETSc*/
	PetscFinalize(); 

	return 0;
}

PetscErrorCode PressurePoissonCreate(  )
{
  PetscErrorCode ierr;
  
  
  PetscFunctionBegin;
  PetscLogEventBegin(EVENT_PressurePoissonCreate,0,0,0,0);
//  PetscLogEventRegister(&EVENT_PressurePoissonCreate,"PressurePoissonCreate", 0);
  ierr = KSPCreate(PETSC_COMM_SELF, &ksp); CHKERRQ(ierr);
  KSPGetPC(ksp, &pc);
//  KSPSetType(ksp, KSPCG);
//  PCSetType(pc, PCICC);
//  PCFactorSetLevels(pc, 0);
//  KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
  KSPSetType(ksp, KSPPREONLY);
  PCSetType(pc, PCCHOLESKY);
  PCFactorSetMatOrderingType(pc, MATORDERING_ND);
  KSPSetFromOptions(ksp);
  KSPSetOperators(ksp, m, m, SAME_PRECONDITIONER);
  PCSetUp(pc);
  
  ierr = IIMCreate(&iim, 12); CHKERRQ(ierr);
//  IIMSetForceComponents(iim, ForceComponentNormalSimple, ForceComponentTangentialSimple);
  CreateGrid2D(d1, d2, &rhsC);
  CreateGrid2D(d1, d2, &p);
  CreateGrid2D(d1, d2, &px);
  CreateGrid2D(d1, d2, &py);
  CreateGrid2D(d1, d2, &u);
  CreateGrid2D(d1, d2, &v);
  ierr = CreateLevelSet2D(d1,d2,&ls); CHKERRQ(ierr);  
  CreateLevelSet2D(d1, d2, &lstemp);
  LevelSetInitializeToStar(ls,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
//  LevelSetInitializeToCircle(ls, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);
  
  PetscLogEventEnd(EVENT_PressurePoissonCreate,0,0,0,0);
  PetscFunctionReturn(0);
}

void PETSc::Solve_withPureNeumann(void)
{
	ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
  	ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
  	
  	ierr = VecAssemblyBegin(x);
  	ierr = VecAssemblyEnd(x);
  	
  	ierr = VecAssemblyBegin(b);
  	ierr = VecAssemblyEnd(b);
  	
	
	MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,PETSC_NULL,&nullsp);
        KSPSetNullSpace(ksp,nullsp);
	
        KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
        
	//KSPSetType(ksp,KSPMINRES);
	//KSPSetType(ksp,KSPGMRES);
	//KSPSetType(ksp,KSPBCGS);
	KSPSetType(ksp,KSPBCGSL);
	KSPBCGSLSetEll(ksp,2);

	//KSPGetPC(ksp, &pc);
	//PCSetType(pc, PCASM);
	//PCSetType(pc, PCMG);
	//PCMGSetLevels(pc, 3, &PETSC_COMM_WORLD);
	//PCMGSetType(pc,PC_MG_MULTIPLICATIVE);
	//PCMGSetCycleType(pc,PC_MG_CYCLE_V);
	//
        KSPSetFromOptions(ksp);
        KSPSetUp(ksp);

	//start_clock("Before Petsc Solve in pure neumann solver");
        KSPSolve(ksp,b,x);
	//stop_clock("After Petsc Solve in pure neumann solver");
}

int main(int argc, char** argv)
{
  Mat            Q;
  Vec            v, a, se;
  KSP            QRsolver;
  PC             pc;
  PetscErrorCode ierr;

  ierr = PetscInitialize(&argc, &argv, NULL, NULL);if (ierr) return ierr;

  ierr = VecCreate(PETSC_COMM_WORLD, &v);CHKERRQ(ierr);
  ierr = MatCreate(PETSC_COMM_WORLD, &Q);CHKERRQ(ierr);
  ierr = MatSetType(Q, MATDENSE);CHKERRQ(ierr);
  ierr = fill(Q, v);CHKERRQ(ierr);

  ierr = MatCreateVecs(Q, &a, NULL);CHKERRQ(ierr);
  ierr = KSPCreate(PETSC_COMM_WORLD, &QRsolver);CHKERRQ(ierr);
  ierr = KSPGetPC(QRsolver, &pc);CHKERRQ(ierr);
  ierr = PCSetType(pc, PCNONE);CHKERRQ(ierr);
  ierr = KSPSetType(QRsolver, KSPLSQR);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(QRsolver);CHKERRQ(ierr);
  ierr = KSPSetOperators(QRsolver, Q, Q);CHKERRQ(ierr);
  ierr = MatViewFromOptions(Q, NULL, "-sys_view");CHKERRQ(ierr);
  ierr = VecViewFromOptions(a, NULL, "-rhs_view");CHKERRQ(ierr);
  ierr = KSPSolve(QRsolver, v, a);CHKERRQ(ierr);
  ierr = KSPLSQRGetStandardErrorVec(QRsolver, &se);CHKERRQ(ierr);
  if (se) {
    ierr = VecViewFromOptions(se, NULL, "-se_view");CHKERRQ(ierr);
  }
  ierr = KSPDestroy(&QRsolver);CHKERRQ(ierr);
  ierr = VecDestroy(&a);CHKERRQ(ierr);
  ierr = VecDestroy(&v);CHKERRQ(ierr);
  ierr = MatDestroy(&Q);CHKERRQ(ierr);

  ierr = PetscFinalize();
  return ierr;
}

PetscErrorCode Petsc_KSPSolve(AppCtx *obj)
{
  PetscErrorCode ierr;
  KSP            ksp;
  PC             pc;

  /*create the ksp context and set the operators,that is, associate the system matrix with it*/
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp,obj->Amat,obj->Amat);CHKERRQ(ierr);

  /*get the preconditioner context, set its type and the tolerances*/
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
  ierr = PCSetType(pc,PCLU);CHKERRQ(ierr);
  ierr = KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);

  /*get the command line options if there are any and set them*/
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

  /*get the linear system (ksp) solve*/
  ierr = KSPSolve(ksp,obj->ksp_rhs,obj->ksp_sol);CHKERRQ(ierr);

  KSPDestroy(&ksp);
  return 0;
}

void Field_solver::create_solver_and_preconditioner( KSP *ksp, PC *pc, Mat *A )
{
    PetscReal rtol = 1.e-12;
    // Default.     
    // Possible to specify from command line using '-ksp_rtol' option.
    
    PetscErrorCode ierr;
    ierr = KSPCreate( PETSC_COMM_WORLD, ksp ); CHKERRXX(ierr);
    ierr = KSPSetOperators( *ksp, *A, *A, DIFFERENT_NONZERO_PATTERN ); CHKERRXX(ierr);
    //ierr = KSPSetOperators( *ksp, *A, *A ); CHKERRXX(ierr);
    ierr = KSPGetPC( *ksp, pc ); CHKERRXX(ierr);
    ierr = PCSetType( *pc, PCGAMG ); CHKERRXX(ierr);
    ierr = KSPSetType( *ksp, KSPGMRES ); CHKERRXX(ierr);
    ierr = KSPSetTolerances( *ksp, rtol,
			     PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT); CHKERRXX(ierr);
    ierr = KSPSetFromOptions( *ksp ); CHKERRXX(ierr);    
    ierr = KSPSetInitialGuessNonzero( *ksp, PETSC_TRUE ); CHKERRXX( ierr );

    // For test purposes
    //ierr = KSPSetInitialGuessNonzero( *ksp, PETSC_FALSE ); CHKERRXX( ierr );

    ierr = KSPSetUp( *ksp ); CHKERRXX(ierr);
    return;
}

int main(int argc,char **args)
{
  Vec            x,b,u;      /* approx solution, RHS, exact solution */
  Mat            A;            /* linear system matrix */
  KSP            ksp;         /* KSP context */
  KSP            *subksp;     /* array of local KSP contexts on this processor */
  PC             pc;           /* PC context */
  PC             subpc;        /* PC context for subdomain */
  PetscReal      norm;         /* norm of solution error */
  PetscErrorCode ierr;
  PetscInt       i,j,Ii,J,*blks,m = 4,n;
  PetscMPIInt    rank,size;
  PetscInt       its,nlocal,first,Istart,Iend;
  PetscScalar    v,one = 1.0,none = -1.0;
  PetscBool      isbjacobi;

  PetscInitialize(&argc,&args,(char*)0,help);
  ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  n    = m+2;

  /* -------------------------------------------------------------------
         Compute the matrix and right-hand-side vector that define
         the linear system, Ax = b.
     ------------------------------------------------------------------- */

  /*
     Create and assemble parallel matrix
  */
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);CHKERRQ(ierr);
  ierr = MatSeqAIJSetPreallocation(A,5,NULL);CHKERRQ(ierr);
  ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
  for (Ii=Istart; Ii<Iend; Ii++) {
    v = -1.0; i = Ii/n; j = Ii - i*n;
    if (i>0)   {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
    if (i<m-1) {J = Ii + n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
    if (j>0)   {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
    if (j<n-1) {J = Ii + 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
    v = 4.0; ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);CHKERRQ(ierr);
  }
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  /*
     Create parallel vectors
  */
  ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
  ierr = VecSetSizes(u,PETSC_DECIDE,m*n);CHKERRQ(ierr);
  ierr = VecSetFromOptions(u);CHKERRQ(ierr);
  ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
  ierr = VecDuplicate(b,&x);CHKERRQ(ierr);

  /*
     Set exact solution; then compute right-hand-side vector.
  */
  ierr = VecSet(u,one);CHKERRQ(ierr);
  ierr = MatMult(A,u,b);CHKERRQ(ierr);

  /*
     Create linear solver context
  */
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);

  /*
     Set operators. Here the matrix that defines the linear system
     also serves as the preconditioning matrix.
  */
  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);

  /*
     Set default preconditioner for this program to be block Jacobi.
     This choice can be overridden at runtime with the option
        -pc_type <type>
  */
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
  ierr = PCSetType(pc,PCBJACOBI);CHKERRQ(ierr);


  /* -------------------------------------------------------------------
                   Define the problem decomposition
     ------------------------------------------------------------------- */

  /*
     Call PCBJacobiSetTotalBlocks() to set individually the size of
     each block in the preconditioner.  This could also be done with
     the runtime option
         -pc_bjacobi_blocks <blocks>
     Also, see the command PCBJacobiSetLocalBlocks() to set the
     local blocks.

      Note: The default decomposition is 1 block per processor.
  */
  ierr = PetscMalloc1(m,&blks);CHKERRQ(ierr);
  for (i=0; i<m; i++) blks[i] = n;
  ierr = PCBJacobiSetTotalBlocks(pc,m,blks);CHKERRQ(ierr);
  ierr = PetscFree(blks);CHKERRQ(ierr);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  PetscReal      hx, hy, area;
  const PetscInt *necon;
  PetscInt       j, k, e, ne, nc, mx, my;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = PetscObjectGetComm((PetscObject) da, &comm);CHKERRQ(ierr);
  ierr = DMSetMatType(da,MATAIJ);CHKERRQ(ierr);
  ierr = DMCreateMatrix(da, &mat);CHKERRQ(ierr);
  ierr = MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
  ierr = DMGetGlobalVector(da, &rhs_u);CHKERRQ(ierr);
  ierr = DMGetGlobalVector(da, &rhs_v);CHKERRQ(ierr);
  ierr = KSPCreate(comm, &ksp);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

  ierr = DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
  hx   = 1.0 / (PetscReal)(mx-1);
  hy   = 1.0 / (PetscReal)(my-1);
  area = 0.5*hx*hy;
  ierr = VecGetArray(user->sol_n.u,       &u_n);CHKERRQ(ierr);
  ierr = VecGetArray(user->sol_n.v,       &v_n);CHKERRQ(ierr);
  ierr = VecGetArray(user->mu,            &mu_n);CHKERRQ(ierr);
  ierr = VecGetArray(user->sol_phi.u,     &u_phi);CHKERRQ(ierr);
  ierr = VecGetArray(user->sol_phi.v,     &v_phi);CHKERRQ(ierr);
  ierr = VecGetArray(user->sol_phi.rho_u, &rho_u_phi);CHKERRQ(ierr);
  ierr = VecGetArray(user->sol_phi.rho_v, &rho_v_phi);CHKERRQ(ierr);
  ierr = DMDAGetElements(da, &ne, &nc, &necon);CHKERRQ(ierr);
  for (e = 0; e < ne; e++) {
    for (j = 0; j < 3; j++) {
      idx[j]      = necon[3*e+j];
      values_u[j] = 0.0;
      values_v[j] = 0.0;
    }
    /* Get basis function deriatives (we need the orientation of the element here) */
    if (idx[1] > idx[0]) {
      psi_x[0] = -hy; psi_x[1] =  hy; psi_x[2] = 0.0;
      psi_y[0] = -hx; psi_y[1] = 0.0; psi_y[2] =  hx;
    } else {
      psi_x[0] =  hy; psi_x[1] = -hy; psi_x[2] = 0.0;
      psi_y[0] =  hx; psi_y[1] = 0.0; psi_y[2] = -hx;
    }
    /*  <\nabla\psi, F^{n+\phi}_e>: Divergence of the element-averaged convective fluxes */
    for (j = 0; j < 3; j++) {
      values_u[j] += psi_x[j]*rho_u_phi[e]*u_phi[e] + psi_y[j]*rho_u_phi[e]*v_phi[e];
      values_v[j] += psi_x[j]*rho_v_phi[e]*u_phi[e] + psi_y[j]*rho_v_phi[e]*v_phi[e];
    }
    /*  -<\nabla\psi, F^n_v>: Divergence of the viscous fluxes */
    for (j = 0; j < 3; j++) {
      /* \tau_{xx} = 2/3 \mu(T) (2 {\partial u\over\partial x} - {\partial v\over\partial y}) */
      /* \tau_{xy} =     \mu(T) (  {\partial u\over\partial y} + {\partial v\over\partial x}) */
      /* \tau_{yy} = 2/3 \mu(T) (2 {\partial v\over\partial y} - {\partial u\over\partial x}) */
      mu     = 0.0;
      tau_xx = 0.0;
      tau_xy = 0.0;
      tau_yy = 0.0;
      for (k = 0; k < 3; k++) {
        mu     += mu_n[idx[k]];
        tau_xx += 2.0*psi_x[k]*u_n[idx[k]] - psi_y[k]*v_n[idx[k]];
        tau_xy +=     psi_y[k]*u_n[idx[k]] + psi_x[k]*v_n[idx[k]];
        tau_yy += 2.0*psi_y[k]*v_n[idx[k]] - psi_x[k]*u_n[idx[k]];
      }
      mu          /= 3.0;
      tau_xx      *= (2.0/3.0)*mu;
      tau_xy      *= mu;
      tau_yy      *= (2.0/3.0)*mu;
      values_u[j] -= area*(psi_x[j]*tau_xx + psi_y[j]*tau_xy);
      values_v[j] -= area*(psi_x[j]*tau_xy + psi_y[j]*tau_yy);
    }
    /* Accumulate to global structures */
    ierr = VecSetValuesLocal(rhs_u, 3, idx, values_u, ADD_VALUES);CHKERRQ(ierr);
    ierr = VecSetValuesLocal(rhs_v, 3, idx, values_v, ADD_VALUES);CHKERRQ(ierr);
    ierr = MatSetValuesLocal(mat, 3, idx, 3, idx, identity, ADD_VALUES);CHKERRQ(ierr);
  }
  ierr = DMDARestoreElements(da, &ne, &nc, &necon);CHKERRQ(ierr);
  ierr = VecRestoreArray(user->sol_n.u,       &u_n);CHKERRQ(ierr);
  ierr = VecRestoreArray(user->sol_n.v,       &v_n);CHKERRQ(ierr);
  ierr = VecRestoreArray(user->mu,            &mu_n);CHKERRQ(ierr);
  ierr = VecRestoreArray(user->sol_phi.u,     &u_phi);CHKERRQ(ierr);
  ierr = VecRestoreArray(user->sol_phi.v,     &v_phi);CHKERRQ(ierr);
  ierr = VecRestoreArray(user->sol_phi.rho_u, &rho_u_phi);CHKERRQ(ierr);
  ierr = VecRestoreArray(user->sol_phi.rho_v, &rho_v_phi);CHKERRQ(ierr);

  ierr = VecAssemblyBegin(rhs_u);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(rhs_v);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(rhs_u);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(rhs_v);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = VecScale(rhs_u,user->dt);CHKERRQ(ierr);
  ierr = VecScale(rhs_v,user->dt);CHKERRQ(ierr);

  ierr = KSPSetOperators(ksp, mat, mat, DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPSolve(ksp, rhs_u, user->sol_np1.rho_u);CHKERRQ(ierr);
  ierr = KSPSolve(ksp, rhs_v, user->sol_np1.rho_v);CHKERRQ(ierr);
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  ierr = MatDestroy(&mat);CHKERRQ(ierr);
  ierr = DMRestoreGlobalVector(da, &rhs_u);CHKERRQ(ierr);
  ierr = DMRestoreGlobalVector(da, &rhs_v);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode PCSetUp_MG(PC pc)
{
  PC_MG          *mg        = (PC_MG*)pc->data;
  PC_MG_Levels   **mglevels = mg->levels;
  PetscErrorCode ierr;
  PetscInt       i,n = mglevels[0]->levels;
  PC             cpc;
  PetscBool      dump = PETSC_FALSE,opsset,use_amat,missinginterpolate = PETSC_FALSE;
  Mat            dA,dB;
  Vec            tvec;
  DM             *dms;
  PetscViewer    viewer = 0;

  PetscFunctionBegin;
  /* FIX: Move this to PCSetFromOptions_MG? */
  if (mg->usedmfornumberoflevels) {
    PetscInt levels;
    ierr = DMGetRefineLevel(pc->dm,&levels);CHKERRQ(ierr);
    levels++;
    if (levels > n) { /* the problem is now being solved on a finer grid */
      ierr     = PCMGSetLevels(pc,levels,NULL);CHKERRQ(ierr);
      n        = levels;
      ierr     = PCSetFromOptions(pc);CHKERRQ(ierr); /* it is bad to call this here, but otherwise will never be called for the new hierarchy */
      mglevels =  mg->levels;
    }
  }
  ierr = KSPGetPC(mglevels[0]->smoothd,&cpc);CHKERRQ(ierr);


  /* If user did not provide fine grid operators OR operator was not updated since last global KSPSetOperators() */
  /* so use those from global PC */
  /* Is this what we always want? What if user wants to keep old one? */
  ierr = KSPGetOperatorsSet(mglevels[n-1]->smoothd,NULL,&opsset);CHKERRQ(ierr);
  if (opsset) {
    Mat mmat;
    ierr = KSPGetOperators(mglevels[n-1]->smoothd,NULL,&mmat);CHKERRQ(ierr);
    if (mmat == pc->pmat) opsset = PETSC_FALSE;
  }

  if (!opsset) {
    ierr = PCGetUseAmat(pc,&use_amat);CHKERRQ(ierr);
    if(use_amat){
      ierr = PetscInfo(pc,"Using outer operators to define finest grid operator \n  because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");CHKERRQ(ierr);
      ierr = KSPSetOperators(mglevels[n-1]->smoothd,pc->mat,pc->pmat);CHKERRQ(ierr);
    }
    else {
      ierr = PetscInfo(pc,"Using matrix (pmat) operators to define finest grid operator \n  because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");CHKERRQ(ierr);
      ierr = KSPSetOperators(mglevels[n-1]->smoothd,pc->pmat,pc->pmat);CHKERRQ(ierr);
    }
  }

  for (i=n-1; i>0; i--) {
    if (!(mglevels[i]->interpolate || mglevels[i]->restrct)) {
      missinginterpolate = PETSC_TRUE;
      continue;
    }
  }
  /*
   Skipping if user has provided all interpolation/restriction needed (since DM might not be able to produce them (when coming from SNES/TS)
   Skipping for galerkin==2 (externally managed hierarchy such as ML and GAMG). Cleaner logic here would be great. Wrap ML/GAMG as DMs?
  */
  if (missinginterpolate && pc->dm && mg->galerkin != 2 && !pc->setupcalled) {
    /* construct the interpolation from the DMs */
    Mat p;
    Vec rscale;
    ierr     = PetscMalloc1(n,&dms);CHKERRQ(ierr);
    dms[n-1] = pc->dm;
    /* Separately create them so we do not get DMKSP interference between levels */
    for (i=n-2; i>-1; i--) {ierr = DMCoarsen(dms[i+1],MPI_COMM_NULL,&dms[i]);CHKERRQ(ierr);}
    for (i=n-2; i>-1; i--) {
      DMKSP     kdm;
      PetscBool dmhasrestrict;
      ierr = KSPSetDM(mglevels[i]->smoothd,dms[i]);CHKERRQ(ierr);
      if (mg->galerkin) {ierr = KSPSetDMActive(mglevels[i]->smoothd,PETSC_FALSE);CHKERRQ(ierr);}
      ierr = DMGetDMKSPWrite(dms[i],&kdm);CHKERRQ(ierr);
      /* Ugly hack so that the next KSPSetUp() will use the RHS that we set. A better fix is to change dmActive to take
       * a bitwise OR of computing the matrix, RHS, and initial iterate. */
      kdm->ops->computerhs = NULL;
      kdm->rhsctx          = NULL;
      if (!mglevels[i+1]->interpolate) {
        ierr = DMCreateInterpolation(dms[i],dms[i+1],&p,&rscale);CHKERRQ(ierr);
        ierr = PCMGSetInterpolation(pc,i+1,p);CHKERRQ(ierr);
        if (rscale) {ierr = PCMGSetRScale(pc,i+1,rscale);CHKERRQ(ierr);}
        ierr = VecDestroy(&rscale);CHKERRQ(ierr);
        ierr = MatDestroy(&p);CHKERRQ(ierr);
      }
      ierr = DMHasCreateRestriction(dms[i],&dmhasrestrict);CHKERRQ(ierr);
      if (dmhasrestrict && !mglevels[i+1]->restrct){
        ierr = DMCreateRestriction(dms[i],dms[i+1],&p);CHKERRQ(ierr);
        ierr = PCMGSetRestriction(pc,i+1,p);CHKERRQ(ierr);
        ierr = MatDestroy(&p);CHKERRQ(ierr);
      }
    }

    for (i=n-2; i>-1; i--) {ierr = DMDestroy(&dms[i]);CHKERRQ(ierr);}
    ierr = PetscFree(dms);CHKERRQ(ierr);
  }

  if (pc->dm && !pc->setupcalled) {
    /* finest smoother also gets DM but it is not active, independent of whether galerkin==2 */
//.........这里部分代码省略.........

int main(int argc,char **args)
{
    Mat            C;
    PetscErrorCode ierr;
    PetscInt       N = 2,rowidx,colidx;
    Vec            u,b,r;
    KSP            ksp;
    PetscReal      norm;
    PetscMPIInt    rank,size;
    PetscScalar    v;

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

    /* create stiffness matrix C = [1 2; 2 3] */
    ierr = MatCreate(PETSC_COMM_WORLD,&C);
    CHKERRQ(ierr);
    ierr = MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N);
    CHKERRQ(ierr);
    ierr = MatSetFromOptions(C);
    CHKERRQ(ierr);
    ierr = MatSetUp(C);
    CHKERRQ(ierr);
    if (!rank) {
        rowidx = 0;
        colidx = 0;
        v = 1.0;
        ierr   = MatSetValues(C,1,&rowidx,1,&colidx,&v,INSERT_VALUES);
        CHKERRQ(ierr);
        rowidx = 0;
        colidx = 1;
        v = 2.0;
        ierr   = MatSetValues(C,1,&rowidx,1,&colidx,&v,INSERT_VALUES);
        CHKERRQ(ierr);

        rowidx = 1;
        colidx = 0;
        v = 2.0;
        ierr   = MatSetValues(C,1,&rowidx,1,&colidx,&v,INSERT_VALUES);
        CHKERRQ(ierr);
        rowidx = 1;
        colidx = 1;
        v = 3.0;
        ierr   = MatSetValues(C,1,&rowidx,1,&colidx,&v,INSERT_VALUES);
        CHKERRQ(ierr);
    }
    ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
    CHKERRQ(ierr);
    ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
    CHKERRQ(ierr);

    /* create right hand side and solution */
    ierr = VecCreate(PETSC_COMM_WORLD,&u);
    CHKERRQ(ierr);
    ierr = VecSetSizes(u,PETSC_DECIDE,N);
    CHKERRQ(ierr);
    ierr = VecSetFromOptions(u);
    CHKERRQ(ierr);
    ierr = VecDuplicate(u,&b);
    CHKERRQ(ierr);
    ierr = VecDuplicate(u,&r);
    CHKERRQ(ierr);
    ierr = VecSet(u,0.0);
    CHKERRQ(ierr);
    ierr = VecSet(b,1.0);
    CHKERRQ(ierr);

    /* solve linear system C*u = b */
    ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);
    CHKERRQ(ierr);
    ierr = KSPSetOperators(ksp,C,C);
    CHKERRQ(ierr);
    ierr = KSPSetFromOptions(ksp);
    CHKERRQ(ierr);
    ierr = KSPSolve(ksp,b,u);
    CHKERRQ(ierr);

    /* check residual r = C*u - b */
    ierr = MatMult(C,u,r);
    CHKERRQ(ierr);
    ierr = VecAXPY(r,-1.0,b);
    CHKERRQ(ierr);
    ierr = VecNorm(r,NORM_2,&norm);
    CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"|| C*u - b|| = %g\n",(double)norm);
    CHKERRQ(ierr);

    /* solve C^T*u = b twice */
    ierr = KSPSolveTranspose(ksp,b,u);
    CHKERRQ(ierr);
    /* check residual r = C^T*u - b */
    ierr = MatMultTranspose(C,u,r);
    CHKERRQ(ierr);
    ierr = VecAXPY(r,-1.0,b);
    CHKERRQ(ierr);
    ierr = VecNorm(r,NORM_2,&norm);
    CHKERRQ(ierr);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    ierr = MatAssemblyEnd(C1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

    /*
       Indicate same nonzero structure of successive linear system matrices
    */
    ierr = MatSetOption(C1,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);CHKERRQ(ierr);

    /*
       Compute right-hand-side vector
    */
    ierr = MatMult(C1,u,b1);CHKERRQ(ierr);

    /*
       Set operators. Here the matrix that defines the linear system
       also serves as the preconditioning matrix.
        - The flag SAME_NONZERO_PATTERN indicates that the
          preconditioning matrix has identical nonzero structure
          as during the last linear solve (although the values of
          the entries have changed). Thus, we can save some
          work in setting up the preconditioner (e.g., no need to
          redo symbolic factorization for ILU/ICC preconditioners).
        - If the nonzero structure of the matrix is different during
          the second linear solve, then the flag DIFFERENT_NONZERO_PATTERN
          must be used instead.  If you are unsure whether the
          matrix structure has changed or not, use the flag
          DIFFERENT_NONZERO_PATTERN.
        - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
          believes your assertion and does not check the structure
          of the matrix.  If you erroneously claim that the structure
          is the same when it actually is not, the new preconditioner
          will not function correctly.  Thus, use this optimization
          feature with caution!
    */
    ierr = KSPSetOperators(ksp1,C1,C1,SAME_NONZERO_PATTERN);CHKERRQ(ierr);

    /*
       Use the previous solution of linear system #1 as the initial
       guess for the next solve of linear system #1.  The user MUST
       call KSPSetInitialGuessNonzero() in indicate use of an initial
       guess vector; otherwise, an initial guess of zero is used.
    */
    if (t>0) {
      ierr = KSPSetInitialGuessNonzero(ksp1,PETSC_TRUE);CHKERRQ(ierr);
    }

    /*
       Solve the first linear system.  Here we explicitly call
       KSPSetUp() for more detailed performance monitoring of
       certain preconditioners, such as ICC and ILU.  This call
       is optional, ase KSPSetUp() will automatically be called
       within KSPSolve() if it hasn't been called already.
    */
    ierr = KSPSetUp(ksp1);CHKERRQ(ierr);
    ierr = KSPSolve(ksp1,b1,x1);CHKERRQ(ierr);
    ierr = KSPGetIterationNumber(ksp1,&its);CHKERRQ(ierr);

    /*
       Check error of solution to first linear system
    */
    ierr = CheckError(u,x1,b1,its,1.e-4,CHECK_ERROR);CHKERRQ(ierr);

    /* - - - - - - - - - - - - Stage 2: - - - - - - - - - - - - - -
                 Assemble and solve second linear system
       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    /*

  long internal_solve()
  {

    /* 
       Create linear solver context
    */
    ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
    ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
    ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
    ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);
    ierr = KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
    ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
 
    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
       Solve the linear system
       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 

    PetscScalar*   solution_vector;
    VecGetArray(x,&solution_vector);
    for (long i =0; i < n; ++i)
      {
	sol_vec.push_back(solution_vector[i]);
      }

    int reason=0;
    KSPGetConvergedReason(ksp,(KSPConvergedReason*)&reason);
    switch(reason)
      {
      case 2:
	std::cout << "### CONVERGED ### " << std::endl;
	std::cout << "### KSP ..KSP_CONVERGED_RTOL .. " << std::endl;
	break;
      case 3:
	std::cout << "### CONVERGED ### " << std::endl;
	std::cout << "### KSP ..KSP_CONVERGED_ATOL .. " << std::endl;
	break;
      case 4:
	std::cout << "### CONVERGED ### " << std::endl;
	std::cout << "### KSP ..KSP_CONVERGED_ITS .. " << std::endl;
	break;
      case 5:
	std::cout << "### CONVERGED ### " << std::endl;
	std::cout << "### KSP ..KSP_CONVERGED_QCG_NEG_CURVE .. " << std::endl;
	break;
      case 6:
	std::cout << "### CONVERGED ### " << std::endl;
	std::cout << "### KSP ..KSP_CONVERGED_QCG_CONSTRAINED .. " << std::endl;
	break;
      case 7:
	std::cout << "### CONVERGED ### " << std::endl;
	std::cout << "### KSP ..KSP_CONVERGED_STEP_LENGTH .. " << std::endl;
	break;
      case -3:
	std::cout << "### !! N O T  !! CONVERGED ### " << std::endl;
	std::cout << "### KSP ..KSP_DIVERGED_ITS .. " << std::endl;
	break;
      case -4:
	std::cout << "### !! N O T  !! CONVERGED ### " << std::endl;
	std::cout << "### KSP ..KSP_DIVERGED_DTOL .. " << std::endl;
	break;
      case -5:
	std::cout << "### !! N O T  !! CONVERGED ### " << std::endl;
	std::cout << "### KSP ..KSP_DIVERGED_BREAKDOWN  .. " << std::endl;
	break;
      case -6:
	std::cout << "### !! N O T  !! CONVERGED ### " << std::endl;
	std::cout << "### KSP ..KSP_DIVERGED_BREAKDOWN_BICG .. " << std::endl;
	break;
      case -7:
	std::cout << "### !! N O T  !! CONVERGED ### " << std::endl;
	std::cout << "### KSP ..KSP_DIVERGED_NONSYMMETRIC .. " << std::endl;
	break;
      case -8:
	std::cout << "### !! N O T  !! CONVERGED ### " << std::endl;
	std::cout << "### KSP ..KSP_DIVERGED_INDEFINITE_PC .. " << std::endl;
	break;

      default:
	std::cout << "### KSP .. no problem detected.. " << std::endl;
	break;

      }

    

    /* 
    ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);



       Check the error
    ierr = VecAXPY(x,neg_one,u);CHKERRQ(ierr);
    ierr  = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
    ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Iterations %D\n",
		       norm,its);CHKERRQ(ierr);
    */
    return 0;
  }