1 static char help[] = "Tests solving linear system with KSPFGMRES + PCSOR (omega != 1) on a matrix obtained from MatTransposeMatMult.\n\n";
2
3 #include <petscksp.h>
4
main(int argc,char ** args)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];
14
15 PetscFunctionBeginUser;
16 PetscCall(PetscInitialize(&argc, &args, NULL, 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 }
54
55 /*TEST
56
57 test:
58 nsize: 1
59 suffix: 1
60
61 test:
62 nsize: 1
63 suffix: 1_mpiaij
64 args: -mtype mpiaij
65
66 test:
67 nsize: 3
68 suffix: 2
69
70 TEST*/
71