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