int
SparseGp_logLikeGrad (SparseGp *gp, HyperParam hp, int lengthInd, double *logLikeGrad)
{
    PetscErrorCode ierr;
    (void) ierr;
    /* compute t' inv(K) (dKdt inv(K) t) */
    /* 1. solve inv(K) t */
    PetscInt N = gp->trainLabels->size;
    Vec invKt;
    ierr = petsc_util_createVec (&invKt, gp->nlocal, N);
    SparseGp_solve (gp, gp->_trainLabels, &invKt);

    /* 2. multiply dKdt by invKt */
    Vec dKdtInvKt;
    ierr = petsc_util_createVec (&dKdtInvKt, gp->nlocal, N);
    SparseGp_KGradient (gp, hp, lengthInd);
    MatMult (gp->_KGradient, invKt, dKdtInvKt);

    /* 3. solve invK (vector from step 2.) */
    Vec invKDkdtInvKt;
    ierr = petsc_util_createVec (&invKDkdtInvKt, gp->nlocal, N);
    SparseGp_solve (gp, dKdtInvKt, &invKDkdtInvKt);

    /* 4. compute inner product */
    double dotProd;
    VecDot (gp->_trainLabels, invKDkdtInvKt, &dotProd);

    /* compute trace (invK dKdt) */
    /* 1. for each column in dKdt, solve */
    double myTrace = 0;

    Vec col, solution;
    petsc_util_createVec (&col, gp->nlocal, N);
    petsc_util_createVec (&solution, gp->nlocal, N);

    for (int i=0; i<N; i++) {
        /* IOU_ROOT_PRINT ("%d --- %d\n", i, N); */
        int row = i;

        ierr = MatGetColumnVector (gp->_KGradient, col, row);

        /* IOU_ROOT_PRINT ("solving...\n"); */
        SparseGp_solve (gp, col, &solution);
        /* IOU_ROOT_PRINT ("solved\n"); */

        if (row < gp->rend && row >= gp->rstart) {
            double ii;
            VecGetValues (solution, 1, &row, &ii);
            myTrace += ii;
        }

    }

    /* gather all trace terms */
    int numProcs;
    MPI_Comm_size (PETSC_COMM_WORLD, &numProcs);
    double diag[numProcs];

    int ret = MPI_Allgather (&myTrace, 1, MPI_DOUBLE, diag, 1, MPI_DOUBLE, PETSC_COMM_WORLD);
    (void) ret;

    double trace = 0;
    for (int i=0; i<numProcs; i++)
        trace += diag[i];

    /* IOU_ROOT_PRINT ("Cleaning up...\n"); */
    /* clean up */
    VecDestroy (&invKt);
    VecDestroy (&dKdtInvKt);
    VecDestroy (&invKDkdtInvKt);
    VecDestroy (&col);
    VecDestroy (&solution);
    /* IOU_ROOT_PRINT ("Done cleaning up...\n"); */

    *logLikeGrad = 0.5*dotProd + 0.5*trace;

    return EXIT_SUCCESS;
}

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  KSP            ksp;
  PC             pc;
  Vec            x,b;
  DM             da;
  Mat            A,Atrans;
  PetscInt       dof=1,M=-8;
  PetscBool      flg,trans=PETSC_FALSE;

  PetscInitialize(&argc,&argv,(char *)0,help);
  ierr = PetscOptionsGetInt(PETSC_NULL,"-dof",&dof,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetBool(PETSC_NULL,"-trans",&trans,PETSC_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,M,M,M);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,PETSC_NULL,PETSC_NULL,PETSC_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 = ComputeRHS(da,b);CHKERRQ(ierr);
  ierr = DMCreateMatrix(da,MATBAIJ,&A);CHKERRQ(ierr);
  ierr = ComputeMatrix(da,A);CHKERRQ(ierr);


  /* A is non-symmetric. Make A = 0.5*(A + Atrans) symmetric for testing icc and cholesky */
  ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&Atrans);CHKERRQ(ierr);
  ierr = MatAXPY(A,1.0,Atrans,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = MatScale(A,0.5);CHKERRQ(ierr);
  ierr = MatDestroy(&Atrans);CHKERRQ(ierr);

  /* Test sbaij matrix */
  flg  = PETSC_FALSE;
  ierr = PetscOptionsGetBool(PETSC_NULL, "-test_sbaij1", &flg,PETSC_NULL);CHKERRQ(ierr);
  if (flg){
    Mat sA;
    PetscBool issymm;
    ierr = MatIsTranspose(A,A,0.0,&issymm);CHKERRQ(ierr);
    if (issymm) {
      ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
    } else {
      printf("Warning: A is non-symmetric\n");
    }
    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);
 
  if (trans) {
    ierr = KSPSolveTranspose(ksp,b,x);CHKERRQ(ierr);
  } else {
    ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
  }

  /* check final residual */
  flg  = PETSC_FALSE;
  ierr = PetscOptionsGetBool(PETSC_NULL, "-check_final_residual", &flg,PETSC_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;
}

//.........这里部分代码省略.........
    ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr);
  } 
  else if (function == CONSTANT) {
    ierr = VecSet(x, 1.0);CHKERRQ(ierr);
  } 
  else if (function == TANH) {
    PetscScalar *a;
    ierr = VecGetArray(x, &a);CHKERRQ(ierr);
    PetscInt j,k = 0;
    for(i = 0; i < 3; ++i) {
      for(j = 0; j < dim[i]; ++j) {
        a[k] = tanh((j - dim[i]/2.0)*(10.0/dim[i]));
        ++k;
      }
    }
    ierr = VecRestoreArray(x, &a);CHKERRQ(ierr);
  }
  if(view_x) {
    ierr = VecView(x, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  }
  ierr = VecCopy(x,xx);CHKERRQ(ierr);
  // Split xx
  ierr = VecStrideGatherAll(xx,xxsplit, INSERT_VALUES);CHKERRQ(ierr); //YES! 'Gather' means 'split' (or maybe 'scatter'?)! 

  ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_SELF, "|x|_2 = %g\n",norm);CHKERRQ(ierr);
  
  /* create USFFT object */
  ierr = MatCreateSeqUSFFT(da,da,&A);CHKERRQ(ierr);
  /* create FFTW object */
  ierr = MatCreateSeqFFTW(PETSC_COMM_SELF,3,dim,&AA);CHKERRQ(ierr);
  
  /* apply USFFT and FFTW FORWARD "preemptively", so the fftw_plans can be reused on different vectors */
  ierr = MatMult(A,x,z);CHKERRQ(ierr);
  for(int ii = 0; ii < 3; ++ii) {
    ierr = MatMult(AA,xxsplit[ii],zzsplit[ii]);CHKERRQ(ierr);
  }
  // Now apply USFFT and FFTW forward several (3) times
  for (i=0; i<3; ++i){
    ierr = MatMult(A,x,y);CHKERRQ(ierr); 
    for(int ii = 0; ii < 3; ++ii) {
      ierr = MatMult(AA,xxsplit[ii],yysplit[ii]);CHKERRQ(ierr);
    }
    ierr = MatMultTranspose(A,y,z);CHKERRQ(ierr);
    for(int ii = 0; ii < 3; ++ii) {
      ierr = MatMult(AA,yysplit[ii],zzsplit[ii]);CHKERRQ(ierr);
    }
  }
  // Unsplit yy
  ierr = VecStrideScatterAll(yysplit, yy, INSERT_VALUES);CHKERRQ(ierr); //YES! 'Scatter' means 'collect' (or maybe 'gather'?)! 
  // Unsplit zz
  ierr = VecStrideScatterAll(zzsplit, zz, INSERT_VALUES);CHKERRQ(ierr); //YES! 'Scatter' means 'collect' (or maybe 'gather'?)! 

  if(view_y) {
    ierr = PetscPrintf(PETSC_COMM_WORLD, "y = \n");CHKERRQ(ierr);
    ierr = VecView(y, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "yy = \n");CHKERRQ(ierr);
    ierr = VecView(yy, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  }
  
  if(view_z) {
    ierr = PetscPrintf(PETSC_COMM_WORLD, "z = \n");CHKERRQ(ierr);
    ierr = VecView(z, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "zz = \n");CHKERRQ(ierr);
    ierr = VecView(zz, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  }

void bsscr_summary(KSP_BSSCR * bsscrp_self, KSP ksp_S, KSP ksp_inner,
		   Mat K,Mat K2,Mat D,Mat G,Mat C,Vec u,Vec p,Vec f,Vec h,Vec t,
		   double penaltyNumber,PetscTruth KisJustK,double mgSetupTime,double scrSolveTime,double a11SingleSolveTime){
    PetscTruth flg, found;
    PetscInt   uSize, pSize, lmax, lmin, iterations;
    PetscReal  rNorm, fNorm, uNorm, uNormInf, pNorm, pNormInf, p_sum, min, max;
    Vec q, qq, t2, t3;
    double solutionAnalysisTime;

      PetscPrintf( PETSC_COMM_WORLD,  "\n\nSCR Solver Summary:\n\n");
      if(bsscrp_self->mg)
	    PetscPrintf( PETSC_COMM_WORLD, "  Multigrid setup:        = %.4g secs \n", mgSetupTime);

      KSPGetIterationNumber( ksp_S, &iterations);
      bsscrp_self->solver->stats.pressure_its = iterations;
      PetscPrintf( PETSC_COMM_WORLD,     "  Pressure Solve:         = %.4g secs / %d its\n", scrSolveTime, iterations);
      KSPGetIterationNumber( ksp_inner, &iterations);
      bsscrp_self->solver->stats.velocity_backsolve_its = iterations;
      PetscPrintf( PETSC_COMM_WORLD,     "  Final V Solve:          = %.4g secs / %d its\n\n", a11SingleSolveTime, iterations);

      /***************************************************************************************************************/

      flg = PETSC_FALSE; /* Off by default */
	  PetscOptionsGetTruth( PETSC_NULL, "-scr_ksp_solution_summary", &flg, &found );

      if(flg) {
	    PetscScalar KuNorm;
	    solutionAnalysisTime = MPI_Wtime();
	    VecGetSize( u, &uSize );
	    VecGetSize( p, &pSize );
	    VecDuplicate( u, &t2 );
	    VecDuplicate( u, &t3 );
	    MatMult( K, u, t3); VecNorm( t3, NORM_2, &KuNorm );
	    double angle, kdot;
	    if(penaltyNumber > 1e-10){/* should change this to ifK2built maybe */
          MatMult( K2, u, t2); VecNorm( t2, NORM_2, &rNorm );
          VecDot(t2,t3,&kdot);
          angle = (kdot/(rNorm*KuNorm));
          PetscPrintf( PETSC_COMM_WORLD,  "  <K u, K2 u>/(|K u| |K2 u|)    = %.6e\n", angle);
	    }
	    VecNorm( t, NORM_2, &rNorm ); /* t = f- G p  should be the formal residual vector, calculated on line 267 in auglag-driver-DGTGD.c */
	    VecDot(t3,t,&kdot);
	    angle = (kdot/(rNorm*KuNorm));
            PetscPrintf( PETSC_COMM_WORLD,  "  <K u, (f-G p)>/(|K u| |f- G p|)    = %.6e\n\n", angle);

	    MatMult( K, u, t2); VecNorm(t2, NORM_2, &KuNorm);
	    VecAYPX( t2, -1.0, t ); /* t2 <- -t2 + t  : t = f- G p  should be the formal residual vector, calculated on line 267 in auglag-driver-DGTGD.c*/
	    VecNorm( t2, NORM_2, &rNorm );
	    VecNorm( f,  NORM_2, &fNorm );
	    if(KisJustK){
          PetscPrintf( PETSC_COMM_WORLD,"Velocity back-solve with original K matrix\n");
          PetscPrintf( PETSC_COMM_WORLD,"Solved    K u = G p -f\n");
          PetscPrintf( PETSC_COMM_WORLD,"Residual with original K matrix\n");
          PetscPrintf( PETSC_COMM_WORLD,  "  |f - K u - G p|                       = %.12e\n", rNorm);
          PetscPrintf( PETSC_COMM_WORLD,  "  |f - K u - G p|/|f|                   = %.12e\n", rNorm/fNorm);
          if(penaltyNumber > 1e-10){/* means the restore_K flag was used */
            //if(K2 && f2){
            MatAXPY(K,penaltyNumber,K2,DIFFERENT_NONZERO_PATTERN);/* Computes K = penaltyNumber*K2 + K */
            //VecAXPY(f,penaltyNumber,f2); /* f = penaltyNumber*f2 + f */
            KisJustK=PETSC_FALSE;
            MatMult( K, u, t2);
            MatMult( G, p, t);
            VecAYPX( t, -1.0, f ); /* t <- -t + f */
            VecAYPX( t2, -1.0, t ); /* t2 <- -t2 + t */
            VecNorm( t2, NORM_2, &rNorm );
            PetscPrintf( PETSC_COMM_WORLD,"Residual with K+K2 matrix and f rhs vector\n");
            PetscPrintf( PETSC_COMM_WORLD,  "  |(f) - (K + K2) u - G p|         = %.12e\n", rNorm);
            //}
          }
	    }
	    else{
          PetscPrintf( PETSC_COMM_WORLD,"Velocity back-solve with K+K2 matrix\n");
          PetscPrintf( PETSC_COMM_WORLD,"Solved    (K + K2) u = G p - (f)\n");
          PetscPrintf( PETSC_COMM_WORLD,"Residual with K+K2 matrix and f rhs vector\n");
          PetscPrintf( PETSC_COMM_WORLD,  "  |(f) - (K + K2) u - G p|         = %.12e\n", rNorm);
          PetscReal KK2Norm,KK2Normf;
          MatNorm(K,NORM_1,&KK2Norm);
          MatNorm(K,NORM_FROBENIUS,&KK2Normf);
          penaltyNumber = -penaltyNumber;
          MatAXPY(K,penaltyNumber,K2,DIFFERENT_NONZERO_PATTERN);/* Computes K = penaltyNumber*K2 + K */
          //VecAXPY(f,penaltyNumber,f2); /* f = penaltyNumber*f2 + f */
          KisJustK=PETSC_FALSE;
          MatMult( K, u, t2);    /* t2 = K*u  */
          MatMult( G, p, t);     /* t  = G*p  */
          VecAYPX( t, -1.0, f ); /* t <- f - t ; t = f - G*p  */
          VecAYPX( t2, -1.0, t ); /* t2 <- t - t2; t2 = f - G*p - K*u  */
          VecNorm( t2, NORM_2, &rNorm );
          PetscPrintf( PETSC_COMM_WORLD,"Residual with original K matrix\n");
          PetscPrintf( PETSC_COMM_WORLD,  "  |f - K u - G p|                       = %.12e\n", rNorm);
          PetscPrintf( PETSC_COMM_WORLD,  "  |f - K u - G p|/|f|                   = %.12e\n", rNorm/fNorm);
          PetscReal KNorm, K2Norm;
          MatNorm(K,NORM_1,&KNorm);	  MatNorm(K2,NORM_1,&K2Norm);
          PetscPrintf( PETSC_COMM_WORLD,"K and K2 norm_1    %.12e %.12e   ratio %.12e\n",KNorm,K2Norm,K2Norm/KNorm);
          MatNorm(K,NORM_INFINITY,&KNorm);  MatNorm(K2,NORM_INFINITY,&K2Norm);
          PetscPrintf( PETSC_COMM_WORLD,"K and K2 norm_inf  %.12e %.12e   ratio %.12e\n",KNorm,K2Norm,K2Norm/KNorm);
          MatNorm(K,NORM_FROBENIUS,&KNorm); MatNorm(K2,NORM_FROBENIUS,&K2Norm);
          PetscPrintf( PETSC_COMM_WORLD,"K and K2 norm_frob %.12e %.12e   ratio %.12e\n",KNorm,K2Norm,K2Norm/KNorm);
          PetscPrintf( PETSC_COMM_WORLD,"K+r*K2 norm_1    %.12e\n",KK2Norm);
          PetscPrintf( PETSC_COMM_WORLD,"K+r*K2 norm_frob %.12e\n",KK2Normf);
          penaltyNumber = -penaltyNumber;
//.........这里部分代码省略.........

int main(int argc,char **args)
{
    Mat            mat;          /* matrix */
    Vec            b,ustar,u;  /* vectors (RHS, exact solution, approx solution) */
    PC             pc;           /* PC context */
    KSP            ksp;          /* KSP context */
    PetscErrorCode ierr;
    PetscInt       n = 10,i,its,col[3];
    PetscScalar    value[3];
    PCType         pcname;
    KSPType        kspname;
    PetscReal      norm,tol=1.e-14;

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

    /* Create and initialize vectors */
    ierr = VecCreateSeq(PETSC_COMM_SELF,n,&b);
    CHKERRQ(ierr);
    ierr = VecCreateSeq(PETSC_COMM_SELF,n,&ustar);
    CHKERRQ(ierr);
    ierr = VecCreateSeq(PETSC_COMM_SELF,n,&u);
    CHKERRQ(ierr);
    ierr = VecSet(ustar,1.0);
    CHKERRQ(ierr);
    ierr = VecSet(u,0.0);
    CHKERRQ(ierr);

    /* Create and assemble matrix */
    ierr     = MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,&mat);
    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(mat,1,&i,3,col,value,INSERT_VALUES);
        CHKERRQ(ierr);
    }
    i    = n - 1;
    col[0] = n - 2;
    col[1] = n - 1;
    ierr = MatSetValues(mat,1,&i,2,col,value,INSERT_VALUES);
    CHKERRQ(ierr);
    i    = 0;
    col[0] = 0;
    col[1] = 1;
    value[0] = 2.0;
    value[1] = -1.0;
    ierr = MatSetValues(mat,1,&i,2,col,value,INSERT_VALUES);
    CHKERRQ(ierr);
    ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
    CHKERRQ(ierr);
    ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
    CHKERRQ(ierr);

    /* Compute right-hand-side vector */
    ierr = MatMult(mat,ustar,b);
    CHKERRQ(ierr);

    /* Create PC context and set up data structures */
    ierr = PCCreate(PETSC_COMM_WORLD,&pc);
    CHKERRQ(ierr);
    ierr = PCSetType(pc,PCNONE);
    CHKERRQ(ierr);
    ierr = PCSetFromOptions(pc);
    CHKERRQ(ierr);
    ierr = PCSetOperators(pc,mat,mat);
    CHKERRQ(ierr);
    ierr = PCSetUp(pc);
    CHKERRQ(ierr);

    /* Create KSP context and set up data structures */
    ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);
    CHKERRQ(ierr);
    ierr = KSPSetType(ksp,KSPRICHARDSON);
    CHKERRQ(ierr);
    ierr = KSPSetFromOptions(ksp);
    CHKERRQ(ierr);
    ierr = PCSetOperators(pc,mat,mat);
    CHKERRQ(ierr);
    ierr = KSPSetPC(ksp,pc);
    CHKERRQ(ierr);
    ierr = KSPSetUp(ksp);
    CHKERRQ(ierr);

    /* Solve the problem */
    ierr = KSPGetType(ksp,&kspname);
    CHKERRQ(ierr);
    ierr = PCGetType(pc,&pcname);
    CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_SELF,"Running %s with %s preconditioning\n",kspname,pcname);
    CHKERRQ(ierr);
    ierr = KSPSolve(ksp,b,u);
    CHKERRQ(ierr);
    ierr = VecAXPY(u,-1.0,ustar);
    CHKERRQ(ierr);
    ierr = VecNorm(u,NORM_2,&norm);
    ierr = KSPGetIterationNumber(ksp,&its);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  ierr = ISDestroy(&is_schur);CHKERRQ(ierr);
  if (use_lu) {
    ierr = MatLUFactorSymbolic(F,A,NULL,NULL,NULL);CHKERRQ(ierr);
  } else {
    ierr = MatCholeskyFactorSymbolic(F,A,NULL,NULL);CHKERRQ(ierr);
  }

  for (nfact = 0; nfact < 3; nfact++) {
    Mat AD;

    if (!nfact) {
      ierr = VecSetRandom(x,rand);CHKERRQ(ierr);
      if (symm && herm) {
        ierr = VecAbs(x);CHKERRQ(ierr);
      }
      ierr = MatDiagonalSet(A,x,ADD_VALUES);CHKERRQ(ierr);
    }
    if (use_lu) {
      ierr = MatLUFactorNumeric(F,A,NULL);CHKERRQ(ierr);
    } else {
      ierr = MatCholeskyFactorNumeric(F,A,NULL);CHKERRQ(ierr);
    }
    ierr = MatFactorCreateSchurComplement(F,&S);CHKERRQ(ierr);
    ierr = MatCreateVecs(S,&xschur,&bschur);CHKERRQ(ierr);
    ierr = VecDuplicate(xschur,&uschur);CHKERRQ(ierr);
    if (nfact == 1) {
      ierr = MatFactorInvertSchurComplement(F);CHKERRQ(ierr);
    }
    for (nsolve = 0; nsolve < 2; nsolve++) {
      ierr = VecSetRandom(x,rand);CHKERRQ(ierr);
      ierr = VecCopy(x,u);CHKERRQ(ierr);

      if (nsolve) {
        ierr = MatMult(A,x,b);CHKERRQ(ierr);
        ierr = MatSolve(F,b,x);CHKERRQ(ierr);
      } else {
        ierr = MatMultTranspose(A,x,b);CHKERRQ(ierr);
        ierr = MatSolveTranspose(F,b,x);CHKERRQ(ierr);
      }
      /* Check the error */
      ierr = VecAXPY(u,-1.0,x);CHKERRQ(ierr);  /* u <- (-1.0)x + u */
      ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr);
      if (norm > tol) {
        PetscReal resi;
        if (nsolve) {
          ierr = MatMult(A,x,u);CHKERRQ(ierr); /* u = A*x */
        } else {
          ierr = MatMultTranspose(A,x,u);CHKERRQ(ierr); /* u = A*x */
        }
        ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr);  /* u <- (-1.0)b + u */
        ierr = VecNorm(u,NORM_2,&resi);CHKERRQ(ierr);
        if (nsolve) {
          ierr = PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatSolve error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);CHKERRQ(ierr);
        } else {
          ierr = PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatSolveTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);CHKERRQ(ierr);
        }
      }
      ierr = VecSetRandom(xschur,rand);CHKERRQ(ierr);
      ierr = VecCopy(xschur,uschur);CHKERRQ(ierr);
      if (nsolve) {
        ierr = MatMult(S,xschur,bschur);CHKERRQ(ierr);
        ierr = MatFactorSolveSchurComplement(F,bschur,xschur);CHKERRQ(ierr);
      } else {
        ierr = MatMultTranspose(S,xschur,bschur);CHKERRQ(ierr);
        ierr = MatFactorSolveSchurComplementTranspose(F,bschur,xschur);CHKERRQ(ierr);
      }

