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;

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

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

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

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

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

 66:   PetscFree2(tetar, tetai);
 67:   VecDestroy(&x);
 68:   VecDestroy(&b);
 69:   VecDestroyVecs(N, &S);
 70:   KSPDestroy(&ksp);
 71:   MatDestroy(&A);
 72:   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*/