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; col=0; val=1.0; PetscCall(MatSetValues(A,1,&row,1,&col,&val,ADD_VALUES)); 29 row = 1; col=3; val=2.0; PetscCall(MatSetValues(A,1,&row,1,&col,&val,ADD_VALUES)); 30 row = 2; col=2; val=3.0; PetscCall(MatSetValues(A,1,&row,1,&col,&val,ADD_VALUES)); 31 row = 3; col=0; val=4.0; PetscCall(MatSetValues(A,1,&row,1,&col,&val,ADD_VALUES)); 32 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 33 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 34 PetscCall(MatSetOptionsPrefix(A,"A_")); 35 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_WORLD)); 36 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n")); 37 38 /* Create seqaij R */ 39 PetscCall(MatCreate(PETSC_COMM_SELF,&R)); 40 PetscCall(MatSetSizes(R,2,4,2,4)); 41 PetscCall(MatSetType(R,MATSEQAIJ)); 42 PetscCall(MatSetFromOptions(R)); 43 PetscCall(MatSetUp(R)); 44 row = 0; col=0; PetscCall(MatSetValues(R,1,&row,1,&col,&one,ADD_VALUES)); 45 row = 0; col=1; PetscCall(MatSetValues(R,1,&row,1,&col,&one,ADD_VALUES)); 46 47 row = 1; col=1; PetscCall(MatSetValues(R,1,&row,1,&col,&one,ADD_VALUES)); 48 row = 1; col=2; PetscCall(MatSetValues(R,1,&row,1,&col,&one,ADD_VALUES)); 49 row = 1; col=3; PetscCall(MatSetValues(R,1,&row,1,&col,&one,ADD_VALUES)); 50 PetscCall(MatAssemblyBegin(R,MAT_FINAL_ASSEMBLY)); 51 PetscCall(MatAssemblyEnd(R,MAT_FINAL_ASSEMBLY)); 52 PetscCall(MatSetOptionsPrefix(R,"R_")); 53 PetscCall(MatView(R,PETSC_VIEWER_STDOUT_WORLD)); 54 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n")); 55 56 /* C = A*R^T */ 57 PetscCall(MatMatTransposeMult(A,R,MAT_INITIAL_MATRIX,2.0,&C)); 58 PetscCall(MatSetOptionsPrefix(C,"ARt_")); 59 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 60 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n")); 61 62 /* Create MatTransposeColoring from symbolic C=A*R^T */ 63 PetscCall(MatColoringCreate(C,&mc)); 64 PetscCall(MatColoringSetDistance(mc,2)); 65 /* PetscCall(MatColoringSetType(mc,MATCOLORINGSL)); */ 66 PetscCall(MatColoringSetFromOptions(mc)); 67 PetscCall(MatColoringApply(mc,&iscoloring)); 68 PetscCall(MatColoringDestroy(&mc)); 69 PetscCall(MatTransposeColoringCreate(C,iscoloring,&matcoloring)); 70 PetscCall(ISColoringDestroy(&iscoloring)); 71 72 /* Create Rt_dense */ 73 PetscCall(MatCreate(PETSC_COMM_WORLD,&Rt_dense)); 74 PetscCall(MatSetSizes(Rt_dense,4,matcoloring->ncolors,PETSC_DECIDE,PETSC_DECIDE)); 75 PetscCall(MatSetType(Rt_dense,MATDENSE)); 76 PetscCall(MatSeqDenseSetPreallocation(Rt_dense,NULL)); 77 PetscCall(MatAssemblyBegin(Rt_dense,MAT_FINAL_ASSEMBLY)); 78 PetscCall(MatAssemblyEnd(Rt_dense,MAT_FINAL_ASSEMBLY)); 79 PetscCall(MatGetLocalSize(Rt_dense,&m,&n)); 80 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Rt_dense: %" PetscInt_FMT ",%" PetscInt_FMT "\n",m,n)); 81 82 /* Get Rt_dense by Apply MatTransposeColoring to R */ 83 PetscCall(MatTransColoringApplySpToDen(matcoloring,R,Rt_dense)); 84 85 /* C_dense = A*Rt_dense */ 86 PetscCall(MatMatMult(A,Rt_dense,MAT_INITIAL_MATRIX,2.0,&C_dense)); 87 PetscCall(MatSetOptionsPrefix(C_dense,"ARt_dense_")); 88 /*PetscCall(MatView(C_dense,PETSC_VIEWER_STDOUT_WORLD)); */ 89 /*PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n")); */ 90 91 /* Recover C from C_dense */ 92 PetscCall(MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&C_sparse)); 93 PetscCall(MatTransColoringApplyDenToSp(matcoloring,C_dense,C_sparse)); 94 PetscCall(MatSetOptionsPrefix(C_sparse,"ARt_color_")); 95 PetscCall(MatView(C_sparse,PETSC_VIEWER_STDOUT_WORLD)); 96 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n")); 97 98 PetscCall(MatDestroy(&C_dense)); 99 PetscCall(MatDestroy(&C_sparse)); 100 PetscCall(MatDestroy(&Rt_dense)); 101 PetscCall(MatTransposeColoringDestroy(&matcoloring)); 102 PetscCall(MatDestroy(&C)); 103 104 /* Test PtAP = P^T*A*P, P = R^T */ 105 PetscCall(MatTranspose(R,MAT_INITIAL_MATRIX,&P)); 106 PetscCall(MatPtAP(A,P,MAT_INITIAL_MATRIX,2.0,&PtAP)); 107 PetscCall(MatSetOptionsPrefix(PtAP,"PtAP_")); 108 PetscCall(MatView(PtAP,PETSC_VIEWER_STDOUT_WORLD)); 109 PetscCall(MatDestroy(&P)); 110 111 /* Test C = RARt */ 112 PetscCall(MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&C)); 113 PetscCall(MatRARt(A,R,MAT_REUSE_MATRIX,2.0,&C)); 114 PetscCall(MatEqual(C,PtAP,&equal)); 115 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: PtAP != RARt"); 116 117 /* Free spaces */ 118 PetscCall(MatDestroy(&C)); 119 PetscCall(MatDestroy(&A)); 120 PetscCall(MatDestroy(&R)); 121 PetscCall(MatDestroy(&PtAP)); 122 PetscCall(PetscFinalize()); 123 return 0; 124 } 125 126 /*TEST 127 128 test: 129 output_file: output/ex161.out 130 test: 131 suffix: 2 132 args: -matmattransmult_via color 133 output_file: output/ex161.out 134 135 test: 136 suffix: 3 137 args: -matmattransmult_via color -matden2sp_brows 3 138 output_file: output/ex161.out 139 140 test: 141 suffix: 4 142 args: -matmattransmult_via color -matrart_via r*art 143 output_file: output/ex161.out 144 145 test: 146 suffix: 5 147 args: -matmattransmult_via color -matrart_via coloring_rart 148 output_file: output/ex161.out 149 150 TEST*/ 151