1 static char help[] = "Example of using MatPreallocator\n\n"; 2 3 #include <petscmat.h> 4 5 PetscErrorCode ex1_nonsquare_bs1(void) 6 { 7 Mat A, preallocator; 8 PetscInt M, N, m, n, bs; 9 10 /* 11 Create the Jacobian matrix 12 */ 13 PetscFunctionBegin; 14 M = 10; 15 N = 8; 16 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 17 PetscCall(MatSetType(A, MATAIJ)); 18 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N)); 19 PetscCall(MatSetBlockSize(A, 1)); 20 PetscCall(MatSetFromOptions(A)); 21 22 /* 23 Get the sizes of the jacobian. 24 */ 25 PetscCall(MatGetLocalSize(A, &m, &n)); 26 PetscCall(MatGetBlockSize(A, &bs)); 27 28 /* 29 Create a preallocator matrix with sizes (local and global) matching the jacobian A. 30 */ 31 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &preallocator)); 32 PetscCall(MatSetType(preallocator, MATPREALLOCATOR)); 33 PetscCall(MatSetSizes(preallocator, m, n, M, N)); 34 PetscCall(MatSetBlockSize(preallocator, bs)); 35 PetscCall(MatSetUp(preallocator)); 36 37 /* 38 Insert non-zero pattern (e.g. perform a sweep over the grid). 39 You can use MatSetValues(), MatSetValuesBlocked() or MatSetValue(). 40 */ 41 { 42 PetscInt ii, jj; 43 PetscScalar vv = 0.0; 44 45 ii = 3; 46 jj = 3; 47 PetscCall(MatSetValues(preallocator, 1, &ii, 1, &jj, &vv, INSERT_VALUES)); 48 49 ii = 7; 50 jj = 4; 51 PetscCall(MatSetValues(preallocator, 1, &ii, 1, &jj, &vv, INSERT_VALUES)); 52 53 ii = 9; 54 jj = 7; 55 PetscCall(MatSetValue(preallocator, ii, jj, vv, INSERT_VALUES)); 56 } 57 PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY)); 58 PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY)); 59 60 /* 61 Push the non-zero pattern defined within preallocator into A. 62 Internally this will call MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE). 63 The arg fill = PETSC_TRUE will insert zeros in the matrix A automatically. 64 */ 65 PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_TRUE, A)); 66 67 /* 68 We no longer require the preallocator object so destroy it. 69 */ 70 PetscCall(MatDestroy(&preallocator)); 71 72 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 73 74 /* 75 Insert non-zero values into A. 76 */ 77 { 78 PetscInt ii, jj; 79 PetscScalar vv; 80 81 ii = 3; 82 jj = 3; 83 vv = 0.3; 84 PetscCall(MatSetValue(A, ii, jj, vv, INSERT_VALUES)); 85 86 ii = 7; 87 jj = 4; 88 vv = 3.3; 89 PetscCall(MatSetValues(A, 1, &ii, 1, &jj, &vv, INSERT_VALUES)); 90 91 ii = 9; 92 jj = 7; 93 vv = 4.3; 94 PetscCall(MatSetValues(A, 1, &ii, 1, &jj, &vv, INSERT_VALUES)); 95 } 96 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 97 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 98 99 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 100 101 PetscCall(MatDestroy(&A)); 102 PetscFunctionReturn(0); 103 } 104 105 PetscErrorCode ex2_square_bsvariable(void) 106 { 107 Mat A, preallocator; 108 PetscInt M, N, m, n, bs = 1; 109 110 PetscFunctionBegin; 111 PetscCall(PetscOptionsGetInt(NULL, NULL, "-block_size", &bs, NULL)); 112 113 /* 114 Create the Jacobian matrix. 115 */ 116 M = 10 * bs; 117 N = 10 * bs; 118 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 119 PetscCall(MatSetType(A, MATAIJ)); 120 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N)); 121 PetscCall(MatSetBlockSize(A, bs)); 122 PetscCall(MatSetFromOptions(A)); 123 124 /* 125 Get the sizes of the jacobian. 126 */ 127 PetscCall(MatGetLocalSize(A, &m, &n)); 128 PetscCall(MatGetBlockSize(A, &bs)); 129 130 /* 131 Create a preallocator matrix with dimensions matching the jacobian A. 132 */ 133 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &preallocator)); 134 PetscCall(MatSetType(preallocator, MATPREALLOCATOR)); 135 PetscCall(MatSetSizes(preallocator, m, n, M, N)); 136 PetscCall(MatSetBlockSize(preallocator, bs)); 137 PetscCall(MatSetUp(preallocator)); 138 139 /* 140 Insert non-zero pattern (e.g. perform a sweep over the grid). 141 You can use MatSetValues(), MatSetValuesBlocked() or MatSetValue(). 142 */ 143 { 144 PetscInt ii, jj; 145 PetscScalar *vv; 146 147 PetscCall(PetscCalloc1(bs * bs, &vv)); 148 149 ii = 0; 150 jj = 9; 151 PetscCall(MatSetValue(preallocator, ii, jj, vv[0], INSERT_VALUES)); 152 153 ii = 0; 154 jj = 0; 155 PetscCall(MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES)); 156 157 ii = 2; 158 jj = 4; 159 PetscCall(MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES)); 160 161 ii = 4; 162 jj = 2; 163 PetscCall(MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES)); 164 165 ii = 4; 166 jj = 4; 167 PetscCall(MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES)); 168 169 ii = 9; 170 jj = 9; 171 PetscCall(MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES)); 172 173 PetscCall(PetscFree(vv)); 174 } 175 PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY)); 176 PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY)); 177 178 /* 179 Push non-zero pattern defined from preallocator into A. 180 Internally this will call MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE). 181 The arg fill = PETSC_TRUE will insert zeros in the matrix A automatically. 182 */ 183 PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_TRUE, A)); 184 185 /* 186 We no longer require the preallocator object so destroy it. 187 */ 188 PetscCall(MatDestroy(&preallocator)); 189 190 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 191 192 { 193 PetscInt ii, jj; 194 PetscScalar *vv; 195 196 PetscCall(PetscCalloc1(bs * bs, &vv)); 197 for (ii = 0; ii < bs * bs; ii++) vv[ii] = (PetscReal)(ii + 1); 198 199 ii = 0; 200 jj = 9; 201 PetscCall(MatSetValue(A, ii, jj, 33.3, INSERT_VALUES)); 202 203 ii = 0; 204 jj = 0; 205 PetscCall(MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES)); 206 207 ii = 2; 208 jj = 4; 209 PetscCall(MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES)); 210 211 ii = 4; 212 jj = 2; 213 PetscCall(MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES)); 214 215 ii = 4; 216 jj = 4; 217 PetscCall(MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES)); 218 219 ii = 9; 220 jj = 9; 221 PetscCall(MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES)); 222 223 PetscCall(PetscFree(vv)); 224 } 225 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 226 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 227 228 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 229 230 PetscCall(MatDestroy(&A)); 231 PetscFunctionReturn(0); 232 } 233 234 int main(int argc, char **args) 235 { 236 PetscInt testid = 0; 237 PetscFunctionBeginUser; 238 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 239 PetscCall(PetscOptionsGetInt(NULL, NULL, "-test_id", &testid, NULL)); 240 switch (testid) { 241 case 0: 242 PetscCall(ex1_nonsquare_bs1()); 243 break; 244 case 1: 245 PetscCall(ex2_square_bsvariable()); 246 break; 247 default: 248 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Invalid value for -test_id. Must be {0,1}"); 249 } 250 PetscCall(PetscFinalize()); 251 return 0; 252 } 253 254 /*TEST 255 256 test: 257 suffix: t0_a_aij 258 nsize: 1 259 args: -test_id 0 -mat_type aij 260 261 test: 262 suffix: t0_b_aij 263 nsize: 6 264 args: -test_id 0 -mat_type aij 265 266 test: 267 suffix: t1_a_aij 268 nsize: 1 269 args: -test_id 1 -mat_type aij 270 271 test: 272 suffix: t2_a_baij 273 nsize: 1 274 args: -test_id 1 -mat_type baij 275 276 test: 277 suffix: t3_a_sbaij 278 nsize: 1 279 args: -test_id 1 -mat_type sbaij 280 281 test: 282 suffix: t4_a_aij_bs3 283 nsize: 1 284 args: -test_id 1 -mat_type aij -block_size 3 285 286 test: 287 suffix: t5_a_baij_bs3 288 nsize: 1 289 args: -test_id 1 -mat_type baij -block_size 3 290 291 test: 292 suffix: t6_a_sbaij_bs3 293 nsize: 1 294 args: -test_id 1 -mat_type sbaij -block_size 3 295 296 test: 297 suffix: t4_b_aij_bs3 298 nsize: 6 299 args: -test_id 1 -mat_type aij -block_size 3 300 301 test: 302 suffix: t5_b_baij_bs3 303 nsize: 6 304 args: -test_id 1 -mat_type baij -block_size 3 305 filter: grep -v Mat_ 306 307 test: 308 suffix: t6_b_sbaij_bs3 309 nsize: 6 310 args: -test_id 1 -mat_type sbaij -block_size 3 311 filter: grep -v Mat_ 312 313 TEST*/ 314