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