xref: /petsc/src/mat/tests/ex221.c (revision 327415f76d85372a4417cf1aaa14db707d4d6c04)
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));
4708401ef6SPierre 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 
79*327415f7SBarry Smith   PetscFunctionBeginUser;
809566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc,&args,(char*)0,help));
819566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
829566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
839566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL));
849566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL));
859566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-ml",&m,NULL));
869566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-nl",&n,NULL));
879566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-square_nc",&nc,NULL));
889566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-rows_only",&ronl,NULL));
899566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-randomize",&randomize,NULL));
909566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-submat",&submat,NULL));
919566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_zerorows",&testzerorows,NULL));
929566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_diagscale",&testdiagscale,NULL));
939566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_getdiag",&testgetdiag,NULL));
949566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_shift",&testshift,NULL));
959566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_scale",&testscale,NULL));
969566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_dup",&testdup,NULL));
979566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_reset",&testreset,NULL));
989566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_submat",&testsubmat,NULL));
999566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_axpy",&testaxpy,NULL));
1009566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_axpy_different",&testaxpyd,NULL));
1019566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_axpy_error",&testaxpyerr,NULL));
1029566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-loop",&ntest,NULL));
1039566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL));
1049566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-diag",&diag,NULL));
1059566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-keep",&keep,NULL));
106c4762a1bSJed Brown   /* This tests square matrices with different row/col layout */
107c4762a1bSJed Brown   if (nc && size > 1) {
108c4762a1bSJed Brown     M = PetscMax(PetscMax(N,M),1);
109c4762a1bSJed Brown     N = M;
110c4762a1bSJed Brown     m = n = 0;
111c4762a1bSJed Brown     if (rank == 0) { m = M-1; n = 1; }
112c4762a1bSJed Brown     else if (rank == 1) { m = 1; n = N-1; }
113c4762a1bSJed Brown   }
1149566063dSJacob Faibussowitsch   PetscCall(MatCreateDense(PETSC_COMM_WORLD,m,n,M,N,NULL,&A));
1159566063dSJacob Faibussowitsch   PetscCall(MatGetLocalSize(A,&m,&n));
1169566063dSJacob Faibussowitsch   PetscCall(MatGetSize(A,&M,&N));
1179566063dSJacob Faibussowitsch   PetscCall(MatGetOwnershipRange(A,&s1,NULL));
118c4762a1bSJed Brown   s2   = 1;
119c4762a1bSJed Brown   while (s2 < M) s2 *= 10;
1209566063dSJacob Faibussowitsch   PetscCall(MatDenseGetArray(A,&data));
121c4762a1bSJed Brown   for (j = 0; j < N; j++) {
122c4762a1bSJed Brown     for (i = 0; i < m; i++) {
123c4762a1bSJed Brown       data[j*m + i] = s2*j + i + s1 + 1;
124c4762a1bSJed Brown     }
125c4762a1bSJed Brown   }
1269566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
1279566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
128c4762a1bSJed Brown 
129c4762a1bSJed Brown   if (submat) {
130c4762a1bSJed Brown     Mat      A2;
131c4762a1bSJed Brown     IS       r,c;
132c4762a1bSJed Brown     PetscInt rst,ren,cst,cen;
133c4762a1bSJed Brown 
1349566063dSJacob Faibussowitsch     PetscCall(MatGetOwnershipRange(A,&rst,&ren));
1359566063dSJacob Faibussowitsch     PetscCall(MatGetOwnershipRangeColumn(A,&cst,&cen));
1369566063dSJacob Faibussowitsch     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r));
1379566063dSJacob Faibussowitsch     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c));
1389566063dSJacob Faibussowitsch     PetscCall(MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&A2));
1399566063dSJacob Faibussowitsch     PetscCall(ISDestroy(&r));
1409566063dSJacob Faibussowitsch     PetscCall(ISDestroy(&c));
1419566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&A));
142c4762a1bSJed Brown     A = A2;
143c4762a1bSJed Brown   }
144c4762a1bSJed Brown 
1459566063dSJacob Faibussowitsch   PetscCall(MatGetSize(A,&M,&N));
1469566063dSJacob Faibussowitsch   PetscCall(MatGetLocalSize(A,&m,&n));
1479566063dSJacob Faibussowitsch   PetscCall(MatHasCongruentLayouts(A,&cong));
148c4762a1bSJed Brown 
1499566063dSJacob Faibussowitsch   PetscCall(MatConvert(A,MATAIJ,MAT_INPLACE_MATRIX,&A));
1509566063dSJacob Faibussowitsch   PetscCall(MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,keep));
1519566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)A,"initial"));
1529566063dSJacob Faibussowitsch   PetscCall(MatViewFromOptions(A,NULL,"-view_mat"));
153c4762a1bSJed Brown 
1549566063dSJacob Faibussowitsch   PetscCall(PetscNew(&user));
1559566063dSJacob Faibussowitsch   PetscCall(MatCreateShell(PETSC_COMM_WORLD,m,n,M,N,user,&S));
1569566063dSJacob Faibussowitsch   PetscCall(MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User));
1579566063dSJacob Faibussowitsch   PetscCall(MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User));
1581baa6e33SBarry Smith   if (cong) PetscCall(MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User));
1599566063dSJacob Faibussowitsch   PetscCall(MatShellSetOperation(S,MATOP_COPY,(void (*)(void))MatCopy_User));
1609566063dSJacob Faibussowitsch   PetscCall(MatShellSetOperation(S,MATOP_DESTROY,(void (*)(void))MatDestroy_User));
1619566063dSJacob Faibussowitsch   PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&user->B));
162c4762a1bSJed Brown 
163c4762a1bSJed Brown   /* Square and rows only scaling */
164c4762a1bSJed Brown   ronl = cong ? ronl : PETSC_TRUE;
165c4762a1bSJed Brown 
166c4762a1bSJed Brown   for (test = 0; test < ntest; test++) {
167c4762a1bSJed Brown     PetscReal err;
168c4762a1bSJed Brown 
1699566063dSJacob Faibussowitsch     PetscCall(MatMultAddEqual(A,S,10,&flg));
170c4762a1bSJed Brown     if (!flg) {
1719566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mult add\n",test));
172c4762a1bSJed Brown     }
1739566063dSJacob Faibussowitsch     PetscCall(MatMultTransposeAddEqual(A,S,10,&flg));
174c4762a1bSJed Brown     if (!flg) {
1759566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mult add (T)\n",test));
176c4762a1bSJed Brown     }
177c4762a1bSJed Brown     if (testzerorows) {
178c4762a1bSJed Brown       Mat       ST,B,C,BT,BTT;
179c4762a1bSJed Brown       IS        zr;
180c4762a1bSJed Brown       Vec       x = NULL, b1 = NULL, b2 = NULL;
181c4762a1bSJed Brown       PetscInt  *idxs = NULL, nr = 0;
182c4762a1bSJed Brown 
183c4762a1bSJed Brown       if (rank == (test%size)) {
184c4762a1bSJed Brown         nr = 1;
1859566063dSJacob Faibussowitsch         PetscCall(PetscMalloc1(nr,&idxs));
186c4762a1bSJed Brown         if (test%2) {
187c4762a1bSJed Brown           idxs[0] = (2*M - 1 - test/2)%M;
188c4762a1bSJed Brown         } else {
189c4762a1bSJed Brown           idxs[0] = (test/2)%M;
190c4762a1bSJed Brown         }
191c4762a1bSJed Brown         idxs[0] = PetscMax(idxs[0],0);
192c4762a1bSJed Brown       }
1939566063dSJacob Faibussowitsch       PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr));
1949566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)zr,"ZR"));
1959566063dSJacob Faibussowitsch       PetscCall(ISViewFromOptions(zr,NULL,"-view_is"));
1969566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(A,&x,&b1));
197c4762a1bSJed Brown       if (randomize) {
1989566063dSJacob Faibussowitsch         PetscCall(VecSetRandom(x,NULL));
1999566063dSJacob Faibussowitsch         PetscCall(VecSetRandom(b1,NULL));
200c4762a1bSJed Brown       } else {
2019566063dSJacob Faibussowitsch         PetscCall(VecSet(x,11.4));
2029566063dSJacob Faibussowitsch         PetscCall(VecSet(b1,-14.2));
203c4762a1bSJed Brown       }
2049566063dSJacob Faibussowitsch       PetscCall(VecDuplicate(b1,&b2));
2059566063dSJacob Faibussowitsch       PetscCall(VecCopy(b1,b2));
2069566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)b1,"A_B1"));
2079566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)b2,"A_B2"));
208c4762a1bSJed Brown       if (size > 1 && !cong) { /* MATMPIAIJ ZeroRows and ZeroRowsColumns are buggy in this case */
2099566063dSJacob Faibussowitsch         PetscCall(VecDestroy(&b1));
210c4762a1bSJed Brown       }
211c4762a1bSJed Brown       if (ronl) {
2129566063dSJacob Faibussowitsch         PetscCall(MatZeroRowsIS(A,zr,diag,x,b1));
2139566063dSJacob Faibussowitsch         PetscCall(MatZeroRowsIS(S,zr,diag,x,b2));
214c4762a1bSJed Brown       } else {
2159566063dSJacob Faibussowitsch         PetscCall(MatZeroRowsColumnsIS(A,zr,diag,x,b1));
2169566063dSJacob Faibussowitsch         PetscCall(MatZeroRowsColumnsIS(S,zr,diag,x,b2));
2179566063dSJacob Faibussowitsch         PetscCall(ISDestroy(&zr));
218c4762a1bSJed Brown         /* Mix zerorows and zerorowscols */
219c4762a1bSJed Brown         nr   = 0;
220c4762a1bSJed Brown         idxs = NULL;
221dd400576SPatrick Sanan         if (rank == 0) {
222c4762a1bSJed Brown           nr   = 1;
2239566063dSJacob Faibussowitsch           PetscCall(PetscMalloc1(nr,&idxs));
224c4762a1bSJed Brown           if (test%2) {
225c4762a1bSJed Brown             idxs[0] = (3*M - 2 - test/2)%M;
226c4762a1bSJed Brown           } else {
227c4762a1bSJed Brown             idxs[0] = (test/2+1)%M;
228c4762a1bSJed Brown           }
229c4762a1bSJed Brown           idxs[0] = PetscMax(idxs[0],0);
230c4762a1bSJed Brown         }
2319566063dSJacob Faibussowitsch         PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr));
2329566063dSJacob Faibussowitsch         PetscCall(PetscObjectSetName((PetscObject)zr,"ZR2"));
2339566063dSJacob Faibussowitsch         PetscCall(ISViewFromOptions(zr,NULL,"-view_is"));
2349566063dSJacob Faibussowitsch         PetscCall(MatZeroRowsIS(A,zr,diag*2.0+PETSC_SMALL,NULL,NULL));
2359566063dSJacob Faibussowitsch         PetscCall(MatZeroRowsIS(S,zr,diag*2.0+PETSC_SMALL,NULL,NULL));
236c4762a1bSJed Brown       }
2379566063dSJacob Faibussowitsch       PetscCall(ISDestroy(&zr));
238c4762a1bSJed Brown 
239c4762a1bSJed Brown       if (b1) {
240c4762a1bSJed Brown         Vec b;
241c4762a1bSJed Brown 
2429566063dSJacob Faibussowitsch         PetscCall(VecViewFromOptions(b1,NULL,"-view_b"));
2439566063dSJacob Faibussowitsch         PetscCall(VecViewFromOptions(b2,NULL,"-view_b"));
2449566063dSJacob Faibussowitsch         PetscCall(VecDuplicate(b1,&b));
2459566063dSJacob Faibussowitsch         PetscCall(VecCopy(b1,b));
2469566063dSJacob Faibussowitsch         PetscCall(VecAXPY(b,-1.0,b2));
2479566063dSJacob Faibussowitsch         PetscCall(VecNorm(b,NORM_INFINITY,&err));
248c4762a1bSJed Brown         if (err >= tol) {
2499566063dSJacob Faibussowitsch           PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error b %g\n",test,(double)err));
250c4762a1bSJed Brown         }
2519566063dSJacob Faibussowitsch         PetscCall(VecDestroy(&b));
252c4762a1bSJed Brown       }
2539566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&b1));
2549566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&b2));
2559566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&x));
2569566063dSJacob Faibussowitsch       PetscCall(MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&B));
257c4762a1bSJed Brown 
2589566063dSJacob Faibussowitsch       PetscCall(MatCreateTranspose(S,&ST));
2599566063dSJacob Faibussowitsch       PetscCall(MatComputeOperator(ST,MATDENSE,&BT));
2609566063dSJacob Faibussowitsch       PetscCall(MatTranspose(BT,MAT_INITIAL_MATRIX,&BTT));
2619566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)B,"S"));
2629566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)BTT,"STT"));
2639566063dSJacob Faibussowitsch       PetscCall(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&C));
2649566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)C,"A"));
265c4762a1bSJed Brown 
2669566063dSJacob Faibussowitsch       PetscCall(MatViewFromOptions(C,NULL,"-view_mat"));
2679566063dSJacob Faibussowitsch       PetscCall(MatViewFromOptions(B,NULL,"-view_mat"));
2689566063dSJacob Faibussowitsch       PetscCall(MatViewFromOptions(BTT,NULL,"-view_mat"));
269c4762a1bSJed Brown 
2709566063dSJacob Faibussowitsch       PetscCall(MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN));
2719566063dSJacob Faibussowitsch       PetscCall(MatNorm(C,NORM_FROBENIUS,&err));
272c4762a1bSJed Brown       if (err >= tol) {
2739566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat mult after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err));
274c4762a1bSJed Brown       }
275c4762a1bSJed Brown 
2769566063dSJacob Faibussowitsch       PetscCall(MatConvert(A,MATDENSE,MAT_REUSE_MATRIX,&C));
2779566063dSJacob Faibussowitsch       PetscCall(MatAXPY(C,-1.0,BTT,SAME_NONZERO_PATTERN));
2789566063dSJacob Faibussowitsch       PetscCall(MatNorm(C,NORM_FROBENIUS,&err));
279c4762a1bSJed Brown       if (err >= tol) {
2809566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat mult transpose after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err));
281c4762a1bSJed Brown       }
282c4762a1bSJed Brown 
2839566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&ST));
2849566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&BTT));
2859566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&BT));
2869566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&B));
2879566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C));
288c4762a1bSJed Brown     }
289c4762a1bSJed Brown     if (testdiagscale) { /* MatDiagonalScale() */
290c4762a1bSJed Brown       Vec vr,vl;
291c4762a1bSJed Brown 
2929566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(A,&vr,&vl));
293c4762a1bSJed Brown       if (randomize) {
2949566063dSJacob Faibussowitsch         PetscCall(VecSetRandom(vr,NULL));
2959566063dSJacob Faibussowitsch         PetscCall(VecSetRandom(vl,NULL));
296c4762a1bSJed Brown       } else {
2979566063dSJacob Faibussowitsch         PetscCall(VecSet(vr,test%2 ? 0.15 : 1.0/0.15));
2989566063dSJacob Faibussowitsch         PetscCall(VecSet(vl,test%2 ? -1.2 : 1.0/-1.2));
299c4762a1bSJed Brown       }
3009566063dSJacob Faibussowitsch       PetscCall(MatDiagonalScale(A,vl,vr));
3019566063dSJacob Faibussowitsch       PetscCall(MatDiagonalScale(S,vl,vr));
3029566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&vr));
3039566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&vl));
304c4762a1bSJed Brown     }
305c4762a1bSJed Brown 
306c4762a1bSJed Brown     if (testscale) { /* MatScale() */
3079566063dSJacob Faibussowitsch       PetscCall(MatScale(A,test%2 ? 1.4 : 1.0/1.4));
3089566063dSJacob Faibussowitsch       PetscCall(MatScale(S,test%2 ? 1.4 : 1.0/1.4));
309c4762a1bSJed Brown     }
310c4762a1bSJed Brown 
311c4762a1bSJed Brown     if (testshift && cong) { /* MatShift() : MATSHELL shift is broken when row/cols layout are not congruent and left/right scaling have been applied */
3129566063dSJacob Faibussowitsch       PetscCall(MatShift(A,test%2 ? -77.5 : 77.5));
3139566063dSJacob Faibussowitsch       PetscCall(MatShift(S,test%2 ? -77.5 : 77.5));
314c4762a1bSJed Brown     }
315c4762a1bSJed Brown 
316c4762a1bSJed Brown     if (testgetdiag && cong) { /* MatGetDiagonal() */
317c4762a1bSJed Brown       Vec dA,dS;
318c4762a1bSJed Brown 
3199566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(A,&dA,NULL));
3209566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(S,&dS,NULL));
3219566063dSJacob Faibussowitsch       PetscCall(MatGetDiagonal(A,dA));
3229566063dSJacob Faibussowitsch       PetscCall(MatGetDiagonal(S,dS));
3239566063dSJacob Faibussowitsch       PetscCall(VecAXPY(dA,-1.0,dS));
3249566063dSJacob Faibussowitsch       PetscCall(VecNorm(dA,NORM_INFINITY,&err));
325c4762a1bSJed Brown       if (err >= tol) {
3269566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error diag %g\n",test,(double)err));
327c4762a1bSJed Brown       }
3289566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&dA));
3299566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&dS));
330c4762a1bSJed Brown     }
331c4762a1bSJed Brown 
332c4762a1bSJed Brown     if (testdup && !test) {
333c4762a1bSJed Brown       Mat A2, S2;
334c4762a1bSJed Brown 
3359566063dSJacob Faibussowitsch       PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2));
3369566063dSJacob Faibussowitsch       PetscCall(MatDuplicate(S,MAT_COPY_VALUES,&S2));
3379566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&A));
3389566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&S));
339c4762a1bSJed Brown       A = A2;
340c4762a1bSJed Brown       S = S2;
341c4762a1bSJed Brown     }
342c4762a1bSJed Brown 
343c4762a1bSJed Brown     if (testsubmat) {
344c4762a1bSJed Brown       Mat      sA,sS,dA,dS,At,St;
345c4762a1bSJed Brown       IS       r,c;
346c4762a1bSJed Brown       PetscInt rst,ren,cst,cen;
347c4762a1bSJed Brown 
3489566063dSJacob Faibussowitsch       PetscCall(MatGetOwnershipRange(A,&rst,&ren));
3499566063dSJacob Faibussowitsch       PetscCall(MatGetOwnershipRangeColumn(A,&cst,&cen));
3509566063dSJacob Faibussowitsch       PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r));
3519566063dSJacob Faibussowitsch       PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c));
3529566063dSJacob Faibussowitsch       PetscCall(MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&sA));
3539566063dSJacob Faibussowitsch       PetscCall(MatCreateSubMatrix(S,r,c,MAT_INITIAL_MATRIX,&sS));
3549566063dSJacob Faibussowitsch       PetscCall(MatMultAddEqual(sA,sS,10,&flg));
355c4762a1bSJed Brown       if (!flg) {
3569566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix mult add\n",test));
357c4762a1bSJed Brown       }
3589566063dSJacob Faibussowitsch       PetscCall(MatMultTransposeAddEqual(sA,sS,10,&flg));
359c4762a1bSJed Brown       if (!flg) {
3609566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix mult add (T)\n",test));
361c4762a1bSJed Brown       }
3629566063dSJacob Faibussowitsch       PetscCall(MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA));
3639566063dSJacob Faibussowitsch       PetscCall(MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS));
3649566063dSJacob Faibussowitsch       PetscCall(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN));
3659566063dSJacob Faibussowitsch       PetscCall(MatNorm(dA,NORM_FROBENIUS,&err));
366c4762a1bSJed Brown       if (err >= tol) {
3679566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix %g\n",test,(double)err));
368c4762a1bSJed Brown       }
3699566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&sA));
3709566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&sS));
3719566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&dA));
3729566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&dS));
3739566063dSJacob Faibussowitsch       PetscCall(MatCreateTranspose(A,&At));
3749566063dSJacob Faibussowitsch       PetscCall(MatCreateTranspose(S,&St));
3759566063dSJacob Faibussowitsch       PetscCall(MatCreateSubMatrix(At,c,r,MAT_INITIAL_MATRIX,&sA));
3769566063dSJacob Faibussowitsch       PetscCall(MatCreateSubMatrix(St,c,r,MAT_INITIAL_MATRIX,&sS));
3779566063dSJacob Faibussowitsch       PetscCall(MatMultAddEqual(sA,sS,10,&flg));
378c4762a1bSJed Brown       if (!flg) {
3799566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix (T) mult add\n",test));
380c4762a1bSJed Brown       }
3819566063dSJacob Faibussowitsch       PetscCall(MatMultTransposeAddEqual(sA,sS,10,&flg));
382c4762a1bSJed Brown       if (!flg) {
3839566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix (T) mult add (T)\n",test));
384c4762a1bSJed Brown       }
3859566063dSJacob Faibussowitsch       PetscCall(MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA));
3869566063dSJacob Faibussowitsch       PetscCall(MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS));
3879566063dSJacob Faibussowitsch       PetscCall(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN));
3889566063dSJacob Faibussowitsch       PetscCall(MatNorm(dA,NORM_FROBENIUS,&err));
389c4762a1bSJed Brown       if (err >= tol) {
3909566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix (T) %g\n",test,(double)err));
391c4762a1bSJed Brown       }
3929566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&sA));
3939566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&sS));
3949566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&dA));
3959566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&dS));
3969566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&At));
3979566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&St));
3989566063dSJacob Faibussowitsch       PetscCall(ISDestroy(&r));
3999566063dSJacob Faibussowitsch       PetscCall(ISDestroy(&c));
400c4762a1bSJed Brown     }
401c4762a1bSJed Brown 
402c4762a1bSJed Brown     if (testaxpy) {
403c4762a1bSJed Brown       Mat          tA,tS,dA,dS;
404c4762a1bSJed Brown       MatStructure str[3] = { SAME_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN };
405c4762a1bSJed Brown 
4069566063dSJacob Faibussowitsch       PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&tA));
407c4762a1bSJed Brown       if (testaxpyd && !(test%2)) {
4089566063dSJacob Faibussowitsch         PetscCall(PetscObjectReference((PetscObject)tA));
409c4762a1bSJed Brown         tS   = tA;
410c4762a1bSJed Brown       } else {
4119566063dSJacob Faibussowitsch         PetscCall(PetscObjectReference((PetscObject)S));
412c4762a1bSJed Brown         tS   = S;
413c4762a1bSJed Brown       }
4149566063dSJacob Faibussowitsch       PetscCall(MatAXPY(A,0.5,tA,str[test%3]));
4159566063dSJacob Faibussowitsch       PetscCall(MatAXPY(S,0.5,tS,str[test%3]));
416c4762a1bSJed Brown       /* this will trigger an error the next MatMult or MatMultTranspose call for S */
4179566063dSJacob Faibussowitsch       if (testaxpyerr) PetscCall(MatScale(tA,0));
4189566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&tA));
4199566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&tS));
4209566063dSJacob Faibussowitsch       PetscCall(MatMultAddEqual(A,S,10,&flg));
421c4762a1bSJed Brown       if (!flg) {
4229566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error axpy mult add\n",test));
423c4762a1bSJed Brown       }
4249566063dSJacob Faibussowitsch       PetscCall(MatMultTransposeAddEqual(A,S,10,&flg));
425c4762a1bSJed Brown       if (!flg) {
4269566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error axpy mult add (T)\n",test));
427c4762a1bSJed Brown       }
4289566063dSJacob Faibussowitsch       PetscCall(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&dA));
4299566063dSJacob Faibussowitsch       PetscCall(MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&dS));
4309566063dSJacob Faibussowitsch       PetscCall(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN));
4319566063dSJacob Faibussowitsch       PetscCall(MatNorm(dA,NORM_FROBENIUS,&err));
432c4762a1bSJed Brown       if (err >= tol) {
4339566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix %g\n",test,(double)err));
434c4762a1bSJed Brown       }
4359566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&dA));
4369566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&dS));
437c4762a1bSJed Brown     }
438c4762a1bSJed Brown 
439c4762a1bSJed Brown     if (testreset && (ntest == 1 || test == ntest-2)) {
440c4762a1bSJed Brown       /* reset MATSHELL */
4419566063dSJacob Faibussowitsch       PetscCall(MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY));
4429566063dSJacob Faibussowitsch       PetscCall(MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY));
443c4762a1bSJed Brown       /* reset A */
4449566063dSJacob Faibussowitsch       PetscCall(MatCopy(user->B,A,DIFFERENT_NONZERO_PATTERN));
445c4762a1bSJed Brown     }
446c4762a1bSJed Brown   }
447c4762a1bSJed Brown 
4489566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&A));
4499566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&S));
4509566063dSJacob Faibussowitsch   PetscCall(PetscFree(user));
4519566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
452b122ec5aSJacob Faibussowitsch   return 0;
453c4762a1bSJed Brown }
454c4762a1bSJed Brown 
455c4762a1bSJed Brown /*TEST
456c4762a1bSJed Brown 
457c4762a1bSJed Brown    testset:
458c4762a1bSJed Brown      suffix: rect
459c4762a1bSJed Brown      requires: !single
460c4762a1bSJed Brown      output_file: output/ex221_1.out
461c4762a1bSJed Brown      nsize: {{1 3}}
462c4762a1bSJed Brown      args: -loop 3 -keep {{0 1}} -M {{12 19}} -N {{19 12}} -submat {{0 1}} -test_axpy_different {{0 1}}
463c4762a1bSJed Brown 
464c4762a1bSJed Brown    testset:
465c4762a1bSJed Brown      suffix: square
466c4762a1bSJed Brown      requires: !single
467c4762a1bSJed Brown      output_file: output/ex221_1.out
468c4762a1bSJed Brown      nsize: {{1 3}}
469c4762a1bSJed Brown      args: -M 21 -N 21 -loop 4 -rows_only {{0 1}} -keep {{0 1}} -submat {{0 1}} -test_axpy_different {{0 1}}
470c4762a1bSJed Brown TEST*/
471