1 2 #include <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 MatCreateMPIAIJ() 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 MatCreateMPIAIJ() 36 . -mat_type seqdense - dense type, uses MatCreateSeqDense() 37 . -mat_type mpidense - dense type, uses MatCreateMPIDense() 38 . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ() 39 - -mat_type mpibaij - block AIJ type, uses MatCreateMPIBAIJ() 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(), MatCreateMPIAIJ(), 56 MatCreateSeqDense(), MatCreateMPIDense(), 57 MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), 58 MatCreateSeqSBAIJ(), MatCreateMPISBAIJ(), 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 #ifndef 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",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 Level: beginner 109 110 .seealso: MatGetSize(), PetscSplitOwnership() 111 @*/ 112 PetscErrorCode MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N) 113 { 114 PetscErrorCode ierr; 115 116 PetscFunctionBegin; 117 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 118 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); 119 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); 120 if (A->ops->setsizes) { 121 /* Since this will not be set until the type has been set, this will NOT be called on the initial 122 call of MatSetSizes() (which must be called BEFORE MatSetType() */ 123 ierr = (*A->ops->setsizes)(A,m,n,M,N);CHKERRQ(ierr); 124 } else { 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 } 128 A->rmap->n = m; 129 A->cmap->n = n; 130 A->rmap->N = M; 131 A->cmap->N = N; 132 133 PetscFunctionReturn(0); 134 } 135 136 #undef __FUNCT__ 137 #define __FUNCT__ "MatSetFromOptions" 138 /*@ 139 MatSetFromOptions - Creates a matrix where the type is determined 140 from the options database. Generates a parallel MPI matrix if the 141 communicator has more than one processor. The default matrix type is 142 AIJ, using the routines MatCreateSeqAIJ() and MatCreateMPIAIJ() if 143 you do not select a type in the options database. 144 145 Collective on Mat 146 147 Input Parameter: 148 . A - the matrix 149 150 Options Database Keys: 151 + -mat_type seqaij - AIJ type, uses MatCreateSeqAIJ() 152 . -mat_type mpiaij - AIJ type, uses MatCreateMPIAIJ() 153 . -mat_type seqdense - dense type, uses MatCreateSeqDense() 154 . -mat_type mpidense - dense type, uses MatCreateMPIDense() 155 . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ() 156 - -mat_type mpibaij - block AIJ type, uses MatCreateMPIBAIJ() 157 158 Even More Options Database Keys: 159 See the manpages for particular formats (e.g., MatCreateSeqAIJ()) 160 for additional format-specific options. 161 162 Level: beginner 163 164 .keywords: matrix, create 165 166 .seealso: MatCreateSeqAIJ((), MatCreateMPIAIJ(), 167 MatCreateSeqDense(), MatCreateMPIDense(), 168 MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), 169 MatCreateSeqSBAIJ(), MatCreateMPISBAIJ(), 170 MatConvert() 171 @*/ 172 PetscErrorCode MatSetFromOptions(Mat B) 173 { 174 PetscErrorCode ierr; 175 const char *deft = MATAIJ; 176 char type[256]; 177 PetscBool flg,set; 178 179 PetscFunctionBegin; 180 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 181 182 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Matrix options","Mat");CHKERRQ(ierr); 183 ierr = PetscOptionsList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);CHKERRQ(ierr); 184 if (flg) { 185 ierr = MatSetType(B,type);CHKERRQ(ierr); 186 } else if (!((PetscObject)B)->type_name) { 187 ierr = MatSetType(B,deft);CHKERRQ(ierr); 188 } 189 190 if (B->ops->setfromoptions) { 191 ierr = (*B->ops->setfromoptions)(B);CHKERRQ(ierr); 192 } 193 194 flg = PETSC_FALSE; 195 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); 196 if (set) {ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,flg);CHKERRQ(ierr);} 197 flg = PETSC_FALSE; 198 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); 199 if (set) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,flg);CHKERRQ(ierr);} 200 201 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 202 ierr = PetscObjectProcessOptionsHandlers((PetscObject)B);CHKERRQ(ierr); 203 ierr = PetscOptionsEnd();CHKERRQ(ierr); 204 205 PetscFunctionReturn(0); 206 } 207 208 #undef __FUNCT__ 209 #define __FUNCT__ "MatSetUpPreallocation" 210 /*@ 211 MatSetUpPreallocation 212 213 Collective on Mat 214 215 Input Parameter: 216 . A - the matrix 217 218 Level: beginner 219 220 .keywords: matrix, create 221 222 .seealso: MatCreateSeqAIJ((), MatCreateMPIAIJ(), 223 MatCreateSeqDense(), MatCreateMPIDense(), 224 MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), 225 MatCreateSeqSBAIJ(), MatCreateMPISBAIJ(), 226 MatConvert() 227 @*/ 228 PetscErrorCode MatSetUpPreallocation(Mat B) 229 { 230 PetscErrorCode ierr; 231 232 PetscFunctionBegin; 233 if (!B->preallocated && B->ops->setuppreallocation) { 234 ierr = PetscInfo(B,"Warning not preallocating matrix storage\n");CHKERRQ(ierr); 235 ierr = (*B->ops->setuppreallocation)(B);CHKERRQ(ierr); 236 } 237 B->preallocated = PETSC_TRUE; 238 PetscFunctionReturn(0); 239 } 240 241 /* 242 Merges some information from Cs header to A; the C object is then destroyed 243 244 This is somewhat different from MatHeaderReplace() it would be nice to merge the code 245 */ 246 #undef __FUNCT__ 247 #define __FUNCT__ "MatHeaderMerge" 248 PetscErrorCode MatHeaderMerge(Mat A,Mat C) 249 { 250 PetscErrorCode ierr; 251 PetscInt refct; 252 PetscOps *Abops; 253 MatOps Aops; 254 char *mtype,*mname; 255 void *spptr; 256 ISLocalToGlobalMapping rmapping,cmapping,rbmapping,cbmapping; 257 258 PetscFunctionBegin; 259 /* save the parts of A we need */ 260 Abops = ((PetscObject)A)->bops; 261 Aops = A->ops; 262 refct = ((PetscObject)A)->refct; 263 mtype = ((PetscObject)A)->type_name; 264 mname = ((PetscObject)A)->name; 265 spptr = A->spptr; 266 rmapping = A->rmapping; 267 cmapping = A->cmapping; 268 rbmapping = A->rbmapping; 269 cbmapping = A->cbmapping; 270 271 /* zero these so the destroy below does not free them */ 272 ((PetscObject)A)->type_name = 0; 273 ((PetscObject)A)->name = 0; 274 275 /* free all the interior data structures from mat */ 276 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 277 278 ierr = PetscFree(C->spptr);CHKERRQ(ierr); 279 ierr = ISLocalToGlobalMappingDestroy(&C->rmapping);CHKERRQ(ierr); 280 ierr = ISLocalToGlobalMappingDestroy(&C->cmapping);CHKERRQ(ierr); 281 ierr = ISLocalToGlobalMappingDestroy(&C->rbmapping);CHKERRQ(ierr); 282 ierr = ISLocalToGlobalMappingDestroy(&C->cbmapping);CHKERRQ(ierr); 283 284 ierr = PetscLayoutDestroy(&A->rmap);CHKERRQ(ierr); 285 ierr = PetscLayoutDestroy(&A->cmap);CHKERRQ(ierr); 286 ierr = PetscFListDestroy(&((PetscObject)A)->qlist);CHKERRQ(ierr); 287 ierr = PetscOListDestroy(&((PetscObject)A)->olist);CHKERRQ(ierr); 288 289 /* copy C over to A */ 290 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 291 292 /* return the parts of A we saved */ 293 ((PetscObject)A)->bops = Abops; 294 A->ops = Aops; 295 ((PetscObject)A)->refct = refct; 296 ((PetscObject)A)->type_name = mtype; 297 ((PetscObject)A)->name = mname; 298 A->spptr = spptr; 299 A->rmapping = rmapping; 300 A->cmapping = cmapping; 301 A->rbmapping = rbmapping; 302 A->cbmapping = cbmapping; 303 304 /* since these two are copied into A we do not want them destroyed in C */ 305 ((PetscObject)C)->qlist = 0; 306 ((PetscObject)C)->olist = 0; 307 ierr = PetscHeaderDestroy(&C);CHKERRQ(ierr); 308 PetscFunctionReturn(0); 309 } 310 /* 311 Replace A's header with that of C; the C object is then destroyed 312 313 This is essentially code moved from MatDestroy() 314 315 This is somewhat different from MatHeaderMerge() it would be nice to merge the code 316 */ 317 #undef __FUNCT__ 318 #define __FUNCT__ "MatHeaderReplace" 319 PetscErrorCode MatHeaderReplace(Mat A,Mat C) 320 { 321 PetscErrorCode ierr; 322 PetscInt refct; 323 324 PetscFunctionBegin; 325 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 326 PetscValidHeaderSpecific(C,MAT_CLASSID,2); 327 if (A == C) PetscFunctionReturn(0); 328 PetscCheckSameComm(A,1,C,2); 329 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); 330 331 /* free all the interior data structures from mat */ 332 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 333 ierr = PetscHeaderDestroy_Private((PetscObject)A);CHKERRQ(ierr); 334 ierr = PetscFree(A->ops);CHKERRQ(ierr); 335 ierr = PetscLayoutDestroy(&A->rmap);CHKERRQ(ierr); 336 ierr = PetscLayoutDestroy(&A->cmap);CHKERRQ(ierr); 337 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 338 ierr = ISLocalToGlobalMappingDestroy(&A->rmapping);CHKERRQ(ierr); 339 ierr = ISLocalToGlobalMappingDestroy(&A->cmapping);CHKERRQ(ierr); 340 ierr = ISLocalToGlobalMappingDestroy(&A->rbmapping);CHKERRQ(ierr); 341 ierr = ISLocalToGlobalMappingDestroy(&A->cbmapping);CHKERRQ(ierr); 342 343 /* copy C over to A */ 344 refct = ((PetscObject)A)->refct; 345 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 346 ((PetscObject)A)->refct = refct; 347 ierr = PetscLogObjectDestroy((PetscObject)C);CHKERRQ(ierr); 348 ierr = PetscFree(C);CHKERRQ(ierr); 349 PetscFunctionReturn(0); 350 } 351