1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 PETSC_INTERN PetscErrorCode MatSetBlockSizes_Default(Mat mat,PetscInt rbs, PetscInt cbs) 5 { 6 PetscFunctionBegin; 7 if (!mat->preallocated) PetscFunctionReturn(0); 8 if (mat->rmap->bs > 0 && mat->rmap->bs != rbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot change row block size %D to %D\n",mat->rmap->bs,rbs); 9 if (mat->cmap->bs > 0 && mat->cmap->bs != cbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot change column block size %D to %D\n",mat->cmap->bs,cbs); 10 PetscFunctionReturn(0); 11 } 12 13 PETSC_INTERN PetscErrorCode MatShift_Basic(Mat Y,PetscScalar a) 14 { 15 PetscErrorCode ierr; 16 PetscInt i,start,end; 17 PetscScalar alpha = a; 18 PetscBool prevoption; 19 20 PetscFunctionBegin; 21 ierr = MatGetOption(Y,MAT_NO_OFF_PROC_ENTRIES,&prevoption);CHKERRQ(ierr); 22 ierr = MatSetOption(Y,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 23 ierr = MatGetOwnershipRange(Y,&start,&end);CHKERRQ(ierr); 24 for (i=start; i<end; i++) { 25 if (i < Y->cmap->N) { 26 ierr = MatSetValues(Y,1,&i,1,&i,&alpha,ADD_VALUES);CHKERRQ(ierr); 27 } 28 } 29 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 30 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 31 ierr = MatSetOption(Y,MAT_NO_OFF_PROC_ENTRIES,prevoption);CHKERRQ(ierr); 32 PetscFunctionReturn(0); 33 } 34 35 /*@ 36 MatCreate - Creates a matrix where the type is determined 37 from either a call to MatSetType() or from the options database 38 with a call to MatSetFromOptions(). The default matrix type is 39 AIJ, using the routines MatCreateSeqAIJ() or MatCreateAIJ() 40 if you do not set a type in the options database. If you never 41 call MatSetType() or MatSetFromOptions() it will generate an 42 error when you try to use the matrix. 43 44 Collective 45 46 Input Parameter: 47 . comm - MPI communicator 48 49 Output Parameter: 50 . A - the matrix 51 52 Options Database Keys: 53 + -mat_type seqaij - AIJ type, uses MatCreateSeqAIJ() 54 . -mat_type mpiaij - AIJ type, uses MatCreateAIJ() 55 . -mat_type seqdense - dense type, uses MatCreateSeqDense() 56 . -mat_type mpidense - dense type, uses MatCreateDense() 57 . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ() 58 - -mat_type mpibaij - block AIJ type, uses MatCreateBAIJ() 59 60 Even More Options Database Keys: 61 See the manpages for particular formats (e.g., MatCreateSeqAIJ()) 62 for additional format-specific options. 63 64 Level: beginner 65 66 .seealso: MatCreateSeqAIJ(), MatCreateAIJ(), 67 MatCreateSeqDense(), MatCreateDense(), 68 MatCreateSeqBAIJ(), MatCreateBAIJ(), 69 MatCreateSeqSBAIJ(), MatCreateSBAIJ(), 70 MatConvert() 71 @*/ 72 PetscErrorCode MatCreate(MPI_Comm comm,Mat *A) 73 { 74 Mat B; 75 PetscErrorCode ierr; 76 77 PetscFunctionBegin; 78 PetscValidPointer(A,2); 79 80 *A = NULL; 81 ierr = MatInitializePackage();CHKERRQ(ierr); 82 83 ierr = PetscHeaderCreate(B,MAT_CLASSID,"Mat","Matrix","Mat",comm,MatDestroy,MatView);CHKERRQ(ierr); 84 ierr = PetscLayoutCreate(comm,&B->rmap);CHKERRQ(ierr); 85 ierr = PetscLayoutCreate(comm,&B->cmap);CHKERRQ(ierr); 86 ierr = PetscStrallocpy(VECSTANDARD,&B->defaultvectype);CHKERRQ(ierr); 87 88 B->congruentlayouts = PETSC_DECIDE; 89 B->preallocated = PETSC_FALSE; 90 *A = B; 91 PetscFunctionReturn(0); 92 } 93 94 /*@ 95 MatSetErrorIfFailure - Causes Mat to generate an error, for example a zero pivot, is detected. 96 97 Logically Collective on Mat 98 99 Input Parameters: 100 + mat - matrix obtained from MatCreate() 101 - flg - PETSC_TRUE indicates you want the error generated 102 103 Level: advanced 104 105 .seealso: PCSetErrorIfFailure() 106 @*/ 107 PetscErrorCode MatSetErrorIfFailure(Mat mat,PetscBool flg) 108 { 109 PetscFunctionBegin; 110 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 111 PetscValidLogicalCollectiveBool(mat,flg,2); 112 mat->erroriffailure = flg; 113 PetscFunctionReturn(0); 114 } 115 116 /*@ 117 MatSetSizes - Sets the local and global sizes, and checks to determine compatibility 118 119 Collective on Mat 120 121 Input Parameters: 122 + A - the matrix 123 . m - number of local rows (or PETSC_DECIDE) 124 . n - number of local columns (or PETSC_DECIDE) 125 . M - number of global rows (or PETSC_DETERMINE) 126 - N - number of global columns (or PETSC_DETERMINE) 127 128 Notes: 129 m (n) and M (N) cannot be both PETSC_DECIDE 130 If one processor calls this with M (N) of PETSC_DECIDE then all processors must, otherwise the program will hang. 131 132 If PETSC_DECIDE is not used for the arguments 'm' and 'n', then the 133 user must ensure that they are chosen to be compatible with the 134 vectors. To do this, one first considers the matrix-vector product 135 'y = A x'. The 'm' that is used in the above routine must match the 136 local size used in the vector creation routine VecCreateMPI() for 'y'. 137 Likewise, the 'n' used must match that used as the local size in 138 VecCreateMPI() for 'x'. 139 140 You cannot change the sizes once they have been set. 141 142 The sizes must be set before MatSetUp() or MatXXXSetPreallocation() is called. 143 144 Level: beginner 145 146 .seealso: MatGetSize(), PetscSplitOwnership() 147 @*/ 148 PetscErrorCode MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N) 149 { 150 PetscFunctionBegin; 151 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 152 PetscValidLogicalCollectiveInt(A,M,4); 153 PetscValidLogicalCollectiveInt(A,N,5); 154 if (M > 0 && m > M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local row size %D cannot be larger than global row size %D",m,M); 155 if (N > 0 && n > N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local column size %D cannot be larger than global column size %D",n,N); 156 if ((A->rmap->n >= 0 && A->rmap->N >= 0) && (A->rmap->n != m || (M > 0 && 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); 157 if ((A->cmap->n >= 0 && A->cmap->N >= 0) && (A->cmap->n != n || (N > 0 && 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); 158 A->rmap->n = m; 159 A->cmap->n = n; 160 A->rmap->N = M > -1 ? M : A->rmap->N; 161 A->cmap->N = N > -1 ? N : A->cmap->N; 162 PetscFunctionReturn(0); 163 } 164 165 /*@ 166 MatSetFromOptions - Creates a matrix where the type is determined 167 from the options database. Generates a parallel MPI matrix if the 168 communicator has more than one processor. The default matrix type is 169 AIJ, using the routines MatCreateSeqAIJ() and MatCreateAIJ() if 170 you do not select a type in the options database. 171 172 Collective on Mat 173 174 Input Parameter: 175 . A - the matrix 176 177 Options Database Keys: 178 + -mat_type seqaij - AIJ type, uses MatCreateSeqAIJ() 179 . -mat_type mpiaij - AIJ type, uses MatCreateAIJ() 180 . -mat_type seqdense - dense type, uses MatCreateSeqDense() 181 . -mat_type mpidense - dense type, uses MatCreateDense() 182 . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ() 183 - -mat_type mpibaij - block AIJ type, uses MatCreateBAIJ() 184 185 Even More Options Database Keys: 186 See the manpages for particular formats (e.g., MatCreateSeqAIJ()) 187 for additional format-specific options. 188 189 Level: beginner 190 191 .seealso: MatCreateSeqAIJ((), MatCreateAIJ(), 192 MatCreateSeqDense(), MatCreateDense(), 193 MatCreateSeqBAIJ(), MatCreateBAIJ(), 194 MatCreateSeqSBAIJ(), MatCreateSBAIJ(), 195 MatConvert() 196 @*/ 197 PetscErrorCode MatSetFromOptions(Mat B) 198 { 199 PetscErrorCode ierr; 200 const char *deft = MATAIJ; 201 char type[256]; 202 PetscBool flg,set; 203 204 PetscFunctionBegin; 205 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 206 207 ierr = PetscObjectOptionsBegin((PetscObject)B);CHKERRQ(ierr); 208 209 if (B->rmap->bs < 0) { 210 PetscInt newbs = -1; 211 ierr = PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSetBlockSize",newbs,&newbs,&flg);CHKERRQ(ierr); 212 if (flg) { 213 ierr = PetscLayoutSetBlockSize(B->rmap,newbs);CHKERRQ(ierr); 214 ierr = PetscLayoutSetBlockSize(B->cmap,newbs);CHKERRQ(ierr); 215 } 216 } 217 218 ierr = PetscOptionsFList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);CHKERRQ(ierr); 219 if (flg) { 220 ierr = MatSetType(B,type);CHKERRQ(ierr); 221 } else if (!((PetscObject)B)->type_name) { 222 ierr = MatSetType(B,deft);CHKERRQ(ierr); 223 } 224 225 ierr = PetscOptionsName("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",&B->checksymmetryonassembly);CHKERRQ(ierr); 226 ierr = PetscOptionsReal("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",B->checksymmetrytol,&B->checksymmetrytol,NULL);CHKERRQ(ierr); 227 ierr = PetscOptionsBool("-mat_null_space_test","Checks if provided null space is correct in MatAssemblyEnd()","MatSetNullSpaceTest",B->checknullspaceonassembly,&B->checknullspaceonassembly,NULL);CHKERRQ(ierr); 228 ierr = PetscOptionsBool("-mat_error_if_failure","Generate an error if an error occurs when factoring the matrix","MatSetErrorIfFailure",B->erroriffailure,&B->erroriffailure,NULL);CHKERRQ(ierr); 229 230 if (B->ops->setfromoptions) { 231 ierr = (*B->ops->setfromoptions)(PetscOptionsObject,B);CHKERRQ(ierr); 232 } 233 234 flg = PETSC_FALSE; 235 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); 236 if (set) {ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,flg);CHKERRQ(ierr);} 237 flg = PETSC_FALSE; 238 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); 239 if (set) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,flg);CHKERRQ(ierr);} 240 241 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 242 ierr = PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)B);CHKERRQ(ierr); 243 ierr = PetscOptionsEnd();CHKERRQ(ierr); 244 PetscFunctionReturn(0); 245 } 246 247 /*@C 248 MatXAIJSetPreallocation - set preallocation for serial and parallel AIJ, BAIJ, and SBAIJ matrices and their unassembled versions. 249 250 Collective on Mat 251 252 Input Arguments: 253 + A - matrix being preallocated 254 . bs - block size 255 . dnnz - number of nonzero column blocks per block row of diagonal part of parallel matrix 256 . onnz - number of nonzero column blocks per block row of off-diagonal part of parallel matrix 257 . dnnzu - number of nonzero column blocks per block row of upper-triangular part of diagonal part of parallel matrix 258 - onnzu - number of nonzero column blocks per block row of upper-triangular part of off-diagonal part of parallel matrix 259 260 Level: beginner 261 262 .seealso: MatSeqAIJSetPreallocation(), MatMPIAIJSetPreallocation(), MatSeqBAIJSetPreallocation(), MatMPIBAIJSetPreallocation(), MatSeqSBAIJSetPreallocation(), MatMPISBAIJSetPreallocation(), 263 PetscSplitOwnership() 264 @*/ 265 PetscErrorCode MatXAIJSetPreallocation(Mat A,PetscInt bs,const PetscInt dnnz[],const PetscInt onnz[],const PetscInt dnnzu[],const PetscInt onnzu[]) 266 { 267 PetscErrorCode ierr; 268 PetscInt cbs; 269 void (*aij)(void); 270 void (*is)(void); 271 void (*hyp)(void) = NULL; 272 273 PetscFunctionBegin; 274 if (bs != PETSC_DECIDE) { /* don't mess with an already set block size */ 275 ierr = MatSetBlockSize(A,bs);CHKERRQ(ierr); 276 } 277 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 278 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 279 ierr = MatGetBlockSizes(A,&bs,&cbs);CHKERRQ(ierr); 280 /* these routines assumes bs == cbs, this should be checked somehow */ 281 ierr = MatSeqBAIJSetPreallocation(A,bs,0,dnnz);CHKERRQ(ierr); 282 ierr = MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);CHKERRQ(ierr); 283 ierr = MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);CHKERRQ(ierr); 284 ierr = MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);CHKERRQ(ierr); 285 /* 286 In general, we have to do extra work to preallocate for scalar (AIJ) or unassembled (IS) matrices so we check whether it will do any 287 good before going on with it. 288 */ 289 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 290 ierr = PetscObjectQueryFunction((PetscObject)A,"MatISSetPreallocation_C",&is);CHKERRQ(ierr); 291 #if defined(PETSC_HAVE_HYPRE) 292 ierr = PetscObjectQueryFunction((PetscObject)A,"MatHYPRESetPreallocation_C",&hyp);CHKERRQ(ierr); 293 #endif 294 if (!aij && !is && !hyp) { 295 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 296 } 297 if (aij || is || hyp) { 298 if (bs == cbs && bs == 1) { 299 ierr = MatSeqAIJSetPreallocation(A,0,dnnz);CHKERRQ(ierr); 300 ierr = MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);CHKERRQ(ierr); 301 ierr = MatISSetPreallocation(A,0,dnnz,0,onnz);CHKERRQ(ierr); 302 #if defined(PETSC_HAVE_HYPRE) 303 ierr = MatHYPRESetPreallocation(A,0,dnnz,0,onnz);CHKERRQ(ierr); 304 #endif 305 } else { /* Convert block-row precallocation to scalar-row */ 306 PetscInt i,m,*sdnnz,*sonnz; 307 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 308 ierr = PetscMalloc2((!!dnnz)*m,&sdnnz,(!!onnz)*m,&sonnz);CHKERRQ(ierr); 309 for (i=0; i<m; i++) { 310 if (dnnz) sdnnz[i] = dnnz[i/bs] * cbs; 311 if (onnz) sonnz[i] = onnz[i/bs] * cbs; 312 } 313 ierr = MatSeqAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL);CHKERRQ(ierr); 314 ierr = MatMPIAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);CHKERRQ(ierr); 315 ierr = MatISSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);CHKERRQ(ierr); 316 #if defined(PETSC_HAVE_HYPRE) 317 ierr = MatHYPRESetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);CHKERRQ(ierr); 318 #endif 319 ierr = PetscFree2(sdnnz,sonnz);CHKERRQ(ierr); 320 } 321 } 322 PetscFunctionReturn(0); 323 } 324 325 /* 326 Merges some information from Cs header to A; the C object is then destroyed 327 328 This is somewhat different from MatHeaderReplace() it would be nice to merge the code 329 */ 330 PetscErrorCode MatHeaderMerge(Mat A,Mat *C) 331 { 332 PetscErrorCode ierr; 333 PetscInt refct; 334 PetscOps Abops; 335 struct _MatOps Aops; 336 char *mtype,*mname,*mprefix; 337 Mat_Product *product; 338 339 PetscFunctionBegin; 340 /* save the parts of A we need */ 341 Abops = ((PetscObject)A)->bops[0]; 342 Aops = A->ops[0]; 343 refct = ((PetscObject)A)->refct; 344 mtype = ((PetscObject)A)->type_name; 345 mname = ((PetscObject)A)->name; 346 mprefix = ((PetscObject)A)->prefix; 347 product = A->product; 348 349 /* zero these so the destroy below does not free them */ 350 ((PetscObject)A)->type_name = 0; 351 ((PetscObject)A)->name = 0; 352 353 /* free all the interior data structures from mat */ 354 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 355 356 ierr = PetscFree(A->defaultvectype);CHKERRQ(ierr); 357 ierr = PetscLayoutDestroy(&A->rmap);CHKERRQ(ierr); 358 ierr = PetscLayoutDestroy(&A->cmap);CHKERRQ(ierr); 359 ierr = PetscFunctionListDestroy(&((PetscObject)A)->qlist);CHKERRQ(ierr); 360 ierr = PetscObjectListDestroy(&((PetscObject)A)->olist);CHKERRQ(ierr); 361 362 /* copy C over to A */ 363 ierr = PetscMemcpy(A,*C,sizeof(struct _p_Mat));CHKERRQ(ierr); 364 365 /* return the parts of A we saved */ 366 ((PetscObject)A)->bops[0] = Abops; 367 A->ops[0] = Aops; 368 ((PetscObject)A)->refct = refct; 369 ((PetscObject)A)->type_name = mtype; 370 ((PetscObject)A)->name = mname; 371 ((PetscObject)A)->prefix = mprefix; 372 A->product = product; 373 374 /* since these two are copied into A we do not want them destroyed in C */ 375 ((PetscObject)*C)->qlist = 0; 376 ((PetscObject)*C)->olist = 0; 377 378 ierr = PetscHeaderDestroy(C);CHKERRQ(ierr); 379 PetscFunctionReturn(0); 380 } 381 /* 382 Replace A's header with that of C; the C object is then destroyed 383 384 This is essentially code moved from MatDestroy() 385 386 This is somewhat different from MatHeaderMerge() it would be nice to merge the code 387 388 Used in DM hence is declared PETSC_EXTERN 389 */ 390 PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A,Mat *C) 391 { 392 PetscErrorCode ierr; 393 PetscInt refct; 394 PetscObjectState state; 395 struct _p_Mat buffer; 396 397 PetscFunctionBegin; 398 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 399 PetscValidHeaderSpecific(*C,MAT_CLASSID,2); 400 if (A == *C) PetscFunctionReturn(0); 401 PetscCheckSameComm(A,1,*C,2); 402 if (((PetscObject)*C)->refct != 1) SETERRQ1(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Object C has refct %D > 1, would leave hanging reference",((PetscObject)*C)->refct); 403 404 /* swap C and A */ 405 refct = ((PetscObject)A)->refct; 406 state = ((PetscObject)A)->state; 407 ierr = PetscMemcpy(&buffer,A,sizeof(struct _p_Mat));CHKERRQ(ierr); 408 ierr = PetscMemcpy(A,*C,sizeof(struct _p_Mat));CHKERRQ(ierr); 409 ierr = PetscMemcpy(*C,&buffer,sizeof(struct _p_Mat));CHKERRQ(ierr); 410 ((PetscObject)A)->refct = refct; 411 ((PetscObject)A)->state = state + 1; 412 413 ((PetscObject)*C)->refct = 1; 414 ierr = MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))NULL);CHKERRQ(ierr); 415 ierr = MatDestroy(C);CHKERRQ(ierr); 416 PetscFunctionReturn(0); 417 } 418 419 /*@ 420 MatBindToCPU - marks a matrix to temporarily stay on the CPU and perform computations on the CPU 421 422 Input Parameters: 423 + A - the matrix 424 - flg - bind to the CPU if value of PETSC_TRUE 425 426 Level: intermediate 427 @*/ 428 PetscErrorCode MatBindToCPU(Mat A,PetscBool flg) 429 { 430 #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) 431 PetscErrorCode ierr; 432 433 PetscFunctionBegin; 434 if (A->boundtocpu == flg) PetscFunctionReturn(0); 435 A->boundtocpu = flg; 436 if (A->ops->bindtocpu) { 437 ierr = (*A->ops->bindtocpu)(A,flg);CHKERRQ(ierr); 438 } 439 PetscFunctionReturn(0); 440 #else 441 return 0; 442 #endif 443 } 444