1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "MatShift_Basic" 6 PETSC_INTERN PetscErrorCode MatShift_Basic(Mat Y,PetscScalar a) 7 { 8 PetscErrorCode ierr; 9 PetscInt i,start,end; 10 PetscScalar alpha = a; 11 PetscBool prevoption; 12 13 PetscFunctionBegin; 14 ierr = MatGetOption(Y,MAT_NO_OFF_PROC_ENTRIES,&prevoption);CHKERRQ(ierr); 15 ierr = MatSetOption(Y,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 16 ierr = MatGetOwnershipRange(Y,&start,&end);CHKERRQ(ierr); 17 for (i=start; i<end; i++) { 18 ierr = MatSetValues(Y,1,&i,1,&i,&alpha,ADD_VALUES);CHKERRQ(ierr); 19 } 20 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 21 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 22 ierr = MatSetOption(Y,MAT_NO_OFF_PROC_ENTRIES,prevoption);CHKERRQ(ierr); 23 PetscFunctionReturn(0); 24 } 25 26 #undef __FUNCT__ 27 #define __FUNCT__ "MatCreate" 28 /*@ 29 MatCreate - Creates a matrix where the type is determined 30 from either a call to MatSetType() or from the options database 31 with a call to MatSetFromOptions(). The default matrix type is 32 AIJ, using the routines MatCreateSeqAIJ() or MatCreateAIJ() 33 if you do not set a type in the options database. If you never 34 call MatSetType() or MatSetFromOptions() it will generate an 35 error when you try to use the matrix. 36 37 Collective on MPI_Comm 38 39 Input Parameter: 40 . comm - MPI communicator 41 42 Output Parameter: 43 . A - the matrix 44 45 Options Database Keys: 46 + -mat_type seqaij - AIJ type, uses MatCreateSeqAIJ() 47 . -mat_type mpiaij - AIJ type, uses MatCreateAIJ() 48 . -mat_type seqdense - dense type, uses MatCreateSeqDense() 49 . -mat_type mpidense - dense type, uses MatCreateDense() 50 . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ() 51 - -mat_type mpibaij - block AIJ type, uses MatCreateBAIJ() 52 53 Even More Options Database Keys: 54 See the manpages for particular formats (e.g., MatCreateSeqAIJ()) 55 for additional format-specific options. 56 57 Notes: 58 59 Level: beginner 60 61 User manual sections: 62 + sec_matcreate 63 - chapter_matrices 64 65 .keywords: matrix, create 66 67 .seealso: MatCreateSeqAIJ(), MatCreateAIJ(), 68 MatCreateSeqDense(), MatCreateDense(), 69 MatCreateSeqBAIJ(), MatCreateBAIJ(), 70 MatCreateSeqSBAIJ(), MatCreateSBAIJ(), 71 MatConvert() 72 @*/ 73 PetscErrorCode MatCreate(MPI_Comm comm,Mat *A) 74 { 75 Mat B; 76 PetscErrorCode ierr; 77 78 PetscFunctionBegin; 79 PetscValidPointer(A,2); 80 81 *A = NULL; 82 ierr = MatInitializePackage();CHKERRQ(ierr); 83 84 ierr = PetscHeaderCreate(B,MAT_CLASSID,"Mat","Matrix","Mat",comm,MatDestroy,MatView);CHKERRQ(ierr); 85 ierr = PetscLayoutCreate(comm,&B->rmap);CHKERRQ(ierr); 86 ierr = PetscLayoutCreate(comm,&B->cmap);CHKERRQ(ierr); 87 88 B->preallocated = PETSC_FALSE; 89 *A = B; 90 PetscFunctionReturn(0); 91 } 92 93 #undef __FUNCT__ 94 #define __FUNCT__ "MatSetErrorIfFailure" 95 /*@ 96 MatSetErrorIfFailure - Causes Mat to generate an error, for example a zero pivot, is detected. 97 98 Logically Collective on Mat 99 100 Input Parameters: 101 + mat - matrix obtained from MatCreate() 102 - flg - PETSC_TRUE indicates you want the error generated 103 104 Level: advanced 105 106 .keywords: Mat, set, initial guess, nonzero 107 108 .seealso: PCSetErrorIfFailure() 109 @*/ 110 PetscErrorCode MatSetErrorIfFailure(Mat mat,PetscBool flg) 111 { 112 PetscFunctionBegin; 113 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 114 PetscValidLogicalCollectiveBool(mat,flg,2); 115 mat->erroriffailure = flg; 116 PetscFunctionReturn(0); 117 } 118 119 #undef __FUNCT__ 120 #define __FUNCT__ "MatSetSizes" 121 /*@ 122 MatSetSizes - Sets the local and global sizes, and checks to determine compatibility 123 124 Collective on Mat 125 126 Input Parameters: 127 + A - the matrix 128 . m - number of local rows (or PETSC_DECIDE) 129 . n - number of local columns (or PETSC_DECIDE) 130 . M - number of global rows (or PETSC_DETERMINE) 131 - N - number of global columns (or PETSC_DETERMINE) 132 133 Notes: 134 m (n) and M (N) cannot be both PETSC_DECIDE 135 If one processor calls this with M (N) of PETSC_DECIDE then all processors must, otherwise the program will hang. 136 137 If PETSC_DECIDE is not used for the arguments 'm' and 'n', then the 138 user must ensure that they are chosen to be compatible with the 139 vectors. To do this, one first considers the matrix-vector product 140 'y = A x'. The 'm' that is used in the above routine must match the 141 local size used in the vector creation routine VecCreateMPI() for 'y'. 142 Likewise, the 'n' used must match that used as the local size in 143 VecCreateMPI() for 'x'. 144 145 You cannot change the sizes once they have been set. 146 147 The sizes must be set before MatSetUp() or MatXXXSetPreallocation() is called. 148 149 Level: beginner 150 151 .seealso: MatGetSize(), PetscSplitOwnership() 152 @*/ 153 PetscErrorCode MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N) 154 { 155 PetscFunctionBegin; 156 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 157 if (M > 0) PetscValidLogicalCollectiveInt(A,M,4); 158 if (N > 0) PetscValidLogicalCollectiveInt(A,N,5); 159 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); 160 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); 161 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); 162 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); 163 A->rmap->n = m; 164 A->cmap->n = n; 165 A->rmap->N = M; 166 A->cmap->N = N; 167 PetscFunctionReturn(0); 168 } 169 170 #undef __FUNCT__ 171 #define __FUNCT__ "MatSetFromOptions" 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 #undef __FUNCT__ 256 #define __FUNCT__ "MatXAIJSetPreallocation" 257 /*@ 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 #undef __FUNCT__ 327 #define __FUNCT__ "MatHeaderMerge" 328 PetscErrorCode MatHeaderMerge(Mat A,Mat *C) 329 { 330 PetscErrorCode ierr; 331 PetscInt refct; 332 PetscOps Abops; 333 struct _MatOps Aops; 334 char *mtype,*mname; 335 void *spptr; 336 337 PetscFunctionBegin; 338 /* save the parts of A we need */ 339 Abops = ((PetscObject)A)->bops[0]; 340 Aops = A->ops[0]; 341 refct = ((PetscObject)A)->refct; 342 mtype = ((PetscObject)A)->type_name; 343 mname = ((PetscObject)A)->name; 344 spptr = A->spptr; 345 346 /* zero these so the destroy below does not free them */ 347 ((PetscObject)A)->type_name = 0; 348 ((PetscObject)A)->name = 0; 349 350 /* free all the interior data structures from mat */ 351 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 352 353 ierr = PetscFree((*C)->spptr);CHKERRQ(ierr); 354 355 ierr = PetscLayoutDestroy(&A->rmap);CHKERRQ(ierr); 356 ierr = PetscLayoutDestroy(&A->cmap);CHKERRQ(ierr); 357 ierr = PetscFunctionListDestroy(&((PetscObject)A)->qlist);CHKERRQ(ierr); 358 ierr = PetscObjectListDestroy(&((PetscObject)A)->olist);CHKERRQ(ierr); 359 360 /* copy C over to A */ 361 ierr = PetscMemcpy(A,*C,sizeof(struct _p_Mat));CHKERRQ(ierr); 362 363 /* return the parts of A we saved */ 364 ((PetscObject)A)->bops[0] = Abops; 365 A->ops[0] = Aops; 366 ((PetscObject)A)->refct = refct; 367 ((PetscObject)A)->type_name = mtype; 368 ((PetscObject)A)->name = mname; 369 A->spptr = spptr; 370 371 /* since these two are copied into A we do not want them destroyed in C */ 372 ((PetscObject)*C)->qlist = 0; 373 ((PetscObject)*C)->olist = 0; 374 375 ierr = PetscHeaderDestroy(C);CHKERRQ(ierr); 376 PetscFunctionReturn(0); 377 } 378 /* 379 Replace A's header with that of C; the C object is then destroyed 380 381 This is essentially code moved from MatDestroy() 382 383 This is somewhat different from MatHeaderMerge() it would be nice to merge the code 384 385 Used in DM hence is declared PETSC_EXTERN 386 */ 387 #undef __FUNCT__ 388 #define __FUNCT__ "MatHeaderReplace" 389 PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A,Mat *C) 390 { 391 PetscErrorCode ierr; 392 PetscInt refct; 393 PetscObjectState state; 394 struct _p_Mat buffer; 395 396 PetscFunctionBegin; 397 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 398 PetscValidHeaderSpecific(*C,MAT_CLASSID,2); 399 if (A == *C) PetscFunctionReturn(0); 400 PetscCheckSameComm(A,1,*C,2); 401 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); 402 403 /* swap C and A */ 404 refct = ((PetscObject)A)->refct; 405 state = ((PetscObject)A)->state; 406 ierr = PetscMemcpy(&buffer,A,sizeof(struct _p_Mat));CHKERRQ(ierr); 407 ierr = PetscMemcpy(A,*C,sizeof(struct _p_Mat));CHKERRQ(ierr); 408 ierr = PetscMemcpy(*C,&buffer,sizeof(struct _p_Mat));CHKERRQ(ierr); 409 ((PetscObject)A)->refct = refct; 410 ((PetscObject)A)->state = state + 1; 411 412 ((PetscObject)*C)->refct = 1; 413 ierr = MatDestroy(C);CHKERRQ(ierr); 414 PetscFunctionReturn(0); 415 } 416