C++ MatDestroy函数代码示例

int main(int argc,char **args)
{
  Mat            A,B;
  PetscErrorCode ierr;
  char           file[PETSC_MAX_PATH_LEN];
  PetscBool      flg;
  PetscViewer    fd;
  MatType        type = MATMPIBAIJ;

  PetscInitialize(&argc,&args,(char*)0,help);
  ierr = PetscOptionsHasName(NULL,"-aij",&flg);CHKERRQ(ierr);
  if (flg) type = MATMPIAIJ;

  ierr = PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");

  /*
     Open binary file.  Note that we use FILE_MODE_READ to indicate
     reading from this file.
  */
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr);

  /*
     Load the matrix; then destroy the viewer.
  */
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetType(A,type);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatLoad(A,fd);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);

  /*
     Open another binary file.  Note that we use FILE_MODE_WRITE to indicate
     reading from this file.
  */
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"fileoutput",FILE_MODE_WRITE,&fd);CHKERRQ(ierr);
  ierr = PetscViewerBinarySetFlowControl(fd,3);CHKERRQ(ierr);
  /*
     Save the matrix and vector; then destroy the viewer.
  */
  ierr = MatView(A,fd);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);

  /* load the new matrix */
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"fileoutput",FILE_MODE_READ,&fd);CHKERRQ(ierr);
  ierr = MatCreate(PETSC_COMM_WORLD,&B);CHKERRQ(ierr);
  ierr = MatSetType(B,type);CHKERRQ(ierr);
  ierr = MatSetFromOptions(B);CHKERRQ(ierr);
  ierr = MatLoad(B,fd);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);

  ierr = MatEqual(A,B,&flg);CHKERRQ(ierr);
  if (flg) {
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Matrices are equal\n");CHKERRQ(ierr);
  } else {
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Matrices are not equal\n");CHKERRQ(ierr);
  }

  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = MatDestroy(&B);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}

