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 /* TODO mat_tests-ex37_nsize-1_mat_type-baij_mat_block_size-2 fails with MatHeaderMerge */ 220 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 221 } 222 PetscFunctionReturn(0); 223 } 224 225 PetscErrorCode MatAXPY_BasicWithPreallocation(Mat B,Mat Y,PetscScalar a,Mat X,MatStructure str) 226 { 227 PetscInt i,start,end,j,ncols,m,n; 228 PetscErrorCode ierr; 229 const PetscInt *row; 230 PetscScalar *val; 231 const PetscScalar *vals; 232 233 PetscFunctionBegin; 234 ierr = MatGetSize(X,&m,&n);CHKERRQ(ierr); 235 ierr = MatGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 236 ierr = MatGetRowUpperTriangular(Y);CHKERRQ(ierr); 237 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 238 if (a == 1.0) { 239 for (i = start; i < end; i++) { 240 ierr = MatGetRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 241 ierr = MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 242 ierr = MatRestoreRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 243 244 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 245 ierr = MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 246 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 247 } 248 } else { 249 PetscInt vs = 100; 250 /* realloc if needed, as this function may be used in parallel */ 251 ierr = PetscMalloc1(vs,&val);CHKERRQ(ierr); 252 for (i=start; i<end; i++) { 253 ierr = MatGetRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 254 ierr = MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 255 ierr = MatRestoreRow(Y,i,&ncols,&row,&vals);CHKERRQ(ierr); 256 257 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 258 if (vs < ncols) { 259 vs = PetscMin(2*ncols,n); 260 ierr = PetscRealloc(vs*sizeof(*val),&val);CHKERRQ(ierr); 261 } 262 for (j=0; j<ncols; j++) val[j] = a*vals[j]; 263 ierr = MatSetValues(B,1,&i,ncols,row,val,ADD_VALUES);CHKERRQ(ierr); 264 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 265 } 266 ierr = PetscFree(val);CHKERRQ(ierr); 267 } 268 ierr = MatRestoreRowUpperTriangular(Y);CHKERRQ(ierr); 269 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 270 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 271 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 272 PetscFunctionReturn(0); 273 } 274 275 /*@ 276 MatShift - Computes Y = Y + a I, where a is a PetscScalar and I is the identity matrix. 277 278 Neighbor-wise Collective on Mat 279 280 Input Parameters: 281 + Y - the matrices 282 - a - the PetscScalar 283 284 Level: intermediate 285 286 Notes: 287 If the matrix Y is missing some diagonal entries this routine can be very slow. To make it fast one should initially 288 fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an 289 entry). No operation is performed when a is zero. 290 291 To form Y = Y + diag(V) use MatDiagonalSet() 292 293 .seealso: MatDiagonalSet(), MatScale(), MatDiagonalScale() 294 @*/ 295 PetscErrorCode MatShift(Mat Y,PetscScalar a) 296 { 297 PetscErrorCode ierr; 298 299 PetscFunctionBegin; 300 PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 301 if (!Y->assembled) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 302 if (Y->factortype) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 303 MatCheckPreallocated(Y,1); 304 if (a == 0.0) PetscFunctionReturn(0); 305 306 if (Y->ops->shift) { 307 ierr = (*Y->ops->shift)(Y,a);CHKERRQ(ierr); 308 } else { 309 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 310 } 311 312 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 313 PetscFunctionReturn(0); 314 } 315 316 PetscErrorCode MatDiagonalSet_Default(Mat Y,Vec D,InsertMode is) 317 { 318 PetscErrorCode ierr; 319 PetscInt i,start,end; 320 const PetscScalar *v; 321 322 PetscFunctionBegin; 323 ierr = MatGetOwnershipRange(Y,&start,&end);CHKERRQ(ierr); 324 ierr = VecGetArrayRead(D,&v);CHKERRQ(ierr); 325 for (i=start; i<end; i++) { 326 ierr = MatSetValues(Y,1,&i,1,&i,v+i-start,is);CHKERRQ(ierr); 327 } 328 ierr = VecRestoreArrayRead(D,&v);CHKERRQ(ierr); 329 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 330 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 331 PetscFunctionReturn(0); 332 } 333 334 /*@ 335 MatDiagonalSet - Computes Y = Y + D, where D is a diagonal matrix 336 that is represented as a vector. Or Y[i,i] = D[i] if InsertMode is 337 INSERT_VALUES. 338 339 Neighbor-wise Collective on Mat 340 341 Input Parameters: 342 + Y - the input matrix 343 . D - the diagonal matrix, represented as a vector 344 - i - INSERT_VALUES or ADD_VALUES 345 346 Notes: 347 If the matrix Y is missing some diagonal entries this routine can be very slow. To make it fast one should initially 348 fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an 349 entry). 350 351 Level: intermediate 352 353 .seealso: MatShift(), MatScale(), MatDiagonalScale() 354 @*/ 355 PetscErrorCode MatDiagonalSet(Mat Y,Vec D,InsertMode is) 356 { 357 PetscErrorCode ierr; 358 PetscInt matlocal,veclocal; 359 360 PetscFunctionBegin; 361 PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 362 PetscValidHeaderSpecific(D,VEC_CLASSID,2); 363 ierr = MatGetLocalSize(Y,&matlocal,NULL);CHKERRQ(ierr); 364 ierr = VecGetLocalSize(D,&veclocal);CHKERRQ(ierr); 365 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); 366 if (Y->ops->diagonalset) { 367 ierr = (*Y->ops->diagonalset)(Y,D,is);CHKERRQ(ierr); 368 } else { 369 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 370 } 371 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 372 PetscFunctionReturn(0); 373 } 374 375 /*@ 376 MatAYPX - Computes Y = a*Y + X. 377 378 Logically on Mat 379 380 Input Parameters: 381 + a - the PetscScalar multiplier 382 . Y - the first matrix 383 . X - the second matrix 384 - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN, UNKNOWN_NONZERO_PATTERN, or SUBSET_NONZERO_PATTERN (nonzeros of X is a subset of Y's) 385 386 Level: intermediate 387 388 .seealso: MatAXPY() 389 @*/ 390 PetscErrorCode MatAYPX(Mat Y,PetscScalar a,Mat X,MatStructure str) 391 { 392 PetscErrorCode ierr; 393 394 PetscFunctionBegin; 395 ierr = MatScale(Y,a);CHKERRQ(ierr); 396 ierr = MatAXPY(Y,1.0,X,str);CHKERRQ(ierr); 397 PetscFunctionReturn(0); 398 } 399 400 /*@ 401 MatComputeOperator - Computes the explicit matrix 402 403 Collective on Mat 404 405 Input Parameter: 406 + inmat - the matrix 407 - mattype - the matrix type for the explicit operator 408 409 Output Parameter: 410 . mat - the explict operator 411 412 Notes: 413 This computation is done by applying the operators to columns of the identity matrix. 414 This routine is costly in general, and is recommended for use only with relatively small systems. 415 Currently, this routine uses a dense matrix format if mattype == NULL. 416 417 Level: advanced 418 419 @*/ 420 PetscErrorCode MatComputeOperator(Mat inmat,MatType mattype,Mat *mat) 421 { 422 PetscErrorCode ierr; 423 424 PetscFunctionBegin; 425 PetscValidHeaderSpecific(inmat,MAT_CLASSID,1); 426 PetscValidPointer(mat,3); 427 ierr = MatConvert_Shell(inmat,mattype ? mattype : MATDENSE,MAT_INITIAL_MATRIX,mat);CHKERRQ(ierr); 428 PetscFunctionReturn(0); 429 } 430 431 /*@ 432 MatComputeOperatorTranspose - Computes the explicit matrix representation of 433 a give matrix that can apply MatMultTranspose() 434 435 Collective on Mat 436 437 Input Parameter: 438 . inmat - the matrix 439 440 Output Parameter: 441 . mat - the explict 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