Actual source code: ex82.c

  1: static const char help[] = "Uses KSPComputeRitz() on a matrix loaded from disk\n";

  3: #include <petscksp.h>

  5: int main(int argc, char **argv)
  6: {
  7:   Mat         A;
  8:   KSP         ksp;
  9:   char        file[PETSC_MAX_PATH_LEN];
 10:   PetscReal  *tetar, *tetai;
 11:   Vec         b, x, *S;
 12:   PetscInt    i, N = 10, Na = N;
 13:   PetscViewer fd;
 14:   PC          pc;
 15:   PetscBool   harmonic = PETSC_FALSE;

 17:   PetscFunctionBeginUser;
 18:   PetscCall(PetscInitialize(&argc, &argv, 0, help));

 20:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), NULL));
 21:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
 22:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 23:   PetscCall(MatLoad(A, fd));
 24:   PetscCall(PetscViewerDestroy(&fd));
 25:   PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));

 27:   PetscCall(MatCreateVecs(A, &b, &x));
 28:   PetscCall(VecSetRandom(b, NULL));
 29:   PetscCall(VecDuplicateVecs(b, N, &S));
 30:   PetscCall(PetscMalloc2(N, &tetar, N, &tetai));

 32:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 33:   PetscCall(KSPSetType(ksp, KSPGMRES));
 34:   PetscCall(KSPSetTolerances(ksp, 10000 * PETSC_MACHINE_EPSILON, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
 35:   PetscCall(KSPGetPC(ksp, &pc));
 36:   PetscCall(PCSetType(pc, PCNONE));
 37:   PetscCall(KSPSetComputeRitz(ksp, PETSC_TRUE));
 38:   PetscCall(KSPSetOperators(ksp, A, A));
 39:   PetscCall(KSPSetFromOptions(ksp));
 40:   PetscCall(KSPSolve(ksp, b, x));
 41:   PetscCall(PetscOptionsHasName(NULL, NULL, "-harmonic", &harmonic));
 42:   PetscCall(KSPComputeRitz(ksp, harmonic ? PETSC_FALSE : PETSC_TRUE, PETSC_TRUE, &Na, S, tetar, tetai));

 44:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%% Number of Ritz pairs %" PetscInt_FMT "\n", Na));
 45:   for (i = 0; i < Na; i++) {
 46:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%% Eigenvalue(s)  %g ", (double)tetar[i]));
 47:     if (tetai[i]) {
 48:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%+gi", (double)tetai[i]));
 49: #if !defined(PETSC_USE_COMPLEX)
 50:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  %g ", (double)tetar[i]));
 51:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%+gi", (double)-tetai[i]));
 52: #endif
 53:     }
 54:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n"));
 55:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%% Eigenvector\n"));
 56:     PetscCall(VecView(S[i], PETSC_VIEWER_STDOUT_WORLD));
 57: #if !defined(PETSC_USE_COMPLEX)
 58:     if (tetai[i]) {
 59:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%% Imaginary part of Eigenvector\n"));
 60:       PetscCall(VecView(S[i + 1], PETSC_VIEWER_STDOUT_WORLD));
 61:       i++;
 62:     }
 63: #endif
 64:   }

 66:   PetscCall(PetscFree2(tetar, tetai));
 67:   PetscCall(VecDestroy(&x));
 68:   PetscCall(VecDestroy(&b));
 69:   PetscCall(VecDestroyVecs(N, &S));
 70:   PetscCall(KSPDestroy(&ksp));
 71:   PetscCall(MatDestroy(&A));
 72:   PetscCall(PetscFinalize());
 73:   return 0;
 74: }

 76: /*TEST

 78:     test:
 79:       requires: !defined(PETSC_USE_64BIT_INDICES) !complex double
 80:       args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/ritz2 -ksp_monitor

 82:     test:
 83:       suffix: 2
 84:       requires: !defined(PETSC_USE_64BIT_INDICES) !complex double
 85:       args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/ritz5 -ksp_monitor

 87:     test:
 88:       suffix: harmonic
 89:       requires: !defined(PETSC_USE_64BIT_INDICES) !complex double
 90:       args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/ritz5 -ksp_monitor -harmonic

 92: TEST*/