//.........这里部分代码省略.........
  ierr = PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);CHKERRQ(ierr);
  n = m;
  ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);

  ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr);
  ierr = MatSetSizes(C,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
  ierr = MatSetType(C,MATELEMENTAL);CHKERRQ(ierr);
  ierr = MatSetFromOptions(C);CHKERRQ(ierr);
  ierr = MatSetUp(C);CHKERRQ(ierr);

  ierr = PetscOptionsHasName(PETSC_NULL,"-row_oriented",&flg);CHKERRQ(ierr);
  if (flg) {ierr = MatSetOption(C,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr);}
  ierr = MatGetOwnershipIS(C,&isrows,&iscols);CHKERRQ(ierr);
  ierr = PetscOptionsHasName(PETSC_NULL,"-Cexp_view_ownership",&flg);CHKERRQ(ierr);
  if (flg) { // View ownership of explicit C
    IS tmp;
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Ownership of explicit C:\n");CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Row index set:\n");CHKERRQ(ierr);
    ierr = ISOnComm(isrows,PETSC_COMM_WORLD,PETSC_USE_POINTER,&tmp);CHKERRQ(ierr);
    ierr = ISView(tmp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = ISDestroy(&tmp);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Column index set:\n");CHKERRQ(ierr);
    ierr = ISOnComm(iscols,PETSC_COMM_WORLD,PETSC_USE_POINTER,&tmp);CHKERRQ(ierr);
    ierr = ISView(tmp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = ISDestroy(&tmp);CHKERRQ(ierr);
  }

  // Set local matrix entries
  ierr = ISGetLocalSize(isrows,&nrows);CHKERRQ(ierr);
  ierr = ISGetIndices(isrows,&rows);CHKERRQ(ierr);
  ierr = ISGetLocalSize(iscols,&ncols);CHKERRQ(ierr);
  ierr = ISGetIndices(iscols,&cols);CHKERRQ(ierr);
  ierr = PetscMalloc(nrows*ncols*sizeof(*v),&v);CHKERRQ(ierr);
  for (i=0; i<nrows; i++) {
    for (j=0; j<ncols; j++) {
      //v[i*ncols+j] = (PetscReal)(rank);
      v[i*ncols+j] = (PetscReal)(rank*10000+100*rows[i]+cols[j]);
    }
  }
  ierr = MatSetValues(C,nrows,rows,ncols,cols,v,INSERT_VALUES);CHKERRQ(ierr);
  ierr = ISRestoreIndices(isrows,&rows);CHKERRQ(ierr);
  ierr = ISRestoreIndices(iscols,&cols);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  // Test MatView()
  if (mats_view){
    ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  }

  // Set unowned matrix entries - add subdiagonals and diagonals from proc[0]
  if (rank == 0) {
    PetscInt M,N,cols[2];
    ierr = MatGetSize(C,&M,&N);CHKERRQ(ierr);
    for (i=0; i<M; i++){
      cols[0] = i;   v[0] = i + 0.5;
      cols[1] = i-1; v[1] = 0.5;
      if (i) {
        ierr = MatSetValues(C,1,&i,2,cols,v,ADD_VALUES);CHKERRQ(ierr);
      } else {
        ierr = MatSetValues(C,1,&i,1,&i,v,ADD_VALUES);CHKERRQ(ierr);
      }
    }
  }
  ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  // Test MatMult()
  ierr = MatComputeExplicitOperator(C,&Caij);CHKERRQ(ierr);
  ierr = MatMultEqual(C,Caij,5,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultEqual() fails");
  ierr = MatMultTransposeEqual(C,Caij,5,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultTransposeEqual() fails");

  // Test MatMultAdd() and MatMultTransposeAddEqual()
  ierr = MatMultAddEqual(C,Caij,5,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultAddEqual() fails");
  ierr = MatMultTransposeAddEqual(C,Caij,5,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultTransposeAddEqual() fails");

  // Test MatMatMult()
  ierr = PetscOptionsHasName(PETSC_NULL,"-test_matmatmult",&Test_MatMatMult);CHKERRQ(ierr);
  if (Test_MatMatMult){
    Mat CCelem,CCaij;
    ierr = MatMatMult(C,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CCelem);CHKERRQ(ierr);
    ierr = MatMatMult(Caij,Caij,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CCaij);CHKERRQ(ierr);
    ierr = MatMultEqual(CCelem,CCaij,5,&flg);CHKERRQ(ierr);
    if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"CCelem != CCaij. MatMatMult() fails");
    ierr = MatDestroy(&CCaij);CHKERRQ(ierr);
    ierr = MatDestroy(&CCelem);CHKERRQ(ierr);
  }

  ierr = PetscFree(v);CHKERRQ(ierr);
  ierr = ISDestroy(&isrows);CHKERRQ(ierr);
  ierr = ISDestroy(&iscols);CHKERRQ(ierr);
  ierr = MatDestroy(&C);CHKERRQ(ierr);
  ierr = MatDestroy(&Caij);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}

PetscErrorCode test1_DAInjection3d(PetscInt mx, PetscInt my, PetscInt mz)
{
  PetscErrorCode   ierr;
  DM               dac,daf;
  PetscViewer      vv;
  Vec              ac,af;
  PetscInt         periodicity;
  DMBoundaryType   bx,by,bz;

  PetscFunctionBeginUser;
  bx = DM_BOUNDARY_NONE;
  by = DM_BOUNDARY_NONE;
  bz = DM_BOUNDARY_NONE;

  periodicity = 0;

  ierr = PetscOptionsGetInt(NULL,"-periodic", &periodicity, NULL);CHKERRQ(ierr);
  if (periodicity==1) {
    bx = DM_BOUNDARY_PERIODIC;
  } else if (periodicity==2) {
    by = DM_BOUNDARY_PERIODIC;
  } else if (periodicity==3) {
    bz = DM_BOUNDARY_PERIODIC;
  }

  ierr = DMDACreate3d(PETSC_COMM_WORLD, bx,by,bz, DMDA_STENCIL_BOX,
                      mx+1, my+1,mz+1,
                      PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE,
                      1, /* 1 dof */
                      1, /* stencil = 1 */
                      NULL,NULL,NULL,
                      &daf);CHKERRQ(ierr);

  ierr = DMSetFromOptions(daf);CHKERRQ(ierr);

  ierr = DMCoarsen(daf,MPI_COMM_NULL,&dac);CHKERRQ(ierr);

  ierr = DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, -1.0,1.0);CHKERRQ(ierr);
  ierr = DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, -1.0,1.0);CHKERRQ(ierr);

  {
    DM         cdaf,cdac;
    Vec        coordsc,coordsf,coordsf2;
    VecScatter inject;
    Mat        interp;
    PetscReal  norm;

    ierr = DMGetCoordinateDM(dac,&cdac);CHKERRQ(ierr);
    ierr = DMGetCoordinateDM(daf,&cdaf);CHKERRQ(ierr);

    ierr = DMGetCoordinates(dac,&coordsc);CHKERRQ(ierr);
    ierr = DMGetCoordinates(daf,&coordsf);CHKERRQ(ierr);

    ierr = DMCreateInjection(cdac,cdaf,&inject);CHKERRQ(ierr);

    ierr = VecScatterBegin(inject,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = VecScatterEnd(inject  ,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = VecScatterDestroy(&inject);CHKERRQ(ierr);

    ierr = DMCreateInterpolation(cdac,cdaf,&interp,NULL);CHKERRQ(ierr);
    ierr = VecDuplicate(coordsf,&coordsf2);CHKERRQ(ierr);
    ierr = MatInterpolate(interp,coordsc,coordsf2);CHKERRQ(ierr);
    ierr = VecAXPY(coordsf2,-1.0,coordsf);CHKERRQ(ierr);
    ierr = VecNorm(coordsf2,NORM_MAX,&norm);CHKERRQ(ierr);
    /* The fine coordinates are only reproduced in certain cases */
    if (!bx && !by && !bz && norm > 1.e-10) {ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm %g\n",(double)norm);CHKERRQ(ierr);}
    ierr = VecDestroy(&coordsf2);CHKERRQ(ierr);
    ierr = MatDestroy(&interp);CHKERRQ(ierr);
  }

  if (0) {
    ierr = DMCreateGlobalVector(dac,&ac);CHKERRQ(ierr);
    ierr = VecZeroEntries(ac);CHKERRQ(ierr);

    ierr = DMCreateGlobalVector(daf,&af);CHKERRQ(ierr);
    ierr = VecZeroEntries(af);CHKERRQ(ierr);

    ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_7.vtk", &vv);CHKERRQ(ierr);
    ierr = PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
    ierr = DMView(dac, vv);CHKERRQ(ierr);
    ierr = VecView(ac, vv);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&vv);CHKERRQ(ierr);

    ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_7.vtk", &vv);CHKERRQ(ierr);
    ierr = PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
    ierr = DMView(daf, vv);CHKERRQ(ierr);
    ierr = VecView(af, vv);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&vv);CHKERRQ(ierr);
    ierr = VecDestroy(&ac);CHKERRQ(ierr);
    ierr = VecDestroy(&af);CHKERRQ(ierr);
  }
  ierr = DMDestroy(&dac);CHKERRQ(ierr);
  ierr = DMDestroy(&daf);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode port_lsd_bfbt(void)
{
  Mat            A;
  Vec            x,b;
  KSP            ksp_A;
  PC             pc_A;
  IS             isu,isp;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = LoadTestMatrices(&A,&x,&b,&isu,&isp);CHKERRQ(ierr);

  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp_A);CHKERRQ(ierr);
  ierr = KSPSetOptionsPrefix(ksp_A,"fc_");CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp_A,A,A);CHKERRQ(ierr);

  ierr = KSPGetPC(ksp_A,&pc_A);CHKERRQ(ierr);
  ierr = PCSetType(pc_A,PCFIELDSPLIT);CHKERRQ(ierr);
  ierr = PCFieldSplitSetBlockSize(pc_A,2);CHKERRQ(ierr);
  ierr = PCFieldSplitSetIS(pc_A,"velocity",isu);CHKERRQ(ierr);
  ierr = PCFieldSplitSetIS(pc_A,"pressure",isp);CHKERRQ(ierr);

  ierr = KSPSetFromOptions(ksp_A);CHKERRQ(ierr);
  ierr = KSPSolve(ksp_A,b,x);CHKERRQ(ierr);

  /* Pull u,p out of x */
  {
    PetscInt    loc;
    PetscReal   max,norm;
    PetscScalar sum;
    Vec         uvec,pvec;
    VecScatter  uscat,pscat;
    Mat         A11,A22;

    /* grab matrices and create the compatable u,p vectors */
    ierr = MatCreateSubMatrix(A,isu,isu,MAT_INITIAL_MATRIX,&A11);CHKERRQ(ierr);
    ierr = MatCreateSubMatrix(A,isp,isp,MAT_INITIAL_MATRIX,&A22);CHKERRQ(ierr);

    ierr = MatCreateVecs(A11,&uvec,NULL);CHKERRQ(ierr);
    ierr = MatCreateVecs(A22,&pvec,NULL);CHKERRQ(ierr);

    /* perform the scatter from x -> (u,p) */
    ierr = VecScatterCreateWithData(x,isu,uvec,NULL,&uscat);CHKERRQ(ierr);
    ierr = VecScatterBegin(uscat,x,uvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = VecScatterEnd(uscat,x,uvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);

    ierr = VecScatterCreateWithData(x,isp,pvec,NULL,&pscat);CHKERRQ(ierr);
    ierr = VecScatterBegin(pscat,x,pvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = VecScatterEnd(pscat,x,pvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);

    PetscPrintf(PETSC_COMM_WORLD,"-- vector vector values --\n");
    ierr = VecMin(uvec,&loc,&max);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"  Min(u)  = %1.6f [loc=%D]\n",(double)max,loc);CHKERRQ(ierr);
    ierr = VecMax(uvec,&loc,&max);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"  Max(u)  = %1.6f [loc=%D]\n",(double)max,loc);CHKERRQ(ierr);
    ierr = VecNorm(uvec,NORM_2,&norm);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"  Norm(u) = %1.6f \n",(double)norm);CHKERRQ(ierr);
    ierr = VecSum(uvec,&sum);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"  Sum(u)  = %1.6f \n",(double)PetscRealPart(sum));CHKERRQ(ierr);

    PetscPrintf(PETSC_COMM_WORLD,"-- pressure vector values --\n");
    ierr = VecMin(pvec,&loc,&max);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"  Min(p)  = %1.6f [loc=%D]\n",(double)max,loc);CHKERRQ(ierr);
    ierr = VecMax(pvec,&loc,&max);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"  Max(p)  = %1.6f [loc=%D]\n",(double)max,loc);CHKERRQ(ierr);
    ierr = VecNorm(pvec,NORM_2,&norm);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"  Norm(p) = %1.6f \n",(double)norm);CHKERRQ(ierr);
    ierr = VecSum(pvec,&sum);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"  Sum(p)  = %1.6f \n",(double)PetscRealPart(sum));CHKERRQ(ierr);

    PetscPrintf(PETSC_COMM_WORLD,"-- Full vector values --\n");
    ierr = VecMin(x,&loc,&max);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"  Min(u,p)  = %1.6f [loc=%D]\n",(double)max,loc);CHKERRQ(ierr);
    ierr = VecMax(x,&loc,&max);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"  Max(u,p)  = %1.6f [loc=%D]\n",(double)max,loc);CHKERRQ(ierr);
    ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"  Norm(u,p) = %1.6f \n",(double)norm);CHKERRQ(ierr);
    ierr = VecSum(x,&sum);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"  Sum(u,p)  = %1.6f \n",(double)PetscRealPart(sum));CHKERRQ(ierr);

    ierr = VecScatterDestroy(&uscat);CHKERRQ(ierr);
    ierr = VecScatterDestroy(&pscat);CHKERRQ(ierr);
    ierr = VecDestroy(&uvec);CHKERRQ(ierr);
    ierr = VecDestroy(&pvec);CHKERRQ(ierr);
    ierr = MatDestroy(&A11);CHKERRQ(ierr);
    ierr = MatDestroy(&A22);CHKERRQ(ierr);
  }

  ierr = KSPDestroy(&ksp_A);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr);
  ierr = ISDestroy(&isu);CHKERRQ(ierr);
  ierr = ISDestroy(&isp);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

