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 PetscErrorCode ierr; 15 16 /* 17 Create the Jacobian matrix 18 */ 19 M = 10; 20 N = 8; 21 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 22 ierr = MatSetType(A,MATAIJ);CHKERRQ(ierr); 23 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 24 ierr = MatSetBlockSize(A,1);CHKERRQ(ierr); 25 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 26 27 /* 28 Get the sizes of the jacobian. 29 */ 30 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 31 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 32 33 /* 34 Create a preallocator matrix with sizes (local and global) matching the jacobian A. 35 */ 36 ierr = MatCreate(PetscObjectComm((PetscObject)A),&preallocator);CHKERRQ(ierr); 37 ierr = MatSetType(preallocator,MATPREALLOCATOR);CHKERRQ(ierr); 38 ierr = MatSetSizes(preallocator,m,n,M,N);CHKERRQ(ierr); 39 ierr = MatSetBlockSize(preallocator,bs);CHKERRQ(ierr); 40 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 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 ierr = MatSetValues(preallocator,1,&ii,1,&jj,&vv,INSERT_VALUES);CHKERRQ(ierr); 52 53 ii = 7; jj = 4; 54 ierr = MatSetValues(preallocator,1,&ii,1,&jj,&vv,INSERT_VALUES);CHKERRQ(ierr); 55 56 ii = 9; jj = 7; 57 ierr = MatSetValue(preallocator,ii,jj,vv,INSERT_VALUES);CHKERRQ(ierr); 58 } 59 ierr = MatAssemblyBegin(preallocator,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 60 ierr = MatAssemblyEnd(preallocator,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 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 ierr = MatPreallocatorPreallocate(preallocator,PETSC_TRUE,A);CHKERRQ(ierr); 68 69 /* 70 We no longer require the preallocator object so destroy it. 71 */ 72 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 73 74 ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 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 ierr = MatSetValue(A,ii,jj,vv,INSERT_VALUES);CHKERRQ(ierr); 85 86 ii = 7; jj = 4; vv = 3.3; 87 ierr = MatSetValues(A,1,&ii,1,&jj,&vv,INSERT_VALUES);CHKERRQ(ierr); 88 89 ii = 9; jj = 7; vv = 4.3; 90 ierr = MatSetValues(A,1,&ii,1,&jj,&vv,INSERT_VALUES);CHKERRQ(ierr); 91 } 92 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 93 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 94 95 ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 96 97 ierr = MatDestroy(&A);CHKERRQ(ierr); 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 PetscErrorCode ierr; 106 107 ierr = PetscOptionsGetInt(NULL,NULL,"-block_size",&bs,NULL);CHKERRQ(ierr); 108 109 /* 110 Create the Jacobian matrix. 111 */ 112 M = 10 * bs; 113 N = 10 * bs; 114 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 115 ierr = MatSetType(A,MATAIJ);CHKERRQ(ierr); 116 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 117 ierr = MatSetBlockSize(A,bs);CHKERRQ(ierr); 118 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 119 120 /* 121 Get the sizes of the jacobian. 122 */ 123 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 124 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 125 126 /* 127 Create a preallocator matrix with dimensions matching the jacobian A. 128 */ 129 ierr = MatCreate(PetscObjectComm((PetscObject)A),&preallocator);CHKERRQ(ierr); 130 ierr = MatSetType(preallocator,MATPREALLOCATOR);CHKERRQ(ierr); 131 ierr = MatSetSizes(preallocator,m,n,M,N);CHKERRQ(ierr); 132 ierr = MatSetBlockSize(preallocator,bs);CHKERRQ(ierr); 133 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 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 ierr = PetscCalloc1(bs*bs,&vv);CHKERRQ(ierr); 144 145 ii = 0; jj = 9; 146 ierr = MatSetValue(preallocator,ii,jj,vv[0],INSERT_VALUES);CHKERRQ(ierr); 147 148 ii = 0; jj = 0; 149 ierr = MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);CHKERRQ(ierr); 150 151 ii = 2; jj = 4; 152 ierr = MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);CHKERRQ(ierr); 153 154 ii = 4; jj = 2; 155 ierr = MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);CHKERRQ(ierr); 156 157 ii = 4; jj = 4; 158 ierr = MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);CHKERRQ(ierr); 159 160 ii = 9; jj = 9; 161 ierr = MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);CHKERRQ(ierr); 162 163 ierr = PetscFree(vv);CHKERRQ(ierr); 164 } 165 ierr = MatAssemblyBegin(preallocator,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 166 ierr = MatAssemblyEnd(preallocator,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 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 ierr = MatPreallocatorPreallocate(preallocator,PETSC_TRUE,A);CHKERRQ(ierr); 174 175 /* 176 We no longer require the preallocator object so destroy it. 177 */ 178 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 179 180 ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 181 182 { 183 PetscInt ii,jj; 184 PetscScalar *vv; 185 186 ierr = PetscCalloc1(bs*bs,&vv);CHKERRQ(ierr); 187 for (ii=0; ii<bs*bs; ii++) { 188 vv[ii] = (PetscReal)(ii+1); 189 } 190 191 ii = 0; jj = 9; 192 ierr = MatSetValue(A,ii,jj,33.3,INSERT_VALUES);CHKERRQ(ierr); 193 194 ii = 0; jj = 0; 195 ierr = MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);CHKERRQ(ierr); 196 197 ii = 2; jj = 4; 198 ierr = MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);CHKERRQ(ierr); 199 200 ii = 4; jj = 2; 201 ierr = MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);CHKERRQ(ierr); 202 203 ii = 4; jj = 4; 204 ierr = MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);CHKERRQ(ierr); 205 206 ii = 9; jj = 9; 207 ierr = MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);CHKERRQ(ierr); 208 209 ierr = PetscFree(vv);CHKERRQ(ierr); 210 } 211 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 212 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 213 214 ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 215 216 ierr = MatDestroy(&A);CHKERRQ(ierr); 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 ierr = PetscOptionsGetInt(NULL,NULL,"-test_id",&testid,NULL);CHKERRQ(ierr); 226 switch (testid) { 227 case 0: 228 ierr = ex1_nonsquare_bs1();CHKERRQ(ierr); 229 break; 230 case 1: 231 ierr = ex2_square_bsvariable();CHKERRQ(ierr); 232 break; 233 default: 234 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Invalid value for -test_id. Must be {0,1}"); 235 break; 236 } 237 ierr = PetscFinalize(); 238 return ierr; 239 } 240 241 /*TEST 242 243 test: 244 suffix: t0_a_aij 245 nsize: 1 246 args: -test_id 0 -mat_type aij 247 248 test: 249 suffix: t0_b_aij 250 nsize: 6 251 args: -test_id 0 -mat_type aij 252 253 test: 254 suffix: t1_a_aij 255 nsize: 1 256 args: -test_id 1 -mat_type aij 257 258 test: 259 suffix: t2_a_baij 260 nsize: 1 261 args: -test_id 1 -mat_type baij 262 263 test: 264 suffix: t3_a_sbaij 265 nsize: 1 266 args: -test_id 1 -mat_type sbaij 267 268 test: 269 suffix: t4_a_aij_bs3 270 nsize: 1 271 args: -test_id 1 -mat_type aij -block_size 3 272 273 test: 274 suffix: t5_a_baij_bs3 275 nsize: 1 276 args: -test_id 1 -mat_type baij -block_size 3 277 278 test: 279 suffix: t6_a_sbaij_bs3 280 nsize: 1 281 args: -test_id 1 -mat_type sbaij -block_size 3 282 283 test: 284 suffix: t4_b_aij_bs3 285 nsize: 6 286 args: -test_id 1 -mat_type aij -block_size 3 287 288 test: 289 suffix: t5_b_baij_bs3 290 nsize: 6 291 args: -test_id 1 -mat_type baij -block_size 3 292 filter: grep -v Mat_ 293 294 test: 295 suffix: t6_b_sbaij_bs3 296 nsize: 6 297 args: -test_id 1 -mat_type sbaij -block_size 3 298 filter: grep -v Mat_ 299 300 TEST*/ 301