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