Actual source code: ex5.c


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

  4: #include <petscksp.h>

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

 13:   PetscFunctionBeginUser;
 14:   PetscCall(PetscInitialize(&argc, &args, (char *)0, 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:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
 29:   PetscCall(PetscRandomSetType(rand, PETSCRANDER48));

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

 54: /*TEST

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

 60: TEST*/