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 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, (char *)0, 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