Actual source code: ex66.c

  1: static char help[] = "Tests MATH2OPUS\n\n";

  3: #include <petscmat.h>
  4: #include <petscsf.h>

  6: static PetscScalar GenEntry_Symm(PetscInt sdim, PetscReal x[], PetscReal y[], PetscCtx ctx)
  7: {
  8:   PetscReal clength = sdim == 3 ? 0.2 : 0.1;
  9:   PetscReal dist, diff = 0.0;

 11:   for (PetscInt d = 0; d < sdim; d++) diff += (x[d] - y[d]) * (x[d] - y[d]);
 12:   dist = PetscSqrtReal(diff);
 13:   return PetscExpReal(-dist / clength);
 14: }

 16: static PetscScalar GenEntry_Unsymm(PetscInt sdim, PetscReal x[], PetscReal y[], PetscCtx ctx)
 17: {
 18:   PetscReal clength = sdim == 3 ? 0.2 : 0.1;
 19:   PetscReal dist, diff = 0.0, nx = 0.0, ny = 0.0;

 21:   for (PetscInt d = 0; d < sdim; d++) nx += x[d] * x[d];
 22:   for (PetscInt d = 0; d < sdim; d++) ny += y[d] * y[d];
 23:   for (PetscInt d = 0; d < sdim; d++) diff += (x[d] - y[d]) * (x[d] - y[d]);
 24:   dist = PetscSqrtReal(diff);
 25:   return nx > ny ? PetscExpReal(-dist / clength) : PetscExpReal(-dist / clength) + 1.;
 26: }

 28: int main(int argc, char **argv)
 29: {
 30:   Mat               A, B, C, D;
 31:   Vec               v, x, y, Ax, Ay, Bx, By;
 32:   PetscRandom       r;
 33:   PetscLayout       map;
 34:   PetscScalar      *Adata = NULL, *Cdata = NULL, scale = 1.0;
 35:   PetscReal        *coords, nA, nD, nB, err, nX, norms[3];
 36:   PetscInt          N, n = 64, dim = 1, i, j, nrhs = 11, lda = 0, ldc = 0, ldu = 0, nlr = 7, nt, ntrials = 2;
 37:   PetscMPIInt       size, rank;
 38:   PetscBool         testlayout = PETSC_FALSE, flg, symm = PETSC_FALSE, Asymm = PETSC_TRUE, kernel = PETSC_TRUE;
 39:   PetscBool         checkexpl = PETSC_FALSE, agpu = PETSC_FALSE, bgpu = PETSC_FALSE, cgpu = PETSC_FALSE, flgglob;
 40:   PetscBool         testtrans, testnorm, randommat = PETSC_TRUE, testorthog, testcompress, testhlru;
 41:   PetscErrorCodeFn *approxnormfunc;
 42:   PetscErrorCodeFn *Anormfunc;

 44: #if PetscDefined(HAVE_MPI_INIT_THREAD)
 45:   PETSC_MPI_THREAD_REQUIRED = MPI_THREAD_MULTIPLE;
 46: #endif
 47:   PetscFunctionBeginUser;
 48:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 49:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ng", &N, &flgglob));
 50:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 51:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
 52:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL));
 53:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-lda", &lda, NULL));
 54:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ldc", &ldc, NULL));
 55:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nlr", &nlr, NULL));
 56:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ldu", &ldu, NULL));
 57:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-matmattrials", &ntrials, NULL));
 58:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-randommat", &randommat, NULL));
 59:   if (!flgglob) PetscCall(PetscOptionsGetBool(NULL, NULL, "-testlayout", &testlayout, NULL));
 60:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-Asymm", &Asymm, NULL));
 61:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-symm", &symm, NULL));
 62:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-kernel", &kernel, NULL));
 63:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-checkexpl", &checkexpl, NULL));
 64:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-agpu", &agpu, NULL));
 65:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-bgpu", &bgpu, NULL));
 66:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-cgpu", &cgpu, NULL));
 67:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-scale", &scale, NULL));
 68:   if (!Asymm) symm = PETSC_FALSE;

 70:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));

 72:   /* Disable tests for unimplemented variants */
 73:   testtrans    = (PetscBool)(size == 1 || symm);
 74:   testnorm     = (PetscBool)(size == 1 || symm);
 75:   testorthog   = (PetscBool)(size == 1 || symm);
 76:   testcompress = (PetscBool)(size == 1 || symm);
 77:   testhlru     = (PetscBool)(size == 1);

 79:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 80:   PetscCall(PetscLayoutCreate(PETSC_COMM_WORLD, &map));
 81:   if (testlayout) {
 82:     if (rank % 2) n = PetscMax(2 * n - 5 * rank, 0);
 83:     else n = 2 * n + rank;
 84:   }
 85:   if (!flgglob) {
 86:     PetscCall(PetscLayoutSetLocalSize(map, n));
 87:     PetscCall(PetscLayoutSetUp(map));
 88:     PetscCall(PetscLayoutGetSize(map, &N));
 89:   } else {
 90:     PetscCall(PetscLayoutSetSize(map, N));
 91:     PetscCall(PetscLayoutSetUp(map));
 92:     PetscCall(PetscLayoutGetLocalSize(map, &n));
 93:   }
 94:   PetscCall(PetscLayoutDestroy(&map));

 96:   if (lda) PetscCall(PetscMalloc1(N * (n + lda), &Adata));
 97:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, n, n, N, N, Adata, &A));
 98:   PetscCall(MatDenseSetLDA(A, n + lda));

