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