xref: /petsc/src/mat/tests/ex62.c (revision dcfd994d081f5f6b19844c8e8bc67e93edfc619a)
1c20d7725SJed Brown 
2c20d7725SJed Brown static char help[] = "Test Matrix products for AIJ matrices\n\
3c20d7725SJed Brown Input arguments are:\n\
4c20d7725SJed Brown   -fA <input_file> -fB <input_file> -fC <input_file>: file to load\n\n";
5c20d7725SJed Brown /* Example of usage:
6c20d7725SJed Brown    ./ex62 -fA <A_binary> -fB <B_binary>
7c20d7725SJed Brown    mpiexec -n 3 ./ex62 -fA medium -fB medium
8c20d7725SJed Brown */
9c20d7725SJed Brown 
10c20d7725SJed Brown #include <petscmat.h>
11c20d7725SJed Brown 
12c20d7725SJed Brown /*
13c20d7725SJed Brown      B = A - B
14c20d7725SJed Brown      norm = norm(B)
15c20d7725SJed Brown */
169371c9d4SSatish Balay PetscErrorCode MatNormDifference(Mat A, Mat B, PetscReal *norm) {
17c20d7725SJed Brown   PetscFunctionBegin;
189566063dSJacob Faibussowitsch   PetscCall(MatAXPY(B, -1.0, A, DIFFERENT_NONZERO_PATTERN));
199566063dSJacob Faibussowitsch   PetscCall(MatNorm(B, NORM_FROBENIUS, norm));
20c20d7725SJed Brown   PetscFunctionReturn(0);
21c20d7725SJed Brown }
22c20d7725SJed Brown 
239371c9d4SSatish Balay int main(int argc, char **args) {
24c20d7725SJed Brown   Mat          A, A_save, B, C, P, C1, R;
25c20d7725SJed Brown   PetscViewer  viewer;
26c20d7725SJed Brown   PetscMPIInt  size, rank;
2720b3374bSStefano Zampini   PetscInt     i, j, *idxn, PM, PN = PETSC_DECIDE, rstart, rend;
28c20d7725SJed Brown   PetscReal    norm;
29c20d7725SJed Brown   PetscRandom  rdm;
3020b3374bSStefano Zampini   char         file[2][PETSC_MAX_PATH_LEN] = {"", ""};
31c20d7725SJed Brown   PetscScalar *a, rval, alpha;
32c20d7725SJed Brown   PetscBool    Test_MatMatMult = PETSC_TRUE, Test_MatTrMat = PETSC_TRUE, Test_MatMatTr = PETSC_TRUE;
3320b3374bSStefano Zampini   PetscBool    Test_MatPtAP = PETSC_TRUE, Test_MatRARt = PETSC_TRUE, flg, seqaij, flgA, flgB;
34c20d7725SJed Brown   MatInfo      info;
3520b3374bSStefano Zampini   PetscInt     nzp = 5; /* num of nonzeros in each row of P */
36c20d7725SJed Brown   MatType      mattype;
3720b3374bSStefano Zampini   const char  *deft = MATAIJ;
3820b3374bSStefano Zampini   char         A_mattype[256], B_mattype[256];
3920b3374bSStefano Zampini   PetscInt     mcheck = 10;
40c20d7725SJed Brown 
41327415f7SBarry Smith   PetscFunctionBeginUser;
429566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
439566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
449566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
45c20d7725SJed Brown 
46c20d7725SJed Brown   /*  Load the matrices A_save and B */
47d0609cedSBarry Smith   PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "");
489566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-test_rart", "Test MatRARt", "", Test_MatRARt, &Test_MatRARt, NULL));
499566063dSJacob Faibussowitsch   PetscCall(PetscOptionsInt("-PN", "Number of columns of P", "", PN, &PN, NULL));
509566063dSJacob Faibussowitsch   PetscCall(PetscOptionsInt("-mcheck", "Number of matmult checks", "", mcheck, &mcheck, NULL));
519566063dSJacob Faibussowitsch   PetscCall(PetscOptionsString("-fA", "Path for matrix A", "", file[0], file[0], sizeof(file[0]), &flg));
5228b400f6SJacob Faibussowitsch   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Must indicate a file name for matrix A with the -fA option.");
539566063dSJacob Faibussowitsch   PetscCall(PetscOptionsString("-fB", "Path for matrix B", "", file[1], file[1], sizeof(file[1]), &flg));
549566063dSJacob Faibussowitsch   PetscCall(PetscOptionsFList("-A_mat_type", "Matrix type", "MatSetType", MatList, deft, A_mattype, 256, &flgA));
559566063dSJacob Faibussowitsch   PetscCall(PetscOptionsFList("-B_mat_type", "Matrix type", "MatSetType", MatList, deft, B_mattype, 256, &flgB));
56d0609cedSBarry Smith   PetscOptionsEnd();
57c20d7725SJed Brown 
589566063dSJacob Faibussowitsch   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[0], FILE_MODE_READ, &viewer));
599566063dSJacob Faibussowitsch   PetscCall(MatCreate(PETSC_COMM_WORLD, &A_save));
609566063dSJacob Faibussowitsch   PetscCall(MatLoad(A_save, viewer));
619566063dSJacob Faibussowitsch   PetscCall(PetscViewerDestroy(&viewer));
62c20d7725SJed Brown 
6320b3374bSStefano Zampini   if (flg) {
649566063dSJacob Faibussowitsch     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[1], FILE_MODE_READ, &viewer));
659566063dSJacob Faibussowitsch     PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
669566063dSJacob Faibussowitsch     PetscCall(MatLoad(B, viewer));
679566063dSJacob Faibussowitsch     PetscCall(PetscViewerDestroy(&viewer));
6820b3374bSStefano Zampini   } else {
699566063dSJacob Faibussowitsch     PetscCall(PetscObjectReference((PetscObject)A_save));
7020b3374bSStefano Zampini     B = A_save;
7120b3374bSStefano Zampini   }
7220b3374bSStefano Zampini 
739371c9d4SSatish Balay   if (flgA) { PetscCall(MatConvert(A_save, A_mattype, MAT_INPLACE_MATRIX, &A_save)); }
749371c9d4SSatish Balay   if (flgB) { PetscCall(MatConvert(B, B_mattype, MAT_INPLACE_MATRIX, &B)); }
759566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(A_save));
769566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(B));
77c20d7725SJed Brown 
789566063dSJacob Faibussowitsch   PetscCall(MatGetType(B, &mattype));
79c20d7725SJed Brown 
809566063dSJacob Faibussowitsch   PetscCall(PetscMalloc(nzp * (sizeof(PetscInt) + sizeof(PetscScalar)), &idxn));
81c20d7725SJed Brown   a = (PetscScalar *)(idxn + nzp);
82c20d7725SJed Brown 
839566063dSJacob Faibussowitsch   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
849566063dSJacob Faibussowitsch   PetscCall(PetscRandomSetFromOptions(rdm));
85c20d7725SJed Brown 
86c20d7725SJed Brown   /* 1) MatMatMult() */
87c20d7725SJed Brown   /* ----------------*/
88c20d7725SJed Brown   if (Test_MatMatMult) {
899566063dSJacob Faibussowitsch     PetscCall(MatDuplicate(A_save, MAT_COPY_VALUES, &A));
90c20d7725SJed Brown 
91c20d7725SJed Brown     /* (1.1) Test developer API */
929566063dSJacob Faibussowitsch     PetscCall(MatProductCreate(A, B, NULL, &C));
939566063dSJacob Faibussowitsch     PetscCall(MatSetOptionsPrefix(C, "AB_"));
949566063dSJacob Faibussowitsch     PetscCall(MatProductSetType(C, MATPRODUCT_AB));
959566063dSJacob Faibussowitsch     PetscCall(MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT));
969566063dSJacob Faibussowitsch     PetscCall(MatProductSetFill(C, PETSC_DEFAULT));
979566063dSJacob Faibussowitsch     PetscCall(MatProductSetFromOptions(C));
98910fa13fSStefano Zampini     /* we can inquire about MATOP_PRODUCTSYMBOLIC even if the destination matrix type has not been set yet */
999566063dSJacob Faibussowitsch     PetscCall(MatHasOperation(C, MATOP_PRODUCTSYMBOLIC, &flg));
1009566063dSJacob Faibussowitsch     PetscCall(MatProductSymbolic(C));
1019566063dSJacob Faibussowitsch     PetscCall(MatProductNumeric(C));
1029566063dSJacob Faibussowitsch     PetscCall(MatMatMultEqual(A, B, C, mcheck, &flg));
10328b400f6SJacob Faibussowitsch     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in C=A*B");
104c20d7725SJed Brown 
105c20d7725SJed Brown     /* Test reuse symbolic C */
106c20d7725SJed Brown     alpha = 0.9;
1079566063dSJacob Faibussowitsch     PetscCall(MatScale(A, alpha));
1089566063dSJacob Faibussowitsch     PetscCall(MatProductNumeric(C));
109c20d7725SJed Brown 
1109566063dSJacob Faibussowitsch     PetscCall(MatMatMultEqual(A, B, C, mcheck, &flg));
11128b400f6SJacob Faibussowitsch     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in C=A*B");
1129566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
113c20d7725SJed Brown 
114c20d7725SJed Brown     /* (1.2) Test user driver */
1159566063dSJacob Faibussowitsch     PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C));
116c20d7725SJed Brown 
117c20d7725SJed Brown     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
118c20d7725SJed Brown     alpha = 1.0;
119c20d7725SJed Brown     for (i = 0; i < 2; i++) {
120c20d7725SJed Brown       alpha -= 0.1;
1219566063dSJacob Faibussowitsch       PetscCall(MatScale(A, alpha));
1229566063dSJacob Faibussowitsch       PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C));
123c20d7725SJed Brown     }
1249566063dSJacob Faibussowitsch     PetscCall(MatMatMultEqual(A, B, C, mcheck, &flg));
12528b400f6SJacob Faibussowitsch     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult()");
1269566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&A));
1274417c5e8SHong Zhang 
1284417c5e8SHong Zhang     /* Test MatProductClear() */
1299566063dSJacob Faibussowitsch     PetscCall(MatProductClear(C));
1309566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
131544a5e07SHong Zhang 
132544a5e07SHong Zhang     /* Test MatMatMult() for dense and aij matrices */
1339566063dSJacob Faibussowitsch     PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &flg, MATSEQAIJ, MATMPIAIJ, ""));
13420b3374bSStefano Zampini     if (flg) {
1359566063dSJacob Faibussowitsch       PetscCall(MatConvert(A_save, MATDENSE, MAT_INITIAL_MATRIX, &A));
1369566063dSJacob Faibussowitsch       PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C));
1379566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C));
1389566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&A));
13920b3374bSStefano Zampini     }
140c20d7725SJed Brown   }
141c20d7725SJed Brown 
142c20d7725SJed Brown   /* Create P and R = P^T  */
143c20d7725SJed Brown   /* --------------------- */
1449566063dSJacob Faibussowitsch   PetscCall(MatGetSize(B, &PM, NULL));
14520b3374bSStefano Zampini   if (PN < 0) PN = PM / 2;
1469566063dSJacob Faibussowitsch   PetscCall(MatCreate(PETSC_COMM_WORLD, &P));
1479566063dSJacob Faibussowitsch   PetscCall(MatSetSizes(P, PETSC_DECIDE, PETSC_DECIDE, PM, PN));
1489566063dSJacob Faibussowitsch   PetscCall(MatSetType(P, MATAIJ));
1499566063dSJacob Faibussowitsch   PetscCall(MatSeqAIJSetPreallocation(P, nzp, NULL));
1509566063dSJacob Faibussowitsch   PetscCall(MatMPIAIJSetPreallocation(P, nzp, NULL, nzp, NULL));
1519566063dSJacob Faibussowitsch   PetscCall(MatGetOwnershipRange(P, &rstart, &rend));
1529371c9d4SSatish Balay   for (i = 0; i < nzp; i++) { PetscCall(PetscRandomGetValue(rdm, &a[i])); }
153c20d7725SJed Brown   for (i = rstart; i < rend; i++) {
154c20d7725SJed Brown     for (j = 0; j < nzp; j++) {
1559566063dSJacob Faibussowitsch       PetscCall(PetscRandomGetValue(rdm, &rval));
156c20d7725SJed Brown       idxn[j] = (PetscInt)(PetscRealPart(rval) * PN);
157c20d7725SJed Brown     }
1589566063dSJacob Faibussowitsch     PetscCall(MatSetValues(P, 1, &i, nzp, idxn, a, ADD_VALUES));
159c20d7725SJed Brown   }
1609566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
1619566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
162c20d7725SJed Brown 
1639566063dSJacob Faibussowitsch   PetscCall(MatTranspose(P, MAT_INITIAL_MATRIX, &R));
1649566063dSJacob Faibussowitsch   PetscCall(MatConvert(P, mattype, MAT_INPLACE_MATRIX, &P));
1659566063dSJacob Faibussowitsch   PetscCall(MatConvert(R, mattype, MAT_INPLACE_MATRIX, &R));
1669566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(P));
1679566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(R));
168c20d7725SJed Brown 
169c20d7725SJed Brown   /* 2) MatTransposeMatMult() */
170c20d7725SJed Brown   /* ------------------------ */
171c20d7725SJed Brown   if (Test_MatTrMat) {
172c20d7725SJed Brown     /* (2.1) Test developer driver C = P^T*B */
1739566063dSJacob Faibussowitsch     PetscCall(MatProductCreate(P, B, NULL, &C));
1749566063dSJacob Faibussowitsch     PetscCall(MatSetOptionsPrefix(C, "AtB_"));
1759566063dSJacob Faibussowitsch     PetscCall(MatProductSetType(C, MATPRODUCT_AtB));
1769566063dSJacob Faibussowitsch     PetscCall(MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT));
1779566063dSJacob Faibussowitsch     PetscCall(MatProductSetFill(C, PETSC_DEFAULT));
1789566063dSJacob Faibussowitsch     PetscCall(MatProductSetFromOptions(C));
1799566063dSJacob Faibussowitsch     PetscCall(MatHasOperation(C, MATOP_PRODUCTSYMBOLIC, &flg));
180263f2b91SStefano Zampini     if (flg) {                                                 /* run tests if supported */
1819566063dSJacob Faibussowitsch       PetscCall(MatProductSymbolic(C));                        /* equivalent to MatSetUp() */
1829566063dSJacob Faibussowitsch       PetscCall(MatSetOption(C, MAT_USE_INODES, PETSC_FALSE)); /* illustrate how to call MatSetOption() */
1839566063dSJacob Faibussowitsch       PetscCall(MatProductNumeric(C));
1849566063dSJacob Faibussowitsch       PetscCall(MatProductNumeric(C)); /* test reuse symbolic C */
18567b3012eSStefano Zampini 
1869566063dSJacob Faibussowitsch       PetscCall(MatTransposeMatMultEqual(P, B, C, mcheck, &flg));
18728b400f6SJacob Faibussowitsch       PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error: developer driver C = P^T*B");
1889566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C));
189c20d7725SJed Brown 
190c20d7725SJed Brown       /* (2.2) Test user driver C = P^T*B */
1919566063dSJacob Faibussowitsch       PetscCall(MatTransposeMatMult(P, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C));
1929566063dSJacob Faibussowitsch       PetscCall(MatTransposeMatMult(P, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C));
1939566063dSJacob Faibussowitsch       PetscCall(MatGetInfo(C, MAT_GLOBAL_SUM, &info));
1949566063dSJacob Faibussowitsch       PetscCall(MatProductClear(C));
195c20d7725SJed Brown 
196c20d7725SJed Brown       /* Compare P^T*B and R*B */
1979566063dSJacob Faibussowitsch       PetscCall(MatMatMult(R, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C1));
1989566063dSJacob Faibussowitsch       PetscCall(MatNormDifference(C, C1, &norm));
19908401ef6SPierre Jolivet       PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatTransposeMatMult(): %g", (double)norm);
2009566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C1));
201c20d7725SJed Brown 
202c20d7725SJed Brown       /* Test MatDuplicate() of C=P^T*B */
2039566063dSJacob Faibussowitsch       PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &C1));
2049566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C1));
205263f2b91SStefano Zampini     } else {
2069566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatTransposeMatMult not supported\n"));
207263f2b91SStefano Zampini     }
2089566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
209c20d7725SJed Brown   }
210c20d7725SJed Brown 
21167b3012eSStefano Zampini   /* 3) MatMatTransposeMult() */
212c20d7725SJed Brown   /* ------------------------ */
213c20d7725SJed Brown   if (Test_MatMatTr) {
214c20d7725SJed Brown     /* C = B*R^T */
2159566063dSJacob Faibussowitsch     PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATSEQAIJ, &seqaij));
21620b3374bSStefano Zampini     if (seqaij) {
2179566063dSJacob Faibussowitsch       PetscCall(MatMatTransposeMult(B, R, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C));
2189566063dSJacob Faibussowitsch       PetscCall(MatSetOptionsPrefix(C, "ABt_")); /* enable '-ABt_' for matrix C */
2199566063dSJacob Faibussowitsch       PetscCall(MatGetInfo(C, MAT_GLOBAL_SUM, &info));
220c20d7725SJed Brown 
221c20d7725SJed Brown       /* Test MAT_REUSE_MATRIX - reuse symbolic C */
2229566063dSJacob Faibussowitsch       PetscCall(MatMatTransposeMult(B, R, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C));
223c20d7725SJed Brown 
224c20d7725SJed Brown       /* Check */
2259566063dSJacob Faibussowitsch       PetscCall(MatMatMult(B, P, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C1));
2269566063dSJacob Faibussowitsch       PetscCall(MatNormDifference(C, C1, &norm));
22708401ef6SPierre Jolivet       PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatMatTransposeMult() %g", (double)norm);
2289566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C1));
2299566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C));
230c20d7725SJed Brown     }
231c20d7725SJed Brown   }
232c20d7725SJed Brown 
233c20d7725SJed Brown   /* 4) Test MatPtAP() */
234c20d7725SJed Brown   /*-------------------*/
235c20d7725SJed Brown   if (Test_MatPtAP) {
2369566063dSJacob Faibussowitsch     PetscCall(MatDuplicate(A_save, MAT_COPY_VALUES, &A));
237c20d7725SJed Brown 
238c20d7725SJed Brown     /* (4.1) Test developer API */
2399566063dSJacob Faibussowitsch     PetscCall(MatProductCreate(A, P, NULL, &C));
2409566063dSJacob Faibussowitsch     PetscCall(MatSetOptionsPrefix(C, "PtAP_"));
2419566063dSJacob Faibussowitsch     PetscCall(MatProductSetType(C, MATPRODUCT_PtAP));
2429566063dSJacob Faibussowitsch     PetscCall(MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT));
2439566063dSJacob Faibussowitsch     PetscCall(MatProductSetFill(C, PETSC_DEFAULT));
2449566063dSJacob Faibussowitsch     PetscCall(MatProductSetFromOptions(C));
2459566063dSJacob Faibussowitsch     PetscCall(MatProductSymbolic(C));
2469566063dSJacob Faibussowitsch     PetscCall(MatProductNumeric(C));
2479566063dSJacob Faibussowitsch     PetscCall(MatPtAPMultEqual(A, P, C, mcheck, &flg));
24828b400f6SJacob Faibussowitsch     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatProduct_PtAP");
2499566063dSJacob Faibussowitsch     PetscCall(MatProductNumeric(C)); /* reuse symbolic C */
250c20d7725SJed Brown 
2519566063dSJacob Faibussowitsch     PetscCall(MatPtAPMultEqual(A, P, C, mcheck, &flg));
25228b400f6SJacob Faibussowitsch     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatProduct_PtAP");
2539566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
254c20d7725SJed Brown 
255c20d7725SJed Brown     /* (4.2) Test user driver */
2569566063dSJacob Faibussowitsch     PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C));
257c20d7725SJed Brown 
258c20d7725SJed Brown     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
259c20d7725SJed Brown     alpha = 1.0;
260c20d7725SJed Brown     for (i = 0; i < 2; i++) {
261c20d7725SJed Brown       alpha -= 0.1;
2629566063dSJacob Faibussowitsch       PetscCall(MatScale(A, alpha));
2639566063dSJacob Faibussowitsch       PetscCall(MatPtAP(A, P, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C));
264c20d7725SJed Brown     }
2659566063dSJacob Faibussowitsch     PetscCall(MatPtAPMultEqual(A, P, C, mcheck, &flg));
26628b400f6SJacob Faibussowitsch     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatPtAP");
267c20d7725SJed Brown 
268c20d7725SJed Brown     /* 5) Test MatRARt() */
269c20d7725SJed Brown     /* ----------------- */
270c20d7725SJed Brown     if (Test_MatRARt) {
271c20d7725SJed Brown       Mat RARt;
27267b3012eSStefano Zampini 
27367b3012eSStefano Zampini       /* (5.1) Test developer driver RARt = R*A*Rt */
2749566063dSJacob Faibussowitsch       PetscCall(MatProductCreate(A, R, NULL, &RARt));
2759566063dSJacob Faibussowitsch       PetscCall(MatSetOptionsPrefix(RARt, "RARt_"));
2769566063dSJacob Faibussowitsch       PetscCall(MatProductSetType(RARt, MATPRODUCT_RARt));
2779566063dSJacob Faibussowitsch       PetscCall(MatProductSetAlgorithm(RARt, MATPRODUCTALGORITHMDEFAULT));
2789566063dSJacob Faibussowitsch       PetscCall(MatProductSetFill(RARt, PETSC_DEFAULT));
2799566063dSJacob Faibussowitsch       PetscCall(MatProductSetFromOptions(RARt));
2809566063dSJacob Faibussowitsch       PetscCall(MatHasOperation(RARt, MATOP_PRODUCTSYMBOLIC, &flg));
281263f2b91SStefano Zampini       if (flg) {
2829566063dSJacob Faibussowitsch         PetscCall(MatProductSymbolic(RARt));                        /* equivalent to MatSetUp() */
2839566063dSJacob Faibussowitsch         PetscCall(MatSetOption(RARt, MAT_USE_INODES, PETSC_FALSE)); /* illustrate how to call MatSetOption() */
2849566063dSJacob Faibussowitsch         PetscCall(MatProductNumeric(RARt));
2859566063dSJacob Faibussowitsch         PetscCall(MatProductNumeric(RARt)); /* test reuse symbolic RARt */
2869566063dSJacob Faibussowitsch         PetscCall(MatDestroy(&RARt));
28767b3012eSStefano Zampini 
28867b3012eSStefano Zampini         /* (2.2) Test user driver RARt = R*A*Rt */
2899566063dSJacob Faibussowitsch         PetscCall(MatRARt(A, R, MAT_INITIAL_MATRIX, 2.0, &RARt));
2909566063dSJacob Faibussowitsch         PetscCall(MatRARt(A, R, MAT_REUSE_MATRIX, 2.0, &RARt));
29167b3012eSStefano Zampini 
2929566063dSJacob Faibussowitsch         PetscCall(MatNormDifference(C, RARt, &norm));
29308401ef6SPierre Jolivet         PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "|PtAP - RARt| = %g", (double)norm);
294263f2b91SStefano Zampini       } else {
2959566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatRARt not supported\n"));
296263f2b91SStefano Zampini       }
2979566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&RARt));
298c20d7725SJed Brown     }
299c20d7725SJed Brown 
3009566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&A));
3019566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
302c20d7725SJed Brown   }
303c20d7725SJed Brown 
304c20d7725SJed Brown   /* Destroy objects */
3059566063dSJacob Faibussowitsch   PetscCall(PetscRandomDestroy(&rdm));
3069566063dSJacob Faibussowitsch   PetscCall(PetscFree(idxn));
307c20d7725SJed Brown 
3089566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&A_save));
3099566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&B));
3109566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&P));
3119566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&R));
312c20d7725SJed Brown 
3139566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
314b122ec5aSJacob Faibussowitsch   return 0;
315c20d7725SJed Brown }
316c20d7725SJed Brown 
317c20d7725SJed Brown /*TEST
318c20d7725SJed Brown    test:
319c20d7725SJed Brown      suffix: 1
320dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
321c20d7725SJed Brown      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium
322c20d7725SJed Brown      output_file: output/ex62_1.out
323c20d7725SJed Brown 
324c20d7725SJed Brown    test:
325c20d7725SJed Brown      suffix: 2_ab_scalable
326dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3273e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm scalable -matmatmult_via scalable -AtB_mat_product_algorithm outerproduct -mattransposematmult_via outerproduct
328c20d7725SJed Brown      output_file: output/ex62_1.out
329c20d7725SJed Brown 
330c20d7725SJed Brown    test:
331c20d7725SJed Brown      suffix: 3_ab_scalable_fast
332dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3333e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm scalable_fast -matmatmult_via scalable_fast -matmattransmult_via color
334c20d7725SJed Brown      output_file: output/ex62_1.out
335c20d7725SJed Brown 
336c20d7725SJed Brown    test:
337c20d7725SJed Brown      suffix: 4_ab_heap
338dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3393e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm heap -matmatmult_via heap -PtAP_mat_product_algorithm rap -matptap_via rap
340c20d7725SJed Brown      output_file: output/ex62_1.out
341c20d7725SJed Brown 
342c20d7725SJed Brown    test:
343c20d7725SJed Brown      suffix: 5_ab_btheap
344dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3453e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm btheap -matmatmult_via btheap -matrart_via r*art
346c20d7725SJed Brown      output_file: output/ex62_1.out
347c20d7725SJed Brown 
348c20d7725SJed Brown    test:
349c20d7725SJed Brown      suffix: 6_ab_llcondensed
350dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3513e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm llcondensed -matmatmult_via llcondensed -matrart_via coloring_rart
352c20d7725SJed Brown      output_file: output/ex62_1.out
353c20d7725SJed Brown 
354c20d7725SJed Brown    test:
355c20d7725SJed Brown      suffix: 7_ab_rowmerge
356dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3573e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm rowmerge -matmatmult_via rowmerge
358c20d7725SJed Brown      output_file: output/ex62_1.out
359c20d7725SJed Brown 
360c20d7725SJed Brown    test:
361c20d7725SJed Brown      suffix: 8_ab_hypre
362dfd57a17SPierre Jolivet      requires: hypre datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3633e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm hypre -matmatmult_via hypre -PtAP_mat_product_algorithm hypre -matptap_via hypre
364c20d7725SJed Brown      output_file: output/ex62_1.out
365c20d7725SJed Brown 
366c20d7725SJed Brown    test:
367263f2b91SStefano Zampini      suffix: hypre_medium
368263f2b91SStefano Zampini      nsize: {{1 3}}
369263f2b91SStefano Zampini      requires: hypre datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
370263f2b91SStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -A_mat_type hypre -B_mat_type hypre -test_rart 0
371263f2b91SStefano Zampini      output_file: output/ex62_hypre.out
372263f2b91SStefano Zampini 
373263f2b91SStefano Zampini    test:
374263f2b91SStefano Zampini      suffix: hypre_tiny
375263f2b91SStefano Zampini      nsize: {{1 3}}
376263f2b91SStefano Zampini      requires: hypre !complex double !defined(PETSC_USE_64BIT_INDICES)
377263f2b91SStefano Zampini      args: -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -A_mat_type hypre -B_mat_type hypre -test_rart 0
378263f2b91SStefano Zampini      output_file: output/ex62_hypre.out
379263f2b91SStefano Zampini 
380263f2b91SStefano Zampini    test:
381cec0a6c6SStefano Zampini      suffix: 9_mkl
38215df37b9SStefano Zampini      TODO: broken MatScale?
383b88df2e7SBarry Smith      requires: mkl_sparse datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
384cec0a6c6SStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -A_mat_type aijmkl -B_mat_type aijmkl
385cec0a6c6SStefano Zampini      output_file: output/ex62_1.out
386cec0a6c6SStefano Zampini 
387cec0a6c6SStefano Zampini    test:
388c20d7725SJed Brown      suffix: 10
389dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
390c20d7725SJed Brown      nsize: 3
391c20d7725SJed Brown      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium
392c20d7725SJed Brown      output_file: output/ex62_1.out
393c20d7725SJed Brown 
394c20d7725SJed Brown    test:
3953cea4f0aSStefano Zampini      suffix: 10_backend
396dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3973cea4f0aSStefano Zampini      nsize: 3
3983e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm backend -matmatmult_via backend -AtB_mat_product_algorithm backend -mattransposematmult_via backend -PtAP_mat_product_algorithm backend -matptap_via backend
3993cea4f0aSStefano Zampini      output_file: output/ex62_1.out
4003cea4f0aSStefano Zampini 
4013cea4f0aSStefano Zampini    test:
402c20d7725SJed Brown      suffix: 11_ab_scalable
403dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
404c20d7725SJed Brown      nsize: 3
4053e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm scalable -matmatmult_via scalable -AtB_mat_product_algorithm scalable -mattransposematmult_via scalable
406c20d7725SJed Brown      output_file: output/ex62_1.out
407c20d7725SJed Brown 
408c20d7725SJed Brown    test:
409c20d7725SJed Brown      suffix: 12_ab_seqmpi
410dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
411c20d7725SJed Brown      nsize: 3
4123e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm seqmpi -matmatmult_via seqmpi -AtB_mat_product_algorithm at*b -mattransposematmult_via at*b
413c20d7725SJed Brown      output_file: output/ex62_1.out
414c20d7725SJed Brown 
415c20d7725SJed Brown    test:
416c20d7725SJed Brown      suffix: 13_ab_hypre
417dfd57a17SPierre Jolivet      requires: hypre datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
418c20d7725SJed Brown      nsize: 3
4193e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm hypre -matmatmult_via hypre -PtAP_mat_product_algorithm hypre -matptap_via hypre
420c20d7725SJed Brown      output_file: output/ex62_1.out
421c20d7725SJed Brown 
42220b3374bSStefano Zampini    test:
42320b3374bSStefano Zampini      suffix: 14_seqaij
424dfd57a17SPierre Jolivet      requires: !complex double !defined(PETSC_USE_64BIT_INDICES)
42520b3374bSStefano Zampini      args: -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
42620b3374bSStefano Zampini      output_file: output/ex62_1.out
42720b3374bSStefano Zampini 
42820b3374bSStefano Zampini    test:
42920b3374bSStefano Zampini      suffix: 14_seqaijcusparse
430dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4311a2c6b5cSJunchao Zhang      args: -A_mat_type aijcusparse -B_mat_type aijcusparse -mat_form_explicit_transpose -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
43220b3374bSStefano Zampini      output_file: output/ex62_1.out
43320b3374bSStefano Zampini 
43420b3374bSStefano Zampini    test:
43565e4b4d4SStefano Zampini      suffix: 14_seqaijcusparse_cpu
436dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4373e662e0bSHong Zhang      args: -A_mat_type aijcusparse -B_mat_type aijcusparse -mat_form_explicit_transpose -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -AB_mat_product_algorithm_backend_cpu -matmatmult_backend_cpu -PtAP_mat_product_algorithm_backend_cpu -matptap_backend_cpu -RARt_mat_product_algorithm_backend_cpu -matrart_backend_cpu
43865e4b4d4SStefano Zampini      output_file: output/ex62_1.out
43965e4b4d4SStefano Zampini 
44065e4b4d4SStefano Zampini    test:
4415cb69cabSStefano Zampini      suffix: 14_mpiaijcusparse_seq
4425cb69cabSStefano Zampini      nsize: 1
443dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4441a2c6b5cSJunchao Zhang      args: -A_mat_type mpiaijcusparse -B_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
4455cb69cabSStefano Zampini      output_file: output/ex62_1.out
4465cb69cabSStefano Zampini 
4475cb69cabSStefano Zampini    test:
44865e4b4d4SStefano Zampini      suffix: 14_mpiaijcusparse_seq_cpu
44965e4b4d4SStefano Zampini      nsize: 1
450dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4514005f4d6SPierre Jolivet      args: -A_mat_type mpiaijcusparse -B_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -AB_mat_product_algorithm_backend_cpu -matmatmult_backend_cpu -PtAP_mat_product_algorithm_backend_cpu -matptap_backend_cpu -test_rart 0
45265e4b4d4SStefano Zampini      output_file: output/ex62_1.out
45365e4b4d4SStefano Zampini 
45465e4b4d4SStefano Zampini    test:
4555cb69cabSStefano Zampini      suffix: 14_mpiaijcusparse
4565cb69cabSStefano Zampini      nsize: 3
457dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4581a2c6b5cSJunchao Zhang      args: -A_mat_type mpiaijcusparse -B_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
4595cb69cabSStefano Zampini      output_file: output/ex62_1.out
4605cb69cabSStefano Zampini 
4615cb69cabSStefano Zampini    test:
46265e4b4d4SStefano Zampini      suffix: 14_mpiaijcusparse_cpu
46365e4b4d4SStefano Zampini      nsize: 3
464dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4654005f4d6SPierre Jolivet      args: -A_mat_type mpiaijcusparse -B_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -AB_mat_product_algorithm_backend_cpu -matmatmult_backend_cpu -PtAP_mat_product_algorithm_backend_cpu -matptap_backend_cpu -test_rart 0
46665e4b4d4SStefano Zampini      output_file: output/ex62_1.out
46765e4b4d4SStefano Zampini 
46865e4b4d4SStefano Zampini    test:
469076ba34aSJunchao Zhang      nsize: {{1 3}}
470076ba34aSJunchao Zhang      suffix: 14_aijkokkos
471*dcfd994dSJunchao Zhang      requires: kokkos_kernels !complex double !defined(PETSC_USE_64BIT_INDICES)
47220b3374bSStefano Zampini      args: -A_mat_type aijkokkos -B_mat_type aijkokkos -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
47320b3374bSStefano Zampini      output_file: output/ex62_1.out
47420b3374bSStefano Zampini 
4755cb69cabSStefano Zampini    # these tests use matrices with many zero rows
47620b3374bSStefano Zampini    test:
47720b3374bSStefano Zampini      suffix: 15_seqaijcusparse
478dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES) datafilespath
4791a2c6b5cSJunchao Zhang      args: -A_mat_type aijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
48020b3374bSStefano Zampini      output_file: output/ex62_1.out
4815cb69cabSStefano Zampini 
4825cb69cabSStefano Zampini    test:
4835cb69cabSStefano Zampini      suffix: 15_mpiaijcusparse_seq
4845cb69cabSStefano Zampini      nsize: 1
485dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES) datafilespath
4861a2c6b5cSJunchao Zhang      args: -A_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
4875cb69cabSStefano Zampini      output_file: output/ex62_1.out
4885cb69cabSStefano Zampini 
4895cb69cabSStefano Zampini    test:
4905cb69cabSStefano Zampini      nsize: 3
4915cb69cabSStefano Zampini      suffix: 15_mpiaijcusparse
492dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES) datafilespath
4931a2c6b5cSJunchao Zhang      args: -A_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
4945cb69cabSStefano Zampini      output_file: output/ex62_1.out
4955cb69cabSStefano Zampini 
496c20d7725SJed Brown TEST*/
497