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