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