//.........这里部分代码省略.........
        col = row+4; ierr = MatSetValues(A,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr);
      }
      v    = 4.0;
      ierr = MatSetValues(A,1,&row,1,&row,&v,INSERT_VALUES);CHKERRQ(ierr);
    }
  }
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD, "Original matrix\n");CHKERRQ(ierr);
  ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  if (size > 1) {
    nsub = 1; /* one subdomain per rank */
  }
  else {
    nsub = 2; /* both subdomains on rank 0 */
  }
  if (rank) {
    jlow = Jlow+1; jhigh = Jhigh+1;
  }
  else {
    jlow = Jlow; jhigh = Jhigh;
  }
  sort_rows = PETSC_FALSE;
  ierr      = PetscOptionsGetBool(NULL,NULL, "-sort_rows", &sort_rows, NULL);CHKERRQ(ierr);
  sort_cols = PETSC_FALSE;
  ierr      = PetscOptionsGetBool(NULL,NULL, "-sort_cols", &sort_cols, NULL);CHKERRQ(ierr);
  for (l = 0; l < nsub; ++l) {
    ierr = PetscMalloc1(12, &subindices);CHKERRQ(ierr);
    k    = 0;
    for (i = 0; i < 4; ++i) {
      for (j = jlow[l]; j < jhigh[l]; ++j) {
        subindices[k] = j*4+i;
        k++;
      }
    }
    ierr = ISCreateGeneral(PETSC_COMM_SELF, 12, subindices, PETSC_OWN_POINTER, rowis+l);CHKERRQ(ierr);
    if ((sort_rows && !sort_cols) || (!sort_rows && sort_cols)) {
      ierr = ISDuplicate(rowis[l],colis+l);CHKERRQ(ierr);
    } else {
      ierr = PetscObjectReference((PetscObject)rowis[l]);CHKERRQ(ierr);
      colis[l] = rowis[l];
    }
    if (sort_rows) {
      ierr = ISSort(rowis[l]);CHKERRQ(ierr);
    }
    if (sort_cols) {
      ierr = ISSort(colis[l]);CHKERRQ(ierr);
    }
  }
  ierr = PetscMalloc1(nsub, &S);CHKERRQ(ierr);
  ierr = MatGetSubMatrices(A,nsub,rowis,colis,MAT_INITIAL_MATRIX, &S);CHKERRQ(ierr);

  show_inversions = PETSC_FALSE;

  ierr = PetscOptionsGetBool(NULL,NULL, "-show_inversions", &show_inversions, NULL);CHKERRQ(ierr);

  inversions = 0;
  for (p = 0; p < size; ++p) {
    if (p == rank) {
      ierr = PetscPrintf(PETSC_COMM_SELF, "[%D:%D]: Number of subdomains: %D:\n", rank, size, nsub);CHKERRQ(ierr);
      for (l = 0; l < nsub; ++l) {
        PetscInt i0, i1;
        ierr = PetscPrintf(PETSC_COMM_SELF, "[%D:%D]: Subdomain row IS %D:\n", rank, size, l);CHKERRQ(ierr);
        ierr = ISView(rowis[l],PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, "[%D:%D]: Subdomain col IS %D:\n", rank, size, l);CHKERRQ(ierr);
        ierr = ISView(colis[l],PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, "[%D:%D]: Submatrix %D:\n", rank, size, l);CHKERRQ(ierr);
        ierr = MatView(S[l],PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
        if (show_inversions) {
          ierr = MatGetOwnershipRange(S[l], &i0,&i1);CHKERRQ(ierr);
          for (i = i0; i < i1; ++i) {
            ierr = MatGetRow(S[l], i, &ncols, &cols, NULL);CHKERRQ(ierr);
            for (j = 1; j < ncols; ++j) {
              if (cols[j] < cols[j-1]) {
                ierr = PetscPrintf(PETSC_COMM_SELF, "***Inversion in row %D: col[%D] = %D < %D = col[%D]\n", i, j, cols[j], cols[j-1], j-1);CHKERRQ(ierr);
                inversions++;
              }
            }
            ierr = MatRestoreRow(S[l], i, &ncols, &cols, NULL);CHKERRQ(ierr);
          }
        }
      }
    }
    ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRQ(ierr);
  }
  if (show_inversions) {
    ierr = MPI_Reduce(&inversions,&total_inversions,1,MPIU_INT, MPI_SUM,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "*Total inversions: %D\n", total_inversions);CHKERRQ(ierr);
  }
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  for (l = 0; l < nsub; ++l) {
    ierr = MatDestroy(&(S[l]));CHKERRQ(ierr);
    ierr = ISDestroy(&(rowis[l]));CHKERRQ(ierr);
    ierr = ISDestroy(&(colis[l]));CHKERRQ(ierr);
  }
  ierr = PetscFree(S);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return ierr;
}