PetscErrorCode ResidualFunction(SNES snes,Vec X, Vec F, Userctx *user)
{
  PetscErrorCode ierr;
  Vec            Xgen,Xnet,Fgen,Fnet;
  PetscScalar    *xgen,*xnet,*fgen,*fnet;
  PetscInt       i,idx=0;
  PetscScalar    Vr,Vi,Vm,Vm2;
  PetscScalar    Eqp,Edp,delta,w; /* Generator variables */
  PetscScalar    Efd,RF,VR; /* Exciter variables */
  PetscScalar    Id,Iq;  /* Generator dq axis currents */
  PetscScalar    Vd,Vq,SE;
  PetscScalar    IGr,IGi,IDr,IDi;
  PetscScalar    Zdq_inv[4],det;
  PetscScalar    PD,QD,Vm0,*v0;
  PetscInt       k;

  PetscFunctionBegin;
  ierr = VecZeroEntries(F);CHKERRQ(ierr);
  ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr);
  ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Fgen,&Fnet);CHKERRQ(ierr);
  ierr = DMCompositeScatter(user->dmpgrid,X,Xgen,Xnet);CHKERRQ(ierr);
  ierr = DMCompositeScatter(user->dmpgrid,F,Fgen,Fnet);CHKERRQ(ierr);

  /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
     The generator current injection, IG, and load current injection, ID are added later
  */
  /* Note that the values in Ybus are stored assuming the imaginary current balance
     equation is ordered first followed by real current balance equation for each bus.
     Thus imaginary current contribution goes in location 2*i, and
     real current contribution in 2*i+1
  */
  ierr = MatMult(user->Ybus,Xnet,Fnet);CHKERRQ(ierr);

  ierr = VecGetArray(Xgen,&xgen);CHKERRQ(ierr);
  ierr = VecGetArray(Xnet,&xnet);CHKERRQ(ierr);
  ierr = VecGetArray(Fgen,&fgen);CHKERRQ(ierr);
  ierr = VecGetArray(Fnet,&fnet);CHKERRQ(ierr);

  /* Generator subsystem */
  for (i=0; i < ngen; i++) {
    Eqp   = xgen[idx];
    Edp   = xgen[idx+1];
    delta = xgen[idx+2];
    w     = xgen[idx+3];
    Id    = xgen[idx+4];
    Iq    = xgen[idx+5];
    Efd   = xgen[idx+6];
    RF    = xgen[idx+7];
    VR    = xgen[idx+8];

    /* Generator differential equations */
    fgen[idx]   = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i];
    fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i];
    fgen[idx+2] = -w + w_s;
    fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i];

    Vr = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */
    Vi = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */

    ierr = ri2dq(Vr,Vi,delta,&Vd,&Vq);CHKERRQ(ierr);
    /* Algebraic equations for stator currents */

    det = Rs[i]*Rs[i] + Xdp[i]*Xqp[i];

    Zdq_inv[0] = Rs[i]/det;
    Zdq_inv[1] = Xqp[i]/det;
    Zdq_inv[2] = -Xdp[i]/det;
    Zdq_inv[3] = Rs[i]/det;

    fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id;
    fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq;

    /* Add generator current injection to network */
    ierr = dq2ri(Id,Iq,delta,&IGr,&IGi);CHKERRQ(ierr);

    fnet[2*gbus[i]]   -= IGi;
    fnet[2*gbus[i]+1] -= IGr;

    Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq);

    SE = k1[i]*PetscExpScalar(k2[i]*Efd);

    /* Exciter differential equations */
    fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i];
    fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i];
    fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i];

    idx = idx + 9;
  }

  ierr = VecGetArray(user->V0,&v0);CHKERRQ(ierr);
  for (i=0; i < nload; i++) {
    Vr  = xnet[2*lbus[i]]; /* Real part of load bus voltage */
    Vi  = xnet[2*lbus[i]+1]; /* Imaginary part of the load bus voltage */
    Vm  = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2 = Vm*Vm;
    Vm0 = PetscSqrtScalar(v0[2*lbus[i]]*v0[2*lbus[i]] + v0[2*lbus[i]+1]*v0[2*lbus[i]+1]);
    PD  = QD = 0.0;
    for (k=0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k]*PD0[i]*PetscPowScalar((Vm/Vm0),ld_betap[k]);
    for (k=0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k]*QD0[i]*PetscPowScalar((Vm/Vm0),ld_betaq[k]);

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

int main(int argc,char **argv)
{
  DM             da;            /* distributed array */
  Vec            x,b,u;         /* approx solution, RHS, exact solution */
  Mat            A;             /* linear system matrix */
  KSP            ksp;           /* linear solver context */
  PetscRandom    rctx;          /* random number generator context */
  PetscReal      norm;          /* norm of solution error */
  PetscInt       i,j,its;
  PetscErrorCode ierr;
  PetscBool      flg = PETSC_FALSE;
  PetscLogStage  stage;
  DMDALocalInfo  info;

  ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
  /*
     Create distributed array to handle parallel distribution.
     The problem size will default to 8 by 7, but this can be
     changed using -da_grid_x M -da_grid_y N
  */
  ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-8,-7,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
         Compute the matrix and right-hand-side vector that define
         the linear system, Ax = b.
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /*
     Create parallel matrix preallocated according to the DMDA, format AIJ by default.
     To use symmetric storage, run with -dm_mat_type sbaij -mat_ignore_lower_triangular
  */
  ierr = DMCreateMatrix(da,&A);CHKERRQ(ierr);

  /*
     Set matrix elements for the 2-D, five-point stencil in parallel.
      - Each processor needs to insert only elements that it owns
        locally (but any non-local elements will be sent to the
        appropriate processor during matrix assembly).
      - Rows and columns are specified by the stencil
      - Entries are normalized for a domain [0,1]x[0,1]
   */
  ierr = PetscLogStageRegister("Assembly", &stage);CHKERRQ(ierr);
  ierr = PetscLogStagePush(stage);CHKERRQ(ierr);
  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
  for (j=info.ys; j<info.ys+info.ym; j++) {
    for (i=info.xs; i<info.xs+info.xm; i++) {
      PetscReal   hx  = 1./info.mx,hy = 1./info.my;
      MatStencil  row = {0},col[5] = {{0}};
      PetscScalar v[5];
      PetscInt    ncols = 0;
      row.j        = j; row.i = i;
      col[ncols].j = j; col[ncols].i = i; v[ncols++] = 2*(hx/hy + hy/hx);
      /* boundaries */
      if (i>0)         {col[ncols].j = j;   col[ncols].i = i-1; v[ncols++] = -hy/hx;}
      if (i<info.mx-1) {col[ncols].j = j;   col[ncols].i = i+1; v[ncols++] = -hy/hx;}
      if (j>0)         {col[ncols].j = j-1; col[ncols].i = i;   v[ncols++] = -hx/hy;}
      if (j<info.my-1) {col[ncols].j = j+1; col[ncols].i = i;   v[ncols++] = -hx/hy;}
      ierr = MatSetValuesStencil(A,1,&row,ncols,col,v,INSERT_VALUES);CHKERRQ(ierr);
    }
  }

  /*
     Assemble matrix, using the 2-step process:
       MatAssemblyBegin(), MatAssemblyEnd()
     Computations can be done while messages are in transition
     by placing code between these two statements.
  */
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = PetscLogStagePop();CHKERRQ(ierr);

  /*
     Create parallel vectors compatible with the DMDA.
  */
  ierr = DMCreateGlobalVector(da,&u);CHKERRQ(ierr);
  ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
  ierr = VecDuplicate(u,&x);CHKERRQ(ierr);

  /*
     Set exact solution; then compute right-hand-side vector.
     By default we use an exact solution of a vector with all
     elements of 1.0;  Alternatively, using the runtime option
     -random_sol forms a solution vector with random components.
  */
  ierr = PetscOptionsGetBool(NULL,NULL,"-random_exact_sol",&flg,NULL);CHKERRQ(ierr);
  if (flg) {
    ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
    ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
    ierr = VecSetRandom(u,rctx);CHKERRQ(ierr);
    ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);
  } else {
    ierr = VecSet(u,1.);CHKERRQ(ierr);
  }
  ierr = MatMult(A,u,b);CHKERRQ(ierr);

  /*
     View the exact solution vector if desired
  */
  flg  = PETSC_FALSE;
  ierr = PetscOptionsGetBool(NULL,NULL,"-view_exact_sol",&flg,NULL);CHKERRQ(ierr);
  if (flg) {ierr = VecView(u,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);}
//.........这里部分代码省略.........

PetscErrorCode KSPFGMRESCycle(PetscInt *itcount,KSP ksp)
{

  KSP_FGMRES     *fgmres = (KSP_FGMRES *)(ksp->data);
  PetscReal      res_norm;
  PetscReal      hapbnd,tt;
  PetscBool      hapend = PETSC_FALSE;  /* indicates happy breakdown ending */
  PetscErrorCode ierr;
  PetscInt       loc_it;                /* local count of # of dir. in Krylov space */
  PetscInt       max_k = fgmres->max_k; /* max # of directions Krylov space */
  Mat            Amat,Pmat;
  MatStructure   pflag;

  PetscFunctionBegin;

  /* Number of pseudo iterations since last restart is the number
     of prestart directions */
  loc_it = 0;

  /* note: (fgmres->it) is always set one less than (loc_it) It is used in
     KSPBUILDSolution_FGMRES, where it is passed to KSPFGMRESBuildSoln.
     Note that when KSPFGMRESBuildSoln is called from this function,
     (loc_it -1) is passed, so the two are equivalent */
  fgmres->it = (loc_it - 1);

  /* initial residual is in VEC_VV(0)  - compute its norm*/
  ierr   = VecNorm(VEC_VV(0),NORM_2,&res_norm);CHKERRQ(ierr);

  /* first entry in right-hand-side of hessenberg system is just
     the initial residual norm */
  *RS(0) = res_norm;

  ksp->rnorm = res_norm;
  KSPLogResidualHistory(ksp,res_norm);

  /* check for the convergence - maybe the current guess is good enough */
  ierr = (*ksp->converged)(ksp,ksp->its,res_norm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  if (ksp->reason) {
    if (itcount) *itcount = 0;
    PetscFunctionReturn(0);
  }

  /* scale VEC_VV (the initial residual) */
  ierr = VecScale(VEC_VV(0),1.0/res_norm);CHKERRQ(ierr);

  /* MAIN ITERATION LOOP BEGINNING*/
  /* keep iterating until we have converged OR generated the max number
     of directions OR reached the max number of iterations for the method */
  while (!ksp->reason && loc_it < max_k && ksp->its < ksp->max_it) {
    if (loc_it) KSPLogResidualHistory(ksp,res_norm);
    fgmres->it = (loc_it - 1);
    ierr = KSPMonitor(ksp,ksp->its,res_norm);CHKERRQ(ierr);

    /* see if more space is needed for work vectors */
    if (fgmres->vv_allocated <= loc_it + VEC_OFFSET + 1) {
      ierr = KSPFGMRESGetNewVectors(ksp,loc_it+1);CHKERRQ(ierr);
      /* (loc_it+1) is passed in as number of the first vector that should
         be allocated */
    }

    /* CHANGE THE PRECONDITIONER? */
    /* ModifyPC is the callback function that can be used to
       change the PC or its attributes before its applied */
    (*fgmres->modifypc)(ksp,ksp->its,loc_it,res_norm,fgmres->modifyctx);


    /* apply PRECONDITIONER to direction vector and store with
       preconditioned vectors in prevec */
    ierr = KSP_PCApply(ksp,VEC_VV(loc_it),PREVEC(loc_it));CHKERRQ(ierr);

    ierr = PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr);
    /* Multiply preconditioned vector by operator - put in VEC_VV(loc_it+1) */
    ierr = MatMult(Amat,PREVEC(loc_it),VEC_VV(1+loc_it));CHKERRQ(ierr);


    /* update hessenberg matrix and do Gram-Schmidt - new direction is in
       VEC_VV(1+loc_it)*/
    ierr = (*fgmres->orthog)(ksp,loc_it);CHKERRQ(ierr);

    /* new entry in hessenburg is the 2-norm of our new direction */
    ierr = VecNorm(VEC_VV(loc_it+1),NORM_2,&tt);CHKERRQ(ierr);
    *HH(loc_it+1,loc_it)   = tt;
    *HES(loc_it+1,loc_it)  = tt;

    /* Happy Breakdown Check */
    hapbnd  = PetscAbsScalar((tt) / *RS(loc_it));
    /* RS(loc_it) contains the res_norm from the last iteration  */
    hapbnd = PetscMin(fgmres->haptol,hapbnd);
    if (tt > hapbnd) {
        /* scale new direction by its norm */
        ierr = VecScale(VEC_VV(loc_it+1),1.0/tt);CHKERRQ(ierr);
    } else {
        /* This happens when the solution is exactly reached. */
        /* So there is no new direction... */
          ierr   = VecSet(VEC_TEMP,0.0);CHKERRQ(ierr); /* set VEC_TEMP to 0 */
          hapend = PETSC_TRUE;
    }
    /* note that for FGMRES we could get HES(loc_it+1, loc_it)  = 0 and the
       current solution would not be exact if HES was singular.  Note that
       HH non-singular implies that HES is no singular, and HES is guaranteed
//.........这里部分代码省略.........

PetscErrorCode PCBDDCNullSpaceAssembleCorrection(PC pc,IS local_dofs)
{
  PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
  PC_IS          *pcis = (PC_IS*)pc->data;
  Mat_IS*        matis = (Mat_IS*)pc->pmat->data;
  KSP            *local_ksp;
  PC             newpc;
  NullSpaceCorrection_ctx  shell_ctx;
  Mat            local_mat,local_pmat,small_mat,inv_small_mat;
  MatStructure   local_mat_struct;
  Vec            work1,work2;
  const Vec      *nullvecs;
  VecScatter     scatter_ctx;
  IS             is_aux;
  MatFactorInfo  matinfo;
  PetscScalar    *basis_mat,*Kbasis_mat,*array,*array_mat;
  PetscScalar    one = 1.0,zero = 0.0, m_one = -1.0;
  PetscInt       basis_dofs,basis_size,nnsp_size,i,k,n_I,n_R;
  PetscBool      nnsp_has_cnst;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  /* Infer the local solver */
  ierr = ISGetSize(local_dofs,&basis_dofs);CHKERRQ(ierr);
  ierr = VecGetSize(pcis->vec1_D,&n_I);CHKERRQ(ierr);
  ierr = VecGetSize(pcbddc->vec1_R,&n_R);CHKERRQ(ierr);
  if (basis_dofs == n_I) {
    /* Dirichlet solver */
    local_ksp = &pcbddc->ksp_D;
  } else if (basis_dofs == n_R) {
    /* Neumann solver */
    local_ksp = &pcbddc->ksp_R;
  } else {
    SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in %s: unknown local IS size %d. n_I=%d, n_R=%d)\n",__FUNCT__,basis_dofs,n_I,n_R);
  }
  ierr = KSPGetOperators(*local_ksp,&local_mat,&local_pmat,&local_mat_struct);CHKERRQ(ierr);

  /* Get null space vecs */
  ierr = MatNullSpaceGetVecs(pcbddc->NullSpace,&nnsp_has_cnst,&nnsp_size,&nullvecs);CHKERRQ(ierr);
  basis_size = nnsp_size;
  if (nnsp_has_cnst) {
    basis_size++;
  }

  /* Create shell ctx */
  ierr = PetscMalloc(sizeof(*shell_ctx),&shell_ctx);CHKERRQ(ierr);

  /* Create work vectors in shell context */
  ierr = VecCreate(PETSC_COMM_SELF,&shell_ctx->work_small_1);CHKERRQ(ierr);
  ierr = VecSetSizes(shell_ctx->work_small_1,basis_size,basis_size);CHKERRQ(ierr);
  ierr = VecSetType(shell_ctx->work_small_1,VECSEQ);CHKERRQ(ierr);
  ierr = VecDuplicate(shell_ctx->work_small_1,&shell_ctx->work_small_2);CHKERRQ(ierr);
  ierr = VecCreate(PETSC_COMM_SELF,&shell_ctx->work_full_1);CHKERRQ(ierr);
  ierr = VecSetSizes(shell_ctx->work_full_1,basis_dofs,basis_dofs);CHKERRQ(ierr);
  ierr = VecSetType(shell_ctx->work_full_1,VECSEQ);CHKERRQ(ierr);
  ierr = VecDuplicate(shell_ctx->work_full_1,&shell_ctx->work_full_2);CHKERRQ(ierr);

  /* Allocate workspace */
  ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_dofs,basis_size,NULL,&shell_ctx->basis_mat );CHKERRQ(ierr);
  ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_dofs,basis_size,NULL,&shell_ctx->Kbasis_mat);CHKERRQ(ierr);
  ierr = MatDenseGetArray(shell_ctx->basis_mat,&basis_mat);CHKERRQ(ierr);
  ierr = MatDenseGetArray(shell_ctx->Kbasis_mat,&Kbasis_mat);CHKERRQ(ierr);

  /* Restrict local null space on selected dofs (Dirichlet or Neumann)
     and compute matrices N and K*N */
  ierr = VecDuplicate(shell_ctx->work_full_1,&work1);CHKERRQ(ierr);
  ierr = VecDuplicate(shell_ctx->work_full_1,&work2);CHKERRQ(ierr);
  ierr = VecScatterCreate(pcis->vec1_N,local_dofs,work1,(IS)0,&scatter_ctx);CHKERRQ(ierr);
  for (k=0;k<nnsp_size;k++) {
    ierr = VecScatterBegin(matis->ctx,nullvecs[k],pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = VecScatterEnd(matis->ctx,nullvecs[k],pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = VecPlaceArray(work1,(const PetscScalar*)&basis_mat[k*basis_dofs]);CHKERRQ(ierr);
    ierr = VecScatterBegin(scatter_ctx,pcis->vec1_N,work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = VecScatterEnd(scatter_ctx,pcis->vec1_N,work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = VecPlaceArray(work2,(const PetscScalar*)&Kbasis_mat[k*basis_dofs]);CHKERRQ(ierr);
    ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr);
    ierr = VecResetArray(work1);CHKERRQ(ierr);
    ierr = VecResetArray(work2);CHKERRQ(ierr);
  }
  if (nnsp_has_cnst) {
    ierr = VecPlaceArray(work1,(const PetscScalar*)&basis_mat[k*basis_dofs]);CHKERRQ(ierr);
    ierr = VecSet(work1,one);CHKERRQ(ierr);
    ierr = VecPlaceArray(work2,(const PetscScalar*)&Kbasis_mat[k*basis_dofs]);CHKERRQ(ierr);
    ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr);
    ierr = VecResetArray(work1);CHKERRQ(ierr);
    ierr = VecResetArray(work2);CHKERRQ(ierr);
  }
  ierr = VecDestroy(&work1);CHKERRQ(ierr);
  ierr = VecDestroy(&work2);CHKERRQ(ierr);
  ierr = VecScatterDestroy(&scatter_ctx);CHKERRQ(ierr);
  ierr = MatDenseRestoreArray(shell_ctx->basis_mat,&basis_mat);CHKERRQ(ierr);
  ierr = MatDenseRestoreArray(shell_ctx->Kbasis_mat,&Kbasis_mat);CHKERRQ(ierr);

  /* Assemble another Mat object in shell context */
  ierr = MatTransposeMatMult(shell_ctx->basis_mat,shell_ctx->Kbasis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&small_mat);CHKERRQ(ierr);
  ierr = MatFactorInfoInitialize(&matinfo);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,basis_size,0,1,&is_aux);CHKERRQ(ierr);
  ierr = MatLUFactor(small_mat,is_aux,is_aux,&matinfo);CHKERRQ(ierr);
  ierr = ISDestroy(&is_aux);CHKERRQ(ierr);
  ierr = PetscMalloc(basis_size*basis_size*sizeof(PetscScalar),&array_mat);CHKERRQ(ierr);
//.........这里部分代码省略.........

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

  /* MatShift() */
  rval = 10*s1norm;
  ierr = MatShift(A,rval);CHKERRQ(ierr);
  ierr = MatShift(B,rval);CHKERRQ(ierr);

  /* Test MatTranspose() */
  ierr = MatTranspose(A,MAT_REUSE_MATRIX,&A);CHKERRQ(ierr);
  ierr = MatTranspose(B,MAT_REUSE_MATRIX,&B);CHKERRQ(ierr);

  /* Now do MatGetValues()  */
  for (i=0; i<30; i++) {
    ierr = PetscRandomGetValue(rdm,&rval);CHKERRQ(ierr);
    cols[0] = (PetscInt)(PetscRealPart(rval)*M);
    ierr = PetscRandomGetValue(rdm,&rval);CHKERRQ(ierr);
    cols[1] = (PetscInt)(PetscRealPart(rval)*M);
    ierr = PetscRandomGetValue(rdm,&rval);CHKERRQ(ierr);
    rows[0] = (PetscInt)(PetscRealPart(rval)*M);
    ierr = PetscRandomGetValue(rdm,&rval);CHKERRQ(ierr);
    rows[1] = (PetscInt)(PetscRealPart(rval)*M);
    ierr = MatGetValues(A,2,rows,2,cols,vals1);CHKERRQ(ierr);
    ierr = MatGetValues(B,2,rows,2,cols,vals2);CHKERRQ(ierr);
    ierr = PetscMemcmp(vals1,vals2,4*sizeof(PetscScalar),&flg);CHKERRQ(ierr);
    if (!flg) {
      ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatGetValues bs = %D\n",bs);CHKERRQ(ierr);
    }
  }

  /* Test MatMult(), MatMultAdd() */
  for (i=0; i<40; i++) {
    ierr = VecSetRandom(xx,rdm);CHKERRQ(ierr);
    ierr = VecSet(s2,0.0);CHKERRQ(ierr);
    ierr = MatMult(A,xx,s1);CHKERRQ(ierr);
    ierr = MatMultAdd(A,xx,s2,s2);CHKERRQ(ierr);
    ierr = VecNorm(s1,NORM_2,&s1norm);CHKERRQ(ierr);
    ierr = VecNorm(s2,NORM_2,&s2norm);CHKERRQ(ierr);
    rnorm = s2norm-s1norm;
    if (rnorm<-tol || rnorm>tol) {
      ierr = PetscPrintf(PETSC_COMM_SELF,"MatMult not equal to MatMultAdd Norm1=%e Norm2=%e bs = %D\n",s1norm,s2norm,bs);CHKERRQ(ierr);
    }
  }

  /* Test MatMult() */
  ierr = MatMultEqual(A,B,10,&flg);CHKERRQ(ierr);
  if (!flg){
    ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatMult()\n");CHKERRQ(ierr);
  }

  /* Test MatMultAdd() */
  ierr = MatMultAddEqual(A,B,10,&flg);CHKERRQ(ierr);
  if (!flg){
    ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatMultAdd()\n");CHKERRQ(ierr);
  }

  /* Test MatMultTranspose() */
  ierr = MatMultTransposeEqual(A,B,10,&flg);CHKERRQ(ierr);
  if (!flg){
    ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatMultTranspose()\n");CHKERRQ(ierr);
  }

  /* Test MatMultTransposeAdd() */
  ierr = MatMultTransposeAddEqual(A,B,10,&flg);CHKERRQ(ierr);
  if (!flg){
    ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatMultTransposeAdd()\n");CHKERRQ(ierr);
  }

PetscErrorCode PCBDDCNullSpaceAdaptGlobal(PC pc)
{
  PC_IS*         pcis = (PC_IS*)  (pc->data);
  PC_BDDC*       pcbddc = (PC_BDDC*)(pc->data);
  KSP            inv_change;
  PC             pc_change;
  const Vec      *nsp_vecs;
  Vec            *new_nsp_vecs;
  PetscInt       i,nsp_size,new_nsp_size,start_new;
  PetscBool      nsp_has_cnst;
  MatNullSpace   new_nsp;
  PetscErrorCode ierr;
  MPI_Comm       comm;

  PetscFunctionBegin;
  /* create KSP for change of basis */
  ierr = KSPCreate(PETSC_COMM_SELF,&inv_change);CHKERRQ(ierr);
  ierr = KSPSetOperators(inv_change,pcbddc->ChangeOfBasisMatrix,pcbddc->ChangeOfBasisMatrix,SAME_PRECONDITIONER);CHKERRQ(ierr);
  ierr = KSPSetType(inv_change,KSPPREONLY);CHKERRQ(ierr);
  ierr = KSPGetPC(inv_change,&pc_change);CHKERRQ(ierr);
  ierr = PCSetType(pc_change,PCLU);CHKERRQ(ierr);
  ierr = KSPSetUp(inv_change);CHKERRQ(ierr);
  /* get nullspace and transform it */
  ierr = MatNullSpaceGetVecs(pcbddc->NullSpace,&nsp_has_cnst,&nsp_size,&nsp_vecs);CHKERRQ(ierr);
  new_nsp_size = nsp_size;
  if (nsp_has_cnst) {
    new_nsp_size++;
  }
  ierr = PetscMalloc(new_nsp_size*sizeof(Vec),&new_nsp_vecs);CHKERRQ(ierr);
  for (i=0;i<new_nsp_size;i++) {
    ierr = VecDuplicate(pcis->vec1_global,&new_nsp_vecs[i]);CHKERRQ(ierr);
  }
  start_new = 0;
  if (nsp_has_cnst) {
    start_new = 1;
    ierr = VecSet(new_nsp_vecs[0],1.0);CHKERRQ(ierr);
    ierr = VecSet(pcis->vec1_B,1.0);CHKERRQ(ierr);
    ierr = KSPSolve(inv_change,pcis->vec1_B,pcis->vec1_B);
    ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,new_nsp_vecs[0],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
    ierr = VecScatterEnd  (pcis->global_to_B,pcis->vec1_B,new_nsp_vecs[0],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
  }
  for (i=0;i<nsp_size;i++) {
    ierr = VecCopy(nsp_vecs[i],new_nsp_vecs[i+start_new]);CHKERRQ(ierr);
    ierr = VecScatterBegin(pcis->global_to_B,nsp_vecs[i],pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = VecScatterEnd  (pcis->global_to_B,nsp_vecs[i],pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = KSPSolve(inv_change,pcis->vec1_B,pcis->vec1_B);
    ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,new_nsp_vecs[i+start_new],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
    ierr = VecScatterEnd  (pcis->global_to_B,pcis->vec1_B,new_nsp_vecs[i+start_new],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
  }
  ierr = PCBDDCOrthonormalizeVecs(new_nsp_size,new_nsp_vecs);CHKERRQ(ierr);
#if 0
  PetscBool nsp_t=PETSC_FALSE;
  ierr = MatNullSpaceTest(pcbddc->NullSpace,pc->pmat,&nsp_t);CHKERRQ(ierr);
  printf("Original Null Space test: %d\n",nsp_t);
  Mat temp_mat;
  Mat_IS* matis = (Mat_IS*)pc->pmat->data;
    temp_mat = matis->A;
    matis->A = pcbddc->local_mat;
    pcbddc->local_mat = temp_mat;
  ierr = MatNullSpaceTest(pcbddc->NullSpace,pc->pmat,&nsp_t);CHKERRQ(ierr);
  printf("Original Null Space, mat changed test: %d\n",nsp_t);
  {
    PetscReal test_norm;
    for (i=0;i<new_nsp_size;i++) {
      ierr = MatMult(pc->pmat,new_nsp_vecs[i],pcis->vec1_global);CHKERRQ(ierr);
      ierr = VecNorm(pcis->vec1_global,NORM_2,&test_norm);CHKERRQ(ierr);
      if (test_norm > 1.e-12) {
        printf("------------ERROR VEC %d------------------\n",i);
        ierr = VecView(pcis->vec1_global,PETSC_VIEWER_STDOUT_WORLD);
        printf("------------------------------------------\n");
      }
    }
  }
#endif

  ierr = KSPDestroy(&inv_change);CHKERRQ(ierr);
  ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
  ierr = MatNullSpaceCreate(comm,PETSC_FALSE,new_nsp_size,new_nsp_vecs,&new_nsp);CHKERRQ(ierr);
  ierr = PCBDDCSetNullSpace(pc,new_nsp);CHKERRQ(ierr);
  ierr = MatNullSpaceDestroy(&new_nsp);CHKERRQ(ierr);
#if 0
  ierr = MatNullSpaceTest(pcbddc->NullSpace,pc->pmat,&nsp_t);CHKERRQ(ierr);
  printf("New Null Space, mat changed: %d\n",nsp_t);
    temp_mat = matis->A;
    matis->A = pcbddc->local_mat;
    pcbddc->local_mat = temp_mat;
  ierr = MatNullSpaceTest(pcbddc->NullSpace,pc->pmat,&nsp_t);CHKERRQ(ierr);
  printf("New Null Space, mat original: %d\n",nsp_t);
#endif

  for (i=0;i<new_nsp_size;i++) {
    ierr = VecDestroy(&new_nsp_vecs[i]);CHKERRQ(ierr);
  }
  ierr = PetscFree(new_nsp_vecs);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode PCBDDCNullSpaceAssembleCoarse(PC pc, MatNullSpace* CoarseNullSpace)
{
  PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
  Mat_IS         *matis = (Mat_IS*)pc->pmat->data;
  MatNullSpace   tempCoarseNullSpace;
  const Vec      *nsp_vecs;
  Vec            *coarse_nsp_vecs,