xref: /petsc/src/mat/tests/ex62.c (revision d0609ced746bc51b019815ca91d747429db24893)
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   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 
439566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc,&args,(char*)0,help));
449566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
459566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
46c20d7725SJed Brown 
47c20d7725SJed Brown   /*  Load the matrices A_save and B */
48*d0609cedSBarry Smith   PetscOptionsBegin(PETSC_COMM_WORLD,"","","");
499566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-test_rart","Test MatRARt","",Test_MatRARt,&Test_MatRARt,NULL));
509566063dSJacob Faibussowitsch   PetscCall(PetscOptionsInt("-PN","Number of columns of P","",PN,&PN,NULL));
519566063dSJacob Faibussowitsch   PetscCall(PetscOptionsInt("-mcheck","Number of matmult checks","",mcheck,&mcheck,NULL));
529566063dSJacob Faibussowitsch   PetscCall(PetscOptionsString("-fA","Path for matrix A","",file[0],file[0],sizeof(file[0]),&flg));
5328b400f6SJacob Faibussowitsch   PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name for matrix A with the -fA option.");
549566063dSJacob Faibussowitsch   PetscCall(PetscOptionsString("-fB","Path for matrix B","",file[1],file[1],sizeof(file[1]),&flg));
559566063dSJacob Faibussowitsch   PetscCall(PetscOptionsFList("-A_mat_type","Matrix type","MatSetType",MatList,deft,A_mattype,256,&flgA));
569566063dSJacob Faibussowitsch   PetscCall(PetscOptionsFList("-B_mat_type","Matrix type","MatSetType",MatList,deft,B_mattype,256,&flgB));
57*d0609cedSBarry Smith   PetscOptionsEnd();
58c20d7725SJed Brown 
599566063dSJacob Faibussowitsch   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&viewer));
609566063dSJacob Faibussowitsch   PetscCall(MatCreate(PETSC_COMM_WORLD,&A_save));
619566063dSJacob Faibussowitsch   PetscCall(MatLoad(A_save,viewer));
629566063dSJacob Faibussowitsch   PetscCall(PetscViewerDestroy(&viewer));
63c20d7725SJed Brown 
6420b3374bSStefano Zampini   if (flg) {
659566063dSJacob Faibussowitsch     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&viewer));
669566063dSJacob Faibussowitsch     PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
679566063dSJacob Faibussowitsch     PetscCall(MatLoad(B,viewer));
689566063dSJacob Faibussowitsch     PetscCall(PetscViewerDestroy(&viewer));
6920b3374bSStefano Zampini   } else {
709566063dSJacob Faibussowitsch     PetscCall(PetscObjectReference((PetscObject)A_save));
7120b3374bSStefano Zampini     B = A_save;
7220b3374bSStefano Zampini   }
7320b3374bSStefano Zampini 
7420b3374bSStefano Zampini   if (flgA) {
759566063dSJacob Faibussowitsch     PetscCall(MatConvert(A_save,A_mattype,MAT_INPLACE_MATRIX,&A_save));
7620b3374bSStefano Zampini   }
7720b3374bSStefano Zampini   if (flgB) {
789566063dSJacob Faibussowitsch     PetscCall(MatConvert(B,B_mattype,MAT_INPLACE_MATRIX,&B));
7920b3374bSStefano Zampini   }
809566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(A_save));
819566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(B));
82c20d7725SJed Brown 
839566063dSJacob Faibussowitsch   PetscCall(MatGetType(B,&mattype));
84c20d7725SJed Brown 
859566063dSJacob Faibussowitsch   PetscCall(PetscMalloc(nzp*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn));
86c20d7725SJed Brown   a    = (PetscScalar*)(idxn + nzp);
87c20d7725SJed Brown 
889566063dSJacob Faibussowitsch   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rdm));
899566063dSJacob Faibussowitsch   PetscCall(PetscRandomSetFromOptions(rdm));
90c20d7725SJed Brown 
91c20d7725SJed Brown   /* 1) MatMatMult() */
92c20d7725SJed Brown   /* ----------------*/
93c20d7725SJed Brown   if (Test_MatMatMult) {
949566063dSJacob Faibussowitsch     PetscCall(MatDuplicate(A_save,MAT_COPY_VALUES,&A));
95c20d7725SJed Brown 
96c20d7725SJed Brown     /* (1.1) Test developer API */
979566063dSJacob Faibussowitsch     PetscCall(MatProductCreate(A,B,NULL,&C));
989566063dSJacob Faibussowitsch     PetscCall(MatSetOptionsPrefix(C,"AB_"));
999566063dSJacob Faibussowitsch     PetscCall(MatProductSetType(C,MATPRODUCT_AB));
1009566063dSJacob Faibussowitsch     PetscCall(MatProductSetAlgorithm(C,MATPRODUCTALGORITHMDEFAULT));
1019566063dSJacob Faibussowitsch     PetscCall(MatProductSetFill(C,PETSC_DEFAULT));
1029566063dSJacob Faibussowitsch     PetscCall(MatProductSetFromOptions(C));
103910fa13fSStefano Zampini     /* we can inquire about MATOP_PRODUCTSYMBOLIC even if the destination matrix type has not been set yet */
1049566063dSJacob Faibussowitsch     PetscCall(MatHasOperation(C,MATOP_PRODUCTSYMBOLIC,&flg));
1059566063dSJacob Faibussowitsch     PetscCall(MatProductSymbolic(C));
1069566063dSJacob Faibussowitsch     PetscCall(MatProductNumeric(C));
1079566063dSJacob Faibussowitsch     PetscCall(MatMatMultEqual(A,B,C,mcheck,&flg));
10828b400f6SJacob Faibussowitsch     PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in C=A*B");
109c20d7725SJed Brown 
110c20d7725SJed Brown     /* Test reuse symbolic C */
111c20d7725SJed Brown     alpha = 0.9;
1129566063dSJacob Faibussowitsch     PetscCall(MatScale(A,alpha));
1139566063dSJacob Faibussowitsch     PetscCall(MatProductNumeric(C));
114c20d7725SJed Brown 
1159566063dSJacob Faibussowitsch     PetscCall(MatMatMultEqual(A,B,C,mcheck,&flg));
11628b400f6SJacob Faibussowitsch     PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in C=A*B");
1179566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
118c20d7725SJed Brown 
119c20d7725SJed Brown     /* (1.2) Test user driver */
1209566063dSJacob Faibussowitsch     PetscCall(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C));
121c20d7725SJed Brown 
122c20d7725SJed Brown     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
123c20d7725SJed Brown     alpha = 1.0;
124c20d7725SJed Brown     for (i=0; i<2; i++) {
125c20d7725SJed Brown       alpha -= 0.1;
1269566063dSJacob Faibussowitsch       PetscCall(MatScale(A,alpha));
1279566063dSJacob Faibussowitsch       PetscCall(MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C));
128c20d7725SJed Brown     }
1299566063dSJacob Faibussowitsch     PetscCall(MatMatMultEqual(A,B,C,mcheck,&flg));
13028b400f6SJacob Faibussowitsch     PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult()");
1319566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&A));
1324417c5e8SHong Zhang 
1334417c5e8SHong Zhang     /* Test MatProductClear() */
1349566063dSJacob Faibussowitsch     PetscCall(MatProductClear(C));
1359566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
136544a5e07SHong Zhang 
137544a5e07SHong Zhang     /* Test MatMatMult() for dense and aij matrices */
1389566063dSJacob Faibussowitsch     PetscCall(PetscObjectTypeCompareAny((PetscObject)A,&flg,MATSEQAIJ,MATMPIAIJ,""));
13920b3374bSStefano Zampini     if (flg) {
1409566063dSJacob Faibussowitsch       PetscCall(MatConvert(A_save,MATDENSE,MAT_INITIAL_MATRIX,&A));
1419566063dSJacob Faibussowitsch       PetscCall(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C));
1429566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C));
1439566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&A));
14420b3374bSStefano Zampini     }
145c20d7725SJed Brown   }
146c20d7725SJed Brown 
147c20d7725SJed Brown   /* Create P and R = P^T  */
148c20d7725SJed Brown   /* --------------------- */
1499566063dSJacob Faibussowitsch   PetscCall(MatGetSize(B,&PM,NULL));
15020b3374bSStefano Zampini   if (PN < 0) PN = PM/2;
1519566063dSJacob Faibussowitsch   PetscCall(MatCreate(PETSC_COMM_WORLD,&P));
1529566063dSJacob Faibussowitsch   PetscCall(MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,PM,PN));
1539566063dSJacob Faibussowitsch   PetscCall(MatSetType(P,MATAIJ));
1549566063dSJacob Faibussowitsch   PetscCall(MatSeqAIJSetPreallocation(P,nzp,NULL));
1559566063dSJacob Faibussowitsch   PetscCall(MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL));
1569566063dSJacob Faibussowitsch   PetscCall(MatGetOwnershipRange(P,&rstart,&rend));
157c20d7725SJed Brown   for (i=0; i<nzp; i++) {
1589566063dSJacob Faibussowitsch     PetscCall(PetscRandomGetValue(rdm,&a[i]));
159c20d7725SJed Brown   }
160c20d7725SJed Brown   for (i=rstart; i<rend; i++) {
161c20d7725SJed Brown     for (j=0; j<nzp; j++) {
1629566063dSJacob Faibussowitsch       PetscCall(PetscRandomGetValue(rdm,&rval));
163c20d7725SJed Brown       idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
164c20d7725SJed Brown     }
1659566063dSJacob Faibussowitsch     PetscCall(MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES));
166c20d7725SJed Brown   }
1679566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY));
1689566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY));
169c20d7725SJed Brown 
1709566063dSJacob Faibussowitsch   PetscCall(MatTranspose(P,MAT_INITIAL_MATRIX,&R));
1719566063dSJacob Faibussowitsch   PetscCall(MatConvert(P,mattype,MAT_INPLACE_MATRIX,&P));
1729566063dSJacob Faibussowitsch   PetscCall(MatConvert(R,mattype,MAT_INPLACE_MATRIX,&R));
1739566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(P));
1749566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(R));
175c20d7725SJed Brown 
176c20d7725SJed Brown   /* 2) MatTransposeMatMult() */
177c20d7725SJed Brown   /* ------------------------ */
178c20d7725SJed Brown   if (Test_MatTrMat) {
179c20d7725SJed Brown     /* (2.1) Test developer driver C = P^T*B */
1809566063dSJacob Faibussowitsch     PetscCall(MatProductCreate(P,B,NULL,&C));
1819566063dSJacob Faibussowitsch     PetscCall(MatSetOptionsPrefix(C,"AtB_"));
1829566063dSJacob Faibussowitsch     PetscCall(MatProductSetType(C,MATPRODUCT_AtB));
1839566063dSJacob Faibussowitsch     PetscCall(MatProductSetAlgorithm(C,MATPRODUCTALGORITHMDEFAULT));
1849566063dSJacob Faibussowitsch     PetscCall(MatProductSetFill(C,PETSC_DEFAULT));
1859566063dSJacob Faibussowitsch     PetscCall(MatProductSetFromOptions(C));
1869566063dSJacob Faibussowitsch     PetscCall(MatHasOperation(C,MATOP_PRODUCTSYMBOLIC,&flg));
187263f2b91SStefano Zampini     if (flg) { /* run tests if supported */
1889566063dSJacob Faibussowitsch       PetscCall(MatProductSymbolic(C)); /* equivalent to MatSetUp() */
1899566063dSJacob Faibussowitsch       PetscCall(MatSetOption(C,MAT_USE_INODES,PETSC_FALSE)); /* illustrate how to call MatSetOption() */
1909566063dSJacob Faibussowitsch       PetscCall(MatProductNumeric(C));
1919566063dSJacob Faibussowitsch       PetscCall(MatProductNumeric(C)); /* test reuse symbolic C */
19267b3012eSStefano Zampini 
1939566063dSJacob Faibussowitsch       PetscCall(MatTransposeMatMultEqual(P,B,C,mcheck,&flg));
19428b400f6SJacob Faibussowitsch       PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error: developer driver C = P^T*B");
1959566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C));
196c20d7725SJed Brown 
197c20d7725SJed Brown       /* (2.2) Test user driver C = P^T*B */
1989566063dSJacob Faibussowitsch       PetscCall(MatTransposeMatMult(P,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C));
1999566063dSJacob Faibussowitsch       PetscCall(MatTransposeMatMult(P,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C));
2009566063dSJacob Faibussowitsch       PetscCall(MatGetInfo(C,MAT_GLOBAL_SUM,&info));
2019566063dSJacob Faibussowitsch       PetscCall(MatProductClear(C));
202c20d7725SJed Brown 
203c20d7725SJed Brown       /* Compare P^T*B and R*B */
2049566063dSJacob Faibussowitsch       PetscCall(MatMatMult(R,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C1));
2059566063dSJacob Faibussowitsch       PetscCall(MatNormDifference(C,C1,&norm));
20608401ef6SPierre Jolivet       PetscCheck(norm <= PETSC_SMALL,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatTransposeMatMult(): %g",(double)norm);
2079566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C1));
208c20d7725SJed Brown 
209c20d7725SJed Brown       /* Test MatDuplicate() of C=P^T*B */
2109566063dSJacob Faibussowitsch       PetscCall(MatDuplicate(C,MAT_COPY_VALUES,&C1));
2119566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C1));
212263f2b91SStefano Zampini     } else {
2139566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatTransposeMatMult not supported\n"));
214263f2b91SStefano Zampini     }
2159566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
216c20d7725SJed Brown   }
217c20d7725SJed Brown 
21867b3012eSStefano Zampini   /* 3) MatMatTransposeMult() */
219c20d7725SJed Brown   /* ------------------------ */
220c20d7725SJed Brown   if (Test_MatMatTr) {
221c20d7725SJed Brown     /* C = B*R^T */
2229566063dSJacob Faibussowitsch     PetscCall(PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij));
22320b3374bSStefano Zampini     if (seqaij) {
2249566063dSJacob Faibussowitsch       PetscCall(MatMatTransposeMult(B,R,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C));
2259566063dSJacob Faibussowitsch       PetscCall(MatSetOptionsPrefix(C,"ABt_")); /* enable '-ABt_' for matrix C */
2269566063dSJacob Faibussowitsch       PetscCall(MatGetInfo(C,MAT_GLOBAL_SUM,&info));
227c20d7725SJed Brown 
228c20d7725SJed Brown       /* Test MAT_REUSE_MATRIX - reuse symbolic C */
2299566063dSJacob Faibussowitsch       PetscCall(MatMatTransposeMult(B,R,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C));
230c20d7725SJed Brown 
231c20d7725SJed Brown       /* Check */
2329566063dSJacob Faibussowitsch       PetscCall(MatMatMult(B,P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C1));
2339566063dSJacob Faibussowitsch       PetscCall(MatNormDifference(C,C1,&norm));
23408401ef6SPierre Jolivet       PetscCheck(norm <= PETSC_SMALL,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatMatTransposeMult() %g",(double)norm);
2359566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C1));
2369566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&C));
237c20d7725SJed Brown     }
238c20d7725SJed Brown   }
239c20d7725SJed Brown 
240c20d7725SJed Brown   /* 4) Test MatPtAP() */
241c20d7725SJed Brown   /*-------------------*/
242c20d7725SJed Brown   if (Test_MatPtAP) {
2439566063dSJacob Faibussowitsch     PetscCall(MatDuplicate(A_save,MAT_COPY_VALUES,&A));
244c20d7725SJed Brown 
245c20d7725SJed Brown     /* (4.1) Test developer API */
2469566063dSJacob Faibussowitsch     PetscCall(MatProductCreate(A,P,NULL,&C));
2479566063dSJacob Faibussowitsch     PetscCall(MatSetOptionsPrefix(C,"PtAP_"));
2489566063dSJacob Faibussowitsch     PetscCall(MatProductSetType(C,MATPRODUCT_PtAP));
2499566063dSJacob Faibussowitsch     PetscCall(MatProductSetAlgorithm(C,MATPRODUCTALGORITHMDEFAULT));
2509566063dSJacob Faibussowitsch     PetscCall(MatProductSetFill(C,PETSC_DEFAULT));
2519566063dSJacob Faibussowitsch     PetscCall(MatProductSetFromOptions(C));
2529566063dSJacob Faibussowitsch     PetscCall(MatProductSymbolic(C));
2539566063dSJacob Faibussowitsch     PetscCall(MatProductNumeric(C));
2549566063dSJacob Faibussowitsch     PetscCall(MatPtAPMultEqual(A,P,C,mcheck,&flg));
25528b400f6SJacob Faibussowitsch     PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatProduct_PtAP");
2569566063dSJacob Faibussowitsch     PetscCall(MatProductNumeric(C)); /* reuse symbolic C */
257c20d7725SJed Brown 
2589566063dSJacob Faibussowitsch     PetscCall(MatPtAPMultEqual(A,P,C,mcheck,&flg));
25928b400f6SJacob Faibussowitsch     PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatProduct_PtAP");
2609566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
261c20d7725SJed Brown 
262c20d7725SJed Brown     /* (4.2) Test user driver */
2639566063dSJacob Faibussowitsch     PetscCall(MatPtAP(A,P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C));
264c20d7725SJed Brown 
265c20d7725SJed Brown     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
266c20d7725SJed Brown     alpha = 1.0;
267c20d7725SJed Brown     for (i=0; i<2; i++) {
268c20d7725SJed Brown       alpha -= 0.1;
2699566063dSJacob Faibussowitsch       PetscCall(MatScale(A,alpha));
2709566063dSJacob Faibussowitsch       PetscCall(MatPtAP(A,P,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C));
271c20d7725SJed Brown     }
2729566063dSJacob Faibussowitsch     PetscCall(MatPtAPMultEqual(A,P,C,mcheck,&flg));
27328b400f6SJacob Faibussowitsch     PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatPtAP");
274c20d7725SJed Brown 
275c20d7725SJed Brown     /* 5) Test MatRARt() */
276c20d7725SJed Brown     /* ----------------- */
277c20d7725SJed Brown     if (Test_MatRARt) {
278c20d7725SJed Brown       Mat RARt;
27967b3012eSStefano Zampini 
28067b3012eSStefano Zampini       /* (5.1) Test developer driver RARt = R*A*Rt */
2819566063dSJacob Faibussowitsch       PetscCall(MatProductCreate(A,R,NULL,&RARt));
2829566063dSJacob Faibussowitsch       PetscCall(MatSetOptionsPrefix(RARt,"RARt_"));
2839566063dSJacob Faibussowitsch       PetscCall(MatProductSetType(RARt,MATPRODUCT_RARt));
2849566063dSJacob Faibussowitsch       PetscCall(MatProductSetAlgorithm(RARt,MATPRODUCTALGORITHMDEFAULT));
2859566063dSJacob Faibussowitsch       PetscCall(MatProductSetFill(RARt,PETSC_DEFAULT));
2869566063dSJacob Faibussowitsch       PetscCall(MatProductSetFromOptions(RARt));
2879566063dSJacob Faibussowitsch       PetscCall(MatHasOperation(RARt,MATOP_PRODUCTSYMBOLIC,&flg));
288263f2b91SStefano Zampini       if (flg) {
2899566063dSJacob Faibussowitsch         PetscCall(MatProductSymbolic(RARt)); /* equivalent to MatSetUp() */
2909566063dSJacob Faibussowitsch         PetscCall(MatSetOption(RARt,MAT_USE_INODES,PETSC_FALSE)); /* illustrate how to call MatSetOption() */
2919566063dSJacob Faibussowitsch         PetscCall(MatProductNumeric(RARt));
2929566063dSJacob Faibussowitsch         PetscCall(MatProductNumeric(RARt)); /* test reuse symbolic RARt */
2939566063dSJacob Faibussowitsch         PetscCall(MatDestroy(&RARt));
29467b3012eSStefano Zampini 
29567b3012eSStefano Zampini         /* (2.2) Test user driver RARt = R*A*Rt */
2969566063dSJacob Faibussowitsch         PetscCall(MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&RARt));
2979566063dSJacob Faibussowitsch         PetscCall(MatRARt(A,R,MAT_REUSE_MATRIX,2.0,&RARt));
29867b3012eSStefano Zampini 
2999566063dSJacob Faibussowitsch         PetscCall(MatNormDifference(C,RARt,&norm));
30008401ef6SPierre Jolivet         PetscCheck(norm <= PETSC_SMALL,PETSC_COMM_SELF,PETSC_ERR_PLIB,"|PtAP - RARt| = %g",(double)norm);
301263f2b91SStefano Zampini       } else {
3029566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatRARt not supported\n"));
303263f2b91SStefano Zampini       }
3049566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&RARt));
305c20d7725SJed Brown     }
306c20d7725SJed Brown 
3079566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&A));
3089566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&C));
309c20d7725SJed Brown   }
310c20d7725SJed Brown 
311c20d7725SJed Brown   /* Destroy objects */
3129566063dSJacob Faibussowitsch   PetscCall(PetscRandomDestroy(&rdm));
3139566063dSJacob Faibussowitsch   PetscCall(PetscFree(idxn));
314c20d7725SJed Brown 
3159566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&A_save));
3169566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&B));
3179566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&P));
3189566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&R));
319c20d7725SJed Brown 
3209566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
321b122ec5aSJacob Faibussowitsch   return 0;
322c20d7725SJed Brown }
323c20d7725SJed Brown 
324c20d7725SJed Brown /*TEST
325c20d7725SJed Brown    test:
326c20d7725SJed Brown      suffix: 1
327dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
328c20d7725SJed Brown      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium
329c20d7725SJed Brown      output_file: output/ex62_1.out
330c20d7725SJed Brown 
331c20d7725SJed Brown    test:
332c20d7725SJed Brown      suffix: 2_ab_scalable
333dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3343e662e0bSHong 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
335c20d7725SJed Brown      output_file: output/ex62_1.out
336c20d7725SJed Brown 
337c20d7725SJed Brown    test:
338c20d7725SJed Brown      suffix: 3_ab_scalable_fast
339dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3403e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm scalable_fast -matmatmult_via scalable_fast -matmattransmult_via color
341c20d7725SJed Brown      output_file: output/ex62_1.out
342c20d7725SJed Brown 
343c20d7725SJed Brown    test:
344c20d7725SJed Brown      suffix: 4_ab_heap
345dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3463e662e0bSHong 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
347c20d7725SJed Brown      output_file: output/ex62_1.out
348c20d7725SJed Brown 
349c20d7725SJed Brown    test:
350c20d7725SJed Brown      suffix: 5_ab_btheap
351dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3523e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm btheap -matmatmult_via btheap -matrart_via r*art
353c20d7725SJed Brown      output_file: output/ex62_1.out
354c20d7725SJed Brown 
355c20d7725SJed Brown    test:
356c20d7725SJed Brown      suffix: 6_ab_llcondensed
357dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3583e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm llcondensed -matmatmult_via llcondensed -matrart_via coloring_rart
359c20d7725SJed Brown      output_file: output/ex62_1.out
360c20d7725SJed Brown 
361c20d7725SJed Brown    test:
362c20d7725SJed Brown      suffix: 7_ab_rowmerge
363dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3643e662e0bSHong Zhang      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_mat_product_algorithm rowmerge -matmatmult_via rowmerge
365c20d7725SJed Brown      output_file: output/ex62_1.out
366c20d7725SJed Brown 
367c20d7725SJed Brown    test:
368c20d7725SJed Brown      suffix: 8_ab_hypre
369dfd57a17SPierre Jolivet      requires: hypre datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3703e662e0bSHong 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
371c20d7725SJed Brown      output_file: output/ex62_1.out
372c20d7725SJed Brown 
373c20d7725SJed Brown    test:
374263f2b91SStefano Zampini      suffix: hypre_medium
375263f2b91SStefano Zampini      nsize: {{1 3}}
376263f2b91SStefano Zampini      requires: hypre datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
377263f2b91SStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -A_mat_type hypre -B_mat_type hypre -test_rart 0
378263f2b91SStefano Zampini      output_file: output/ex62_hypre.out
379263f2b91SStefano Zampini 
380263f2b91SStefano Zampini    test:
381263f2b91SStefano Zampini      suffix: hypre_tiny
382263f2b91SStefano Zampini      nsize: {{1 3}}
383263f2b91SStefano Zampini      requires: hypre !complex double !defined(PETSC_USE_64BIT_INDICES)
384263f2b91SStefano 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
385263f2b91SStefano Zampini      output_file: output/ex62_hypre.out
386263f2b91SStefano Zampini 
387263f2b91SStefano Zampini    test:
388cec0a6c6SStefano Zampini      suffix: 9_mkl
38915df37b9SStefano Zampini      TODO: broken MatScale?
390b88df2e7SBarry Smith      requires: mkl_sparse datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
391cec0a6c6SStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -A_mat_type aijmkl -B_mat_type aijmkl
392cec0a6c6SStefano Zampini      output_file: output/ex62_1.out
393cec0a6c6SStefano Zampini 
394cec0a6c6SStefano Zampini    test:
395c20d7725SJed Brown      suffix: 10
396dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
397c20d7725SJed Brown      nsize: 3
398c20d7725SJed Brown      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium
399c20d7725SJed Brown      output_file: output/ex62_1.out
400c20d7725SJed Brown 
401c20d7725SJed Brown    test:
4023cea4f0aSStefano Zampini      suffix: 10_backend
403dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
4043cea4f0aSStefano Zampini      nsize: 3
4053e662e0bSHong 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
4063cea4f0aSStefano Zampini      output_file: output/ex62_1.out
4073cea4f0aSStefano Zampini 
4083cea4f0aSStefano Zampini    test:
409c20d7725SJed Brown      suffix: 11_ab_scalable
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 scalable -matmatmult_via scalable -AtB_mat_product_algorithm scalable -mattransposematmult_via scalable
413c20d7725SJed Brown      output_file: output/ex62_1.out
414c20d7725SJed Brown 
415c20d7725SJed Brown    test:
416c20d7725SJed Brown      suffix: 12_ab_seqmpi
417dfd57a17SPierre Jolivet      requires: 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 seqmpi -matmatmult_via seqmpi -AtB_mat_product_algorithm at*b -mattransposematmult_via at*b
420c20d7725SJed Brown      output_file: output/ex62_1.out
421c20d7725SJed Brown 
422c20d7725SJed Brown    test:
423c20d7725SJed Brown      suffix: 13_ab_hypre
424dfd57a17SPierre Jolivet      requires: hypre datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
425c20d7725SJed Brown      nsize: 3
4263e662e0bSHong 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
427c20d7725SJed Brown      output_file: output/ex62_1.out
428c20d7725SJed Brown 
42920b3374bSStefano Zampini    test:
43020b3374bSStefano Zampini      suffix: 14_seqaij
431dfd57a17SPierre Jolivet      requires: !complex double !defined(PETSC_USE_64BIT_INDICES)
43220b3374bSStefano Zampini      args: -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
43320b3374bSStefano Zampini      output_file: output/ex62_1.out
43420b3374bSStefano Zampini 
43520b3374bSStefano Zampini    test:
43620b3374bSStefano Zampini      suffix: 14_seqaijcusparse
437dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4381a2c6b5cSJunchao 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
43920b3374bSStefano Zampini      output_file: output/ex62_1.out
44020b3374bSStefano Zampini 
44120b3374bSStefano Zampini    test:
44265e4b4d4SStefano Zampini      suffix: 14_seqaijcusparse_cpu
443dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4443e662e0bSHong 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
44565e4b4d4SStefano Zampini      output_file: output/ex62_1.out
44665e4b4d4SStefano Zampini 
44765e4b4d4SStefano Zampini    test:
4485cb69cabSStefano Zampini      suffix: 14_mpiaijcusparse_seq
4495cb69cabSStefano Zampini      nsize: 1
450dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4511a2c6b5cSJunchao 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
4525cb69cabSStefano Zampini      output_file: output/ex62_1.out
4535cb69cabSStefano Zampini 
4545cb69cabSStefano Zampini    test:
45565e4b4d4SStefano Zampini      suffix: 14_mpiaijcusparse_seq_cpu
45665e4b4d4SStefano Zampini      nsize: 1
457dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4584005f4d6SPierre 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
45965e4b4d4SStefano Zampini      output_file: output/ex62_1.out
46065e4b4d4SStefano Zampini 
46165e4b4d4SStefano Zampini    test:
4625cb69cabSStefano Zampini      suffix: 14_mpiaijcusparse
4635cb69cabSStefano Zampini      nsize: 3
464dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4651a2c6b5cSJunchao 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
4665cb69cabSStefano Zampini      output_file: output/ex62_1.out
4675cb69cabSStefano Zampini 
4685cb69cabSStefano Zampini    test:
46965e4b4d4SStefano Zampini      suffix: 14_mpiaijcusparse_cpu
47065e4b4d4SStefano Zampini      nsize: 3
471dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4724005f4d6SPierre 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
47365e4b4d4SStefano Zampini      output_file: output/ex62_1.out
47465e4b4d4SStefano Zampini 
47565e4b4d4SStefano Zampini    test:
476076ba34aSJunchao Zhang      nsize: {{1 3}}
477076ba34aSJunchao Zhang      suffix: 14_aijkokkos
4783078479eSJunchao Zhang      requires: !sycl kokkos_kernels !complex double !defined(PETSC_USE_64BIT_INDICES)
47920b3374bSStefano 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
48020b3374bSStefano Zampini      output_file: output/ex62_1.out
48120b3374bSStefano Zampini 
4825cb69cabSStefano Zampini    # these tests use matrices with many zero rows
48320b3374bSStefano Zampini    test:
48420b3374bSStefano Zampini      suffix: 15_seqaijcusparse
485dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES) datafilespath
4861a2c6b5cSJunchao Zhang      args: -A_mat_type aijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
48720b3374bSStefano Zampini      output_file: output/ex62_1.out
4885cb69cabSStefano Zampini 
4895cb69cabSStefano Zampini    test:
4905cb69cabSStefano Zampini      suffix: 15_mpiaijcusparse_seq
4915cb69cabSStefano Zampini      nsize: 1
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 
4965cb69cabSStefano Zampini    test:
4975cb69cabSStefano Zampini      nsize: 3
4985cb69cabSStefano Zampini      suffix: 15_mpiaijcusparse
499dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES) datafilespath
5001a2c6b5cSJunchao Zhang      args: -A_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
5015cb69cabSStefano Zampini      output_file: output/ex62_1.out
5025cb69cabSStefano Zampini 
503c20d7725SJed Brown TEST*/
504