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 = PetscViewerDestroy(&B->viewonassembly);CHKERRQ(ierr); 187 ierr = PetscOptionsViewer("-mat_view","Display mat with the viewer on MatAssemblyEnd()","MatView",&B->viewonassembly,&B->viewformatonassembly,NULL);CHKERRQ(ierr); 188 ierr = PetscOptionsName("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",&B->checksymmetryonassembly);CHKERRQ(ierr); 189 ierr = PetscOptionsReal("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",0.0,&B->checksymmetrytol,NULL);CHKERRQ(ierr); 190 ierr = PetscOptionsBool("-mat_null_space_test","Checks if provided null space is correct in MatAssemblyEnd()","MatSetNullSpaceTest",PETSC_FALSE,&B->checknullspaceonassembly,NULL);CHKERRQ(ierr); 191 192 if (B->ops->setfromoptions) { 193 ierr = (*B->ops->setfromoptions)(B);CHKERRQ(ierr); 194 } 195 196 flg = PETSC_FALSE; 197 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); 198 if (set) {ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,flg);CHKERRQ(ierr);} 199 flg = PETSC_FALSE; 200 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); 201 if (set) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,flg);CHKERRQ(ierr);} 202 203 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 204 ierr = PetscObjectProcessOptionsHandlers((PetscObject)B);CHKERRQ(ierr); 205 ierr = PetscOptionsEnd();CHKERRQ(ierr); 206 PetscFunctionReturn(0); 207 } 208 209 #undef __FUNCT__ 210 #define __FUNCT__ "MatXAIJSetPreallocation" 211 /*@ 212 MatXAIJSetPreallocation - set preallocation for serial and parallel AIJ, BAIJ, and SBAIJ matrices 213 214 Collective on Mat 215 216 Input Arguments: 217 + A - matrix being preallocated 218 . bs - block size 219 . dnnz - number of nonzero blocks per block row of diagonal part of parallel matrix 220 . onnz - number of nonzero blocks per block row of off-diagonal part of parallel matrix 221 . dnnzu - number of nonzero blocks per block row of upper-triangular part of diagonal part of parallel matrix 222 - onnzu - number of nonzero blocks per block row of upper-triangular part of off-diagonal part of parallel matrix 223 224 Level: beginner 225 226 .seealso: MatSeqAIJSetPreallocation(), MatMPIAIJSetPreallocation(), MatSeqBAIJSetPreallocation(), MatMPIBAIJSetPreallocation(), MatSeqSBAIJSetPreallocation(), MatMPISBAIJSetPreallocation(), 227 PetscSplitOwnership() 228 @*/ 229 PetscErrorCode MatXAIJSetPreallocation(Mat A,PetscInt bs,const PetscInt dnnz[],const PetscInt onnz[],const PetscInt dnnzu[],const PetscInt onnzu[]) 230 { 231 PetscErrorCode ierr; 232 void (*aij)(void); 233 234 PetscFunctionBegin; 235 ierr = MatSetBlockSize(A,bs);CHKERRQ(ierr); 236 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 237 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 238 ierr = MatSeqBAIJSetPreallocation(A,bs,0,dnnz);CHKERRQ(ierr); 239 ierr = MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);CHKERRQ(ierr); 240 ierr = MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);CHKERRQ(ierr); 241 ierr = MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);CHKERRQ(ierr); 242 /* 243 In general, we have to do extra work to preallocate for scalar (AIJ) matrices so we check whether it will do any 244 good before going on with it. 245 */ 246 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 247 if (!aij) { 248 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 249 } 250 if (aij) { 251 if (bs == 1) { 252 ierr = MatSeqAIJSetPreallocation(A,0,dnnz);CHKERRQ(ierr); 253 ierr = MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);CHKERRQ(ierr); 254 } else { /* Convert block-row precallocation to scalar-row */ 255 PetscInt i,m,*sdnnz,*sonnz; 256 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 257 ierr = PetscMalloc2((!!dnnz)*m,&sdnnz,(!!onnz)*m,&sonnz);CHKERRQ(ierr); 258 for (i=0; i<m; i++) { 259 if (dnnz) sdnnz[i] = dnnz[i/bs] * bs; 260 if (onnz) sonnz[i] = onnz[i/bs] * bs; 261 } 262 ierr = MatSeqAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL);CHKERRQ(ierr); 263 ierr = MatMPIAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);CHKERRQ(ierr); 264 ierr = PetscFree2(sdnnz,sonnz);CHKERRQ(ierr); 265 } 266 } 267 PetscFunctionReturn(0); 268 } 269 270 /* 271 Merges some information from Cs header to A; the C object is then destroyed 272 273 This is somewhat different from MatHeaderReplace() it would be nice to merge the code 274 */ 275 #undef __FUNCT__ 276 #define __FUNCT__ "MatHeaderMerge" 277 PetscErrorCode MatHeaderMerge(Mat A,Mat C) 278 { 279 PetscErrorCode ierr; 280 PetscInt refct; 281 PetscOps *Abops; 282 MatOps Aops; 283 char *mtype,*mname; 284 void *spptr; 285 286 PetscFunctionBegin; 287 /* save the parts of A we need */ 288 Abops = ((PetscObject)A)->bops; 289 Aops = A->ops; 290 refct = ((PetscObject)A)->refct; 291 mtype = ((PetscObject)A)->type_name; 292 mname = ((PetscObject)A)->name; 293 spptr = A->spptr; 294 295 /* zero these so the destroy below does not free them */ 296 ((PetscObject)A)->type_name = 0; 297 ((PetscObject)A)->name = 0; 298 299 /* free all the interior data structures from mat */ 300 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 301 302 ierr = PetscFree(C->spptr);CHKERRQ(ierr); 303 304 ierr = PetscLayoutDestroy(&A->rmap);CHKERRQ(ierr); 305 ierr = PetscLayoutDestroy(&A->cmap);CHKERRQ(ierr); 306 ierr = PetscFunctionListDestroy(&((PetscObject)A)->qlist);CHKERRQ(ierr); 307 ierr = PetscObjectListDestroy(&((PetscObject)A)->olist);CHKERRQ(ierr); 308 309 /* copy C over to A */ 310 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 311 312 /* return the parts of A we saved */ 313 ((PetscObject)A)->bops = Abops; 314 A->ops = Aops; 315 ((PetscObject)A)->refct = refct; 316 ((PetscObject)A)->type_name = mtype; 317 ((PetscObject)A)->name = mname; 318 A->spptr = spptr; 319 320 /* since these two are copied into A we do not want them destroyed in C */ 321 ((PetscObject)C)->qlist = 0; 322 ((PetscObject)C)->olist = 0; 323 324 ierr = PetscHeaderDestroy(&C);CHKERRQ(ierr); 325 PetscFunctionReturn(0); 326 } 327 /* 328 Replace A's header with that of C; the C object is then destroyed 329 330 This is essentially code moved from MatDestroy() 331 332 This is somewhat different from MatHeaderMerge() it would be nice to merge the code 333 334 Used in DM hence is declared PETSC_EXTERN 335 */ 336 #undef __FUNCT__ 337 #define __FUNCT__ "MatHeaderReplace" 338 PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A,Mat C) 339 { 340 PetscErrorCode ierr; 341 PetscInt refct; 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 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 361 362 ((PetscObject)A)->refct = refct; 363 364 ierr = PetscFree(C);CHKERRQ(ierr); 365 PetscFunctionReturn(0); 366 } 367