static char help[] = "Test file for the PCFactorSetShiftType()\n"; /* * Test file for the PCFactorSetShiftType() routine or -pc_factor_shift_type POSITIVE_DEFINITE option. * The test matrix is the example from Kershaw's paper [J.Comp.Phys 1978] * of a positive definite matrix for which ILU(0) will give a negative pivot. * This means that the CG method will break down; the Manteuffel shift * [Math. Comp. 1980] repairs this. * * Run the executable twice: * 1/ without options: the iterative method diverges because of an * indefinite preconditioner * 2/ with -pc_factor_shift_positive_definite option (or comment in the PCFactorSetShiftType() line below): * the method will now successfully converge. */ #include int main(int argc,char **argv) { KSP ksp; PC pc; Mat A,M; Vec X,B,D; MPI_Comm comm; PetscScalar v; KSPConvergedReason reason; PetscInt i,j,its; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr; comm = MPI_COMM_SELF; /* * Construct the Kershaw matrix * and a suitable rhs / initial guess */ ierr = MatCreateSeqAIJ(comm,4,4,4,0,&A);CHKERRQ(ierr); ierr = VecCreateSeq(comm,4,&B);CHKERRQ(ierr); ierr = VecDuplicate(B,&X);CHKERRQ(ierr); for (i=0; i<4; i++) { v = 3; ierr = MatSetValues(A,1,&i,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); v = 1; ierr = VecSetValues(B,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); ierr = VecSetValues(X,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); } i =0; v=0; ierr = VecSetValues(X,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); for (i=0; i<3; i++) { v = -2; j=i+1; ierr = MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); ierr = MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); } i=0; j=3; v=2; ierr = MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); ierr = MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = VecAssemblyBegin(B);CHKERRQ(ierr); ierr = VecAssemblyEnd(B);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"\nThe Kershaw matrix:\n\n");CHKERRQ(ierr); ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); /* * A Conjugate Gradient method * with ILU(0) preconditioning */ ierr = KSPCreate(comm,&ksp);CHKERRQ(ierr); ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr); ierr = KSPSetType(ksp,KSPCG);CHKERRQ(ierr); ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); /* * ILU preconditioner; * The iterative method will break down unless you comment in the SetShift * line below, or use the -pc_factor_shift_positive_definite option. * Run the code twice: once as given to see the negative pivot and the * divergence behaviour, then comment in the Shift line, or add the * command line option, and see that the pivots are all positive and * the method converges. */ ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); ierr = PCSetType(pc,PCICC);CHKERRQ(ierr); /* ierr = PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE);CHKERRQ(ierr); */ ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); ierr = KSPSetUp(ksp);CHKERRQ(ierr); /* * Now that the factorisation is done, show the pivots; * note that the last one is negative. This in itself is not an error, * but it will make the iterative method diverge. */ ierr = PCFactorGetMatrix(pc,&M);CHKERRQ(ierr); ierr = VecDuplicate(B,&D);CHKERRQ(ierr); ierr = MatGetDiagonal(M,D);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"\nPivots:\n\n");CHKERRQ(ierr); ierr = VecView(D,0);CHKERRQ(ierr); /* * Solve the system; * without the shift this will diverge with * an indefinite preconditioner */ ierr = KSPSolve(ksp,B,X);CHKERRQ(ierr); ierr = KSPGetConvergedReason(ksp,&reason);CHKERRQ(ierr); if (reason==KSP_DIVERGED_INDEFINITE_PC) { ierr = PetscPrintf(PETSC_COMM_WORLD,"\nDivergence because of indefinite preconditioner;\n");CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Run the executable again but with -pc_factor_shift_positive_definite option.\n");CHKERRQ(ierr); } else if (reason<0) { ierr = PetscPrintf(PETSC_COMM_WORLD,"\nOther kind of divergence: this should not happen.\n");CHKERRQ(ierr); } else { ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"\nConvergence in %d iterations.\n",(int)its);CHKERRQ(ierr); } ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr); ierr = KSPDestroy(&ksp);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr); ierr = VecDestroy(&B);CHKERRQ(ierr); ierr = VecDestroy(&X);CHKERRQ(ierr); ierr = VecDestroy(&D);CHKERRQ(ierr); ierr = PetscFinalize(); return ierr; } /*TEST test: filter: sed -e "s/in 5 iterations/in 4 iterations/g" TEST*/