16f79c3a4SBarry Smith 2af0996ceSBarry Smith #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 36f79c3a4SBarry Smith 4d54c938cSPierre Jolivet static PetscErrorCode MatTransposeAXPY_Private(Mat Y,PetscScalar a,Mat X,MatStructure str,Mat T) 5d54c938cSPierre Jolivet { 6d54c938cSPierre Jolivet Mat A,F; 75f80ce2aSJacob Faibussowitsch PetscErrorCode (*f)(Mat,Mat*); 8d54c938cSPierre Jolivet 9d54c938cSPierre Jolivet PetscFunctionBegin; 109566063dSJacob Faibussowitsch PetscCall(PetscObjectQueryFunction((PetscObject)T,"MatTransposeGetMat_C",&f)); 11d54c938cSPierre Jolivet if (f) { 129566063dSJacob Faibussowitsch PetscCall(MatTransposeGetMat(T,&A)); 13d54c938cSPierre Jolivet if (T == X) { 149566063dSJacob Faibussowitsch PetscCall(PetscInfo(NULL,"Explicitly transposing X of type MATTRANSPOSEMAT to perform MatAXPY()\n")); 159566063dSJacob Faibussowitsch PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX,&F)); 16d54c938cSPierre Jolivet A = Y; 17d54c938cSPierre Jolivet } else { 189566063dSJacob Faibussowitsch PetscCall(PetscInfo(NULL,"Transposing X because Y of type MATTRANSPOSEMAT to perform MatAXPY()\n")); 199566063dSJacob Faibussowitsch PetscCall(MatTranspose(X,MAT_INITIAL_MATRIX,&F)); 20d54c938cSPierre Jolivet } 21d54c938cSPierre Jolivet } else { 229566063dSJacob Faibussowitsch PetscCall(MatHermitianTransposeGetMat(T,&A)); 23d54c938cSPierre Jolivet if (T == X) { 249566063dSJacob Faibussowitsch PetscCall(PetscInfo(NULL,"Explicitly Hermitian transposing X of type MATTRANSPOSEMAT to perform MatAXPY()\n")); 259566063dSJacob Faibussowitsch PetscCall(MatHermitianTranspose(A,MAT_INITIAL_MATRIX,&F)); 26d54c938cSPierre Jolivet A = Y; 27d54c938cSPierre Jolivet } else { 289566063dSJacob Faibussowitsch PetscCall(PetscInfo(NULL,"Hermitian transposing X because Y of type MATTRANSPOSEMAT to perform MatAXPY()\n")); 299566063dSJacob Faibussowitsch PetscCall(MatHermitianTranspose(X,MAT_INITIAL_MATRIX,&F)); 30d54c938cSPierre Jolivet } 31d54c938cSPierre Jolivet } 329566063dSJacob Faibussowitsch PetscCall(MatAXPY(A,a,F,str)); 339566063dSJacob Faibussowitsch PetscCall(MatDestroy(&F)); 34d54c938cSPierre Jolivet PetscFunctionReturn(0); 35d54c938cSPierre Jolivet } 36d54c938cSPierre Jolivet 3706be10caSBarry Smith /*@ 3821c89e3eSBarry Smith MatAXPY - Computes Y = a*X + Y. 396f79c3a4SBarry Smith 403f9fe445SBarry Smith Logically Collective on Mat 41fee21e36SBarry Smith 4298a79cdbSBarry Smith Input Parameters: 43607cd303SBarry Smith + a - the scalar multiplier 44607cd303SBarry Smith . X - the first matrix 45607cd303SBarry Smith . Y - the second matrix 46531d0c6aSBarry Smith - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN, UNKNOWN_NONZERO_PATTERN, or SUBSET_NONZERO_PATTERN (nonzeros of X is a subset of Y's) 4798a79cdbSBarry Smith 48edf9a46cSStefano Zampini Notes: No operation is performed when a is zero. 49edf9a46cSStefano Zampini 502860a424SLois Curfman McInnes Level: intermediate 512860a424SLois Curfman McInnes 52db781477SPatrick Sanan .seealso: `MatAYPX()` 5306be10caSBarry Smith @*/ 547087cfbeSBarry Smith PetscErrorCode MatAXPY(Mat Y,PetscScalar a,Mat X,MatStructure str) 556f79c3a4SBarry Smith { 56646531bbSStefano Zampini PetscInt M1,M2,N1,N2; 57c1ac3661SBarry Smith PetscInt m1,m2,n1,n2; 58a683ea4aSPierre Jolivet MatType t1,t2; 59a683ea4aSPierre Jolivet PetscBool sametype,transpose; 606f79c3a4SBarry Smith 613a40ed3dSBarry Smith PetscFunctionBegin; 620700a824SBarry Smith PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 63c5eb9154SBarry Smith PetscValidLogicalCollectiveScalar(Y,a,2); 64646531bbSStefano Zampini PetscValidHeaderSpecific(X,MAT_CLASSID,3); 65646531bbSStefano Zampini PetscCheckSameComm(Y,1,X,3); 669566063dSJacob Faibussowitsch PetscCall(MatGetSize(X,&M1,&N1)); 679566063dSJacob Faibussowitsch PetscCall(MatGetSize(Y,&M2,&N2)); 689566063dSJacob Faibussowitsch PetscCall(MatGetLocalSize(X,&m1,&n1)); 699566063dSJacob Faibussowitsch PetscCall(MatGetLocalSize(Y,&m2,&n2)); 702c71b3e2SJacob Faibussowitsch PetscCheck(M1 == M2 && N1 == N2,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_SIZ,"Non conforming matrix add: global sizes X %" PetscInt_FMT " x %" PetscInt_FMT ", Y %" PetscInt_FMT " x %" PetscInt_FMT,M1,N1,M2,N2); 712c71b3e2SJacob Faibussowitsch PetscCheck(m1 == m2 && n1 == n2,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Non conforming matrix add: local sizes X %" PetscInt_FMT " x %" PetscInt_FMT ", Y %" PetscInt_FMT " x %" PetscInt_FMT,m1,n1,m2,n2); 722c71b3e2SJacob Faibussowitsch PetscCheck(Y->assembled,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix (Y)"); 732c71b3e2SJacob Faibussowitsch PetscCheck(X->assembled,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix (X)"); 74edf9a46cSStefano Zampini if (a == (PetscScalar)0.0) PetscFunctionReturn(0); 75edf9a46cSStefano Zampini if (Y == X) { 769566063dSJacob Faibussowitsch PetscCall(MatScale(Y,1.0 + a)); 77edf9a46cSStefano Zampini PetscFunctionReturn(0); 78edf9a46cSStefano Zampini } 799566063dSJacob Faibussowitsch PetscCall(MatGetType(X,&t1)); 809566063dSJacob Faibussowitsch PetscCall(MatGetType(Y,&t2)); 819566063dSJacob Faibussowitsch PetscCall(PetscStrcmp(t1,t2,&sametype)); 829566063dSJacob Faibussowitsch PetscCall(PetscLogEventBegin(MAT_AXPY,Y,0,0,0)); 8393c87e65SStefano Zampini if (Y->ops->axpy && (sametype || X->ops->axpy == Y->ops->axpy)) { 84*dbbe0bcdSBarry Smith PetscUseTypeMethod(Y,axpy,a,X,str); 85d4bb536fSBarry Smith } else { 869566063dSJacob Faibussowitsch PetscCall(PetscStrcmp(t1,MATTRANSPOSEMAT,&transpose)); 87a683ea4aSPierre Jolivet if (transpose) { 889566063dSJacob Faibussowitsch PetscCall(MatTransposeAXPY_Private(Y,a,X,str,X)); 89a683ea4aSPierre Jolivet } else { 909566063dSJacob Faibussowitsch PetscCall(PetscStrcmp(t2,MATTRANSPOSEMAT,&transpose)); 91a683ea4aSPierre Jolivet if (transpose) { 929566063dSJacob Faibussowitsch PetscCall(MatTransposeAXPY_Private(Y,a,X,str,Y)); 93a683ea4aSPierre Jolivet } else { 949566063dSJacob Faibussowitsch PetscCall(MatAXPY_Basic(Y,a,X,str)); 95607cd303SBarry Smith } 96a683ea4aSPierre Jolivet } 97a683ea4aSPierre Jolivet } 989566063dSJacob Faibussowitsch PetscCall(PetscLogEventEnd(MAT_AXPY,Y,0,0,0)); 99607cd303SBarry Smith PetscFunctionReturn(0); 100607cd303SBarry Smith } 101607cd303SBarry Smith 102646531bbSStefano Zampini PetscErrorCode MatAXPY_Basic_Preallocate(Mat Y, Mat X, Mat *B) 103646531bbSStefano Zampini { 10427afe9e8SStefano Zampini PetscErrorCode (*preall)(Mat,Mat,Mat*) = NULL; 105646531bbSStefano Zampini 106646531bbSStefano Zampini PetscFunctionBegin; 107646531bbSStefano Zampini /* look for any available faster alternative to the general preallocator */ 1089566063dSJacob Faibussowitsch PetscCall(PetscObjectQueryFunction((PetscObject)Y,"MatAXPYGetPreallocation_C",&preall)); 109646531bbSStefano Zampini if (preall) { 1109566063dSJacob Faibussowitsch PetscCall((*preall)(Y,X,B)); 111646531bbSStefano Zampini } else { /* Use MatPrellocator, assumes same row-col distribution */ 112646531bbSStefano Zampini Mat preallocator; 113646531bbSStefano Zampini PetscInt r,rstart,rend; 114646531bbSStefano Zampini PetscInt m,n,M,N; 115646531bbSStefano Zampini 1169566063dSJacob Faibussowitsch PetscCall(MatGetRowUpperTriangular(Y)); 1179566063dSJacob Faibussowitsch PetscCall(MatGetRowUpperTriangular(X)); 1189566063dSJacob Faibussowitsch PetscCall(MatGetSize(Y,&M,&N)); 1199566063dSJacob Faibussowitsch PetscCall(MatGetLocalSize(Y,&m,&n)); 1209566063dSJacob Faibussowitsch PetscCall(MatCreate(PetscObjectComm((PetscObject)Y),&preallocator)); 1219566063dSJacob Faibussowitsch PetscCall(MatSetType(preallocator,MATPREALLOCATOR)); 1229566063dSJacob Faibussowitsch PetscCall(MatSetLayouts(preallocator,Y->rmap,Y->cmap)); 1239566063dSJacob Faibussowitsch PetscCall(MatSetUp(preallocator)); 1249566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRange(preallocator,&rstart,&rend)); 125646531bbSStefano Zampini for (r = rstart; r < rend; ++r) { 126646531bbSStefano Zampini PetscInt ncols; 127646531bbSStefano Zampini const PetscInt *row; 128646531bbSStefano Zampini const PetscScalar *vals; 129646531bbSStefano Zampini 1309566063dSJacob Faibussowitsch PetscCall(MatGetRow(Y,r,&ncols,&row,&vals)); 1319566063dSJacob Faibussowitsch PetscCall(MatSetValues(preallocator,1,&r,ncols,row,vals,INSERT_VALUES)); 1329566063dSJacob Faibussowitsch PetscCall(MatRestoreRow(Y,r,&ncols,&row,&vals)); 1339566063dSJacob Faibussowitsch PetscCall(MatGetRow(X,r,&ncols,&row,&vals)); 1349566063dSJacob Faibussowitsch PetscCall(MatSetValues(preallocator,1,&r,ncols,row,vals,INSERT_VALUES)); 1359566063dSJacob Faibussowitsch PetscCall(MatRestoreRow(X,r,&ncols,&row,&vals)); 136646531bbSStefano Zampini } 1379566063dSJacob Faibussowitsch PetscCall(MatSetOption(preallocator,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE)); 1389566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(preallocator,MAT_FINAL_ASSEMBLY)); 1399566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(preallocator,MAT_FINAL_ASSEMBLY)); 1409566063dSJacob Faibussowitsch PetscCall(MatRestoreRowUpperTriangular(Y)); 1419566063dSJacob Faibussowitsch PetscCall(MatRestoreRowUpperTriangular(X)); 142646531bbSStefano Zampini 1439566063dSJacob Faibussowitsch PetscCall(MatCreate(PetscObjectComm((PetscObject)Y),B)); 1449566063dSJacob Faibussowitsch PetscCall(MatSetType(*B,((PetscObject)Y)->type_name)); 1459566063dSJacob Faibussowitsch PetscCall(MatSetLayouts(*B,Y->rmap,Y->cmap)); 1469566063dSJacob Faibussowitsch PetscCall(MatPreallocatorPreallocate(preallocator,PETSC_FALSE,*B)); 1479566063dSJacob Faibussowitsch PetscCall(MatDestroy(&preallocator)); 148646531bbSStefano Zampini } 149646531bbSStefano Zampini PetscFunctionReturn(0); 150646531bbSStefano Zampini } 151646531bbSStefano Zampini 152f4df32b1SMatthew Knepley PetscErrorCode MatAXPY_Basic(Mat Y,PetscScalar a,Mat X,MatStructure str) 153607cd303SBarry Smith { 154be705e3aSPierre Jolivet PetscBool isshell,isdense,isnest; 155edf9a46cSStefano Zampini 156edf9a46cSStefano Zampini PetscFunctionBegin; 1579566063dSJacob Faibussowitsch PetscCall(MatIsShell(Y,&isshell)); 158edf9a46cSStefano Zampini if (isshell) { /* MatShell has special support for AXPY */ 159edf9a46cSStefano Zampini PetscErrorCode (*f)(Mat,PetscScalar,Mat,MatStructure); 160edf9a46cSStefano Zampini 1619566063dSJacob Faibussowitsch PetscCall(MatGetOperation(Y,MATOP_AXPY,(void (**)(void))&f)); 162edf9a46cSStefano Zampini if (f) { 1639566063dSJacob Faibussowitsch PetscCall((*f)(Y,a,X,str)); 164edf9a46cSStefano Zampini PetscFunctionReturn(0); 165edf9a46cSStefano Zampini } 166edf9a46cSStefano Zampini } 16793c87e65SStefano Zampini /* no need to preallocate if Y is dense */ 1689566063dSJacob Faibussowitsch PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)Y,&isdense,MATSEQDENSE,MATMPIDENSE,"")); 169be705e3aSPierre Jolivet if (isdense) { 1709566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)X,MATNEST,&isnest)); 171be705e3aSPierre Jolivet if (isnest) { 1729566063dSJacob Faibussowitsch PetscCall(MatAXPY_Dense_Nest(Y,a,X)); 173be705e3aSPierre Jolivet PetscFunctionReturn(0); 174be705e3aSPierre Jolivet } 175be705e3aSPierre Jolivet if (str == DIFFERENT_NONZERO_PATTERN || str == UNKNOWN_NONZERO_PATTERN) str = SUBSET_NONZERO_PATTERN; 176be705e3aSPierre Jolivet } 177a6ab3590SBarry Smith if (str != DIFFERENT_NONZERO_PATTERN && str != UNKNOWN_NONZERO_PATTERN) { 178edf9a46cSStefano Zampini PetscInt i,start,end,j,ncols,m,n; 17938baddfdSBarry Smith const PetscInt *row; 180b3cc6726SBarry Smith PetscScalar *val; 181b3cc6726SBarry Smith const PetscScalar *vals; 182607cd303SBarry Smith 1839566063dSJacob Faibussowitsch PetscCall(MatGetSize(X,&m,&n)); 1849566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRange(X,&start,&end)); 1859566063dSJacob Faibussowitsch PetscCall(MatGetRowUpperTriangular(X)); 186f4df32b1SMatthew Knepley if (a == 1.0) { 187d4bb536fSBarry Smith for (i = start; i < end; i++) { 1889566063dSJacob Faibussowitsch PetscCall(MatGetRow(X,i,&ncols,&row,&vals)); 1899566063dSJacob Faibussowitsch PetscCall(MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES)); 1909566063dSJacob Faibussowitsch PetscCall(MatRestoreRow(X,i,&ncols,&row,&vals)); 191d4bb536fSBarry Smith } 192d4bb536fSBarry Smith } else { 19321a3365eSStefano Zampini PetscInt vs = 100; 19421a3365eSStefano Zampini /* realloc if needed, as this function may be used in parallel */ 1959566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(vs,&val)); 19606be10caSBarry Smith for (i=start; i<end; i++) { 1979566063dSJacob Faibussowitsch PetscCall(MatGetRow(X,i,&ncols,&row,&vals)); 19821a3365eSStefano Zampini if (vs < ncols) { 19921a3365eSStefano Zampini vs = PetscMin(2*ncols,n); 2009566063dSJacob Faibussowitsch PetscCall(PetscRealloc(vs*sizeof(*val),&val)); 2016f79c3a4SBarry Smith } 20221a3365eSStefano Zampini for (j=0; j<ncols; j++) val[j] = a*vals[j]; 2039566063dSJacob Faibussowitsch PetscCall(MatSetValues(Y,1,&i,ncols,row,val,ADD_VALUES)); 2049566063dSJacob Faibussowitsch PetscCall(MatRestoreRow(X,i,&ncols,&row,&vals)); 2056f79c3a4SBarry Smith } 2069566063dSJacob Faibussowitsch PetscCall(PetscFree(val)); 207d4bb536fSBarry Smith } 2089566063dSJacob Faibussowitsch PetscCall(MatRestoreRowUpperTriangular(X)); 2099566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY)); 2109566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY)); 211edf9a46cSStefano Zampini } else { 212edf9a46cSStefano Zampini Mat B; 213edf9a46cSStefano Zampini 2149566063dSJacob Faibussowitsch PetscCall(MatAXPY_Basic_Preallocate(Y,X,&B)); 2159566063dSJacob Faibussowitsch PetscCall(MatAXPY_BasicWithPreallocation(B,Y,a,X,str)); 2169566063dSJacob Faibussowitsch PetscCall(MatHeaderMerge(Y,&B)); 217edf9a46cSStefano Zampini } 2183a40ed3dSBarry Smith PetscFunctionReturn(0); 2196f79c3a4SBarry Smith } 220052efed2SBarry Smith 221ec7775f6SShri Abhyankar PetscErrorCode MatAXPY_BasicWithPreallocation(Mat B,Mat Y,PetscScalar a,Mat X,MatStructure str) 222ec7775f6SShri Abhyankar { 223ec7775f6SShri Abhyankar PetscInt i,start,end,j,ncols,m,n; 224ec7775f6SShri Abhyankar const PetscInt *row; 225ec7775f6SShri Abhyankar PetscScalar *val; 226ec7775f6SShri Abhyankar const PetscScalar *vals; 227ec7775f6SShri Abhyankar 228ec7775f6SShri Abhyankar PetscFunctionBegin; 2299566063dSJacob Faibussowitsch PetscCall(MatGetSize(X,&m,&n)); 2309566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRange(X,&start,&end)); 2319566063dSJacob Faibussowitsch PetscCall(MatGetRowUpperTriangular(Y)); 2329566063dSJacob Faibussowitsch PetscCall(MatGetRowUpperTriangular(X)); 233ec7775f6SShri Abhyankar if (a == 1.0) { 234ec7775f6SShri Abhyankar for (i = start; i < end; i++) { 2359566063dSJacob Faibussowitsch PetscCall(MatGetRow(Y,i,&ncols,&row,&vals)); 2369566063dSJacob Faibussowitsch PetscCall(MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES)); 2379566063dSJacob Faibussowitsch PetscCall(MatRestoreRow(Y,i,&ncols,&row,&vals)); 238ec7775f6SShri Abhyankar 2399566063dSJacob Faibussowitsch PetscCall(MatGetRow(X,i,&ncols,&row,&vals)); 2409566063dSJacob Faibussowitsch PetscCall(MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES)); 2419566063dSJacob Faibussowitsch PetscCall(MatRestoreRow(X,i,&ncols,&row,&vals)); 242ec7775f6SShri Abhyankar } 243ec7775f6SShri Abhyankar } else { 24421a3365eSStefano Zampini PetscInt vs = 100; 24521a3365eSStefano Zampini /* realloc if needed, as this function may be used in parallel */ 2469566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(vs,&val)); 247ec7775f6SShri Abhyankar for (i=start; i<end; i++) { 2489566063dSJacob Faibussowitsch PetscCall(MatGetRow(Y,i,&ncols,&row,&vals)); 2499566063dSJacob Faibussowitsch PetscCall(MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES)); 2509566063dSJacob Faibussowitsch PetscCall(MatRestoreRow(Y,i,&ncols,&row,&vals)); 251ec7775f6SShri Abhyankar 2529566063dSJacob Faibussowitsch PetscCall(MatGetRow(X,i,&ncols,&row,&vals)); 25321a3365eSStefano Zampini if (vs < ncols) { 25421a3365eSStefano Zampini vs = PetscMin(2*ncols,n); 2559566063dSJacob Faibussowitsch PetscCall(PetscRealloc(vs*sizeof(*val),&val)); 256ec7775f6SShri Abhyankar } 25721a3365eSStefano Zampini for (j=0; j<ncols; j++) val[j] = a*vals[j]; 2589566063dSJacob Faibussowitsch PetscCall(MatSetValues(B,1,&i,ncols,row,val,ADD_VALUES)); 2599566063dSJacob Faibussowitsch PetscCall(MatRestoreRow(X,i,&ncols,&row,&vals)); 260ec7775f6SShri Abhyankar } 2619566063dSJacob Faibussowitsch PetscCall(PetscFree(val)); 262ec7775f6SShri Abhyankar } 2639566063dSJacob Faibussowitsch PetscCall(MatRestoreRowUpperTriangular(Y)); 2649566063dSJacob Faibussowitsch PetscCall(MatRestoreRowUpperTriangular(X)); 2659566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 2669566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 267ec7775f6SShri Abhyankar PetscFunctionReturn(0); 268ec7775f6SShri Abhyankar } 269ec7775f6SShri Abhyankar 270052efed2SBarry Smith /*@ 27187828ca2SBarry Smith MatShift - Computes Y = Y + a I, where a is a PetscScalar and I is the identity matrix. 272052efed2SBarry Smith 2733f9fe445SBarry Smith Neighbor-wise Collective on Mat 274fee21e36SBarry Smith 27598a79cdbSBarry Smith Input Parameters: 27698a79cdbSBarry Smith + Y - the matrices 27787828ca2SBarry Smith - a - the PetscScalar 27898a79cdbSBarry Smith 2792860a424SLois Curfman McInnes Level: intermediate 2802860a424SLois Curfman McInnes 28195452b02SPatrick Sanan Notes: 28284e19e3eSJunchao Zhang If Y is a rectangular matrix, the shift is done on the main diagonal Y_{ii} of the matrix (https://en.wikipedia.org/wiki/Main_diagonal) 28384e19e3eSJunchao Zhang 28495452b02SPatrick Sanan If the matrix Y is missing some diagonal entries this routine can be very slow. To make it fast one should initially 2856f33a894SBarry Smith fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an 286edf9a46cSStefano Zampini entry). No operation is performed when a is zero. 2876f33a894SBarry Smith 2880c0fd78eSBarry Smith To form Y = Y + diag(V) use MatDiagonalSet() 2890c0fd78eSBarry Smith 290db781477SPatrick Sanan .seealso: `MatDiagonalSet()`, `MatScale()`, `MatDiagonalScale()` 291052efed2SBarry Smith @*/ 2927087cfbeSBarry Smith PetscErrorCode MatShift(Mat Y,PetscScalar a) 293052efed2SBarry Smith { 2943a40ed3dSBarry Smith PetscFunctionBegin; 2950700a824SBarry Smith PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 2965f80ce2aSJacob Faibussowitsch PetscCheck(Y->assembled,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2975f80ce2aSJacob Faibussowitsch PetscCheck(!Y->factortype,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2984994cf47SJed Brown MatCheckPreallocated(Y,1); 2991b71c1a7SBarry Smith if (a == 0.0) PetscFunctionReturn(0); 300b50b34bdSBarry Smith 301*dbbe0bcdSBarry Smith if (Y->ops->shift) PetscUseTypeMethod(Y,shift,a); 3029566063dSJacob Faibussowitsch else PetscCall(MatShift_Basic(Y,a)); 3037d68702bSBarry Smith 3049566063dSJacob Faibussowitsch PetscCall(PetscObjectStateIncrease((PetscObject)Y)); 3053a40ed3dSBarry Smith PetscFunctionReturn(0); 306052efed2SBarry Smith } 3076d84be18SBarry Smith 3087087cfbeSBarry Smith PetscErrorCode MatDiagonalSet_Default(Mat Y,Vec D,InsertMode is) 30909f38230SBarry Smith { 31067576b8bSBarry Smith PetscInt i,start,end; 311fa2e578bSStefano Zampini const PetscScalar *v; 31209f38230SBarry Smith 31309f38230SBarry Smith PetscFunctionBegin; 3149566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRange(Y,&start,&end)); 3159566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(D,&v)); 31609f38230SBarry Smith for (i=start; i<end; i++) { 3179566063dSJacob Faibussowitsch PetscCall(MatSetValues(Y,1,&i,1,&i,v+i-start,is)); 31809f38230SBarry Smith } 3199566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(D,&v)); 3209566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY)); 3219566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY)); 32209f38230SBarry Smith PetscFunctionReturn(0); 32309f38230SBarry Smith } 32409f38230SBarry Smith 3256d84be18SBarry Smith /*@ 326f56f2b3fSBarry Smith MatDiagonalSet - Computes Y = Y + D, where D is a diagonal matrix 327f56f2b3fSBarry Smith that is represented as a vector. Or Y[i,i] = D[i] if InsertMode is 328f56f2b3fSBarry Smith INSERT_VALUES. 3296d84be18SBarry Smith 33048e586daSJose E. Roman Neighbor-wise Collective on Mat 33148e586daSJose E. Roman 3326d84be18SBarry Smith Input Parameters: 33398a79cdbSBarry Smith + Y - the input matrix 334f56f2b3fSBarry Smith . D - the diagonal matrix, represented as a vector 335f56f2b3fSBarry Smith - i - INSERT_VALUES or ADD_VALUES 3366d84be18SBarry Smith 33795452b02SPatrick Sanan Notes: 33895452b02SPatrick Sanan If the matrix Y is missing some diagonal entries this routine can be very slow. To make it fast one should initially 3396f33a894SBarry Smith fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an 3406f33a894SBarry Smith entry). 3416f33a894SBarry Smith 3422860a424SLois Curfman McInnes Level: intermediate 3432860a424SLois Curfman McInnes 344db781477SPatrick Sanan .seealso: `MatShift()`, `MatScale()`, `MatDiagonalScale()` 3456d84be18SBarry Smith @*/ 3467087cfbeSBarry Smith PetscErrorCode MatDiagonalSet(Mat Y,Vec D,InsertMode is) 3476d84be18SBarry Smith { 34867576b8bSBarry Smith PetscInt matlocal,veclocal; 3496d84be18SBarry Smith 3503a40ed3dSBarry Smith PetscFunctionBegin; 3510700a824SBarry Smith PetscValidHeaderSpecific(Y,MAT_CLASSID,1); 3520700a824SBarry Smith PetscValidHeaderSpecific(D,VEC_CLASSID,2); 3539566063dSJacob Faibussowitsch PetscCall(MatGetLocalSize(Y,&matlocal,NULL)); 3549566063dSJacob Faibussowitsch PetscCall(VecGetLocalSize(D,&veclocal)); 3555f80ce2aSJacob Faibussowitsch PetscCheck(matlocal == veclocal,PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number local rows of matrix %" PetscInt_FMT " does not match that of vector for diagonal %" PetscInt_FMT,matlocal,veclocal); 356*dbbe0bcdSBarry Smith if (Y->ops->diagonalset) PetscUseTypeMethod(Y,diagonalset,D,is); 357*dbbe0bcdSBarry Smith else PetscCall(MatDiagonalSet_Default(Y,D,is)); 3589566063dSJacob Faibussowitsch PetscCall(PetscObjectStateIncrease((PetscObject)Y)); 3593a40ed3dSBarry Smith PetscFunctionReturn(0); 3606d84be18SBarry Smith } 361d4bb536fSBarry Smith 362d4bb536fSBarry Smith /*@ 36304aac2b0SHong Zhang MatAYPX - Computes Y = a*Y + X. 364d4bb536fSBarry Smith 3653f9fe445SBarry Smith Logically on Mat 366fee21e36SBarry Smith 36798a79cdbSBarry Smith Input Parameters: 36804aac2b0SHong Zhang + a - the PetscScalar multiplier 36904aac2b0SHong Zhang . Y - the first matrix 37004aac2b0SHong Zhang . X - the second matrix 371531d0c6aSBarry Smith - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN, UNKNOWN_NONZERO_PATTERN, or SUBSET_NONZERO_PATTERN (nonzeros of X is a subset of Y's) 372d4bb536fSBarry Smith 3732860a424SLois Curfman McInnes Level: intermediate 3742860a424SLois Curfman McInnes 375db781477SPatrick Sanan .seealso: `MatAXPY()` 376d4bb536fSBarry Smith @*/ 3777087cfbeSBarry Smith PetscErrorCode MatAYPX(Mat Y,PetscScalar a,Mat X,MatStructure str) 378d4bb536fSBarry Smith { 3793a40ed3dSBarry Smith PetscFunctionBegin; 3809566063dSJacob Faibussowitsch PetscCall(MatScale(Y,a)); 3819566063dSJacob Faibussowitsch PetscCall(MatAXPY(Y,1.0,X,str)); 3823a40ed3dSBarry Smith PetscFunctionReturn(0); 383d4bb536fSBarry Smith } 384b0a32e0cSBarry Smith 385b0a32e0cSBarry Smith /*@ 3860bacdadaSStefano Zampini MatComputeOperator - Computes the explicit matrix 387b0a32e0cSBarry Smith 388b0a32e0cSBarry Smith Collective on Mat 389b0a32e0cSBarry Smith 390d8d19677SJose E. Roman Input Parameters: 391186a3e20SStefano Zampini + inmat - the matrix 392186a3e20SStefano Zampini - mattype - the matrix type for the explicit operator 393b0a32e0cSBarry Smith 394b0a32e0cSBarry Smith Output Parameter: 395a5b23f4aSJose E. Roman . mat - the explicit operator 396b0a32e0cSBarry Smith 397b0a32e0cSBarry Smith Notes: 398186a3e20SStefano Zampini This computation is done by applying the operators to columns of the identity matrix. 399186a3e20SStefano Zampini This routine is costly in general, and is recommended for use only with relatively small systems. 400186a3e20SStefano Zampini Currently, this routine uses a dense matrix format if mattype == NULL. 401b0a32e0cSBarry Smith 402b0a32e0cSBarry Smith Level: advanced 403b0a32e0cSBarry Smith 404b0a32e0cSBarry Smith @*/ 4050bacdadaSStefano Zampini PetscErrorCode MatComputeOperator(Mat inmat,MatType mattype,Mat *mat) 406b0a32e0cSBarry Smith { 407b0a32e0cSBarry Smith PetscFunctionBegin; 4080700a824SBarry Smith PetscValidHeaderSpecific(inmat,MAT_CLASSID,1); 409186a3e20SStefano Zampini PetscValidPointer(mat,3); 4109566063dSJacob Faibussowitsch PetscCall(MatConvert_Shell(inmat,mattype ? mattype : MATDENSE,MAT_INITIAL_MATRIX,mat)); 41124f910e3SHong Zhang PetscFunctionReturn(0); 41224f910e3SHong Zhang } 4134325cce7SMatthew G Knepley 4144325cce7SMatthew G Knepley /*@ 4150bacdadaSStefano Zampini MatComputeOperatorTranspose - Computes the explicit matrix representation of 416f3b1f45cSBarry Smith a give matrix that can apply MatMultTranspose() 417f3b1f45cSBarry Smith 418f3b1f45cSBarry Smith Collective on Mat 419f3b1f45cSBarry Smith 4206b867d5aSJose E. Roman Input Parameters: 4216b867d5aSJose E. Roman + inmat - the matrix 4226b867d5aSJose E. Roman - mattype - the matrix type for the explicit operator 423f3b1f45cSBarry Smith 424f3b1f45cSBarry Smith Output Parameter: 425a5b23f4aSJose E. Roman . mat - the explicit operator transposed 426f3b1f45cSBarry Smith 427f3b1f45cSBarry Smith Notes: 428186a3e20SStefano Zampini This computation is done by applying the transpose of the operator to columns of the identity matrix. 429186a3e20SStefano Zampini This routine is costly in general, and is recommended for use only with relatively small systems. 430186a3e20SStefano Zampini Currently, this routine uses a dense matrix format if mattype == NULL. 431f3b1f45cSBarry Smith 432f3b1f45cSBarry Smith Level: advanced 433f3b1f45cSBarry Smith 434f3b1f45cSBarry Smith @*/ 4350bacdadaSStefano Zampini PetscErrorCode MatComputeOperatorTranspose(Mat inmat,MatType mattype,Mat *mat) 436f3b1f45cSBarry Smith { 437186a3e20SStefano Zampini Mat A; 438f3b1f45cSBarry Smith 439f3b1f45cSBarry Smith PetscFunctionBegin; 440f3b1f45cSBarry Smith PetscValidHeaderSpecific(inmat,MAT_CLASSID,1); 441186a3e20SStefano Zampini PetscValidPointer(mat,3); 4429566063dSJacob Faibussowitsch PetscCall(MatCreateTranspose(inmat,&A)); 4439566063dSJacob Faibussowitsch PetscCall(MatConvert_Shell(A,mattype ? mattype : MATDENSE,MAT_INITIAL_MATRIX,mat)); 4449566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A)); 445f3b1f45cSBarry Smith PetscFunctionReturn(0); 446f3b1f45cSBarry Smith } 447f3b1f45cSBarry Smith 448f3b1f45cSBarry Smith /*@ 4494325cce7SMatthew G Knepley MatChop - Set all values in the matrix less than the tolerance to zero 4504325cce7SMatthew G Knepley 4514325cce7SMatthew G Knepley Input Parameters: 4524325cce7SMatthew G Knepley + A - The matrix 4534325cce7SMatthew G Knepley - tol - The zero tolerance 4544325cce7SMatthew G Knepley 4554325cce7SMatthew G Knepley Output Parameters: 4564325cce7SMatthew G Knepley . A - The chopped matrix 4574325cce7SMatthew G Knepley 4584325cce7SMatthew G Knepley Level: intermediate 4594325cce7SMatthew G Knepley 460db781477SPatrick Sanan .seealso: `MatCreate()`, `MatZeroEntries()` 4613fc99919SSatish Balay @*/ 4624325cce7SMatthew G Knepley PetscErrorCode MatChop(Mat A, PetscReal tol) 4634325cce7SMatthew G Knepley { 464038df967SPierre Jolivet Mat a; 4654325cce7SMatthew G Knepley PetscScalar *newVals; 466cf5d3338SPierre Jolivet PetscInt *newCols, rStart, rEnd, numRows, maxRows, r, colMax = 0; 467038df967SPierre Jolivet PetscBool flg; 4684325cce7SMatthew G Knepley 4694325cce7SMatthew G Knepley PetscFunctionBegin; 4709566063dSJacob Faibussowitsch PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)A, &flg, MATSEQDENSE, MATMPIDENSE, "")); 471038df967SPierre Jolivet if (flg) { 4729566063dSJacob Faibussowitsch PetscCall(MatDenseGetLocalMatrix(A, &a)); 4739566063dSJacob Faibussowitsch PetscCall(MatDenseGetLDA(a, &r)); 4749566063dSJacob Faibussowitsch PetscCall(MatGetSize(a, &rStart, &rEnd)); 4759566063dSJacob Faibussowitsch PetscCall(MatDenseGetArray(a, &newVals)); 476038df967SPierre Jolivet for (; colMax < rEnd; ++colMax) { 477038df967SPierre Jolivet for (maxRows = 0; maxRows < rStart; ++maxRows) { 478038df967SPierre Jolivet newVals[maxRows + colMax * r] = PetscAbsScalar(newVals[maxRows + colMax * r]) < tol ? 0.0 : newVals[maxRows + colMax * r]; 479038df967SPierre Jolivet } 480038df967SPierre Jolivet } 4819566063dSJacob Faibussowitsch PetscCall(MatDenseRestoreArray(a, &newVals)); 482038df967SPierre Jolivet } else { 4839566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRange(A, &rStart, &rEnd)); 4849566063dSJacob Faibussowitsch PetscCall(MatGetRowUpperTriangular(A)); 4854325cce7SMatthew G Knepley for (r = rStart; r < rEnd; ++r) { 4864325cce7SMatthew G Knepley PetscInt ncols; 4874325cce7SMatthew G Knepley 4889566063dSJacob Faibussowitsch PetscCall(MatGetRow(A, r, &ncols, NULL, NULL)); 4895f80ce2aSJacob Faibussowitsch colMax = PetscMax(colMax, ncols); 4909566063dSJacob Faibussowitsch PetscCall(MatRestoreRow(A, r, &ncols, NULL, NULL)); 4914325cce7SMatthew G Knepley } 4924325cce7SMatthew G Knepley numRows = rEnd - rStart; 4931c2dc1cbSBarry Smith PetscCall(MPIU_Allreduce(&numRows, &maxRows, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)A))); 4949566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(colMax, &newCols, colMax, &newVals)); 4959566063dSJacob Faibussowitsch PetscCall(MatGetOption(A, MAT_NO_OFF_PROC_ENTRIES, &flg)); /* cache user-defined value */ 4969566063dSJacob Faibussowitsch PetscCall(MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 497cf5d3338SPierre Jolivet /* short-circuit code in MatAssemblyBegin() and MatAssemblyEnd() */ 498cf5d3338SPierre Jolivet /* that are potentially called many times depending on the distribution of A */ 4994325cce7SMatthew G Knepley for (r = rStart; r < rStart+maxRows; ++r) { 5004325cce7SMatthew G Knepley const PetscScalar *vals; 5014325cce7SMatthew G Knepley const PetscInt *cols; 502fad45679SMatthew G. Knepley PetscInt ncols, newcols, c; 5034325cce7SMatthew G Knepley 5044325cce7SMatthew G Knepley if (r < rEnd) { 5059566063dSJacob Faibussowitsch PetscCall(MatGetRow(A, r, &ncols, &cols, &vals)); 5064325cce7SMatthew G Knepley for (c = 0; c < ncols; ++c) { 5074325cce7SMatthew G Knepley newCols[c] = cols[c]; 5084325cce7SMatthew G Knepley newVals[c] = PetscAbsScalar(vals[c]) < tol ? 0.0 : vals[c]; 5094325cce7SMatthew G Knepley } 510fad45679SMatthew G. Knepley newcols = ncols; 5119566063dSJacob Faibussowitsch PetscCall(MatRestoreRow(A, r, &ncols, &cols, &vals)); 5129566063dSJacob Faibussowitsch PetscCall(MatSetValues(A, 1, &r, newcols, newCols, newVals, INSERT_VALUES)); 5134325cce7SMatthew G Knepley } 5149566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 5159566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 5164325cce7SMatthew G Knepley } 5179566063dSJacob Faibussowitsch PetscCall(MatRestoreRowUpperTriangular(A)); 5189566063dSJacob Faibussowitsch PetscCall(PetscFree2(newCols, newVals)); 5199566063dSJacob Faibussowitsch PetscCall(MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, flg)); /* reset option to its user-defined value */ 520038df967SPierre Jolivet } 5214325cce7SMatthew G Knepley PetscFunctionReturn(0); 5224325cce7SMatthew G Knepley } 523