100:   /* Create random points; these are replicated in order to populate a dense matrix and to compare sequential and dense runs
101:      The constructor for MATH2OPUS can take as input the distributed coordinates and replicates them internally in case
102:      a kernel construction is requested */
103:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &r));
104:   PetscCall(PetscRandomSetFromOptions(r));
105:   PetscCall(PetscRandomSetSeed(r, 123456));
106:   PetscCall(PetscRandomSeed(r));
107:   PetscCall(PetscMalloc1(N * dim, &coords));
108:   PetscCall(PetscRandomGetValuesReal(r, N * dim, coords));
109:   PetscCall(PetscRandomDestroy(&r));

111:   if (kernel || !randommat) {
112:     MatH2OpusKernelFn *k = Asymm ? GenEntry_Symm : GenEntry_Unsymm;
113:     PetscInt           ist, ien;

115:     PetscCall(MatGetOwnershipRange(A, &ist, &ien));
116:     for (i = ist; i < ien; i++) {
117:       for (j = 0; j < N; j++) PetscCall(MatSetValue(A, i, j, (*k)(dim, coords + i * dim, coords + j * dim, NULL), INSERT_VALUES));
118:     }
119:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
120:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
121:     if (kernel) {
122:       PetscCall(MatCreateH2OpusFromKernel(PETSC_COMM_WORLD, n, n, N, N, dim, coords + ist * dim, PETSC_TRUE, k, NULL, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, &B));
123:     } else {
124:       PetscCall(MatCreateH2OpusFromMat(A, dim, coords + ist * dim, PETSC_TRUE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, &B));
125:     }
126:   } else {
127:     PetscInt ist;

129:     PetscCall(MatGetOwnershipRange(A, &ist, NULL));
130:     PetscCall(MatSetRandom(A, NULL));
131:     if (Asymm) {
132:       PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B));
133:       PetscCall(MatAXPY(A, 1.0, B, SAME_NONZERO_PATTERN));
134:       PetscCall(MatDestroy(&B));
135:       PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
136:     }
137:     PetscCall(MatCreateH2OpusFromMat(A, dim, coords + ist * dim, PETSC_TRUE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, &B));
138:   }
139:   PetscCall(PetscFree(coords));
140:   if (agpu) PetscCall(MatConvert(A, MATDENSECUDA, MAT_INPLACE_MATRIX, &A));
141:   PetscCall(MatViewFromOptions(A, NULL, "-A_view"));

143:   PetscCall(MatSetOption(B, MAT_SYMMETRIC, symm));

