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, UNKNOWN_NONZERO_PATTERN, or SUBSET_NONZERO_PATTERN (nonzeros of X is a subset of Y's) 47 48 Notes: No operation is performed when a is zero. 49 50 Level: intermediate 51 52 .seealso: MatAYPX() 53 @*/ 54 PetscErrorCode MatAXPY(Mat Y,PetscScalar a,Mat X,MatStructure str) 55 { 56 PetscErrorCode ierr; 57 PetscInt M1,M2,N1,N2; 58 PetscInt m1,m2,n1,n2; 59 MatType t1,t2; 60 PetscBool sametype,transpose; 61 62 PetscFunctionBegin; 63 PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 64 PetscValidLogicalCollectiveScalar(Y,a,2); 65 PetscValidHeaderSpecific(X,MAT_CLASSID,3); 66 PetscCheckSameComm(Y,1,X,3); 67 ierr = MatGetSize(X,&M1,&N1);CHKERRQ(ierr); 68 ierr = MatGetSize(Y,&M2,&N2);CHKERRQ(ierr); 69 ierr = MatGetLocalSize(X,&m1,&n1);CHKERRQ(ierr); 70 ierr = MatGetLocalSize(Y,&m2,&n2);CHKERRQ(ierr); 71 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); 72 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); 73 if (!Y->assembled) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix (Y)"); 74 if (!X->assembled) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix (X)"); 75 if (a == (PetscScalar)0.0) PetscFunctionReturn(0); 76 if (Y == X) { 77 ierr = MatScale(Y,1.0 + a);CHKERRQ(ierr); 78 PetscFunctionReturn(0); 79 } 80 ierr = MatGetType(X,&t1);CHKERRQ(ierr); 81 ierr = MatGetType(Y,&t2);CHKERRQ(ierr); 82 ierr = PetscStrcmp(t1,t2,&sametype);CHKERRQ(ierr); 83 ierr = PetscLogEventBegin(MAT_AXPY,Y,0,0,0);CHKERRQ(ierr); 84 if (Y->ops->axpy && (sametype || X->ops->axpy == Y->ops->axpy)) { 85 ierr = (*Y->ops->axpy)(Y,a,X,str);CHKERRQ(ierr); 86 } else { 87 ierr = PetscStrcmp(t1,MATTRANSPOSEMAT,&transpose);CHKERRQ(ierr); 88 if (transpose) { 89 ierr = MatTransposeAXPY_Private(Y,a,X,str,X);CHKERRQ(ierr); 90 } else { 91 ierr = PetscStrcmp(t2,MATTRANSPOSEMAT,&transpose);CHKERRQ(ierr); 92 if (transpose) { 93 ierr = MatTransposeAXPY_Private(Y,a,X,str,Y);CHKERRQ(ierr); 94 } else { 95 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 96 } 97 } 98 } 99 ierr = PetscLogEventEnd(MAT_AXPY,Y,0,0,0);CHKERRQ(ierr); 100 PetscFunctionReturn(0); 101 } 102 103 PetscErrorCode MatAXPY_Basic_Preallocate(Mat Y, Mat X, Mat *B) 104 { 105 PetscErrorCode ierr; 106 PetscErrorCode (*preall)(Mat,Mat,Mat*) = NULL; 107 108 PetscFunctionBegin; 109 /* look for any available faster alternative to the general preallocator */ 110 ierr = PetscObjectQueryFunction((PetscObject)Y,"MatAXPYGetPreallocation_C",&preall);CHKERRQ(ierr); 111 if (preall) { 112 ierr = (*preall)(Y,X,B);CHKERRQ(ierr); 113 } else { /* Use MatPrellocator, assumes same row-col distribution */ 114 Mat preallocator; 115 PetscInt r,rstart,rend; 116 PetscInt m,n,M,N; 117 118 ierr = MatGetRowUpperTriangular(Y);CHKERRQ(ierr); 119 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 120 ierr = MatGetSize(Y,&M,&N);CHKERRQ(ierr); 121 ierr = MatGetLocalSize(Y,&m,&n);CHKERRQ(ierr); 122 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&preallocator);CHKERRQ(ierr); 123 ierr = MatSetType(preallocator,MATPREALLOCATOR);CHKERRQ(ierr); 124 ierr = MatSetLayouts(preallocator,Y->rmap,Y->cmap);CHKERRQ(ierr); 125 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 126 ierr = MatGetOwnershipRange(preallocator,&rstart,&rend);CHKERRQ(ierr); 127 for (r = rstart; r < rend; ++r) { 128 PetscInt ncols; 129 const PetscInt *row; 130 const PetscScalar *vals; 131 132 ierr = MatGetRow(Y,r,&ncols,&row,&vals);CHKERRQ(ierr); 133 ierr = MatSetValues(preallocator,1,&r,ncols,row,vals,INSERT_VALUES);CHKERRQ(ierr); 134 ierr = MatRestoreRow(Y,r,&ncols,&row,&vals);CHKERRQ(ierr); 135 ierr = MatGetRow(X,r,&ncols,&row,&vals);CHKERRQ(ierr); 136 ierr = MatSetValues(preallocator,1,&r,ncols,row,vals,INSERT_VALUES);CHKERRQ(ierr); 137 ierr = MatRestoreRow(X,r,&ncols,&row,&vals);CHKERRQ(ierr); 138 } 139 ierr = MatSetOption(preallocator,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 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,isnest; 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) { 173 ierr = PetscObjectTypeCompare((PetscObject)X,MATNEST,&isnest);CHKERRQ(ierr); 174 if (isnest) { 175 ierr = MatAXPY_Dense_Nest(Y,a,X);CHKERRQ(ierr); 176 PetscFunctionReturn(0); 177 } 178 if (str == DIFFERENT_NONZERO_PATTERN || str == UNKNOWN_NONZERO_PATTERN) str = SUBSET_NONZERO_PATTERN; 179 } 180 if (str != DIFFERENT_NONZERO_PATTERN && str != UNKNOWN_NONZERO_PATTERN) { 181 PetscInt i,start,end,j,ncols,m,n; 182 const PetscInt *row; 183 PetscScalar *val; 184 const PetscScalar *vals; 185 186 ierr = MatGetSize(X,&m,&n);CHKERRQ(ierr); 187 ierr = MatGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 188 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 189 if (a == 1.0) { 190 for (i = start; i < end; i++) { 191 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 192 ierr = MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 193 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 194 } 195 } else { 196 PetscInt vs = 100; 197 /* realloc if needed, as this function may be used in parallel */ 198 ierr = PetscMalloc1(vs,&val);CHKERRQ(ierr); 199 for (i=start; i<end; i++) { 200 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 201 if (vs < ncols) { 202 vs = PetscMin(2*ncols,n); 203 ierr = PetscRealloc(vs*sizeof(*val),&val);CHKERRQ(ierr); 204 } 205 for (j=0; j<ncols; j++) val[j] = a*vals[j]; 206 ierr = MatSetValues(Y,1,&i,ncols,row,val,ADD_VALUES);CHKERRQ(ierr); 207 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 208 } 209 ierr = PetscFree(val);CHKERRQ(ierr); 210 } 211 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 212 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 213 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 214 } else { 215 Mat B; 216 217 ierr = MatAXPY_Basic_Preallocate(Y,X,&B);CHKERRQ(ierr); 218 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 219 ierr = MatHeaderMerge(Y,&B);CHKERRQ(ierr); 220 } 221 PetscFunctionReturn(0); 222 } 223 224 PetscErrorCode MatAXPY_BasicWithPreallocation(Mat B,Mat Y,PetscScalar a,Mat X,MatStructure str) 225 { 226 PetscInt i,start,end,j,ncols,m,n; 227 PetscErrorCode ierr; 228 const PetscInt *row; 229 PetscScalar *val; 230 const PetscScalar *vals; 231 232 PetscFunctionBegin; 233 ierr = MatGetSize(X,&m,&n);CHKERRQ(ierr); 234 ierr = MatGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 235 ierr = MatGetRowUpperTriangular(Y);CHKERRQ(ierr); 236 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 237 if (a == 1.0) { 238 for (i = start; i < end; i++) { 239 ierr = MatGetRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 240 ierr = MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 241 ierr = MatRestoreRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 242 243 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 244 ierr = MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 245 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 246 } 247 } else { 248 PetscInt vs = 100; 249 /* realloc if needed, as this function may be used in parallel */ 250 ierr = PetscMalloc1(vs,&val);CHKERRQ(ierr); 251 for (i=start; i<end; i++) { 252 ierr = MatGetRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 253 ierr = MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 254 ierr = MatRestoreRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 255 256 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 257 if (vs < ncols) { 258 vs = PetscMin(2*ncols,n); 259 ierr = PetscRealloc(vs*sizeof(*val),&val);CHKERRQ(ierr); 260 } 261 for (j=0; j<ncols; j++) val[j] = a*vals[j]; 262 ierr = MatSetValues(B,1,&i,ncols,row,val,ADD_VALUES);CHKERRQ(ierr); 263 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 264 } 265 ierr = PetscFree(val);CHKERRQ(ierr); 266 } 267 ierr = MatRestoreRowUpperTriangular(Y);CHKERRQ(ierr); 268 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 269 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 270 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 271 PetscFunctionReturn(0); 272 } 273 274 /*@ 275 MatShift - Computes Y = Y + a I, where a is a PetscScalar and I is the identity matrix. 276 277 Neighbor-wise Collective on Mat 278 279 Input Parameters: 280 + Y - the matrices 281 - a - the PetscScalar 282 283 Level: intermediate 284 285 Notes: 286 If the matrix Y is missing some diagonal entries this routine can be very slow. To make it fast one should initially 287 fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an 288 entry). No operation is performed when a is zero. 289 290 To form Y = Y + diag(V) use MatDiagonalSet() 291 292 .seealso: MatDiagonalSet(), MatScale(), MatDiagonalScale() 293 @*/ 294 PetscErrorCode MatShift(Mat Y,PetscScalar a) 295 { 296 PetscErrorCode ierr; 297 298 PetscFunctionBegin; 299 PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 300 if (!Y->assembled) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 301 if (Y->factortype) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 302 MatCheckPreallocated(Y,1); 303 if (a == 0.0) PetscFunctionReturn(0); 304 305 if (Y->ops->shift) { 306 ierr = (*Y->ops->shift)(Y,a);CHKERRQ(ierr); 307 } else { 308 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 309 } 310 311 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 312 PetscFunctionReturn(0); 313 } 314 315 PetscErrorCode MatDiagonalSet_Default(Mat Y,Vec D,InsertMode is) 316 { 317 PetscErrorCode ierr; 318 PetscInt i,start,end; 319 const PetscScalar *v; 320 321 PetscFunctionBegin; 322 ierr = MatGetOwnershipRange(Y,&start,&end);CHKERRQ(ierr); 323 ierr = VecGetArrayRead(D,&v);CHKERRQ(ierr); 324 for (i=start; i<end; i++) { 325 ierr = MatSetValues(Y,1,&i,1,&i,v+i-start,is);CHKERRQ(ierr); 326 } 327 ierr = VecRestoreArrayRead(D,&v);CHKERRQ(ierr); 328 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 329 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 330 PetscFunctionReturn(0); 331 } 332 333 /*@ 334 MatDiagonalSet - Computes Y = Y + D, where D is a diagonal matrix 335 that is represented as a vector. Or Y[i,i] = D[i] if InsertMode is 336 INSERT_VALUES. 337 338 Neighbor-wise Collective on Mat 339 340 Input Parameters: 341 + Y - the input matrix 342 . D - the diagonal matrix, represented as a vector 343 - i - INSERT_VALUES or ADD_VALUES 344 345 Notes: 346 If the matrix Y is missing some diagonal entries this routine can be very slow. To make it fast one should initially 347 fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an 348 entry). 349 350 Level: intermediate 351 352 .seealso: MatShift(), MatScale(), MatDiagonalScale() 353 @*/ 354 PetscErrorCode MatDiagonalSet(Mat Y,Vec D,InsertMode is) 355 { 356 PetscErrorCode ierr; 357 PetscInt matlocal,veclocal; 358 359 PetscFunctionBegin; 360 PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 361 PetscValidHeaderSpecific(D,VEC_CLASSID,2); 362 ierr = MatGetLocalSize(Y,&matlocal,NULL);CHKERRQ(ierr); 363 ierr = VecGetLocalSize(D,&veclocal);CHKERRQ(ierr); 364 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); 365 if (Y->ops->diagonalset) { 366 ierr = (*Y->ops->diagonalset)(Y,D,is);CHKERRQ(ierr); 367 } else { 368 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 369 } 370 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 371 PetscFunctionReturn(0); 372 } 373 374 /*@ 375 MatAYPX - Computes Y = a*Y + X. 376 377 Logically on Mat 378 379 Input Parameters: 380 + a - the PetscScalar multiplier 381 . Y - the first matrix 382 . X - the second matrix 383 - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN, UNKNOWN_NONZERO_PATTERN, or SUBSET_NONZERO_PATTERN (nonzeros of X is a subset of Y's) 384 385 Level: intermediate 386 387 .seealso: MatAXPY() 388 @*/ 389 PetscErrorCode MatAYPX(Mat Y,PetscScalar a,Mat X,MatStructure str) 390 { 391 PetscErrorCode ierr; 392 393 PetscFunctionBegin; 394 ierr = MatScale(Y,a);CHKERRQ(ierr); 395 ierr = MatAXPY(Y,1.0,X,str);CHKERRQ(ierr); 396 PetscFunctionReturn(0); 397 } 398 399 /*@ 400 MatComputeOperator - Computes the explicit matrix 401 402 Collective on Mat 403 404 Input Parameters: 405 + inmat - the matrix 406 - mattype - the matrix type for the explicit operator 407 408 Output Parameter: 409 . mat - the explicit operator 410 411 Notes: 412 This computation is done by applying the operators to columns of the identity matrix. 413 This routine is costly in general, and is recommended for use only with relatively small systems. 414 Currently, this routine uses a dense matrix format if mattype == NULL. 415 416 Level: advanced 417 418 @*/ 419 PetscErrorCode MatComputeOperator(Mat inmat,MatType mattype,Mat *mat) 420 { 421 PetscErrorCode ierr; 422 423 PetscFunctionBegin; 424 PetscValidHeaderSpecific(inmat,MAT_CLASSID,1); 425 PetscValidPointer(mat,3); 426 ierr = MatConvert_Shell(inmat,mattype ? mattype : MATDENSE,MAT_INITIAL_MATRIX,mat);CHKERRQ(ierr); 427 PetscFunctionReturn(0); 428 } 429 430 /*@ 431 MatComputeOperatorTranspose - Computes the explicit matrix representation of 432 a give matrix that can apply MatMultTranspose() 433 434 Collective on Mat 435 436 Input Parameters: 437 + inmat - the matrix 438 - mattype - the matrix type for the explicit operator 439 440 Output Parameter: 441 . mat - the explicit operator transposed 442 443 Notes: 444 This computation is done by applying the transpose of the operator to columns of the identity matrix. 445 This routine is costly in general, and is recommended for use only with relatively small systems. 446 Currently, this routine uses a dense matrix format if mattype == NULL. 447 448 Level: advanced 449 450 @*/ 451 PetscErrorCode MatComputeOperatorTranspose(Mat inmat,MatType mattype,Mat *mat) 452 { 453 Mat A; 454 PetscErrorCode ierr; 455 456 PetscFunctionBegin; 457 PetscValidHeaderSpecific(inmat,MAT_CLASSID,1); 458 PetscValidPointer(mat,3); 459 ierr = MatCreateTranspose(inmat,&A);CHKERRQ(ierr); 460 ierr = MatConvert_Shell(A,mattype ? mattype : MATDENSE,MAT_INITIAL_MATRIX,mat);CHKERRQ(ierr); 461 ierr = MatDestroy(&A);CHKERRQ(ierr); 462 PetscFunctionReturn(0); 463 } 464 465 /*@ 466 MatChop - Set all values in the matrix less than the tolerance to zero 467 468 Input Parameters: 469 + A - The matrix 470 - tol - The zero tolerance 471 472 Output Parameters: 473 . A - The chopped matrix 474 475 Level: intermediate 476 477 .seealso: MatCreate(), MatZeroEntries() 478 @*/ 479 PetscErrorCode MatChop(Mat A, PetscReal tol) 480 { 481 Mat a; 482 PetscScalar *newVals; 483 PetscInt *newCols, rStart, rEnd, numRows, maxRows, r, colMax = 0; 484 PetscBool flg; 485 PetscErrorCode ierr; 486 487 PetscFunctionBegin; 488 ierr = PetscObjectBaseTypeCompareAny((PetscObject)A, &flg, MATSEQDENSE, MATMPIDENSE, "");CHKERRQ(ierr); 489 if (flg) { 490 ierr = MatDenseGetLocalMatrix(A, &a);CHKERRQ(ierr); 491 ierr = MatDenseGetLDA(a, &r);CHKERRQ(ierr); 492 ierr = MatGetSize(a, &rStart, &rEnd);CHKERRQ(ierr); 493 ierr = MatDenseGetArray(a, &newVals);CHKERRQ(ierr); 494 for (; colMax < rEnd; ++colMax) { 495 for (maxRows = 0; maxRows < rStart; ++maxRows) { 496 newVals[maxRows + colMax * r] = PetscAbsScalar(newVals[maxRows + colMax * r]) < tol ? 0.0 : newVals[maxRows + colMax * r]; 497 } 498 } 499 ierr = MatDenseRestoreArray(a, &newVals);CHKERRQ(ierr); 500 } else { 501 ierr = MatGetOwnershipRange(A, &rStart, &rEnd);CHKERRQ(ierr); 502 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 503 for (r = rStart; r < rEnd; ++r) { 504 PetscInt ncols; 505 506 ierr = MatGetRow(A, r, &ncols, NULL, NULL);CHKERRQ(ierr); 507 colMax = PetscMax(colMax, ncols);CHKERRQ(ierr); 508 ierr = MatRestoreRow(A, r, &ncols, NULL, NULL);CHKERRQ(ierr); 509 } 510 numRows = rEnd - rStart; 511 ierr = MPIU_Allreduce(&numRows, &maxRows, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)A));CHKERRMPI(ierr); 512 ierr = PetscMalloc2(colMax, &newCols, colMax, &newVals);CHKERRQ(ierr); 513 ierr = MatGetOption(A, MAT_NO_OFF_PROC_ENTRIES, &flg);CHKERRQ(ierr); /* cache user-defined value */ 514 ierr = MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);CHKERRQ(ierr); 515 /* short-circuit code in MatAssemblyBegin() and MatAssemblyEnd() */ 516 /* that are potentially called many times depending on the distribution of A */ 517 for (r = rStart; r < rStart+maxRows; ++r) { 518 const PetscScalar *vals; 519 const PetscInt *cols; 520 PetscInt ncols, newcols, c; 521 522 if (r < rEnd) { 523 ierr = MatGetRow(A, r, &ncols, &cols, &vals);CHKERRQ(ierr); 524 for (c = 0; c < ncols; ++c) { 525 newCols[c] = cols[c]; 526 newVals[c] = PetscAbsScalar(vals[c]) < tol ? 0.0 : vals[c]; 527 } 528 newcols = ncols; 529 ierr = MatRestoreRow(A, r, &ncols, &cols, &vals);CHKERRQ(ierr); 530 ierr = MatSetValues(A, 1, &r, newcols, newCols, newVals, INSERT_VALUES);CHKERRQ(ierr); 531 } 532 ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 533 ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 534 } 535 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 536 ierr = PetscFree2(newCols, newVals);CHKERRQ(ierr); 537 ierr = MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, flg);CHKERRQ(ierr); /* reset option to its user-defined value */ 538 } 539 PetscFunctionReturn(0); 540 } 541