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; 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 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 195 ierr = PetscObjectProcessOptionsHandlers((PetscObject)B);CHKERRQ(ierr); 196 ierr = PetscOptionsEnd();CHKERRQ(ierr); 197 198 PetscFunctionReturn(0); 199 } 200 201 #undef __FUNCT__ 202 #define __FUNCT__ "MatSetUpPreallocation" 203 /*@ 204 MatSetUpPreallocation 205 206 Collective on Mat 207 208 Input Parameter: 209 . A - the matrix 210 211 Level: beginner 212 213 .keywords: matrix, create 214 215 .seealso: MatCreateSeqAIJ((), MatCreateMPIAIJ(), 216 MatCreateSeqDense(), MatCreateMPIDense(), 217 MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), 218 MatCreateSeqSBAIJ(), MatCreateMPISBAIJ(), 219 MatConvert() 220 @*/ 221 PetscErrorCode MatSetUpPreallocation(Mat B) 222 { 223 PetscErrorCode ierr; 224 225 PetscFunctionBegin; 226 if (!B->preallocated && B->ops->setuppreallocation) { 227 ierr = PetscInfo(B,"Warning not preallocating matrix storage\n");CHKERRQ(ierr); 228 ierr = (*B->ops->setuppreallocation)(B);CHKERRQ(ierr); 229 } 230 B->preallocated = PETSC_TRUE; 231 PetscFunctionReturn(0); 232 } 233 234 /* 235 Merges some information from Cs header to A; the C object is then destroyed 236 237 This is somewhat different from MatHeaderReplace() it would be nice to merge the code 238 */ 239 #undef __FUNCT__ 240 #define __FUNCT__ "MatHeaderMerge" 241 PetscErrorCode MatHeaderMerge(Mat A,Mat C) 242 { 243 PetscErrorCode ierr; 244 PetscInt refct; 245 PetscOps *Abops; 246 MatOps Aops; 247 char *mtype,*mname; 248 void *spptr; 249 ISLocalToGlobalMapping rmapping,cmapping,rbmapping,cbmapping; 250 251 PetscFunctionBegin; 252 /* save the parts of A we need */ 253 Abops = ((PetscObject)A)->bops; 254 Aops = A->ops; 255 refct = ((PetscObject)A)->refct; 256 mtype = ((PetscObject)A)->type_name; 257 mname = ((PetscObject)A)->name; 258 spptr = A->spptr; 259 rmapping = A->rmapping; 260 cmapping = A->cmapping; 261 rbmapping = A->rbmapping; 262 cbmapping = A->cbmapping; 263 264 /* zero these so the destroy below does not free them */ 265 ((PetscObject)A)->type_name = 0; 266 ((PetscObject)A)->name = 0; 267 268 /* free all the interior data structures from mat */ 269 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 270 271 ierr = PetscFree(C->spptr);CHKERRQ(ierr); 272 if (C->rmapping) {ierr = ISLocalToGlobalMappingDestroy(C->rmapping);CHKERRQ(ierr);} 273 if (C->cmapping) {ierr = ISLocalToGlobalMappingDestroy(C->cmapping);CHKERRQ(ierr);} 274 if (C->rbmapping) {ierr = ISLocalToGlobalMappingDestroy(C->rbmapping);CHKERRQ(ierr);} 275 if (C->cbmapping) {ierr = ISLocalToGlobalMappingDestroy(C->cbmapping);CHKERRQ(ierr);} 276 277 ierr = PetscLayoutDestroy(A->rmap);CHKERRQ(ierr); 278 ierr = PetscLayoutDestroy(A->cmap);CHKERRQ(ierr); 279 ierr = PetscFListDestroy(&((PetscObject)A)->qlist);CHKERRQ(ierr); 280 ierr = PetscOListDestroy(((PetscObject)A)->olist);CHKERRQ(ierr); 281 282 /* copy C over to A */ 283 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 284 285 /* return the parts of A we saved */ 286 ((PetscObject)A)->bops = Abops; 287 A->ops = Aops; 288 ((PetscObject)A)->refct = refct; 289 ((PetscObject)A)->type_name = mtype; 290 ((PetscObject)A)->name = mname; 291 A->spptr = spptr; 292 A->rmapping = rmapping; 293 A->cmapping = cmapping; 294 A->rbmapping = rbmapping; 295 A->cbmapping = cbmapping; 296 297 /* since these two are copied into A we do not want them destroyed in C */ 298 ((PetscObject)C)->qlist = 0; 299 ((PetscObject)C)->olist = 0; 300 ierr = PetscHeaderDestroy(C);CHKERRQ(ierr); 301 PetscFunctionReturn(0); 302 } 303 /* 304 Replace A's header with that of C; the C object is then destroyed 305 306 This is essentially code moved from MatDestroy() 307 308 This is somewhat different from MatHeaderMerge() it would be nice to merge the code 309 */ 310 #undef __FUNCT__ 311 #define __FUNCT__ "MatHeaderReplace" 312 PetscErrorCode MatHeaderReplace(Mat A,Mat C) 313 { 314 PetscErrorCode ierr; 315 PetscInt refct; 316 317 PetscFunctionBegin; 318 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 319 PetscValidHeaderSpecific(C,MAT_CLASSID,2); 320 if (A == C) PetscFunctionReturn(0); 321 PetscCheckSameComm(A,1,C,2); 322 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); 323 324 /* free all the interior data structures from mat */ 325 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 326 ierr = PetscHeaderDestroy_Private((PetscObject)A);CHKERRQ(ierr); 327 ierr = PetscFree(A->ops);CHKERRQ(ierr); 328 ierr = PetscLayoutDestroy(A->rmap);CHKERRQ(ierr); 329 ierr = PetscLayoutDestroy(A->cmap);CHKERRQ(ierr); 330 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 331 if (A->rmapping) {ierr = ISLocalToGlobalMappingDestroy(A->rmapping);CHKERRQ(ierr);} 332 if (A->cmapping) {ierr = ISLocalToGlobalMappingDestroy(A->cmapping);CHKERRQ(ierr);} 333 if (A->rbmapping) {ierr = ISLocalToGlobalMappingDestroy(A->rbmapping);CHKERRQ(ierr);} 334 if (A->cbmapping) {ierr = ISLocalToGlobalMappingDestroy(A->cbmapping);CHKERRQ(ierr);} 335 336 /* copy C over to A */ 337 refct = ((PetscObject)A)->refct; 338 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 339 ((PetscObject)A)->refct = refct; 340 ierr = PetscLogObjectDestroy((PetscObject)C);CHKERRQ(ierr); 341 ierr = PetscFree(C);CHKERRQ(ierr); 342 PetscFunctionReturn(0); 343 } 344