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 = MatSetOption(preallocator,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 141 ierr = MatAssemblyBegin(preallocator,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 142 ierr = MatAssemblyEnd(preallocator,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 143 ierr = MatRestoreRowUpperTriangular(Y);CHKERRQ(ierr); 144 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 145 146 ierr = MatCreate(PetscObjectComm((PetscObject)Y),B);CHKERRQ(ierr); 147 ierr = MatSetType(*B,((PetscObject)Y)->type_name);CHKERRQ(ierr); 148 ierr = MatSetLayouts(*B,Y->rmap,Y->cmap);CHKERRQ(ierr); 149 ierr = MatPreallocatorPreallocate(preallocator,PETSC_FALSE,*B);CHKERRQ(ierr); 150 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 151 } 152 PetscFunctionReturn(0); 153 } 154 155 PetscErrorCode MatAXPY_Basic(Mat Y,PetscScalar a,Mat X,MatStructure str) 156 { 157 PetscErrorCode ierr; 158 PetscBool isshell,isdense,isnest; 159 160 PetscFunctionBegin; 161 ierr = MatIsShell(Y,&isshell);CHKERRQ(ierr); 162 if (isshell) { /* MatShell has special support for AXPY */ 163 PetscErrorCode (*f)(Mat,PetscScalar,Mat,MatStructure); 164 165 ierr = MatGetOperation(Y,MATOP_AXPY,(void (**)(void))&f);CHKERRQ(ierr); 166 if (f) { 167 ierr = (*f)(Y,a,X,str);CHKERRQ(ierr); 168 PetscFunctionReturn(0); 169 } 170 } 171 /* no need to preallocate if Y is dense */ 172 ierr = PetscObjectBaseTypeCompareAny((PetscObject)Y,&isdense,MATSEQDENSE,MATMPIDENSE,"");CHKERRQ(ierr); 173 if (isdense) { 174 ierr = PetscObjectTypeCompare((PetscObject)X,MATNEST,&isnest);CHKERRQ(ierr); 175 if (isnest) { 176 ierr = MatAXPY_Dense_Nest(Y,a,X);CHKERRQ(ierr); 177 PetscFunctionReturn(0); 178 } 179 if (str == DIFFERENT_NONZERO_PATTERN || str == UNKNOWN_NONZERO_PATTERN) str = SUBSET_NONZERO_PATTERN; 180 } 181 if (str != DIFFERENT_NONZERO_PATTERN && str != UNKNOWN_NONZERO_PATTERN) { 182 PetscInt i,start,end,j,ncols,m,n; 183 const PetscInt *row; 184 PetscScalar *val; 185 const PetscScalar *vals; 186 187 ierr = MatGetSize(X,&m,&n);CHKERRQ(ierr); 188 ierr = MatGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 189 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 190 if (a == 1.0) { 191 for (i = start; i < end; i++) { 192 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 193 ierr = MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 194 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 195 } 196 } else { 197 PetscInt vs = 100; 198 /* realloc if needed, as this function may be used in parallel */ 199 ierr = PetscMalloc1(vs,&val);CHKERRQ(ierr); 200 for (i=start; i<end; i++) { 201 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 202 if (vs < ncols) { 203 vs = PetscMin(2*ncols,n); 204 ierr = PetscRealloc(vs*sizeof(*val),&val);CHKERRQ(ierr); 205 } 206 for (j=0; j<ncols; j++) val[j] = a*vals[j]; 207 ierr = MatSetValues(Y,1,&i,ncols,row,val,ADD_VALUES);CHKERRQ(ierr); 208 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 209 } 210 ierr = PetscFree(val);CHKERRQ(ierr); 211 } 212 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 213 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 214 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 215 } else { 216 Mat B; 217 218 ierr = MatAXPY_Basic_Preallocate(Y,X,&B);CHKERRQ(ierr); 219 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 220 /* TODO mat_tests-ex37_nsize-1_mat_type-baij_mat_block_size-2 fails with MatHeaderMerge */ 221 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 222 } 223 PetscFunctionReturn(0); 224 } 225 226 PetscErrorCode MatAXPY_BasicWithPreallocation(Mat B,Mat Y,PetscScalar a,Mat X,MatStructure str) 227 { 228 PetscInt i,start,end,j,ncols,m,n; 229 PetscErrorCode ierr; 230 const PetscInt *row; 231 PetscScalar *val; 232 const PetscScalar *vals; 233 234 PetscFunctionBegin; 235 ierr = MatGetSize(X,&m,&n);CHKERRQ(ierr); 236 ierr = MatGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 237 ierr = MatGetRowUpperTriangular(Y);CHKERRQ(ierr); 238 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 239 if (a == 1.0) { 240 for (i = start; i < end; i++) { 241 ierr = MatGetRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 242 ierr = MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 243 ierr = MatRestoreRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 244 245 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 246 ierr = MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 247 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 248 } 249 } else { 250 PetscInt vs = 100; 251 /* realloc if needed, as this function may be used in parallel */ 252 ierr = PetscMalloc1(vs,&val);CHKERRQ(ierr); 253 for (i=start; i<end; i++) { 254 ierr = MatGetRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 255 ierr = MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 256 ierr = MatRestoreRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 257 258 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 259 if (vs < ncols) { 260 vs = PetscMin(2*ncols,n); 261 ierr = PetscRealloc(vs*sizeof(*val),&val);CHKERRQ(ierr); 262 } 263 for (j=0; j<ncols; j++) val[j] = a*vals[j]; 264 ierr = MatSetValues(B,1,&i,ncols,row,val,ADD_VALUES);CHKERRQ(ierr); 265 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 266 } 267 ierr = PetscFree(val);CHKERRQ(ierr); 268 } 269 ierr = MatRestoreRowUpperTriangular(Y);CHKERRQ(ierr); 270 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 271 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 272 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 273 PetscFunctionReturn(0); 274 } 275 276 /*@ 277 MatShift - Computes Y = Y + a I, where a is a PetscScalar and I is the identity matrix. 278 279 Neighbor-wise Collective on Mat 280 281 Input Parameters: 282 + Y - the matrices 283 - a - the PetscScalar 284 285 Level: intermediate 286 287 Notes: 288 If the matrix Y is missing some diagonal entries this routine can be very slow. To make it fast one should initially 289 fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an 290 entry). No operation is performed when a is zero. 291 292 To form Y = Y + diag(V) use MatDiagonalSet() 293 294 .seealso: MatDiagonalSet(), MatScale(), MatDiagonalScale() 295 @*/ 296 PetscErrorCode MatShift(Mat Y,PetscScalar a) 297 { 298 PetscErrorCode ierr; 299 300 PetscFunctionBegin; 301 PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 302 if (!Y->assembled) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 303 if (Y->factortype) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 304 MatCheckPreallocated(Y,1); 305 if (a == 0.0) PetscFunctionReturn(0); 306 307 if (Y->ops->shift) { 308 ierr = (*Y->ops->shift)(Y,a);CHKERRQ(ierr); 309 } else { 310 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 311 } 312 313 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 314 PetscFunctionReturn(0); 315 } 316 317 PetscErrorCode MatDiagonalSet_Default(Mat Y,Vec D,InsertMode is) 318 { 319 PetscErrorCode ierr; 320 PetscInt i,start,end; 321 const PetscScalar *v; 322 323 PetscFunctionBegin; 324 ierr = MatGetOwnershipRange(Y,&start,&end);CHKERRQ(ierr); 325 ierr = VecGetArrayRead(D,&v);CHKERRQ(ierr); 326 for (i=start; i<end; i++) { 327 ierr = MatSetValues(Y,1,&i,1,&i,v+i-start,is);CHKERRQ(ierr); 328 } 329 ierr = VecRestoreArrayRead(D,&v);CHKERRQ(ierr); 330 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 331 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 332 PetscFunctionReturn(0); 333 } 334 335 /*@ 336 MatDiagonalSet - Computes Y = Y + D, where D is a diagonal matrix 337 that is represented as a vector. Or Y[i,i] = D[i] if InsertMode is 338 INSERT_VALUES. 339 340 Neighbor-wise Collective on Mat 341 342 Input Parameters: 343 + Y - the input matrix 344 . D - the diagonal matrix, represented as a vector 345 - i - INSERT_VALUES or ADD_VALUES 346 347 Notes: 348 If the matrix Y is missing some diagonal entries this routine can be very slow. To make it fast one should initially 349 fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an 350 entry). 351 352 Level: intermediate 353 354 .seealso: MatShift(), MatScale(), MatDiagonalScale() 355 @*/ 356 PetscErrorCode MatDiagonalSet(Mat Y,Vec D,InsertMode is) 357 { 358 PetscErrorCode ierr; 359 PetscInt matlocal,veclocal; 360 361 PetscFunctionBegin; 362 PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 363 PetscValidHeaderSpecific(D,VEC_CLASSID,2); 364 ierr = MatGetLocalSize(Y,&matlocal,NULL);CHKERRQ(ierr); 365 ierr = VecGetLocalSize(D,&veclocal);CHKERRQ(ierr); 366 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); 367 if (Y->ops->diagonalset) { 368 ierr = (*Y->ops->diagonalset)(Y,D,is);CHKERRQ(ierr); 369 } else { 370 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 371 } 372 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 373 PetscFunctionReturn(0); 374 } 375 376 /*@ 377 MatAYPX - Computes Y = a*Y + X. 378 379 Logically on Mat 380 381 Input Parameters: 382 + a - the PetscScalar multiplier 383 . Y - the first matrix 384 . X - the second matrix 385 - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN or SUBSET_NONZERO_PATTERN 386 387 Level: intermediate 388 389 .seealso: MatAXPY() 390 @*/ 391 PetscErrorCode MatAYPX(Mat Y,PetscScalar a,Mat X,MatStructure str) 392 { 393 PetscErrorCode ierr; 394 395 PetscFunctionBegin; 396 ierr = MatScale(Y,a);CHKERRQ(ierr); 397 ierr = MatAXPY(Y,1.0,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 Mat a; 483 PetscScalar *newVals; 484 PetscInt *newCols; 485 PetscInt rStart, rEnd, numRows, maxRows, r, colMax = 0; 486 PetscBool flg; 487 PetscErrorCode ierr; 488 489 PetscFunctionBegin; 490 ierr = PetscObjectBaseTypeCompareAny((PetscObject)A, &flg, MATSEQDENSE, MATMPIDENSE, "");CHKERRQ(ierr); 491 if (flg) { 492 ierr = MatDenseGetLocalMatrix(A, &a);CHKERRQ(ierr); 493 ierr = MatDenseGetLDA(a, &r);CHKERRQ(ierr); 494 ierr = MatGetSize(a, &rStart, &rEnd);CHKERRQ(ierr); 495 ierr = MatDenseGetArray(a, &newVals);CHKERRQ(ierr); 496 for (; colMax < rEnd; ++colMax) { 497 for (maxRows = 0; maxRows < rStart; ++maxRows) { 498 newVals[maxRows + colMax * r] = PetscAbsScalar(newVals[maxRows + colMax * r]) < tol ? 0.0 : newVals[maxRows + colMax * r]; 499 } 500 } 501 ierr = MatDenseRestoreArray(a, &newVals);CHKERRQ(ierr); 502 } else { 503 ierr = MatGetOwnershipRange(A, &rStart, &rEnd);CHKERRQ(ierr); 504 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 505 for (r = rStart; r < rEnd; ++r) { 506 PetscInt ncols; 507 508 ierr = MatGetRow(A, r, &ncols, NULL, NULL);CHKERRQ(ierr); 509 colMax = PetscMax(colMax, ncols);CHKERRQ(ierr); 510 ierr = MatRestoreRow(A, r, &ncols, NULL, NULL);CHKERRQ(ierr); 511 } 512 numRows = rEnd - rStart; 513 ierr = MPIU_Allreduce(&numRows, &maxRows, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 514 ierr = PetscMalloc2(colMax,&newCols,colMax,&newVals);CHKERRQ(ierr); 515 ierr = MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);CHKERRQ(ierr); 516 for (r = rStart; r < rStart+maxRows; ++r) { 517 const PetscScalar *vals; 518 const PetscInt *cols; 519 PetscInt ncols, newcols, c; 520 521 if (r < rEnd) { 522 ierr = MatGetRow(A, r, &ncols, &cols, &vals);CHKERRQ(ierr); 523 for (c = 0; c < ncols; ++c) { 524 newCols[c] = cols[c]; 525 newVals[c] = PetscAbsScalar(vals[c]) < tol ? 0.0 : vals[c]; 526 } 527 newcols = ncols; 528 ierr = MatRestoreRow(A, r, &ncols, &cols, &vals);CHKERRQ(ierr); 529 ierr = MatSetValues(A, 1, &r, newcols, newCols, newVals, INSERT_VALUES);CHKERRQ(ierr); 530 } 531 ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 532 ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 533 } 534 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 535 ierr = PetscFree2(newCols,newVals);CHKERRQ(ierr); 536 } 537 PetscFunctionReturn(0); 538 } 539