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