xref: /petsc/src/mat/tests/ex258.c (revision 609caa7c8c030312b00807b4f015fd827bb80932)
1 static char help[] = "Test MatProductReplaceMats() \n\
2 Modified from the code contributed by Pierre Jolivet \n\n";
3 
4 #include <petscmat.h>
5 
main(int argc,char ** args)6 int main(int argc, char **args)
7 {
8   PetscInt  n = 2, convert;
9   Mat       A, B, Bdense, Conjugate;
10   PetscBool conjugate = PETSC_FALSE, equal, flg;
11 
12   PetscFunctionBeginUser;
13   PetscCall(PetscInitialize(&argc, &args, NULL, help));
14 
15   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
16   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
17   PetscCall(MatSetType(A, MATDENSE));
18   PetscCall(MatSetFromOptions(A));
19   PetscCall(MatSeqDenseSetPreallocation(A, NULL));
20   PetscCall(MatMPIDenseSetPreallocation(A, NULL));
21   PetscCall(MatSetRandom(A, NULL));
22   PetscCall(MatViewFromOptions(A, NULL, "-A_view"));
23   PetscCall(PetscOptionsGetBool(NULL, NULL, "-conjugate", &conjugate, NULL));
24 
25   for (convert = 0; convert < 2; convert++) {
26     /* convert dense matrix A to aij format */
27     if (convert) PetscCall(MatConvert(A, MATAIJ, MAT_INPLACE_MATRIX, &A));
28 
29     /* compute B = A^T * A or  B = A^H * A */
30     PetscCall(MatProductCreate(A, A, NULL, &B));
31 
32     flg = PETSC_FALSE;
33     PetscCall(PetscOptionsGetBool(NULL, NULL, "-atb", &flg, NULL));
34     if (flg) {
35       PetscCall(MatProductSetType(B, MATPRODUCT_AtB));
36     } else {
37       PetscCall(PetscOptionsGetBool(NULL, NULL, "-ptap", &flg, NULL));
38       if (flg) {
39         PetscCall(MatProductSetType(B, MATPRODUCT_PtAP));
40       } else {
41         PetscCall(PetscOptionsGetBool(NULL, NULL, "-abt", &flg, NULL));
42         if (flg) {
43           PetscCall(MatProductSetType(B, MATPRODUCT_ABt));
44         } else {
45           PetscCall(MatProductSetType(B, MATPRODUCT_AB));
46         }
47       }
48     }
49     PetscCall(MatProductSetFromOptions(B));
50     PetscCall(MatProductSymbolic(B));
51 
52     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &Conjugate));
53     if (conjugate) PetscCall(MatConjugate(Conjugate));
54 
55     /* replace input A by Conjugate */
56     PetscCall(MatProductReplaceMats(Conjugate, NULL, NULL, B));
57 
58     PetscCall(MatProductNumeric(B));
59     PetscCall(MatViewFromOptions(B, NULL, "-product_view"));
60 
61     PetscCall(MatDestroy(&Conjugate));
62     if (!convert) {
63       Bdense = B;
64       B      = NULL;
65     }
66   }
67 
68   /* Compare Bdense and B */
69   PetscCall(MatMultEqual(Bdense, B, 10, &equal));
70   PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Bdense != B");
71 
72   PetscCall(MatDestroy(&Bdense));
73   PetscCall(MatDestroy(&B));
74   PetscCall(MatDestroy(&A));
75   PetscCall(PetscFinalize());
76   return 0;
77 }
78 
79 /*TEST
80 
81    test:
82       suffix: 1
83       args: -conjugate false -atb
84       output_file: output/empty.out
85 
86    test:
87       suffix: 2
88       args: -conjugate true -atb
89       output_file: output/empty.out
90 
91    test:
92       suffix: 3
93       args: -conjugate false
94       output_file: output/empty.out
95 
96    test:
97       suffix: 4
98       args: -ptap
99       output_file: output/empty.out
100 
101    test:
102       suffix: 5
103       args: -abt
104       output_file: output/empty.out
105 
106    test:
107       suffix: 6
108       nsize: 2
109       args: -conjugate false -atb
110       output_file: output/empty.out
111 
112    test:
113       suffix: 7
114       nsize: 2
115       args: -conjugate true -atb
116       output_file: output/empty.out
117 
118    test:
119       suffix: 8
120       nsize: 2
121       args: -conjugate false
122       output_file: output/empty.out
123 
124    test:
125       suffix: 9
126       nsize: 2
127       args: -ptap
128       output_file: output/empty.out
129 
130    test:
131       suffix: 10
132       nsize: 2
133       args: -abt
134       output_file: output/empty.out
135 
136    test:
137       suffix: 11
138       nsize: 2
139       args: -conjugate true -atb -mat_product_algorithm backend
140       TODO: bug: MatProductReplaceMats() does not change the product for this test
141 
142 TEST*/
143