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