145:   /* assemble the H-matrix */
146:   PetscCall(MatBindToCPU(B, (PetscBool)!bgpu));
147:   PetscCall(MatSetFromOptions(B));
148:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
149:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
150:   PetscCall(MatViewFromOptions(B, NULL, "-B_view"));

152:   /* Test MatScale */
153:   PetscCall(MatScale(A, scale));
154:   PetscCall(MatScale(B, scale));

156:   /* Test MatMult */
157:   PetscCall(MatCreateVecs(A, &Ax, &Ay));
158:   PetscCall(MatCreateVecs(B, &Bx, &By));
159:   PetscCall(VecSetRandom(Ax, NULL));
160:   PetscCall(VecCopy(Ax, Bx));
161:   PetscCall(MatMult(A, Ax, Ay));
162:   PetscCall(MatMult(B, Bx, By));
163:   PetscCall(VecViewFromOptions(Ay, NULL, "-mult_vec_view"));
164:   PetscCall(VecViewFromOptions(By, NULL, "-mult_vec_view"));
165:   PetscCall(VecNorm(Ay, NORM_INFINITY, &nX));
166:   PetscCall(VecAXPY(Ay, -1.0, By));
167:   PetscCall(VecViewFromOptions(Ay, NULL, "-mult_vec_view"));
168:   PetscCall(VecNorm(Ay, NORM_INFINITY, &err));
169:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMult err %g\n", err / nX));
170:   PetscCall(VecScale(By, -1.0));
171:   PetscCall(MatMultAdd(B, Bx, By, By));
172:   PetscCall(VecNorm(By, NORM_INFINITY, &err));
173:   PetscCall(VecViewFromOptions(By, NULL, "-mult_vec_view"));
174:   if (err > 10. * PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMultAdd err %g\n", err));

176:   /* Test MatNorm */
177:   PetscCall(MatNorm(A, NORM_INFINITY, &norms[0]));
178:   PetscCall(MatNorm(A, NORM_1, &norms[1]));
179:   norms[2] = -1.; /* NORM_2 not supported */
180:   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]));
181:   PetscCall(MatGetOperation(A, MATOP_NORM, &Anormfunc));
182:   PetscCall(MatGetOperation(B, MATOP_NORM, &approxnormfunc));
183:   PetscCall(MatSetOperation(A, MATOP_NORM, approxnormfunc));
184:   PetscCall(MatNorm(A, NORM_INFINITY, &norms[0]));
185:   PetscCall(MatNorm(A, NORM_1, &norms[1]));
186:   PetscCall(MatNorm(A, NORM_2, &norms[2]));
187:   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]));
188:   if (testnorm) {
189:     PetscCall(MatNorm(B, NORM_INFINITY, &norms[0]));
190:     PetscCall(MatNorm(B, NORM_1, &norms[1]));
191:     PetscCall(MatNorm(B, NORM_2, &norms[2]));
192:   } else {
193:     norms[0] = -1.;
194:     norms[1] = -1.;
195:     norms[2] = -1.;
196:   }
197:   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]));
198:   PetscCall(MatSetOperation(A, MATOP_NORM, Anormfunc));

200:   /* Test MatDuplicate */
201:   PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &D));
202:   PetscCall(MatSetOption(D, MAT_SYMMETRIC, symm));
203:   PetscCall(MatMultEqual(B, D, 10, &flg));
204:   if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMult error after MatDuplicate\n"));
205:   if (testtrans) {
206:     PetscCall(MatMultTransposeEqual(B, D, 10, &flg));
207:     if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMultTranspose error after MatDuplicate\n"));
208:   }
209:   PetscCall(MatDestroy(&D));

