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