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