211:   if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */
212:     PetscCall(VecSetRandom(Ay, NULL));
213:     PetscCall(VecCopy(Ay, By));
214:     PetscCall(MatMultTranspose(A, Ay, Ax));
215:     PetscCall(MatMultTranspose(B, By, Bx));
216:     PetscCall(VecViewFromOptions(Ax, NULL, "-multtrans_vec_view"));
217:     PetscCall(VecViewFromOptions(Bx, NULL, "-multtrans_vec_view"));
218:     PetscCall(VecNorm(Ax, NORM_INFINITY, &nX));
219:     PetscCall(VecAXPY(Ax, -1.0, Bx));
220:     PetscCall(VecViewFromOptions(Ax, NULL, "-multtrans_vec_view"));
221:     PetscCall(VecNorm(Ax, NORM_INFINITY, &err));
222:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMultTranspose err %g\n", err / nX));
223:     PetscCall(VecScale(Bx, -1.0));
224:     PetscCall(MatMultTransposeAdd(B, By, Bx, Bx));
225:     PetscCall(VecNorm(Bx, NORM_INFINITY, &err));
226:     if (err > 10. * PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMultTransposeAdd err %g\n", err));
227:   }
228:   PetscCall(VecDestroy(&Ax));
229:   PetscCall(VecDestroy(&Ay));
230:   PetscCall(VecDestroy(&Bx));
231:   PetscCall(VecDestroy(&By));

233:   /* Test MatMatMult */
234:   if (ldc) PetscCall(PetscMalloc1(nrhs * (n + ldc), &Cdata));
235:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, n, PETSC_DECIDE, N, nrhs, Cdata, &C));
236:   PetscCall(MatDenseSetLDA(C, n + ldc));
237:   PetscCall(MatSetRandom(C, NULL));
238:   if (cgpu) PetscCall(MatConvert(C, MATDENSECUDA, MAT_INPLACE_MATRIX, &C));
239:   for (nt = 0; nt < ntrials; nt++) {
240:     PetscCall(MatMatMult(B, C, nt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, PETSC_DETERMINE, &D));
241:     PetscCall(MatViewFromOptions(D, NULL, "-bc_view"));
242:     PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)D, &flg, MATSEQDENSE, MATMPIDENSE, ""));
243:     if (flg) {
244:       PetscCall(MatCreateVecs(B, &x, &y));
245:       PetscCall(MatCreateVecs(D, NULL, &v));
246:       for (i = 0; i < nrhs; i++) {
247:         PetscCall(MatGetColumnVector(D, v, i));
248:         PetscCall(MatGetColumnVector(C, x, i));
249:         PetscCall(MatMult(B, x, y));
250:         PetscCall(VecAXPY(y, -1.0, v));
251:         PetscCall(VecNorm(y, NORM_INFINITY, &err));
252:         if (err > 10. * PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMat err %" PetscInt_FMT " %g\n", i, err));
253:       }
254:       PetscCall(VecDestroy(&y));
255:       PetscCall(VecDestroy(&x));
256:       PetscCall(VecDestroy(&v));
257:     }
258:   }
259:   PetscCall(MatDestroy(&D));

261:   /* Test MatTransposeMatMult */
262:   if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */
263:     for (nt = 0; nt < ntrials; nt++) {
264:       PetscCall(MatTransposeMatMult(B, C, nt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, PETSC_DETERMINE, &D));
265:       PetscCall(MatViewFromOptions(D, NULL, "-btc_view"));
266:       PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)D, &flg, MATSEQDENSE, MATMPIDENSE, ""));
267:       if (flg) {
268:         PetscCall(MatCreateVecs(B, &y, &x));
269:         PetscCall(MatCreateVecs(D, NULL, &v));
270:         for (i = 0; i < nrhs; i++) {
271:           PetscCall(MatGetColumnVector(D, v, i));
272:           PetscCall(MatGetColumnVector(C, x, i));
273:           PetscCall(MatMultTranspose(B, x, y));
274:           PetscCall(VecAXPY(y, -1.0, v));
275:           PetscCall(VecNorm(y, NORM_INFINITY, &err));
276:           if (err > 10. * PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatTransMat err %" PetscInt_FMT " %g\n", i, err));
277:         }
278:         PetscCall(VecDestroy(&y));
279:         PetscCall(VecDestroy(&x));
280:         PetscCall(VecDestroy(&v));
281:       }
282:     }
283:     PetscCall(MatDestroy(&D));
284:   }

