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