//.........这里部分代码省略.........
      const PetscInt    *cone;
      PetscInt           coneSize, cl, i, j, d, r;

      ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);CHKERRQ(ierr);
      /* Only works for simplices */
      for (i = 0, r = 0; i < dim+1; ++i) {
        for (j = 0; j < i; ++j, ++r) {
          for (d = 0; d < dim; ++d) e[d] = cellCoords[i*dim+d] - cellCoords[j*dim+d];
          /* FORTRAN ORDERING */
          if (dim == 2) {
            eqns[0*size+r] = PetscSqr(e[0]);
            eqns[1*size+r] = 2.0*e[0]*e[1];
            eqns[2*size+r] = PetscSqr(e[1]);
          } else {
            eqns[0*size+r] = PetscSqr(e[0]);
            eqns[1*size+r] = 2.0*e[0]*e[1];
            eqns[2*size+r] = 2.0*e[0]*e[2];
            eqns[3*size+r] = PetscSqr(e[1]);
            eqns[4*size+r] = 2.0*e[1]*e[2];
            eqns[5*size+r] = PetscSqr(e[2]);
          }
        }
      }
      ierr = MatSetUnfactored(A);CHKERRQ(ierr);
      ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);CHKERRQ(ierr);
      ierr = MatLUFactor(A, NULL, NULL, NULL);CHKERRQ(ierr);
      ierr = MatSolve(A, mb, mx);CHKERRQ(ierr);
      ierr = VecGetArrayRead(mx, &sol);CHKERRQ(ierr);
      ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);CHKERRQ(ierr);
      ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr);
      ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr);
      for (cl = 0; cl < coneSize; ++cl) {
        const PetscInt v = cone[cl] - vStart;

        if (dim == 2) {
          metric[v*4+0] += vol*coarseRatio*sol[0];
          metric[v*4+1] += vol*coarseRatio*sol[1];
          metric[v*4+2] += vol*coarseRatio*sol[1];
          metric[v*4+3] += vol*coarseRatio*sol[2];
        } else {
          metric[v*9+0] += vol*coarseRatio*sol[0];
          metric[v*9+1] += vol*coarseRatio*sol[1];
          metric[v*9+3] += vol*coarseRatio*sol[1];
          metric[v*9+2] += vol*coarseRatio*sol[2];
          metric[v*9+6] += vol*coarseRatio*sol[2];
          metric[v*9+4] += vol*coarseRatio*sol[3];
          metric[v*9+5] += vol*coarseRatio*sol[4];
          metric[v*9+7] += vol*coarseRatio*sol[4];
          metric[v*9+8] += vol*coarseRatio*sol[5];
        }
      }
      ierr = VecRestoreArrayRead(mx, &sol);CHKERRQ(ierr);
    }
    for (v = 0; v < numVertices; ++v) {
      const PetscInt *support;
      PetscInt        supportSize, s;
      PetscReal       vol, totVol = 0.0;

      ierr = DMPlexGetSupport(udm, v+vStart, &support);CHKERRQ(ierr);
      ierr = DMPlexGetSupportSize(udm, v+vStart, &supportSize);CHKERRQ(ierr);
      for (s = 0; s < supportSize; ++s) {ierr = DMPlexComputeCellGeometryFVM(dm, support[s], &vol, NULL, NULL);CHKERRQ(ierr); totVol += vol;}
      for (s = 0; s < PetscSqr(dim); ++s) metric[v*PetscSqr(dim)+s] /= totVol;
    }
    ierr = VecDestroy(&mx);CHKERRQ(ierr);
    ierr = VecDestroy(&mb);CHKERRQ(ierr);
    ierr = MatDestroy(&A);CHKERRQ(ierr);
    ierr = DMDestroy(&udm);CHKERRQ(ierr);
    ierr = PetscFree(eqns);CHKERRQ(ierr);
    pragmatic_set_metric(metric);
    pragmatic_adapt();
    /* Read out mesh */
    pragmatic_get_info(&numCoarseVertices, &numCoarseCells);
    ierr = PetscMalloc1(numCoarseVertices*dim, &coarseCoords);CHKERRQ(ierr);
    switch (dim) {
    case 2:
      pragmatic_get_coords_2d(x, y);
      numCorners = 3;
      for (v = 0; v < numCoarseVertices; ++v) {coarseCoords[v*2+0] = x[v]; coarseCoords[v*2+1] = y[v];}
      break;
    case 3:
      pragmatic_get_coords_3d(x, y, z);
      numCorners = 4;
      for (v = 0; v < numCoarseVertices; ++v) {coarseCoords[v*3+0] = x[v]; coarseCoords[v*3+1] = y[v]; coarseCoords[v*3+2] = z[v];}
      break;
    default:
      SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic coarsening defined for dimension %d", dim);
    }
    pragmatic_get_elements(cells);
    /* TODO Read out markers for boundary */
    ierr = DMPlexCreateFromCellList(PetscObjectComm((PetscObject) dm), dim, numCoarseCells, numCoarseVertices, numCorners, PETSC_TRUE, cells, dim, coarseCoords, &mesh->coarseMesh);CHKERRQ(ierr);
    pragmatic_finalize();
    ierr = PetscFree5(x, y, z, metric, cells);CHKERRQ(ierr);
    ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr);
    ierr = PetscFree(coarseCoords);CHKERRQ(ierr);
  }
#endif
  ierr = PetscObjectReference((PetscObject) mesh->coarseMesh);CHKERRQ(ierr);
  *dmCoarsened = mesh->coarseMesh;
  PetscFunctionReturn(0);
}

