Actual source code: ex25.c

  1: static char help[] = "Tests CG, MINRES and SYMMLQ on the symmetric indefinite matrices: afiro \n\n";

  3: #include <petscksp.h>

  5: int main(int argc, char **args)
  6: {
  7:   Mat         C;
  8:   PetscScalar none = -1.0;
  9:   PetscMPIInt rank, size;
 10:   PetscInt    its, k;
 11:   PetscReal   err_norm, res_norm;
 12:   Vec         x, b, u, u_tmp;
 13:   PC          pc;
 14:   KSP         ksp;
 15:   PetscViewer view;
 16:   char        filein[128]; /* input file name */

 18:   PetscFunctionBeginUser;
 19:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 20:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 21:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));

 23:   /* Load the binary data file "filein". Set runtime option: -f filein */
 24:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n Load dataset ...\n"));
 25:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filein, sizeof(filein), NULL));
 26:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filein, FILE_MODE_READ, &view));
 27:   PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
 28:   PetscCall(MatSetType(C, MATMPISBAIJ));
 29:   PetscCall(MatLoad(C, view));
 30:   PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
 31:   PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
 32:   PetscCall(VecLoad(b, view));
 33:   PetscCall(VecLoad(u, view));
 34:   PetscCall(PetscViewerDestroy(&view));
 35:   /* PetscCall(VecView(b,VIEWER_STDOUT_WORLD)); */
 36:   /* PetscCall(MatView(C,VIEWER_STDOUT_WORLD)); */

 38:   PetscCall(VecDuplicate(u, &u_tmp));

 40:   /* Check accuracy of the data */
 41:   /*
 42:   PetscCall(MatMult(C,u,u_tmp));
 43:   PetscCall(VecAXPY(u_tmp,none,b));
 44:   PetscCall(VecNorm(u_tmp,NORM_2,&res_norm));
 45:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Accuracy of the loading data: | b - A*u |_2 : %g \n",(double)res_norm));
 46:   */

 48:   /* Setup and solve for system */
 49:   PetscCall(VecDuplicate(b, &x));
 50:   for (k = 0; k < 3; k++) {
 51:     if (k == 0) { /* CG  */
 52:       PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 53:       PetscCall(KSPSetType(ksp, KSPCG));
 54:       PetscCall(KSPSetOperators(ksp, C, C));
 55:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n CG: \n"));
 56:     } else if (k == 1) { /* MINRES */
 57:       PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 58:       PetscCall(KSPSetType(ksp, KSPMINRES));
 59:       PetscCall(KSPSetOperators(ksp, C, C));
 60:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n MINRES: \n"));
 61:     } else { /* SYMMLQ */
 62:       PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 63:       PetscCall(KSPSetOperators(ksp, C, C));
 64:       PetscCall(KSPSetType(ksp, KSPSYMMLQ));
 65:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n SYMMLQ: \n"));
 66:     }

 68:     PetscCall(KSPGetPC(ksp, &pc));
 69:     PetscCall(PCSetType(pc, PCNONE));

 71:     /*
 72:     Set runtime options, e.g.,
 73:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
 74:                          -pc_type jacobi -pc_jacobi_type rowmax
 75:     These options will override those specified above as long as
 76:     KSPSetFromOptions() is called _after_ any other customization routines.
 77:     */
 78:     PetscCall(KSPSetFromOptions(ksp));

 80:     /* Solve linear system; */
 81:     PetscCall(KSPSolve(ksp, b, x));
 82:     PetscCall(KSPGetIterationNumber(ksp, &its));

 84:     /* Check error */
 85:     PetscCall(VecCopy(u, u_tmp));
 86:     PetscCall(VecAXPY(u_tmp, none, x));
 87:     PetscCall(VecNorm(u_tmp, NORM_2, &err_norm));
 88:     PetscCall(MatMult(C, x, u_tmp));
 89:     PetscCall(VecAXPY(u_tmp, none, b));
 90:     PetscCall(VecNorm(u_tmp, NORM_2, &res_norm));

 92:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its));
 93:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm: %g;", (double)res_norm));
 94:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Error norm: %g.\n", (double)err_norm));

 96:     PetscCall(KSPDestroy(&ksp));
 97:   }

 99:   /*
100:        Free work space.  All PETSc objects should be destroyed when they
101:        are no longer needed.
102:   */
103:   PetscCall(VecDestroy(&b));
104:   PetscCall(VecDestroy(&u));
105:   PetscCall(VecDestroy(&x));
106:   PetscCall(VecDestroy(&u_tmp));
107:   PetscCall(MatDestroy(&C));

109:   PetscCall(PetscFinalize());
110:   return 0;
111: }

113: /*TEST

115:     test:
116:       args: -f ${DATAFILESPATH}/matrices/indefinite/afiro -ksp_rtol 1.e-3
117:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)

119: TEST*/