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;

 18:   PetscInitialize(&argc, &argv, (char *)0, help);
 19:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

 22:   /* Create seqaij A */
 23:   MatCreate(PETSC_COMM_SELF, &A);
 24:   MatSetSizes(A, 4, 4, 4, 4);
 25:   MatSetType(A, MATSEQAIJ);
 26:   MatSetFromOptions(A);
 27:   MatSetUp(A);
 28:   row = 0;
 29:   col = 0;
 30:   val = 1.0;
 31:   MatSetValues(A, 1, &row, 1, &col, &val, ADD_VALUES);
 32:   row = 1;
 33:   col = 3;
 34:   val = 2.0;
 35:   MatSetValues(A, 1, &row, 1, &col, &val, ADD_VALUES);
 36:   row = 2;
 37:   col = 2;
 38:   val = 3.0;
 39:   MatSetValues(A, 1, &row, 1, &col, &val, ADD_VALUES);
 40:   row = 3;
 41:   col = 0;
 42:   val = 4.0;
 43:   MatSetValues(A, 1, &row, 1, &col, &val, ADD_VALUES);
 44:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 45:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 46:   MatSetOptionsPrefix(A, "A_");
 47:   MatView(A, PETSC_VIEWER_STDOUT_WORLD);
 48:   PetscPrintf(PETSC_COMM_SELF, "\n");

 50:   /* Create seqaij R */
 51:   MatCreate(PETSC_COMM_SELF, &R);
 52:   MatSetSizes(R, 2, 4, 2, 4);
 53:   MatSetType(R, MATSEQAIJ);
 54:   MatSetFromOptions(R);
 55:   MatSetUp(R);
 56:   row = 0;
 57:   col = 0;
 58:   MatSetValues(R, 1, &row, 1, &col, &one, ADD_VALUES);
 59:   row = 0;
 60:   col = 1;
 61:   MatSetValues(R, 1, &row, 1, &col, &one, ADD_VALUES);

 63:   row = 1;
 64:   col = 1;
 65:   MatSetValues(R, 1, &row, 1, &col, &one, ADD_VALUES);
 66:   row = 1;
 67:   col = 2;
 68:   MatSetValues(R, 1, &row, 1, &col, &one, ADD_VALUES);
 69:   row = 1;
 70:   col = 3;
 71:   MatSetValues(R, 1, &row, 1, &col, &one, ADD_VALUES);
 72:   MatAssemblyBegin(R, MAT_FINAL_ASSEMBLY);
 73:   MatAssemblyEnd(R, MAT_FINAL_ASSEMBLY);
 74:   MatSetOptionsPrefix(R, "R_");
 75:   MatView(R, PETSC_VIEWER_STDOUT_WORLD);
 76:   PetscPrintf(PETSC_COMM_SELF, "\n");

 78:   /* C = A*R^T */
 79:   MatMatTransposeMult(A, R, MAT_INITIAL_MATRIX, 2.0, &C);
 80:   MatSetOptionsPrefix(C, "ARt_");
 81:   MatView(C, PETSC_VIEWER_STDOUT_WORLD);
 82:   PetscPrintf(PETSC_COMM_SELF, "\n");

 84:   /* Create MatTransposeColoring from symbolic C=A*R^T */
 85:   MatColoringCreate(C, &mc);
 86:   MatColoringSetDistance(mc, 2);
 87:   /* MatColoringSetType(mc,MATCOLORINGSL); */
 88:   MatColoringSetFromOptions(mc);
 89:   MatColoringApply(mc, &iscoloring);
 90:   MatColoringDestroy(&mc);
 91:   MatTransposeColoringCreate(C, iscoloring, &matcoloring);
 92:   ISColoringDestroy(&iscoloring);

 94:   /* Create Rt_dense */
 95:   MatCreate(PETSC_COMM_WORLD, &Rt_dense);
 96:   MatSetSizes(Rt_dense, 4, matcoloring->ncolors, PETSC_DECIDE, PETSC_DECIDE);
 97:   MatSetType(Rt_dense, MATDENSE);
 98:   MatSeqDenseSetPreallocation(Rt_dense, NULL);
 99:   MatAssemblyBegin(Rt_dense, MAT_FINAL_ASSEMBLY);
100:   MatAssemblyEnd(Rt_dense, MAT_FINAL_ASSEMBLY);
101:   MatGetLocalSize(Rt_dense, &m, &n);
102:   PetscPrintf(PETSC_COMM_SELF, "Rt_dense: %" PetscInt_FMT ",%" PetscInt_FMT "\n", m, n);

104:   /* Get Rt_dense by Apply MatTransposeColoring to R */
105:   MatTransColoringApplySpToDen(matcoloring, R, Rt_dense);

107:   /* C_dense = A*Rt_dense */
108:   MatMatMult(A, Rt_dense, MAT_INITIAL_MATRIX, 2.0, &C_dense);
109:   MatSetOptionsPrefix(C_dense, "ARt_dense_");
110:   /*MatView(C_dense,PETSC_VIEWER_STDOUT_WORLD); */
111:   /*PetscPrintf(PETSC_COMM_SELF,"\n"); */

113:   /* Recover C from C_dense */
114:   MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &C_sparse);
115:   MatTransColoringApplyDenToSp(matcoloring, C_dense, C_sparse);
116:   MatSetOptionsPrefix(C_sparse, "ARt_color_");
117:   MatView(C_sparse, PETSC_VIEWER_STDOUT_WORLD);
118:   PetscPrintf(PETSC_COMM_SELF, "\n");

120:   MatDestroy(&C_dense);
121:   MatDestroy(&C_sparse);
122:   MatDestroy(&Rt_dense);
123:   MatTransposeColoringDestroy(&matcoloring);
124:   MatDestroy(&C);

126:   /* Test PtAP = P^T*A*P, P = R^T */
127:   MatTranspose(R, MAT_INITIAL_MATRIX, &P);
128:   MatPtAP(A, P, MAT_INITIAL_MATRIX, 2.0, &PtAP);
129:   MatSetOptionsPrefix(PtAP, "PtAP_");
130:   MatView(PtAP, PETSC_VIEWER_STDOUT_WORLD);
131:   MatDestroy(&P);

133:   /* Test C = RARt */
134:   MatRARt(A, R, MAT_INITIAL_MATRIX, 2.0, &C);
135:   MatRARt(A, R, MAT_REUSE_MATRIX, 2.0, &C);
136:   MatEqual(C, PtAP, &equal);

139:   /* Free spaces */
140:   MatDestroy(&C);
141:   MatDestroy(&A);
142:   MatDestroy(&R);
143:   MatDestroy(&PtAP);
144:   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*/