//.........这里部分代码省略.........
        }
        ierr = VecDestroy(&b);CHKERRQ(ierr);
      }
      ierr = VecDestroy(&b1);CHKERRQ(ierr);
      ierr = VecDestroy(&b2);CHKERRQ(ierr);
      ierr = VecDestroy(&x);CHKERRQ(ierr);
      ierr = MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr);

      ierr = MatCreateTranspose(S,&ST);CHKERRQ(ierr);
      ierr = MatComputeOperator(ST,MATDENSE,&BT);CHKERRQ(ierr);
      ierr = MatTranspose(BT,MAT_INITIAL_MATRIX,&BTT);CHKERRQ(ierr);
      ierr = PetscObjectSetName((PetscObject)B,"S");CHKERRQ(ierr);
      ierr = PetscObjectSetName((PetscObject)BTT,"STT");CHKERRQ(ierr);
      ierr = MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);
      ierr = PetscObjectSetName((PetscObject)C,"A");CHKERRQ(ierr);

      ierr = MatViewFromOptions(C,NULL,"-view_mat");CHKERRQ(ierr);
      ierr = MatViewFromOptions(B,NULL,"-view_mat");CHKERRQ(ierr);
      ierr = MatViewFromOptions(BTT,NULL,"-view_mat");CHKERRQ(ierr);

      ierr = MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
      ierr = MatNorm(C,NORM_FROBENIUS,&err);CHKERRQ(ierr);
      if (err >= tol) {
        ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mat mult %g\n",test,(double)err);CHKERRQ(ierr);
      }

      ierr = MatConvert(A,MATDENSE,MAT_REUSE_MATRIX,&C);CHKERRQ(ierr);
      ierr = MatAXPY(C,-1.0,BTT,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
      ierr = MatNorm(C,NORM_FROBENIUS,&err);CHKERRQ(ierr);
      if (err >= tol) {
        ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mat mult transpose %g\n",test,(double)err);CHKERRQ(ierr);
      }

      ierr = MatDestroy(&ST);CHKERRQ(ierr);
      ierr = MatDestroy(&BTT);CHKERRQ(ierr);
      ierr = MatDestroy(&BT);CHKERRQ(ierr);
      ierr = MatDestroy(&B);CHKERRQ(ierr);
      ierr = MatDestroy(&C);CHKERRQ(ierr);
    }
    if (testdiagscale) { /* MatDiagonalScale() */
      Vec vr,vl;

      ierr = MatCreateVecs(A,&vr,&vl);CHKERRQ(ierr);
      if (randomize) {
        ierr = VecSetRandom(vr,NULL);CHKERRQ(ierr);
        ierr = VecSetRandom(vl,NULL);CHKERRQ(ierr);
      } else {
        ierr = VecSet(vr,test%2 ? 0.15 : 1.0/0.15);CHKERRQ(ierr);
        ierr = VecSet(vl,test%2 ? -1.2 : 1.0/-1.2);CHKERRQ(ierr);
      }
      ierr = MatDiagonalScale(A,vl,vr);CHKERRQ(ierr);
      ierr = MatDiagonalScale(S,vl,vr);CHKERRQ(ierr);
      ierr = VecDestroy(&vr);CHKERRQ(ierr);
      ierr = VecDestroy(&vl);CHKERRQ(ierr);
    }

    if (testscale) { /* MatScale() */
      ierr = MatScale(A,test%2 ? 1.4 : 1.0/1.4);CHKERRQ(ierr);
      ierr = MatScale(S,test%2 ? 1.4 : 1.0/1.4);CHKERRQ(ierr);
    }

    if (testshift && cong) { /* MatShift() : MATSHELL shift is broken when row/cols layout are not congruent and left/right scaling have been applied */
      ierr = MatShift(A,test%2 ? -77.5 : 77.5);CHKERRQ(ierr);
      ierr = MatShift(S,test%2 ? -77.5 : 77.5);CHKERRQ(ierr);
    }

