xref: /petsc/src/mat/tests/ex221.c (revision 08401ef684002a709c6d3db98a0c9f54a8bcf1ec)
1c4762a1bSJed Brown static char help[] = "Tests various routines for MATSHELL\n\n";
2c4762a1bSJed Brown 
3c4762a1bSJed Brown #include <petscmat.h>
4c4762a1bSJed Brown 
5c4762a1bSJed Brown typedef struct _n_User *User;
6c4762a1bSJed Brown struct _n_User {
7c4762a1bSJed Brown   Mat B;
8c4762a1bSJed Brown };
9c4762a1bSJed Brown 
10c4762a1bSJed Brown static PetscErrorCode MatGetDiagonal_User(Mat A,Vec X)
11c4762a1bSJed Brown {
12c4762a1bSJed Brown   User           user;
13c4762a1bSJed Brown 
14c4762a1bSJed Brown   PetscFunctionBegin;
159566063dSJacob Faibussowitsch   PetscCall(MatShellGetContext(A,&user));
169566063dSJacob Faibussowitsch   PetscCall(MatGetDiagonal(user->B,X));
17c4762a1bSJed Brown   PetscFunctionReturn(0);
18c4762a1bSJed Brown }
19c4762a1bSJed Brown 
20c4762a1bSJed Brown static PetscErrorCode MatMult_User(Mat A,Vec X,Vec Y)
21c4762a1bSJed Brown {
22c4762a1bSJed Brown   User           user;
23c4762a1bSJed Brown 
24c4762a1bSJed Brown   PetscFunctionBegin;
259566063dSJacob Faibussowitsch   PetscCall(MatShellGetContext(A,&user));
269566063dSJacob Faibussowitsch   PetscCall(MatMult(user->B,X,Y));
27c4762a1bSJed Brown   PetscFunctionReturn(0);
28c4762a1bSJed Brown }
29c4762a1bSJed Brown 
30c4762a1bSJed Brown static PetscErrorCode MatMultTranspose_User(Mat A,Vec X,Vec Y)
31c4762a1bSJed Brown {
32c4762a1bSJed Brown   User           user;
33c4762a1bSJed Brown 
34c4762a1bSJed Brown   PetscFunctionBegin;
359566063dSJacob Faibussowitsch   PetscCall(MatShellGetContext(A,&user));
369566063dSJacob Faibussowitsch   PetscCall(MatMultTranspose(user->B,X,Y));
37c4762a1bSJed Brown   PetscFunctionReturn(0);
38c4762a1bSJed Brown }
39c4762a1bSJed Brown 
40c4762a1bSJed Brown static PetscErrorCode MatCopy_User(Mat A,Mat X,MatStructure str)
41c4762a1bSJed Brown {
42c4762a1bSJed Brown   User           user,userX;
43c4762a1bSJed Brown 
44c4762a1bSJed Brown   PetscFunctionBegin;
459566063dSJacob Faibussowitsch   PetscCall(MatShellGetContext(A,&user));
469566063dSJacob Faibussowitsch   PetscCall(MatShellGetContext(X,&userX));
47*08401ef6SPierre Jolivet   PetscCheck(user == userX,PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen");
489566063dSJacob Faibussowitsch   PetscCall(PetscObjectReference((PetscObject)user->B));
49c4762a1bSJed Brown   PetscFunctionReturn(0);
50c4762a1bSJed Brown }
51c4762a1bSJed Brown 
52c4762a1bSJed Brown static PetscErrorCode MatDestroy_User(Mat A)
53c4762a1bSJed Brown {
54c4762a1bSJed Brown   User           user;
55c4762a1bSJed Brown 
56c4762a1bSJed Brown   PetscFunctionBegin;
579566063dSJacob Faibussowitsch   PetscCall(MatShellGetContext(A,&user));
589566063dSJacob Faibussowitsch   PetscCall(PetscObjectDereference((PetscObject)user->B));
59c4762a1bSJed Brown   PetscFunctionReturn(0);
60c4762a1bSJed Brown }
61c4762a1bSJed Brown 
62c4762a1bSJed Brown int main(int argc,char **args)
63c4762a1bSJed Brown {
64c4762a1bSJed Brown   User           user;
65c4762a1bSJed Brown   Mat            A,S;
66c4762a1bSJed Brown   PetscScalar    *data,diag = 1.3;
67c4762a1bSJed Brown   PetscReal      tol = PETSC_SMALL;
68c4762a1bSJed Brown   PetscInt       i,j,m = PETSC_DECIDE,n = PETSC_DECIDE,M = 17,N = 15,s1,s2;
69c4762a1bSJed Brown   PetscInt       test, ntest = 2;
70c4762a1bSJed Brown   PetscMPIInt    rank,size;
71c4762a1bSJed Brown   PetscBool      nc = PETSC_FALSE, cong, flg;
72c4762a1bSJed Brown   PetscBool      ronl = PETSC_TRUE;
73c4762a1bSJed Brown   PetscBool      randomize = PETSC_FALSE, submat = PETSC_FALSE;
74c4762a1bSJed Brown   PetscBool      keep = PETSC_FALSE;
75c4762a1bSJed Brown   PetscBool      testzerorows = PETSC_TRUE, testdiagscale = PETSC_TRUE, testgetdiag = PETSC_TRUE, testsubmat = PETSC_TRUE;
76c4762a1bSJed Brown   PetscBool      testshift = PETSC_TRUE, testscale = PETSC_TRUE, testdup = PETSC_TRUE, testreset = PETSC_TRUE;
77c4762a1bSJed Brown   PetscBool      testaxpy = PETSC_TRUE, testaxpyd = PETSC_TRUE, testaxpyerr = PETSC_FALSE;
78c4762a1bSJed Brown 
799566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc,&args,(char*)0,help));
809566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
819566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
829566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL));
839566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL));
849566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-ml",&m,NULL));
859566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-nl",&n,NULL));
869566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-square_nc",&nc,NULL));
879566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-rows_only",&ronl,NULL));
889566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-randomize",&randomize,NULL));
899566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-submat",&submat,NULL));
909566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_zerorows",&testzerorows,NULL));
919566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_diagscale",&testdiagscale,NULL));
929566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_getdiag",&testgetdiag,NULL));
939566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_shift",&testshift,NULL));
949566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_scale",&testscale,NULL));
959566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_dup",&testdup,NULL));
969566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_reset",&testreset,NULL));
979566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_submat",&testsubmat,NULL));
989566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_axpy",&testaxpy,NULL));
999566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_axpy_different",&testaxpyd,NULL));
1009566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_axpy_error",&testaxpyerr,NULL));
1019566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-loop",&ntest,NULL));
1029566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL));
1039566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-diag",&diag,NULL));
1049566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-keep",&keep,NULL));
105c4762a1bSJed Brown   /* This tests square matrices with different row/col layout */
106c4762a1bSJed Brown   if (nc && size > 1) {
107c4762a1bSJed Brown     M = PetscMax(PetscMax(N,M),1);
108c4762a1bSJed Brown     N = M;
109c4762a1bSJed Brown     m = n = 0;
110c4762a1bSJed Brown     if (rank == 0) { m = M-1; n = 1; }
111c4762a1bSJed Brown     else if (rank == 1) { m = 1; n = N-1; }
112c4762a1bSJed Brown   }
1139566063dSJacob Faibussowitsch   PetscCall(MatCreateDense(PETSC_COMM_WORLD,m,n,M,N,NULL,&A));
1149566063dSJacob Faibussowitsch   PetscCall(MatGetLocalSize(A,&m,&n));
1159566063dSJacob Faibussowitsch   PetscCall(MatGetSize(A,&M,&N));
1169566063dSJacob Faibussowitsch   PetscCall(MatGetOwnershipRange(A,&s1,NULL));
117c4762a1bSJed Brown   s2   = 1;
118c4762a1bSJed Brown   while (s2 < M) s2 *= 10;
1199566063dSJacob Faibussowitsch   PetscCall(MatDenseGetArray(A,&data));
120c4762a1bSJed Brown   for (j = 0; j < N; j++) {
121c4762a1bSJed Brown     for (i = 0; i < m; i++) {
122c4762a1bSJed Brown       data[j*m + i] = s2*j + i + s1 + 1;
123c4762a1bSJed Brown     }
124c4762a1bSJed Brown   }
1259566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
1269566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
127c4762a1bSJed Brown 
128c4762a1bSJed Brown   if (submat) {
129c4762a1bSJed Brown     Mat      A2;
130c4762a1bSJed Brown     IS       r,c;
131c4762a1bSJed Brown     PetscInt rst,ren,cst,cen;
132c4762a1bSJed Brown 
1339566063dSJacob Faibussowitsch     PetscCall(MatGetOwnershipRange(A,&rst,&ren));
1349566063dSJacob Faibussowitsch     PetscCall(MatGetOwnershipRangeColumn(A,&cst,&cen));
1359566063dSJacob Faibussowitsch     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r));
1369566063dSJacob Faibussowitsch     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c));
1379566063dSJacob Faibussowitsch     PetscCall(MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&A2));
1389566063dSJacob Faibussowitsch     PetscCall(ISDestroy(&r));
1399566063dSJacob Faibussowitsch     PetscCall(ISDestroy(&c));
1409566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&A));
141c4762a1bSJed Brown     A = A2;
142c4762a1bSJed Brown   }
143c4762a1bSJed Brown 
1449566063dSJacob Faibussowitsch   PetscCall(MatGetSize(A,&M,&N));
1459566063dSJacob Faibussowitsch   PetscCall(MatGetLocalSize(A,&m,&n));
1469566063dSJacob Faibussowitsch   PetscCall(MatHasCongruentLayouts(A,&cong));
147c4762a1bSJed Brown 
1489566063dSJacob Faibussowitsch   PetscCall(MatConvert(A,MATAIJ,MAT_INPLACE_MATRIX,&A));
1499566063dSJacob Faibussowitsch   PetscCall(MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,keep));
1509566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)A,"initial"));
1519566063dSJacob Faibussowitsch   PetscCall(MatViewFromOptions(A,NULL,"-view_mat"));
152c4762a1bSJed Brown 
1539566063dSJacob Faibussowitsch   PetscCall(PetscNew(&user));
1549566063dSJacob Faibussowitsch   PetscCall(MatCreateShell(PETSC_COMM_WORLD,m,n,M,N,user,&S));
1559566063dSJacob Faibussowitsch   PetscCall(MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User));
1569566063dSJacob Faibussowitsch   PetscCall(MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User));
157c4762a1bSJed Brown   if (cong) {
1589566063dSJacob Faibussowitsch     PetscCall(MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User));
159c4762a1bSJed Brown   }
1609566063dSJacob Faibussowitsch   PetscCall(MatShellSetOperation(S,MATOP_COPY,(void (*)(void))MatCopy_User));
1619566063dSJacob Faibussowitsch   PetscCall(MatShellSetOperation(S,MATOP_DESTROY,(void (*)(void))MatDestroy_User));
1629566063dSJacob Faibussowitsch   PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&user->B));
163c4762a1bSJed Brown 
164c4762a1bSJed Brown   /* Square and rows only scaling */
165c4762a1bSJed Brown   ronl = cong ? ronl : PETSC_TRUE;
166c4762a1bSJed Brown 
167c4762a1bSJed Brown   for (test = 0; test < ntest; test++) {
168c4762a1bSJed Brown     PetscReal err;
169c4762a1bSJed Brown 
1709566063dSJacob Faibussowitsch     PetscCall(MatMultAddEqual(A,S,10,&flg));
171c4762a1bSJed Brown     if (!flg) {
1729566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mult add\n",test));
173c4762a1bSJed Brown     }
1749566063dSJacob Faibussowitsch     PetscCall(MatMultTransposeAddEqual(A,S,10,&flg));
175c4762a1bSJed Brown     if (!flg) {
1769566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mult add (T)\n",test));
177c4762a1bSJed Brown     }
178c4762a1bSJed Brown     if (testzerorows) {
179c4762a1bSJed Brown       Mat       ST,B,C,BT,BTT;
180c4762a1bSJed Brown       IS        zr;
181c4762a1bSJed Brown       Vec       x = NULL, b1 = NULL, b2 = NULL;
182c4762a1bSJed Brown       PetscInt  *idxs = NULL, nr = 0;
183c4762a1bSJed Brown 
184c4762a1bSJed Brown       if (rank == (test%size)) {
185c4762a1bSJed Brown         nr = 1;
1869566063dSJacob Faibussowitsch         PetscCall(PetscMalloc1(nr,&idxs));
187c4762a1bSJed Brown         if (test%2) {
188c4762a1bSJed Brown           idxs[0] = (2*M - 1 - test/2)%M;
189c4762a1bSJed Brown         } else {
190c4762a1bSJed Brown           idxs[0] = (test/2)%M;
191c4762a1bSJed Brown         }
192c4762a1bSJed Brown         idxs[0] = PetscMax(idxs[0],0);
193c4762a1bSJed Brown       }
1949566063dSJacob Faibussowitsch       PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr));
1959566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)zr,"ZR"));
1969566063dSJacob Faibussowitsch       PetscCall(ISViewFromOptions(zr,NULL,"-view_is"));
1979566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(A,&x,&b1));
198c4762a1bSJed Brown       if (randomize) {
1999566063dSJacob Faibussowitsch         PetscCall(VecSetRandom(x,NULL));
2009566063dSJacob Faibussowitsch         PetscCall(VecSetRandom(b1,NULL));
201c4762a1bSJed Brown       } else {
2029566063dSJacob Faibussowitsch         PetscCall(VecSet(x,11.4));
2039566063dSJacob Faibussowitsch         PetscCall(VecSet(b1,-14.2));
204c4762a1bSJed Brown       }
2059566063dSJacob Faibussowitsch       PetscCall(VecDuplicate(b1,&b2));
2069566063dSJacob Faibussowitsch       PetscCall(VecCopy(b1,b2));
2079566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)b1,"A_B1"));
2089566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)b2,"A_B2"));
209c4762a1bSJed Brown       if (size > 1 && !cong) { /* MATMPIAIJ ZeroRows and ZeroRowsColumns are buggy in this case */
2109566063dSJacob Faibussowitsch         PetscCall(VecDestroy(&b1));
211c4762a1bSJed Brown       }
212c4762a1bSJed Brown       if (ronl) {
2139566063dSJacob Faibussowitsch         PetscCall(MatZeroRowsIS(A,zr,diag,x,b1));
2149566063dSJacob Faibussowitsch         PetscCall(MatZeroRowsIS(S,zr,diag,x,b2));
215c4762a1bSJed Brown       } else {
2169566063dSJacob Faibussowitsch         PetscCall(MatZeroRowsColumnsIS(A,zr,diag,x,b1));
2179566063dSJacob Faibussowitsch         PetscCall(MatZeroRowsColumnsIS(S,zr,diag,x,b2));
2189566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&zr));
219c4762a1bSJed Brown         /* Mix zerorows and zerorowscols */
220c4762a1bSJed Brown         nr   = 0;
221c4762a1bSJed Brown         idxs = NULL;
222dd400576SPatrick Sanan         if (rank == 0) {
223c4762a1bSJed Brown           nr   = 1;
2249566063dSJacob Faibussowitsch           PetscCall(PetscMalloc1(nr,&idxs));
225c4762a1bSJed Brown           if (test%2) {
226c4762a1bSJed Brown             idxs[0] = (3*M - 2 - test/2)%M;
227c4762a1bSJed Brown           } else {
228c4762a1bSJed Brown             idxs[0] = (test/2+1)%M;
229c4762a1bSJed Brown           }
230c4762a1bSJed Brown           idxs[0] = PetscMax(idxs[0],0);
231c4762a1bSJed Brown         }
2329566063dSJacob Faibussowitsch         PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr));
2339566063dSJacob Faibussowitsch         PetscCall(PetscObjectSetName((PetscObject)zr,"ZR2"));
2349566063dSJacob Faibussowitsch         PetscCall(ISViewFromOptions(zr,NULL,"-view_is"));
2359566063dSJacob Faibussowitsch         PetscCall(MatZeroRowsIS(A,zr,diag*2.0+PETSC_SMALL,NULL,NULL));
2369566063dSJacob Faibussowitsch         PetscCall(MatZeroRowsIS(S,zr,diag*2.0+PETSC_SMALL,NULL,NULL));
237c4762a1bSJed Brown       }
2389566063dSJacob Faibussowitsch       PetscCall(ISDestroy(&zr));
239c4762a1bSJed Brown 
240c4762a1bSJed Brown       if (b1) {
241c4762a1bSJed Brown         Vec b;
242c4762a1bSJed Brown 
2439566063dSJacob Faibussowitsch         PetscCall(VecViewFromOptions(b1,NULL,"-view_b"));
2449566063dSJacob Faibussowitsch         PetscCall(VecViewFromOptions(b2,NULL,"-view_b"));
2459566063dSJacob Faibussowitsch         PetscCall(VecDuplicate(b1,&b));
2469566063dSJacob Faibussowitsch         PetscCall(VecCopy(b1,b));
2479566063dSJacob Faibussowitsch         PetscCall(VecAXPY(b,-1.0,b2));
2489566063dSJacob Faibussowitsch         PetscCall(VecNorm(b,NORM_INFINITY,&err));
249c4762a1bSJed Brown         if (err >= tol) {
2509566063dSJacob Faibussowitsch           PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error b %g\n",test,(double)err));
251c4762a1bSJed Brown         }
2529566063dSJacob Faibussowitsch         PetscCall(VecDestroy(&b));
253c4762a1bSJed Brown       }
2549566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&b1));
2559566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&b2));
2569566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&x));
2579566063dSJacob Faibussowitsch       PetscCall(MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&B));
258c4762a1bSJed Brown 
2599566063dSJacob Faibussowitsch       PetscCall(MatCreateTranspose(S,&ST));
2609566063dSJacob Faibussowitsch       PetscCall(MatComputeOperator(ST,MATDENSE,&BT));
2619566063dSJacob Faibussowitsch       PetscCall(MatTranspose(BT,MAT_INITIAL_MATRIX,&BTT));
2629566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)B,"S"));
2639566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)BTT,"STT"));
2649566063dSJacob Faibussowitsch       PetscCall(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&C));
2659566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)C,"A"));
266c4762a1bSJed Brown 
2679566063dSJacob Faibussowitsch       PetscCall(MatViewFromOptions(C,NULL,"-view_mat"));
2689566063dSJacob Faibussowitsch       PetscCall(MatViewFromOptions(B,NULL,"-view_mat"));
2699566063dSJacob Faibussowitsch       PetscCall(MatViewFromOptions(BTT,NULL,"-view_mat"));
270c4762a1bSJed Brown 
2719566063dSJacob Faibussowitsch       PetscCall(MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN));
2729566063dSJacob Faibussowitsch       PetscCall(MatNorm(C,NORM_FROBENIUS,&err));
273c4762a1bSJed Brown       if (err >= tol) {
2749566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat mult after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err));
275c4762a1bSJed Brown       }
276c4762a1bSJed Brown 
2779566063dSJacob Faibussowitsch       PetscCall(MatConvert(A,MATDENSE,MAT_REUSE_MATRIX,&C));
2789566063dSJacob Faibussowitsch       PetscCall(MatAXPY(C,-1.0,BTT,SAME_NONZERO_PATTERN));
2799566063dSJacob Faibussowitsch       PetscCall(MatNorm(C,NORM_FROBENIUS,&err));
280c4762a1bSJed Brown       if (err >= tol) {
2819566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat mult transpose after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err));
282c4762a1bSJed Brown       }
283c4762a1bSJed Brown 
2849566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&ST));
2859566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&BTT));
2869566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&BT));
2879566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&B));
2889566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C));
289c4762a1bSJed Brown     }
290c4762a1bSJed Brown     if (testdiagscale) { /* MatDiagonalScale() */
291c4762a1bSJed Brown       Vec vr,vl;
292c4762a1bSJed Brown 
2939566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(A,&vr,&vl));
294c4762a1bSJed Brown       if (randomize) {
2959566063dSJacob Faibussowitsch         PetscCall(VecSetRandom(vr,NULL));
2969566063dSJacob Faibussowitsch         PetscCall(VecSetRandom(vl,NULL));
297c4762a1bSJed Brown       } else {
2989566063dSJacob Faibussowitsch         PetscCall(VecSet(vr,test%2 ? 0.15 : 1.0/0.15));
2999566063dSJacob Faibussowitsch         PetscCall(VecSet(vl,test%2 ? -1.2 : 1.0/-1.2));
300c4762a1bSJed Brown       }
3019566063dSJacob Faibussowitsch       PetscCall(MatDiagonalScale(A,vl,vr));
3029566063dSJacob Faibussowitsch       PetscCall(MatDiagonalScale(S,vl,vr));
3039566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&vr));
3049566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&vl));
305c4762a1bSJed Brown     }
306c4762a1bSJed Brown 
307c4762a1bSJed Brown     if (testscale) { /* MatScale() */
3089566063dSJacob Faibussowitsch       PetscCall(MatScale(A,test%2 ? 1.4 : 1.0/1.4));
3099566063dSJacob Faibussowitsch       PetscCall(MatScale(S,test%2 ? 1.4 : 1.0/1.4));
310c4762a1bSJed Brown     }
311c4762a1bSJed Brown 
312c4762a1bSJed Brown     if (testshift && cong) { /* MatShift() : MATSHELL shift is broken when row/cols layout are not congruent and left/right scaling have been applied */
3139566063dSJacob Faibussowitsch       PetscCall(MatShift(A,test%2 ? -77.5 : 77.5));
3149566063dSJacob Faibussowitsch       PetscCall(MatShift(S,test%2 ? -77.5 : 77.5));
315c4762a1bSJed Brown     }
316c4762a1bSJed Brown 
317c4762a1bSJed Brown     if (testgetdiag && cong) { /* MatGetDiagonal() */
318c4762a1bSJed Brown       Vec dA,dS;
319c4762a1bSJed Brown 
3209566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(A,&dA,NULL));
3219566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(S,&dS,NULL));
3229566063dSJacob Faibussowitsch       PetscCall(MatGetDiagonal(A,dA));
3239566063dSJacob Faibussowitsch       PetscCall(MatGetDiagonal(S,dS));
3249566063dSJacob Faibussowitsch       PetscCall(VecAXPY(dA,-1.0,dS));
3259566063dSJacob Faibussowitsch       PetscCall(VecNorm(dA,NORM_INFINITY,&err));
326c4762a1bSJed Brown       if (err >= tol) {
3279566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error diag %g\n",test,(double)err));
328c4762a1bSJed Brown       }
3299566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&dA));
3309566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&dS));
331c4762a1bSJed Brown     }
332c4762a1bSJed Brown 
333c4762a1bSJed Brown     if (testdup && !test) {
334c4762a1bSJed Brown       Mat A2, S2;
335c4762a1bSJed Brown 
3369566063dSJacob Faibussowitsch       PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2));
3379566063dSJacob Faibussowitsch       PetscCall(MatDuplicate(S,MAT_COPY_VALUES,&S2));
3389566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&A));
3399566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&S));
340c4762a1bSJed Brown       A = A2;
341c4762a1bSJed Brown       S = S2;
342c4762a1bSJed Brown     }
343c4762a1bSJed Brown 
344c4762a1bSJed Brown     if (testsubmat) {
345c4762a1bSJed Brown       Mat      sA,sS,dA,dS,At,St;
346c4762a1bSJed Brown       IS       r,c;
347c4762a1bSJed Brown       PetscInt rst,ren,cst,cen;
348c4762a1bSJed Brown 
3499566063dSJacob Faibussowitsch       PetscCall(MatGetOwnershipRange(A,&rst,&ren));
3509566063dSJacob Faibussowitsch       PetscCall(MatGetOwnershipRangeColumn(A,&cst,&cen));
3519566063dSJacob Faibussowitsch       PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r));
3529566063dSJacob Faibussowitsch       PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c));
3539566063dSJacob Faibussowitsch       PetscCall(MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&sA));
3549566063dSJacob Faibussowitsch       PetscCall(MatCreateSubMatrix(S,r,c,MAT_INITIAL_MATRIX,&sS));
3559566063dSJacob Faibussowitsch       PetscCall(MatMultAddEqual(sA,sS,10,&flg));
356c4762a1bSJed Brown       if (!flg) {
3579566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix mult add\n",test));
358c4762a1bSJed Brown       }
3599566063dSJacob Faibussowitsch       PetscCall(MatMultTransposeAddEqual(sA,sS,10,&flg));
360c4762a1bSJed Brown       if (!flg) {
3619566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix mult add (T)\n",test));
362c4762a1bSJed Brown       }
3639566063dSJacob Faibussowitsch       PetscCall(MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA));
3649566063dSJacob Faibussowitsch       PetscCall(MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS));
3659566063dSJacob Faibussowitsch       PetscCall(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN));
3669566063dSJacob Faibussowitsch       PetscCall(MatNorm(dA,NORM_FROBENIUS,&err));
367c4762a1bSJed Brown       if (err >= tol) {
3689566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix %g\n",test,(double)err));
369c4762a1bSJed Brown       }
3709566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&sA));
3719566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&sS));
3729566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&dA));
3739566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&dS));
3749566063dSJacob Faibussowitsch       PetscCall(MatCreateTranspose(A,&At));
3759566063dSJacob Faibussowitsch       PetscCall(MatCreateTranspose(S,&St));
3769566063dSJacob Faibussowitsch       PetscCall(MatCreateSubMatrix(At,c,r,MAT_INITIAL_MATRIX,&sA));
3779566063dSJacob Faibussowitsch       PetscCall(MatCreateSubMatrix(St,c,r,MAT_INITIAL_MATRIX,&sS));
3789566063dSJacob Faibussowitsch       PetscCall(MatMultAddEqual(sA,sS,10,&flg));
379c4762a1bSJed Brown       if (!flg) {
3809566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix (T) mult add\n",test));
381c4762a1bSJed Brown       }
3829566063dSJacob Faibussowitsch       PetscCall(MatMultTransposeAddEqual(sA,sS,10,&flg));
383c4762a1bSJed Brown       if (!flg) {
3849566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix (T) mult add (T)\n",test));
385c4762a1bSJed Brown       }
3869566063dSJacob Faibussowitsch       PetscCall(MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA));
3879566063dSJacob Faibussowitsch       PetscCall(MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS));
3889566063dSJacob Faibussowitsch       PetscCall(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN));
3899566063dSJacob Faibussowitsch       PetscCall(MatNorm(dA,NORM_FROBENIUS,&err));
390c4762a1bSJed Brown       if (err >= tol) {
3919566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix (T) %g\n",test,(double)err));
392c4762a1bSJed Brown       }
3939566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&sA));
3949566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&sS));
3959566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&dA));
3969566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&dS));
3979566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&At));
3989566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&St));
3999566063dSJacob Faibussowitsch       PetscCall(ISDestroy(&r));
4009566063dSJacob Faibussowitsch       PetscCall(ISDestroy(&c));
401c4762a1bSJed Brown     }
402c4762a1bSJed Brown 
403c4762a1bSJed Brown     if (testaxpy) {
404c4762a1bSJed Brown       Mat          tA,tS,dA,dS;
405c4762a1bSJed Brown       MatStructure str[3] = { SAME_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN };
406c4762a1bSJed Brown 
4079566063dSJacob Faibussowitsch       PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&tA));
408c4762a1bSJed Brown       if (testaxpyd && !(test%2)) {
4099566063dSJacob Faibussowitsch         PetscCall(PetscObjectReference((PetscObject)tA));
410c4762a1bSJed Brown         tS   = tA;
411c4762a1bSJed Brown       } else {
4129566063dSJacob Faibussowitsch         PetscCall(PetscObjectReference((PetscObject)S));
413c4762a1bSJed Brown         tS   = S;
414c4762a1bSJed Brown       }
4159566063dSJacob Faibussowitsch       PetscCall(MatAXPY(A,0.5,tA,str[test%3]));
4169566063dSJacob Faibussowitsch       PetscCall(MatAXPY(S,0.5,tS,str[test%3]));
417c4762a1bSJed Brown       /* this will trigger an error the next MatMult or MatMultTranspose call for S */
4189566063dSJacob Faibussowitsch       if (testaxpyerr) PetscCall(MatScale(tA,0));
4199566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&tA));
4209566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&tS));
4219566063dSJacob Faibussowitsch       PetscCall(MatMultAddEqual(A,S,10,&flg));
422c4762a1bSJed Brown       if (!flg) {
4239566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error axpy mult add\n",test));
424c4762a1bSJed Brown       }
4259566063dSJacob Faibussowitsch       PetscCall(MatMultTransposeAddEqual(A,S,10,&flg));
426c4762a1bSJed Brown       if (!flg) {
4279566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error axpy mult add (T)\n",test));
428c4762a1bSJed Brown       }
4299566063dSJacob Faibussowitsch       PetscCall(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&dA));
4309566063dSJacob Faibussowitsch       PetscCall(MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&dS));
4319566063dSJacob Faibussowitsch       PetscCall(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN));
4329566063dSJacob Faibussowitsch       PetscCall(MatNorm(dA,NORM_FROBENIUS,&err));
433c4762a1bSJed Brown       if (err >= tol) {
4349566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix %g\n",test,(double)err));
435c4762a1bSJed Brown       }
4369566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&dA));
4379566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&dS));
438c4762a1bSJed Brown     }
439c4762a1bSJed Brown 
440c4762a1bSJed Brown     if (testreset && (ntest == 1 || test == ntest-2)) {
441c4762a1bSJed Brown       /* reset MATSHELL */
4429566063dSJacob Faibussowitsch       PetscCall(MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY));
4439566063dSJacob Faibussowitsch       PetscCall(MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY));
444c4762a1bSJed Brown       /* reset A */
4459566063dSJacob Faibussowitsch       PetscCall(MatCopy(user->B,A,DIFFERENT_NONZERO_PATTERN));
446c4762a1bSJed Brown     }
447c4762a1bSJed Brown   }
448c4762a1bSJed Brown 
4499566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&A));
4509566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&S));
4519566063dSJacob Faibussowitsch   PetscCall(PetscFree(user));
4529566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
453b122ec5aSJacob Faibussowitsch   return 0;
454c4762a1bSJed Brown }
455c4762a1bSJed Brown 
456c4762a1bSJed Brown /*TEST
457c4762a1bSJed Brown 
458c4762a1bSJed Brown    testset:
459c4762a1bSJed Brown      suffix: rect
460c4762a1bSJed Brown      requires: !single
461c4762a1bSJed Brown      output_file: output/ex221_1.out
462c4762a1bSJed Brown      nsize: {{1 3}}
463c4762a1bSJed Brown      args: -loop 3 -keep {{0 1}} -M {{12 19}} -N {{19 12}} -submat {{0 1}} -test_axpy_different {{0 1}}
464c4762a1bSJed Brown 
465c4762a1bSJed Brown    testset:
466c4762a1bSJed Brown      suffix: square
467c4762a1bSJed Brown      requires: !single
468c4762a1bSJed Brown      output_file: output/ex221_1.out
469c4762a1bSJed Brown      nsize: {{1 3}}
470c4762a1bSJed Brown      args: -M 21 -N 21 -loop 4 -rows_only {{0 1}} -keep {{0 1}} -submat {{0 1}} -test_axpy_different {{0 1}}
471c4762a1bSJed Brown TEST*/
472