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