xref: /petsc/src/mat/tests/ex66.c (revision 9566063d113dddea24716c546802770db7481bc0)
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*9566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc,&argv,(char*) 0,help));
50*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-ng",&N,&flgglob));
51*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
52*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL));
53*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL));
54*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-lda",&lda,NULL));
55*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-ldc",&ldc,NULL));
56*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-nlr",&nlr,NULL));
57*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-ldu",&ldu,NULL));
58*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL,NULL,"-matmattrials",&ntrials,NULL));
59*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-randommat",&randommat,NULL));
60*9566063dSJacob Faibussowitsch   if (!flgglob) PetscCall(PetscOptionsGetBool(NULL,NULL,"-testlayout",&testlayout,NULL));
61*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-Asymm",&Asymm,NULL));
62*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-symm",&symm,NULL));
63*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-kernel",&kernel,NULL));
64*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-checkexpl",&checkexpl,NULL));
65*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-agpu",&agpu,NULL));
66*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-bgpu",&bgpu,NULL));
67*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,NULL,"-cgpu",&cgpu,NULL));
68*9566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-scale",&scale,NULL));
69c4fb69feSStefano Zampini   if (!Asymm) symm = PETSC_FALSE;
70c4fb69feSStefano Zampini 
71*9566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
72300d917bSStefano Zampini 
73300d917bSStefano Zampini   /* Disable tests for unimplemented variants */
74c4fb69feSStefano Zampini   testtrans = (PetscBool)(size == 1 || symm);
75c4fb69feSStefano Zampini   testnorm = (PetscBool)(size == 1 || symm);
7653022affSStefano Zampini   testorthog = (PetscBool)(size == 1 || symm);
7753022affSStefano Zampini   testcompress = (PetscBool)(size == 1 || symm);
78300d917bSStefano Zampini   testhlru = (PetscBool)(size == 1);
79c4fb69feSStefano Zampini 
80*9566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
81*9566063dSJacob Faibussowitsch   PetscCall(PetscLayoutCreate(PETSC_COMM_WORLD,&map));
82c4fb69feSStefano Zampini   if (testlayout) {
83c4fb69feSStefano Zampini     if (rank%2) n = PetscMax(2*n-5*rank,0);
84c4fb69feSStefano Zampini     else n = 2*n+rank;
85c4fb69feSStefano Zampini   }
8653022affSStefano Zampini   if (!flgglob) {
87*9566063dSJacob Faibussowitsch     PetscCall(PetscLayoutSetLocalSize(map,n));
88*9566063dSJacob Faibussowitsch     PetscCall(PetscLayoutSetUp(map));
89*9566063dSJacob Faibussowitsch     PetscCall(PetscLayoutGetSize(map,&N));
9053022affSStefano Zampini   } else {
91*9566063dSJacob Faibussowitsch     PetscCall(PetscLayoutSetSize(map,N));
92*9566063dSJacob Faibussowitsch     PetscCall(PetscLayoutSetUp(map));
93*9566063dSJacob Faibussowitsch     PetscCall(PetscLayoutGetLocalSize(map,&n));
9453022affSStefano Zampini   }
95*9566063dSJacob Faibussowitsch   PetscCall(PetscLayoutDestroy(&map));
96c4fb69feSStefano Zampini 
97c4fb69feSStefano Zampini   if (lda) {
98*9566063dSJacob Faibussowitsch     PetscCall(PetscMalloc1(N*(n+lda),&Adata));
99c4fb69feSStefano Zampini   }
100*9566063dSJacob Faibussowitsch   PetscCall(MatCreateDense(PETSC_COMM_WORLD,n,n,N,N,Adata,&A));
101*9566063dSJacob Faibussowitsch   PetscCall(MatDenseSetLDA(A,n+lda));
102c4fb69feSStefano Zampini 
10353022affSStefano Zampini   /* Create random points; these are replicated in order to populate a dense matrix and to compare sequential and dense runs
10453022affSStefano Zampini      The constructor for MATH2OPUS can take as input the distributed coordinates and replicates them internally in case
10553022affSStefano Zampini      a kernel construction is requested */
106*9566063dSJacob Faibussowitsch   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&r));
107*9566063dSJacob Faibussowitsch   PetscCall(PetscRandomSetFromOptions(r));
108*9566063dSJacob Faibussowitsch   PetscCall(PetscRandomSetSeed(r,123456));
109*9566063dSJacob Faibussowitsch   PetscCall(PetscRandomSeed(r));
110*9566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(N*dim,&coords));
111*9566063dSJacob Faibussowitsch   PetscCall(PetscRandomGetValuesReal(r,N*dim,coords));
112*9566063dSJacob Faibussowitsch   PetscCall(PetscRandomDestroy(&r));
11353022affSStefano Zampini 
114c4fb69feSStefano Zampini   if (kernel || !randommat) {
11553022affSStefano Zampini     MatH2OpusKernel k = Asymm ? GenEntry_Symm : GenEntry_Unsymm;
116c4fb69feSStefano Zampini     PetscInt        ist,ien;
117c4fb69feSStefano Zampini 
118*9566063dSJacob Faibussowitsch     PetscCall(MatGetOwnershipRange(A,&ist,&ien));
119c4fb69feSStefano Zampini     for (i = ist; i < ien; i++) {
120c4fb69feSStefano Zampini       for (j = 0; j < N; j++) {
121*9566063dSJacob Faibussowitsch         PetscCall(MatSetValue(A,i,j,(*k)(dim,coords + i*dim,coords + j*dim,NULL),INSERT_VALUES));
122c4fb69feSStefano Zampini       }
123c4fb69feSStefano Zampini     }
124*9566063dSJacob Faibussowitsch     PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
125*9566063dSJacob Faibussowitsch     PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
126c4fb69feSStefano Zampini     if (kernel) {
127*9566063dSJacob 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));
128c4fb69feSStefano Zampini     } else {
129*9566063dSJacob Faibussowitsch       PetscCall(MatCreateH2OpusFromMat(A,dim,coords + ist*dim,PETSC_TRUE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,&B));
130c4fb69feSStefano Zampini     }
131c4fb69feSStefano Zampini   } else {
13253022affSStefano Zampini     PetscInt ist;
13353022affSStefano Zampini 
134*9566063dSJacob Faibussowitsch     PetscCall(MatGetOwnershipRange(A,&ist,NULL));
135*9566063dSJacob Faibussowitsch     PetscCall(MatSetRandom(A,NULL));
136c4fb69feSStefano Zampini     if (Asymm) {
137*9566063dSJacob Faibussowitsch       PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX,&B));
138*9566063dSJacob Faibussowitsch       PetscCall(MatAXPY(A,1.0,B,SAME_NONZERO_PATTERN));
139*9566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&B));
140*9566063dSJacob Faibussowitsch       PetscCall(MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE));
141c4fb69feSStefano Zampini     }
142*9566063dSJacob Faibussowitsch     PetscCall(MatCreateH2OpusFromMat(A,dim,coords + ist*dim,PETSC_TRUE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,&B));
143c4fb69feSStefano Zampini   }
144*9566063dSJacob Faibussowitsch   PetscCall(PetscFree(coords));
145c4fb69feSStefano Zampini   if (agpu) {
146*9566063dSJacob Faibussowitsch     PetscCall(MatConvert(A,MATDENSECUDA,MAT_INPLACE_MATRIX,&A));
147c4fb69feSStefano Zampini   }
148*9566063dSJacob Faibussowitsch   PetscCall(MatViewFromOptions(A,NULL,"-A_view"));
149c4fb69feSStefano Zampini 
150*9566063dSJacob Faibussowitsch   PetscCall(MatSetOption(B,MAT_SYMMETRIC,symm));
151c4fb69feSStefano Zampini 
152c4fb69feSStefano Zampini   /* assemble the H-matrix */
153*9566063dSJacob Faibussowitsch   PetscCall(MatBindToCPU(B,(PetscBool)!bgpu));
154*9566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(B));
155*9566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
156*9566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
157*9566063dSJacob Faibussowitsch   PetscCall(MatViewFromOptions(B,NULL,"-B_view"));
158c4fb69feSStefano Zampini 
159c4fb69feSStefano Zampini   /* Test MatScale */
160*9566063dSJacob Faibussowitsch   PetscCall(MatScale(A,scale));
161*9566063dSJacob Faibussowitsch   PetscCall(MatScale(B,scale));
162c4fb69feSStefano Zampini 
163c4fb69feSStefano Zampini   /* Test MatMult */
164*9566063dSJacob Faibussowitsch   PetscCall(MatCreateVecs(A,&Ax,&Ay));
165*9566063dSJacob Faibussowitsch   PetscCall(MatCreateVecs(B,&Bx,&By));
166*9566063dSJacob Faibussowitsch   PetscCall(VecSetRandom(Ax,NULL));
167*9566063dSJacob Faibussowitsch   PetscCall(VecCopy(Ax,Bx));
168*9566063dSJacob Faibussowitsch   PetscCall(MatMult(A,Ax,Ay));
169*9566063dSJacob Faibussowitsch   PetscCall(MatMult(B,Bx,By));
170*9566063dSJacob Faibussowitsch   PetscCall(VecViewFromOptions(Ay,NULL,"-mult_vec_view"));
171*9566063dSJacob Faibussowitsch   PetscCall(VecViewFromOptions(By,NULL,"-mult_vec_view"));
172*9566063dSJacob Faibussowitsch   PetscCall(VecNorm(Ay,NORM_INFINITY,&nX));
173*9566063dSJacob Faibussowitsch   PetscCall(VecAXPY(Ay,-1.0,By));
174*9566063dSJacob Faibussowitsch   PetscCall(VecViewFromOptions(Ay,NULL,"-mult_vec_view"));
175*9566063dSJacob Faibussowitsch   PetscCall(VecNorm(Ay,NORM_INFINITY,&err));
176*9566063dSJacob Faibussowitsch   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMult err %g\n",err/nX));
177*9566063dSJacob Faibussowitsch   PetscCall(VecScale(By,-1.0));
178*9566063dSJacob Faibussowitsch   PetscCall(MatMultAdd(B,Bx,By,By));
179*9566063dSJacob Faibussowitsch   PetscCall(VecNorm(By,NORM_INFINITY,&err));
180*9566063dSJacob Faibussowitsch   PetscCall(VecViewFromOptions(By,NULL,"-mult_vec_view"));
1814c3a5f3cSStefano Zampini   if (err > 10.*PETSC_SMALL) {
182*9566063dSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMultAdd err %g\n",err));
183c4fb69feSStefano Zampini   }
184c4fb69feSStefano Zampini 
185c4fb69feSStefano Zampini   /* Test MatNorm */
186*9566063dSJacob Faibussowitsch   PetscCall(MatNorm(A,NORM_INFINITY,&norms[0]));
187*9566063dSJacob Faibussowitsch   PetscCall(MatNorm(A,NORM_1,&norms[1]));
188c4fb69feSStefano Zampini   norms[2] = -1.; /* NORM_2 not supported */
189*9566063dSJacob 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]));
190*9566063dSJacob Faibussowitsch   PetscCall(MatGetOperation(A,MATOP_NORM,&Anormfunc));
191*9566063dSJacob Faibussowitsch   PetscCall(MatGetOperation(B,MATOP_NORM,&approxnormfunc));
192*9566063dSJacob Faibussowitsch   PetscCall(MatSetOperation(A,MATOP_NORM,approxnormfunc));
193*9566063dSJacob Faibussowitsch   PetscCall(MatNorm(A,NORM_INFINITY,&norms[0]));
194*9566063dSJacob Faibussowitsch   PetscCall(MatNorm(A,NORM_1,&norms[1]));
195*9566063dSJacob Faibussowitsch   PetscCall(MatNorm(A,NORM_2,&norms[2]));
196*9566063dSJacob 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]));
197c4fb69feSStefano Zampini   if (testnorm) {
198*9566063dSJacob Faibussowitsch     PetscCall(MatNorm(B,NORM_INFINITY,&norms[0]));
199*9566063dSJacob Faibussowitsch     PetscCall(MatNorm(B,NORM_1,&norms[1]));
200*9566063dSJacob Faibussowitsch     PetscCall(MatNorm(B,NORM_2,&norms[2]));
201c4fb69feSStefano Zampini   } else {
202c4fb69feSStefano Zampini     norms[0] = -1.;
203c4fb69feSStefano Zampini     norms[1] = -1.;
204c4fb69feSStefano Zampini     norms[2] = -1.;
205c4fb69feSStefano Zampini   }
206*9566063dSJacob 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]));
207*9566063dSJacob Faibussowitsch   PetscCall(MatSetOperation(A,MATOP_NORM,Anormfunc));
208c4fb69feSStefano Zampini 
209c4fb69feSStefano Zampini   /* Test MatDuplicate */
210*9566063dSJacob Faibussowitsch   PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&D));
211*9566063dSJacob Faibussowitsch   PetscCall(MatSetOption(D,MAT_SYMMETRIC,symm));
212*9566063dSJacob Faibussowitsch   PetscCall(MatMultEqual(B,D,10,&flg));
213c4fb69feSStefano Zampini   if (!flg) {
214*9566063dSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMult error after MatDuplicate\n"));
215c4fb69feSStefano Zampini   }
216c4fb69feSStefano Zampini   if (testtrans) {
217*9566063dSJacob Faibussowitsch     PetscCall(MatMultTransposeEqual(B,D,10,&flg));
218c4fb69feSStefano Zampini     if (!flg) {
219*9566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMultTranspose error after MatDuplicate\n"));
220c4fb69feSStefano Zampini     }
221c4fb69feSStefano Zampini   }
222*9566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&D));
223c4fb69feSStefano Zampini 
224c4fb69feSStefano Zampini   if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */
225*9566063dSJacob Faibussowitsch     PetscCall(VecSetRandom(Ay,NULL));
226*9566063dSJacob Faibussowitsch     PetscCall(VecCopy(Ay,By));
227*9566063dSJacob Faibussowitsch     PetscCall(MatMultTranspose(A,Ay,Ax));
228*9566063dSJacob Faibussowitsch     PetscCall(MatMultTranspose(B,By,Bx));
229*9566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(Ax,NULL,"-multtrans_vec_view"));
230*9566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(Bx,NULL,"-multtrans_vec_view"));
231*9566063dSJacob Faibussowitsch     PetscCall(VecNorm(Ax,NORM_INFINITY,&nX));
232*9566063dSJacob Faibussowitsch     PetscCall(VecAXPY(Ax,-1.0,Bx));
233*9566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(Ax,NULL,"-multtrans_vec_view"));
234*9566063dSJacob Faibussowitsch     PetscCall(VecNorm(Ax,NORM_INFINITY,&err));
235*9566063dSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMultTranspose err %g\n",err/nX));
236*9566063dSJacob Faibussowitsch     PetscCall(VecScale(Bx,-1.0));
237*9566063dSJacob Faibussowitsch     PetscCall(MatMultTransposeAdd(B,By,Bx,Bx));
238*9566063dSJacob Faibussowitsch     PetscCall(VecNorm(Bx,NORM_INFINITY,&err));
2394c3a5f3cSStefano Zampini     if (err > 10.*PETSC_SMALL) {
240*9566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMultTransposeAdd err %g\n",err));
241c4fb69feSStefano Zampini     }
242c4fb69feSStefano Zampini   }
243*9566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&Ax));
244*9566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&Ay));
245*9566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&Bx));
246*9566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&By));
247c4fb69feSStefano Zampini 
248c4fb69feSStefano Zampini   /* Test MatMatMult */
249c4fb69feSStefano Zampini   if (ldc) {
250*9566063dSJacob Faibussowitsch     PetscCall(PetscMalloc1(nrhs*(n+ldc),&Cdata));
251c4fb69feSStefano Zampini   }
252*9566063dSJacob Faibussowitsch   PetscCall(MatCreateDense(PETSC_COMM_WORLD,n,PETSC_DECIDE,N,nrhs,Cdata,&C));
253*9566063dSJacob Faibussowitsch   PetscCall(MatDenseSetLDA(C,n+ldc));
254*9566063dSJacob Faibussowitsch   PetscCall(MatSetRandom(C,NULL));
255c4fb69feSStefano Zampini   if (cgpu) {
256*9566063dSJacob Faibussowitsch     PetscCall(MatConvert(C,MATDENSECUDA,MAT_INPLACE_MATRIX,&C));
257c4fb69feSStefano Zampini   }
258c4fb69feSStefano Zampini   for (nt = 0; nt < ntrials; nt++) {
259*9566063dSJacob Faibussowitsch     PetscCall(MatMatMult(B,C,nt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D));
260*9566063dSJacob Faibussowitsch     PetscCall(MatViewFromOptions(D,NULL,"-bc_view"));
261*9566063dSJacob Faibussowitsch     PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)D,&flg,MATSEQDENSE,MATMPIDENSE,""));
262c4fb69feSStefano Zampini     if (flg) {
263*9566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(B,&x,&y));
264*9566063dSJacob Faibussowitsch       PetscCall(MatCreateVecs(D,NULL,&v));
265c4fb69feSStefano Zampini       for (i = 0; i < nrhs; i++) {
266*9566063dSJacob Faibussowitsch         PetscCall(MatGetColumnVector(D,v,i));
267*9566063dSJacob Faibussowitsch         PetscCall(MatGetColumnVector(C,x,i));
268*9566063dSJacob Faibussowitsch         PetscCall(MatMult(B,x,y));
269*9566063dSJacob Faibussowitsch         PetscCall(VecAXPY(y,-1.0,v));
270*9566063dSJacob Faibussowitsch         PetscCall(VecNorm(y,NORM_INFINITY,&err));
271*9566063dSJacob Faibussowitsch         if (err > 10.*PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMat err %" PetscInt_FMT " %g\n",i,err));
272c4fb69feSStefano Zampini       }
273*9566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&y));
274*9566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&x));
275*9566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&v));
276c4fb69feSStefano Zampini     }
277c4fb69feSStefano Zampini   }
278*9566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&D));
279c4fb69feSStefano Zampini 
280c4fb69feSStefano Zampini   /* Test MatTransposeMatMult */
281c4fb69feSStefano Zampini   if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */
282c4fb69feSStefano Zampini     for (nt = 0; nt < ntrials; nt++) {
283*9566063dSJacob Faibussowitsch       PetscCall(MatTransposeMatMult(B,C,nt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D));
284*9566063dSJacob Faibussowitsch       PetscCall(MatViewFromOptions(D,NULL,"-btc_view"));
285*9566063dSJacob Faibussowitsch       PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)D,&flg,MATSEQDENSE,MATMPIDENSE,""));
286c4fb69feSStefano Zampini       if (flg) {
287*9566063dSJacob Faibussowitsch         PetscCall(MatCreateVecs(B,&y,&x));
288*9566063dSJacob Faibussowitsch         PetscCall(MatCreateVecs(D,NULL,&v));
289c4fb69feSStefano Zampini         for (i = 0; i < nrhs; i++) {
290*9566063dSJacob Faibussowitsch           PetscCall(MatGetColumnVector(D,v,i));
291*9566063dSJacob Faibussowitsch           PetscCall(MatGetColumnVector(C,x,i));
292*9566063dSJacob Faibussowitsch           PetscCall(MatMultTranspose(B,x,y));
293*9566063dSJacob Faibussowitsch           PetscCall(VecAXPY(y,-1.0,v));
294*9566063dSJacob Faibussowitsch           PetscCall(VecNorm(y,NORM_INFINITY,&err));
295*9566063dSJacob Faibussowitsch           if (err > 10.*PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatTransMat err %" PetscInt_FMT " %g\n",i,err));
296c4fb69feSStefano Zampini         }
297*9566063dSJacob Faibussowitsch         PetscCall(VecDestroy(&y));
298*9566063dSJacob Faibussowitsch         PetscCall(VecDestroy(&x));
299*9566063dSJacob Faibussowitsch         PetscCall(VecDestroy(&v));
300c4fb69feSStefano Zampini       }
301c4fb69feSStefano Zampini     }
302*9566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&D));
303c4fb69feSStefano Zampini   }
304c4fb69feSStefano Zampini 
305300d917bSStefano Zampini   /* Test basis orthogonalization */
30653022affSStefano Zampini   if (testorthog) {
307*9566063dSJacob Faibussowitsch     PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&D));
308*9566063dSJacob Faibussowitsch     PetscCall(MatSetOption(D,MAT_SYMMETRIC,symm));
309*9566063dSJacob Faibussowitsch     PetscCall(MatH2OpusOrthogonalize(D));
310*9566063dSJacob Faibussowitsch     PetscCall(MatMultEqual(B,D,10,&flg));
31153022affSStefano Zampini     if (!flg) {
312*9566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMult error after basis ortogonalization\n"));
31353022affSStefano Zampini     }
314*9566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&D));
31553022affSStefano Zampini   }
31653022affSStefano Zampini 
317300d917bSStefano Zampini   /* Test matrix compression */
31853022affSStefano Zampini   if (testcompress) {
319*9566063dSJacob Faibussowitsch     PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&D));
320*9566063dSJacob Faibussowitsch     PetscCall(MatSetOption(D,MAT_SYMMETRIC,symm));
321*9566063dSJacob Faibussowitsch     PetscCall(MatH2OpusCompress(D,PETSC_SMALL));
322*9566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&D));
32353022affSStefano Zampini   }
32453022affSStefano Zampini 
325300d917bSStefano Zampini   /* Test low-rank update */
326300d917bSStefano Zampini   if (testhlru) {
327300d917bSStefano Zampini     Mat         U, V;
328300d917bSStefano Zampini     PetscScalar *Udata = NULL, *Vdata = NULL;
329300d917bSStefano Zampini 
330300d917bSStefano Zampini     if (ldu) {
331*9566063dSJacob Faibussowitsch       PetscCall(PetscMalloc1(nlr*(n+ldu),&Udata));
332*9566063dSJacob Faibussowitsch       PetscCall(PetscMalloc1(nlr*(n+ldu+2),&Vdata));
333300d917bSStefano Zampini     }
334*9566063dSJacob Faibussowitsch     PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&D));
335*9566063dSJacob Faibussowitsch     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)D),n,PETSC_DECIDE,N,nlr,Udata,&U));
336*9566063dSJacob Faibussowitsch     PetscCall(MatDenseSetLDA(U,n+ldu));
337*9566063dSJacob Faibussowitsch     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)D),n,PETSC_DECIDE,N,nlr,Vdata,&V));
338*9566063dSJacob Faibussowitsch     if (ldu) PetscCall(MatDenseSetLDA(V,n+ldu+2));
339*9566063dSJacob Faibussowitsch     PetscCall(MatSetRandom(U,NULL));
340*9566063dSJacob Faibussowitsch     PetscCall(MatSetRandom(V,NULL));
341*9566063dSJacob Faibussowitsch     PetscCall(MatH2OpusLowRankUpdate(D,U,V,0.5));
342*9566063dSJacob Faibussowitsch     PetscCall(MatH2OpusLowRankUpdate(D,U,V,-0.5));
343*9566063dSJacob Faibussowitsch     PetscCall(MatMultEqual(B,D,10,&flg));
344300d917bSStefano Zampini     if (!flg) {
345*9566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMult error after low-rank update\n"));
346300d917bSStefano Zampini     }
347*9566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&D));
348*9566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&U));
349*9566063dSJacob Faibussowitsch     PetscCall(PetscFree(Udata));
350*9566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&V));
351*9566063dSJacob Faibussowitsch     PetscCall(PetscFree(Vdata));
352300d917bSStefano Zampini   }
353300d917bSStefano Zampini 
354c4fb69feSStefano Zampini   /* check explicit operator */
355c4fb69feSStefano Zampini   if (checkexpl) {
356c4fb69feSStefano Zampini     Mat Be, Bet;
357c4fb69feSStefano Zampini 
358*9566063dSJacob Faibussowitsch     PetscCall(MatComputeOperator(B,MATDENSE,&D));
359*9566063dSJacob Faibussowitsch     PetscCall(MatDuplicate(D,MAT_COPY_VALUES,&Be));
360*9566063dSJacob Faibussowitsch     PetscCall(MatNorm(D,NORM_FROBENIUS,&nB));
361*9566063dSJacob Faibussowitsch     PetscCall(MatViewFromOptions(D,NULL,"-expl_view"));
362*9566063dSJacob Faibussowitsch     PetscCall(MatAXPY(D,-1.0,A,SAME_NONZERO_PATTERN));
363*9566063dSJacob Faibussowitsch     PetscCall(MatViewFromOptions(D,NULL,"-diff_view"));
364*9566063dSJacob Faibussowitsch     PetscCall(MatNorm(D,NORM_FROBENIUS,&nD));
365*9566063dSJacob Faibussowitsch     PetscCall(MatNorm(A,NORM_FROBENIUS,&nA));
366*9566063dSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Approximation error %g (%g / %g, %g)\n",nD/nA,nD,nA,nB));
367*9566063dSJacob Faibussowitsch     PetscCall(MatDestroy(&D));
368c4fb69feSStefano Zampini 
369c4fb69feSStefano Zampini     if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */
370*9566063dSJacob Faibussowitsch       PetscCall(MatTranspose(A,MAT_INPLACE_MATRIX,&A));
371*9566063dSJacob Faibussowitsch       PetscCall(MatComputeOperatorTranspose(B,MATDENSE,&D));
372*9566063dSJacob Faibussowitsch       PetscCall(MatDuplicate(D,MAT_COPY_VALUES,&Bet));
373*9566063dSJacob Faibussowitsch       PetscCall(MatNorm(D,NORM_FROBENIUS,&nB));
374*9566063dSJacob Faibussowitsch       PetscCall(MatViewFromOptions(D,NULL,"-expl_trans_view"));
375*9566063dSJacob Faibussowitsch       PetscCall(MatAXPY(D,-1.0,A,SAME_NONZERO_PATTERN));
376*9566063dSJacob Faibussowitsch       PetscCall(MatViewFromOptions(D,NULL,"-diff_trans_view"));
377*9566063dSJacob Faibussowitsch       PetscCall(MatNorm(D,NORM_FROBENIUS,&nD));
378*9566063dSJacob Faibussowitsch       PetscCall(MatNorm(A,NORM_FROBENIUS,&nA));
379*9566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Approximation error transpose %g (%g / %g, %g)\n",nD/nA,nD,nA,nB));
380*9566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&D));
381c4fb69feSStefano Zampini 
382*9566063dSJacob Faibussowitsch       PetscCall(MatTranspose(Bet,MAT_INPLACE_MATRIX,&Bet));
383*9566063dSJacob Faibussowitsch       PetscCall(MatAXPY(Be,-1.0,Bet,SAME_NONZERO_PATTERN));
384*9566063dSJacob Faibussowitsch       PetscCall(MatViewFromOptions(Be,NULL,"-diff_expl_view"));
385*9566063dSJacob Faibussowitsch       PetscCall(MatNorm(Be,NORM_FROBENIUS,&nB));
386*9566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Approximation error B - (B^T)^T %g\n",nB));
387*9566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&Be));
388*9566063dSJacob Faibussowitsch       PetscCall(MatDestroy(&Bet));
389c4fb69feSStefano Zampini     }
390c4fb69feSStefano Zampini   }
391*9566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&A));
392*9566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&B));
393*9566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&C));
394*9566063dSJacob Faibussowitsch   PetscCall(PetscFree(Cdata));
395*9566063dSJacob Faibussowitsch   PetscCall(PetscFree(Adata));
396*9566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
397b122ec5aSJacob Faibussowitsch   return 0;
398c4fb69feSStefano Zampini }
399c4fb69feSStefano Zampini 
400c4fb69feSStefano Zampini /*TEST
401c4fb69feSStefano Zampini 
402c4fb69feSStefano Zampini    build:
40353022affSStefano Zampini      requires: h2opus
404c4fb69feSStefano Zampini 
405c4fb69feSStefano Zampini #tests from kernel
406c4fb69feSStefano Zampini    test:
40753022affSStefano Zampini      requires: h2opus
408c4fb69feSStefano Zampini      nsize: 1
409c4fb69feSStefano Zampini      suffix: 1
410c4fb69feSStefano Zampini      args: -n {{17 33}} -kernel 1 -dim {{1 2 3}} -symm {{0 1}} -checkexpl -bgpu 0
411c4fb69feSStefano Zampini 
412c4fb69feSStefano Zampini    test:
41353022affSStefano Zampini      requires: h2opus
414c4fb69feSStefano Zampini      nsize: 1
415c4fb69feSStefano Zampini      suffix: 1_ld
416c4fb69feSStefano Zampini      output_file: output/ex66_1.out
417300d917bSStefano Zampini      args: -n 33 -kernel 1 -dim 1 -lda 13 -ldc 11 -ldu 17 -symm 0 -checkexpl -bgpu 0
418c4fb69feSStefano Zampini 
419c4fb69feSStefano Zampini    test:
42053022affSStefano Zampini      requires: h2opus cuda
421c4fb69feSStefano Zampini      nsize: 1
422c4fb69feSStefano Zampini      suffix: 1_cuda
423c4fb69feSStefano Zampini      output_file: output/ex66_1.out
424c4fb69feSStefano Zampini      args: -n {{17 33}} -kernel 1 -dim {{1 2 3}} -symm {{0 1}} -checkexpl -bgpu 1
425c4fb69feSStefano Zampini 
426c4fb69feSStefano Zampini    test:
42753022affSStefano Zampini      requires: h2opus cuda
428c4fb69feSStefano Zampini      nsize: 1
429c4fb69feSStefano Zampini      suffix: 1_cuda_ld
430c4fb69feSStefano Zampini      output_file: output/ex66_1.out
431300d917bSStefano Zampini      args: -n 33 -kernel 1 -dim 1 -lda 13 -ldc 11 -ldu 17 -symm 0 -checkexpl -bgpu 1
432c4fb69feSStefano Zampini 
433c4fb69feSStefano Zampini    test:
43453022affSStefano Zampini      requires: h2opus
43553022affSStefano Zampini      nsize: {{2 3}}
436c4fb69feSStefano Zampini      suffix: 1_par
43753022affSStefano Zampini      args: -n 64 -symm -kernel 1 -dim 1 -ldc 12 -testlayout {{0 1}} -bgpu 0 -cgpu 0
438c4fb69feSStefano Zampini 
439c4fb69feSStefano Zampini    test:
44053022affSStefano Zampini      requires: h2opus cuda
44153022affSStefano Zampini      nsize: {{2 3}}
442c4fb69feSStefano Zampini      suffix: 1_par_cuda
44353022affSStefano Zampini      args: -n 64 -symm -kernel 1 -dim 1 -ldc 12 -testlayout {{0 1}} -bgpu {{0 1}} -cgpu {{0 1}}
444c4fb69feSStefano Zampini      output_file: output/ex66_1_par.out
445c4fb69feSStefano Zampini 
446c4fb69feSStefano Zampini #tests from matrix sampling (parallel or unsymmetric not supported)
447c4fb69feSStefano Zampini    test:
44853022affSStefano Zampini      requires: h2opus
449c4fb69feSStefano Zampini      nsize: 1
450c4fb69feSStefano Zampini      suffix: 2
451c4fb69feSStefano Zampini      args: -n {{17 33}} -kernel 0 -dim 2 -symm 1 -checkexpl -bgpu 0
452c4fb69feSStefano Zampini 
453c4fb69feSStefano Zampini    test:
45453022affSStefano Zampini      requires: h2opus cuda
455c4fb69feSStefano Zampini      nsize: 1
456c4fb69feSStefano Zampini      suffix: 2_cuda
457c4fb69feSStefano Zampini      output_file: output/ex66_2.out
458c4fb69feSStefano Zampini      args: -n {{17 33}} -kernel 0 -dim 2 -symm 1 -checkexpl -bgpu {{0 1}} -agpu {{0 1}}
459c4fb69feSStefano Zampini 
46053022affSStefano Zampini #tests view operation
46153022affSStefano Zampini    test:
46253022affSStefano Zampini      requires: h2opus !cuda
46353022affSStefano Zampini      filter: grep -v "MPI processes" | grep -v "\[" | grep -v "\]"
46453022affSStefano Zampini      nsize: {{1 2 3}}
46553022affSStefano Zampini      suffix: view
46653022affSStefano 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
46753022affSStefano Zampini 
46853022affSStefano Zampini    test:
46953022affSStefano Zampini      requires: h2opus cuda
47053022affSStefano Zampini      filter: grep -v "MPI processes" | grep -v "\[" | grep -v "\]"
47153022affSStefano Zampini      nsize: {{1 2 3}}
47253022affSStefano Zampini      suffix: view_cuda
47353022affSStefano 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
47453022affSStefano Zampini 
475c4fb69feSStefano Zampini TEST*/
476