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 MatCreateAIJ() 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 MatCreateAIJ() 36 . -mat_type seqdense - dense type, uses MatCreateSeqDense() 37 . -mat_type mpidense - dense type, uses MatCreateDense() 38 . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ() 39 - -mat_type mpibaij - block AIJ type, uses MatCreateBAIJ() 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(), MatCreateAIJ(), 56 MatCreateSeqDense(), MatCreateDense(), 57 MatCreateSeqBAIJ(), MatCreateBAIJ(), 58 MatCreateSeqSBAIJ(), MatCreateSBAIJ(), 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 MatCreateAIJ() 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 MatCreateAIJ() 153 . -mat_type seqdense - dense type, uses MatCreateSeqDense() 154 . -mat_type mpidense - dense type, uses MatCreateDense() 155 . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ() 156 - -mat_type mpibaij - block AIJ type, uses MatCreateBAIJ() 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((), MatCreateAIJ(), 167 MatCreateSeqDense(), MatCreateDense(), 168 MatCreateSeqBAIJ(), MatCreateBAIJ(), 169 MatCreateSeqSBAIJ(), MatCreateSBAIJ(), 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__ "MatXAIJSetPreallocation" 210 /*@ 211 MatXAIJSetPreallocation - set preallocation for serial and parallel AIJ, BAIJ, and SBAIJ matrices 212 213 Collective on Mat 214 215 Input Arguments: 216 + A - matrix being preallocated 217 . bs - block size 218 . dnnz - number of nonzero blocks per block row of diagonal part of parallel matrix 219 . onnz - number of nonzero blocks per block row of off-diagonal part of parallel matrix 220 . dnnzu - number of nonzero blocks per block row of upper-triangular part of diagonal part of parallel matrix 221 - onnzu - number of nonzero blocks per block row of upper-triangular part of off-diagonal part of parallel matrix 222 223 Level: beginner 224 225 .seealso: MatSeqAIJSetPreallocation(), MatMPIAIJSetPreallocation(), MatSeqBAIJSetPreallocation(), MatMPIBAIJSetPreallocation(), MatSeqSBAIJSetPreallocation(), MatMPISBAIJSetPreallocation() 226 @*/ 227 PetscErrorCode MatXAIJSetPreallocation(Mat A,PetscInt bs,const PetscInt *dnnz,const PetscInt *onnz,const PetscInt *dnnzu,const PetscInt *onnzu) 228 { 229 PetscErrorCode ierr; 230 void (*aij)(void); 231 232 PetscFunctionBegin; 233 ierr = MatSeqBAIJSetPreallocation(A,bs,0,dnnz);CHKERRQ(ierr); 234 ierr = MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);CHKERRQ(ierr); 235 ierr = MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);CHKERRQ(ierr); 236 ierr = MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);CHKERRQ(ierr); 237 /* 238 In general, we have to do extra work to preallocate for scalar (AIJ) matrices so we check whether it will do any 239 good before going on with it. 240 */ 241 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 242 if (!aij) { 243 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 244 } 245 if (aij) { 246 if (bs == 1) { 247 ierr = MatSeqAIJSetPreallocation(A,0,dnnz);CHKERRQ(ierr); 248 ierr = MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);CHKERRQ(ierr); 249 } else { /* Convert block-row precallocation to scalar-row */ 250 PetscInt i,m,*sdnnz,*sonnz; 251 ierr = MatGetLocalSize(A,&m,PETSC_NULL);CHKERRQ(ierr); 252 ierr = PetscMalloc2((!!dnnz)*m*bs,PetscInt,&sdnnz,(!!onnz)*m*bs,PetscInt,&sonnz);CHKERRQ(ierr); 253 for (i=0; i<m*bs; i++) { 254 if (dnnz) sdnnz[i] = dnnz[i/bs] * bs; 255 if (onnz) sonnz[i] = onnz[i/bs] * bs; 256 } 257 ierr = MatSeqAIJSetPreallocation(A,0,dnnz?sdnnz:PETSC_NULL);CHKERRQ(ierr); 258 ierr = MatMPIAIJSetPreallocation(A,0,dnnz?sdnnz:PETSC_NULL,0,onnz?sonnz:PETSC_NULL);CHKERRQ(ierr); 259 ierr = PetscFree2(sdnnz,sonnz);CHKERRQ(ierr); 260 } 261 } 262 ierr = MatSetBlockSize(A,bs);CHKERRQ(ierr); 263 PetscFunctionReturn(0); 264 } 265 266 /* 267 Merges some information from Cs header to A; the C object is then destroyed 268 269 This is somewhat different from MatHeaderReplace() it would be nice to merge the code 270 */ 271 #undef __FUNCT__ 272 #define __FUNCT__ "MatHeaderMerge" 273 PetscErrorCode MatHeaderMerge(Mat A,Mat C) 274 { 275 PetscErrorCode ierr; 276 PetscInt refct; 277 PetscOps *Abops; 278 MatOps Aops; 279 char *mtype,*mname; 280 void *spptr; 281 282 PetscFunctionBegin; 283 /* save the parts of A we need */ 284 Abops = ((PetscObject)A)->bops; 285 Aops = A->ops; 286 refct = ((PetscObject)A)->refct; 287 mtype = ((PetscObject)A)->type_name; 288 mname = ((PetscObject)A)->name; 289 spptr = A->spptr; 290 291 /* zero these so the destroy below does not free them */ 292 ((PetscObject)A)->type_name = 0; 293 ((PetscObject)A)->name = 0; 294 295 /* free all the interior data structures from mat */ 296 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 297 298 ierr = PetscFree(C->spptr);CHKERRQ(ierr); 299 300 ierr = PetscLayoutDestroy(&A->rmap);CHKERRQ(ierr); 301 ierr = PetscLayoutDestroy(&A->cmap);CHKERRQ(ierr); 302 ierr = PetscFListDestroy(&((PetscObject)A)->qlist);CHKERRQ(ierr); 303 ierr = PetscOListDestroy(&((PetscObject)A)->olist);CHKERRQ(ierr); 304 305 /* copy C over to A */ 306 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 307 308 /* return the parts of A we saved */ 309 ((PetscObject)A)->bops = Abops; 310 A->ops = Aops; 311 ((PetscObject)A)->refct = refct; 312 ((PetscObject)A)->type_name = mtype; 313 ((PetscObject)A)->name = mname; 314 A->spptr = spptr; 315 316 /* since these two are copied into A we do not want them destroyed in C */ 317 ((PetscObject)C)->qlist = 0; 318 ((PetscObject)C)->olist = 0; 319 ierr = PetscHeaderDestroy(&C);CHKERRQ(ierr); 320 PetscFunctionReturn(0); 321 } 322 /* 323 Replace A's header with that of C; the C object is then destroyed 324 325 This is essentially code moved from MatDestroy() 326 327 This is somewhat different from MatHeaderMerge() it would be nice to merge the code 328 */ 329 #undef __FUNCT__ 330 #define __FUNCT__ "MatHeaderReplace" 331 PetscErrorCode MatHeaderReplace(Mat A,Mat C) 332 { 333 PetscErrorCode ierr; 334 PetscInt refct; 335 336 PetscFunctionBegin; 337 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 338 PetscValidHeaderSpecific(C,MAT_CLASSID,2); 339 if (A == C) PetscFunctionReturn(0); 340 PetscCheckSameComm(A,1,C,2); 341 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); 342 343 /* free all the interior data structures from mat */ 344 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 345 ierr = PetscHeaderDestroy_Private((PetscObject)A);CHKERRQ(ierr); 346 ierr = PetscFree(A->ops);CHKERRQ(ierr); 347 ierr = PetscLayoutDestroy(&A->rmap);CHKERRQ(ierr); 348 ierr = PetscLayoutDestroy(&A->cmap);CHKERRQ(ierr); 349 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 350 351 /* copy C over to A */ 352 refct = ((PetscObject)A)->refct; 353 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 354 ((PetscObject)A)->refct = refct; 355 ierr = PetscLogObjectDestroy((PetscObject)C);CHKERRQ(ierr); 356 ierr = PetscFree(C);CHKERRQ(ierr); 357 PetscFunctionReturn(0); 358 } 359