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