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