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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 232 } 233 234 int main(int argc, char **args) 235 { 236 PetscInt testid = 0; 237 238 PetscFunctionBeginUser; 239 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 240 PetscCall(PetscOptionsGetInt(NULL, NULL, "-test_id", &testid, NULL)); 241 switch (testid) { 242 case 0: 243 PetscCall(ex1_nonsquare_bs1()); 244 break; 245 case 1: 246 PetscCall(ex2_square_bsvariable()); 247 break; 248 default: 249 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Invalid value for -test_id. Must be {0,1}"); 250 } 251 PetscCall(PetscFinalize()); 252 return 0; 253 } 254 255 /*TEST 256 257 test: 258 suffix: t0_a_aij 259 nsize: 1 260 args: -test_id 0 -mat_type aij 261 262 test: 263 suffix: t0_b_aij 264 nsize: 6 265 args: -test_id 0 -mat_type aij 266 267 test: 268 suffix: t1_a_aij 269 nsize: 1 270 args: -test_id 1 -mat_type aij 271 272 test: 273 suffix: t2_a_baij 274 nsize: 1 275 args: -test_id 1 -mat_type baij 276 277 test: 278 suffix: t3_a_sbaij 279 nsize: 1 280 args: -test_id 1 -mat_type sbaij 281 282 test: 283 suffix: t4_a_aij_bs3 284 nsize: 1 285 args: -test_id 1 -mat_type aij -block_size 3 286 287 test: 288 suffix: t5_a_baij_bs3 289 nsize: 1 290 args: -test_id 1 -mat_type baij -block_size 3 291 292 test: 293 suffix: t6_a_sbaij_bs3 294 nsize: 1 295 args: -test_id 1 -mat_type sbaij -block_size 3 296 297 test: 298 suffix: t4_b_aij_bs3 299 nsize: 6 300 args: -test_id 1 -mat_type aij -block_size 3 301 302 test: 303 suffix: t5_b_baij_bs3 304 nsize: 6 305 args: -test_id 1 -mat_type baij -block_size 3 306 307 test: 308 suffix: t6_b_sbaij_bs3 309 nsize: 6 310 args: -test_id 1 -mat_type sbaij -block_size 3 311 312 TEST*/ 313