1 static char help[] = "Tests MatPtAP() for MPIMAIJ and MPIAIJ \n "; 2 3 #include <petscdmda.h> 4 5 int main(int argc, char **argv) 6 { 7 DM coarsedm, finedm; 8 PetscMPIInt size, rank; 9 PetscInt M, N, Z, i, nrows; 10 PetscScalar one = 1.0; 11 PetscReal fill = 2.0; 12 Mat A, P, C; 13 PetscScalar *array, alpha; 14 PetscBool Test_3D = PETSC_FALSE, flg; 15 const PetscInt *ia, *ja; 16 PetscInt dof; 17 MPI_Comm comm; 18 19 PetscFunctionBeginUser; 20 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 21 comm = PETSC_COMM_WORLD; 22 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 23 PetscCallMPI(MPI_Comm_size(comm, &size)); 24 M = 10; 25 N = 10; 26 Z = 10; 27 dof = 10; 28 29 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_3D", &Test_3D, NULL)); 30 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL)); 31 PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL)); 32 PetscCall(PetscOptionsGetInt(NULL, NULL, "-Z", &Z, NULL)); 33 /* Set up distributed array for fine grid */ 34 if (!Test_3D) { 35 PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, 1, NULL, NULL, &coarsedm)); 36 } else { 37 PetscCall(DMDACreate3d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, Z, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, 1, NULL, NULL, NULL, &coarsedm)); 38 } 39 PetscCall(DMSetFromOptions(coarsedm)); 40 PetscCall(DMSetUp(coarsedm)); 41 42 /* This makes sure the coarse DMDA has the same partition as the fine DMDA */ 43 PetscCall(DMRefine(coarsedm, PetscObjectComm((PetscObject)coarsedm), &finedm)); 44 45 /*------------------------------------------------------------*/ 46 PetscCall(DMSetMatType(finedm, MATAIJ)); 47 PetscCall(DMCreateMatrix(finedm, &A)); 48 49 /* set val=one to A */ 50 if (size == 1) { 51 PetscCall(MatGetRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg)); 52 if (flg) { 53 PetscCall(MatSeqAIJGetArray(A, &array)); 54 for (i = 0; i < ia[nrows]; i++) array[i] = one; 55 PetscCall(MatSeqAIJRestoreArray(A, &array)); 56 } 57 PetscCall(MatRestoreRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg)); 58 } else { 59 Mat AA, AB; 60 PetscCall(MatMPIAIJGetSeqAIJ(A, &AA, &AB, NULL)); 61 PetscCall(MatGetRowIJ(AA, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg)); 62 if (flg) { 63 PetscCall(MatSeqAIJGetArray(AA, &array)); 64 for (i = 0; i < ia[nrows]; i++) array[i] = one; 65 PetscCall(MatSeqAIJRestoreArray(AA, &array)); 66 } 67 PetscCall(MatRestoreRowIJ(AA, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg)); 68 PetscCall(MatGetRowIJ(AB, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg)); 69 if (flg) { 70 PetscCall(MatSeqAIJGetArray(AB, &array)); 71 for (i = 0; i < ia[nrows]; i++) array[i] = one; 72 PetscCall(MatSeqAIJRestoreArray(AB, &array)); 73 } 74 PetscCall(MatRestoreRowIJ(AB, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg)); 75 } 76 /* Create interpolation between the fine and coarse grids */ 77 PetscCall(DMCreateInterpolation(coarsedm, finedm, &P, NULL)); 78 79 /* Test P^T * A * P - MatPtAP() */ 80 /*------------------------------*/ 81 /* (1) Developer API */ 82 PetscCall(MatProductCreate(A, P, NULL, &C)); 83 PetscCall(MatProductSetType(C, MATPRODUCT_PtAP)); 84 PetscCall(MatProductSetAlgorithm(C, "allatonce")); 85 PetscCall(MatProductSetFill(C, PETSC_DEFAULT)); 86 PetscCall(MatProductSetFromOptions(C)); 87 PetscCall(MatProductSymbolic(C)); 88 PetscCall(MatProductNumeric(C)); 89 PetscCall(MatProductNumeric(C)); /* Test reuse of symbolic C */ 90 91 { /* Test MatProductView() */ 92 PetscViewer viewer; 93 PetscCall(PetscViewerASCIIOpen(comm, NULL, &viewer)); 94 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO)); 95 PetscCall(MatProductView(C, viewer)); 96 PetscCall(PetscViewerPopFormat(viewer)); 97 PetscCall(PetscViewerDestroy(&viewer)); 98 } 99 100 PetscCall(MatPtAPMultEqual(A, P, C, 10, &flg)); 101 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatProduct_PtAP"); 102 PetscCall(MatDestroy(&C)); 103 104 /* (2) User API */ 105 PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, fill, &C)); 106 /* Test MAT_REUSE_MATRIX - reuse symbolic C */ 107 alpha = 1.0; 108 for (i = 0; i < 1; i++) { 109 alpha -= 0.1; 110 PetscCall(MatScale(A, alpha)); 111 PetscCall(MatPtAP(A, P, MAT_REUSE_MATRIX, fill, &C)); 112 } 113 114 /* Free intermediate data structures created for reuse of C=Pt*A*P */ 115 PetscCall(MatProductClear(C)); 116 117 PetscCall(MatPtAPMultEqual(A, P, C, 10, &flg)); 118 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatPtAP"); 119 120 PetscCall(MatDestroy(&C)); 121 PetscCall(MatDestroy(&A)); 122 PetscCall(MatDestroy(&P)); 123 PetscCall(DMDestroy(&finedm)); 124 PetscCall(DMDestroy(&coarsedm)); 125 PetscCall(PetscFinalize()); 126 return 0; 127 } 128 129 /*TEST 130 131 test: 132 args: -M 10 -N 10 -Z 10 133 output_file: output/empty.out 134 135 test: 136 suffix: allatonce 137 nsize: 4 138 args: -M 10 -N 10 -Z 10 139 output_file: output/ex89_2.out 140 141 test: 142 suffix: allatonce_merged 143 nsize: 4 144 args: -M 10 -M 5 -M 10 -mat_product_algorithm allatonce_merged 145 output_file: output/ex89_3.out 146 147 test: 148 suffix: nonscalable_3D 149 nsize: 4 150 args: -M 10 -M 5 -M 10 -test_3D 1 -mat_product_algorithm nonscalable 151 output_file: output/ex89_4.out 152 153 test: 154 suffix: allatonce_merged_3D 155 nsize: 4 156 args: -M 10 -M 5 -M 10 -test_3D 1 -mat_product_algorithm allatonce_merged 157 output_file: output/ex89_3.out 158 159 test: 160 suffix: nonscalable 161 nsize: 4 162 args: -M 10 -N 10 -Z 10 -mat_product_algorithm nonscalable 163 output_file: output/ex89_5.out 164 165 TEST*/ 166