xref: /petsc/src/mat/tests/ex66.c (revision 327415f76d85372a4417cf1aaa14db707d4d6c04)
153022affSStefano Zampini static char help[] = "Tests MATH2OPUS\n\n";
2c4fb69feSStefano Zampini 
3c4fb69feSStefano Zampini #include <petscmat.h>
4c4fb69feSStefano Zampini #include <petscsf.h>
5c4fb69feSStefano Zampini 
6c4fb69feSStefano Zampini static PetscScalar GenEntry_Symm(PetscInt sdim, PetscReal x[], PetscReal y[], void *ctx)
7c4fb69feSStefano Zampini {
8c4fb69feSStefano Zampini     PetscInt  d;
9c4fb69feSStefano Zampini     PetscReal clength = sdim == 3 ? 0.2 : 0.1;
10c4fb69feSStefano Zampini     PetscReal dist, diff = 0.0;
11c4fb69feSStefano Zampini 
12c4fb69feSStefano Zampini     for (d = 0; d < sdim; d++) { diff += (x[d] - y[d]) * (x[d] - y[d]); }
13c4fb69feSStefano Zampini     dist = PetscSqrtReal(diff);
14c4fb69feSStefano Zampini     return PetscExpReal(-dist / clength);
15c4fb69feSStefano Zampini }
16c4fb69feSStefano Zampini 
17c4fb69feSStefano Zampini static PetscScalar GenEntry_Unsymm(PetscInt sdim, PetscReal x[], PetscReal y[], void *ctx)
18c4fb69feSStefano Zampini {
19c4fb69feSStefano Zampini     PetscInt  d;
20c4fb69feSStefano Zampini     PetscReal clength = sdim == 3 ? 0.2 : 0.1;
21c4fb69feSStefano Zampini     PetscReal dist, diff = 0.0, nx = 0.0, ny = 0.0;
22c4fb69feSStefano Zampini 
23c4fb69feSStefano Zampini     for (d = 0; d < sdim; d++) { nx += x[d]*x[d]; }
24c4fb69feSStefano Zampini     for (d = 0; d < sdim; d++) { ny += y[d]*y[d]; }
25c4fb69feSStefano Zampini     for (d = 0; d < sdim; d++) { diff += (x[d] - y[d]) * (x[d] - y[d]); }
26c4fb69feSStefano Zampini     dist = PetscSqrtReal(diff);
27c4fb69feSStefano Zampini     return nx > ny ? PetscExpReal(-dist / clength) : PetscExpReal(-dist / clength) + 1.;
28c4fb69feSStefano Zampini }
29c4fb69feSStefano Zampini 
30c4fb69feSStefano Zampini int main(int argc,char **argv)
31c4fb69feSStefano Zampini {
32c4fb69feSStefano Zampini   Mat            A,B,C,D;
33c4fb69feSStefano Zampini   Vec            v,x,y,Ax,Ay,Bx,By;
34c4fb69feSStefano Zampini   PetscRandom    r;
35c4fb69feSStefano Zampini   PetscLayout    map;
36c4fb69feSStefano Zampini   PetscScalar    *Adata = NULL, *Cdata = NULL, scale = 1.0;
37c4fb69feSStefano Zampini   PetscReal      *coords,nA,nD,nB,err,nX,norms[3];
38300d917bSStefano Zampini   PetscInt       N, n = 64, dim = 1, i, j, nrhs = 11, lda = 0, ldc = 0, ldu = 0, nlr = 7, nt, ntrials = 2;
39c4fb69feSStefano Zampini   PetscMPIInt    size,rank;
40c4fb69feSStefano Zampini   PetscBool      testlayout = PETSC_FALSE, flg, symm = PETSC_FALSE, Asymm = PETSC_TRUE, kernel = PETSC_TRUE;
4153022affSStefano Zampini   PetscBool      checkexpl = PETSC_FALSE, agpu = PETSC_FALSE, bgpu = PETSC_FALSE, cgpu = PETSC_FALSE, flgglob;
42300d917bSStefano Zampini   PetscBool      testtrans, testnorm, randommat = PETSC_TRUE, testorthog, testcompress, testhlru;
43c4fb69feSStefano Zampini   void           (*approxnormfunc)(void);
44c4fb69feSStefano Zampini   void           (*Anormfunc)(void);
45c4fb69feSStefano Zampini 
4653022affSStefano Zampini #if defined(PETSC_HAVE_MPI_INIT_THREAD)
47c4fb69feSStefano Zampini   PETSC_MPI_THREAD_REQUIRED = MPI_THREAD_MULTIPLE;
4853022affSStefano Zampini #endif
49*327415f7SBarry Smith   PetscFunctionBeginUser;
509566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc,&argv,(char*) 0,help));
519566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-ng",&N,&flgglob));
529566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
539566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL));
549566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL));
559566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-lda",&lda,NULL));
569566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-ldc",&ldc,NULL));
579566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-nlr",&nlr,NULL));
589566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-ldu",&ldu,NULL));
599566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-matmattrials",&ntrials,NULL));
609566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-randommat",&randommat,NULL));
619566063dSJacob Faibussowitsch   if (!flgglob) PetscCall(PetscOptionsGetBool(NULL,NULL,"-testlayout",&testlayout,NULL));
629566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-Asymm",&Asymm,NULL));
639566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-symm",&symm,NULL));
649566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-kernel",&kernel,NULL));
659566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-checkexpl",&checkexpl,NULL));
669566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-agpu",&agpu,NULL));
679566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-bgpu",&bgpu,NULL));
689566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-cgpu",&cgpu,NULL));
699566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-scale",&scale,NULL));
70c4fb69feSStefano Zampini   if (!Asymm) symm = PETSC_FALSE;
71c4fb69feSStefano Zampini 
729566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
73300d917bSStefano Zampini 
74300d917bSStefano Zampini   /* Disable tests for unimplemented variants */
75c4fb69feSStefano Zampini   testtrans = (PetscBool)(size == 1 || symm);
76c4fb69feSStefano Zampini   testnorm = (PetscBool)(size == 1 || symm);
7753022affSStefano Zampini   testorthog = (PetscBool)(size == 1 || symm);
7853022affSStefano Zampini   testcompress = (PetscBool)(size == 1 || symm);
79300d917bSStefano Zampini   testhlru = (PetscBool)(size == 1);
80c4fb69feSStefano Zampini 
819566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
829566063dSJacob Faibussowitsch   PetscCall(PetscLayoutCreate(PETSC_COMM_WORLD,&map));
83c4fb69feSStefano Zampini   if (testlayout) {
84c4fb69feSStefano Zampini     if (rank%2) n = PetscMax(2*n-5*rank,0);
85c4fb69feSStefano Zampini     else n = 2*n+rank;
86c4fb69feSStefano Zampini   }
8753022affSStefano Zampini   if (!flgglob) {
889566063dSJacob Faibussowitsch     PetscCall(PetscLayoutSetLocalSize(map,n));
899566063dSJacob Faibussowitsch     PetscCall(PetscLayoutSetUp(map));
909566063dSJacob Faibussowitsch     PetscCall(PetscLayoutGetSize(map,&N));
9153022affSStefano Zampini   } else {
929566063dSJacob Faibussowitsch     PetscCall(PetscLayoutSetSize(map,N));
939566063dSJacob Faibussowitsch     PetscCall(PetscLayoutSetUp(map));
949566063dSJacob Faibussowitsch     PetscCall(PetscLayoutGetLocalSize(map,&n));
9553022affSStefano Zampini   }
969566063dSJacob Faibussowitsch   PetscCall(PetscLayoutDestroy(&map));
97c4fb69feSStefano Zampini 
98c4fb69feSStefano Zampini   if (lda) {
999566063dSJacob Faibussowitsch     PetscCall(PetscMalloc1(N*(n+lda),&Adata));
100c4fb69feSStefano Zampini   }
1019566063dSJacob Faibussowitsch   PetscCall(MatCreateDense(PETSC_COMM_WORLD,n,n,N,N,Adata,&A));
1029566063dSJacob Faibussowitsch   PetscCall(MatDenseSetLDA(A,n+lda));
103c4fb69feSStefano Zampini 
10453022affSStefano Zampini   /* Create random points; these are replicated in order to populate a dense matrix and to compare sequential and dense runs
10553022affSStefano Zampini      The constructor for MATH2OPUS can take as input the distributed coordinates and replicates them internally in case
10653022affSStefano Zampini      a kernel construction is requested */
1079566063dSJacob Faibussowitsch   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&r));
1089566063dSJacob Faibussowitsch   PetscCall(PetscRandomSetFromOptions(r));
1099566063dSJacob Faibussowitsch   PetscCall(PetscRandomSetSeed(r,123456));
1109566063dSJacob Faibussowitsch   PetscCall(PetscRandomSeed(r));
1119566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(N*dim,&coords));
1129566063dSJacob Faibussowitsch   PetscCall(PetscRandomGetValuesReal(r,N*dim,coords));
1139566063dSJacob Faibussowitsch   PetscCall(PetscRandomDestroy(&r));
11453022affSStefano Zampini 
115c4fb69feSStefano Zampini   if (kernel || !randommat) {
11653022affSStefano Zampini     MatH2OpusKernel k = Asymm ? GenEntry_Symm : GenEntry_Unsymm;
117c4fb69feSStefano Zampini     PetscInt        ist,ien;
118c4fb69feSStefano Zampini 
1199566063dSJacob Faibussowitsch     PetscCall(MatGetOwnershipRange(A,&ist,&ien));
120c4fb69feSStefano Zampini     for (i = ist; i < ien; i++) {
121c4fb69feSStefano Zampini       for (j = 0; j < N; j++) {
1229566063dSJacob Faibussowitsch         PetscCall(MatSetValue(A,i,j,(*k)(dim,coords + i*dim,coords + j*dim,NULL),INSERT_VALUES));
123c4fb69feSStefano Zampini       }
124c4fb69feSStefano Zampini     }
1259566063dSJacob Faibussowitsch     PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
1269566063dSJacob Faibussowitsch     PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
127c4fb69feSStefano Zampini     if (kernel) {
1289566063dSJacob Faibussowitsch       PetscCall(MatCreateH2OpusFromKernel(PETSC_COMM_WORLD,n,n,N,N,dim,coords + ist*dim,PETSC_TRUE,k,NULL,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,&B));
129c4fb69feSStefano Zampini     } else {
1309566063dSJacob Faibussowitsch       PetscCall(MatCreateH2OpusFromMat(A,dim,coords + ist*dim,PETSC_TRUE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,&B));
131c4fb69feSStefano Zampini     }
132c4fb69feSStefano Zampini   } else {
13353022affSStefano Zampini     PetscInt ist;
13453022affSStefano Zampini 
1359566063dSJacob Faibussowitsch     PetscCall(MatGetOwnershipRange(A,&ist,NULL));
1369566063dSJacob Faibussowitsch     PetscCall(MatSetRandom(A,NULL));
137c4fb69feSStefano Zampini     if (Asymm) {
1389566063dSJacob Faibussowitsch       PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX,&B));
1399566063dSJacob Faibussowitsch       PetscCall(MatAXPY(A,1.0,B,SAME_NONZERO_PATTERN));
1409566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&B));
1419566063dSJacob Faibussowitsch       PetscCall(MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE));
142c4fb69feSStefano Zampini     }
1439566063dSJacob Faibussowitsch     PetscCall(MatCreateH2OpusFromMat(A,dim,coords + ist*dim,PETSC_TRUE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,&B));
144c4fb69feSStefano Zampini   }
1459566063dSJacob Faibussowitsch   PetscCall(PetscFree(coords));
146c4fb69feSStefano Zampini   if (agpu) {
1479566063dSJacob Faibussowitsch     PetscCall(MatConvert(A,MATDENSECUDA,MAT_INPLACE_MATRIX,&A));
148c4fb69feSStefano Zampini   }
1499566063dSJacob Faibussowitsch   PetscCall(MatViewFromOptions(A,NULL,"-A_view"));
150c4fb69feSStefano Zampini 
1519566063dSJacob Faibussowitsch   PetscCall(MatSetOption(B,MAT_SYMMETRIC,symm));
152c4fb69feSStefano Zampini 
153c4fb69feSStefano Zampini   /* assemble the H-matrix */
1549566063dSJacob Faibussowitsch   PetscCall(MatBindToCPU(B,(PetscBool)!bgpu));
1559566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(B));
1569566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
1579566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
1589566063dSJacob Faibussowitsch   PetscCall(MatViewFromOptions(B,NULL,"-B_view"));
159c4fb69feSStefano Zampini 
160c4fb69feSStefano Zampini   /* Test MatScale */
1619566063dSJacob Faibussowitsch   PetscCall(MatScale(A,scale));
1629566063dSJacob Faibussowitsch   PetscCall(MatScale(B,scale));
163c4fb69feSStefano Zampini 
164c4fb69feSStefano Zampini   /* Test MatMult */
1659566063dSJacob Faibussowitsch   PetscCall(MatCreateVecs(A,&Ax,&Ay));
1669566063dSJacob Faibussowitsch   PetscCall(MatCreateVecs(B,&Bx,&By));
1679566063dSJacob Faibussowitsch   PetscCall(VecSetRandom(Ax,NULL));
1689566063dSJacob Faibussowitsch   PetscCall(VecCopy(Ax,Bx));
1699566063dSJacob Faibussowitsch   PetscCall(MatMult(A,Ax,Ay));
1709566063dSJacob Faibussowitsch   PetscCall(MatMult(B,Bx,By));
1719566063dSJacob Faibussowitsch   PetscCall(VecViewFromOptions(Ay,NULL,"-mult_vec_view"));
1729566063dSJacob Faibussowitsch   PetscCall(VecViewFromOptions(By,NULL,"-mult_vec_view"));
1739566063dSJacob Faibussowitsch   PetscCall(VecNorm(Ay,NORM_INFINITY,&nX));
1749566063dSJacob Faibussowitsch   PetscCall(VecAXPY(Ay,-1.0,By));
1759566063dSJacob Faibussowitsch   PetscCall(VecViewFromOptions(Ay,NULL,"-mult_vec_view"));
1769566063dSJacob Faibussowitsch   PetscCall(VecNorm(Ay,NORM_INFINITY,&err));
1779566063dSJacob Faibussowitsch   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMult err %g\n",err/nX));
1789566063dSJacob Faibussowitsch   PetscCall(VecScale(By,-1.0));
1799566063dSJacob Faibussowitsch   PetscCall(MatMultAdd(B,Bx,By,By));
1809566063dSJacob Faibussowitsch   PetscCall(VecNorm(By,NORM_INFINITY,&err));
1819566063dSJacob Faibussowitsch   PetscCall(VecViewFromOptions(By,NULL,"-mult_vec_view"));
1824c3a5f3cSStefano Zampini   if (err > 10.*PETSC_SMALL) {
1839566063dSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMultAdd err %g\n",err));
184c4fb69feSStefano Zampini   }
185c4fb69feSStefano Zampini 
186c4fb69feSStefano Zampini   /* Test MatNorm */
1879566063dSJacob Faibussowitsch   PetscCall(MatNorm(A,NORM_INFINITY,&norms[0]));
1889566063dSJacob Faibussowitsch   PetscCall(MatNorm(A,NORM_1,&norms[1]));
189c4fb69feSStefano Zampini   norms[2] = -1.; /* NORM_2 not supported */
1909566063dSJacob Faibussowitsch   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"A Matrix norms:        infty=%g, norm_1=%g, norm_2=%g\n",(double)norms[0],(double)norms[1],(double)norms[2]));
1919566063dSJacob Faibussowitsch   PetscCall(MatGetOperation(A,MATOP_NORM,&Anormfunc));
1929566063dSJacob Faibussowitsch   PetscCall(MatGetOperation(B,MATOP_NORM,&approxnormfunc));
1939566063dSJacob Faibussowitsch   PetscCall(MatSetOperation(A,MATOP_NORM,approxnormfunc));
1949566063dSJacob Faibussowitsch   PetscCall(MatNorm(A,NORM_INFINITY,&norms[0]));
1959566063dSJacob Faibussowitsch   PetscCall(MatNorm(A,NORM_1,&norms[1]));
1969566063dSJacob Faibussowitsch   PetscCall(MatNorm(A,NORM_2,&norms[2]));
1979566063dSJacob Faibussowitsch   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"A Approx Matrix norms: infty=%g, norm_1=%g, norm_2=%g\n",(double)norms[0],(double)norms[1],(double)norms[2]));
198c4fb69feSStefano Zampini   if (testnorm) {
1999566063dSJacob Faibussowitsch     PetscCall(MatNorm(B,NORM_INFINITY,&norms[0]));
2009566063dSJacob Faibussowitsch     PetscCall(MatNorm(B,NORM_1,&norms[1]));
2019566063dSJacob Faibussowitsch     PetscCall(MatNorm(B,NORM_2,&norms[2]));
202c4fb69feSStefano Zampini   } else {
203c4fb69feSStefano Zampini     norms[0] = -1.;
204c4fb69feSStefano Zampini     norms[1] = -1.;
205c4fb69feSStefano Zampini     norms[2] = -1.;
206c4fb69feSStefano Zampini   }
2079566063dSJacob Faibussowitsch   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"B Approx Matrix norms: infty=%g, norm_1=%g, norm_2=%g\n",(double)norms[0],(double)norms[1],(double)norms[2]));
2089566063dSJacob Faibussowitsch   PetscCall(MatSetOperation(A,MATOP_NORM,Anormfunc));
209c4fb69feSStefano Zampini 
210c4fb69feSStefano Zampini   /* Test MatDuplicate */
2119566063dSJacob Faibussowitsch   PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&D));
2129566063dSJacob Faibussowitsch   PetscCall(MatSetOption(D,MAT_SYMMETRIC,symm));
2139566063dSJacob Faibussowitsch   PetscCall(MatMultEqual(B,D,10,&flg));
214c4fb69feSStefano Zampini   if (!flg) {
2159566063dSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMult error after MatDuplicate\n"));
216c4fb69feSStefano Zampini   }
217c4fb69feSStefano Zampini   if (testtrans) {
2189566063dSJacob Faibussowitsch     PetscCall(MatMultTransposeEqual(B,D,10,&flg));
219c4fb69feSStefano Zampini     if (!flg) {
2209566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMultTranspose error after MatDuplicate\n"));
221c4fb69feSStefano Zampini     }
222c4fb69feSStefano Zampini   }
2239566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&D));
224c4fb69feSStefano Zampini 
225c4fb69feSStefano Zampini   if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */
2269566063dSJacob Faibussowitsch     PetscCall(VecSetRandom(Ay,NULL));
2279566063dSJacob Faibussowitsch     PetscCall(VecCopy(Ay,By));
2289566063dSJacob Faibussowitsch     PetscCall(MatMultTranspose(A,Ay,Ax));
2299566063dSJacob Faibussowitsch     PetscCall(MatMultTranspose(B,By,Bx));
2309566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(Ax,NULL,"-multtrans_vec_view"));
2319566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(Bx,NULL,"-multtrans_vec_view"));
2329566063dSJacob Faibussowitsch     PetscCall(VecNorm(Ax,NORM_INFINITY,&nX));
2339566063dSJacob Faibussowitsch     PetscCall(VecAXPY(Ax,-1.0,Bx));
2349566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(Ax,NULL,"-multtrans_vec_view"));
2359566063dSJacob Faibussowitsch     PetscCall(VecNorm(Ax,NORM_INFINITY,&err));
2369566063dSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMultTranspose err %g\n",err/nX));
2379566063dSJacob Faibussowitsch     PetscCall(VecScale(Bx,-1.0));
2389566063dSJacob Faibussowitsch     PetscCall(MatMultTransposeAdd(B,By,Bx,Bx));
2399566063dSJacob Faibussowitsch     PetscCall(VecNorm(Bx,NORM_INFINITY,&err));
2404c3a5f3cSStefano Zampini     if (err > 10.*PETSC_SMALL) {
2419566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMultTransposeAdd err %g\n",err));
242c4fb69feSStefano Zampini     }
243c4fb69feSStefano Zampini   }
2449566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&Ax));
2459566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&Ay));
2469566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&Bx));
2479566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&By));
248c4fb69feSStefano Zampini 
249c4fb69feSStefano Zampini   /* Test MatMatMult */
250c4fb69feSStefano Zampini   if (ldc) {
2519566063dSJacob Faibussowitsch     PetscCall(PetscMalloc1(nrhs*(n+ldc),&Cdata));
252c4fb69feSStefano Zampini   }
2539566063dSJacob Faibussowitsch   PetscCall(MatCreateDense(PETSC_COMM_WORLD,n,PETSC_DECIDE,N,nrhs,Cdata,&C));
2549566063dSJacob Faibussowitsch   PetscCall(MatDenseSetLDA(C,n+ldc));
2559566063dSJacob Faibussowitsch   PetscCall(MatSetRandom(C,NULL));
256c4fb69feSStefano Zampini   if (cgpu) {
2579566063dSJacob Faibussowitsch     PetscCall(MatConvert(C,MATDENSECUDA,MAT_INPLACE_MATRIX,&C));
258c4fb69feSStefano Zampini   }
259c4fb69feSStefano Zampini   for (nt = 0; nt < ntrials; nt++) {
2609566063dSJacob Faibussowitsch     PetscCall(MatMatMult(B,C,nt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D));
2619566063dSJacob Faibussowitsch     PetscCall(MatViewFromOptions(D,NULL,"-bc_view"));
2629566063dSJacob Faibussowitsch     PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)D,&flg,MATSEQDENSE,MATMPIDENSE,""));
263c4fb69feSStefano Zampini     if (flg) {
2649566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(B,&x,&y));
2659566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(D,NULL,&v));
266c4fb69feSStefano Zampini       for (i = 0; i < nrhs; i++) {
2679566063dSJacob Faibussowitsch         PetscCall(MatGetColumnVector(D,v,i));
2689566063dSJacob Faibussowitsch         PetscCall(MatGetColumnVector(C,x,i));
2699566063dSJacob Faibussowitsch         PetscCall(MatMult(B,x,y));
2709566063dSJacob Faibussowitsch         PetscCall(VecAXPY(y,-1.0,v));
2719566063dSJacob Faibussowitsch         PetscCall(VecNorm(y,NORM_INFINITY,&err));
2729566063dSJacob Faibussowitsch         if (err > 10.*PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMat err %" PetscInt_FMT " %g\n",i,err));
273c4fb69feSStefano Zampini       }
2749566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&y));
2759566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&x));
2769566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&v));
277c4fb69feSStefano Zampini     }
278c4fb69feSStefano Zampini   }
2799566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&D));
280c4fb69feSStefano Zampini 
281c4fb69feSStefano Zampini   /* Test MatTransposeMatMult */
282c4fb69feSStefano Zampini   if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */
283c4fb69feSStefano Zampini     for (nt = 0; nt < ntrials; nt++) {
2849566063dSJacob Faibussowitsch       PetscCall(MatTransposeMatMult(B,C,nt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D));
2859566063dSJacob Faibussowitsch       PetscCall(MatViewFromOptions(D,NULL,"-btc_view"));
2869566063dSJacob Faibussowitsch       PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)D,&flg,MATSEQDENSE,MATMPIDENSE,""));
287c4fb69feSStefano Zampini       if (flg) {
2889566063dSJacob Faibussowitsch         PetscCall(MatCreateVecs(B,&y,&x));
2899566063dSJacob Faibussowitsch         PetscCall(MatCreateVecs(D,NULL,&v));
290c4fb69feSStefano Zampini         for (i = 0; i < nrhs; i++) {
2919566063dSJacob Faibussowitsch           PetscCall(MatGetColumnVector(D,v,i));
2929566063dSJacob Faibussowitsch           PetscCall(MatGetColumnVector(C,x,i));
2939566063dSJacob Faibussowitsch           PetscCall(MatMultTranspose(B,x,y));
2949566063dSJacob Faibussowitsch           PetscCall(VecAXPY(y,-1.0,v));
2959566063dSJacob Faibussowitsch           PetscCall(VecNorm(y,NORM_INFINITY,&err));
2969566063dSJacob Faibussowitsch           if (err > 10.*PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatTransMat err %" PetscInt_FMT " %g\n",i,err));
297c4fb69feSStefano Zampini         }
2989566063dSJacob Faibussowitsch         PetscCall(VecDestroy(&y));
2999566063dSJacob Faibussowitsch         PetscCall(VecDestroy(&x));
3009566063dSJacob Faibussowitsch         PetscCall(VecDestroy(&v));
301c4fb69feSStefano Zampini       }
302c4fb69feSStefano Zampini     }
3039566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&D));
304c4fb69feSStefano Zampini   }
305c4fb69feSStefano Zampini 
306300d917bSStefano Zampini   /* Test basis orthogonalization */
30753022affSStefano Zampini   if (testorthog) {
3089566063dSJacob Faibussowitsch     PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&D));
3099566063dSJacob Faibussowitsch     PetscCall(MatSetOption(D,MAT_SYMMETRIC,symm));
3109566063dSJacob Faibussowitsch     PetscCall(MatH2OpusOrthogonalize(D));
3119566063dSJacob Faibussowitsch     PetscCall(MatMultEqual(B,D,10,&flg));
31253022affSStefano Zampini     if (!flg) {
3139566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMult error after basis ortogonalization\n"));
31453022affSStefano Zampini     }
3159566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&D));
31653022affSStefano Zampini   }
31753022affSStefano Zampini 
318300d917bSStefano Zampini   /* Test matrix compression */
31953022affSStefano Zampini   if (testcompress) {
3209566063dSJacob Faibussowitsch     PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&D));
3219566063dSJacob Faibussowitsch     PetscCall(MatSetOption(D,MAT_SYMMETRIC,symm));
3229566063dSJacob Faibussowitsch     PetscCall(MatH2OpusCompress(D,PETSC_SMALL));
3239566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&D));
32453022affSStefano Zampini   }
32553022affSStefano Zampini 
326300d917bSStefano Zampini   /* Test low-rank update */
327300d917bSStefano Zampini   if (testhlru) {
328300d917bSStefano Zampini     Mat         U, V;
329300d917bSStefano Zampini     PetscScalar *Udata = NULL, *Vdata = NULL;
330300d917bSStefano Zampini 
331300d917bSStefano Zampini     if (ldu) {
3329566063dSJacob Faibussowitsch       PetscCall(PetscMalloc1(nlr*(n+ldu),&Udata));
3339566063dSJacob Faibussowitsch       PetscCall(PetscMalloc1(nlr*(n+ldu+2),&Vdata));
334300d917bSStefano Zampini     }
3359566063dSJacob Faibussowitsch     PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&D));
3369566063dSJacob Faibussowitsch     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)D),n,PETSC_DECIDE,N,nlr,Udata,&U));
3379566063dSJacob Faibussowitsch     PetscCall(MatDenseSetLDA(U,n+ldu));
3389566063dSJacob Faibussowitsch     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)D),n,PETSC_DECIDE,N,nlr,Vdata,&V));
3399566063dSJacob Faibussowitsch     if (ldu) PetscCall(MatDenseSetLDA(V,n+ldu+2));
3409566063dSJacob Faibussowitsch     PetscCall(MatSetRandom(U,NULL));
3419566063dSJacob Faibussowitsch     PetscCall(MatSetRandom(V,NULL));
3429566063dSJacob Faibussowitsch     PetscCall(MatH2OpusLowRankUpdate(D,U,V,0.5));
3439566063dSJacob Faibussowitsch     PetscCall(MatH2OpusLowRankUpdate(D,U,V,-0.5));
3449566063dSJacob Faibussowitsch     PetscCall(MatMultEqual(B,D,10,&flg));
345300d917bSStefano Zampini     if (!flg) {
3469566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMult error after low-rank update\n"));
347300d917bSStefano Zampini     }
3489566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&D));
3499566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&U));
3509566063dSJacob Faibussowitsch     PetscCall(PetscFree(Udata));
3519566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&V));
3529566063dSJacob Faibussowitsch     PetscCall(PetscFree(Vdata));
353300d917bSStefano Zampini   }
354300d917bSStefano Zampini 
355c4fb69feSStefano Zampini   /* check explicit operator */
356c4fb69feSStefano Zampini   if (checkexpl) {
357c4fb69feSStefano Zampini     Mat Be, Bet;
358c4fb69feSStefano Zampini 
3599566063dSJacob Faibussowitsch     PetscCall(MatComputeOperator(B,MATDENSE,&D));
3609566063dSJacob Faibussowitsch     PetscCall(MatDuplicate(D,MAT_COPY_VALUES,&Be));
3619566063dSJacob Faibussowitsch     PetscCall(MatNorm(D,NORM_FROBENIUS,&nB));
3629566063dSJacob Faibussowitsch     PetscCall(MatViewFromOptions(D,NULL,"-expl_view"));
3639566063dSJacob Faibussowitsch     PetscCall(MatAXPY(D,-1.0,A,SAME_NONZERO_PATTERN));
3649566063dSJacob Faibussowitsch     PetscCall(MatViewFromOptions(D,NULL,"-diff_view"));
3659566063dSJacob Faibussowitsch     PetscCall(MatNorm(D,NORM_FROBENIUS,&nD));
3669566063dSJacob Faibussowitsch     PetscCall(MatNorm(A,NORM_FROBENIUS,&nA));
3679566063dSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Approximation error %g (%g / %g, %g)\n",nD/nA,nD,nA,nB));
3689566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&D));
369c4fb69feSStefano Zampini 
370c4fb69feSStefano Zampini     if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */
3719566063dSJacob Faibussowitsch       PetscCall(MatTranspose(A,MAT_INPLACE_MATRIX,&A));
3729566063dSJacob Faibussowitsch       PetscCall(MatComputeOperatorTranspose(B,MATDENSE,&D));
3739566063dSJacob Faibussowitsch       PetscCall(MatDuplicate(D,MAT_COPY_VALUES,&Bet));
3749566063dSJacob Faibussowitsch       PetscCall(MatNorm(D,NORM_FROBENIUS,&nB));
3759566063dSJacob Faibussowitsch       PetscCall(MatViewFromOptions(D,NULL,"-expl_trans_view"));
3769566063dSJacob Faibussowitsch       PetscCall(MatAXPY(D,-1.0,A,SAME_NONZERO_PATTERN));
3779566063dSJacob Faibussowitsch       PetscCall(MatViewFromOptions(D,NULL,"-diff_trans_view"));
3789566063dSJacob Faibussowitsch       PetscCall(MatNorm(D,NORM_FROBENIUS,&nD));
3799566063dSJacob Faibussowitsch       PetscCall(MatNorm(A,NORM_FROBENIUS,&nA));
3809566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Approximation error transpose %g (%g / %g, %g)\n",nD/nA,nD,nA,nB));
3819566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&D));
382c4fb69feSStefano Zampini 
3839566063dSJacob Faibussowitsch       PetscCall(MatTranspose(Bet,MAT_INPLACE_MATRIX,&Bet));
3849566063dSJacob Faibussowitsch       PetscCall(MatAXPY(Be,-1.0,Bet,SAME_NONZERO_PATTERN));
3859566063dSJacob Faibussowitsch       PetscCall(MatViewFromOptions(Be,NULL,"-diff_expl_view"));
3869566063dSJacob Faibussowitsch       PetscCall(MatNorm(Be,NORM_FROBENIUS,&nB));
3879566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Approximation error B - (B^T)^T %g\n",nB));
3889566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&Be));
3899566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&Bet));
390c4fb69feSStefano Zampini     }
391c4fb69feSStefano Zampini   }
3929566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&A));
3939566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&B));
3949566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&C));
3959566063dSJacob Faibussowitsch   PetscCall(PetscFree(Cdata));
3969566063dSJacob Faibussowitsch   PetscCall(PetscFree(Adata));
3979566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
398b122ec5aSJacob Faibussowitsch   return 0;
399c4fb69feSStefano Zampini }
400c4fb69feSStefano Zampini 
401c4fb69feSStefano Zampini /*TEST
402c4fb69feSStefano Zampini 
403c4fb69feSStefano Zampini    build:
40453022affSStefano Zampini      requires: h2opus
405c4fb69feSStefano Zampini 
406c4fb69feSStefano Zampini #tests from kernel
407c4fb69feSStefano Zampini    test:
40853022affSStefano Zampini      requires: h2opus
409c4fb69feSStefano Zampini      nsize: 1
410c4fb69feSStefano Zampini      suffix: 1
411c4fb69feSStefano Zampini      args: -n {{17 33}} -kernel 1 -dim {{1 2 3}} -symm {{0 1}} -checkexpl -bgpu 0
412c4fb69feSStefano Zampini 
413c4fb69feSStefano Zampini    test:
41453022affSStefano Zampini      requires: h2opus
415c4fb69feSStefano Zampini      nsize: 1
416c4fb69feSStefano Zampini      suffix: 1_ld
417c4fb69feSStefano Zampini      output_file: output/ex66_1.out
418300d917bSStefano Zampini      args: -n 33 -kernel 1 -dim 1 -lda 13 -ldc 11 -ldu 17 -symm 0 -checkexpl -bgpu 0
419c4fb69feSStefano Zampini 
420c4fb69feSStefano Zampini    test:
42153022affSStefano Zampini      requires: h2opus cuda
422c4fb69feSStefano Zampini      nsize: 1
423c4fb69feSStefano Zampini      suffix: 1_cuda
424c4fb69feSStefano Zampini      output_file: output/ex66_1.out
425c4fb69feSStefano Zampini      args: -n {{17 33}} -kernel 1 -dim {{1 2 3}} -symm {{0 1}} -checkexpl -bgpu 1
426c4fb69feSStefano Zampini 
427c4fb69feSStefano Zampini    test:
42853022affSStefano Zampini      requires: h2opus cuda
429c4fb69feSStefano Zampini      nsize: 1
430c4fb69feSStefano Zampini      suffix: 1_cuda_ld
431c4fb69feSStefano Zampini      output_file: output/ex66_1.out
432300d917bSStefano Zampini      args: -n 33 -kernel 1 -dim 1 -lda 13 -ldc 11 -ldu 17 -symm 0 -checkexpl -bgpu 1
433c4fb69feSStefano Zampini 
434c4fb69feSStefano Zampini    test:
43553022affSStefano Zampini      requires: h2opus
43653022affSStefano Zampini      nsize: {{2 3}}
437c4fb69feSStefano Zampini      suffix: 1_par
43853022affSStefano Zampini      args: -n 64 -symm -kernel 1 -dim 1 -ldc 12 -testlayout {{0 1}} -bgpu 0 -cgpu 0
439c4fb69feSStefano Zampini 
440c4fb69feSStefano Zampini    test:
44153022affSStefano Zampini      requires: h2opus cuda
44253022affSStefano Zampini      nsize: {{2 3}}
443c4fb69feSStefano Zampini      suffix: 1_par_cuda
44453022affSStefano Zampini      args: -n 64 -symm -kernel 1 -dim 1 -ldc 12 -testlayout {{0 1}} -bgpu {{0 1}} -cgpu {{0 1}}
445c4fb69feSStefano Zampini      output_file: output/ex66_1_par.out
446c4fb69feSStefano Zampini 
447c4fb69feSStefano Zampini #tests from matrix sampling (parallel or unsymmetric not supported)
448c4fb69feSStefano Zampini    test:
44953022affSStefano Zampini      requires: h2opus
450c4fb69feSStefano Zampini      nsize: 1
451c4fb69feSStefano Zampini      suffix: 2
452c4fb69feSStefano Zampini      args: -n {{17 33}} -kernel 0 -dim 2 -symm 1 -checkexpl -bgpu 0
453c4fb69feSStefano Zampini 
454c4fb69feSStefano Zampini    test:
45553022affSStefano Zampini      requires: h2opus cuda
456c4fb69feSStefano Zampini      nsize: 1
457c4fb69feSStefano Zampini      suffix: 2_cuda
458c4fb69feSStefano Zampini      output_file: output/ex66_2.out
459c4fb69feSStefano Zampini      args: -n {{17 33}} -kernel 0 -dim 2 -symm 1 -checkexpl -bgpu {{0 1}} -agpu {{0 1}}
460c4fb69feSStefano Zampini 
46153022affSStefano Zampini #tests view operation
46253022affSStefano Zampini    test:
46353022affSStefano Zampini      requires: h2opus !cuda
4648cc725e6SPierre Jolivet      filter: grep -v " MPI process" | grep -v "\[" | grep -v "\]"
46553022affSStefano Zampini      nsize: {{1 2 3}}
46653022affSStefano Zampini      suffix: view
46753022affSStefano Zampini      args: -ng 64 -kernel 1 -dim 2 -symm 1 -checkexpl -B_view -mat_h2opus_leafsize 17 -mat_h2opus_normsamples 13 -mat_h2opus_indexmap_view ::ascii_matlab -mat_approximate_norm_samples 2 -mat_h2opus_normsamples 2
46853022affSStefano Zampini 
46953022affSStefano Zampini    test:
47053022affSStefano Zampini      requires: h2opus cuda
4718cc725e6SPierre Jolivet      filter: grep -v " MPI process" | grep -v "\[" | grep -v "\]"
47253022affSStefano Zampini      nsize: {{1 2 3}}
47353022affSStefano Zampini      suffix: view_cuda
47453022affSStefano Zampini      args: -ng 64 -kernel 1 -dim 2 -symm 1 -checkexpl -bgpu -B_view -mat_h2opus_leafsize 17 -mat_h2opus_normsamples 13 -mat_h2opus_indexmap_view ::ascii_matlab -mat_approximate_norm_samples 2 -mat_h2opus_normsamples 2
47553022affSStefano Zampini 
476c4fb69feSStefano Zampini TEST*/
477