1 #define PETSCMAT_DLL 2 3 #include "private/matimpl.h" /*I "petscmat.h" I*/ 4 5 #if 0 6 #undef __FUNCT__ 7 #define __FUNCT__ "MatPublish_Base" 8 static PetscErrorCode MatPublish_Base(PetscObject obj) 9 { 10 PetscFunctionBegin; 11 PetscFunctionReturn(0); 12 } 13 #endif 14 15 #undef __FUNCT__ 16 #define __FUNCT__ "MatCreate" 17 /*@ 18 MatCreate - Creates a matrix where the type is determined 19 from either a call to MatSetType() or from the options database 20 with a call to MatSetFromOptions(). The default matrix type is 21 AIJ, using the routines MatCreateSeqAIJ() or MatCreateMPIAIJ() 22 if you do not set a type in the options database. If you never 23 call MatSetType() or MatSetFromOptions() it will generate an 24 error when you try to use the matrix. 25 26 Collective on MPI_Comm 27 28 Input Parameter: 29 . comm - MPI communicator 30 31 Output Parameter: 32 . A - the matrix 33 34 Options Database Keys: 35 + -mat_type seqaij - AIJ type, uses MatCreateSeqAIJ() 36 . -mat_type mpiaij - AIJ type, uses MatCreateMPIAIJ() 37 . -mat_type seqdense - dense type, uses MatCreateSeqDense() 38 . -mat_type mpidense - dense type, uses MatCreateMPIDense() 39 . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ() 40 - -mat_type mpibaij - block AIJ type, uses MatCreateMPIBAIJ() 41 42 Even More Options Database Keys: 43 See the manpages for particular formats (e.g., MatCreateSeqAIJ()) 44 for additional format-specific options. 45 46 Notes: 47 48 Level: beginner 49 50 User manual sections: 51 + sec_matcreate 52 - chapter_matrices 53 54 .keywords: matrix, create 55 56 .seealso: MatCreateSeqAIJ(), MatCreateMPIAIJ(), 57 MatCreateSeqDense(), MatCreateMPIDense(), 58 MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), 59 MatCreateSeqSBAIJ(), MatCreateMPISBAIJ(), 60 MatConvert() 61 @*/ 62 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate(MPI_Comm comm,Mat *A) 63 { 64 Mat B; 65 PetscErrorCode ierr; 66 67 PetscFunctionBegin; 68 PetscValidPointer(A,2); 69 70 *A = PETSC_NULL; 71 #ifndef PETSC_USE_DYNAMIC_LIBRARIES 72 ierr = MatInitializePackage(PETSC_NULL);CHKERRQ(ierr); 73 #endif 74 75 ierr = PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_CLASSID,0,"Mat",comm,MatDestroy,MatView);CHKERRQ(ierr); 76 ierr = PetscLayoutCreate(comm,&B->rmap);CHKERRQ(ierr); 77 ierr = PetscLayoutCreate(comm,&B->cmap);CHKERRQ(ierr); 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 Level: beginner 110 111 .seealso: MatGetSize(), PetscSplitOwnership() 112 @*/ 113 PetscErrorCode PETSCMAT_DLLEXPORT MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N) 114 { 115 PetscErrorCode ierr; 116 117 PetscFunctionBegin; 118 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 119 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); 120 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); 121 if (A->ops->setsizes) { 122 /* Since this will not be set until the type has been set, this will NOT be called on the initial 123 call of MatSetSizes() (which must be called BEFORE MatSetType() */ 124 ierr = (*A->ops->setsizes)(A,m,n,M,N);CHKERRQ(ierr); 125 } else { 126 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); 127 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); 128 } 129 A->rmap->n = m; 130 A->cmap->n = n; 131 A->rmap->N = M; 132 A->cmap->N = N; 133 134 PetscFunctionReturn(0); 135 } 136 137 #undef __FUNCT__ 138 #define __FUNCT__ "MatSetFromOptions" 139 /*@ 140 MatSetFromOptions - Creates a matrix where the type is determined 141 from the options database. Generates a parallel MPI matrix if the 142 communicator has more than one processor. The default matrix type is 143 AIJ, using the routines MatCreateSeqAIJ() and MatCreateMPIAIJ() if 144 you do not select a type in the options database. 145 146 Collective on Mat 147 148 Input Parameter: 149 . A - the matrix 150 151 Options Database Keys: 152 + -mat_type seqaij - AIJ type, uses MatCreateSeqAIJ() 153 . -mat_type mpiaij - AIJ type, uses MatCreateMPIAIJ() 154 . -mat_type seqdense - dense type, uses MatCreateSeqDense() 155 . -mat_type mpidense - dense type, uses MatCreateMPIDense() 156 . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ() 157 - -mat_type mpibaij - block AIJ type, uses MatCreateMPIBAIJ() 158 159 Even More Options Database Keys: 160 See the manpages for particular formats (e.g., MatCreateSeqAIJ()) 161 for additional format-specific options. 162 163 Level: beginner 164 165 .keywords: matrix, create 166 167 .seealso: MatCreateSeqAIJ((), MatCreateMPIAIJ(), 168 MatCreateSeqDense(), MatCreateMPIDense(), 169 MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), 170 MatCreateSeqSBAIJ(), MatCreateMPISBAIJ(), 171 MatConvert() 172 @*/ 173 PetscErrorCode PETSCMAT_DLLEXPORT MatSetFromOptions(Mat B) 174 { 175 PetscErrorCode ierr; 176 const char *deft = MATAIJ; 177 char type[256]; 178 PetscTruth flg; 179 180 PetscFunctionBegin; 181 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 182 183 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Matrix options","Mat");CHKERRQ(ierr); 184 ierr = PetscOptionsList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);CHKERRQ(ierr); 185 if (flg) { 186 ierr = MatSetType(B,type);CHKERRQ(ierr); 187 } else if (!((PetscObject)B)->type_name) { 188 ierr = MatSetType(B,deft);CHKERRQ(ierr); 189 } 190 191 if (B->ops->setfromoptions) { 192 ierr = (*B->ops->setfromoptions)(B);CHKERRQ(ierr); 193 } 194 195 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 196 ierr = PetscObjectProcessOptionsHandlers((PetscObject)B);CHKERRQ(ierr); 197 ierr = PetscOptionsEnd();CHKERRQ(ierr); 198 199 PetscFunctionReturn(0); 200 } 201 202 #undef __FUNCT__ 203 #define __FUNCT__ "MatSetUpPreallocation" 204 /*@ 205 MatSetUpPreallocation 206 207 Collective on Mat 208 209 Input Parameter: 210 . A - the matrix 211 212 Level: beginner 213 214 .keywords: matrix, create 215 216 .seealso: MatCreateSeqAIJ((), MatCreateMPIAIJ(), 217 MatCreateSeqDense(), MatCreateMPIDense(), 218 MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), 219 MatCreateSeqSBAIJ(), MatCreateMPISBAIJ(), 220 MatConvert() 221 @*/ 222 PetscErrorCode PETSCMAT_DLLEXPORT MatSetUpPreallocation(Mat B) 223 { 224 PetscErrorCode ierr; 225 226 PetscFunctionBegin; 227 if (!B->preallocated && B->ops->setuppreallocation) { 228 ierr = PetscInfo(B,"Warning not preallocating matrix storage\n");CHKERRQ(ierr); 229 ierr = (*B->ops->setuppreallocation)(B);CHKERRQ(ierr); 230 } 231 B->preallocated = PETSC_TRUE; 232 PetscFunctionReturn(0); 233 } 234 235 /* 236 Merges some information from Cs header to A; the C object is then destroyed 237 238 This is somewhat different from MatHeaderReplace() it would be nice to merge the code 239 */ 240 #undef __FUNCT__ 241 #define __FUNCT__ "MatHeaderMerge" 242 PetscErrorCode MatHeaderMerge(Mat A,Mat C) 243 { 244 PetscErrorCode ierr; 245 PetscInt refct; 246 PetscOps *Abops; 247 MatOps Aops; 248 char *mtype,*mname; 249 void *spptr; 250 ISLocalToGlobalMapping mapping,bmapping; 251 252 PetscFunctionBegin; 253 /* save the parts of A we need */ 254 Abops = ((PetscObject)A)->bops; 255 Aops = A->ops; 256 refct = ((PetscObject)A)->refct; 257 mtype = ((PetscObject)A)->type_name; 258 mname = ((PetscObject)A)->name; 259 spptr = A->spptr; 260 mapping = A->mapping; 261 bmapping = A->bmapping; 262 263 /* zero these so the destroy below does not free them */ 264 ((PetscObject)A)->type_name = 0; 265 ((PetscObject)A)->name = 0; 266 267 /* free all the interior data structures from mat */ 268 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 269 270 ierr = PetscFree(C->spptr);CHKERRQ(ierr); 271 if (C->mapping) { 272 ierr = ISLocalToGlobalMappingDestroy(C->mapping);CHKERRQ(ierr); 273 } 274 if (C->bmapping) { 275 ierr = ISLocalToGlobalMappingDestroy(C->bmapping);CHKERRQ(ierr); 276 } 277 278 ierr = PetscLayoutDestroy(A->rmap);CHKERRQ(ierr); 279 ierr = PetscLayoutDestroy(A->cmap);CHKERRQ(ierr); 280 ierr = PetscFListDestroy(&((PetscObject)A)->qlist);CHKERRQ(ierr); 281 ierr = PetscOListDestroy(((PetscObject)A)->olist);CHKERRQ(ierr); 282 283 /* copy C over to A */ 284 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 285 286 /* return the parts of A we saved */ 287 ((PetscObject)A)->bops = Abops; 288 A->ops = Aops; 289 ((PetscObject)A)->refct = refct; 290 ((PetscObject)A)->type_name = mtype; 291 ((PetscObject)A)->name = mname; 292 A->spptr = spptr; 293 A->mapping = mapping; 294 A->bmapping = bmapping; 295 296 /* since these two are copied into A we do not want them destroyed in C */ 297 ((PetscObject)C)->qlist = 0; 298 ((PetscObject)C)->olist = 0; 299 ierr = PetscHeaderDestroy(C);CHKERRQ(ierr); 300 PetscFunctionReturn(0); 301 } 302 /* 303 Replace A's header with that of C; the C object is then destroyed 304 305 This is essentially code moved from MatDestroy() 306 307 This is somewhat different from MatHeaderMerge() it would be nice to merge the code 308 */ 309 #undef __FUNCT__ 310 #define __FUNCT__ "MatHeaderReplace" 311 PetscErrorCode MatHeaderReplace(Mat A,Mat C) 312 { 313 PetscErrorCode ierr; 314 315 PetscFunctionBegin; 316 if (A == C) PetscFunctionReturn(0); 317 318 /* free all the interior data structures from mat */ 319 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 320 ierr = PetscHeaderDestroy_Private((PetscObject)A);CHKERRQ(ierr); 321 ierr = PetscFree(A->ops);CHKERRQ(ierr); 322 ierr = PetscLayoutDestroy(A->rmap);CHKERRQ(ierr); 323 ierr = PetscLayoutDestroy(A->cmap);CHKERRQ(ierr); 324 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 325 if (A->mapping) { 326 ierr = ISLocalToGlobalMappingDestroy(A->mapping);CHKERRQ(ierr); 327 } 328 if (A->bmapping) { 329 ierr = ISLocalToGlobalMappingDestroy(A->bmapping);CHKERRQ(ierr); 330 } 331 332 /* copy C over to A */ 333 if (C) { 334 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 335 ierr = PetscLogObjectDestroy((PetscObject)C);CHKERRQ(ierr); 336 ierr = PetscFree(C);CHKERRQ(ierr); 337 } 338 PetscFunctionReturn(0); 339 } 340