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 = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 235 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 236 ierr = MatSeqBAIJSetPreallocation(A,bs,0,dnnz);CHKERRQ(ierr); 237 ierr = MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);CHKERRQ(ierr); 238 ierr = MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);CHKERRQ(ierr); 239 ierr = MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);CHKERRQ(ierr); 240 /* 241 In general, we have to do extra work to preallocate for scalar (AIJ) matrices so we check whether it will do any 242 good before going on with it. 243 */ 244 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 245 if (!aij) { 246 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 247 } 248 if (aij) { 249 if (bs == 1) { 250 ierr = MatSeqAIJSetPreallocation(A,0,dnnz);CHKERRQ(ierr); 251 ierr = MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);CHKERRQ(ierr); 252 } else { /* Convert block-row precallocation to scalar-row */ 253 PetscInt i,m,*sdnnz,*sonnz; 254 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 255 ierr = PetscMalloc2((!!dnnz)*m,&sdnnz,(!!onnz)*m,&sonnz);CHKERRQ(ierr); 256 for (i=0; i<m; i++) { 257 if (dnnz) sdnnz[i] = dnnz[i/bs] * bs; 258 if (onnz) sonnz[i] = onnz[i/bs] * bs; 259 } 260 ierr = MatSeqAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL);CHKERRQ(ierr); 261 ierr = MatMPIAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);CHKERRQ(ierr); 262 ierr = PetscFree2(sdnnz,sonnz);CHKERRQ(ierr); 263 } 264 } 265 PetscFunctionReturn(0); 266 } 267 268 /* 269 Merges some information from Cs header to A; the C object is then destroyed 270 271 This is somewhat different from MatHeaderReplace() it would be nice to merge the code 272 */ 273 #undef __FUNCT__ 274 #define __FUNCT__ "MatHeaderMerge" 275 PetscErrorCode MatHeaderMerge(Mat A,Mat C) 276 { 277 PetscErrorCode ierr; 278 PetscInt refct; 279 PetscOps *Abops; 280 MatOps Aops; 281 char *mtype,*mname; 282 void *spptr; 283 284 PetscFunctionBegin; 285 /* save the parts of A we need */ 286 Abops = ((PetscObject)A)->bops; 287 Aops = A->ops; 288 refct = ((PetscObject)A)->refct; 289 mtype = ((PetscObject)A)->type_name; 290 mname = ((PetscObject)A)->name; 291 spptr = A->spptr; 292 293 /* zero these so the destroy below does not free them */ 294 ((PetscObject)A)->type_name = 0; 295 ((PetscObject)A)->name = 0; 296 297 /* free all the interior data structures from mat */ 298 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 299 300 ierr = PetscFree(C->spptr);CHKERRQ(ierr); 301 302 ierr = PetscLayoutDestroy(&A->rmap);CHKERRQ(ierr); 303 ierr = PetscLayoutDestroy(&A->cmap);CHKERRQ(ierr); 304 ierr = PetscFunctionListDestroy(&((PetscObject)A)->qlist);CHKERRQ(ierr); 305 ierr = PetscObjectListDestroy(&((PetscObject)A)->olist);CHKERRQ(ierr); 306 307 /* copy C over to A */ 308 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 309 310 /* return the parts of A we saved */ 311 ((PetscObject)A)->bops = Abops; 312 A->ops = Aops; 313 ((PetscObject)A)->refct = refct; 314 ((PetscObject)A)->type_name = mtype; 315 ((PetscObject)A)->name = mname; 316 A->spptr = spptr; 317 318 /* since these two are copied into A we do not want them destroyed in C */ 319 ((PetscObject)C)->qlist = 0; 320 ((PetscObject)C)->olist = 0; 321 322 ierr = PetscHeaderDestroy(&C);CHKERRQ(ierr); 323 PetscFunctionReturn(0); 324 } 325 /* 326 Replace A's header with that of C; the C object is then destroyed 327 328 This is essentially code moved from MatDestroy() 329 330 This is somewhat different from MatHeaderMerge() it would be nice to merge the code 331 332 Used in DM hence is declared PETSC_EXTERN 333 */ 334 #undef __FUNCT__ 335 #define __FUNCT__ "MatHeaderReplace" 336 PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A,Mat C) 337 { 338 PetscErrorCode ierr; 339 PetscInt refct; 340 341 PetscFunctionBegin; 342 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 343 PetscValidHeaderSpecific(C,MAT_CLASSID,2); 344 if (A == C) PetscFunctionReturn(0); 345 PetscCheckSameComm(A,1,C,2); 346 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); 347 348 /* free all the interior data structures from mat */ 349 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 350 ierr = PetscHeaderDestroy_Private((PetscObject)A);CHKERRQ(ierr); 351 ierr = PetscFree(A->ops);CHKERRQ(ierr); 352 ierr = PetscLayoutDestroy(&A->rmap);CHKERRQ(ierr); 353 ierr = PetscLayoutDestroy(&A->cmap);CHKERRQ(ierr); 354 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 355 356 /* copy C over to A */ 357 refct = ((PetscObject)A)->refct; 358 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 359 360 ((PetscObject)A)->refct = refct; 361 362 ierr = PetscFree(C);CHKERRQ(ierr); 363 PetscFunctionReturn(0); 364 } 365