286:   /* Test basis orthogonalization */
287:   if (testorthog) {
288:     PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &D));
289:     PetscCall(MatSetOption(D, MAT_SYMMETRIC, symm));
290:     PetscCall(MatH2OpusOrthogonalize(D));
291:     PetscCall(MatMultEqual(B, D, 10, &flg));
292:     if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMult error after basis ortogonalization\n"));
293:     PetscCall(MatDestroy(&D));
294:   }

296:   /* Test matrix compression */
297:   if (testcompress) {
298:     PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &D));
299:     PetscCall(MatSetOption(D, MAT_SYMMETRIC, symm));
300:     PetscCall(MatH2OpusCompress(D, PETSC_SMALL));
301:     PetscCall(MatDestroy(&D));
302:   }

304:   /* Test low-rank update */
305:   if (testhlru) {
306:     Mat          U, V;
307:     PetscScalar *Udata = NULL, *Vdata = NULL;

309:     if (ldu) {
310:       PetscCall(PetscMalloc1(nlr * (n + ldu), &Udata));
311:       PetscCall(PetscMalloc1(nlr * (n + ldu + 2), &Vdata));
312:     }
313:     PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &D));
314:     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)D), n, PETSC_DECIDE, N, nlr, Udata, &U));
315:     PetscCall(MatDenseSetLDA(U, n + ldu));
316:     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)D), n, PETSC_DECIDE, N, nlr, Vdata, &V));
317:     if (ldu) PetscCall(MatDenseSetLDA(V, n + ldu + 2));
318:     PetscCall(MatSetRandom(U, NULL));
319:     PetscCall(MatSetRandom(V, NULL));
320:     PetscCall(MatH2OpusLowRankUpdate(D, U, V, 0.5));
321:     PetscCall(MatH2OpusLowRankUpdate(D, U, V, -0.5));
322:     PetscCall(MatMultEqual(B, D, 10, &flg));
323:     if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMult error after low-rank update\n"));
324:     PetscCall(MatDestroy(&D));
325:     PetscCall(MatDestroy(&U));
326:     PetscCall(PetscFree(Udata));
327:     PetscCall(MatDestroy(&V));
328:     PetscCall(PetscFree(Vdata));
329:   }

331:   /* check explicit operator */
332:   if (checkexpl) {
333:     Mat Be, Bet;

335:     PetscCall(MatComputeOperator(B, MATDENSE, &D));
336:     PetscCall(MatDuplicate(D, MAT_COPY_VALUES, &Be));
337:     PetscCall(MatNorm(D, NORM_FROBENIUS, &nB));
338:     PetscCall(MatViewFromOptions(D, NULL, "-expl_view"));
339:     PetscCall(MatAXPY(D, -1.0, A, SAME_NONZERO_PATTERN));
340:     PetscCall(MatViewFromOptions(D, NULL, "-diff_view"));
341:     PetscCall(MatNorm(D, NORM_FROBENIUS, &nD));
342:     PetscCall(MatNorm(A, NORM_FROBENIUS, &nA));
343:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Approximation error %g (%g / %g, %g)\n", nD / nA, nD, nA, nB));
344:     PetscCall(MatDestroy(&D));

346:     if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */
347:       PetscCall(MatTranspose(A, MAT_INPLACE_MATRIX, &A));
348:       PetscCall(MatComputeOperatorTranspose(B, MATDENSE, &D));
349:       PetscCall(MatDuplicate(D, MAT_COPY_VALUES, &Bet));
350:       PetscCall(MatNorm(D, NORM_FROBENIUS, &nB));
351:       PetscCall(MatViewFromOptions(D, NULL, "-expl_trans_view"));
352:       PetscCall(MatAXPY(D, -1.0, A, SAME_NONZERO_PATTERN));
353:       PetscCall(MatViewFromOptions(D, NULL, "-diff_trans_view"));
354:       PetscCall(MatNorm(D, NORM_FROBENIUS, &nD));
355:       PetscCall(MatNorm(A, NORM_FROBENIUS, &nA));
356:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Approximation error transpose %g (%g / %g, %g)\n", nD / nA, nD, nA, nB));
357:       PetscCall(MatDestroy(&D));

359:       PetscCall(MatTranspose(Bet, MAT_INPLACE_MATRIX, &Bet));
360:       PetscCall(MatAXPY(Be, -1.0, Bet, SAME_NONZERO_PATTERN));
361:       PetscCall(MatViewFromOptions(Be, NULL, "-diff_expl_view"));
362:       PetscCall(MatNorm(Be, NORM_FROBENIUS, &nB));
363:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Approximation error B - (B^T)^T %g\n", nB));
364:       PetscCall(MatDestroy(&Be));
365:       PetscCall(MatDestroy(&Bet));
366:     }
367:   }
368:   PetscCall(MatDestroy(&A));
369:   PetscCall(MatDestroy(&B));
370:   PetscCall(MatDestroy(&C));
371:   PetscCall(PetscFree(Cdata));
372:   PetscCall(PetscFree(Adata));
373:   PetscCall(PetscFinalize());
374:   return 0;
375: }

