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:   Mat         A;
  8:   Vec         x, f, u, b;
  9:   PetscInt    n = 4, nup = 10;
 10:   PetscRandom rand;

 12:   PetscFunctionBeginUser;
 13:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));

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

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

 27:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
 28:   PetscCall(PetscRandomSetType(rand, PETSCRANDER48));

 30:   /* create vectors */
 31:   PetscCall(MatCreateVecs(A, &x, &f));
 32:   PetscCall(VecDuplicate(x, &u));
 33:   PetscCall(VecDuplicate(f, &b));
 34:   PetscCall(VecSetRandom(u, rand));
 35:   for (PetscInt i = 0; i < nup; i++) {
 36:     PetscCall(VecSetRandom(x, rand));
 37:     PetscCall(VecSetRandom(f, rand));
 38:     PetscCall(MatLMVMUpdate(A, x, f));
 39:     PetscCall(MatView(A, NULL));
 40:     PetscCall(MatMult(A, u, b));
 41:     PetscCall(MatSolve(A, b, x));
 42:   }
 43:   PetscCall(PetscRandomDestroy(&rand));
 44:   PetscCall(VecDestroy(&u));
 45:   PetscCall(VecDestroy(&b));
 46:   PetscCall(VecDestroy(&x));
 47:   PetscCall(VecDestroy(&f));
 48:   PetscCall(MatDestroy(&A));
 49:   PetscCall(PetscFinalize());
 50:   return 0;
 51: }

 53: /*TEST

 55:     test:
 56:       requires: !complex
 57:       args: -mat_type {{lmvmdfp lmvmbfgs lmvmsr1 lmvmbroyden lmvmbadbroyden lmvmsymbroyden lmvmsymbadbroyden lmvmdiagbroyden}separate output}

 59: TEST*/