1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 static PetscErrorCode MatTransposeAXPY_Private(Mat Y,PetscScalar a,Mat X,MatStructure str,Mat T) 5 { 6 PetscErrorCode ierr,(*f)(Mat,Mat*); 7 Mat A,F; 8 9 PetscFunctionBegin; 10 ierr = PetscObjectQueryFunction((PetscObject)T,"MatTransposeGetMat_C",&f);CHKERRQ(ierr); 11 if (f) { 12 ierr = MatTransposeGetMat(T,&A);CHKERRQ(ierr); 13 if (T == X) { 14 ierr = PetscInfo(NULL,"Explicitly transposing X of type MATTRANSPOSEMAT to perform MatAXPY()\n");CHKERRQ(ierr); 15 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&F);CHKERRQ(ierr); 16 A = Y; 17 } else { 18 ierr = PetscInfo(NULL,"Transposing X because Y of type MATTRANSPOSEMAT to perform MatAXPY()\n");CHKERRQ(ierr); 19 ierr = MatTranspose(X,MAT_INITIAL_MATRIX,&F);CHKERRQ(ierr); 20 } 21 } else { 22 ierr = MatHermitianTransposeGetMat(T,&A);CHKERRQ(ierr); 23 if (T == X) { 24 ierr = PetscInfo(NULL,"Explicitly Hermitian transposing X of type MATTRANSPOSEMAT to perform MatAXPY()\n");CHKERRQ(ierr); 25 ierr = MatHermitianTranspose(A,MAT_INITIAL_MATRIX,&F);CHKERRQ(ierr); 26 A = Y; 27 } else { 28 ierr = PetscInfo(NULL,"Hermitian transposing X because Y of type MATTRANSPOSEMAT to perform MatAXPY()\n");CHKERRQ(ierr); 29 ierr = MatHermitianTranspose(X,MAT_INITIAL_MATRIX,&F);CHKERRQ(ierr); 30 } 31 } 32 ierr = MatAXPY(A,a,F,str);CHKERRQ(ierr); 33 ierr = MatDestroy(&F);CHKERRQ(ierr); 34 PetscFunctionReturn(0); 35 } 36 37 /*@ 38 MatAXPY - Computes Y = a*X + Y. 39 40 Logically Collective on Mat 41 42 Input Parameters: 43 + a - the scalar multiplier 44 . X - the first matrix 45 . Y - the second matrix 46 - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN 47 or SUBSET_NONZERO_PATTERN (nonzeros of X is a subset of Y's) 48 49 Level: intermediate 50 51 .keywords: matrix, add 52 53 .seealso: MatAYPX() 54 @*/ 55 PetscErrorCode MatAXPY(Mat Y,PetscScalar a,Mat X,MatStructure str) 56 { 57 PetscErrorCode ierr; 58 PetscInt M1,M2,N1,N2; 59 PetscInt m1,m2,n1,n2; 60 MatType t1,t2; 61 PetscBool sametype,transpose; 62 63 PetscFunctionBegin; 64 PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 65 PetscValidLogicalCollectiveScalar(Y,a,2); 66 PetscValidHeaderSpecific(X,MAT_CLASSID,3); 67 PetscCheckSameComm(Y,1,X,3); 68 ierr = MatGetSize(X,&M1,&N1);CHKERRQ(ierr); 69 ierr = MatGetSize(Y,&M2,&N2);CHKERRQ(ierr); 70 ierr = MatGetLocalSize(X,&m1,&n1);CHKERRQ(ierr); 71 ierr = MatGetLocalSize(Y,&m2,&n2);CHKERRQ(ierr); 72 if (M1 != M2 || N1 != N2) SETERRQ4(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_SIZ,"Non conforming matrix add: global sizes %D x %D, %D x %D",M1,M2,N1,N2); 73 if (m1 != m2 || n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Non conforming matrix add: local sizes %D x %D, %D x %D",m1,m2,n1,n2); 74 75 ierr = MatGetType(X,&t1);CHKERRQ(ierr); 76 ierr = MatGetType(Y,&t2);CHKERRQ(ierr); 77 ierr = PetscStrcmp(t1,t2,&sametype);CHKERRQ(ierr); 78 ierr = PetscLogEventBegin(MAT_AXPY,Y,0,0,0);CHKERRQ(ierr); 79 if (Y->ops->axpy && sametype) { 80 ierr = (*Y->ops->axpy)(Y,a,X,str);CHKERRQ(ierr); 81 } else { 82 ierr = PetscStrcmp(t1,MATTRANSPOSEMAT,&transpose);CHKERRQ(ierr); 83 if (transpose) { 84 ierr = MatTransposeAXPY_Private(Y,a,X,str,X);CHKERRQ(ierr); 85 } else { 86 ierr = PetscStrcmp(t2,MATTRANSPOSEMAT,&transpose);CHKERRQ(ierr); 87 if (transpose) { 88 ierr = MatTransposeAXPY_Private(Y,a,X,str,Y);CHKERRQ(ierr); 89 } else { 90 if (str != DIFFERENT_NONZERO_PATTERN) { 91 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 92 } else { 93 Mat B; 94 95 ierr = MatAXPY_Basic_Preallocate(Y,X,&B);CHKERRQ(ierr); 96 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 97 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 98 } 99 } 100 } 101 } 102 ierr = PetscLogEventEnd(MAT_AXPY,Y,0,0,0);CHKERRQ(ierr); 103 #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) 104 if (Y->valid_GPU_matrix != PETSC_OFFLOAD_UNALLOCATED) { 105 Y->valid_GPU_matrix = PETSC_OFFLOAD_CPU; 106 } 107 #endif 108 PetscFunctionReturn(0); 109 } 110 111 PetscErrorCode MatAXPY_Basic_Preallocate(Mat Y, Mat X, Mat *B) 112 { 113 PetscErrorCode ierr; 114 PetscErrorCode (*preall)(Mat,Mat,Mat*) = NULL; 115 116 PetscFunctionBegin; 117 /* look for any available faster alternative to the general preallocator */ 118 ierr = PetscObjectQueryFunction((PetscObject)Y,"MatAXPYGetPreallocation_C",&preall);CHKERRQ(ierr); 119 if (preall) { 120 ierr = (*preall)(Y,X,B);CHKERRQ(ierr); 121 } else { /* Use MatPrellocator, assumes same row-col distribution */ 122 Mat preallocator; 123 PetscInt r,rstart,rend; 124 PetscInt m,n,M,N; 125 126 ierr = MatGetRowUpperTriangular(Y);CHKERRQ(ierr); 127 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 128 ierr = MatGetSize(Y,&M,&N);CHKERRQ(ierr); 129 ierr = MatGetLocalSize(Y,&m,&n);CHKERRQ(ierr); 130 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&preallocator);CHKERRQ(ierr); 131 ierr = MatSetType(preallocator,MATPREALLOCATOR);CHKERRQ(ierr); 132 ierr = MatSetSizes(preallocator,m,n,M,N);CHKERRQ(ierr); 133 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 134 ierr = MatGetOwnershipRange(preallocator,&rstart,&rend);CHKERRQ(ierr); 135 for (r = rstart; r < rend; ++r) { 136 PetscInt ncols; 137 const PetscInt *row; 138 const PetscScalar *vals; 139 140 ierr = MatGetRow(Y,r,&ncols,&row,&vals);CHKERRQ(ierr); 141 ierr = MatSetValues(preallocator,1,&r,ncols,row,vals,INSERT_VALUES);CHKERRQ(ierr); 142 ierr = MatRestoreRow(Y,r,&ncols,&row,&vals);CHKERRQ(ierr); 143 ierr = MatGetRow(X,r,&ncols,&row,&vals);CHKERRQ(ierr); 144 ierr = MatSetValues(preallocator,1,&r,ncols,row,vals,INSERT_VALUES);CHKERRQ(ierr); 145 ierr = MatRestoreRow(X,r,&ncols,&row,&vals);CHKERRQ(ierr); 146 } 147 ierr = MatAssemblyBegin(preallocator,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 148 ierr = MatAssemblyEnd(preallocator,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 149 ierr = MatRestoreRowUpperTriangular(Y);CHKERRQ(ierr); 150 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 151 152 ierr = MatCreate(PetscObjectComm((PetscObject)Y),B);CHKERRQ(ierr); 153 ierr = MatSetType(*B,((PetscObject)Y)->type_name);CHKERRQ(ierr); 154 ierr = MatSetSizes(*B,m,n,M,N);CHKERRQ(ierr); 155 ierr = MatPreallocatorPreallocate(preallocator,PETSC_FALSE,*B);CHKERRQ(ierr); 156 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 157 } 158 PetscFunctionReturn(0); 159 } 160 161 PetscErrorCode MatAXPY_Basic(Mat Y,PetscScalar a,Mat X,MatStructure str) 162 { 163 PetscInt i,start,end,j,ncols,m,n; 164 PetscErrorCode ierr; 165 const PetscInt *row; 166 PetscScalar *val; 167 const PetscScalar *vals; 168 169 PetscFunctionBegin; 170 ierr = MatGetSize(X,&m,&n);CHKERRQ(ierr); 171 ierr = MatGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 172 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 173 if (a == 1.0) { 174 for (i = start; i < end; i++) { 175 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 176 ierr = MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 177 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 178 } 179 } else { 180 PetscInt vs = 100; 181 /* realloc if needed, as this function may be used in parallel */ 182 ierr = PetscMalloc1(vs,&val);CHKERRQ(ierr); 183 for (i=start; i<end; i++) { 184 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 185 if (vs < ncols) { 186 vs = PetscMin(2*ncols,n); 187 ierr = PetscRealloc(vs*sizeof(*val),&val);CHKERRQ(ierr); 188 } 189 for (j=0; j<ncols; j++) val[j] = a*vals[j]; 190 ierr = MatSetValues(Y,1,&i,ncols,row,val,ADD_VALUES);CHKERRQ(ierr); 191 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 192 } 193 ierr = PetscFree(val);CHKERRQ(ierr); 194 } 195 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 196 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 197 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 198 PetscFunctionReturn(0); 199 } 200 201 PetscErrorCode MatAXPY_BasicWithPreallocation(Mat B,Mat Y,PetscScalar a,Mat X,MatStructure str) 202 { 203 PetscInt i,start,end,j,ncols,m,n; 204 PetscErrorCode ierr; 205 const PetscInt *row; 206 PetscScalar *val; 207 const PetscScalar *vals; 208 209 PetscFunctionBegin; 210 ierr = MatGetSize(X,&m,&n);CHKERRQ(ierr); 211 ierr = MatGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 212 ierr = MatGetRowUpperTriangular(Y);CHKERRQ(ierr); 213 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 214 if (a == 1.0) { 215 for (i = start; i < end; i++) { 216 ierr = MatGetRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 217 ierr = MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 218 ierr = MatRestoreRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 219 220 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 221 ierr = MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 222 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 223 } 224 } else { 225 PetscInt vs = 100; 226 /* realloc if needed, as this function may be used in parallel */ 227 ierr = PetscMalloc1(vs,&val);CHKERRQ(ierr); 228 for (i=start; i<end; i++) { 229 ierr = MatGetRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 230 ierr = MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 231 ierr = MatRestoreRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 232 233 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 234 if (vs < ncols) { 235 vs = PetscMin(2*ncols,n); 236 ierr = PetscRealloc(vs*sizeof(*val),&val);CHKERRQ(ierr); 237 } 238 for (j=0; j<ncols; j++) val[j] = a*vals[j]; 239 ierr = MatSetValues(B,1,&i,ncols,row,val,ADD_VALUES);CHKERRQ(ierr); 240 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 241 } 242 ierr = PetscFree(val);CHKERRQ(ierr); 243 } 244 ierr = MatRestoreRowUpperTriangular(Y);CHKERRQ(ierr); 245 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 246 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 247 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 248 PetscFunctionReturn(0); 249 } 250 251 /*@ 252 MatShift - Computes Y = Y + a I, where a is a PetscScalar and I is the identity matrix. 253 254 Neighbor-wise Collective on Mat 255 256 Input Parameters: 257 + Y - the matrices 258 - a - the PetscScalar 259 260 Level: intermediate 261 262 Notes: 263 If the matrix Y is missing some diagonal entries this routine can be very slow. To make it fast one should initially 264 fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an 265 entry). 266 267 To form Y = Y + diag(V) use MatDiagonalSet() 268 269 Developers Note: If the local "diagonal part" of the matrix Y has no entries then the local diagonal part is 270 preallocated with 1 nonzero per row for the to be added values. This allows for fast shifting of an empty matrix. 271 272 .keywords: matrix, add, shift 273 274 .seealso: MatDiagonalSet(), MatScale(), MatDiagonalScale() 275 @*/ 276 PetscErrorCode MatShift(Mat Y,PetscScalar a) 277 { 278 PetscErrorCode ierr; 279 280 PetscFunctionBegin; 281 PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 282 if (!Y->assembled) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 283 if (Y->factortype) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 284 MatCheckPreallocated(Y,1); 285 if (a == 0.0) PetscFunctionReturn(0); 286 287 if (Y->ops->shift) { 288 ierr = (*Y->ops->shift)(Y,a);CHKERRQ(ierr); 289 } else { 290 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 291 } 292 293 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 294 #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) 295 if (Y->valid_GPU_matrix != PETSC_OFFLOAD_UNALLOCATED) { 296 Y->valid_GPU_matrix = PETSC_OFFLOAD_CPU; 297 } 298 #endif 299 PetscFunctionReturn(0); 300 } 301 302 PetscErrorCode MatDiagonalSet_Default(Mat Y,Vec D,InsertMode is) 303 { 304 PetscErrorCode ierr; 305 PetscInt i,start,end; 306 PetscScalar *v; 307 308 PetscFunctionBegin; 309 ierr = MatGetOwnershipRange(Y,&start,&end);CHKERRQ(ierr); 310 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 311 for (i=start; i<end; i++) { 312 ierr = MatSetValues(Y,1,&i,1,&i,v+i-start,is);CHKERRQ(ierr); 313 } 314 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 315 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 316 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 317 PetscFunctionReturn(0); 318 } 319 320 /*@ 321 MatDiagonalSet - Computes Y = Y + D, where D is a diagonal matrix 322 that is represented as a vector. Or Y[i,i] = D[i] if InsertMode is 323 INSERT_VALUES. 324 325 Input Parameters: 326 + Y - the input matrix 327 . D - the diagonal matrix, represented as a vector 328 - i - INSERT_VALUES or ADD_VALUES 329 330 Neighbor-wise Collective on Mat and Vec 331 332 Notes: 333 If the matrix Y is missing some diagonal entries this routine can be very slow. To make it fast one should initially 334 fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an 335 entry). 336 337 Level: intermediate 338 339 .keywords: matrix, add, shift, diagonal 340 341 .seealso: MatShift(), MatScale(), MatDiagonalScale() 342 @*/ 343 PetscErrorCode MatDiagonalSet(Mat Y,Vec D,InsertMode is) 344 { 345 PetscErrorCode ierr; 346 PetscInt matlocal,veclocal; 347 348 PetscFunctionBegin; 349 PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 350 PetscValidHeaderSpecific(D,VEC_CLASSID,2); 351 ierr = MatGetLocalSize(Y,&matlocal,NULL);CHKERRQ(ierr); 352 ierr = VecGetLocalSize(D,&veclocal);CHKERRQ(ierr); 353 if (matlocal != veclocal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number local rows of matrix %D does not match that of vector for diagonal %D",matlocal,veclocal); 354 if (Y->ops->diagonalset) { 355 ierr = (*Y->ops->diagonalset)(Y,D,is);CHKERRQ(ierr); 356 } else { 357 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 358 } 359 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 360 PetscFunctionReturn(0); 361 } 362 363 /*@ 364 MatAYPX - Computes Y = a*Y + X. 365 366 Logically on Mat 367 368 Input Parameters: 369 + a - the PetscScalar multiplier 370 . Y - the first matrix 371 . X - the second matrix 372 - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN or SUBSET_NONZERO_PATTERN 373 374 Level: intermediate 375 376 .keywords: matrix, add 377 378 .seealso: MatAXPY() 379 @*/ 380 PetscErrorCode MatAYPX(Mat Y,PetscScalar a,Mat X,MatStructure str) 381 { 382 PetscScalar one = 1.0; 383 PetscErrorCode ierr; 384 PetscInt mX,mY,nX,nY; 385 386 PetscFunctionBegin; 387 PetscValidHeaderSpecific(X,MAT_CLASSID,3); 388 PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 389 PetscValidLogicalCollectiveScalar(Y,a,2); 390 ierr = MatGetSize(X,&mX,&nX);CHKERRQ(ierr); 391 ierr = MatGetSize(X,&mY,&nY);CHKERRQ(ierr); 392 if (mX != mY || nX != nY) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Non conforming matrices: %D %D first %D %D second",mX,mY,nX,nY); 393 ierr = MatScale(Y,a);CHKERRQ(ierr); 394 ierr = MatAXPY(Y,one,X,str);CHKERRQ(ierr); 395 PetscFunctionReturn(0); 396 } 397 398 /*@ 399 MatComputeOperator - Computes the explicit matrix 400 401 Collective on Mat 402 403 Input Parameter: 404 + inmat - the matrix 405 - mattype - the matrix type for the explicit operator 406 407 Output Parameter: 408 . mat - the explict operator 409 410 Notes: 411 This computation is done by applying the operators to columns of the identity matrix. 412 This routine is costly in general, and is recommended for use only with relatively small systems. 413 Currently, this routine uses a dense matrix format if mattype == NULL. 414 415 Level: advanced 416 417 .keywords: Mat, compute, explicit, operator 418 @*/ 419 PetscErrorCode MatComputeOperator(Mat inmat,MatType mattype,Mat *mat) 420 { 421 PetscErrorCode ierr; 422 423 PetscFunctionBegin; 424 PetscValidHeaderSpecific(inmat,MAT_CLASSID,1); 425 PetscValidPointer(mat,3); 426 ierr = MatConvert_Shell(inmat,mattype ? mattype : MATDENSE,MAT_INITIAL_MATRIX,mat);CHKERRQ(ierr); 427 PetscFunctionReturn(0); 428 } 429 430 /*@ 431 MatComputeOperatorTranspose - Computes the explicit matrix representation of 432 a give matrix that can apply MatMultTranspose() 433 434 Collective on Mat 435 436 Input Parameter: 437 . inmat - the matrix 438 439 Output Parameter: 440 . mat - the explict operator transposed 441 442 Notes: 443 This computation is done by applying the transpose of the operator to columns of the identity matrix. 444 This routine is costly in general, and is recommended for use only with relatively small systems. 445 Currently, this routine uses a dense matrix format if mattype == NULL. 446 447 Level: advanced 448 449 .keywords: Mat, compute, explicit, operator 450 @*/ 451 PetscErrorCode MatComputeOperatorTranspose(Mat inmat,MatType mattype,Mat *mat) 452 { 453 Mat A; 454 PetscErrorCode ierr; 455 456 PetscFunctionBegin; 457 PetscValidHeaderSpecific(inmat,MAT_CLASSID,1); 458 PetscValidPointer(mat,3); 459 ierr = MatCreateTranspose(inmat,&A);CHKERRQ(ierr); 460 ierr = MatConvert_Shell(A,mattype ? mattype : MATDENSE,MAT_INITIAL_MATRIX,mat);CHKERRQ(ierr); 461 ierr = MatDestroy(&A);CHKERRQ(ierr); 462 PetscFunctionReturn(0); 463 } 464 465 /*@ 466 MatChop - Set all values in the matrix less than the tolerance to zero 467 468 Input Parameters: 469 + A - The matrix 470 - tol - The zero tolerance 471 472 Output Parameters: 473 . A - The chopped matrix 474 475 Level: intermediate 476 477 .seealso: MatCreate(), MatZeroEntries() 478 @*/ 479 PetscErrorCode MatChop(Mat A, PetscReal tol) 480 { 481 PetscScalar *newVals; 482 PetscInt *newCols; 483 PetscInt rStart, rEnd, numRows, maxRows, r, colMax = 0; 484 PetscErrorCode ierr; 485 486 PetscFunctionBegin; 487 ierr = MatGetOwnershipRange(A, &rStart, &rEnd);CHKERRQ(ierr); 488 for (r = rStart; r < rEnd; ++r) { 489 PetscInt ncols; 490 491 ierr = MatGetRow(A, r, &ncols, NULL, NULL);CHKERRQ(ierr); 492 colMax = PetscMax(colMax, ncols);CHKERRQ(ierr); 493 ierr = MatRestoreRow(A, r, &ncols, NULL, NULL);CHKERRQ(ierr); 494 } 495 numRows = rEnd - rStart; 496 ierr = MPIU_Allreduce(&numRows, &maxRows, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 497 ierr = PetscMalloc2(colMax,&newCols,colMax,&newVals);CHKERRQ(ierr); 498 for (r = rStart; r < rStart+maxRows; ++r) { 499 const PetscScalar *vals; 500 const PetscInt *cols; 501 PetscInt ncols, newcols, c; 502 503 if (r < rEnd) { 504 ierr = MatGetRow(A, r, &ncols, &cols, &vals);CHKERRQ(ierr); 505 for (c = 0; c < ncols; ++c) { 506 newCols[c] = cols[c]; 507 newVals[c] = PetscAbsScalar(vals[c]) < tol ? 0.0 : vals[c]; 508 } 509 newcols = ncols; 510 ierr = MatRestoreRow(A, r, &ncols, &cols, &vals);CHKERRQ(ierr); 511 ierr = MatSetValues(A, 1, &r, newcols, newCols, newVals, INSERT_VALUES);CHKERRQ(ierr); 512 } 513 ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 514 ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 515 } 516 ierr = PetscFree2(newCols,newVals);CHKERRQ(ierr); 517 PetscFunctionReturn(0); 518 } 519