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","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 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 = PetscObjectOptionsBegin((PetscObject)B);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 - If the user has not set preallocation for this matrix then a default preallocation that is likely to be inefficient is used. 212 213 Collective on Mat 214 215 Input Parameter: 216 . A - the matrix 217 218 Level: advanced 219 220 Notes: See the Performance chapter of the PETSc users manual for how to preallocate matrices 221 222 .keywords: matrix, create 223 224 .seealso: MatCreateSeqAIJ((), MatCreateMPIAIJ(), 225 MatCreateSeqDense(), MatCreateMPIDense(), 226 MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), 227 MatCreateSeqSBAIJ(), MatCreateMPISBAIJ(), 228 MatConvert() 229 @*/ 230 PetscErrorCode MatSetUpPreallocation(Mat B) 231 { 232 PetscErrorCode ierr; 233 234 PetscFunctionBegin; 235 if (!B->preallocated && B->ops->setuppreallocation) { 236 ierr = PetscInfo(B,"Warning not preallocating matrix storage\n");CHKERRQ(ierr); 237 ierr = (*B->ops->setuppreallocation)(B);CHKERRQ(ierr); 238 } 239 B->preallocated = PETSC_TRUE; 240 PetscFunctionReturn(0); 241 } 242 243 /* 244 Merges some information from Cs header to A; the C object is then destroyed 245 246 This is somewhat different from MatHeaderReplace() it would be nice to merge the code 247 */ 248 #undef __FUNCT__ 249 #define __FUNCT__ "MatHeaderMerge" 250 PetscErrorCode MatHeaderMerge(Mat A,Mat C) 251 { 252 PetscErrorCode ierr; 253 PetscInt refct; 254 PetscOps *Abops; 255 MatOps Aops; 256 char *mtype,*mname; 257 void *spptr; 258 259 PetscFunctionBegin; 260 /* save the parts of A we need */ 261 Abops = ((PetscObject)A)->bops; 262 Aops = A->ops; 263 refct = ((PetscObject)A)->refct; 264 mtype = ((PetscObject)A)->type_name; 265 mname = ((PetscObject)A)->name; 266 spptr = A->spptr; 267 268 /* zero these so the destroy below does not free them */ 269 ((PetscObject)A)->type_name = 0; 270 ((PetscObject)A)->name = 0; 271 272 /* free all the interior data structures from mat */ 273 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 274 275 ierr = PetscFree(C->spptr);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 293 /* since these two are copied into A we do not want them destroyed in C */ 294 ((PetscObject)C)->qlist = 0; 295 ((PetscObject)C)->olist = 0; 296 ierr = PetscHeaderDestroy(&C);CHKERRQ(ierr); 297 PetscFunctionReturn(0); 298 } 299 /* 300 Replace A's header with that of C; the C object is then destroyed 301 302 This is essentially code moved from MatDestroy() 303 304 This is somewhat different from MatHeaderMerge() it would be nice to merge the code 305 */ 306 #undef __FUNCT__ 307 #define __FUNCT__ "MatHeaderReplace" 308 PetscErrorCode MatHeaderReplace(Mat A,Mat C) 309 { 310 PetscErrorCode ierr; 311 PetscInt refct; 312 313 PetscFunctionBegin; 314 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 315 PetscValidHeaderSpecific(C,MAT_CLASSID,2); 316 if (A == C) PetscFunctionReturn(0); 317 PetscCheckSameComm(A,1,C,2); 318 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); 319 320 /* free all the interior data structures from mat */ 321 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 322 ierr = PetscHeaderDestroy_Private((PetscObject)A);CHKERRQ(ierr); 323 ierr = PetscFree(A->ops);CHKERRQ(ierr); 324 ierr = PetscLayoutDestroy(&A->rmap);CHKERRQ(ierr); 325 ierr = PetscLayoutDestroy(&A->cmap);CHKERRQ(ierr); 326 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 327 328 /* copy C over to A */ 329 refct = ((PetscObject)A)->refct; 330 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 331 ((PetscObject)A)->refct = refct; 332 ierr = PetscLogObjectDestroy((PetscObject)C);CHKERRQ(ierr); 333 ierr = PetscFree(C);CHKERRQ(ierr); 334 PetscFunctionReturn(0); 335 } 336