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