xref: /petsc/src/mat/tests/ex62.c (revision dfd57a172ac9fa6c7b5fe6de6ab5df85cefc2996)
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   PetscErrorCode ierr;
19c20d7725SJed Brown 
20c20d7725SJed Brown   PetscFunctionBegin;
21c20d7725SJed Brown   ierr = MatAXPY(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
22c20d7725SJed Brown   ierr = MatNorm(B,NORM_FROBENIUS,norm);CHKERRQ(ierr);
23c20d7725SJed Brown   PetscFunctionReturn(0);
24c20d7725SJed Brown }
25c20d7725SJed Brown 
26c20d7725SJed Brown int main(int argc,char **args)
27c20d7725SJed Brown {
28c20d7725SJed Brown   Mat            A,A_save,B,C,P,C1,R;
29c20d7725SJed Brown   PetscViewer    viewer;
30c20d7725SJed Brown   PetscErrorCode ierr;
31c20d7725SJed Brown   PetscMPIInt    size,rank;
3220b3374bSStefano Zampini   PetscInt       i,j,*idxn,PM,PN = PETSC_DECIDE,rstart,rend;
33c20d7725SJed Brown   PetscReal      norm;
34c20d7725SJed Brown   PetscRandom    rdm;
3520b3374bSStefano Zampini   char           file[2][PETSC_MAX_PATH_LEN] = { "", ""};
36c20d7725SJed Brown   PetscScalar    *a,rval,alpha;
37c20d7725SJed Brown   PetscBool      Test_MatMatMult=PETSC_TRUE,Test_MatTrMat=PETSC_TRUE,Test_MatMatTr=PETSC_TRUE;
3820b3374bSStefano Zampini   PetscBool      Test_MatPtAP=PETSC_TRUE,Test_MatRARt=PETSC_TRUE,flg,seqaij,flgA,flgB;
39c20d7725SJed Brown   MatInfo        info;
4020b3374bSStefano Zampini   PetscInt       nzp  = 5; /* num of nonzeros in each row of P */
41c20d7725SJed Brown   MatType        mattype;
4220b3374bSStefano Zampini   const char     *deft = MATAIJ;
4320b3374bSStefano Zampini   char           A_mattype[256], B_mattype[256];
4420b3374bSStefano Zampini   PetscInt       mcheck = 10;
45c20d7725SJed Brown 
46c20d7725SJed Brown   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
47ffc4695bSBarry Smith   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr);
48ffc4695bSBarry Smith   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr);
49c20d7725SJed Brown 
50c20d7725SJed Brown   /*  Load the matrices A_save and B */
5120b3374bSStefano Zampini   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","","");CHKERRQ(ierr);
5220b3374bSStefano Zampini   ierr = PetscOptionsInt("-PN","Number of columns of P","",PN,&PN,NULL);CHKERRQ(ierr);
5320b3374bSStefano Zampini   ierr = PetscOptionsInt("-mcheck","Number of matmult checks","",mcheck,&mcheck,NULL);CHKERRQ(ierr);
5420b3374bSStefano Zampini   ierr = PetscOptionsString("-fA","Path for matrix A","",file[0],file[0],sizeof(file[0]),&flg);CHKERRQ(ierr);
55d8185827SBarry Smith   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name for matrix A with the -fA option.");
5620b3374bSStefano Zampini   ierr = PetscOptionsString("-fB","Path for matrix B","",file[1],file[1],sizeof(file[1]),&flg);CHKERRQ(ierr);
5720b3374bSStefano Zampini   ierr = PetscOptionsFList("-A_mat_type","Matrix type","MatSetType",MatList,deft,A_mattype,256,&flgA);CHKERRQ(ierr);
5820b3374bSStefano Zampini   ierr = PetscOptionsFList("-B_mat_type","Matrix type","MatSetType",MatList,deft,B_mattype,256,&flgB);CHKERRQ(ierr);
5920b3374bSStefano Zampini   ierr = PetscOptionsEnd();CHKERRQ(ierr);
60c20d7725SJed Brown 
61c20d7725SJed Brown   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&viewer);CHKERRQ(ierr);
62c20d7725SJed Brown   ierr = MatCreate(PETSC_COMM_WORLD,&A_save);CHKERRQ(ierr);
63c20d7725SJed Brown   ierr = MatLoad(A_save,viewer);CHKERRQ(ierr);
64c20d7725SJed Brown   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
65c20d7725SJed Brown 
6620b3374bSStefano Zampini   if (flg) {
67c20d7725SJed Brown     ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&viewer);CHKERRQ(ierr);
68c20d7725SJed Brown     ierr = MatCreate(PETSC_COMM_WORLD,&B);CHKERRQ(ierr);
69c20d7725SJed Brown     ierr = MatLoad(B,viewer);CHKERRQ(ierr);
70c20d7725SJed Brown     ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
7120b3374bSStefano Zampini   } else {
7220b3374bSStefano Zampini     ierr = PetscObjectReference((PetscObject)A_save);CHKERRQ(ierr);
7320b3374bSStefano Zampini     B = A_save;
7420b3374bSStefano Zampini   }
7520b3374bSStefano Zampini 
7620b3374bSStefano Zampini   if (flgA) {
7720b3374bSStefano Zampini     ierr = MatConvert(A_save,A_mattype,MAT_INPLACE_MATRIX,&A_save);CHKERRQ(ierr);
7820b3374bSStefano Zampini   }
7920b3374bSStefano Zampini   if (flgB) {
8020b3374bSStefano Zampini     ierr = MatConvert(B,B_mattype,MAT_INPLACE_MATRIX,&B);CHKERRQ(ierr);
8120b3374bSStefano Zampini   }
8220b3374bSStefano Zampini   ierr = MatSetFromOptions(A_save);CHKERRQ(ierr);
8320b3374bSStefano Zampini   ierr = MatSetFromOptions(B);CHKERRQ(ierr);
84c20d7725SJed Brown 
85c20d7725SJed Brown   ierr = MatGetType(B,&mattype);CHKERRQ(ierr);
86c20d7725SJed Brown 
8720b3374bSStefano Zampini   ierr = PetscMalloc(nzp*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn);CHKERRQ(ierr);
88c20d7725SJed Brown   a    = (PetscScalar*)(idxn + nzp);
89c20d7725SJed Brown 
90c20d7725SJed Brown   ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rdm);CHKERRQ(ierr);
91c20d7725SJed Brown   ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr);
92c20d7725SJed Brown 
93c20d7725SJed Brown   /* 1) MatMatMult() */
94c20d7725SJed Brown   /* ----------------*/
95c20d7725SJed Brown   if (Test_MatMatMult) {
96c20d7725SJed Brown     ierr = MatDuplicate(A_save,MAT_COPY_VALUES,&A);CHKERRQ(ierr);
97c20d7725SJed Brown 
98c20d7725SJed Brown     /* (1.1) Test developer API */
99c20d7725SJed Brown     ierr = MatProductCreate(A,B,NULL,&C);CHKERRQ(ierr);
10067b3012eSStefano Zampini     ierr = MatSetOptionsPrefix(C,"AB_");CHKERRQ(ierr);
101c20d7725SJed Brown     ierr = MatProductSetType(C,MATPRODUCT_AB);CHKERRQ(ierr);
102910fa13fSStefano Zampini     ierr = MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr);
103c20d7725SJed Brown     ierr = MatProductSetFill(C,PETSC_DEFAULT);CHKERRQ(ierr);
104c20d7725SJed Brown     ierr = MatProductSetFromOptions(C);CHKERRQ(ierr);
105910fa13fSStefano Zampini     /* we can inquire about MATOP_PRODUCTSYMBOLIC even if the destination matrix type has not been set yet */
106910fa13fSStefano Zampini     ierr = MatHasOperation(C,MATOP_PRODUCTSYMBOLIC,&flg);CHKERRQ(ierr);
107c20d7725SJed Brown     ierr = MatProductSymbolic(C);CHKERRQ(ierr);
108c20d7725SJed Brown     ierr = MatProductNumeric(C);CHKERRQ(ierr);
10920b3374bSStefano Zampini     ierr = MatMatMultEqual(A,B,C,mcheck,&flg);CHKERRQ(ierr);
11020b3374bSStefano Zampini     if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in C=A*B");
111c20d7725SJed Brown 
112c20d7725SJed Brown     /* Test reuse symbolic C */
113c20d7725SJed Brown     alpha = 0.9;
114c20d7725SJed Brown     ierr = MatScale(A,alpha);CHKERRQ(ierr);
115c20d7725SJed Brown     ierr = MatProductNumeric(C);CHKERRQ(ierr);
116c20d7725SJed Brown 
11720b3374bSStefano Zampini     ierr = MatMatMultEqual(A,B,C,mcheck,&flg);CHKERRQ(ierr);
118c20d7725SJed Brown     if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in C=A*B");
119c20d7725SJed Brown     ierr = MatDestroy(&C);CHKERRQ(ierr);
120c20d7725SJed Brown 
121c20d7725SJed Brown     /* (1.2) Test user driver */
122c20d7725SJed Brown     ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
123c20d7725SJed Brown 
124c20d7725SJed Brown     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
125c20d7725SJed Brown     alpha = 1.0;
126c20d7725SJed Brown     for (i=0; i<2; i++) {
127c20d7725SJed Brown       alpha -= 0.1;
128c20d7725SJed Brown       ierr   = MatScale(A,alpha);CHKERRQ(ierr);
129c20d7725SJed Brown       ierr   = MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
130c20d7725SJed Brown     }
13120b3374bSStefano Zampini     ierr = MatMatMultEqual(A,B,C,mcheck,&flg);CHKERRQ(ierr);
132c20d7725SJed Brown     if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult()");
133c20d7725SJed Brown     ierr = MatDestroy(&A);CHKERRQ(ierr);
1344417c5e8SHong Zhang 
1354417c5e8SHong Zhang     /* Test MatProductClear() */
1364417c5e8SHong Zhang     ierr = MatProductClear(C);CHKERRQ(ierr);
137c20d7725SJed Brown     ierr = MatDestroy(&C);CHKERRQ(ierr);
138544a5e07SHong Zhang 
139544a5e07SHong Zhang     /* Test MatMatMult() for dense and aij matrices */
14020b3374bSStefano Zampini     ierr = PetscObjectTypeCompareAny((PetscObject)A,&flg,MATSEQAIJ,MATMPIAIJ,"");CHKERRQ(ierr);
14120b3374bSStefano Zampini     if (flg) {
142544a5e07SHong Zhang       ierr = MatConvert(A_save,MATDENSE,MAT_INITIAL_MATRIX,&A);CHKERRQ(ierr);
143544a5e07SHong Zhang       ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
144544a5e07SHong Zhang       ierr = MatDestroy(&C);CHKERRQ(ierr);
145544a5e07SHong Zhang       ierr = MatDestroy(&A);CHKERRQ(ierr);
14620b3374bSStefano Zampini     }
147c20d7725SJed Brown   }
148c20d7725SJed Brown 
149c20d7725SJed Brown   /* Create P and R = P^T  */
150c20d7725SJed Brown   /* --------------------- */
15120b3374bSStefano Zampini   ierr = MatGetSize(B,&PM,NULL);CHKERRQ(ierr);
15220b3374bSStefano Zampini   if (PN < 0) PN = PM/2;
153c20d7725SJed Brown   ierr = MatCreate(PETSC_COMM_WORLD,&P);CHKERRQ(ierr);
15420b3374bSStefano Zampini   ierr = MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,PM,PN);CHKERRQ(ierr);
155c20d7725SJed Brown   ierr = MatSetType(P,mattype);CHKERRQ(ierr);
156c20d7725SJed Brown   ierr = MatSeqAIJSetPreallocation(P,nzp,NULL);CHKERRQ(ierr);
157c20d7725SJed Brown   ierr = MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL);CHKERRQ(ierr);
158c20d7725SJed Brown   ierr = MatGetOwnershipRange(P,&rstart,&rend);CHKERRQ(ierr);
159c20d7725SJed Brown   for (i=0; i<nzp; i++) {
160c20d7725SJed Brown     ierr = PetscRandomGetValue(rdm,&a[i]);CHKERRQ(ierr);
161c20d7725SJed Brown   }
162c20d7725SJed Brown   for (i=rstart; i<rend; i++) {
163c20d7725SJed Brown     for (j=0; j<nzp; j++) {
164c20d7725SJed Brown       ierr    = PetscRandomGetValue(rdm,&rval);CHKERRQ(ierr);
165c20d7725SJed Brown       idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
166c20d7725SJed Brown     }
167c20d7725SJed Brown     ierr = MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);CHKERRQ(ierr);
168c20d7725SJed Brown   }
169c20d7725SJed Brown   ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
170c20d7725SJed Brown   ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
17120b3374bSStefano Zampini   ierr = MatSetFromOptions(P);CHKERRQ(ierr);
172c20d7725SJed Brown 
173c20d7725SJed Brown   ierr = MatTranspose(P,MAT_INITIAL_MATRIX,&R);CHKERRQ(ierr);
17420b3374bSStefano Zampini   ierr = MatSetFromOptions(R);CHKERRQ(ierr);
175c20d7725SJed Brown 
176c20d7725SJed Brown   /* 2) MatTransposeMatMult() */
177c20d7725SJed Brown   /* ------------------------ */
178c20d7725SJed Brown   if (Test_MatTrMat) {
179c20d7725SJed Brown     /* (2.1) Test developer driver C = P^T*B */
180c20d7725SJed Brown     ierr = MatProductCreate(P,B,NULL,&C);CHKERRQ(ierr);
18167b3012eSStefano Zampini     ierr = MatSetOptionsPrefix(C,"AtB_");CHKERRQ(ierr);
182c20d7725SJed Brown     ierr = MatProductSetType(C,MATPRODUCT_AtB);CHKERRQ(ierr);
183910fa13fSStefano Zampini     ierr = MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr);
184c20d7725SJed Brown     ierr = MatProductSetFill(C,PETSC_DEFAULT);CHKERRQ(ierr);
185c20d7725SJed Brown     ierr = MatProductSetFromOptions(C);CHKERRQ(ierr);
1860e1ed72eSHong Zhang     ierr = MatProductSymbolic(C);CHKERRQ(ierr); /* equivalent to MatSetUp() */
1870ad02fcaSStefano Zampini     ierr = MatSetOption(C,MAT_USE_INODES,PETSC_FALSE);CHKERRQ(ierr); /* illustrate how to call MatSetOption() */
188c20d7725SJed Brown     ierr = MatProductNumeric(C);CHKERRQ(ierr);
18967b3012eSStefano Zampini     ierr = MatProductNumeric(C);CHKERRQ(ierr); /* test reuse symbolic C */
19067b3012eSStefano Zampini 
19120b3374bSStefano Zampini     ierr = MatTransposeMatMultEqual(P,B,C,mcheck,&flg);CHKERRQ(ierr);
192c20d7725SJed Brown     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error: developer driver C = P^T*B");
193c20d7725SJed Brown     ierr = MatDestroy(&C);CHKERRQ(ierr);
194c20d7725SJed Brown 
195c20d7725SJed Brown     /* (2.2) Test user driver C = P^T*B */
196c20d7725SJed Brown     ierr = MatTransposeMatMult(P,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
197c20d7725SJed Brown     ierr = MatTransposeMatMult(P,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
198c20d7725SJed Brown     ierr = MatGetInfo(C,MAT_GLOBAL_SUM,&info);CHKERRQ(ierr);
1996718818eSStefano Zampini     ierr = MatProductClear(C);CHKERRQ(ierr);
200c20d7725SJed Brown 
201c20d7725SJed Brown     /* Compare P^T*B and R*B */
202c20d7725SJed Brown     ierr = MatMatMult(R,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C1);CHKERRQ(ierr);
203c20d7725SJed Brown     ierr = MatNormDifference(C,C1,&norm);CHKERRQ(ierr);
204c20d7725SJed Brown     if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatTransposeMatMult(): %g",(double)norm);
205c20d7725SJed Brown     ierr = MatDestroy(&C1);CHKERRQ(ierr);
206c20d7725SJed Brown 
207c20d7725SJed Brown     /* Test MatDuplicate() of C=P^T*B */
208c20d7725SJed Brown     ierr = MatDuplicate(C,MAT_COPY_VALUES,&C1);CHKERRQ(ierr);
209c20d7725SJed Brown     ierr = MatDestroy(&C1);CHKERRQ(ierr);
210c20d7725SJed Brown     ierr = MatDestroy(&C);CHKERRQ(ierr);
211c20d7725SJed Brown   }
212c20d7725SJed Brown 
21367b3012eSStefano Zampini   /* 3) MatMatTransposeMult() */
214c20d7725SJed Brown   /* ------------------------ */
215c20d7725SJed Brown   if (Test_MatMatTr) {
216c20d7725SJed Brown     /* C = B*R^T */
21720b3374bSStefano Zampini     ierr = PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);CHKERRQ(ierr);
21820b3374bSStefano Zampini     if (seqaij) {
219c20d7725SJed Brown       ierr = MatMatTransposeMult(B,R,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
22067b3012eSStefano Zampini       ierr = MatSetOptionsPrefix(C,"ABt_");CHKERRQ(ierr); /* enable '-ABt_' for matrix C */
221c20d7725SJed Brown       ierr = MatGetInfo(C,MAT_GLOBAL_SUM,&info);CHKERRQ(ierr);
222c20d7725SJed Brown 
223c20d7725SJed Brown       /* Test MAT_REUSE_MATRIX - reuse symbolic C */
224c20d7725SJed Brown       ierr = MatMatTransposeMult(B,R,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
225c20d7725SJed Brown 
226c20d7725SJed Brown       /* Check */
227c20d7725SJed Brown       ierr = MatMatMult(B,P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C1);CHKERRQ(ierr);
228c20d7725SJed Brown       ierr = MatNormDifference(C,C1,&norm);CHKERRQ(ierr);
229c20d7725SJed Brown       if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatMatTransposeMult() %g",(double)norm);
230c20d7725SJed Brown       ierr = MatDestroy(&C1);CHKERRQ(ierr);
231c20d7725SJed Brown       ierr = MatDestroy(&C);CHKERRQ(ierr);
232c20d7725SJed Brown     }
233c20d7725SJed Brown   }
234c20d7725SJed Brown 
235c20d7725SJed Brown   /* 4) Test MatPtAP() */
236c20d7725SJed Brown   /*-------------------*/
237c20d7725SJed Brown   if (Test_MatPtAP) {
238c20d7725SJed Brown     ierr = MatDuplicate(A_save,MAT_COPY_VALUES,&A);CHKERRQ(ierr);
239c20d7725SJed Brown 
240c20d7725SJed Brown     /* (4.1) Test developer API */
241c20d7725SJed Brown     ierr = MatProductCreate(A,P,NULL,&C);CHKERRQ(ierr);
24267b3012eSStefano Zampini     ierr = MatSetOptionsPrefix(C,"PtAP_");CHKERRQ(ierr);
243c20d7725SJed Brown     ierr = MatProductSetType(C,MATPRODUCT_PtAP);CHKERRQ(ierr);
244910fa13fSStefano Zampini     ierr = MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr);
245c20d7725SJed Brown     ierr = MatProductSetFill(C,PETSC_DEFAULT);CHKERRQ(ierr);
246c20d7725SJed Brown     ierr = MatProductSetFromOptions(C);CHKERRQ(ierr);
247c20d7725SJed Brown     ierr = MatProductSymbolic(C);CHKERRQ(ierr);
248c20d7725SJed Brown     ierr = MatProductNumeric(C);CHKERRQ(ierr);
24920b3374bSStefano Zampini     ierr = MatPtAPMultEqual(A,P,C,mcheck,&flg);CHKERRQ(ierr);
25020b3374bSStefano Zampini     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatProduct_PtAP");
251c20d7725SJed Brown     ierr = MatProductNumeric(C);CHKERRQ(ierr); /* reuse symbolic C */
252c20d7725SJed Brown 
25320b3374bSStefano Zampini     ierr = MatPtAPMultEqual(A,P,C,mcheck,&flg);CHKERRQ(ierr);
254c20d7725SJed Brown     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatProduct_PtAP");
255c20d7725SJed Brown     ierr = MatDestroy(&C);CHKERRQ(ierr);
256c20d7725SJed Brown 
257c20d7725SJed Brown     /* (4.2) Test user driver */
258c20d7725SJed Brown     ierr = MatPtAP(A,P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
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;
264c20d7725SJed Brown       ierr   = MatScale(A,alpha);CHKERRQ(ierr);
265c20d7725SJed Brown       ierr   = MatPtAP(A,P,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
266c20d7725SJed Brown     }
26720b3374bSStefano Zampini     ierr = MatPtAPMultEqual(A,P,C,mcheck,&flg);CHKERRQ(ierr);
268c20d7725SJed Brown     if (!flg) SETERRQ(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;
274c20d7725SJed Brown       ierr = MatTranspose(P,MAT_REUSE_MATRIX,&R);CHKERRQ(ierr);
27567b3012eSStefano Zampini 
27667b3012eSStefano Zampini       /* (5.1) Test developer driver RARt = R*A*Rt */
27767b3012eSStefano Zampini       ierr = MatProductCreate(A,R,NULL,&RARt);CHKERRQ(ierr);
27867b3012eSStefano Zampini       ierr = MatSetOptionsPrefix(RARt,"RARt_");CHKERRQ(ierr);
27967b3012eSStefano Zampini       ierr = MatProductSetType(RARt,MATPRODUCT_RARt);CHKERRQ(ierr);
280910fa13fSStefano Zampini       ierr = MatProductSetAlgorithm(RARt,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr);
28167b3012eSStefano Zampini       ierr = MatProductSetFill(RARt,PETSC_DEFAULT);CHKERRQ(ierr);
28267b3012eSStefano Zampini       ierr = MatProductSetFromOptions(RARt);CHKERRQ(ierr);
28367b3012eSStefano Zampini       ierr = MatProductSymbolic(RARt);CHKERRQ(ierr); /* equivalent to MatSetUp() */
28467b3012eSStefano Zampini       ierr = MatSetOption(RARt,MAT_USE_INODES,PETSC_FALSE);CHKERRQ(ierr); /* illustrate how to call MatSetOption() */
28567b3012eSStefano Zampini       ierr = MatProductNumeric(RARt);CHKERRQ(ierr);
28667b3012eSStefano Zampini       ierr = MatProductNumeric(RARt);CHKERRQ(ierr); /* test reuse symbolic RARt */
28767b3012eSStefano Zampini       ierr = MatDestroy(&RARt);CHKERRQ(ierr);
28867b3012eSStefano Zampini 
28967b3012eSStefano Zampini       /* (2.2) Test user driver RARt = R*A*Rt */
290c20d7725SJed Brown       ierr = MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&RARt);CHKERRQ(ierr);
291c20d7725SJed Brown       ierr = MatRARt(A,R,MAT_REUSE_MATRIX,2.0,&RARt);CHKERRQ(ierr);
29267b3012eSStefano Zampini 
293c20d7725SJed Brown       ierr = MatNormDifference(C,RARt,&norm);CHKERRQ(ierr);
294c20d7725SJed Brown       if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"|PtAP - RARt| = %g",(double)norm);
295c20d7725SJed Brown       ierr = MatDestroy(&RARt);CHKERRQ(ierr);
296c20d7725SJed Brown     }
297c20d7725SJed Brown 
298c20d7725SJed Brown     ierr = MatDestroy(&A);CHKERRQ(ierr);
299c20d7725SJed Brown     ierr = MatDestroy(&C);CHKERRQ(ierr);
300c20d7725SJed Brown   }
301c20d7725SJed Brown 
302c20d7725SJed Brown   /* Destroy objects */
303c20d7725SJed Brown   ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr);
304c20d7725SJed Brown   ierr = PetscFree(idxn);CHKERRQ(ierr);
305c20d7725SJed Brown 
306c20d7725SJed Brown   ierr = MatDestroy(&A_save);CHKERRQ(ierr);
307c20d7725SJed Brown   ierr = MatDestroy(&B);CHKERRQ(ierr);
308c20d7725SJed Brown   ierr = MatDestroy(&P);CHKERRQ(ierr);
309c20d7725SJed Brown   ierr = MatDestroy(&R);CHKERRQ(ierr);
310c20d7725SJed Brown 
311c20d7725SJed Brown   PetscFinalize();
312c20d7725SJed Brown   return ierr;
313c20d7725SJed Brown }
314c20d7725SJed Brown 
315c20d7725SJed Brown /*TEST
316c20d7725SJed Brown    test:
317c20d7725SJed Brown      suffix: 1
318*dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
319c20d7725SJed Brown      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium
320c20d7725SJed Brown      output_file: output/ex62_1.out
321c20d7725SJed Brown 
322c20d7725SJed Brown    test:
323c20d7725SJed Brown      suffix: 2_ab_scalable
324*dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
32567b3012eSStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via scalable -matmatmult_via scalable -AtB_matproduct_atb_via outerproduct -mattransposematmult_via outerproduct
326c20d7725SJed Brown      output_file: output/ex62_1.out
327c20d7725SJed Brown 
328c20d7725SJed Brown    test:
329c20d7725SJed Brown      suffix: 3_ab_scalable_fast
330*dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
33167b3012eSStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via scalable_fast -matmatmult_via scalable_fast -matmattransmult_via color
332c20d7725SJed Brown      output_file: output/ex62_1.out
333c20d7725SJed Brown 
334c20d7725SJed Brown    test:
335c20d7725SJed Brown      suffix: 4_ab_heap
336*dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
33767b3012eSStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via heap -matmatmult_via heap -PtAP_matproduct_ptap_via rap -matptap_via rap
338c20d7725SJed Brown      output_file: output/ex62_1.out
339c20d7725SJed Brown 
340c20d7725SJed Brown    test:
341c20d7725SJed Brown      suffix: 5_ab_btheap
342*dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
34367b3012eSStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via btheap -matmatmult_via btheap -matrart_via r*art
344c20d7725SJed Brown      output_file: output/ex62_1.out
345c20d7725SJed Brown 
346c20d7725SJed Brown    test:
347c20d7725SJed Brown      suffix: 6_ab_llcondensed
348*dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
34967b3012eSStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via llcondensed -matmatmult_via llcondensed -matrart_via coloring_rart
350c20d7725SJed Brown      output_file: output/ex62_1.out
351c20d7725SJed Brown 
352c20d7725SJed Brown    test:
353c20d7725SJed Brown      suffix: 7_ab_rowmerge
354*dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
35567b3012eSStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via rowmerge -matmatmult_via rowmerge
356c20d7725SJed Brown      output_file: output/ex62_1.out
357c20d7725SJed Brown 
358c20d7725SJed Brown    test:
359c20d7725SJed Brown      suffix: 8_ab_hypre
360*dfd57a17SPierre Jolivet      requires: hypre datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
36167b3012eSStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via hypre -matmatmult_via hypre -PtAP_matproduct_ptap_via hypre -matptap_via hypre
362c20d7725SJed Brown      output_file: output/ex62_1.out
363c20d7725SJed Brown 
364c20d7725SJed Brown    test:
365cec0a6c6SStefano Zampini      suffix: 9_mkl
36615df37b9SStefano Zampini      TODO: broken MatScale?
367*dfd57a17SPierre Jolivet      requires: mkl datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
368cec0a6c6SStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -A_mat_type aijmkl -B_mat_type aijmkl
369cec0a6c6SStefano Zampini      output_file: output/ex62_1.out
370cec0a6c6SStefano Zampini 
371cec0a6c6SStefano Zampini    test:
372c20d7725SJed Brown      suffix: 10
373*dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
374c20d7725SJed Brown      nsize: 3
375c20d7725SJed Brown      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium
376c20d7725SJed Brown      output_file: output/ex62_1.out
377c20d7725SJed Brown 
378c20d7725SJed Brown    test:
3793cea4f0aSStefano Zampini      suffix: 10_backend
380*dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
3813cea4f0aSStefano Zampini      nsize: 3
3823cea4f0aSStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via backend -matmatmult_via backend -AtB_matproduct_atb_via backend -mattransposematmult_via backend -PtAP_matproduct_ptap_via backend -matptap_via backend
3833cea4f0aSStefano Zampini      output_file: output/ex62_1.out
3843cea4f0aSStefano Zampini 
3853cea4f0aSStefano Zampini    test:
386c20d7725SJed Brown      suffix: 11_ab_scalable
387*dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
388c20d7725SJed Brown      nsize: 3
38967b3012eSStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via scalable -matmatmult_via scalable -AtB_matproduct_atb_via scalable -mattransposematmult_via scalable
390c20d7725SJed Brown      output_file: output/ex62_1.out
391c20d7725SJed Brown 
392c20d7725SJed Brown    test:
393c20d7725SJed Brown      suffix: 12_ab_seqmpi
394*dfd57a17SPierre Jolivet      requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
395c20d7725SJed Brown      nsize: 3
39667b3012eSStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via seqmpi -matmatmult_via seqmpi -AtB_matproduct_atb_via at*b -mattransposematmult_via at*b
397c20d7725SJed Brown      output_file: output/ex62_1.out
398c20d7725SJed Brown 
399c20d7725SJed Brown    test:
400c20d7725SJed Brown      suffix: 13_ab_hypre
401*dfd57a17SPierre Jolivet      requires: hypre datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
402c20d7725SJed Brown      nsize: 3
40367b3012eSStefano Zampini      args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via hypre -matmatmult_via hypre -PtAP_matproduct_ptap_via hypre -matptap_via hypre
404c20d7725SJed Brown      output_file: output/ex62_1.out
405c20d7725SJed Brown 
40620b3374bSStefano Zampini    test:
40720b3374bSStefano Zampini      suffix: 14_seqaij
408*dfd57a17SPierre Jolivet      requires: !complex double !defined(PETSC_USE_64BIT_INDICES)
40920b3374bSStefano Zampini      args: -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
41020b3374bSStefano Zampini      output_file: output/ex62_1.out
41120b3374bSStefano Zampini 
41220b3374bSStefano Zampini    test:
41320b3374bSStefano Zampini      suffix: 14_seqaijcusparse
414*dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4151a2c6b5cSJunchao 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
41620b3374bSStefano Zampini      output_file: output/ex62_1.out
41720b3374bSStefano Zampini 
41820b3374bSStefano Zampini    test:
41965e4b4d4SStefano Zampini      suffix: 14_seqaijcusparse_cpu
420*dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
42165e4b4d4SStefano Zampini      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_matproduct_ab_backend_cpu -matmatmult_backend_cpu -PtAP_matproduct_ptap_backend_cpu -matptap_backend_cpu -RARt_matproduct_rart_backend_cpu -matrart_backend_cpu
42265e4b4d4SStefano Zampini      output_file: output/ex62_1.out
42365e4b4d4SStefano Zampini 
42465e4b4d4SStefano Zampini    test:
4255cb69cabSStefano Zampini      suffix: 14_mpiaijcusparse_seq
4265cb69cabSStefano Zampini      nsize: 1
427*dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4281a2c6b5cSJunchao 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
4295cb69cabSStefano Zampini      output_file: output/ex62_1.out
4305cb69cabSStefano Zampini 
4315cb69cabSStefano Zampini    test:
43265e4b4d4SStefano Zampini      suffix: 14_mpiaijcusparse_seq_cpu
43365e4b4d4SStefano Zampini      nsize: 1
434*dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
43565e4b4d4SStefano Zampini      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_matproduct_ab_backend_cpu -matmatmult_backend_cpu -PtAP_matproduct_ptap_backend_cpu -matptap_backend_cpu
43665e4b4d4SStefano Zampini      output_file: output/ex62_1.out
43765e4b4d4SStefano Zampini 
43865e4b4d4SStefano Zampini    test:
4395cb69cabSStefano Zampini      suffix: 14_mpiaijcusparse
4405cb69cabSStefano Zampini      nsize: 3
441*dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
4421a2c6b5cSJunchao 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
4435cb69cabSStefano Zampini      output_file: output/ex62_1.out
4445cb69cabSStefano Zampini 
4455cb69cabSStefano Zampini    test:
44665e4b4d4SStefano Zampini      suffix: 14_mpiaijcusparse_cpu
44765e4b4d4SStefano Zampini      nsize: 3
448*dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES)
44965e4b4d4SStefano Zampini      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_matproduct_ab_backend_cpu -matmatmult_backend_cpu -PtAP_matproduct_ptap_backend_cpu -matptap_backend_cpu
45065e4b4d4SStefano Zampini      output_file: output/ex62_1.out
45165e4b4d4SStefano Zampini 
45265e4b4d4SStefano Zampini    test:
45320b3374bSStefano Zampini      suffix: 14_seqaijkokkos
454*dfd57a17SPierre Jolivet      requires: kokkos_kernels !complex double !defined(PETSC_USE_64BIT_INDICES)
45520b3374bSStefano 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
45620b3374bSStefano Zampini      output_file: output/ex62_1.out
45720b3374bSStefano Zampini 
4585cb69cabSStefano Zampini    # these tests use matrices with many zero rows
45920b3374bSStefano Zampini    test:
46020b3374bSStefano Zampini      suffix: 15_seqaijcusparse
461*dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES) datafilespath
4621a2c6b5cSJunchao Zhang      args: -A_mat_type aijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
46320b3374bSStefano Zampini      output_file: output/ex62_1.out
4645cb69cabSStefano Zampini 
4655cb69cabSStefano Zampini    test:
4665cb69cabSStefano Zampini      suffix: 15_mpiaijcusparse_seq
4675cb69cabSStefano Zampini      nsize: 1
468*dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES) datafilespath
4691a2c6b5cSJunchao Zhang      args: -A_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
4705cb69cabSStefano Zampini      output_file: output/ex62_1.out
4715cb69cabSStefano Zampini 
4725cb69cabSStefano Zampini    test:
4735cb69cabSStefano Zampini      nsize: 3
4745cb69cabSStefano Zampini      suffix: 15_mpiaijcusparse
475*dfd57a17SPierre Jolivet      requires: cuda !complex double !defined(PETSC_USE_64BIT_INDICES) datafilespath
4761a2c6b5cSJunchao Zhang      args: -A_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
4775cb69cabSStefano Zampini      output_file: output/ex62_1.out
4785cb69cabSStefano Zampini 
479c20d7725SJed Brown TEST*/
480