Actual source code: ex34.c
2: static char help[] = "Tests solving linear system with KSPFGMRES + PCSOR (omega != 1) on a matrix obtained from MatTransposeMatMult.\n\n";
4: #include <petscksp.h>
6: int main(int argc, char **args)
7: {
8: Mat A, Ad, B;
9: PetscInt N = 10, M = 3;
10: PetscBool no_inodes = PETSC_TRUE, flg;
11: KSP ksp;
12: PC pc;
13: Vec x, y;
14: char mtype[256];
16: PetscFunctionBeginUser;
17: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
18: PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
19: PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
20: PetscCall(PetscOptionsGetString(NULL, NULL, "-mtype", mtype, sizeof(mtype), &flg));
21: PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_inodes", &no_inodes, NULL));
22: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, NULL, &Ad));
23: PetscCall(MatSetRandom(Ad, NULL));
24: PetscCall(MatConvert(Ad, flg ? mtype : MATAIJ, MAT_INITIAL_MATRIX, &A));
25: PetscCall(MatProductCreate(A, A, NULL, &B));
26: PetscCall(MatProductSetType(B, MATPRODUCT_AtB));
27: PetscCall(MatProductSetAlgorithm(B, "default"));
28: PetscCall(MatProductSetFill(B, PETSC_DEFAULT));
29: PetscCall(MatProductSetFromOptions(B));
30: PetscCall(MatProductSymbolic(B));
31: if (no_inodes) PetscCall(MatSetOption(B, MAT_USE_INODES, PETSC_FALSE));
32: PetscCall(MatProductNumeric(B));
33: PetscCall(MatTransposeMatMultEqual(A, A, B, 10, &flg));
34: if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Wrong MatTransposeMat"));
35: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
36: PetscCall(KSPSetOperators(ksp, B, B));
37: PetscCall(KSPGetPC(ksp, &pc));
38: PetscCall(PCSetType(pc, PCSOR));
39: PetscCall(PCSORSetOmega(pc, 1.1));
40: PetscCall(KSPSetUp(ksp));
41: PetscCall(KSPView(ksp, NULL));
42: PetscCall(KSPGetPC(ksp, &pc));
43: PetscCall(MatCreateVecs(B, &y, &x));
44: PetscCall(VecSetRandom(x, NULL));
45: PetscCall(PCApply(pc, x, y));
46: PetscCall(KSPDestroy(&ksp));
47: PetscCall(VecDestroy(&x));
48: PetscCall(VecDestroy(&y));
49: PetscCall(MatDestroy(&B));
50: PetscCall(MatDestroy(&Ad));
51: PetscCall(MatDestroy(&A));
52: PetscCall(PetscFinalize());
53: return 0;
54: }
56: /*TEST
58: test:
59: nsize: 1
60: suffix: 1
62: test:
63: nsize: 1
64: suffix: 1_mpiaij
65: args: -mtype mpiaij
67: test:
68: nsize: 3
69: suffix: 2
71: TEST*/