int main(int argc,char **args)
{
  MatType         mtype = MATMPIAIJ; /* matrix format */
  Mat             A,B;               /* matrix */
  PetscViewer     fd;                /* viewer */
  char            file[PETSC_MAX_PATH_LEN];         /* input file name */
  PetscBool       flg,viewMats,viewIS,viewVecs;
  PetscInt        ierr,*nlocal,m,n;
  PetscMPIInt     rank,size;
  MatPartitioning part;
  IS              is,isn;
  Vec             xin, xout;
  VecScatter      scat;

  PetscInitialize(&argc,&args,(char*)0,help);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
  ierr = PetscOptionsHasName(NULL, "-view_mats", &viewMats);CHKERRQ(ierr);
  ierr = PetscOptionsHasName(NULL, "-view_is", &viewIS);CHKERRQ(ierr);
  ierr = PetscOptionsHasName(NULL, "-view_vecs", &viewVecs);CHKERRQ(ierr);

  /*
     Determine file from which we read the matrix
  */
  ierr = PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);

  /*
       Open binary file.  Note that we use FILE_MODE_READ to indicate
       reading from this file.
  */
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr);

  /*
      Load the matrix and vector; then destroy the viewer.
  */
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetType(A,mtype);CHKERRQ(ierr);
  ierr = MatLoad(A,fd);CHKERRQ(ierr);
  ierr = VecCreate(PETSC_COMM_WORLD,&xin);CHKERRQ(ierr);
  ierr = VecLoad(xin,fd);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
  if (viewMats) {
    if (!rank) printf("Original matrix:\n");
    ierr = MatView(A,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
  }
  if (viewVecs) {
    if (!rank) printf("Original vector:\n");
    ierr = VecView(xin,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  }

  /* Partition the graph of the matrix */
  ierr = MatPartitioningCreate(PETSC_COMM_WORLD,&part);CHKERRQ(ierr);
  ierr = MatPartitioningSetAdjacency(part,A);CHKERRQ(ierr);
  ierr = MatPartitioningSetFromOptions(part);CHKERRQ(ierr);

  /* get new processor owner number of each vertex */
  ierr = MatPartitioningApply(part,&is);CHKERRQ(ierr);
  if (viewIS) {
    if (!rank) printf("IS1 - new processor ownership:\n");
    ierr = ISView(is,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  }

  /* get new global number of each old global number */
  ierr = ISPartitioningToNumbering(is,&isn);CHKERRQ(ierr);
  if (viewIS) {
    if (!rank) printf("IS2 - new global numbering:\n");
    ierr = ISView(isn,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  }

  /* get number of new vertices for each processor */
  ierr = PetscMalloc(size*sizeof(PetscInt),&nlocal);CHKERRQ(ierr);
  ierr = ISPartitioningCount(is,size,nlocal);CHKERRQ(ierr);
  ierr = ISDestroy(&is);CHKERRQ(ierr);

  /* get old global number of each new global number */
  ierr = ISInvertPermutation(isn,nlocal[rank],&is);CHKERRQ(ierr);
  ierr = PetscFree(nlocal);CHKERRQ(ierr);
  ierr = ISDestroy(&isn);CHKERRQ(ierr);
  ierr = MatPartitioningDestroy(&part);CHKERRQ(ierr);
  if (viewIS) {
    if (!rank) printf("IS3=inv(IS2) - old global number of each new global number:\n");
    ierr = ISView(is,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  }

  /* move the matrix rows to the new processes they have been assigned to by the permutation */
  ierr = ISSort(is);CHKERRQ(ierr);
  ierr = MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);

  /* move the vector rows to the new processes they have been assigned to */
  ierr = MatGetLocalSize(B,&m,&n);CHKERRQ(ierr);
  ierr = VecCreateMPI(PETSC_COMM_WORLD,m,PETSC_DECIDE,&xout);CHKERRQ(ierr);
  ierr = VecScatterCreate(xin,is,xout,NULL,&scat);CHKERRQ(ierr);
  ierr = VecScatterBegin(scat,xin,xout,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterEnd(scat,xin,xout,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterDestroy(&scat);CHKERRQ(ierr);
  ierr = ISDestroy(&is);CHKERRQ(ierr);
  if (viewMats) {
    if (!rank) printf("Partitioned matrix:\n");
    ierr = MatView(B,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
//.........这里部分代码省略.........

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      preonly,lu,redundant,cholesky,svd,dump = PETSC_FALSE,opsset,use_amat;
  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);
    }
  }

  /* Skipping this for galerkin==2 (externally managed hierarchy such as ML and GAMG). Cleaner logic here would be great. Wrap ML/GAMG as DMs? */
  if (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;
    for (i=n-2; i>-1; i--) {
      DMKSP kdm;
      ierr = DMCoarsen(dms[i+1],MPI_COMM_NULL,&dms[i]);CHKERRQ(ierr);
      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);
      }
    }

    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 */
    ierr = KSPSetDM(mglevels[n-1]->smoothd,pc->dm);CHKERRQ(ierr);
    ierr = KSPSetDMActive(mglevels[n-1]->smoothd,PETSC_FALSE);CHKERRQ(ierr);
  }

  if (mg->galerkin == 1) {
    Mat B;
    /* currently only handle case where mat and pmat are the same on coarser levels */
    ierr = KSPGetOperators(mglevels[n-1]->smoothd,&dA,&dB);CHKERRQ(ierr);
    if (!pc->setupcalled) {
      for (i=n-2; i>-1; i--) {
        if (mglevels[i+1]->interpolate == mglevels[i+1]->restrct || !mglevels[i+1]->restrct) {
          ierr = MatPtAP(dB,mglevels[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);CHKERRQ(ierr);
        } else {
          ierr = MatMatMatMult(mglevels[i+1]->restrct,dB,mglevels[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);CHKERRQ(ierr);
        }
//.........这里部分代码省略.........

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  KSP            ksp;
  PC             pc;
  Vec            x,b;
  DM             da;
  Mat            A;
  PetscInt       dof=1;
  PetscBool      flg;
  PetscScalar    zero=0.0;

  PetscInitialize(&argc,&argv,(char*)0,help);
  ierr = PetscOptionsGetInt(NULL,"-dof",&dof,NULL);CHKERRQ(ierr);

  ierr = DMDACreate(PETSC_COMM_WORLD,&da);CHKERRQ(ierr);
  ierr = DMDASetDim(da,3);CHKERRQ(ierr);
  ierr = DMDASetBoundaryType(da,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE);CHKERRQ(ierr);
  ierr = DMDASetStencilType(da,DMDA_STENCIL_STAR);CHKERRQ(ierr);
  ierr = DMDASetSizes(da,3,3,3);CHKERRQ(ierr);
  ierr = DMDASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
  ierr = DMDASetDof(da,dof);CHKERRQ(ierr);
  ierr = DMDASetStencilWidth(da,1);CHKERRQ(ierr);
  ierr = DMDASetOwnershipRanges(da,NULL,NULL,NULL);CHKERRQ(ierr);
  ierr = DMSetFromOptions(da);CHKERRQ(ierr);
  ierr = DMSetUp(da);CHKERRQ(ierr);

  ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr);
  ierr = DMCreateGlobalVector(da,&b);CHKERRQ(ierr);
  ierr = DMCreateMatrix(da,MATAIJ,&A);CHKERRQ(ierr);
  ierr = VecSet(b,zero);CHKERRQ(ierr);

  /* Test sbaij matrix */
  flg  = PETSC_FALSE;
  ierr = PetscOptionsGetBool(NULL,"-test_sbaij",&flg,NULL);CHKERRQ(ierr);
  if (flg) {
    Mat sA;
    ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
    ierr = MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);CHKERRQ(ierr);
    ierr = MatDestroy(&A);CHKERRQ(ierr);
    A    = sA;
  }

  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
  ierr = PCSetDM(pc,(DM)da);CHKERRQ(ierr);

  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);

  /* check final residual */
  flg  = PETSC_FALSE;
  ierr = PetscOptionsGetBool(NULL, "-check_final_residual", &flg,NULL);CHKERRQ(ierr);
  if (flg) {
    Vec       b1;
    PetscReal norm;
    ierr = KSPGetSolution(ksp,&x);CHKERRQ(ierr);
    ierr = VecDuplicate(b,&b1);CHKERRQ(ierr);
    ierr = MatMult(A,x,b1);CHKERRQ(ierr);
    ierr = VecAXPY(b1,-1.0,b);CHKERRQ(ierr);
    ierr = VecNorm(b1,NORM_2,&norm);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Final residual %g\n",norm);CHKERRQ(ierr);
    ierr = VecDestroy(&b1);CHKERRQ(ierr);
  }

  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = DMDestroy(&da);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}

int main(int argc,char **args)
{
  Vec            x,b,u;      /* approx solution, RHS, exact solution */
  Mat            A;            /* linear system matrix */
  KSP            ksp;         /* KSP context */
  PetscErrorCode ierr;
  PetscInt       i,n = 10,col[3],its,i1,i2;
  PetscScalar    none = -1.0,value[3],avalue;
  PetscReal      norm;
  PC             pc;

  ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
  ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr);

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

  /* create a solution that is orthogonal to the constants */
  ierr = VecGetOwnershipRange(u,&i1,&i2);CHKERRQ(ierr);
  for (i=i1; i<i2; i++) {
    avalue = i;
    VecSetValues(u,1,&i,&avalue,INSERT_VALUES);
  }
  ierr   = VecAssemblyBegin(u);CHKERRQ(ierr);
  ierr   = VecAssemblyEnd(u);CHKERRQ(ierr);
  ierr   = VecSum(u,&avalue);CHKERRQ(ierr);
  avalue = -avalue/(PetscReal)n;
  ierr   = VecShift(u,avalue);CHKERRQ(ierr);

  /* Create and assemble matrix */
  ierr     = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr     = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
  ierr     = MatSetFromOptions(A);CHKERRQ(ierr);
  value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
  for (i=1; i<n-1; i++) {
    col[0] = i-1; col[1] = i; col[2] = i+1;
    ierr   = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
  }
  i    = n - 1; col[0] = n - 2; col[1] = n - 1; value[1] = 1.0;
  ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
  i    = 0; col[0] = 0; col[1] = 1; value[0] = 1.0; value[1] = -1.0;
  ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatMult(A,u,b);CHKERRQ(ierr);

  /* Create KSP context; set operators and options; solve linear system */
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);

  /* Insure that preconditioner has same null space as matrix */
  /* currently does not do anything */
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);

  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
  /* ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */

  /* Check error */
  ierr = VecAXPY(x,none,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 %g, Iterations %D\n",(double)norm,its);CHKERRQ(ierr);

  /* Free work space */
  ierr = VecDestroy(&x);CHKERRQ(ierr);ierr = VecDestroy(&u);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr);ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return ierr;
}

int main(int argc,char **args)
{
  Mat              *A,B;           /* matrix */
  PetscErrorCode   ierr;
  Vec              x,y,v,v2,z;
  PetscReal        rnorm;
  PetscInt         n = 20;         /* size of the matrix */
  PetscInt         nmat = 3;       /* number of matrices */
  PetscInt         i;
  PetscRandom      rctx;
  MatCompositeType type;

  ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
  ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(NULL,NULL,"-nmat",&nmat,NULL);CHKERRQ(ierr);

  /*
     Create random matrices
  */
  ierr = PetscMalloc1(nmat+3,&A);CHKERRQ(ierr);
  ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
  ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n/2,3,NULL,3,NULL,&A[0]);CHKERRQ(ierr);
  for (i = 1; i < nmat+1; i++) {
    ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,3,NULL,3,NULL,&A[i]);CHKERRQ(ierr);
  }
  ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n/2,n,3,NULL,3,NULL,&A[nmat+1]);CHKERRQ(ierr);
  for (i = 0; i < nmat+2; i++) {
    ierr = MatSetRandom(A[i],rctx);CHKERRQ(ierr);
  }

  ierr = MatCreateVecs(A[1],&x,&y);CHKERRQ(ierr);
  ierr = VecDuplicate(y,&z);CHKERRQ(ierr);
  ierr = MatCreateVecs(A[0],&v,NULL);CHKERRQ(ierr);
  ierr = VecDuplicate(v,&v2);CHKERRQ(ierr);

  ierr = VecSet(x,1.0);CHKERRQ(ierr);
  ierr = VecSet(y,0.0);CHKERRQ(ierr);
  ierr = MatMult(A[1],x,z);CHKERRQ(ierr);
  for (i = 2; i < nmat+1; i++) {
    ierr = MatMultAdd(A[i],x,z,z);CHKERRQ(ierr);
  }

  ierr = MatCreateComposite(PETSC_COMM_WORLD,nmat,A+1,&B);CHKERRQ(ierr);
  ierr = MatMultAdd(B,x,y,y);CHKERRQ(ierr);
  ierr = VecAXPY(y,-1.0,z);CHKERRQ(ierr);
  ierr = VecNorm(y,NORM_2,&rnorm);CHKERRQ(ierr);
  if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with composite add %g\n",(double)rnorm);CHKERRQ(ierr);
  }

  ierr = MatCompositeSetMatStructure(B,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); /* default */
  ierr = MatCompositeMerge(B);CHKERRQ(ierr);
  ierr = MatMult(B,x,y);CHKERRQ(ierr);
  ierr = MatDestroy(&B);CHKERRQ(ierr);
  ierr = VecAXPY(y,-1.0,z);CHKERRQ(ierr);
  ierr = VecNorm(y,NORM_2,&rnorm);CHKERRQ(ierr);
  if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with composite add after merge %g\n",(double)rnorm);CHKERRQ(ierr);
  }

  /*
     Test n x n/2 multiplicative composite
  */
  ierr = VecSet(v,1.0);CHKERRQ(ierr);
  ierr = MatMult(A[0],v,z);CHKERRQ(ierr);
  for (i = 1; i < nmat; i++) {
    ierr = MatMult(A[i],z,y);CHKERRQ(ierr);
    ierr = VecCopy(y,z);CHKERRQ(ierr);
  }

  ierr = MatCreateComposite(PETSC_COMM_WORLD,nmat,A,&B);CHKERRQ(ierr);
  ierr = MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr);
  ierr = MatCompositeSetMergeType(B,MAT_COMPOSITE_MERGE_LEFT);CHKERRQ(ierr);
  ierr = MatSetFromOptions(B);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); /* do MatCompositeMerge() if -mat_composite_merge 1 */
  ierr = MatMult(B,v,y);CHKERRQ(ierr);
  ierr = MatDestroy(&B);CHKERRQ(ierr);
  ierr = VecAXPY(y,-1.0,z);CHKERRQ(ierr);
  ierr = VecNorm(y,NORM_2,&rnorm);CHKERRQ(ierr);
  if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative %g\n",(double)rnorm);CHKERRQ(ierr);
  }

  /*
     Test n/2 x n multiplicative composite
  */
  ierr = VecSet(x,1.0);CHKERRQ(ierr);
  ierr = MatMult(A[2],x,z);CHKERRQ(ierr);
  for (i = 3; i < nmat+1; i++) {
    ierr = MatMult(A[i],z,y);CHKERRQ(ierr);
    ierr = VecCopy(y,z);CHKERRQ(ierr);
  }
  ierr = MatMult(A[nmat+1],z,v);CHKERRQ(ierr);

  ierr = MatCreateComposite(PETSC_COMM_WORLD,nmat,A+2,&B);CHKERRQ(ierr);
  ierr = MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr);
  ierr = MatSetFromOptions(B);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); /* do MatCompositeMerge() if -mat_composite_merge 1 */
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
  PetscErrorCode ierr;
  Mat            A,B,Asub,Bsub;
  PetscInt       ms,idxrow[3],idxcol[4];
  Vec            left,right,X,Y,X1,Y1;
  IS             isrow,iscol;
  PetscBool      random = PETSC_TRUE;

  ierr = PetscInitialize(&argc,&argv,PETSC_NULL,help);CHKERRQ(ierr);
  ierr = AssembleMatrix(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = AssembleMatrix(PETSC_COMM_WORLD,&B);CHKERRQ(ierr);
  ierr = MatShellSetOperation(B,MATOP_GET_SUBMATRIX,PETSC_NULL);CHKERRQ(ierr);
  ierr = MatShellSetOperation(B,MATOP_GET_SUBMATRICES,PETSC_NULL);CHKERRQ(ierr);
  ierr = MatGetOwnershipRange(A,&ms,PETSC_NULL);CHKERRQ(ierr);

  idxrow[0] = ms+1;
  idxrow[1] = ms+2;
  idxrow[2] = ms+4;
  ierr = ISCreateGeneral(PETSC_COMM_WORLD,3,idxrow,PETSC_USE_POINTER,&isrow);CHKERRQ(ierr);

  idxcol[0] = ms+1;
  idxcol[1] = ms+2;
  idxcol[2] = ms+4;
  idxcol[3] = ms+5;
  ierr = ISCreateGeneral(PETSC_COMM_WORLD,4,idxcol,PETSC_USE_POINTER,&iscol);CHKERRQ(ierr);

  ierr = MatGetSubMatrix(A,isrow,iscol,MAT_INITIAL_MATRIX,&Asub);CHKERRQ(ierr);
  ierr = MatGetSubMatrix(B,isrow,iscol,MAT_INITIAL_MATRIX,&Bsub);CHKERRQ(ierr);

  ierr = MatGetVecs(Asub,&right,&left);CHKERRQ(ierr);
  ierr = VecDuplicate(right,&X);CHKERRQ(ierr);
  ierr = VecDuplicate(right,&X1);CHKERRQ(ierr);
  ierr = VecDuplicate(left,&Y);CHKERRQ(ierr);
  ierr = VecDuplicate(left,&Y1);CHKERRQ(ierr);

  ierr = PetscOptionsGetBool(PETSC_NULL,"-random",&random,PETSC_NULL);CHKERRQ(ierr);
  if (random) {
    ierr = VecSetRandom(right,PETSC_NULL);CHKERRQ(ierr);
    ierr = VecSetRandom(left,PETSC_NULL);CHKERRQ(ierr);
    ierr = VecSetRandom(X,PETSC_NULL);CHKERRQ(ierr);
    ierr = VecSetRandom(Y,PETSC_NULL);CHKERRQ(ierr);
    ierr = VecSetRandom(X1,PETSC_NULL);CHKERRQ(ierr);
    ierr = VecSetRandom(Y1,PETSC_NULL);CHKERRQ(ierr);
  } else {
    ierr = VecSet(right,1.0);CHKERRQ(ierr);
    ierr = VecSet(left,2.0);CHKERRQ(ierr);
    ierr = VecSet(X,3.0);CHKERRQ(ierr);
    ierr = VecSet(Y,4.0);CHKERRQ(ierr);
    ierr = VecSet(X1,3.0);CHKERRQ(ierr);
    ierr = VecSet(Y1,4.0);CHKERRQ(ierr);
  }
  ierr = CheckMatrices(Asub,Bsub,left,right,X,Y,X1,Y1);CHKERRQ(ierr);
  ierr = ISDestroy(&isrow);CHKERRQ(ierr);
  ierr = ISDestroy(&iscol);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = MatDestroy(&B);CHKERRQ(ierr);
  ierr = MatDestroy(&Asub);CHKERRQ(ierr);
  ierr = MatDestroy(&Bsub);CHKERRQ(ierr);
  ierr = VecDestroy(&left);CHKERRQ(ierr);
  ierr = VecDestroy(&right);CHKERRQ(ierr);
  ierr = VecDestroy(&X);CHKERRQ(ierr);
  ierr = VecDestroy(&Y);CHKERRQ(ierr);
  ierr = VecDestroy(&X1);CHKERRQ(ierr);
  ierr = VecDestroy(&Y1);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}

