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 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 219 PetscCall(PetscOptionsGetInt(NULL,NULL,"-test_id",&testid,NULL)); 220 switch (testid) { 221 case 0: 222 PetscCall(ex1_nonsquare_bs1()); 223 break; 224 case 1: 225 PetscCall(ex2_square_bsvariable()); 226 break; 227 default: 228 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Invalid value for -test_id. Must be {0,1}"); 229 } 230 PetscCall(PetscFinalize()); 231 return 0; 232 } 233 234 /*TEST 235 236 test: 237 suffix: t0_a_aij 238 nsize: 1 239 args: -test_id 0 -mat_type aij 240 241 test: 242 suffix: t0_b_aij 243 nsize: 6 244 args: -test_id 0 -mat_type aij 245 246 test: 247 suffix: t1_a_aij 248 nsize: 1 249 args: -test_id 1 -mat_type aij 250 251 test: 252 suffix: t2_a_baij 253 nsize: 1 254 args: -test_id 1 -mat_type baij 255 256 test: 257 suffix: t3_a_sbaij 258 nsize: 1 259 args: -test_id 1 -mat_type sbaij 260 261 test: 262 suffix: t4_a_aij_bs3 263 nsize: 1 264 args: -test_id 1 -mat_type aij -block_size 3 265 266 test: 267 suffix: t5_a_baij_bs3 268 nsize: 1 269 args: -test_id 1 -mat_type baij -block_size 3 270 271 test: 272 suffix: t6_a_sbaij_bs3 273 nsize: 1 274 args: -test_id 1 -mat_type sbaij -block_size 3 275 276 test: 277 suffix: t4_b_aij_bs3 278 nsize: 6 279 args: -test_id 1 -mat_type aij -block_size 3 280 281 test: 282 suffix: t5_b_baij_bs3 283 nsize: 6 284 args: -test_id 1 -mat_type baij -block_size 3 285 filter: grep -v Mat_ 286 287 test: 288 suffix: t6_b_sbaij_bs3 289 nsize: 6 290 args: -test_id 1 -mat_type sbaij -block_size 3 291 filter: grep -v Mat_ 292 293 TEST*/ 294