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