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*/