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