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