1 2 #include <petsc-private/matimpl.h> /*I "petscmat.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "MatCreate" 6 /*@ 7 MatCreate - Creates a matrix where the type is determined 8 from either a call to MatSetType() or from the options database 9 with a call to MatSetFromOptions(). The default matrix type is 10 AIJ, using the routines MatCreateSeqAIJ() or MatCreateAIJ() 11 if you do not set a type in the options database. If you never 12 call MatSetType() or MatSetFromOptions() it will generate an 13 error when you try to use the matrix. 14 15 Collective on MPI_Comm 16 17 Input Parameter: 18 . comm - MPI communicator 19 20 Output Parameter: 21 . A - the matrix 22 23 Options Database Keys: 24 + -mat_type seqaij - AIJ type, uses MatCreateSeqAIJ() 25 . -mat_type mpiaij - AIJ type, uses MatCreateAIJ() 26 . -mat_type seqdense - dense type, uses MatCreateSeqDense() 27 . -mat_type mpidense - dense type, uses MatCreateDense() 28 . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ() 29 - -mat_type mpibaij - block AIJ type, uses MatCreateBAIJ() 30 31 Even More Options Database Keys: 32 See the manpages for particular formats (e.g., MatCreateSeqAIJ()) 33 for additional format-specific options. 34 35 Notes: 36 37 Level: beginner 38 39 User manual sections: 40 + sec_matcreate 41 - chapter_matrices 42 43 .keywords: matrix, create 44 45 .seealso: MatCreateSeqAIJ(), MatCreateAIJ(), 46 MatCreateSeqDense(), MatCreateDense(), 47 MatCreateSeqBAIJ(), MatCreateBAIJ(), 48 MatCreateSeqSBAIJ(), MatCreateSBAIJ(), 49 MatConvert() 50 @*/ 51 PetscErrorCode MatCreate(MPI_Comm comm,Mat *A) 52 { 53 Mat B; 54 PetscErrorCode ierr; 55 56 PetscFunctionBegin; 57 PetscValidPointer(A,2); 58 59 *A = NULL; 60 ierr = MatInitializePackage();CHKERRQ(ierr); 61 62 ierr = PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_CLASSID,"Mat","Matrix","Mat",comm,MatDestroy,MatView);CHKERRQ(ierr); 63 ierr = PetscLayoutCreate(comm,&B->rmap);CHKERRQ(ierr); 64 ierr = PetscLayoutCreate(comm,&B->cmap);CHKERRQ(ierr); 65 66 B->preallocated = PETSC_FALSE; 67 *A = B; 68 PetscFunctionReturn(0); 69 } 70 71 #undef __FUNCT__ 72 #define __FUNCT__ "MatSetSizes" 73 /*@ 74 MatSetSizes - Sets the local and global sizes, and checks to determine compatibility 75 76 Collective on Mat 77 78 Input Parameters: 79 + A - the matrix 80 . m - number of local rows (or PETSC_DECIDE) 81 . n - number of local columns (or PETSC_DECIDE) 82 . M - number of global rows (or PETSC_DETERMINE) 83 - N - number of global columns (or PETSC_DETERMINE) 84 85 Notes: 86 m (n) and M (N) cannot be both PETSC_DECIDE 87 If one processor calls this with M (N) of PETSC_DECIDE then all processors must, otherwise the program will hang. 88 89 If PETSC_DECIDE is not used for the arguments 'm' and 'n', then the 90 user must ensure that they are chosen to be compatible with the 91 vectors. To do this, one first considers the matrix-vector product 92 'y = A x'. The 'm' that is used in the above routine must match the 93 local size used in the vector creation routine VecCreateMPI() for 'y'. 94 Likewise, the 'n' used must match that used as the local size in 95 VecCreateMPI() for 'x'. 96 97 You cannot change the sizes once they have been set. 98 99 The sizes must be set before MatSetUp() or MatXXXSetPreallocation() is called. 100 101 Level: beginner 102 103 .seealso: MatGetSize(), PetscSplitOwnership() 104 @*/ 105 PetscErrorCode MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N) 106 { 107 PetscFunctionBegin; 108 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 109 if (M > 0) PetscValidLogicalCollectiveInt(A,M,4); 110 if (N > 0) PetscValidLogicalCollectiveInt(A,N,5); 111 if (M > 0 && m > M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local column size %D cannot be larger than global column size %D",m,M); 112 if (N > 0 && n > N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local row size %D cannot be larger than global row size %D",n,N); 113 if ((A->rmap->n >= 0 || A->rmap->N >= 0) && (A->rmap->n != m || A->rmap->N != M)) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot change/reset row sizes to %D local %D global after previously setting them to %D local %D global",m,M,A->rmap->n,A->rmap->N); 114 if ((A->cmap->n >= 0 || A->cmap->N >= 0) && (A->cmap->n != n || A->cmap->N != N)) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot change/reset column sizes to %D local %D global after previously setting them to %D local %D global",n,N,A->cmap->n,A->cmap->N); 115 A->rmap->n = m; 116 A->cmap->n = n; 117 A->rmap->N = M; 118 A->cmap->N = N; 119 PetscFunctionReturn(0); 120 } 121 122 #undef __FUNCT__ 123 #define __FUNCT__ "MatSetFromOptions" 124 /*@ 125 MatSetFromOptions - Creates a matrix where the type is determined 126 from the options database. Generates a parallel MPI matrix if the 127 communicator has more than one processor. The default matrix type is 128 AIJ, using the routines MatCreateSeqAIJ() and MatCreateAIJ() if 129 you do not select a type in the options database. 130 131 Collective on Mat 132 133 Input Parameter: 134 . A - the matrix 135 136 Options Database Keys: 137 + -mat_type seqaij - AIJ type, uses MatCreateSeqAIJ() 138 . -mat_type mpiaij - AIJ type, uses MatCreateAIJ() 139 . -mat_type seqdense - dense type, uses MatCreateSeqDense() 140 . -mat_type mpidense - dense type, uses MatCreateDense() 141 . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ() 142 - -mat_type mpibaij - block AIJ type, uses MatCreateBAIJ() 143 144 Even More Options Database Keys: 145 See the manpages for particular formats (e.g., MatCreateSeqAIJ()) 146 for additional format-specific options. 147 148 Level: beginner 149 150 .keywords: matrix, create 151 152 .seealso: MatCreateSeqAIJ((), MatCreateAIJ(), 153 MatCreateSeqDense(), MatCreateDense(), 154 MatCreateSeqBAIJ(), MatCreateBAIJ(), 155 MatCreateSeqSBAIJ(), MatCreateSBAIJ(), 156 MatConvert() 157 @*/ 158 PetscErrorCode MatSetFromOptions(Mat B) 159 { 160 PetscErrorCode ierr; 161 const char *deft = MATAIJ; 162 char type[256]; 163 PetscBool flg,set; 164 165 PetscFunctionBegin; 166 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 167 168 ierr = PetscObjectOptionsBegin((PetscObject)B);CHKERRQ(ierr); 169 170 if (B->rmap->bs < 0) { 171 PetscInt newbs = -1; 172 ierr = PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSetBlockSize",newbs,&newbs,&flg);CHKERRQ(ierr); 173 if (flg) { 174 ierr = PetscLayoutSetBlockSize(B->rmap,newbs);CHKERRQ(ierr); 175 ierr = PetscLayoutSetBlockSize(B->cmap,newbs);CHKERRQ(ierr); 176 } 177 } 178 179 ierr = PetscOptionsFList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);CHKERRQ(ierr); 180 if (flg) { 181 ierr = MatSetType(B,type);CHKERRQ(ierr); 182 } else if (!((PetscObject)B)->type_name) { 183 ierr = MatSetType(B,deft);CHKERRQ(ierr); 184 } 185 186 ierr = PetscOptionsName("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",&B->checksymmetryonassembly);CHKERRQ(ierr); 187 ierr = PetscOptionsReal("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",0.0,&B->checksymmetrytol,NULL);CHKERRQ(ierr); 188 ierr = PetscOptionsBool("-mat_null_space_test","Checks if provided null space is correct in MatAssemblyEnd()","MatSetNullSpaceTest",PETSC_FALSE,&B->checknullspaceonassembly,NULL);CHKERRQ(ierr); 189 190 if (B->ops->setfromoptions) { 191 ierr = (*B->ops->setfromoptions)(B);CHKERRQ(ierr); 192 } 193 194 flg = PETSC_FALSE; 195 ierr = PetscOptionsBool("-mat_new_nonzero_location_err","Generate an error if new nonzeros are created in the matrix structure (useful to test preallocation)","MatSetOption",flg,&flg,&set);CHKERRQ(ierr); 196 if (set) {ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,flg);CHKERRQ(ierr);} 197 flg = PETSC_FALSE; 198 ierr = PetscOptionsBool("-mat_new_nonzero_allocation_err","Generate an error if new nonzeros are allocated in the matrix structure (useful to test preallocation)","MatSetOption",flg,&flg,&set);CHKERRQ(ierr); 199 if (set) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,flg);CHKERRQ(ierr);} 200 201 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 202 ierr = PetscObjectProcessOptionsHandlers((PetscObject)B);CHKERRQ(ierr); 203 ierr = PetscOptionsEnd();CHKERRQ(ierr); 204 PetscFunctionReturn(0); 205 } 206 207 #undef __FUNCT__ 208 #define __FUNCT__ "MatXAIJSetPreallocation" 209 /*@ 210 MatXAIJSetPreallocation - set preallocation for serial and parallel AIJ, BAIJ, and SBAIJ matrices 211 212 Collective on Mat 213 214 Input Arguments: 215 + A - matrix being preallocated 216 . bs - block size 217 . dnnz - number of nonzero blocks per block row of diagonal part of parallel matrix 218 . onnz - number of nonzero blocks per block row of off-diagonal part of parallel matrix 219 . dnnzu - number of nonzero blocks per block row of upper-triangular part of diagonal part of parallel matrix 220 - onnzu - number of nonzero blocks per block row of upper-triangular part of off-diagonal part of parallel matrix 221 222 Level: beginner 223 224 .seealso: MatSeqAIJSetPreallocation(), MatMPIAIJSetPreallocation(), MatSeqBAIJSetPreallocation(), MatMPIBAIJSetPreallocation(), MatSeqSBAIJSetPreallocation(), MatMPISBAIJSetPreallocation(), 225 PetscSplitOwnership() 226 @*/ 227 PetscErrorCode MatXAIJSetPreallocation(Mat A,PetscInt bs,const PetscInt dnnz[],const PetscInt onnz[],const PetscInt dnnzu[],const PetscInt onnzu[]) 228 { 229 PetscErrorCode ierr; 230 void (*aij)(void); 231 232 PetscFunctionBegin; 233 ierr = MatSetBlockSize(A,bs);CHKERRQ(ierr); 234 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 235 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 236 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 237 ierr = MatSeqBAIJSetPreallocation(A,bs,0,dnnz);CHKERRQ(ierr); 238 ierr = MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);CHKERRQ(ierr); 239 ierr = MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);CHKERRQ(ierr); 240 ierr = MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);CHKERRQ(ierr); 241 /* 242 In general, we have to do extra work to preallocate for scalar (AIJ) matrices so we check whether it will do any 243 good before going on with it. 244 */ 245 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 246 if (!aij) { 247 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 248 } 249 if (aij) { 250 if (bs == 1) { 251 ierr = MatSeqAIJSetPreallocation(A,0,dnnz);CHKERRQ(ierr); 252 ierr = MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);CHKERRQ(ierr); 253 } else { /* Convert block-row precallocation to scalar-row */ 254 PetscInt i,m,*sdnnz,*sonnz; 255 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 256 ierr = PetscMalloc2((!!dnnz)*m,&sdnnz,(!!onnz)*m,&sonnz);CHKERRQ(ierr); 257 for (i=0; i<m; i++) { 258 if (dnnz) sdnnz[i] = dnnz[i/bs] * bs; 259 if (onnz) sonnz[i] = onnz[i/bs] * bs; 260 } 261 ierr = MatSeqAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL);CHKERRQ(ierr); 262 ierr = MatMPIAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);CHKERRQ(ierr); 263 ierr = PetscFree2(sdnnz,sonnz);CHKERRQ(ierr); 264 } 265 } 266 PetscFunctionReturn(0); 267 } 268 269 /* 270 Merges some information from Cs header to A; the C object is then destroyed 271 272 This is somewhat different from MatHeaderReplace() it would be nice to merge the code 273 */ 274 #undef __FUNCT__ 275 #define __FUNCT__ "MatHeaderMerge" 276 PetscErrorCode MatHeaderMerge(Mat A,Mat C) 277 { 278 PetscErrorCode ierr; 279 PetscInt refct; 280 PetscOps *Abops; 281 MatOps Aops; 282 char *mtype,*mname; 283 void *spptr; 284 285 PetscFunctionBegin; 286 /* save the parts of A we need */ 287 Abops = ((PetscObject)A)->bops; 288 Aops = A->ops; 289 refct = ((PetscObject)A)->refct; 290 mtype = ((PetscObject)A)->type_name; 291 mname = ((PetscObject)A)->name; 292 spptr = A->spptr; 293 294 /* zero these so the destroy below does not free them */ 295 ((PetscObject)A)->type_name = 0; 296 ((PetscObject)A)->name = 0; 297 298 /* free all the interior data structures from mat */ 299 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 300 301 ierr = PetscFree(C->spptr);CHKERRQ(ierr); 302 303 ierr = PetscLayoutDestroy(&A->rmap);CHKERRQ(ierr); 304 ierr = PetscLayoutDestroy(&A->cmap);CHKERRQ(ierr); 305 ierr = PetscFunctionListDestroy(&((PetscObject)A)->qlist);CHKERRQ(ierr); 306 ierr = PetscObjectListDestroy(&((PetscObject)A)->olist);CHKERRQ(ierr); 307 308 /* copy C over to A */ 309 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 310 311 /* return the parts of A we saved */ 312 ((PetscObject)A)->bops = Abops; 313 A->ops = Aops; 314 ((PetscObject)A)->refct = refct; 315 ((PetscObject)A)->type_name = mtype; 316 ((PetscObject)A)->name = mname; 317 A->spptr = spptr; 318 319 /* since these two are copied into A we do not want them destroyed in C */ 320 ((PetscObject)C)->qlist = 0; 321 ((PetscObject)C)->olist = 0; 322 323 ierr = PetscHeaderDestroy(&C);CHKERRQ(ierr); 324 PetscFunctionReturn(0); 325 } 326 /* 327 Replace A's header with that of C; the C object is then destroyed 328 329 This is essentially code moved from MatDestroy() 330 331 This is somewhat different from MatHeaderMerge() it would be nice to merge the code 332 333 Used in DM hence is declared PETSC_EXTERN 334 */ 335 #undef __FUNCT__ 336 #define __FUNCT__ "MatHeaderReplace" 337 PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A,Mat C) 338 { 339 PetscErrorCode ierr; 340 PetscInt refct; 341 PetscObjectState state; 342 343 PetscFunctionBegin; 344 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 345 PetscValidHeaderSpecific(C,MAT_CLASSID,2); 346 if (A == C) PetscFunctionReturn(0); 347 PetscCheckSameComm(A,1,C,2); 348 if (((PetscObject)C)->refct != 1) SETERRQ1(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Object C has refct %D > 1, would leave hanging reference",((PetscObject)C)->refct); 349 350 /* free all the interior data structures from mat */ 351 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 352 ierr = PetscHeaderDestroy_Private((PetscObject)A);CHKERRQ(ierr); 353 ierr = PetscFree(A->ops);CHKERRQ(ierr); 354 ierr = PetscLayoutDestroy(&A->rmap);CHKERRQ(ierr); 355 ierr = PetscLayoutDestroy(&A->cmap);CHKERRQ(ierr); 356 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 357 358 /* copy C over to A */ 359 refct = ((PetscObject)A)->refct; 360 state = ((PetscObject)A)->state; 361 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 362 363 ((PetscObject)A)->refct = refct; 364 ((PetscObject)A)->state = state + 1; 365 366 ierr = PetscFree(C);CHKERRQ(ierr); 367 PetscFunctionReturn(0); 368 } 369