xref: /petsc/src/mat/tests/ex62.c (revision 08401ef684002a709c6d3db98a0c9f54a8bcf1ec)
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 */
16c20d7725SJed Brown PetscErrorCode MatNormDifference(Mat A,Mat B,PetscReal *norm)
17c20d7725SJed Brown {
18c20d7725SJed Brown   PetscFunctionBegin;
199566063dSJacob Faibussowitsch   PetscCall(MatAXPY(B,-1.0,A,DIFFERENT_NONZERO_PATTERN));
209566063dSJacob Faibussowitsch   PetscCall(MatNorm(B,NORM_FROBENIUS,norm));
21c20d7725SJed Brown   PetscFunctionReturn(0);
22c20d7725SJed Brown }
23c20d7725SJed Brown 
24c20d7725SJed Brown int main(int argc,char **args)
25c20d7725SJed Brown {
26c20d7725SJed Brown   Mat            A,A_save,B,C,P,C1,R;
27c20d7725SJed Brown   PetscViewer    viewer;
28c20d7725SJed Brown   PetscErrorCode ierr;
29c20d7725SJed Brown   PetscMPIInt    size,rank;
3020b3374bSStefano Zampini   PetscInt       i,j,*idxn,PM,PN = PETSC_DECIDE,rstart,rend;
31c20d7725SJed Brown   PetscReal      norm;
32c20d7725SJed Brown   PetscRandom    rdm;
3320b3374bSStefano Zampini   char           file[2][PETSC_MAX_PATH_LEN] = { "", ""};
34c20d7725SJed Brown   PetscScalar    *a,rval,alpha;
35c20d7725SJed Brown   PetscBool      Test_MatMatMult=PETSC_TRUE,Test_MatTrMat=PETSC_TRUE,Test_MatMatTr=PETSC_TRUE;
3620b3374bSStefano Zampini   PetscBool      Test_MatPtAP=PETSC_TRUE,Test_MatRARt=PETSC_TRUE,flg,seqaij,flgA,flgB;
37c20d7725SJed Brown   MatInfo        info;
3820b3374bSStefano Zampini   PetscInt       nzp  = 5; /* num of nonzeros in each row of P */
39c20d7725SJed Brown   MatType        mattype;
4020b3374bSStefano Zampini   const char     *deft = MATAIJ;
4120b3374bSStefano Zampini   char           A_mattype[256], B_mattype[256];
4220b3374bSStefano Zampini   PetscInt       mcheck = 10;
43c20d7725SJed Brown 
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 */
499566063dSJacob Faibussowitsch   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","","");PetscCall(ierr);
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));
589566063dSJacob Faibussowitsch   ierr = PetscOptionsEnd();PetscCall(ierr);
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 
7520b3374bSStefano Zampini   if (flgA) {
769566063dSJacob Faibussowitsch     PetscCall(MatConvert(A_save,A_mattype,MAT_INPLACE_MATRIX,&A_save));
7720b3374bSStefano Zampini   }
7820b3374bSStefano Zampini   if (flgB) {
799566063dSJacob Faibussowitsch     PetscCall(MatConvert(B,B_mattype,MAT_INPLACE_MATRIX,&B));
8020b3374bSStefano Zampini   }
819566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(A_save));
829566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(B));
83c20d7725SJed Brown 
849566063dSJacob Faibussowitsch   PetscCall(MatGetType(B,&mattype));
85c20d7725SJed Brown 
869566063dSJacob Faibussowitsch   PetscCall(PetscMalloc(nzp*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn));
87c20d7725SJed Brown   a    = (PetscScalar*)(idxn + nzp);
88c20d7725SJed Brown 
899566063dSJacob Faibussowitsch   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rdm));
909566063dSJacob Faibussowitsch   PetscCall(PetscRandomSetFromOptions(rdm));
91c20d7725SJed Brown 
92c20d7725SJed Brown   /* 1) MatMatMult() */
93c20d7725SJed Brown   /* ----------------*/
94c20d7725SJed Brown   if (Test_MatMatMult) {
959566063dSJacob Faibussowitsch     PetscCall(MatDuplicate(A_save,MAT_COPY_VALUES,&A));
96c20d7725SJed Brown 
97c20d7725SJed Brown     /* (1.1) Test developer API */
989566063dSJacob Faibussowitsch     PetscCall(MatProductCreate(A,B,NULL,&C));
999566063dSJacob Faibussowitsch     PetscCall(MatSetOptionsPrefix(C,"AB_"));
1009566063dSJacob Faibussowitsch     PetscCall(MatProductSetType(C,MATPRODUCT_AB));
1019566063dSJacob Faibussowitsch     PetscCall(MatProductSetAlgorithm(C,MATPRODUCTALGORITHMDEFAULT));
1029566063dSJacob Faibussowitsch     PetscCall(MatProductSetFill(C,PETSC_DEFAULT));
1039566063dSJacob Faibussowitsch     PetscCall(MatProductSetFromOptions(C));
104910fa13fSStefano Zampini     /* we can inquire about MATOP_PRODUCTSYMBOLIC even if the destination matrix type has not been set yet */
1059566063dSJacob Faibussowitsch     PetscCall(MatHasOperation(C,MATOP_PRODUCTSYMBOLIC,&flg));
1069566063dSJacob Faibussowitsch     PetscCall(MatProductSymbolic(C));
1079566063dSJacob Faibussowitsch     PetscCall(MatProductNumeric(C));
1089566063dSJacob Faibussowitsch     PetscCall(MatMatMultEqual(A,B,C,mcheck,&flg));
10928b400f6SJacob Faibussowitsch     PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in C=A*B");
110c20d7725SJed Brown 
111c20d7725SJed Brown     /* Test reuse symbolic C */
112c20d7725SJed Brown     alpha = 0.9;
1139566063dSJacob Faibussowitsch     PetscCall(MatScale(A,alpha));
1149566063dSJacob Faibussowitsch     PetscCall(MatProductNumeric(C));
115c20d7725SJed Brown 
1169566063dSJacob Faibussowitsch     PetscCall(MatMatMultEqual(A,B,C,mcheck,&flg));
11728b400f6SJacob Faibussowitsch     PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in C=A*B");
1189566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
119c20d7725SJed Brown 
120c20d7725SJed Brown     /* (1.2) Test user driver */
1219566063dSJacob Faibussowitsch     PetscCall(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C));
122c20d7725SJed Brown 
123c20d7725SJed Brown     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
124c20d7725SJed Brown     alpha = 1.0;
125c20d7725SJed Brown     for (i=0; i<2; i++) {
126c20d7725SJed Brown       alpha -= 0.1;
1279566063dSJacob Faibussowitsch       PetscCall(MatScale(A,alpha));
1289566063dSJacob Faibussowitsch       PetscCall(MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C));
129c20d7725SJed Brown     }
1309566063dSJacob Faibussowitsch     PetscCall(MatMatMultEqual(A,B,C,mcheck,&flg));
13128b400f6SJacob Faibussowitsch     PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult()");
1329566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&A));
1334417c5e8SHong Zhang 
1344417c5e8SHong Zhang     /* Test MatProductClear() */
1359566063dSJacob Faibussowitsch     PetscCall(MatProductClear(C));
1369566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
137544a5e07SHong Zhang 
138544a5e07SHong Zhang     /* Test MatMatMult() for dense and aij matrices */
1399566063dSJacob Faibussowitsch     PetscCall(PetscObjectTypeCompareAny((PetscObject)A,&flg,MATSEQAIJ,MATMPIAIJ,""));
14020b3374bSStefano Zampini     if (flg) {
1419566063dSJacob Faibussowitsch       PetscCall(MatConvert(A_save,MATDENSE,MAT_INITIAL_MATRIX,&A));
1429566063dSJacob Faibussowitsch       PetscCall(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C));
1439566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C));
1449566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&A));
14520b3374bSStefano Zampini     }
146c20d7725SJed Brown   }
147c20d7725SJed Brown 
148c20d7725SJed Brown   /* Create P and R = P^T  */
149c20d7725SJed Brown   /* --------------------- */
1509566063dSJacob Faibussowitsch   PetscCall(MatGetSize(B,&PM,NULL));
15120b3374bSStefano Zampini   if (PN < 0) PN = PM/2;
1529566063dSJacob Faibussowitsch   PetscCall(MatCreate(PETSC_COMM_WORLD,&P));
1539566063dSJacob Faibussowitsch   PetscCall(MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,PM,PN));
1549566063dSJacob Faibussowitsch   PetscCall(MatSetType(P,MATAIJ));
1559566063dSJacob Faibussowitsch   PetscCall(MatSeqAIJSetPreallocation(P,nzp,NULL));
1569566063dSJacob Faibussowitsch   PetscCall(MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL));
1579566063dSJacob Faibussowitsch   PetscCall(MatGetOwnershipRange(P,&rstart,&rend));
158c20d7725SJed Brown   for (i=0; i<nzp; i++) {
1599566063dSJacob Faibussowitsch     PetscCall(PetscRandomGetValue(rdm,&a[i]));
160c20d7725SJed Brown   }
161c20d7725SJed Brown   for (i=rstart; i<rend; i++) {
162c20d7725SJed Brown     for (j=0; j<nzp; j++) {
1639566063dSJacob Faibussowitsch       PetscCall(PetscRandomGetValue(rdm,&rval));
164c20d7725SJed Brown       idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
165c20d7725SJed Brown     }
1669566063dSJacob Faibussowitsch     PetscCall(MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES));
167c20d7725SJed Brown   }
1689566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY));
1699566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY));
170c20d7725SJed Brown 
1719566063dSJacob Faibussowitsch   PetscCall(MatTranspose(P,MAT_INITIAL_MATRIX,&R));
1729566063dSJacob Faibussowitsch   PetscCall(MatConvert(P,mattype,MAT_INPLACE_MATRIX,&P));
1739566063dSJacob Faibussowitsch   PetscCall(MatConvert(R,mattype,MAT_INPLACE_MATRIX,&R));
1749566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(P));
1759566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(R));
176c20d7725SJed Brown 
177c20d7725SJed Brown   /* 2) MatTransposeMatMult() */
178c20d7725SJed Brown   /* ------------------------ */
179c20d7725SJed Brown   if (Test_MatTrMat) {
180c20d7725SJed Brown     /* (2.1) Test developer driver C = P^T*B */
1819566063dSJacob Faibussowitsch     PetscCall(MatProductCreate(P,B,NULL,&C));
1829566063dSJacob Faibussowitsch     PetscCall(MatSetOptionsPrefix(C,"AtB_"));
1839566063dSJacob Faibussowitsch     PetscCall(MatProductSetType(C,MATPRODUCT_AtB));
1849566063dSJacob Faibussowitsch     PetscCall(MatProductSetAlgorithm(C,MATPRODUCTALGORITHMDEFAULT));
1859566063dSJacob Faibussowitsch     PetscCall(MatProductSetFill(C,PETSC_DEFAULT));
1869566063dSJacob Faibussowitsch     PetscCall(MatProductSetFromOptions(C));
1879566063dSJacob Faibussowitsch     PetscCall(MatHasOperation(C,MATOP_PRODUCTSYMBOLIC,&flg));
188263f2b91SStefano Zampini     if (flg) { /* run tests if supported */
1899566063dSJacob Faibussowitsch       PetscCall(MatProductSymbolic(C)); /* equivalent to MatSetUp() */
1909566063dSJacob Faibussowitsch       PetscCall(MatSetOption(C,MAT_USE_INODES,PETSC_FALSE)); /* illustrate how to call MatSetOption() */
1919566063dSJacob Faibussowitsch       PetscCall(MatProductNumeric(C));
1929566063dSJacob Faibussowitsch       PetscCall(MatProductNumeric(C)); /* test reuse symbolic C */
19367b3012eSStefano Zampini 
1949566063dSJacob Faibussowitsch       PetscCall(MatTransposeMatMultEqual(P,B,C,mcheck,&flg));
19528b400f6SJacob Faibussowitsch       PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error: developer driver C = P^T*B");
1969566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C));
197c20d7725SJed Brown 
198c20d7725SJed Brown       /* (2.2) Test user driver C = P^T*B */
1999566063dSJacob Faibussowitsch       PetscCall(MatTransposeMatMult(P,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C));
2009566063dSJacob Faibussowitsch       PetscCall(MatTransposeMatMult(P,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C));
2019566063dSJacob Faibussowitsch       PetscCall(MatGetInfo(C,MAT_GLOBAL_SUM,&info));
2029566063dSJacob Faibussowitsch       PetscCall(MatProductClear(C));
203c20d7725SJed Brown 
204c20d7725SJed Brown       /* Compare P^T*B and R*B */
2059566063dSJacob Faibussowitsch       PetscCall(MatMatMult(R,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C1));
2069566063dSJacob Faibussowitsch       PetscCall(MatNormDifference(C,C1,&norm));
207*08401ef6SPierre Jolivet       PetscCheck(norm <= PETSC_SMALL,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatTransposeMatMult(): %g",(double)norm);
2089566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C1));
209c20d7725SJed Brown 
210c20d7725SJed Brown       /* Test MatDuplicate() of C=P^T*B */
2119566063dSJacob Faibussowitsch       PetscCall(MatDuplicate(C,MAT_COPY_VALUES,&C1));
2129566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C1));
213263f2b91SStefano Zampini     } else {
2149566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatTransposeMatMult not supported\n"));
215263f2b91SStefano Zampini     }
2169566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
217c20d7725SJed Brown   }
218c20d7725SJed Brown 
21967b3012eSStefano Zampini   /* 3) MatMatTransposeMult() */
220c20d7725SJed Brown   /* ------------------------ */
221c20d7725SJed Brown   if (Test_MatMatTr) {
222c20d7725SJed Brown     /* C = B*R^T */
2239566063dSJacob Faibussowitsch     PetscCall(PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij));
22420b3374bSStefano Zampini     if (seqaij) {
2259566063dSJacob Faibussowitsch       PetscCall(MatMatTransposeMult(B,R,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C));
2269566063dSJacob Faibussowitsch       PetscCall(MatSetOptionsPrefix(C,"ABt_")); /* enable '-ABt_' for matrix C */
2279566063dSJacob Faibussowitsch       PetscCall(MatGetInfo(C,MAT_GLOBAL_SUM,&info));
228c20d7725SJed Brown 
229c20d7725SJed Brown       /* Test MAT_REUSE_MATRIX - reuse symbolic C */
2309566063dSJacob Faibussowitsch       PetscCall(MatMatTransposeMult(B,R,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C));
231c20d7725SJed Brown 
232c20d7725SJed Brown       /* Check */
2339566063dSJacob Faibussowitsch       PetscCall(MatMatMult(B,P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C1));
2349566063dSJacob Faibussowitsch       PetscCall(MatNormDifference(C,C1,&norm));
235*08401ef6SPierre Jolivet       PetscCheck(norm <= PETSC_SMALL,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatMatTransposeMult() %g",(double)norm);
2369566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C1));
2379566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C));
238c20d7725SJed Brown     }
239c20d7725SJed Brown   }
240c20d7725SJed Brown 
241c20d7725SJed Brown   /* 4) Test MatPtAP() */
242c20d7725SJed Brown   /*-------------------*/
243c20d7725SJed Brown   if (Test_MatPtAP) {
2449566063dSJacob Faibussowitsch     PetscCall(MatDuplicate(A_save,MAT_COPY_VALUES,&A));
245c20d7725SJed Brown 
246c20d7725SJed Brown     /* (4.1) Test developer API */
2479566063dSJacob Faibussowitsch     PetscCall(MatProductCreate(A,P,NULL,&C));
2489566063dSJacob Faibussowitsch     PetscCall(MatSetOptionsPrefix(C,"PtAP_"));
2499566063dSJacob Faibussowitsch     PetscCall(MatProductSetType(C,MATPRODUCT_PtAP));
2509566063dSJacob Faibussowitsch     PetscCall(MatProductSetAlgorithm(C,MATPRODUCTALGORITHMDEFAULT));
2519566063dSJacob Faibussowitsch     PetscCall(MatProductSetFill(C,PETSC_DEFAULT));
2529566063dSJacob Faibussowitsch     PetscCall(MatProductSetFromOptions(C));
2539566063dSJacob Faibussowitsch     PetscCall(MatProductSymbolic(C));
2549566063dSJacob Faibussowitsch     PetscCall(MatProductNumeric(C));
2559566063dSJacob Faibussowitsch     PetscCall(MatPtAPMultEqual(A,P,C,mcheck,&flg));
25628b400f6SJacob Faibussowitsch     PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatProduct_PtAP");
2579566063dSJacob Faibussowitsch     PetscCall(MatProductNumeric(C)); /* reuse symbolic C */
258c20d7725SJed Brown 
2599566063dSJacob Faibussowitsch     PetscCall(MatPtAPMultEqual(A,P,C,mcheck,&flg));
26028b400f6SJacob Faibussowitsch     PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatProduct_PtAP");
2619566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
262c20d7725SJed Brown 
263c20d7725SJed Brown     /* (4.2) Test user driver */
2649566063dSJacob Faibussowitsch     PetscCall(MatPtAP(A,P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C));
265c20d7725SJed Brown 
266c20d7725SJed Brown     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
267c20d7725SJed Brown     alpha = 1.0;
268c20d7725SJed Brown     for (i=0; i<2; i++) {
269c20d7725SJed Brown       alpha -= 0.1;
2709566063dSJacob Faibussowitsch       PetscCall(MatScale(A,alpha));
2719566063dSJacob Faibussowitsch       PetscCall(MatPtAP(A,P,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C));
272c20d7725SJed Brown     }
2739566063dSJacob Faibussowitsch     PetscCall(MatPtAPMultEqual(A,P,C,mcheck,&flg));
27428b400f6SJacob Faibussowitsch     PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatPtAP");
275c20d7725SJed Brown 
276c20d7725SJed Brown     /* 5) Test MatRARt() */
277c20d7725SJed Brown     /* ----------------- */
278c20d7725SJed Brown     if (Test_MatRARt) {
279c20d7725SJed Brown       Mat RARt;
28067b3012eSStefano Zampini 
28167b3012eSStefano Zampini       /* (5.1) Test developer driver RARt = R*A*Rt */
2829566063dSJacob Faibussowitsch       PetscCall(MatProductCreate(A,R,NULL,&RARt));
2839566063dSJacob Faibussowitsch       PetscCall(MatSetOptionsPrefix(RARt,"RARt_"));
2849566063dSJacob Faibussowitsch       PetscCall(MatProductSetType(RARt,MATPRODUCT_RARt));
2859566063dSJacob Faibussowitsch       PetscCall(MatProductSetAlgorithm(RARt,MATPRODUCTALGORITHMDEFAULT));
2869566063dSJacob Faibussowitsch       PetscCall(MatProductSetFill(RARt,PETSC_DEFAULT));
2879566063dSJacob Faibussowitsch       PetscCall(MatProductSetFromOptions(RARt));
2889566063dSJacob Faibussowitsch       PetscCall(MatHasOperation(RARt,MATOP_PRODUCTSYMBOLIC,&flg));
289263f2b91SStefano Zampini       if (flg) {
2909566063dSJacob Faibussowitsch         PetscCall(MatProductSymbolic(RARt)); /* equivalent to MatSetUp() */
2919566063dSJacob Faibussowitsch         PetscCall(MatSetOption(RARt,MAT_USE_INODES,PETSC_FALSE)); /* illustrate how to call MatSetOption() */
2929566063dSJacob Faibussowitsch         PetscCall(MatProductNumeric(RARt));
2939566063dSJacob Faibussowitsch         PetscCall(MatProductNumeric(RARt)); /* test reuse symbolic RARt */
2949566063dSJacob Faibussowitsch         PetscCall(MatDestroy(&RARt));
29567b3012eSStefano Zampini 
29667b3012eSStefano Zampini         /* (2.2) Test user driver RARt = R*A*Rt */
2979566063dSJacob Faibussowitsch         PetscCall(MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&RARt));
2989566063dSJacob Faibussowitsch         PetscCall(MatRARt(A,R,MAT_REUSE_MATRIX,2.0,&RARt));
29967b3012eSStefano Zampini 
3009566063dSJacob Faibussowitsch         PetscCall(MatNormDifference(C,RARt,&norm));
301*08401ef6SPierre Jolivet         PetscCheck(norm <= PETSC_SMALL,PETSC_COMM_SELF,PETSC_ERR_PLIB,"|PtAP - RARt| = %g",(double)norm);
302263f2b91SStefano Zampini       } else {
3039566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatRARt not supported\n"));
304263f2b91SStefano Zampini       }
3059566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&RARt));
306c20d7725SJed Brown     }
307c20d7725SJed Brown 
3089566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&A));
3099566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
310c20d7725SJed Brown   }
311c20d7725SJed Brown 
312c20d7725SJed Brown   /* Destroy objects */
3139566063dSJacob Faibussowitsch   PetscCall(PetscRandomDestroy(&rdm));
3149566063dSJacob Faibussowitsch   PetscCall(PetscFree(idxn));
315c20d7725SJed Brown 
3169566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&A_save));
3179566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&B));
3189566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&P));
3199566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&R));
320c20d7725SJed Brown 
3219566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
322b122ec5aSJacob Faibussowitsch   return 0;
323c20d7725SJed Brown }
324c20d7725SJed Brown 
325c20d7725SJed Brown /*TEST
326c20d7725SJed Brown    test:
327c20d7725SJed Brown      suffix: 1
328dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
329c20d7725SJed Brown      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium
330c20d7725SJed Brown      output_file: output/ex62_1.out
331c20d7725SJed Brown 
332c20d7725SJed Brown    test:
333c20d7725SJed Brown      suffix: 2_ab_scalable
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 -matmatmult_via scalable -AtB_mat_product_algorithm outerproduct -mattransposematmult_via outerproduct
336c20d7725SJed Brown      output_file: output/ex62_1.out
337c20d7725SJed Brown 
338c20d7725SJed Brown    test:
339c20d7725SJed Brown      suffix: 3_ab_scalable_fast
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 scalable_fast -matmatmult_via scalable_fast -matmattransmult_via color
342c20d7725SJed Brown      output_file: output/ex62_1.out
343c20d7725SJed Brown 
344c20d7725SJed Brown    test:
345c20d7725SJed Brown      suffix: 4_ab_heap
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 heap -matmatmult_via heap -PtAP_mat_product_algorithm rap -matptap_via rap
348c20d7725SJed Brown      output_file: output/ex62_1.out
349c20d7725SJed Brown 
350c20d7725SJed Brown    test:
351c20d7725SJed Brown      suffix: 5_ab_btheap
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 btheap -matmatmult_via btheap -matrart_via r*art
354c20d7725SJed Brown      output_file: output/ex62_1.out
355c20d7725SJed Brown 
356c20d7725SJed Brown    test:
357c20d7725SJed Brown      suffix: 6_ab_llcondensed
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 llcondensed -matmatmult_via llcondensed -matrart_via coloring_rart
360c20d7725SJed Brown      output_file: output/ex62_1.out
361c20d7725SJed Brown 
362c20d7725SJed Brown    test:
363c20d7725SJed Brown      suffix: 7_ab_rowmerge
364dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3653e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm rowmerge -matmatmult_via rowmerge
366c20d7725SJed Brown      output_file: output/ex62_1.out
367c20d7725SJed Brown 
368c20d7725SJed Brown    test:
369c20d7725SJed Brown      suffix: 8_ab_hypre
370dfd57a17SPierre Jolivet      requires: hypre datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3713e662e0bSHong 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
372c20d7725SJed Brown      output_file: output/ex62_1.out
373c20d7725SJed Brown 
374c20d7725SJed Brown    test:
375263f2b91SStefano Zampini      suffix: hypre_medium
376263f2b91SStefano Zampini      nsize: {{1 3}}
377263f2b91SStefano Zampini      requires: hypre datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
378263f2b91SStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -A_mat_type hypre -B_mat_type hypre -test_rart 0
379263f2b91SStefano Zampini      output_file: output/ex62_hypre.out
380263f2b91SStefano Zampini 
381263f2b91SStefano Zampini    test:
382263f2b91SStefano Zampini      suffix: hypre_tiny
383263f2b91SStefano Zampini      nsize: {{1 3}}
384263f2b91SStefano Zampini      requires: hypre !complex double !defined(PETSC_USE_64BIT_INDICES)
385263f2b91SStefano 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
386263f2b91SStefano Zampini      output_file: output/ex62_hypre.out
387263f2b91SStefano Zampini 
388263f2b91SStefano Zampini    test:
389cec0a6c6SStefano Zampini      suffix: 9_mkl
39015df37b9SStefano Zampini      TODO: broken MatScale?
391b88df2e7SBarry Smith      requires: mkl_sparse datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
392cec0a6c6SStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -A_mat_type aijmkl -B_mat_type aijmkl
393cec0a6c6SStefano Zampini      output_file: output/ex62_1.out
394cec0a6c6SStefano Zampini 
395cec0a6c6SStefano Zampini    test:
396c20d7725SJed Brown      suffix: 10
397dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
398c20d7725SJed Brown      nsize: 3
399c20d7725SJed Brown      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium
400c20d7725SJed Brown      output_file: output/ex62_1.out
401c20d7725SJed Brown 
402c20d7725SJed Brown    test:
4033cea4f0aSStefano Zampini      suffix: 10_backend
404dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
4053cea4f0aSStefano Zampini      nsize: 3
4063e662e0bSHong 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
4073cea4f0aSStefano Zampini      output_file: output/ex62_1.out
4083cea4f0aSStefano Zampini 
4093cea4f0aSStefano Zampini    test:
410c20d7725SJed Brown      suffix: 11_ab_scalable
411dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
412c20d7725SJed Brown      nsize: 3
4133e662e0bSHong 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
414c20d7725SJed Brown      output_file: output/ex62_1.out
415c20d7725SJed Brown 
416c20d7725SJed Brown    test:
417c20d7725SJed Brown      suffix: 12_ab_seqmpi
418dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
419c20d7725SJed Brown      nsize: 3
4203e662e0bSHong 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
421c20d7725SJed Brown      output_file: output/ex62_1.out
422c20d7725SJed Brown 
423c20d7725SJed Brown    test:
424c20d7725SJed Brown      suffix: 13_ab_hypre
425dfd57a17SPierre Jolivet      requires: hypre datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
426c20d7725SJed Brown      nsize: 3
4273e662e0bSHong 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
428c20d7725SJed Brown      output_file: output/ex62_1.out
429c20d7725SJed Brown 
43020b3374bSStefano Zampini    test:
43120b3374bSStefano Zampini      suffix: 14_seqaij
432dfd57a17SPierre Jolivet      requires: !complex double !defined(PETSC_USE_64BIT_INDICES)
43320b3374bSStefano Zampini      args: -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:
43720b3374bSStefano Zampini      suffix: 14_seqaijcusparse
438dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4391a2c6b5cSJunchao 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
44020b3374bSStefano Zampini      output_file: output/ex62_1.out
44120b3374bSStefano Zampini 
44220b3374bSStefano Zampini    test:
44365e4b4d4SStefano Zampini      suffix: 14_seqaijcusparse_cpu
444dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4453e662e0bSHong 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
44665e4b4d4SStefano Zampini      output_file: output/ex62_1.out
44765e4b4d4SStefano Zampini 
44865e4b4d4SStefano Zampini    test:
4495cb69cabSStefano Zampini      suffix: 14_mpiaijcusparse_seq
4505cb69cabSStefano Zampini      nsize: 1
451dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4521a2c6b5cSJunchao 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
4535cb69cabSStefano Zampini      output_file: output/ex62_1.out
4545cb69cabSStefano Zampini 
4555cb69cabSStefano Zampini    test:
45665e4b4d4SStefano Zampini      suffix: 14_mpiaijcusparse_seq_cpu
45765e4b4d4SStefano Zampini      nsize: 1
458dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4594005f4d6SPierre 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
46065e4b4d4SStefano Zampini      output_file: output/ex62_1.out
46165e4b4d4SStefano Zampini 
46265e4b4d4SStefano Zampini    test:
4635cb69cabSStefano Zampini      suffix: 14_mpiaijcusparse
4645cb69cabSStefano Zampini      nsize: 3
465dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4661a2c6b5cSJunchao 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
4675cb69cabSStefano Zampini      output_file: output/ex62_1.out
4685cb69cabSStefano Zampini 
4695cb69cabSStefano Zampini    test:
47065e4b4d4SStefano Zampini      suffix: 14_mpiaijcusparse_cpu
47165e4b4d4SStefano Zampini      nsize: 3
472dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4734005f4d6SPierre 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
47465e4b4d4SStefano Zampini      output_file: output/ex62_1.out
47565e4b4d4SStefano Zampini 
47665e4b4d4SStefano Zampini    test:
477076ba34aSJunchao Zhang      nsize: {{1 3}}
478076ba34aSJunchao Zhang      suffix: 14_aijkokkos
4793078479eSJunchao Zhang      requires: !sycl kokkos_kernels !complex double !defined(PETSC_USE_64BIT_INDICES)
48020b3374bSStefano 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
48120b3374bSStefano Zampini      output_file: output/ex62_1.out
48220b3374bSStefano Zampini 
4835cb69cabSStefano Zampini    # these tests use matrices with many zero rows
48420b3374bSStefano Zampini    test:
48520b3374bSStefano Zampini      suffix: 15_seqaijcusparse
486dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES) datafilespath
4871a2c6b5cSJunchao Zhang      args: -A_mat_type aijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
48820b3374bSStefano Zampini      output_file: output/ex62_1.out
4895cb69cabSStefano Zampini 
4905cb69cabSStefano Zampini    test:
4915cb69cabSStefano Zampini      suffix: 15_mpiaijcusparse_seq
4925cb69cabSStefano Zampini      nsize: 1
493dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES) datafilespath
4941a2c6b5cSJunchao Zhang      args: -A_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
4955cb69cabSStefano Zampini      output_file: output/ex62_1.out
4965cb69cabSStefano Zampini 
4975cb69cabSStefano Zampini    test:
4985cb69cabSStefano Zampini      nsize: 3
4995cb69cabSStefano Zampini      suffix: 15_mpiaijcusparse
500dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES) datafilespath
5011a2c6b5cSJunchao Zhang      args: -A_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
5025cb69cabSStefano Zampini      output_file: output/ex62_1.out
5035cb69cabSStefano Zampini 
504c20d7725SJed Brown TEST*/
505