1 static char help[] = "Benchmark for MatMatMult() of AIJ matrices using different 2d finite-difference stencils.\n\n"; 2 3 #include <petscmat.h> 4 5 /* Converts 3d grid coordinates (i,j,k) for a grid of size m \times n to global indexing. Pass k = 0 for a 2d grid. */ 6 int global_index(PetscInt i,PetscInt j,PetscInt k, PetscInt m, PetscInt n) { return i + j * m + k * m * n; } 7 8 int main(int argc,char **argv) 9 { 10 Mat A,B,C,PtAP,PtAP_copy,PtAP_squared; 11 PetscInt i,M,N,Istart,Iend,n=7,j,J,Ii,m=8,k,o=1; 12 PetscScalar v; 13 PetscBool equal=PETSC_FALSE,mat_view=PETSC_FALSE; 14 char stencil[PETSC_MAX_PATH_LEN]; 15 #if defined(PETSC_USE_LOG) 16 PetscLogStage fullMatMatMultStage; 17 #endif 18 19 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 20 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 21 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 22 PetscCall(PetscOptionsGetInt(NULL,NULL,"-o",&o,NULL)); 23 PetscCall(PetscOptionsHasName(NULL,NULL,"-result_view",&mat_view)); 24 PetscCall(PetscOptionsGetString(NULL,NULL,"-stencil",stencil,sizeof(stencil),NULL)); 25 26 /* Create a aij matrix A */ 27 M = N = m*n*o; 28 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 29 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N)); 30 PetscCall(MatSetType(A,MATAIJ)); 31 PetscCall(MatSetFromOptions(A)); 32 33 /* Consistency checks */ 34 PetscCheckFalse(o < 1 || m < 1 || n < 1,PETSC_COMM_WORLD,PETSC_ERR_USER,"Dimensions need to be larger than zero!"); 35 36 /************ 2D stencils ***************/ 37 PetscCall(PetscStrcmp(stencil, "2d5point", &equal)); 38 if (equal) { /* 5-point stencil, 2D */ 39 PetscCall(MatMPIAIJSetPreallocation(A,5,NULL,5,NULL)); 40 PetscCall(MatSeqAIJSetPreallocation(A,5,NULL)); 41 PetscCall(MatGetOwnershipRange(A,&Istart,&Iend)); 42 for (Ii=Istart; Ii<Iend; Ii++) { 43 v = -1.0; k = Ii / (m*n); j = (Ii - k * m * n) / m; i = (Ii - k * m * n - j * m); 44 if (i>0) {J = global_index(i-1,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 45 if (i<m-1) {J = global_index(i+1,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 46 if (j>0) {J = global_index(i,j-1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 47 if (j<n-1) {J = global_index(i,j+1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 48 v = 4.0; PetscCall(MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES)); 49 } 50 } 51 PetscCall(PetscStrcmp(stencil, "2d9point", &equal)); 52 if (equal) { /* 9-point stencil, 2D */ 53 PetscCall(MatMPIAIJSetPreallocation(A,9,NULL,9,NULL)); 54 PetscCall(MatSeqAIJSetPreallocation(A,9,NULL)); 55 PetscCall(MatGetOwnershipRange(A,&Istart,&Iend)); 56 for (Ii=Istart; Ii<Iend; Ii++) { 57 v = -1.0; k = Ii / (m*n); j = (Ii - k * m * n) / m; i = (Ii - k * m * n - j * m); 58 if (i>0) {J = global_index(i-1,j, k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 59 if (i>0 && j>0) {J = global_index(i-1,j-1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 60 if (j>0) {J = global_index(i, j-1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 61 if (i<m-1 && j>0) {J = global_index(i+1,j-1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 62 if (i<m-1) {J = global_index(i+1,j, k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 63 if (i<m-1 && j<n-1) {J = global_index(i+1,j+1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 64 if (j<n-1) {J = global_index(i, j+1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 65 if (i>0 && j<n-1) {J = global_index(i-1,j+1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 66 v = 8.0; PetscCall(MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES)); 67 } 68 } 69 PetscCall(PetscStrcmp(stencil, "2d9point2", &equal)); 70 if (equal) { /* 9-point Cartesian stencil (width 2 per coordinate), 2D */ 71 PetscCall(MatMPIAIJSetPreallocation(A,9,NULL,9,NULL)); 72 PetscCall(MatSeqAIJSetPreallocation(A,9,NULL)); 73 PetscCall(MatGetOwnershipRange(A,&Istart,&Iend)); 74 for (Ii=Istart; Ii<Iend; Ii++) { 75 v = -1.0; k = Ii / (m*n); j = (Ii - k * m * n) / m; i = (Ii - k * m * n - j * m); 76 if (i>0) {J = global_index(i-1,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 77 if (i>1) {J = global_index(i-2,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 78 if (i<m-1) {J = global_index(i+1,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 79 if (i<m-2) {J = global_index(i+2,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 80 if (j>0) {J = global_index(i,j-1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 81 if (j>1) {J = global_index(i,j-2,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 82 if (j<n-1) {J = global_index(i,j+1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 83 if (j<n-2) {J = global_index(i,j+2,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 84 v = 8.0; PetscCall(MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES)); 85 } 86 } 87 PetscCall(PetscStrcmp(stencil, "2d13point", &equal)); 88 if (equal) { /* 13-point Cartesian stencil (width 3 per coordinate), 2D */ 89 PetscCall(MatMPIAIJSetPreallocation(A,13,NULL,13,NULL)); 90 PetscCall(MatSeqAIJSetPreallocation(A,13,NULL)); 91 PetscCall(MatGetOwnershipRange(A,&Istart,&Iend)); 92 for (Ii=Istart; Ii<Iend; Ii++) { 93 v = -1.0; k = Ii / (m*n); j = (Ii - k * m * n) / m; i = (Ii - k * m * n - j * m); 94 if (i>0) {J = global_index(i-1,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 95 if (i>1) {J = global_index(i-2,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 96 if (i>2) {J = global_index(i-3,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 97 if (i<m-1) {J = global_index(i+1,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 98 if (i<m-2) {J = global_index(i+2,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 99 if (i<m-3) {J = global_index(i+3,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 100 if (j>0) {J = global_index(i,j-1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 101 if (j>1) {J = global_index(i,j-2,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 102 if (j>2) {J = global_index(i,j-3,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 103 if (j<n-1) {J = global_index(i,j+1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 104 if (j<n-2) {J = global_index(i,j+2,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 105 if (j<n-3) {J = global_index(i,j+3,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 106 v = 12.0; PetscCall(MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES)); 107 } 108 } 109 /************ 3D stencils ***************/ 110 PetscCall(PetscStrcmp(stencil, "3d7point", &equal)); 111 if (equal) { /* 7-point stencil, 3D */ 112 PetscCall(MatMPIAIJSetPreallocation(A,7,NULL,7,NULL)); 113 PetscCall(MatSeqAIJSetPreallocation(A,7,NULL)); 114 PetscCall(MatGetOwnershipRange(A,&Istart,&Iend)); 115 for (Ii=Istart; Ii<Iend; Ii++) { 116 v = -1.0; k = Ii / (m*n); j = (Ii - k * m * n) / m; i = (Ii - k * m * n - j * m); 117 if (i>0) {J = global_index(i-1,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 118 if (i<m-1) {J = global_index(i+1,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 119 if (j>0) {J = global_index(i,j-1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 120 if (j<n-1) {J = global_index(i,j+1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 121 if (k>0) {J = global_index(i,j,k-1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 122 if (k<o-1) {J = global_index(i,j,k+1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 123 v = 6.0; PetscCall(MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES)); 124 } 125 } 126 PetscCall(PetscStrcmp(stencil, "3d13point", &equal)); 127 if (equal) { /* 13-point stencil, 3D */ 128 PetscCall(MatMPIAIJSetPreallocation(A,13,NULL,13,NULL)); 129 PetscCall(MatSeqAIJSetPreallocation(A,13,NULL)); 130 PetscCall(MatGetOwnershipRange(A,&Istart,&Iend)); 131 for (Ii=Istart; Ii<Iend; Ii++) { 132 v = -1.0; k = Ii / (m*n); j = (Ii - k * m * n) / m; i = (Ii - k * m * n - j * m); 133 if (i>0) {J = global_index(i-1,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 134 if (i>1) {J = global_index(i-2,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 135 if (i<m-1) {J = global_index(i+1,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 136 if (i<m-2) {J = global_index(i+2,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 137 if (j>0) {J = global_index(i,j-1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 138 if (j>1) {J = global_index(i,j-2,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 139 if (j<n-1) {J = global_index(i,j+1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 140 if (j<n-2) {J = global_index(i,j+2,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 141 if (k>0) {J = global_index(i,j,k-1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 142 if (k>1) {J = global_index(i,j,k-2,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 143 if (k<o-1) {J = global_index(i,j,k+1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 144 if (k<o-2) {J = global_index(i,j,k+2,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 145 v = 12.0; PetscCall(MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES)); 146 } 147 } 148 PetscCall(PetscStrcmp(stencil, "3d19point", &equal)); 149 if (equal) { /* 19-point stencil, 3D */ 150 PetscCall(MatMPIAIJSetPreallocation(A,19,NULL,19,NULL)); 151 PetscCall(MatSeqAIJSetPreallocation(A,19,NULL)); 152 PetscCall(MatGetOwnershipRange(A,&Istart,&Iend)); 153 for (Ii=Istart; Ii<Iend; Ii++) { 154 v = -1.0; k = Ii / (m*n); j = (Ii - k * m * n) / m; i = (Ii - k * m * n - j * m); 155 /* one hop */ 156 if (i>0) {J = global_index(i-1,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 157 if (i<m-1) {J = global_index(i+1,j,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 158 if (j>0) {J = global_index(i,j-1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 159 if (j<n-1) {J = global_index(i,j+1,k,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 160 if (k>0) {J = global_index(i,j,k-1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 161 if (k<o-1) {J = global_index(i,j,k+1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 162 /* two hops */ 163 if (i>0 && j>0) {J = global_index(i-1,j-1,k ,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 164 if (i>0 && k>0) {J = global_index(i-1,j, k-1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 165 if (i>0 && j<n-1) {J = global_index(i-1,j+1,k ,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 166 if (i>0 && k<o-1) {J = global_index(i-1,j, k+1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 167 if (i<m-1 && j>0) {J = global_index(i+1,j-1,k ,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 168 if (i<m-1 && k>0) {J = global_index(i+1,j, k-1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 169 if (i<m-1 && j<n-1) {J = global_index(i+1,j+1,k ,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 170 if (i<m-1 && k<o-1) {J = global_index(i+1,j, k+1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 171 if (j>0 && k>0) {J = global_index(i, j-1,k-1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 172 if (j>0 && k<o-1) {J = global_index(i, j-1,k+1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 173 if (j<n-1 && k>0) {J = global_index(i, j+1,k-1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 174 if (j<n-1 && k<o-1) {J = global_index(i, j+1,k+1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 175 v = 18.0; PetscCall(MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES)); 176 } 177 } 178 PetscCall(PetscStrcmp(stencil, "3d27point", &equal)); 179 if (equal) { /* 27-point stencil, 3D */ 180 PetscCall(MatMPIAIJSetPreallocation(A,27,NULL,27,NULL)); 181 PetscCall(MatSeqAIJSetPreallocation(A,27,NULL)); 182 PetscCall(MatGetOwnershipRange(A,&Istart,&Iend)); 183 for (Ii=Istart; Ii<Iend; Ii++) { 184 v = -1.0; k = Ii / (m*n); j = (Ii - k * m * n) / m; i = (Ii - k * m * n - j * m); 185 if (k>0) { 186 if (j>0) { 187 if (i>0) {J = global_index(i-1,j-1,k-1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 188 J = global_index(i, j-1,k-1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES)); 189 if (i<m-1) {J = global_index(i+1,j-1,k-1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 190 } 191 { 192 if (i>0) {J = global_index(i-1,j, k-1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 193 J = global_index(i, j, k-1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES)); 194 if (i<m-1) {J = global_index(i+1,j, k-1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 195 } 196 if (j<n-1) { 197 if (i>0) {J = global_index(i-1,j+1,k-1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 198 J = global_index(i, j+1,k-1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES)); 199 if (i<m-1) {J = global_index(i+1,j+1,k-1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 200 } 201 } 202 { 203 if (j>0) { 204 if (i>0) {J = global_index(i-1,j-1,k ,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 205 J = global_index(i, j-1,k ,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES)); 206 if (i<m-1) {J = global_index(i+1,j-1,k ,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 207 } 208 { 209 if (i>0) {J = global_index(i-1,j, k ,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 210 J = global_index(i, j, k ,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES)); 211 if (i<m-1) {J = global_index(i+1,j, k ,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 212 } 213 if (j<n-1) { 214 if (i>0) {J = global_index(i-1,j+1,k ,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 215 J = global_index(i, j+1,k ,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES)); 216 if (i<m-1) {J = global_index(i+1,j+1,k ,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 217 } 218 } 219 if (k<o-1) { 220 if (j>0) { 221 if (i>0) {J = global_index(i-1,j-1,k+1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 222 J = global_index(i, j-1,k+1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES)); 223 if (i<m-1) {J = global_index(i+1,j-1,k+1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 224 } 225 { 226 if (i>0) {J = global_index(i-1,j, k+1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 227 J = global_index(i, j, k+1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES)); 228 if (i<m-1) {J = global_index(i+1,j, k+1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 229 } 230 if (j<n-1) { 231 if (i>0) {J = global_index(i-1,j+1,k+1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 232 J = global_index(i, j+1,k+1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES)); 233 if (i<m-1) {J = global_index(i+1,j+1,k+1,m,n); PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES));} 234 } 235 } 236 v = 26.0; PetscCall(MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES)); 237 } 238 } 239 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 240 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 241 242 /* Copy A into B in order to have a more representative benchmark (A*A has more cache hits than A*B) */ 243 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&B)); 244 245 PetscCall(PetscLogStageRegister("Full MatMatMult",&fullMatMatMultStage)); 246 247 /* Test C = A*B */ 248 PetscCall(PetscLogStagePush(fullMatMatMultStage)); 249 PetscCall(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C)); 250 251 /* Test PtAP_squared = PtAP(C,C)*PtAP(C,C) */ 252 PetscCall(MatPtAP(C,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&PtAP)); 253 PetscCall(MatDuplicate(PtAP,MAT_COPY_VALUES,&PtAP_copy)); 254 PetscCall(MatMatMult(PtAP,PtAP_copy,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&PtAP_squared)); 255 256 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 257 PetscCall(MatView(PtAP_squared,PETSC_VIEWER_STDOUT_WORLD)); 258 259 PetscCall(MatDestroy(&PtAP_squared)); 260 PetscCall(MatDestroy(&PtAP_copy)); 261 PetscCall(MatDestroy(&PtAP)); 262 PetscCall(MatDestroy(&C)); 263 PetscCall(MatDestroy(&B)); 264 PetscCall(MatDestroy(&A)); 265 PetscCall(PetscFinalize()); 266 return 0; 267 } 268 269 /*TEST 270 271 test: 272 suffix: 1 273 nsize: 1 274 args: -m 8 -n 8 -stencil 2d5point -matmatmult_via sorted 275 276 test: 277 suffix: 2 278 nsize: 1 279 args: -m 5 -n 5 -o 5 -stencil 3d27point -matmatmult_via rowmerge 280 281 test: 282 suffix: 3 283 nsize: 4 284 args: -m 6 -n 6 -stencil 2d5point -matmatmult_via seqmpi 285 286 TEST*/ 287