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->congruentlayouts = -1; 89 B->preallocated = PETSC_FALSE; 90 *A = B; 91 PetscFunctionReturn(0); 92 } 93 94 #undef __FUNCT__ 95 #define __FUNCT__ "MatSetErrorIfFailure" 96 /*@ 97 MatSetErrorIfFailure - Causes Mat to generate an error, for example a zero pivot, is detected. 98 99 Logically Collective on Mat 100 101 Input Parameters: 102 + mat - matrix obtained from MatCreate() 103 - flg - PETSC_TRUE indicates you want the error generated 104 105 Level: advanced 106 107 .keywords: Mat, set, initial guess, nonzero 108 109 .seealso: PCSetErrorIfFailure() 110 @*/ 111 PetscErrorCode MatSetErrorIfFailure(Mat mat,PetscBool flg) 112 { 113 PetscFunctionBegin; 114 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 115 PetscValidLogicalCollectiveBool(mat,flg,2); 116 mat->erroriffailure = flg; 117 PetscFunctionReturn(0); 118 } 119 120 #undef __FUNCT__ 121 #define __FUNCT__ "MatSetSizes" 122 /*@ 123 MatSetSizes - Sets the local and global sizes, and checks to determine compatibility 124 125 Collective on Mat 126 127 Input Parameters: 128 + A - the matrix 129 . m - number of local rows (or PETSC_DECIDE) 130 . n - number of local columns (or PETSC_DECIDE) 131 . M - number of global rows (or PETSC_DETERMINE) 132 - N - number of global columns (or PETSC_DETERMINE) 133 134 Notes: 135 m (n) and M (N) cannot be both PETSC_DECIDE 136 If one processor calls this with M (N) of PETSC_DECIDE then all processors must, otherwise the program will hang. 137 138 If PETSC_DECIDE is not used for the arguments 'm' and 'n', then the 139 user must ensure that they are chosen to be compatible with the 140 vectors. To do this, one first considers the matrix-vector product 141 'y = A x'. The 'm' that is used in the above routine must match the 142 local size used in the vector creation routine VecCreateMPI() for 'y'. 143 Likewise, the 'n' used must match that used as the local size in 144 VecCreateMPI() for 'x'. 145 146 You cannot change the sizes once they have been set. 147 148 The sizes must be set before MatSetUp() or MatXXXSetPreallocation() is called. 149 150 Level: beginner 151 152 .seealso: MatGetSize(), PetscSplitOwnership() 153 @*/ 154 PetscErrorCode MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N) 155 { 156 PetscFunctionBegin; 157 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 158 if (M > 0) PetscValidLogicalCollectiveInt(A,M,4); 159 if (N > 0) PetscValidLogicalCollectiveInt(A,N,5); 160 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); 161 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); 162 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); 163 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); 164 A->rmap->n = m; 165 A->cmap->n = n; 166 A->rmap->N = M; 167 A->cmap->N = N; 168 PetscFunctionReturn(0); 169 } 170 171 #undef __FUNCT__ 172 #define __FUNCT__ "MatSetFromOptions" 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 239 if (B->ops->setfromoptions) { 240 ierr = (*B->ops->setfromoptions)(PetscOptionsObject,B);CHKERRQ(ierr); 241 } 242 243 flg = PETSC_FALSE; 244 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); 245 if (set) {ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,flg);CHKERRQ(ierr);} 246 flg = PETSC_FALSE; 247 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); 248 if (set) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,flg);CHKERRQ(ierr);} 249 250 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 251 ierr = PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)B);CHKERRQ(ierr); 252 ierr = PetscOptionsEnd();CHKERRQ(ierr); 253 PetscFunctionReturn(0); 254 } 255 256 #undef __FUNCT__ 257 #define __FUNCT__ "MatXAIJSetPreallocation" 258 /*@ 259 MatXAIJSetPreallocation - set preallocation for serial and parallel AIJ, BAIJ, and SBAIJ matrices and their unassembled versions. 260 261 Collective on Mat 262 263 Input Arguments: 264 + A - matrix being preallocated 265 . bs - block size 266 . dnnz - number of nonzero blocks per block row of diagonal part of parallel matrix 267 . onnz - number of nonzero blocks per block row of off-diagonal part of parallel matrix 268 . dnnzu - number of nonzero blocks per block row of upper-triangular part of diagonal part of parallel matrix 269 - onnzu - number of nonzero blocks per block row of upper-triangular part of off-diagonal part of parallel matrix 270 271 Level: beginner 272 273 .seealso: MatSeqAIJSetPreallocation(), MatMPIAIJSetPreallocation(), MatSeqBAIJSetPreallocation(), MatMPIBAIJSetPreallocation(), MatSeqSBAIJSetPreallocation(), MatMPISBAIJSetPreallocation(), 274 PetscSplitOwnership() 275 @*/ 276 PetscErrorCode MatXAIJSetPreallocation(Mat A,PetscInt bs,const PetscInt dnnz[],const PetscInt onnz[],const PetscInt dnnzu[],const PetscInt onnzu[]) 277 { 278 PetscErrorCode ierr; 279 void (*aij)(void); 280 void (*is)(void); 281 282 PetscFunctionBegin; 283 ierr = MatSetBlockSize(A,bs);CHKERRQ(ierr); 284 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 285 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 286 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 287 ierr = MatSeqBAIJSetPreallocation(A,bs,0,dnnz);CHKERRQ(ierr); 288 ierr = MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);CHKERRQ(ierr); 289 ierr = MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);CHKERRQ(ierr); 290 ierr = MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);CHKERRQ(ierr); 291 /* 292 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 293 good before going on with it. 294 */ 295 ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 296 ierr = PetscObjectQueryFunction((PetscObject)A,"MatISSetPreallocation_C",&is);CHKERRQ(ierr); 297 if (!aij && !is) { 298 ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 299 } 300 if (aij || is) { 301 if (bs == 1) { 302 ierr = MatSeqAIJSetPreallocation(A,0,dnnz);CHKERRQ(ierr); 303 ierr = MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);CHKERRQ(ierr); 304 ierr = MatISSetPreallocation(A,0,dnnz,0,onnz);CHKERRQ(ierr); 305 } else { /* Convert block-row precallocation to scalar-row */ 306 PetscInt i,m,*sdnnz,*sonnz; 307 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 308 ierr = PetscMalloc2((!!dnnz)*m,&sdnnz,(!!onnz)*m,&sonnz);CHKERRQ(ierr); 309 for (i=0; i<m; i++) { 310 if (dnnz) sdnnz[i] = dnnz[i/bs] * bs; 311 if (onnz) sonnz[i] = onnz[i/bs] * bs; 312 } 313 ierr = MatSeqAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL);CHKERRQ(ierr); 314 ierr = MatMPIAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);CHKERRQ(ierr); 315 ierr = MatISSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);CHKERRQ(ierr); 316 ierr = PetscFree2(sdnnz,sonnz);CHKERRQ(ierr); 317 } 318 } 319 PetscFunctionReturn(0); 320 } 321 322 /* 323 Merges some information from Cs header to A; the C object is then destroyed 324 325 This is somewhat different from MatHeaderReplace() it would be nice to merge the code 326 */ 327 #undef __FUNCT__ 328 #define __FUNCT__ "MatHeaderMerge" 329 PetscErrorCode MatHeaderMerge(Mat A,Mat *C) 330 { 331 PetscErrorCode ierr; 332 PetscInt refct; 333 PetscOps Abops; 334 struct _MatOps Aops; 335 char *mtype,*mname; 336 void *spptr; 337 338 PetscFunctionBegin; 339 /* save the parts of A we need */ 340 Abops = ((PetscObject)A)->bops[0]; 341 Aops = A->ops[0]; 342 refct = ((PetscObject)A)->refct; 343 mtype = ((PetscObject)A)->type_name; 344 mname = ((PetscObject)A)->name; 345 spptr = A->spptr; 346 347 /* zero these so the destroy below does not free them */ 348 ((PetscObject)A)->type_name = 0; 349 ((PetscObject)A)->name = 0; 350 351 /* free all the interior data structures from mat */ 352 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 353 354 ierr = PetscFree((*C)->spptr);CHKERRQ(ierr); 355 356 ierr = PetscLayoutDestroy(&A->rmap);CHKERRQ(ierr); 357 ierr = PetscLayoutDestroy(&A->cmap);CHKERRQ(ierr); 358 ierr = PetscFunctionListDestroy(&((PetscObject)A)->qlist);CHKERRQ(ierr); 359 ierr = PetscObjectListDestroy(&((PetscObject)A)->olist);CHKERRQ(ierr); 360 361 /* copy C over to A */ 362 ierr = PetscMemcpy(A,*C,sizeof(struct _p_Mat));CHKERRQ(ierr); 363 364 /* return the parts of A we saved */ 365 ((PetscObject)A)->bops[0] = Abops; 366 A->ops[0] = Aops; 367 ((PetscObject)A)->refct = refct; 368 ((PetscObject)A)->type_name = mtype; 369 ((PetscObject)A)->name = mname; 370 A->spptr = spptr; 371 372 /* since these two are copied into A we do not want them destroyed in C */ 373 ((PetscObject)*C)->qlist = 0; 374 ((PetscObject)*C)->olist = 0; 375 376 ierr = PetscHeaderDestroy(C);CHKERRQ(ierr); 377 PetscFunctionReturn(0); 378 } 379 /* 380 Replace A's header with that of C; the C object is then destroyed 381 382 This is essentially code moved from MatDestroy() 383 384 This is somewhat different from MatHeaderMerge() it would be nice to merge the code 385 386 Used in DM hence is declared PETSC_EXTERN 387 */ 388 #undef __FUNCT__ 389 #define __FUNCT__ "MatHeaderReplace" 390 PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A,Mat *C) 391 { 392 PetscErrorCode ierr; 393 PetscInt refct; 394 PetscObjectState state; 395 struct _p_Mat buffer; 396 397 PetscFunctionBegin; 398 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 399 PetscValidHeaderSpecific(*C,MAT_CLASSID,2); 400 if (A == *C) PetscFunctionReturn(0); 401 PetscCheckSameComm(A,1,*C,2); 402 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); 403 404 /* swap C and A */ 405 refct = ((PetscObject)A)->refct; 406 state = ((PetscObject)A)->state; 407 ierr = PetscMemcpy(&buffer,A,sizeof(struct _p_Mat));CHKERRQ(ierr); 408 ierr = PetscMemcpy(A,*C,sizeof(struct _p_Mat));CHKERRQ(ierr); 409 ierr = PetscMemcpy(*C,&buffer,sizeof(struct _p_Mat));CHKERRQ(ierr); 410 ((PetscObject)A)->refct = refct; 411 ((PetscObject)A)->state = state + 1; 412 413 ((PetscObject)*C)->refct = 1; 414 ierr = MatDestroy(C);CHKERRQ(ierr); 415 PetscFunctionReturn(0); 416 } 417