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