Actual source code: ex161.c
1: static char help[] = "Test MatTransposeColoring for SeqAIJ matrices. Used for '-matmattransmult_color' on MatMatTransposeMult \n\n";
3: #include <petscmat.h>
4: #include <petsc/private/matimpl.h>
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;
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!");
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"));
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));
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"));
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"));
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));
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));
104: /* Get Rt_dense by Apply MatTransposeColoring to R */
105: PetscCall(MatTransColoringApplySpToDen(matcoloring, R, Rt_dense));
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")); */
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"));
120: PetscCall(MatDestroy(&C_dense));
121: PetscCall(MatDestroy(&C_sparse));
122: PetscCall(MatDestroy(&Rt_dense));
123: PetscCall(MatTransposeColoringDestroy(&matcoloring));
124: PetscCall(MatDestroy(&C));
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));
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");
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: }
148: /*TEST
150: test:
151: output_file: output/ex161.out
152: test:
153: suffix: 2
154: args: -matmattransmult_via color
155: output_file: output/ex161.out
157: test:
158: suffix: 3
159: args: -matmattransmult_via color -matden2sp_brows 3
160: output_file: output/ex161.out
162: test:
163: suffix: 4
164: args: -matmattransmult_via color -matrart_via r*art
165: output_file: output/ex161.out
167: test:
168: suffix: 5
169: args: -matmattransmult_via color -matrart_via coloring_rart
170: output_file: output/ex161.out
172: TEST*/