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 || 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 = 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,isdense; 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 /* 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 = SUBSET_NONZERO_PATTERN; 173 if (str != DIFFERENT_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 PetscScalar *v; 314 315 PetscFunctionBegin; 316 ierr = MatGetOwnershipRange(Y,&start,&end);CHKERRQ(ierr); 317 ierr = VecGetArray(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 = VecRestoreArray(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 PetscScalar one = 1.0; 386 PetscErrorCode ierr; 387 PetscInt mX,mY,nX,nY; 388 389 PetscFunctionBegin; 390 PetscValidHeaderSpecific(X,MAT_CLASSID,3); 391 PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 392 PetscValidLogicalCollectiveScalar(Y,a,2); 393 ierr = MatGetSize(X,&mX,&nX);CHKERRQ(ierr); 394 ierr = MatGetSize(X,&mY,&nY);CHKERRQ(ierr); 395 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); 396 ierr = MatScale(Y,a);CHKERRQ(ierr); 397 ierr = MatAXPY(Y,one,X,str);CHKERRQ(ierr); 398 PetscFunctionReturn(0); 399 } 400 401 /*@ 402 MatComputeOperator - Computes the explicit matrix 403 404 Collective on Mat 405 406 Input Parameter: 407 + inmat - the matrix 408 - mattype - the matrix type for the explicit operator 409 410 Output Parameter: 411 . mat - the explict operator 412 413 Notes: 414 This computation is done by applying the operators to columns of the identity matrix. 415 This routine is costly in general, and is recommended for use only with relatively small systems. 416 Currently, this routine uses a dense matrix format if mattype == NULL. 417 418 Level: advanced 419 420 @*/ 421 PetscErrorCode MatComputeOperator(Mat inmat,MatType mattype,Mat *mat) 422 { 423 PetscErrorCode ierr; 424 425 PetscFunctionBegin; 426 PetscValidHeaderSpecific(inmat,MAT_CLASSID,1); 427 PetscValidPointer(mat,3); 428 ierr = MatConvert_Shell(inmat,mattype ? mattype : MATDENSE,MAT_INITIAL_MATRIX,mat);CHKERRQ(ierr); 429 PetscFunctionReturn(0); 430 } 431 432 /*@ 433 MatComputeOperatorTranspose - Computes the explicit matrix representation of 434 a give matrix that can apply MatMultTranspose() 435 436 Collective on Mat 437 438 Input Parameter: 439 . inmat - the matrix 440 441 Output Parameter: 442 . mat - the explict operator transposed 443 444 Notes: 445 This computation is done by applying the transpose of the operator to columns of the identity matrix. 446 This routine is costly in general, and is recommended for use only with relatively small systems. 447 Currently, this routine uses a dense matrix format if mattype == NULL. 448 449 Level: advanced 450 451 @*/ 452 PetscErrorCode MatComputeOperatorTranspose(Mat inmat,MatType mattype,Mat *mat) 453 { 454 Mat A; 455 PetscErrorCode ierr; 456 457 PetscFunctionBegin; 458 PetscValidHeaderSpecific(inmat,MAT_CLASSID,1); 459 PetscValidPointer(mat,3); 460 ierr = MatCreateTranspose(inmat,&A);CHKERRQ(ierr); 461 ierr = MatConvert_Shell(A,mattype ? mattype : MATDENSE,MAT_INITIAL_MATRIX,mat);CHKERRQ(ierr); 462 ierr = MatDestroy(&A);CHKERRQ(ierr); 463 PetscFunctionReturn(0); 464 } 465 466 /*@ 467 MatChop - Set all values in the matrix less than the tolerance to zero 468 469 Input Parameters: 470 + A - The matrix 471 - tol - The zero tolerance 472 473 Output Parameters: 474 . A - The chopped matrix 475 476 Level: intermediate 477 478 .seealso: MatCreate(), MatZeroEntries() 479 @*/ 480 PetscErrorCode MatChop(Mat A, PetscReal tol) 481 { 482 PetscScalar *newVals; 483 PetscInt *newCols; 484 PetscInt rStart, rEnd, numRows, maxRows, r, colMax = 0; 485 PetscErrorCode ierr; 486 487 PetscFunctionBegin; 488 ierr = MatGetOwnershipRange(A, &rStart, &rEnd);CHKERRQ(ierr); 489 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 490 for (r = rStart; r < rEnd; ++r) { 491 PetscInt ncols; 492 493 ierr = MatGetRow(A, r, &ncols, NULL, NULL);CHKERRQ(ierr); 494 colMax = PetscMax(colMax, ncols);CHKERRQ(ierr); 495 ierr = MatRestoreRow(A, r, &ncols, NULL, NULL);CHKERRQ(ierr); 496 } 497 numRows = rEnd - rStart; 498 ierr = MPIU_Allreduce(&numRows, &maxRows, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 499 ierr = PetscMalloc2(colMax,&newCols,colMax,&newVals);CHKERRQ(ierr); 500 for (r = rStart; r < rStart+maxRows; ++r) { 501 const PetscScalar *vals; 502 const PetscInt *cols; 503 PetscInt ncols, newcols, c; 504 505 if (r < rEnd) { 506 ierr = MatGetRow(A, r, &ncols, &cols, &vals);CHKERRQ(ierr); 507 for (c = 0; c < ncols; ++c) { 508 newCols[c] = cols[c]; 509 newVals[c] = PetscAbsScalar(vals[c]) < tol ? 0.0 : vals[c]; 510 } 511 newcols = ncols; 512 ierr = MatRestoreRow(A, r, &ncols, &cols, &vals);CHKERRQ(ierr); 513 ierr = MatSetValues(A, 1, &r, newcols, newCols, newVals, INSERT_VALUES);CHKERRQ(ierr); 514 } 515 ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 516 ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 517 } 518 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 519 ierr = PetscFree2(newCols,newVals);CHKERRQ(ierr); 520 PetscFunctionReturn(0); 521 } 522