1 static char help[] = "Test MatTransposeColoring for SeqAIJ matrices. Used for '-matmattransmult_color' on MatMatTransposeMult \n\n";
2
3 #include <petscmat.h>
4 #include <petsc/private/matimpl.h> /* Need struct _p_MatTransposeColoring for this test. */
5
main(int argc,char ** argv)6 int main(int argc, char **argv)
7 {
8 Mat A, R, C, C_dense, C_sparse, Rt_dense, P, PtAP;
9 PetscInt row, col, m, n;
10 MatScalar one = 1.0, val;
11 MatColoring mc;
12 MatTransposeColoring matcoloring = 0;
13 ISColoring iscoloring;
14 PetscBool equal;
15 PetscMPIInt size;
16
17 PetscFunctionBeginUser;
18 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
19 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
20 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
21
22 /* Create seqaij A */
23 PetscCall(MatCreate(PETSC_COMM_SELF, &A));
24 PetscCall(MatSetSizes(A, 4, 4, 4, 4));
25 PetscCall(MatSetType(A, MATSEQAIJ));
26 PetscCall(MatSetFromOptions(A));
27 PetscCall(MatSetUp(A));
28 row = 0;
29 col = 0;
30 val = 1.0;
31 PetscCall(MatSetValues(A, 1, &row, 1, &col, &val, ADD_VALUES));
32 row = 1;
33 col = 3;
34 val = 2.0;
35 PetscCall(MatSetValues(A, 1, &row, 1, &col, &val, ADD_VALUES));
36 row = 2;
37 col = 2;
38 val = 3.0;
39 PetscCall(MatSetValues(A, 1, &row, 1, &col, &val, ADD_VALUES));
40 row = 3;
41 col = 0;
42 val = 4.0;
43 PetscCall(MatSetValues(A, 1, &row, 1, &col, &val, ADD_VALUES));
44 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
45 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
46 PetscCall(MatSetOptionsPrefix(A, "A_"));
47 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
48 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
49
50 /* Create seqaij R */
51 PetscCall(MatCreate(PETSC_COMM_SELF, &R));
52 PetscCall(MatSetSizes(R, 2, 4, 2, 4));
53 PetscCall(MatSetType(R, MATSEQAIJ));
54 PetscCall(MatSetFromOptions(R));
55 PetscCall(MatSetUp(R));
56 row = 0;
57 col = 0;
58 PetscCall(MatSetValues(R, 1, &row, 1, &col, &one, ADD_VALUES));
59 row = 0;
60 col = 1;
61 PetscCall(MatSetValues(R, 1, &row, 1, &col, &one, ADD_VALUES));
62
63 row = 1;
64 col = 1;
65 PetscCall(MatSetValues(R, 1, &row, 1, &col, &one, ADD_VALUES));
66 row = 1;
67 col = 2;
68 PetscCall(MatSetValues(R, 1, &row, 1, &col, &one, ADD_VALUES));
69 row = 1;
70 col = 3;
71 PetscCall(MatSetValues(R, 1, &row, 1, &col, &one, ADD_VALUES));
72 PetscCall(MatAssemblyBegin(R, MAT_FINAL_ASSEMBLY));
73 PetscCall(MatAssemblyEnd(R, MAT_FINAL_ASSEMBLY));
74 PetscCall(MatSetOptionsPrefix(R, "R_"));
75 PetscCall(MatView(R, PETSC_VIEWER_STDOUT_WORLD));
76 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
77
78 /* C = A*R^T */
79 PetscCall(MatMatTransposeMult(A, R, MAT_INITIAL_MATRIX, 2.0, &C));
80 PetscCall(MatSetOptionsPrefix(C, "ARt_"));
81 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
82 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
83
84 /* Create MatTransposeColoring from symbolic C=A*R^T */
85 PetscCall(MatColoringCreate(C, &mc));
86 PetscCall(MatColoringSetDistance(mc, 2));
87 /* PetscCall(MatColoringSetType(mc,MATCOLORINGSL)); */
88 PetscCall(MatColoringSetFromOptions(mc));
89 PetscCall(MatColoringApply(mc, &iscoloring));
90 PetscCall(MatColoringDestroy(&mc));
91 PetscCall(MatTransposeColoringCreate(C, iscoloring, &matcoloring));
92 PetscCall(ISColoringDestroy(&iscoloring));
93
94 /* Create Rt_dense */
95 PetscCall(MatCreate(PETSC_COMM_WORLD, &Rt_dense));
96 PetscCall(MatSetSizes(Rt_dense, 4, matcoloring->ncolors, PETSC_DECIDE, PETSC_DECIDE));
97 PetscCall(MatSetType(Rt_dense, MATDENSE));
98 PetscCall(MatSeqDenseSetPreallocation(Rt_dense, NULL));
99 PetscCall(MatAssemblyBegin(Rt_dense, MAT_FINAL_ASSEMBLY));
100 PetscCall(MatAssemblyEnd(Rt_dense, MAT_FINAL_ASSEMBLY));
101 PetscCall(MatGetLocalSize(Rt_dense, &m, &n));
102 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Rt_dense: %" PetscInt_FMT ",%" PetscInt_FMT "\n", m, n));
103
104 /* Get Rt_dense by Apply MatTransposeColoring to R */
105 PetscCall(MatTransColoringApplySpToDen(matcoloring, R, Rt_dense));
106
107 /* C_dense = A*Rt_dense */
108 PetscCall(MatMatMult(A, Rt_dense, MAT_INITIAL_MATRIX, 2.0, &C_dense));
109 PetscCall(MatSetOptionsPrefix(C_dense, "ARt_dense_"));
110 /*PetscCall(MatView(C_dense,PETSC_VIEWER_STDOUT_WORLD)); */
111 /*PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n")); */
112
113 /* Recover C from C_dense */
114 PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &C_sparse));
115 PetscCall(MatTransColoringApplyDenToSp(matcoloring, C_dense, C_sparse));
116 PetscCall(MatSetOptionsPrefix(C_sparse, "ARt_color_"));
117 PetscCall(MatView(C_sparse, PETSC_VIEWER_STDOUT_WORLD));
118 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
119
120 PetscCall(MatDestroy(&C_dense));
121 PetscCall(MatDestroy(&C_sparse));
122 PetscCall(MatDestroy(&Rt_dense));
123 PetscCall(MatTransposeColoringDestroy(&matcoloring));
124 PetscCall(MatDestroy(&C));
125
126 /* Test PtAP = P^T*A*P, P = R^T */
127 PetscCall(MatTranspose(R, MAT_INITIAL_MATRIX, &P));
128 PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, 2.0, &PtAP));
129 PetscCall(MatSetOptionsPrefix(PtAP, "PtAP_"));
130 PetscCall(MatView(PtAP, PETSC_VIEWER_STDOUT_WORLD));
131 PetscCall(MatDestroy(&P));
132
133 /* Test C = RARt */
134 PetscCall(MatRARt(A, R, MAT_INITIAL_MATRIX, 2.0, &C));
135 PetscCall(MatRARt(A, R, MAT_REUSE_MATRIX, 2.0, &C));
136 PetscCall(MatEqual(C, PtAP, &equal));
137 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: PtAP != RARt");
138
139 /* Free spaces */
140 PetscCall(MatDestroy(&C));
141 PetscCall(MatDestroy(&A));
142 PetscCall(MatDestroy(&R));
143 PetscCall(MatDestroy(&PtAP));
144 PetscCall(PetscFinalize());
145 return 0;
146 }
147
148 /*TEST
149
150 test:
151 output_file: output/ex161.out
152 test:
153 suffix: 2
154 args: -matmattransmult_via color
155 output_file: output/ex161.out
156
157 test:
158 suffix: 3
159 args: -matmattransmult_via color -matden2sp_brows 3
160 output_file: output/ex161.out
161
162 test:
163 suffix: 4
164 args: -matmattransmult_via color -matrart_via r*art
165 output_file: output/ex161.out
166
167 test:
168 suffix: 5
169 args: -matmattransmult_via color -matrart_via coloring_rart
170 output_file: output/ex161.out
171
172 TEST*/
173