xref: /petsc/src/mat/tests/ex66.c (revision d5b43468fb8780a8feea140ccd6fa3e6a50411cc)
1 static char help[] = "Tests MATH2OPUS\n\n";
2 
3 #include <petscmat.h>
4 #include <petscsf.h>
5 
6 static PetscScalar GenEntry_Symm(PetscInt sdim, PetscReal x[], PetscReal y[], void *ctx)
7 {
8   PetscInt  d;
9   PetscReal clength = sdim == 3 ? 0.2 : 0.1;
10   PetscReal dist, diff = 0.0;
11 
12   for (d = 0; d < sdim; d++) diff += (x[d] - y[d]) * (x[d] - y[d]);
13   dist = PetscSqrtReal(diff);
14   return PetscExpReal(-dist / clength);
15 }
16 
17 static PetscScalar GenEntry_Unsymm(PetscInt sdim, PetscReal x[], PetscReal y[], void *ctx)
18 {
19   PetscInt  d;
20   PetscReal clength = sdim == 3 ? 0.2 : 0.1;
21   PetscReal dist, diff = 0.0, nx = 0.0, ny = 0.0;
22 
23   for (d = 0; d < sdim; d++) nx += x[d] * x[d];
24   for (d = 0; d < sdim; d++) ny += y[d] * y[d];
25   for (d = 0; d < sdim; d++) diff += (x[d] - y[d]) * (x[d] - y[d]);
26   dist = PetscSqrtReal(diff);
27   return nx > ny ? PetscExpReal(-dist / clength) : PetscExpReal(-dist / clength) + 1.;
28 }
29 
30 int main(int argc, char **argv)
31 {
32   Mat          A, B, C, D;
33   Vec          v, x, y, Ax, Ay, Bx, By;
34   PetscRandom  r;
35   PetscLayout  map;
36   PetscScalar *Adata = NULL, *Cdata = NULL, scale = 1.0;
37   PetscReal   *coords, nA, nD, nB, err, nX, norms[3];
38   PetscInt     N, n = 64, dim = 1, i, j, nrhs = 11, lda = 0, ldc = 0, ldu = 0, nlr = 7, nt, ntrials = 2;
39   PetscMPIInt  size, rank;
40   PetscBool    testlayout = PETSC_FALSE, flg, symm = PETSC_FALSE, Asymm = PETSC_TRUE, kernel = PETSC_TRUE;
41   PetscBool    checkexpl = PETSC_FALSE, agpu = PETSC_FALSE, bgpu = PETSC_FALSE, cgpu = PETSC_FALSE, flgglob;
42   PetscBool    testtrans, testnorm, randommat = PETSC_TRUE, testorthog, testcompress, testhlru;
43   void (*approxnormfunc)(void);
44   void (*Anormfunc)(void);
45 
46 #if defined(PETSC_HAVE_MPI_INIT_THREAD)
47   PETSC_MPI_THREAD_REQUIRED = MPI_THREAD_MULTIPLE;
48 #endif
49   PetscFunctionBeginUser;
50   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
51   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ng", &N, &flgglob));
52   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
53   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
54   PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL));
55   PetscCall(PetscOptionsGetInt(NULL, NULL, "-lda", &lda, NULL));
56   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ldc", &ldc, NULL));
57   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nlr", &nlr, NULL));
58   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ldu", &ldu, NULL));
59   PetscCall(PetscOptionsGetInt(NULL, NULL, "-matmattrials", &ntrials, NULL));
60   PetscCall(PetscOptionsGetBool(NULL, NULL, "-randommat", &randommat, NULL));
61   if (!flgglob) PetscCall(PetscOptionsGetBool(NULL, NULL, "-testlayout", &testlayout, NULL));
62   PetscCall(PetscOptionsGetBool(NULL, NULL, "-Asymm", &Asymm, NULL));
63   PetscCall(PetscOptionsGetBool(NULL, NULL, "-symm", &symm, NULL));
64   PetscCall(PetscOptionsGetBool(NULL, NULL, "-kernel", &kernel, NULL));
65   PetscCall(PetscOptionsGetBool(NULL, NULL, "-checkexpl", &checkexpl, NULL));
66   PetscCall(PetscOptionsGetBool(NULL, NULL, "-agpu", &agpu, NULL));
67   PetscCall(PetscOptionsGetBool(NULL, NULL, "-bgpu", &bgpu, NULL));
68   PetscCall(PetscOptionsGetBool(NULL, NULL, "-cgpu", &cgpu, NULL));
69   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-scale", &scale, NULL));
70   if (!Asymm) symm = PETSC_FALSE;
71 
72   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
73 
74   /* Disable tests for unimplemented variants */
75   testtrans    = (PetscBool)(size == 1 || symm);
76   testnorm     = (PetscBool)(size == 1 || symm);
77   testorthog   = (PetscBool)(size == 1 || symm);
78   testcompress = (PetscBool)(size == 1 || symm);
79   testhlru     = (PetscBool)(size == 1);
80 
81   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
82   PetscCall(PetscLayoutCreate(PETSC_COMM_WORLD, &map));
83   if (testlayout) {
84     if (rank % 2) n = PetscMax(2 * n - 5 * rank, 0);
85     else n = 2 * n + rank;
86   }
87   if (!flgglob) {
88     PetscCall(PetscLayoutSetLocalSize(map, n));
89     PetscCall(PetscLayoutSetUp(map));
90     PetscCall(PetscLayoutGetSize(map, &N));
91   } else {
92     PetscCall(PetscLayoutSetSize(map, N));
93     PetscCall(PetscLayoutSetUp(map));
94     PetscCall(PetscLayoutGetLocalSize(map, &n));
95   }
96   PetscCall(PetscLayoutDestroy(&map));
97 
98   if (lda) PetscCall(PetscMalloc1(N * (n + lda), &Adata));
99   PetscCall(MatCreateDense(PETSC_COMM_WORLD, n, n, N, N, Adata, &A));
100   PetscCall(MatDenseSetLDA(A, n + lda));
101 
102   /* Create random points; these are replicated in order to populate a dense matrix and to compare sequential and dense runs
103      The constructor for MATH2OPUS can take as input the distributed coordinates and replicates them internally in case
104      a kernel construction is requested */
105   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &r));
106   PetscCall(PetscRandomSetFromOptions(r));
107   PetscCall(PetscRandomSetSeed(r, 123456));
108   PetscCall(PetscRandomSeed(r));
109   PetscCall(PetscMalloc1(N * dim, &coords));
110   PetscCall(PetscRandomGetValuesReal(r, N * dim, coords));
111   PetscCall(PetscRandomDestroy(&r));
112 
113   if (kernel || !randommat) {
114     MatH2OpusKernel k = Asymm ? GenEntry_Symm : GenEntry_Unsymm;
115     PetscInt        ist, ien;
116 
117     PetscCall(MatGetOwnershipRange(A, &ist, &ien));
118     for (i = ist; i < ien; i++) {
119       for (j = 0; j < N; j++) PetscCall(MatSetValue(A, i, j, (*k)(dim, coords + i * dim, coords + j * dim, NULL), INSERT_VALUES));
120     }
121     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
122     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
123     if (kernel) {
124       PetscCall(MatCreateH2OpusFromKernel(PETSC_COMM_WORLD, n, n, N, N, dim, coords + ist * dim, PETSC_TRUE, k, NULL, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, &B));
125     } else {
126       PetscCall(MatCreateH2OpusFromMat(A, dim, coords + ist * dim, PETSC_TRUE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, &B));
127     }
128   } else {
129     PetscInt ist;
130 
131     PetscCall(MatGetOwnershipRange(A, &ist, NULL));
132     PetscCall(MatSetRandom(A, NULL));
133     if (Asymm) {
134       PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B));
135       PetscCall(MatAXPY(A, 1.0, B, SAME_NONZERO_PATTERN));
136       PetscCall(MatDestroy(&B));
137       PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
138     }
139     PetscCall(MatCreateH2OpusFromMat(A, dim, coords + ist * dim, PETSC_TRUE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, &B));
140   }
141   PetscCall(PetscFree(coords));
142   if (agpu) PetscCall(MatConvert(A, MATDENSECUDA, MAT_INPLACE_MATRIX, &A));
143   PetscCall(MatViewFromOptions(A, NULL, "-A_view"));
144 
145   PetscCall(MatSetOption(B, MAT_SYMMETRIC, symm));
146 
147   /* assemble the H-matrix */
148   PetscCall(MatBindToCPU(B, (PetscBool)!bgpu));
149   PetscCall(MatSetFromOptions(B));
150   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
151   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
152   PetscCall(MatViewFromOptions(B, NULL, "-B_view"));
153 
154   /* Test MatScale */
155   PetscCall(MatScale(A, scale));
156   PetscCall(MatScale(B, scale));
157 
158   /* Test MatMult */
159   PetscCall(MatCreateVecs(A, &Ax, &Ay));
160   PetscCall(MatCreateVecs(B, &Bx, &By));
161   PetscCall(VecSetRandom(Ax, NULL));
162   PetscCall(VecCopy(Ax, Bx));
163   PetscCall(MatMult(A, Ax, Ay));
164   PetscCall(MatMult(B, Bx, By));
165   PetscCall(VecViewFromOptions(Ay, NULL, "-mult_vec_view"));
166   PetscCall(VecViewFromOptions(By, NULL, "-mult_vec_view"));
167   PetscCall(VecNorm(Ay, NORM_INFINITY, &nX));
168   PetscCall(VecAXPY(Ay, -1.0, By));
169   PetscCall(VecViewFromOptions(Ay, NULL, "-mult_vec_view"));
170   PetscCall(VecNorm(Ay, NORM_INFINITY, &err));
171   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMult err %g\n", err / nX));
172   PetscCall(VecScale(By, -1.0));
173   PetscCall(MatMultAdd(B, Bx, By, By));
174   PetscCall(VecNorm(By, NORM_INFINITY, &err));
175   PetscCall(VecViewFromOptions(By, NULL, "-mult_vec_view"));
176   if (err > 10. * PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMultAdd err %g\n", err));
177 
178   /* Test MatNorm */
179   PetscCall(MatNorm(A, NORM_INFINITY, &norms[0]));
180   PetscCall(MatNorm(A, NORM_1, &norms[1]));
181   norms[2] = -1.; /* NORM_2 not supported */
182   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]));
183   PetscCall(MatGetOperation(A, MATOP_NORM, &Anormfunc));
184   PetscCall(MatGetOperation(B, MATOP_NORM, &approxnormfunc));
185   PetscCall(MatSetOperation(A, MATOP_NORM, approxnormfunc));
186   PetscCall(MatNorm(A, NORM_INFINITY, &norms[0]));
187   PetscCall(MatNorm(A, NORM_1, &norms[1]));
188   PetscCall(MatNorm(A, NORM_2, &norms[2]));
189   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]));
190   if (testnorm) {
191     PetscCall(MatNorm(B, NORM_INFINITY, &norms[0]));
192     PetscCall(MatNorm(B, NORM_1, &norms[1]));
193     PetscCall(MatNorm(B, NORM_2, &norms[2]));
194   } else {
195     norms[0] = -1.;
196     norms[1] = -1.;
197     norms[2] = -1.;
198   }
199   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]));
200   PetscCall(MatSetOperation(A, MATOP_NORM, Anormfunc));
201 
202   /* Test MatDuplicate */
203   PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &D));
204   PetscCall(MatSetOption(D, MAT_SYMMETRIC, symm));
205   PetscCall(MatMultEqual(B, D, 10, &flg));
206   if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMult error after MatDuplicate\n"));
207   if (testtrans) {
208     PetscCall(MatMultTransposeEqual(B, D, 10, &flg));
209     if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMultTranspose error after MatDuplicate\n"));
210   }
211   PetscCall(MatDestroy(&D));
212 
213   if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */
214     PetscCall(VecSetRandom(Ay, NULL));
215     PetscCall(VecCopy(Ay, By));
216     PetscCall(MatMultTranspose(A, Ay, Ax));
217     PetscCall(MatMultTranspose(B, By, Bx));
218     PetscCall(VecViewFromOptions(Ax, NULL, "-multtrans_vec_view"));
219     PetscCall(VecViewFromOptions(Bx, NULL, "-multtrans_vec_view"));
220     PetscCall(VecNorm(Ax, NORM_INFINITY, &nX));
221     PetscCall(VecAXPY(Ax, -1.0, Bx));
222     PetscCall(VecViewFromOptions(Ax, NULL, "-multtrans_vec_view"));
223     PetscCall(VecNorm(Ax, NORM_INFINITY, &err));
224     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMultTranspose err %g\n", err / nX));
225     PetscCall(VecScale(Bx, -1.0));
226     PetscCall(MatMultTransposeAdd(B, By, Bx, Bx));
227     PetscCall(VecNorm(Bx, NORM_INFINITY, &err));
228     if (err > 10. * PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMultTransposeAdd err %g\n", err));
229   }
230   PetscCall(VecDestroy(&Ax));
231   PetscCall(VecDestroy(&Ay));
232   PetscCall(VecDestroy(&Bx));
233   PetscCall(VecDestroy(&By));
234 
235   /* Test MatMatMult */
236   if (ldc) PetscCall(PetscMalloc1(nrhs * (n + ldc), &Cdata));
237   PetscCall(MatCreateDense(PETSC_COMM_WORLD, n, PETSC_DECIDE, N, nrhs, Cdata, &C));
238   PetscCall(MatDenseSetLDA(C, n + ldc));
239   PetscCall(MatSetRandom(C, NULL));
240   if (cgpu) PetscCall(MatConvert(C, MATDENSECUDA, MAT_INPLACE_MATRIX, &C));
241   for (nt = 0; nt < ntrials; nt++) {
242     PetscCall(MatMatMult(B, C, nt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, PETSC_DEFAULT, &D));
243     PetscCall(MatViewFromOptions(D, NULL, "-bc_view"));
244     PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)D, &flg, MATSEQDENSE, MATMPIDENSE, ""));
245     if (flg) {
246       PetscCall(MatCreateVecs(B, &x, &y));
247       PetscCall(MatCreateVecs(D, NULL, &v));
248       for (i = 0; i < nrhs; i++) {
249         PetscCall(MatGetColumnVector(D, v, i));
250         PetscCall(MatGetColumnVector(C, x, i));
251         PetscCall(MatMult(B, x, y));
252         PetscCall(VecAXPY(y, -1.0, v));
253         PetscCall(VecNorm(y, NORM_INFINITY, &err));
254         if (err > 10. * PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMat err %" PetscInt_FMT " %g\n", i, err));
255       }
256       PetscCall(VecDestroy(&y));
257       PetscCall(VecDestroy(&x));
258       PetscCall(VecDestroy(&v));
259     }
260   }
261   PetscCall(MatDestroy(&D));
262 
263   /* Test MatTransposeMatMult */
264   if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */
265     for (nt = 0; nt < ntrials; nt++) {
266       PetscCall(MatTransposeMatMult(B, C, nt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, PETSC_DEFAULT, &D));
267       PetscCall(MatViewFromOptions(D, NULL, "-btc_view"));
268       PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)D, &flg, MATSEQDENSE, MATMPIDENSE, ""));
269       if (flg) {
270         PetscCall(MatCreateVecs(B, &y, &x));
271         PetscCall(MatCreateVecs(D, NULL, &v));
272         for (i = 0; i < nrhs; i++) {
273           PetscCall(MatGetColumnVector(D, v, i));
274           PetscCall(MatGetColumnVector(C, x, i));
275           PetscCall(MatMultTranspose(B, x, y));
276           PetscCall(VecAXPY(y, -1.0, v));
277           PetscCall(VecNorm(y, NORM_INFINITY, &err));
278           if (err > 10. * PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatTransMat err %" PetscInt_FMT " %g\n", i, err));
279         }
280         PetscCall(VecDestroy(&y));
281         PetscCall(VecDestroy(&x));
282         PetscCall(VecDestroy(&v));
283       }
284     }
285     PetscCall(MatDestroy(&D));
286   }
287 
288   /* Test basis orthogonalization */
289   if (testorthog) {
290     PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &D));
291     PetscCall(MatSetOption(D, MAT_SYMMETRIC, symm));
292     PetscCall(MatH2OpusOrthogonalize(D));
293     PetscCall(MatMultEqual(B, D, 10, &flg));
294     if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMult error after basis ortogonalization\n"));
295     PetscCall(MatDestroy(&D));
296   }
297 
298   /* Test matrix compression */
299   if (testcompress) {
300     PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &D));
301     PetscCall(MatSetOption(D, MAT_SYMMETRIC, symm));
302     PetscCall(MatH2OpusCompress(D, PETSC_SMALL));
303     PetscCall(MatDestroy(&D));
304   }
305 
306   /* Test low-rank update */
307   if (testhlru) {
308     Mat          U, V;
309     PetscScalar *Udata = NULL, *Vdata = NULL;
310 
311     if (ldu) {
312       PetscCall(PetscMalloc1(nlr * (n + ldu), &Udata));
313       PetscCall(PetscMalloc1(nlr * (n + ldu + 2), &Vdata));
314     }
315     PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &D));
316     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)D), n, PETSC_DECIDE, N, nlr, Udata, &U));
317     PetscCall(MatDenseSetLDA(U, n + ldu));
318     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)D), n, PETSC_DECIDE, N, nlr, Vdata, &V));
319     if (ldu) PetscCall(MatDenseSetLDA(V, n + ldu + 2));
320     PetscCall(MatSetRandom(U, NULL));
321     PetscCall(MatSetRandom(V, NULL));
322     PetscCall(MatH2OpusLowRankUpdate(D, U, V, 0.5));
323     PetscCall(MatH2OpusLowRankUpdate(D, U, V, -0.5));
324     PetscCall(MatMultEqual(B, D, 10, &flg));
325     if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMult error after low-rank update\n"));
326     PetscCall(MatDestroy(&D));
327     PetscCall(MatDestroy(&U));
328     PetscCall(PetscFree(Udata));
329     PetscCall(MatDestroy(&V));
330     PetscCall(PetscFree(Vdata));
331   }
332 
333   /* check explicit operator */
334   if (checkexpl) {
335     Mat Be, Bet;
336 
337     PetscCall(MatComputeOperator(B, MATDENSE, &D));
338     PetscCall(MatDuplicate(D, MAT_COPY_VALUES, &Be));
339     PetscCall(MatNorm(D, NORM_FROBENIUS, &nB));
340     PetscCall(MatViewFromOptions(D, NULL, "-expl_view"));
341     PetscCall(MatAXPY(D, -1.0, A, SAME_NONZERO_PATTERN));
342     PetscCall(MatViewFromOptions(D, NULL, "-diff_view"));
343     PetscCall(MatNorm(D, NORM_FROBENIUS, &nD));
344     PetscCall(MatNorm(A, NORM_FROBENIUS, &nA));
345     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Approximation error %g (%g / %g, %g)\n", nD / nA, nD, nA, nB));
346     PetscCall(MatDestroy(&D));
347 
348     if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */
349       PetscCall(MatTranspose(A, MAT_INPLACE_MATRIX, &A));
350       PetscCall(MatComputeOperatorTranspose(B, MATDENSE, &D));
351       PetscCall(MatDuplicate(D, MAT_COPY_VALUES, &Bet));
352       PetscCall(MatNorm(D, NORM_FROBENIUS, &nB));
353       PetscCall(MatViewFromOptions(D, NULL, "-expl_trans_view"));
354       PetscCall(MatAXPY(D, -1.0, A, SAME_NONZERO_PATTERN));
355       PetscCall(MatViewFromOptions(D, NULL, "-diff_trans_view"));
356       PetscCall(MatNorm(D, NORM_FROBENIUS, &nD));
357       PetscCall(MatNorm(A, NORM_FROBENIUS, &nA));
358       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Approximation error transpose %g (%g / %g, %g)\n", nD / nA, nD, nA, nB));
359       PetscCall(MatDestroy(&D));
360 
361       PetscCall(MatTranspose(Bet, MAT_INPLACE_MATRIX, &Bet));
362       PetscCall(MatAXPY(Be, -1.0, Bet, SAME_NONZERO_PATTERN));
363       PetscCall(MatViewFromOptions(Be, NULL, "-diff_expl_view"));
364       PetscCall(MatNorm(Be, NORM_FROBENIUS, &nB));
365       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Approximation error B - (B^T)^T %g\n", nB));
366       PetscCall(MatDestroy(&Be));
367       PetscCall(MatDestroy(&Bet));
368     }
369   }
370   PetscCall(MatDestroy(&A));
371   PetscCall(MatDestroy(&B));
372   PetscCall(MatDestroy(&C));
373   PetscCall(PetscFree(Cdata));
374   PetscCall(PetscFree(Adata));
375   PetscCall(PetscFinalize());
376   return 0;
377 }
378 
379 /*TEST
380 
381    build:
382      requires: h2opus
383 
384 #tests from kernel
385    test:
386      requires: h2opus
387      nsize: 1
388      suffix: 1
389      args: -n {{17 33}} -kernel 1 -dim {{1 2 3}} -symm {{0 1}} -checkexpl -bgpu 0
390 
391    test:
392      requires: h2opus
393      nsize: 1
394      suffix: 1_ld
395      output_file: output/ex66_1.out
396      args: -n 33 -kernel 1 -dim 1 -lda 13 -ldc 11 -ldu 17 -symm 0 -checkexpl -bgpu 0
397 
398    test:
399      requires: h2opus cuda
400      nsize: 1
401      suffix: 1_cuda
402      output_file: output/ex66_1.out
403      args: -n {{17 33}} -kernel 1 -dim {{1 2 3}} -symm {{0 1}} -checkexpl -bgpu 1
404 
405    test:
406      requires: h2opus cuda
407      nsize: 1
408      suffix: 1_cuda_ld
409      output_file: output/ex66_1.out
410      args: -n 33 -kernel 1 -dim 1 -lda 13 -ldc 11 -ldu 17 -symm 0 -checkexpl -bgpu 1
411 
412    test:
413      requires: h2opus
414      nsize: {{2 3}}
415      suffix: 1_par
416      args: -n 64 -symm -kernel 1 -dim 1 -ldc 12 -testlayout {{0 1}} -bgpu 0 -cgpu 0
417 
418    test:
419      requires: h2opus cuda
420      nsize: {{2 3}}
421      suffix: 1_par_cuda
422      args: -n 64 -symm -kernel 1 -dim 1 -ldc 12 -testlayout {{0 1}} -bgpu {{0 1}} -cgpu {{0 1}}
423      output_file: output/ex66_1_par.out
424 
425 #tests from matrix sampling (parallel or unsymmetric not supported)
426    test:
427      requires: h2opus
428      nsize: 1
429      suffix: 2
430      args: -n {{17 33}} -kernel 0 -dim 2 -symm 1 -checkexpl -bgpu 0
431 
432    test:
433      requires: h2opus cuda
434      nsize: 1
435      suffix: 2_cuda
436      output_file: output/ex66_2.out
437      args: -n {{17 29}} -kernel 0 -dim 2 -symm 1 -checkexpl -bgpu {{0 1}} -agpu {{0 1}}
438 
439 #tests view operation
440    test:
441      requires: h2opus !cuda
442      filter: grep -v " MPI process" | grep -v "\[" | grep -v "\]"
443      nsize: {{1 2 3}}
444      suffix: view
445      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
446 
447    test:
448      requires: h2opus cuda
449      filter: grep -v " MPI process" | grep -v "\[" | grep -v "\]"
450      nsize: {{1 2 3}}
451      suffix: view_cuda
452      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
453 
454 TEST*/
455