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