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 257 PetscFunctionBegin; 258 /* save the parts of A we need */ 259 Abops = ((PetscObject)A)->bops; 260 Aops = A->ops; 261 refct = ((PetscObject)A)->refct; 262 mtype = ((PetscObject)A)->type_name; 263 mname = ((PetscObject)A)->name; 264 spptr = A->spptr; 265 266 /* zero these so the destroy below does not free them */ 267 ((PetscObject)A)->type_name = 0; 268 ((PetscObject)A)->name = 0; 269 270 /* free all the interior data structures from mat */ 271 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 272 273 ierr = PetscFree(C->spptr);CHKERRQ(ierr); 274 275 ierr = PetscLayoutDestroy(&A->rmap);CHKERRQ(ierr); 276 ierr = PetscLayoutDestroy(&A->cmap);CHKERRQ(ierr); 277 ierr = PetscFListDestroy(&((PetscObject)A)->qlist);CHKERRQ(ierr); 278 ierr = PetscOListDestroy(&((PetscObject)A)->olist);CHKERRQ(ierr); 279 280 /* copy C over to A */ 281 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 282 283 /* return the parts of A we saved */ 284 ((PetscObject)A)->bops = Abops; 285 A->ops = Aops; 286 ((PetscObject)A)->refct = refct; 287 ((PetscObject)A)->type_name = mtype; 288 ((PetscObject)A)->name = mname; 289 A->spptr = spptr; 290 291 /* since these two are copied into A we do not want them destroyed in C */ 292 ((PetscObject)C)->qlist = 0; 293 ((PetscObject)C)->olist = 0; 294 ierr = PetscHeaderDestroy(&C);CHKERRQ(ierr); 295 PetscFunctionReturn(0); 296 } 297 /* 298 Replace A's header with that of C; the C object is then destroyed 299 300 This is essentially code moved from MatDestroy() 301 302 This is somewhat different from MatHeaderMerge() it would be nice to merge the code 303 */ 304 #undef __FUNCT__ 305 #define __FUNCT__ "MatHeaderReplace" 306 PetscErrorCode MatHeaderReplace(Mat A,Mat C) 307 { 308 PetscErrorCode ierr; 309 PetscInt refct; 310 311 PetscFunctionBegin; 312 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 313 PetscValidHeaderSpecific(C,MAT_CLASSID,2); 314 if (A == C) PetscFunctionReturn(0); 315 PetscCheckSameComm(A,1,C,2); 316 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); 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 326 /* copy C over to A */ 327 refct = ((PetscObject)A)->refct; 328 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 329 ((PetscObject)A)->refct = refct; 330 ierr = PetscLogObjectDestroy((PetscObject)C);CHKERRQ(ierr); 331 ierr = PetscFree(C);CHKERRQ(ierr); 332 PetscFunctionReturn(0); 333 } 334