377: /*TEST

379:    build:
380:      requires: h2opus

382: #tests from kernel
383:    test:
384:      requires: h2opus
385:      nsize: 1
386:      suffix: 1
387:      args: -n {{17 33}} -kernel 1 -dim {{1 2 3}} -symm {{0 1}} -checkexpl -bgpu 0

389:    test:
390:      requires: h2opus
391:      nsize: 1
392:      suffix: 1_ld
393:      output_file: output/ex66_1.out
394:      args: -n 33 -kernel 1 -dim 1 -lda 13 -ldc 11 -ldu 17 -symm 0 -checkexpl -bgpu 0

396:    test:
397:      requires: h2opus cuda
398:      nsize: 1
399:      suffix: 1_cuda
400:      output_file: output/ex66_1.out
401:      args: -n {{17 33}} -kernel 1 -dim {{1 2 3}} -symm {{0 1}} -checkexpl -bgpu 1

403:    test:
404:      requires: h2opus cuda
405:      nsize: 1
406:      suffix: 1_cuda_ld
407:      output_file: output/ex66_1.out
408:      args: -n 33 -kernel 1 -dim 1 -lda 13 -ldc 11 -ldu 17 -symm 0 -checkexpl -bgpu 1

410:    test:
411:      requires: h2opus
412:      nsize: {{2 3}}
413:      suffix: 1_par
414:      args: -n 64 -symm -kernel 1 -dim 1 -ldc 12 -testlayout {{0 1}} -bgpu 0 -cgpu 0

416:    test:
417:      requires: h2opus cuda
418:      nsize: {{2 3}}
419:      suffix: 1_par_cuda
420:      args: -n 64 -symm -kernel 1 -dim 1 -ldc 12 -testlayout {{0 1}} -bgpu {{0 1}} -cgpu {{0 1}}
421:      output_file: output/ex66_1_par.out

423: #tests from matrix sampling (parallel or unsymmetric not supported)
424:    test:
425:      requires: h2opus
426:      nsize: 1
427:      suffix: 2
428:      args: -n {{17 33}} -kernel 0 -dim 2 -symm 1 -checkexpl -bgpu 0

430:    test:
431:      requires: h2opus cuda
432:      nsize: 1
433:      suffix: 2_cuda
434:      output_file: output/ex66_2.out
435:      args: -n {{17 29}} -kernel 0 -dim 2 -symm 1 -checkexpl -bgpu {{0 1}} -agpu {{0 1}}

437: #tests view operation
438:    test:
439:      requires: h2opus !cuda
440:      filter: grep -v " MPI process" | grep -v "\[" | grep -v "\]"
441:      nsize: {{1 2 3}}
442:      suffix: view
443:      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

445:    test:
446:      requires: h2opus cuda
447:      filter: grep -v " MPI process" | grep -v "\[" | grep -v "\]"
448:      nsize: {{1 2 3}}
449:      suffix: view_cuda
450:      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

452: TEST*/