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 ierr = MatSetValues(Y,1,&i,1,&i,&alpha,ADD_VALUES);CHKERRQ(ierr); 26 } 27 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 28 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 29 ierr = MatSetOption(Y,MAT_NO_OFF_PROC_ENTRIES,prevoption);CHKERRQ(ierr); 30 PetscFunctionReturn(0); 31 } 32 33 /*@ 34 MatCreate - Creates a matrix where the type is determined 35 from either a call to MatSetType() or from the options database 36 with a call to MatSetFromOptions(). The default matrix type is 37 AIJ, using the routines MatCreateSeqAIJ() or MatCreateAIJ() 38 if you do not set a type in the options database. If you never 39 call MatSetType() or MatSetFromOptions() it will generate an 40 error when you try to use the matrix. 41 42 Collective on MPI_Comm 43 44 Input Parameter: 45 . comm - MPI communicator 46 47 Output Parameter: 48 . A - the matrix 49 50 Options Database Keys: 51 + -mat_type seqaij - AIJ type, uses MatCreateSeqAIJ() 52 . -mat_type mpiaij - AIJ type, uses MatCreateAIJ() 53 . -mat_type seqdense - dense type, uses MatCreateSeqDense() 54 . -mat_type mpidense - dense type, uses MatCreateDense() 55 . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ() 56 - -mat_type mpibaij - block AIJ type, uses MatCreateBAIJ() 57 58 Even More Options Database Keys: 59 See the manpages for particular formats (e.g., MatCreateSeqAIJ()) 60 for additional format-specific options. 61 62 Notes: 63 64 Level: beginner 65 66 User manual sections: 67 + sec_matcreate 68 - chapter_matrices 69 70 .keywords: matrix, create 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 .keywords: Mat, set, initial guess, nonzero 112 113 .seealso: PCSetErrorIfFailure() 114 @*/ 115 PetscErrorCode MatSetErrorIfFailure(Mat mat,PetscBool flg) 116 { 117 PetscFunctionBegin; 118 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 119 PetscValidLogicalCollectiveBool(mat,flg,2); 120 mat->erroriffailure = flg; 121 PetscFunctionReturn(0); 122 } 123 124 /*@ 125 MatSetSizes - Sets the local and global sizes, and checks to determine compatibility 126 127 Collective on Mat 128 129 Input Parameters: 130 + A - the matrix 131 . m - number of local rows (or PETSC_DECIDE) 132 . n - number of local columns (or PETSC_DECIDE) 133 . M - number of global rows (or PETSC_DETERMINE) 134 - N - number of global columns (or PETSC_DETERMINE) 135 136 Notes: 137 m (n) and M (N) cannot be both PETSC_DECIDE 138 If one processor calls this with M (N) of PETSC_DECIDE then all processors must, otherwise the program will hang. 139 140 If PETSC_DECIDE is not used for the arguments 'm' and 'n', then the 141 user must ensure that they are chosen to be compatible with the 142 vectors. To do this, one first considers the matrix-vector product 143 'y = A x'. The 'm' that is used in the above routine must match the 144 local size used in the vector creation routine VecCreateMPI() for 'y'. 145 Likewise, the 'n' used must match that used as the local size in 146 VecCreateMPI() for 'x'. 147 148 You cannot change the sizes once they have been set. 149 150 The sizes must be set before MatSetUp() or MatXXXSetPreallocation() is called. 151 152 Level: beginner 153 154 .seealso: MatGetSize(), PetscSplitOwnership() 155 @*/ 156 PetscErrorCode MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N) 157 { 158 PetscFunctionBegin; 159 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 160 if (M > 0) PetscValidLogicalCollectiveInt(A,M,4); 161 if (N > 0) PetscValidLogicalCollectiveInt(A,N,5); 162 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); 163 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); 164 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); 165 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); 166 A->rmap->n = m; 167 A->cmap->n = n; 168 A->rmap->N = M > -1 ? M : A->rmap->N; 169 A->cmap->N = N > -1 ? N : A->cmap->N; 170 PetscFunctionReturn(0); 171 } 172 173 /*@ 174 MatSetFromOptions - Creates a matrix where the type is determined 175 from the options database. Generates a parallel MPI matrix if the 176 communicator has more than one processor. The default matrix type is 177 AIJ, using the routines MatCreateSeqAIJ() and MatCreateAIJ() if 178 you do not select a type in the options database. 179 180 Collective on Mat 181 182 Input Parameter: 183 . A - the matrix 184 185 Options Database Keys: 186 + -mat_type seqaij - AIJ type, uses MatCreateSeqAIJ() 187 . -mat_type mpiaij - AIJ type, uses MatCreateAIJ() 188 . -mat_type seqdense - dense type, uses MatCreateSeqDense() 189 . -mat_type mpidense - dense type, uses MatCreateDense() 190 . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ() 191 - -mat_type mpibaij - block AIJ type, uses MatCreateBAIJ() 192 193 Even More Options Database Keys: 194 See the manpages for particular formats (e.g., MatCreateSeqAIJ()) 195 for additional format-specific options. 196 197 Level: beginner 198 199 .keywords: matrix, create 200 201 .seealso: MatCreateSeqAIJ((), MatCreateAIJ(), 202 MatCreateSeqDense(), MatCreateDense(), 203 MatCreateSeqBAIJ(), MatCreateBAIJ(), 204 MatCreateSeqSBAIJ(), MatCreateSBAIJ(), 205 MatConvert() 206 @*/ 207 PetscErrorCode MatSetFromOptions(Mat B) 208 { 209 PetscErrorCode ierr; 210 const char *deft = MATAIJ; 211 char type[256]; 212 PetscBool flg,set; 213 214 PetscFunctionBegin; 215 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 216 217 ierr = PetscObjectOptionsBegin((PetscObject)B);CHKERRQ(ierr); 218 219 if (B->rmap->bs < 0) { 220 PetscInt newbs = -1; 221 ierr = PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSetBlockSize",newbs,&newbs,&flg);CHKERRQ(ierr); 222 if (flg) { 223 ierr = PetscLayoutSetBlockSize(B->rmap,newbs);CHKERRQ(ierr); 224 ierr = PetscLayoutSetBlockSize(B->cmap,newbs);CHKERRQ(ierr); 225 } 226 } 227 228 ierr = PetscOptionsFList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);CHKERRQ(ierr); 229 if (flg) { 230 ierr = MatSetType(B,type);CHKERRQ(ierr); 231 } else if (!((PetscObject)B)->type_name) { 232 ierr = MatSetType(B,deft);CHKERRQ(ierr); 233 } 234 235 ierr = PetscOptionsName("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",&B->checksymmetryonassembly);CHKERRQ(ierr); 236 ierr = PetscOptionsReal("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",B->checksymmetrytol,&B->checksymmetrytol,NULL);CHKERRQ(ierr); 237 ierr = PetscOptionsBool("-mat_null_space_test","Checks if provided null space is correct in MatAssemblyEnd()","MatSetNullSpaceTest",B->checknullspaceonassembly,&B->checknullspaceonassembly,NULL);CHKERRQ(ierr); 238 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); 239 240 if (B->ops->setfromoptions) { 241 ierr = (*B->ops->setfromoptions)(PetscOptionsObject,B);CHKERRQ(ierr); 242 } 243 244 flg = PETSC_FALSE; 245 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); 246 if (set) {ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,flg);CHKERRQ(ierr);} 247 flg = PETSC_FALSE; 248 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); 249 if (set) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,flg);CHKERRQ(ierr);} 250 251 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 252 ierr = PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)B);CHKERRQ(ierr); 253 ierr = PetscOptionsEnd();CHKERRQ(ierr); 254 PetscFunctionReturn(0); 255 } 256 257 /*@C 258 MatXAIJSetPreallocation - set preallocation for serial and parallel AIJ, BAIJ, and SBAIJ matrices and their unassembled versions. 259 260 Collective on Mat 261 262 Input Arguments: 263 + A - matrix being preallocated 264 . bs - block size 265 . dnnz - number of nonzero blocks per block row of diagonal part of parallel matrix 266 . onnz - number of nonzero blocks per block row of off-diagonal part of parallel matrix 267 . dnnzu - number of nonzero blocks per block row of upper-triangular part of diagonal part of parallel matrix 268 - onnzu - number of nonzero blocks per block row of upper-triangular part of off-diagonal part of parallel matrix 269 270 Level: beginner 271 272 .seealso: MatSeqAIJSetPreallocation(), MatMPIAIJSetPreallocation(), MatSeqBAIJSetPreallocation(), MatMPIBAIJSetPreallocation(), MatSeqSBAIJSetPreallocation(), MatMPISBAIJSetPreallocation(), 273 PetscSplitOwnership() 274 @*/ 275 PetscErrorCode MatXAIJSetPreallocation(Mat A,PetscInt bs,const PetscInt dnnz[],const PetscInt onnz[],const PetscInt dnnzu[],const PetscInt onnzu[]) 276 { 277 PetscErrorCode ierr; 278 void (*aij)(void); 279 void (*is)(void); 280 281 PetscFunctionBegin; 282 ierr = MatSetBlockSize(A,bs);CHKERRQ(ierr); 283 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 284 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 285 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 286 ierr = MatSeqBAIJSetPreallocation(A,bs,0,dnnz);CHKERRQ(ierr); 287 ierr = MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);CHKERRQ(ierr); 288 ierr = MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);CHKERRQ(ierr); 289 ierr = MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);CHKERRQ(ierr); 290 /* 291 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 292 good before going on with it. 293 */ 294 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 295 ierr = PetscObjectQueryFunction((PetscObject)A,"MatISSetPreallocation_C",&is);CHKERRQ(ierr); 296 if (!aij && !is) { 297 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 298 } 299 if (aij || is) { 300 if (bs == 1) { 301 ierr = MatSeqAIJSetPreallocation(A,0,dnnz);CHKERRQ(ierr); 302 ierr = MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);CHKERRQ(ierr); 303 ierr = MatISSetPreallocation(A,0,dnnz,0,onnz);CHKERRQ(ierr); 304 } else { /* Convert block-row precallocation to scalar-row */ 305 PetscInt i,m,*sdnnz,*sonnz; 306 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 307 ierr = PetscMalloc2((!!dnnz)*m,&sdnnz,(!!onnz)*m,&sonnz);CHKERRQ(ierr); 308 for (i=0; i<m; i++) { 309 if (dnnz) sdnnz[i] = dnnz[i/bs] * bs; 310 if (onnz) sonnz[i] = onnz[i/bs] * bs; 311 } 312 ierr = MatSeqAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL);CHKERRQ(ierr); 313 ierr = MatMPIAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);CHKERRQ(ierr); 314 ierr = MatISSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);CHKERRQ(ierr); 315 ierr = PetscFree2(sdnnz,sonnz);CHKERRQ(ierr); 316 } 317 } 318 PetscFunctionReturn(0); 319 } 320 321 /* 322 Merges some information from Cs header to A; the C object is then destroyed 323 324 This is somewhat different from MatHeaderReplace() it would be nice to merge the code 325 */ 326 PetscErrorCode MatHeaderMerge(Mat A,Mat *C) 327 { 328 PetscErrorCode ierr; 329 PetscInt refct; 330 PetscOps Abops; 331 struct _MatOps Aops; 332 char *mtype,*mname; 333 334 PetscFunctionBegin; 335 /* save the parts of A we need */ 336 Abops = ((PetscObject)A)->bops[0]; 337 Aops = A->ops[0]; 338 refct = ((PetscObject)A)->refct; 339 mtype = ((PetscObject)A)->type_name; 340 mname = ((PetscObject)A)->name; 341 342 /* zero these so the destroy below does not free them */ 343 ((PetscObject)A)->type_name = 0; 344 ((PetscObject)A)->name = 0; 345 346 /* free all the interior data structures from mat */ 347 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 348 349 ierr = PetscFree(A->defaultvectype);CHKERRQ(ierr); 350 ierr = PetscLayoutDestroy(&A->rmap);CHKERRQ(ierr); 351 ierr = PetscLayoutDestroy(&A->cmap);CHKERRQ(ierr); 352 ierr = PetscFunctionListDestroy(&((PetscObject)A)->qlist);CHKERRQ(ierr); 353 ierr = PetscObjectListDestroy(&((PetscObject)A)->olist);CHKERRQ(ierr); 354 355 /* copy C over to A */ 356 ierr = PetscMemcpy(A,*C,sizeof(struct _p_Mat));CHKERRQ(ierr); 357 358 /* return the parts of A we saved */ 359 ((PetscObject)A)->bops[0] = Abops; 360 A->ops[0] = Aops; 361 ((PetscObject)A)->refct = refct; 362 ((PetscObject)A)->type_name = mtype; 363 ((PetscObject)A)->name = mname; 364 365 /* since these two are copied into A we do not want them destroyed in C */ 366 ((PetscObject)*C)->qlist = 0; 367 ((PetscObject)*C)->olist = 0; 368 369 ierr = PetscHeaderDestroy(C);CHKERRQ(ierr); 370 PetscFunctionReturn(0); 371 } 372 /* 373 Replace A's header with that of C; the C object is then destroyed 374 375 This is essentially code moved from MatDestroy() 376 377 This is somewhat different from MatHeaderMerge() it would be nice to merge the code 378 379 Used in DM hence is declared PETSC_EXTERN 380 */ 381 PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A,Mat *C) 382 { 383 PetscErrorCode ierr; 384 PetscInt refct; 385 PetscObjectState state; 386 struct _p_Mat buffer; 387 388 PetscFunctionBegin; 389 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 390 PetscValidHeaderSpecific(*C,MAT_CLASSID,2); 391 if (A == *C) PetscFunctionReturn(0); 392 PetscCheckSameComm(A,1,*C,2); 393 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); 394 395 /* swap C and A */ 396 refct = ((PetscObject)A)->refct; 397 state = ((PetscObject)A)->state; 398 ierr = PetscMemcpy(&buffer,A,sizeof(struct _p_Mat));CHKERRQ(ierr); 399 ierr = PetscMemcpy(A,*C,sizeof(struct _p_Mat));CHKERRQ(ierr); 400 ierr = PetscMemcpy(*C,&buffer,sizeof(struct _p_Mat));CHKERRQ(ierr); 401 ((PetscObject)A)->refct = refct; 402 ((PetscObject)A)->state = state + 1; 403 404 ((PetscObject)*C)->refct = 1; 405 ierr = MatDestroy(C);CHKERRQ(ierr); 406 PetscFunctionReturn(0); 407 } 408