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