Actual source code: ex5.c

  1: static char help[] = "Tests MATLMVM classes.\n\n";

  3: #include <petscksp.h>

  5: int main(int argc, char **args)
  6: {
  7:   PC          P;
  8:   Mat         A;
  9:   Vec         x, x2, f, u, b;
 10:   PetscInt    n = 4, nup = 10;
 11:   PetscRandom rand;

 13:   PetscFunctionBeginUser;
 14:   PetscCall(PetscInitialize(&argc, &args, NULL, help));

 16:   /* make sure LMVM classes are registered */
 17:   PetscCall(KSPInitializePackage());

 19:   /* create matrix */
 20:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 21:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
 22:   PetscCall(MatSetType(A, MATLMVMBFGS));
 23:   PetscCall(MatSetFromOptions(A));
 24:   PetscCall(MatSetUp(A));
 25:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 26:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 28:   /* create preconditioner wrapping MatSolve */
 29:   PetscCall(PCCreate(PETSC_COMM_WORLD, &P));
 30:   PetscCall(PCSetType(P, PCMAT));
 31:   PetscCall(PCSetOperators(P, A, A));
 32:   PetscCall(PCSetUp(P));
 33:   PetscCall(PCView(P, NULL));

 35:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
 36:   PetscCall(PetscRandomSetType(rand, PETSCRANDER48));

 38:   /* create vectors */
 39:   PetscCall(MatCreateVecs(A, &x, &f));
 40:   PetscCall(VecDuplicate(x, &x2));
 41:   PetscCall(VecDuplicate(x, &u));
 42:   PetscCall(VecDuplicate(f, &b));
 43:   PetscCall(VecSetRandom(u, rand));

 45:   /* Test various routines */
 46:   for (PetscInt i = 0; i < nup; i++) {
 47:     PetscReal err;

 49:     PetscCall(VecSetRandom(x, rand));
 50:     PetscCall(VecSetRandom(f, rand));
 51:     PetscCall(MatLMVMUpdate(A, x, f));
 52:     PetscCall(MatView(A, NULL));
 53:     PetscCall(MatMult(A, u, b));
 54:     PetscCall(MatSolve(A, b, x));
 55:     PetscCall(PCApply(P, b, x2));
 56:     PetscCall(VecAXPY(x2, -1.0, x));
 57:     PetscCall(VecNorm(x2, NORM_2, &err));
 58:     if (err > PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error %" PetscInt_FMT ": %g\n", i, (double)err));
 59:   }

 61:   /* cleanup */
 62:   PetscCall(PetscRandomDestroy(&rand));
 63:   PetscCall(VecDestroy(&u));
 64:   PetscCall(VecDestroy(&b));
 65:   PetscCall(VecDestroy(&x));
 66:   PetscCall(VecDestroy(&x2));
 67:   PetscCall(VecDestroy(&f));
 68:   PetscCall(MatDestroy(&A));
 69:   PetscCall(PCDestroy(&P));
 70:   PetscCall(PetscFinalize());
 71:   return 0;
 72: }

 74: /*TEST

 76:     test:
 77:       requires: !complex
 78:       args: -mat_type {{lmvmdfp lmvmbfgs lmvmsr1 lmvmbroyden lmvmbadbroyden lmvmsymbroyden lmvmsymbadbroyden lmvmdiagbroyden}separate output}

 80: TEST*/