xref: /petsc/src/mat/utils/axpy.c (revision dbbe0bcd3f3a8fbab5a45420dc06f8387e5764c6)
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