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