xref: /petsc/src/mat/tests/ex62.c (revision 3ba1676111f5c958fe6c2729b46ca4d523958bb3)
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 */
16d71ae5a4SJacob Faibussowitsch PetscErrorCode MatNormDifference(Mat A, Mat B, PetscReal *norm)
17d71ae5a4SJacob Faibussowitsch {
18c20d7725SJed Brown   PetscFunctionBegin;
199566063dSJacob Faibussowitsch   PetscCall(MatAXPY(B, -1.0, A, DIFFERENT_NONZERO_PATTERN));
209566063dSJacob Faibussowitsch   PetscCall(MatNorm(B, NORM_FROBENIUS, norm));
21*3ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
22c20d7725SJed Brown }
23c20d7725SJed Brown 
24d71ae5a4SJacob Faibussowitsch int main(int argc, char **args)
25d71ae5a4SJacob Faibussowitsch {
26c20d7725SJed Brown   Mat          A, A_save, B, C, P, C1, R;
27c20d7725SJed Brown   PetscViewer  viewer;
28c20d7725SJed Brown   PetscMPIInt  size, rank;
2920b3374bSStefano Zampini   PetscInt     i, j, *idxn, PM, PN = PETSC_DECIDE, rstart, rend;
30c20d7725SJed Brown   PetscReal    norm;
31c20d7725SJed Brown   PetscRandom  rdm;
3220b3374bSStefano Zampini   char         file[2][PETSC_MAX_PATH_LEN] = {"", ""};
33c20d7725SJed Brown   PetscScalar *a, rval, alpha;
34c20d7725SJed Brown   PetscBool    Test_MatMatMult = PETSC_TRUE, Test_MatTrMat = PETSC_TRUE, Test_MatMatTr = PETSC_TRUE;
3520b3374bSStefano Zampini   PetscBool    Test_MatPtAP = PETSC_TRUE, Test_MatRARt = PETSC_TRUE, flg, seqaij, flgA, flgB;
36c20d7725SJed Brown   MatInfo      info;
3720b3374bSStefano Zampini   PetscInt     nzp = 5; /* num of nonzeros in each row of P */
38c20d7725SJed Brown   MatType      mattype;
3920b3374bSStefano Zampini   const char  *deft = MATAIJ;
4020b3374bSStefano Zampini   char         A_mattype[256], B_mattype[256];
4120b3374bSStefano Zampini   PetscInt     mcheck = 10;
42c20d7725SJed Brown 
43327415f7SBarry Smith   PetscFunctionBeginUser;
449566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
459566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
469566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
47c20d7725SJed Brown 
48c20d7725SJed Brown   /*  Load the matrices A_save and B */
49d0609cedSBarry Smith   PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "");
509566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-test_rart", "Test MatRARt", "", Test_MatRARt, &Test_MatRARt, NULL));
519566063dSJacob Faibussowitsch   PetscCall(PetscOptionsInt("-PN", "Number of columns of P", "", PN, &PN, NULL));
529566063dSJacob Faibussowitsch   PetscCall(PetscOptionsInt("-mcheck", "Number of matmult checks", "", mcheck, &mcheck, NULL));
539566063dSJacob Faibussowitsch   PetscCall(PetscOptionsString("-fA", "Path for matrix A", "", file[0], file[0], sizeof(file[0]), &flg));
5428b400f6SJacob Faibussowitsch   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Must indicate a file name for matrix A with the -fA option.");
559566063dSJacob Faibussowitsch   PetscCall(PetscOptionsString("-fB", "Path for matrix B", "", file[1], file[1], sizeof(file[1]), &flg));
569566063dSJacob Faibussowitsch   PetscCall(PetscOptionsFList("-A_mat_type", "Matrix type", "MatSetType", MatList, deft, A_mattype, 256, &flgA));
579566063dSJacob Faibussowitsch   PetscCall(PetscOptionsFList("-B_mat_type", "Matrix type", "MatSetType", MatList, deft, B_mattype, 256, &flgB));
58d0609cedSBarry Smith   PetscOptionsEnd();
59c20d7725SJed Brown 
609566063dSJacob Faibussowitsch   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[0], FILE_MODE_READ, &viewer));
619566063dSJacob Faibussowitsch   PetscCall(MatCreate(PETSC_COMM_WORLD, &A_save));
629566063dSJacob Faibussowitsch   PetscCall(MatLoad(A_save, viewer));
639566063dSJacob Faibussowitsch   PetscCall(PetscViewerDestroy(&viewer));
64c20d7725SJed Brown 
6520b3374bSStefano Zampini   if (flg) {
669566063dSJacob Faibussowitsch     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[1], FILE_MODE_READ, &viewer));
679566063dSJacob Faibussowitsch     PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
689566063dSJacob Faibussowitsch     PetscCall(MatLoad(B, viewer));
699566063dSJacob Faibussowitsch     PetscCall(PetscViewerDestroy(&viewer));
7020b3374bSStefano Zampini   } else {
719566063dSJacob Faibussowitsch     PetscCall(PetscObjectReference((PetscObject)A_save));
7220b3374bSStefano Zampini     B = A_save;
7320b3374bSStefano Zampini   }
7420b3374bSStefano Zampini 
7548a46eb9SPierre Jolivet   if (flgA) PetscCall(MatConvert(A_save, A_mattype, MAT_INPLACE_MATRIX, &A_save));
7648a46eb9SPierre Jolivet   if (flgB) PetscCall(MatConvert(B, B_mattype, MAT_INPLACE_MATRIX, &B));
779566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(A_save));
789566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(B));
79c20d7725SJed Brown 
809566063dSJacob Faibussowitsch   PetscCall(MatGetType(B, &mattype));
81c20d7725SJed Brown 
829566063dSJacob Faibussowitsch   PetscCall(PetscMalloc(nzp * (sizeof(PetscInt) + sizeof(PetscScalar)), &idxn));
83c20d7725SJed Brown   a = (PetscScalar *)(idxn + nzp);
84c20d7725SJed Brown 
859566063dSJacob Faibussowitsch   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
869566063dSJacob Faibussowitsch   PetscCall(PetscRandomSetFromOptions(rdm));
87c20d7725SJed Brown 
88c20d7725SJed Brown   /* 1) MatMatMult() */
89c20d7725SJed Brown   /* ----------------*/
90c20d7725SJed Brown   if (Test_MatMatMult) {
919566063dSJacob Faibussowitsch     PetscCall(MatDuplicate(A_save, MAT_COPY_VALUES, &A));
92c20d7725SJed Brown 
93c20d7725SJed Brown     /* (1.1) Test developer API */
949566063dSJacob Faibussowitsch     PetscCall(MatProductCreate(A, B, NULL, &C));
959566063dSJacob Faibussowitsch     PetscCall(MatSetOptionsPrefix(C, "AB_"));
969566063dSJacob Faibussowitsch     PetscCall(MatProductSetType(C, MATPRODUCT_AB));
979566063dSJacob Faibussowitsch     PetscCall(MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT));
989566063dSJacob Faibussowitsch     PetscCall(MatProductSetFill(C, PETSC_DEFAULT));
999566063dSJacob Faibussowitsch     PetscCall(MatProductSetFromOptions(C));
100910fa13fSStefano Zampini     /* we can inquire about MATOP_PRODUCTSYMBOLIC even if the destination matrix type has not been set yet */
1019566063dSJacob Faibussowitsch     PetscCall(MatHasOperation(C, MATOP_PRODUCTSYMBOLIC, &flg));
1029566063dSJacob Faibussowitsch     PetscCall(MatProductSymbolic(C));
1039566063dSJacob Faibussowitsch     PetscCall(MatProductNumeric(C));
1049566063dSJacob Faibussowitsch     PetscCall(MatMatMultEqual(A, B, C, mcheck, &flg));
10528b400f6SJacob Faibussowitsch     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in C=A*B");
106c20d7725SJed Brown 
107c20d7725SJed Brown     /* Test reuse symbolic C */
108c20d7725SJed Brown     alpha = 0.9;
1099566063dSJacob Faibussowitsch     PetscCall(MatScale(A, alpha));
1109566063dSJacob Faibussowitsch     PetscCall(MatProductNumeric(C));
111c20d7725SJed Brown 
1129566063dSJacob Faibussowitsch     PetscCall(MatMatMultEqual(A, B, C, mcheck, &flg));
11328b400f6SJacob Faibussowitsch     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in C=A*B");
1149566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
115c20d7725SJed Brown 
116c20d7725SJed Brown     /* (1.2) Test user driver */
1179566063dSJacob Faibussowitsch     PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C));
118c20d7725SJed Brown 
119c20d7725SJed Brown     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
120c20d7725SJed Brown     alpha = 1.0;
121c20d7725SJed Brown     for (i = 0; i < 2; i++) {
122c20d7725SJed Brown       alpha -= 0.1;
1239566063dSJacob Faibussowitsch       PetscCall(MatScale(A, alpha));
1249566063dSJacob Faibussowitsch       PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C));
125c20d7725SJed Brown     }
1269566063dSJacob Faibussowitsch     PetscCall(MatMatMultEqual(A, B, C, mcheck, &flg));
12728b400f6SJacob Faibussowitsch     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult()");
1289566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&A));
1294417c5e8SHong Zhang 
1304417c5e8SHong Zhang     /* Test MatProductClear() */
1319566063dSJacob Faibussowitsch     PetscCall(MatProductClear(C));
1329566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
133544a5e07SHong Zhang 
134544a5e07SHong Zhang     /* Test MatMatMult() for dense and aij matrices */
1359566063dSJacob Faibussowitsch     PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &flg, MATSEQAIJ, MATMPIAIJ, ""));
13620b3374bSStefano Zampini     if (flg) {
1379566063dSJacob Faibussowitsch       PetscCall(MatConvert(A_save, MATDENSE, MAT_INITIAL_MATRIX, &A));
1389566063dSJacob Faibussowitsch       PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C));
1399566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C));
1409566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&A));
14120b3374bSStefano Zampini     }
142c20d7725SJed Brown   }
143c20d7725SJed Brown 
144c20d7725SJed Brown   /* Create P and R = P^T  */
145c20d7725SJed Brown   /* --------------------- */
1469566063dSJacob Faibussowitsch   PetscCall(MatGetSize(B, &PM, NULL));
14720b3374bSStefano Zampini   if (PN < 0) PN = PM / 2;
1489566063dSJacob Faibussowitsch   PetscCall(MatCreate(PETSC_COMM_WORLD, &P));
1499566063dSJacob Faibussowitsch   PetscCall(MatSetSizes(P, PETSC_DECIDE, PETSC_DECIDE, PM, PN));
1509566063dSJacob Faibussowitsch   PetscCall(MatSetType(P, MATAIJ));
1519566063dSJacob Faibussowitsch   PetscCall(MatSeqAIJSetPreallocation(P, nzp, NULL));
1529566063dSJacob Faibussowitsch   PetscCall(MatMPIAIJSetPreallocation(P, nzp, NULL, nzp, NULL));
1539566063dSJacob Faibussowitsch   PetscCall(MatGetOwnershipRange(P, &rstart, &rend));
15448a46eb9SPierre Jolivet   for (i = 0; i < nzp; i++) PetscCall(PetscRandomGetValue(rdm, &a[i]));
155c20d7725SJed Brown   for (i = rstart; i < rend; i++) {
156c20d7725SJed Brown     for (j = 0; j < nzp; j++) {
1579566063dSJacob Faibussowitsch       PetscCall(PetscRandomGetValue(rdm, &rval));
158c20d7725SJed Brown       idxn[j] = (PetscInt)(PetscRealPart(rval) * PN);
159c20d7725SJed Brown     }
1609566063dSJacob Faibussowitsch     PetscCall(MatSetValues(P, 1, &i, nzp, idxn, a, ADD_VALUES));
161c20d7725SJed Brown   }
1629566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
1639566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
164c20d7725SJed Brown 
1659566063dSJacob Faibussowitsch   PetscCall(MatTranspose(P, MAT_INITIAL_MATRIX, &R));
1669566063dSJacob Faibussowitsch   PetscCall(MatConvert(P, mattype, MAT_INPLACE_MATRIX, &P));
1679566063dSJacob Faibussowitsch   PetscCall(MatConvert(R, mattype, MAT_INPLACE_MATRIX, &R));
1689566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(P));
1699566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(R));
170c20d7725SJed Brown 
171c20d7725SJed Brown   /* 2) MatTransposeMatMult() */
172c20d7725SJed Brown   /* ------------------------ */
173c20d7725SJed Brown   if (Test_MatTrMat) {
174c20d7725SJed Brown     /* (2.1) Test developer driver C = P^T*B */
1759566063dSJacob Faibussowitsch     PetscCall(MatProductCreate(P, B, NULL, &C));
1769566063dSJacob Faibussowitsch     PetscCall(MatSetOptionsPrefix(C, "AtB_"));
1779566063dSJacob Faibussowitsch     PetscCall(MatProductSetType(C, MATPRODUCT_AtB));
1789566063dSJacob Faibussowitsch     PetscCall(MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT));
1799566063dSJacob Faibussowitsch     PetscCall(MatProductSetFill(C, PETSC_DEFAULT));
1809566063dSJacob Faibussowitsch     PetscCall(MatProductSetFromOptions(C));
1819566063dSJacob Faibussowitsch     PetscCall(MatHasOperation(C, MATOP_PRODUCTSYMBOLIC, &flg));
182263f2b91SStefano Zampini     if (flg) {                                                 /* run tests if supported */
1839566063dSJacob Faibussowitsch       PetscCall(MatProductSymbolic(C));                        /* equivalent to MatSetUp() */
1849566063dSJacob Faibussowitsch       PetscCall(MatSetOption(C, MAT_USE_INODES, PETSC_FALSE)); /* illustrate how to call MatSetOption() */
1859566063dSJacob Faibussowitsch       PetscCall(MatProductNumeric(C));
1869566063dSJacob Faibussowitsch       PetscCall(MatProductNumeric(C)); /* test reuse symbolic C */
18767b3012eSStefano Zampini 
1889566063dSJacob Faibussowitsch       PetscCall(MatTransposeMatMultEqual(P, B, C, mcheck, &flg));
18928b400f6SJacob Faibussowitsch       PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error: developer driver C = P^T*B");
1909566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C));
191c20d7725SJed Brown 
192c20d7725SJed Brown       /* (2.2) Test user driver C = P^T*B */
1939566063dSJacob Faibussowitsch       PetscCall(MatTransposeMatMult(P, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C));
1949566063dSJacob Faibussowitsch       PetscCall(MatTransposeMatMult(P, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C));
1959566063dSJacob Faibussowitsch       PetscCall(MatGetInfo(C, MAT_GLOBAL_SUM, &info));
1969566063dSJacob Faibussowitsch       PetscCall(MatProductClear(C));
197c20d7725SJed Brown 
198c20d7725SJed Brown       /* Compare P^T*B and R*B */
1999566063dSJacob Faibussowitsch       PetscCall(MatMatMult(R, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C1));
2009566063dSJacob Faibussowitsch       PetscCall(MatNormDifference(C, C1, &norm));
20108401ef6SPierre Jolivet       PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatTransposeMatMult(): %g", (double)norm);
2029566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C1));
203c20d7725SJed Brown 
204c20d7725SJed Brown       /* Test MatDuplicate() of C=P^T*B */
2059566063dSJacob Faibussowitsch       PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &C1));
2069566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C1));
207263f2b91SStefano Zampini     } else {
2089566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatTransposeMatMult not supported\n"));
209263f2b91SStefano Zampini     }
2109566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
211c20d7725SJed Brown   }
212c20d7725SJed Brown 
21367b3012eSStefano Zampini   /* 3) MatMatTransposeMult() */
214c20d7725SJed Brown   /* ------------------------ */
215c20d7725SJed Brown   if (Test_MatMatTr) {
216c20d7725SJed Brown     /* C = B*R^T */
2179566063dSJacob Faibussowitsch     PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATSEQAIJ, &seqaij));
21820b3374bSStefano Zampini     if (seqaij) {
2199566063dSJacob Faibussowitsch       PetscCall(MatMatTransposeMult(B, R, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C));
2209566063dSJacob Faibussowitsch       PetscCall(MatSetOptionsPrefix(C, "ABt_")); /* enable '-ABt_' for matrix C */
2219566063dSJacob Faibussowitsch       PetscCall(MatGetInfo(C, MAT_GLOBAL_SUM, &info));
222c20d7725SJed Brown 
223c20d7725SJed Brown       /* Test MAT_REUSE_MATRIX - reuse symbolic C */
2249566063dSJacob Faibussowitsch       PetscCall(MatMatTransposeMult(B, R, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C));
225c20d7725SJed Brown 
226c20d7725SJed Brown       /* Check */
2279566063dSJacob Faibussowitsch       PetscCall(MatMatMult(B, P, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C1));
2289566063dSJacob Faibussowitsch       PetscCall(MatNormDifference(C, C1, &norm));
22908401ef6SPierre Jolivet       PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatMatTransposeMult() %g", (double)norm);
2309566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C1));
2319566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C));
232c20d7725SJed Brown     }
233c20d7725SJed Brown   }
234c20d7725SJed Brown 
235c20d7725SJed Brown   /* 4) Test MatPtAP() */
236c20d7725SJed Brown   /*-------------------*/
237c20d7725SJed Brown   if (Test_MatPtAP) {
2389566063dSJacob Faibussowitsch     PetscCall(MatDuplicate(A_save, MAT_COPY_VALUES, &A));
239c20d7725SJed Brown 
240c20d7725SJed Brown     /* (4.1) Test developer API */
2419566063dSJacob Faibussowitsch     PetscCall(MatProductCreate(A, P, NULL, &C));
2429566063dSJacob Faibussowitsch     PetscCall(MatSetOptionsPrefix(C, "PtAP_"));
2439566063dSJacob Faibussowitsch     PetscCall(MatProductSetType(C, MATPRODUCT_PtAP));
2449566063dSJacob Faibussowitsch     PetscCall(MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT));
2459566063dSJacob Faibussowitsch     PetscCall(MatProductSetFill(C, PETSC_DEFAULT));
2469566063dSJacob Faibussowitsch     PetscCall(MatProductSetFromOptions(C));
2479566063dSJacob Faibussowitsch     PetscCall(MatProductSymbolic(C));
2489566063dSJacob Faibussowitsch     PetscCall(MatProductNumeric(C));
2499566063dSJacob Faibussowitsch     PetscCall(MatPtAPMultEqual(A, P, C, mcheck, &flg));
25028b400f6SJacob Faibussowitsch     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatProduct_PtAP");
2519566063dSJacob Faibussowitsch     PetscCall(MatProductNumeric(C)); /* reuse symbolic C */
252c20d7725SJed Brown 
2539566063dSJacob Faibussowitsch     PetscCall(MatPtAPMultEqual(A, P, C, mcheck, &flg));
25428b400f6SJacob Faibussowitsch     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatProduct_PtAP");
2559566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
256c20d7725SJed Brown 
257c20d7725SJed Brown     /* (4.2) Test user driver */
2589566063dSJacob Faibussowitsch     PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C));
259c20d7725SJed Brown 
260c20d7725SJed Brown     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
261c20d7725SJed Brown     alpha = 1.0;
262c20d7725SJed Brown     for (i = 0; i < 2; i++) {
263c20d7725SJed Brown       alpha -= 0.1;
2649566063dSJacob Faibussowitsch       PetscCall(MatScale(A, alpha));
2659566063dSJacob Faibussowitsch       PetscCall(MatPtAP(A, P, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C));
266c20d7725SJed Brown     }
2679566063dSJacob Faibussowitsch     PetscCall(MatPtAPMultEqual(A, P, C, mcheck, &flg));
26828b400f6SJacob Faibussowitsch     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatPtAP");
269c20d7725SJed Brown 
270c20d7725SJed Brown     /* 5) Test MatRARt() */
271c20d7725SJed Brown     /* ----------------- */
272c20d7725SJed Brown     if (Test_MatRARt) {
273c20d7725SJed Brown       Mat RARt;
27467b3012eSStefano Zampini 
27567b3012eSStefano Zampini       /* (5.1) Test developer driver RARt = R*A*Rt */
2769566063dSJacob Faibussowitsch       PetscCall(MatProductCreate(A, R, NULL, &RARt));
2779566063dSJacob Faibussowitsch       PetscCall(MatSetOptionsPrefix(RARt, "RARt_"));
2789566063dSJacob Faibussowitsch       PetscCall(MatProductSetType(RARt, MATPRODUCT_RARt));
2799566063dSJacob Faibussowitsch       PetscCall(MatProductSetAlgorithm(RARt, MATPRODUCTALGORITHMDEFAULT));
2809566063dSJacob Faibussowitsch       PetscCall(MatProductSetFill(RARt, PETSC_DEFAULT));
2819566063dSJacob Faibussowitsch       PetscCall(MatProductSetFromOptions(RARt));
2829566063dSJacob Faibussowitsch       PetscCall(MatHasOperation(RARt, MATOP_PRODUCTSYMBOLIC, &flg));
283263f2b91SStefano Zampini       if (flg) {
2849566063dSJacob Faibussowitsch         PetscCall(MatProductSymbolic(RARt));                        /* equivalent to MatSetUp() */
2859566063dSJacob Faibussowitsch         PetscCall(MatSetOption(RARt, MAT_USE_INODES, PETSC_FALSE)); /* illustrate how to call MatSetOption() */
2869566063dSJacob Faibussowitsch         PetscCall(MatProductNumeric(RARt));
2879566063dSJacob Faibussowitsch         PetscCall(MatProductNumeric(RARt)); /* test reuse symbolic RARt */
2889566063dSJacob Faibussowitsch         PetscCall(MatDestroy(&RARt));
28967b3012eSStefano Zampini 
29067b3012eSStefano Zampini         /* (2.2) Test user driver RARt = R*A*Rt */
2919566063dSJacob Faibussowitsch         PetscCall(MatRARt(A, R, MAT_INITIAL_MATRIX, 2.0, &RARt));
2929566063dSJacob Faibussowitsch         PetscCall(MatRARt(A, R, MAT_REUSE_MATRIX, 2.0, &RARt));
29367b3012eSStefano Zampini 
2949566063dSJacob Faibussowitsch         PetscCall(MatNormDifference(C, RARt, &norm));
29508401ef6SPierre Jolivet         PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "|PtAP - RARt| = %g", (double)norm);
296263f2b91SStefano Zampini       } else {
2979566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatRARt not supported\n"));
298263f2b91SStefano Zampini       }
2999566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&RARt));
300c20d7725SJed Brown     }
301c20d7725SJed Brown 
3029566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&A));
3039566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
304c20d7725SJed Brown   }
305c20d7725SJed Brown 
306c20d7725SJed Brown   /* Destroy objects */
3079566063dSJacob Faibussowitsch   PetscCall(PetscRandomDestroy(&rdm));
3089566063dSJacob Faibussowitsch   PetscCall(PetscFree(idxn));
309c20d7725SJed Brown 
3109566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&A_save));
3119566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&B));
3129566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&P));
3139566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&R));
314c20d7725SJed Brown 
3159566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
316b122ec5aSJacob Faibussowitsch   return 0;
317c20d7725SJed Brown }
318c20d7725SJed Brown 
319c20d7725SJed Brown /*TEST
320c20d7725SJed Brown    test:
321c20d7725SJed Brown      suffix: 1
322dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
323c20d7725SJed Brown      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium
324c20d7725SJed Brown      output_file: output/ex62_1.out
325c20d7725SJed Brown 
326c20d7725SJed Brown    test:
327c20d7725SJed Brown      suffix: 2_ab_scalable
328dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3293e662e0bSHong 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
330c20d7725SJed Brown      output_file: output/ex62_1.out
331c20d7725SJed Brown 
332c20d7725SJed Brown    test:
333c20d7725SJed Brown      suffix: 3_ab_scalable_fast
334dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3353e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm scalable_fast -matmatmult_via scalable_fast -matmattransmult_via color
336c20d7725SJed Brown      output_file: output/ex62_1.out
337c20d7725SJed Brown 
338c20d7725SJed Brown    test:
339c20d7725SJed Brown      suffix: 4_ab_heap
340dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3413e662e0bSHong 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
342c20d7725SJed Brown      output_file: output/ex62_1.out
343c20d7725SJed Brown 
344c20d7725SJed Brown    test:
345c20d7725SJed Brown      suffix: 5_ab_btheap
346dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3473e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm btheap -matmatmult_via btheap -matrart_via r*art
348c20d7725SJed Brown      output_file: output/ex62_1.out
349c20d7725SJed Brown 
350c20d7725SJed Brown    test:
351c20d7725SJed Brown      suffix: 6_ab_llcondensed
352dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3533e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm llcondensed -matmatmult_via llcondensed -matrart_via coloring_rart
354c20d7725SJed Brown      output_file: output/ex62_1.out
355c20d7725SJed Brown 
356c20d7725SJed Brown    test:
357c20d7725SJed Brown      suffix: 7_ab_rowmerge
358dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3593e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm rowmerge -matmatmult_via rowmerge
360c20d7725SJed Brown      output_file: output/ex62_1.out
361c20d7725SJed Brown 
362c20d7725SJed Brown    test:
363c20d7725SJed Brown      suffix: 8_ab_hypre
364dfd57a17SPierre Jolivet      requires: hypre datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3653e662e0bSHong 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
366c20d7725SJed Brown      output_file: output/ex62_1.out
367c20d7725SJed Brown 
368c20d7725SJed Brown    test:
369263f2b91SStefano Zampini      suffix: hypre_medium
370263f2b91SStefano Zampini      nsize: {{1 3}}
371263f2b91SStefano Zampini      requires: hypre datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
372263f2b91SStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -A_mat_type hypre -B_mat_type hypre -test_rart 0
373263f2b91SStefano Zampini      output_file: output/ex62_hypre.out
374263f2b91SStefano Zampini 
375263f2b91SStefano Zampini    test:
376263f2b91SStefano Zampini      suffix: hypre_tiny
377263f2b91SStefano Zampini      nsize: {{1 3}}
378263f2b91SStefano Zampini      requires: hypre !complex double !defined(PETSC_USE_64BIT_INDICES)
379263f2b91SStefano 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
380263f2b91SStefano Zampini      output_file: output/ex62_hypre.out
381263f2b91SStefano Zampini 
382263f2b91SStefano Zampini    test:
383cec0a6c6SStefano Zampini      suffix: 9_mkl
38415df37b9SStefano Zampini      TODO: broken MatScale?
385b88df2e7SBarry Smith      requires: mkl_sparse datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
386cec0a6c6SStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -A_mat_type aijmkl -B_mat_type aijmkl
387cec0a6c6SStefano Zampini      output_file: output/ex62_1.out
388cec0a6c6SStefano Zampini 
389cec0a6c6SStefano Zampini    test:
390c20d7725SJed Brown      suffix: 10
391dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
392c20d7725SJed Brown      nsize: 3
393c20d7725SJed Brown      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium
394c20d7725SJed Brown      output_file: output/ex62_1.out
395c20d7725SJed Brown 
396c20d7725SJed Brown    test:
3973cea4f0aSStefano Zampini      suffix: 10_backend
398dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3993cea4f0aSStefano Zampini      nsize: 3
4003e662e0bSHong 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
4013cea4f0aSStefano Zampini      output_file: output/ex62_1.out
4023cea4f0aSStefano Zampini 
4033cea4f0aSStefano Zampini    test:
404c20d7725SJed Brown      suffix: 11_ab_scalable
405dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
406c20d7725SJed Brown      nsize: 3
4073e662e0bSHong 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
408c20d7725SJed Brown      output_file: output/ex62_1.out
409c20d7725SJed Brown 
410c20d7725SJed Brown    test:
411c20d7725SJed Brown      suffix: 12_ab_seqmpi
412dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
413c20d7725SJed Brown      nsize: 3
4143e662e0bSHong 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
415c20d7725SJed Brown      output_file: output/ex62_1.out
416c20d7725SJed Brown 
417c20d7725SJed Brown    test:
418c20d7725SJed Brown      suffix: 13_ab_hypre
419dfd57a17SPierre Jolivet      requires: hypre datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
420c20d7725SJed Brown      nsize: 3
4213e662e0bSHong 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
422c20d7725SJed Brown      output_file: output/ex62_1.out
423c20d7725SJed Brown 
42420b3374bSStefano Zampini    test:
42520b3374bSStefano Zampini      suffix: 14_seqaij
426dfd57a17SPierre Jolivet      requires: !complex double !defined(PETSC_USE_64BIT_INDICES)
42720b3374bSStefano Zampini      args: -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
42820b3374bSStefano Zampini      output_file: output/ex62_1.out
42920b3374bSStefano Zampini 
43020b3374bSStefano Zampini    test:
43120b3374bSStefano Zampini      suffix: 14_seqaijcusparse
432dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4331a2c6b5cSJunchao 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
43420b3374bSStefano Zampini      output_file: output/ex62_1.out
43520b3374bSStefano Zampini 
43620b3374bSStefano Zampini    test:
43765e4b4d4SStefano Zampini      suffix: 14_seqaijcusparse_cpu
438dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4393e662e0bSHong 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
44065e4b4d4SStefano Zampini      output_file: output/ex62_1.out
44165e4b4d4SStefano Zampini 
44265e4b4d4SStefano Zampini    test:
4435cb69cabSStefano Zampini      suffix: 14_mpiaijcusparse_seq
4445cb69cabSStefano Zampini      nsize: 1
445dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4461a2c6b5cSJunchao 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
4475cb69cabSStefano Zampini      output_file: output/ex62_1.out
4485cb69cabSStefano Zampini 
4495cb69cabSStefano Zampini    test:
45065e4b4d4SStefano Zampini      suffix: 14_mpiaijcusparse_seq_cpu
45165e4b4d4SStefano Zampini      nsize: 1
452dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4534005f4d6SPierre 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
45465e4b4d4SStefano Zampini      output_file: output/ex62_1.out
45565e4b4d4SStefano Zampini 
45665e4b4d4SStefano Zampini    test:
4575cb69cabSStefano Zampini      suffix: 14_mpiaijcusparse
4585cb69cabSStefano Zampini      nsize: 3
459dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4601a2c6b5cSJunchao 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
4615cb69cabSStefano Zampini      output_file: output/ex62_1.out
4625cb69cabSStefano Zampini 
4635cb69cabSStefano Zampini    test:
46465e4b4d4SStefano Zampini      suffix: 14_mpiaijcusparse_cpu
46565e4b4d4SStefano Zampini      nsize: 3
466dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4674005f4d6SPierre 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
46865e4b4d4SStefano Zampini      output_file: output/ex62_1.out
46965e4b4d4SStefano Zampini 
47065e4b4d4SStefano Zampini    test:
471076ba34aSJunchao Zhang      nsize: {{1 3}}
472076ba34aSJunchao Zhang      suffix: 14_aijkokkos
473dcfd994dSJunchao Zhang      requires: kokkos_kernels !complex double !defined(PETSC_USE_64BIT_INDICES)
47420b3374bSStefano 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
47520b3374bSStefano Zampini      output_file: output/ex62_1.out
47620b3374bSStefano Zampini 
4775cb69cabSStefano Zampini    # these tests use matrices with many zero rows
47820b3374bSStefano Zampini    test:
47920b3374bSStefano Zampini      suffix: 15_seqaijcusparse
480dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES) datafilespath
4811a2c6b5cSJunchao Zhang      args: -A_mat_type aijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
48220b3374bSStefano Zampini      output_file: output/ex62_1.out
4835cb69cabSStefano Zampini 
4845cb69cabSStefano Zampini    test:
4855cb69cabSStefano Zampini      suffix: 15_mpiaijcusparse_seq
4865cb69cabSStefano Zampini      nsize: 1
487dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES) datafilespath
4881a2c6b5cSJunchao Zhang      args: -A_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
4895cb69cabSStefano Zampini      output_file: output/ex62_1.out
4905cb69cabSStefano Zampini 
4915cb69cabSStefano Zampini    test:
4925cb69cabSStefano Zampini      nsize: 3
4935cb69cabSStefano Zampini      suffix: 15_mpiaijcusparse
494dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES) datafilespath
4951a2c6b5cSJunchao Zhang      args: -A_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
4965cb69cabSStefano Zampini      output_file: output/ex62_1.out
4975cb69cabSStefano Zampini 
498c20d7725SJed Brown TEST*/
499