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