int main(int argc, char **argv)
{
  PetscErrorCode ierr;
  PetscInt       n=10000,its,dfid=1;
  Vec            x,b,u;
  Mat            A;
  KSP            ksp;
  PC             pc,pcnoise;
  PCNoise_Ctx    ctx={0,NULL};
  PetscReal      eta=0.1,norm;
  PetscScalar(*diagfunc)(PetscInt,PetscInt);

  ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
  /* Process command line options */
  ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(NULL,NULL,"-eta",&eta,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(NULL,NULL,"-diagfunc",&dfid,NULL);CHKERRQ(ierr);
  switch(dfid){
    case 1:
      diagfunc = diagFunc1;
      break;
    case 2:
      diagfunc = diagFunc2;
      break;
    case 3:
      diagfunc = diagFunc3;
      break;
    default:
      SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unrecognized diagfunc option");
  }

  /* Create a diagonal matrix with a given distribution of diagonal elements */
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatSetUp(A);CHKERRQ(ierr);
  ierr = AssembleDiagonalMatrix(A,diagfunc);CHKERRQ(ierr);

  /* Allocate vectors and manufacture an exact solution and rhs */
  ierr = MatCreateVecs(A,&x,NULL);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject)x,"Computed Solution");CHKERRQ(ierr);
  ierr = MatCreateVecs(A,&b,NULL);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject)b,"RHS");CHKERRQ(ierr);
  ierr = MatCreateVecs(A,&u,NULL);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject)u,"Reference Solution");CHKERRQ(ierr);
  ierr = VecSet(u,1.0);CHKERRQ(ierr);
  ierr = MatMult(A,u,b);CHKERRQ(ierr);

  /* Create a KSP object */
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);

  /* Set up a composite preconditioner */
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
  ierr = PCSetType(pc,PCCOMPOSITE);CHKERRQ(ierr); /* default composite with single Identity PC */
  ierr = PCCompositeSetType(pc,PC_COMPOSITE_ADDITIVE);CHKERRQ(ierr);
  ierr = PCCompositeAddPC(pc,PCNONE);CHKERRQ(ierr);
  if(eta > 0){
    ierr = PCCompositeAddPC(pc,PCSHELL);CHKERRQ(ierr);
    ierr = PCCompositeGetPC(pc,1,&pcnoise);CHKERRQ(ierr);
    ctx.eta = eta;
    ierr = PCShellSetContext(pcnoise,&ctx);CHKERRQ(ierr);
    ierr = PCShellSetApply(pcnoise,PCApply_Noise);CHKERRQ(ierr);
    ierr = PCShellSetSetUp(pcnoise,PCSetup_Noise);CHKERRQ(ierr);
    ierr = PCShellSetDestroy(pcnoise,PCDestroy_Noise);CHKERRQ(ierr);
    ierr = PCShellSetName(pcnoise,"Noise PC");CHKERRQ(ierr);
  }

  /* Set KSP from options (this can override the PC just defined) */
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

  /* Solve */
  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);

  /* Compute error */
  ierr = VecAXPY(x,-1.0,u);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject)x,"Error");CHKERRQ(ierr);
  ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
  ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);CHKERRQ(ierr);

  /* Destroy objects and finalize */
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr);
  ierr = VecDestroy(&u);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return ierr;
}

