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