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