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 Notes: No operation is performed when a is zero. 50 51 Level: intermediate 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 X %D x %D, Y %D x %D",M1,N1,M2,N2); 73 if (m1 != m2 || n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Non conforming matrix add: local sizes X %D x %D, Y %D x %D",m1,n1,m2,n2); 74 if (!Y->assembled) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix (Y)"); 75 if (!X->assembled) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix (X)"); 76 if (a == (PetscScalar)0.0) PetscFunctionReturn(0); 77 if (Y == X) { 78 ierr = MatScale(Y,1.0 + a);CHKERRQ(ierr); 79 PetscFunctionReturn(0); 80 } 81 ierr = MatGetType(X,&t1);CHKERRQ(ierr); 82 ierr = MatGetType(Y,&t2);CHKERRQ(ierr); 83 ierr = PetscStrcmp(t1,t2,&sametype);CHKERRQ(ierr); 84 ierr = PetscLogEventBegin(MAT_AXPY,Y,0,0,0);CHKERRQ(ierr); 85 if (Y->ops->axpy && (sametype || X->ops->axpy == Y->ops->axpy)) { 86 ierr = (*Y->ops->axpy)(Y,a,X,str);CHKERRQ(ierr); 87 } else { 88 ierr = PetscStrcmp(t1,MATTRANSPOSEMAT,&transpose);CHKERRQ(ierr); 89 if (transpose) { 90 ierr = MatTransposeAXPY_Private(Y,a,X,str,X);CHKERRQ(ierr); 91 } else { 92 ierr = PetscStrcmp(t2,MATTRANSPOSEMAT,&transpose);CHKERRQ(ierr); 93 if (transpose) { 94 ierr = MatTransposeAXPY_Private(Y,a,X,str,Y);CHKERRQ(ierr); 95 } else { 96 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 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 = MatSetLayouts(preallocator,Y->rmap,Y->cmap);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 = MatSetLayouts(*B,Y->rmap,Y->cmap);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 PetscErrorCode ierr; 157 PetscBool isshell,isdense; 158 159 PetscFunctionBegin; 160 ierr = MatIsShell(Y,&isshell);CHKERRQ(ierr); 161 if (isshell) { /* MatShell has special support for AXPY */ 162 PetscErrorCode (*f)(Mat,PetscScalar,Mat,MatStructure); 163 164 ierr = MatGetOperation(Y,MATOP_AXPY,(void (**)(void))&f);CHKERRQ(ierr); 165 if (f) { 166 ierr = (*f)(Y,a,X,str);CHKERRQ(ierr); 167 PetscFunctionReturn(0); 168 } 169 } 170 /* no need to preallocate if Y is dense */ 171 ierr = PetscObjectBaseTypeCompareAny((PetscObject)Y,&isdense,MATSEQDENSE,MATMPIDENSE,"");CHKERRQ(ierr); 172 if (isdense && (str == DIFFERENT_NONZERO_PATTERN || str == UNKNOWN_NONZERO_PATTERN)) str = SUBSET_NONZERO_PATTERN; 173 if (str != DIFFERENT_NONZERO_PATTERN && str != UNKNOWN_NONZERO_PATTERN) { 174 PetscInt i,start,end,j,ncols,m,n; 175 const PetscInt *row; 176 PetscScalar *val; 177 const PetscScalar *vals; 178 179 ierr = MatGetSize(X,&m,&n);CHKERRQ(ierr); 180 ierr = MatGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 181 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 182 if (a == 1.0) { 183 for (i = start; i < end; i++) { 184 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 185 ierr = MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 186 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 187 } 188 } else { 189 PetscInt vs = 100; 190 /* realloc if needed, as this function may be used in parallel */ 191 ierr = PetscMalloc1(vs,&val);CHKERRQ(ierr); 192 for (i=start; i<end; i++) { 193 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 194 if (vs < ncols) { 195 vs = PetscMin(2*ncols,n); 196 ierr = PetscRealloc(vs*sizeof(*val),&val);CHKERRQ(ierr); 197 } 198 for (j=0; j<ncols; j++) val[j] = a*vals[j]; 199 ierr = MatSetValues(Y,1,&i,ncols,row,val,ADD_VALUES);CHKERRQ(ierr); 200 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 201 } 202 ierr = PetscFree(val);CHKERRQ(ierr); 203 } 204 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 205 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 206 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 207 } else { 208 Mat B; 209 210 ierr = MatAXPY_Basic_Preallocate(Y,X,&B);CHKERRQ(ierr); 211 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 212 /* TODO mat_tests-ex37_nsize-1_mat_type-baij_mat_block_size-2 fails with MatHeaderMerge */ 213 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 214 } 215 PetscFunctionReturn(0); 216 } 217 218 PetscErrorCode MatAXPY_BasicWithPreallocation(Mat B,Mat Y,PetscScalar a,Mat X,MatStructure str) 219 { 220 PetscInt i,start,end,j,ncols,m,n; 221 PetscErrorCode ierr; 222 const PetscInt *row; 223 PetscScalar *val; 224 const PetscScalar *vals; 225 226 PetscFunctionBegin; 227 ierr = MatGetSize(X,&m,&n);CHKERRQ(ierr); 228 ierr = MatGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 229 ierr = MatGetRowUpperTriangular(Y);CHKERRQ(ierr); 230 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 231 if (a == 1.0) { 232 for (i = start; i < end; i++) { 233 ierr = MatGetRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 234 ierr = MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 235 ierr = MatRestoreRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 236 237 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 238 ierr = MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 239 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 240 } 241 } else { 242 PetscInt vs = 100; 243 /* realloc if needed, as this function may be used in parallel */ 244 ierr = PetscMalloc1(vs,&val);CHKERRQ(ierr); 245 for (i=start; i<end; i++) { 246 ierr = MatGetRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 247 ierr = MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 248 ierr = MatRestoreRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 249 250 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 251 if (vs < ncols) { 252 vs = PetscMin(2*ncols,n); 253 ierr = PetscRealloc(vs*sizeof(*val),&val);CHKERRQ(ierr); 254 } 255 for (j=0; j<ncols; j++) val[j] = a*vals[j]; 256 ierr = MatSetValues(B,1,&i,ncols,row,val,ADD_VALUES);CHKERRQ(ierr); 257 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 258 } 259 ierr = PetscFree(val);CHKERRQ(ierr); 260 } 261 ierr = MatRestoreRowUpperTriangular(Y);CHKERRQ(ierr); 262 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 263 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 264 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 265 PetscFunctionReturn(0); 266 } 267 268 /*@ 269 MatShift - Computes Y = Y + a I, where a is a PetscScalar and I is the identity matrix. 270 271 Neighbor-wise Collective on Mat 272 273 Input Parameters: 274 + Y - the matrices 275 - a - the PetscScalar 276 277 Level: intermediate 278 279 Notes: 280 If the matrix Y is missing some diagonal entries this routine can be very slow. To make it fast one should initially 281 fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an 282 entry). No operation is performed when a is zero. 283 284 To form Y = Y + diag(V) use MatDiagonalSet() 285 286 .seealso: MatDiagonalSet(), MatScale(), MatDiagonalScale() 287 @*/ 288 PetscErrorCode MatShift(Mat Y,PetscScalar a) 289 { 290 PetscErrorCode ierr; 291 292 PetscFunctionBegin; 293 PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 294 if (!Y->assembled) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 295 if (Y->factortype) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 296 MatCheckPreallocated(Y,1); 297 if (a == 0.0) PetscFunctionReturn(0); 298 299 if (Y->ops->shift) { 300 ierr = (*Y->ops->shift)(Y,a);CHKERRQ(ierr); 301 } else { 302 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 303 } 304 305 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 306 PetscFunctionReturn(0); 307 } 308 309 PetscErrorCode MatDiagonalSet_Default(Mat Y,Vec D,InsertMode is) 310 { 311 PetscErrorCode ierr; 312 PetscInt i,start,end; 313 const PetscScalar *v; 314 315 PetscFunctionBegin; 316 ierr = MatGetOwnershipRange(Y,&start,&end);CHKERRQ(ierr); 317 ierr = VecGetArrayRead(D,&v);CHKERRQ(ierr); 318 for (i=start; i<end; i++) { 319 ierr = MatSetValues(Y,1,&i,1,&i,v+i-start,is);CHKERRQ(ierr); 320 } 321 ierr = VecRestoreArrayRead(D,&v);CHKERRQ(ierr); 322 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 323 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 324 PetscFunctionReturn(0); 325 } 326 327 /*@ 328 MatDiagonalSet - Computes Y = Y + D, where D is a diagonal matrix 329 that is represented as a vector. Or Y[i,i] = D[i] if InsertMode is 330 INSERT_VALUES. 331 332 Neighbor-wise Collective on Mat 333 334 Input Parameters: 335 + Y - the input matrix 336 . D - the diagonal matrix, represented as a vector 337 - i - INSERT_VALUES or ADD_VALUES 338 339 Notes: 340 If the matrix Y is missing some diagonal entries this routine can be very slow. To make it fast one should initially 341 fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an 342 entry). 343 344 Level: intermediate 345 346 .seealso: MatShift(), MatScale(), MatDiagonalScale() 347 @*/ 348 PetscErrorCode MatDiagonalSet(Mat Y,Vec D,InsertMode is) 349 { 350 PetscErrorCode ierr; 351 PetscInt matlocal,veclocal; 352 353 PetscFunctionBegin; 354 PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 355 PetscValidHeaderSpecific(D,VEC_CLASSID,2); 356 ierr = MatGetLocalSize(Y,&matlocal,NULL);CHKERRQ(ierr); 357 ierr = VecGetLocalSize(D,&veclocal);CHKERRQ(ierr); 358 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); 359 if (Y->ops->diagonalset) { 360 ierr = (*Y->ops->diagonalset)(Y,D,is);CHKERRQ(ierr); 361 } else { 362 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 363 } 364 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 365 PetscFunctionReturn(0); 366 } 367 368 /*@ 369 MatAYPX - Computes Y = a*Y + X. 370 371 Logically on Mat 372 373 Input Parameters: 374 + a - the PetscScalar multiplier 375 . Y - the first matrix 376 . X - the second matrix 377 - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN or SUBSET_NONZERO_PATTERN 378 379 Level: intermediate 380 381 .seealso: MatAXPY() 382 @*/ 383 PetscErrorCode MatAYPX(Mat Y,PetscScalar a,Mat X,MatStructure str) 384 { 385 PetscErrorCode ierr; 386 387 PetscFunctionBegin; 388 ierr = MatScale(Y,a);CHKERRQ(ierr); 389 ierr = MatAXPY(Y,1.0,X,str);CHKERRQ(ierr); 390 PetscFunctionReturn(0); 391 } 392 393 /*@ 394 MatComputeOperator - Computes the explicit matrix 395 396 Collective on Mat 397 398 Input Parameter: 399 + inmat - the matrix 400 - mattype - the matrix type for the explicit operator 401 402 Output Parameter: 403 . mat - the explict operator 404 405 Notes: 406 This computation is done by applying the operators to columns of the identity matrix. 407 This routine is costly in general, and is recommended for use only with relatively small systems. 408 Currently, this routine uses a dense matrix format if mattype == NULL. 409 410 Level: advanced 411 412 @*/ 413 PetscErrorCode MatComputeOperator(Mat inmat,MatType mattype,Mat *mat) 414 { 415 PetscErrorCode ierr; 416 417 PetscFunctionBegin; 418 PetscValidHeaderSpecific(inmat,MAT_CLASSID,1); 419 PetscValidPointer(mat,3); 420 ierr = MatConvert_Shell(inmat,mattype ? mattype : MATDENSE,MAT_INITIAL_MATRIX,mat);CHKERRQ(ierr); 421 PetscFunctionReturn(0); 422 } 423 424 /*@ 425 MatComputeOperatorTranspose - Computes the explicit matrix representation of 426 a give matrix that can apply MatMultTranspose() 427 428 Collective on Mat 429 430 Input Parameter: 431 . inmat - the matrix 432 433 Output Parameter: 434 . mat - the explict operator transposed 435 436 Notes: 437 This computation is done by applying the transpose of the operator to columns of the identity matrix. 438 This routine is costly in general, and is recommended for use only with relatively small systems. 439 Currently, this routine uses a dense matrix format if mattype == NULL. 440 441 Level: advanced 442 443 @*/ 444 PetscErrorCode MatComputeOperatorTranspose(Mat inmat,MatType mattype,Mat *mat) 445 { 446 Mat A; 447 PetscErrorCode ierr; 448 449 PetscFunctionBegin; 450 PetscValidHeaderSpecific(inmat,MAT_CLASSID,1); 451 PetscValidPointer(mat,3); 452 ierr = MatCreateTranspose(inmat,&A);CHKERRQ(ierr); 453 ierr = MatConvert_Shell(A,mattype ? mattype : MATDENSE,MAT_INITIAL_MATRIX,mat);CHKERRQ(ierr); 454 ierr = MatDestroy(&A);CHKERRQ(ierr); 455 PetscFunctionReturn(0); 456 } 457 458 /*@ 459 MatChop - Set all values in the matrix less than the tolerance to zero 460 461 Input Parameters: 462 + A - The matrix 463 - tol - The zero tolerance 464 465 Output Parameters: 466 . A - The chopped matrix 467 468 Level: intermediate 469 470 .seealso: MatCreate(), MatZeroEntries() 471 @*/ 472 PetscErrorCode MatChop(Mat A, PetscReal tol) 473 { 474 PetscScalar *newVals; 475 PetscInt *newCols; 476 PetscInt rStart, rEnd, numRows, maxRows, r, colMax = 0; 477 PetscErrorCode ierr; 478 479 PetscFunctionBegin; 480 ierr = MatGetOwnershipRange(A, &rStart, &rEnd);CHKERRQ(ierr); 481 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 482 for (r = rStart; r < rEnd; ++r) { 483 PetscInt ncols; 484 485 ierr = MatGetRow(A, r, &ncols, NULL, NULL);CHKERRQ(ierr); 486 colMax = PetscMax(colMax, ncols);CHKERRQ(ierr); 487 ierr = MatRestoreRow(A, r, &ncols, NULL, NULL);CHKERRQ(ierr); 488 } 489 numRows = rEnd - rStart; 490 ierr = MPIU_Allreduce(&numRows, &maxRows, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 491 ierr = PetscMalloc2(colMax,&newCols,colMax,&newVals);CHKERRQ(ierr); 492 for (r = rStart; r < rStart+maxRows; ++r) { 493 const PetscScalar *vals; 494 const PetscInt *cols; 495 PetscInt ncols, newcols, c; 496 497 if (r < rEnd) { 498 ierr = MatGetRow(A, r, &ncols, &cols, &vals);CHKERRQ(ierr); 499 for (c = 0; c < ncols; ++c) { 500 newCols[c] = cols[c]; 501 newVals[c] = PetscAbsScalar(vals[c]) < tol ? 0.0 : vals[c]; 502 } 503 newcols = ncols; 504 ierr = MatRestoreRow(A, r, &ncols, &cols, &vals);CHKERRQ(ierr); 505 ierr = MatSetValues(A, 1, &r, newcols, newCols, newVals, INSERT_VALUES);CHKERRQ(ierr); 506 } 507 ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 508 ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 509 } 510 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 511 ierr = PetscFree2(newCols,newVals);CHKERRQ(ierr); 512 PetscFunctionReturn(0); 513 } 514