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*/