int main(int argc,char **argv)
{
  SNES                snes;    /* nonlinear solver context */
  Vec                 x,r;     /* solution, residual vectors */
  Mat                 J;       /* Jacobian matrix */
  PetscErrorCode      ierr;
  PetscInt            its;
  PetscScalar         *xx;
  SNESConvergedReason reason;

  PetscInitialize(&argc,&argv,(char*)0,help);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Create nonlinear solver context
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Create matrix and vector data structures; set corresponding routines
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /*
     Create vectors for solution and nonlinear function
  */
  ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
  ierr = VecSetSizes(x,PETSC_DECIDE,2);CHKERRQ(ierr);
  ierr = VecSetFromOptions(x);CHKERRQ(ierr);
  ierr = VecDuplicate(x,&r);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);

  /*
     Set function evaluation routine and vector.
  */
  ierr = SNESSetFunction(snes,r,FormFunction1,NULL);CHKERRQ(ierr);

  /*
     Set Jacobian matrix data structure and Jacobian evaluation routine
  */
  ierr = SNESSetJacobian(snes,J,J,FormJacobian1,NULL);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Customize nonlinear solver; set runtime options
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Evaluate initial guess; then solve nonlinear system
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr  = VecGetArray(x,&xx);CHKERRQ(ierr);
  xx[0] = -1.2; xx[1] = 1.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,NULL,x);CHKERRQ(ierr);
  ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
  ierr = SNESGetConvergedReason(snes,&reason);CHKERRQ(ierr);
  /*
     Some systems computes a residual that is identically zero, thus converging
     due to FNORM_ABS, others converge due to FNORM_RELATIVE.  Here, we only
     report the reason if the iteration did not converge so that the tests are
     reproducible.
  */
  ierr = PetscPrintf(PETSC_COMM_WORLD,"%s number of SNES iterations = %D\n\n",reason>0 ? "CONVERGED" : SNESConvergedReasons[reason],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);
  ierr = PetscFinalize();
  return 0;
}