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) { 99 PetscCall(PetscMalloc1(N*(n+lda),&Adata)); 100 } 101 PetscCall(MatCreateDense(PETSC_COMM_WORLD,n,n,N,N,Adata,&A)); 102 PetscCall(MatDenseSetLDA(A,n+lda)); 103 104 /* Create random points; these are replicated in order to populate a dense matrix and to compare sequential and dense runs 105 The constructor for MATH2OPUS can take as input the distributed coordinates and replicates them internally in case 106 a kernel construction is requested */ 107 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&r)); 108 PetscCall(PetscRandomSetFromOptions(r)); 109 PetscCall(PetscRandomSetSeed(r,123456)); 110 PetscCall(PetscRandomSeed(r)); 111 PetscCall(PetscMalloc1(N*dim,&coords)); 112 PetscCall(PetscRandomGetValuesReal(r,N*dim,coords)); 113 PetscCall(PetscRandomDestroy(&r)); 114 115 if (kernel || !randommat) { 116 MatH2OpusKernel k = Asymm ? GenEntry_Symm : GenEntry_Unsymm; 117 PetscInt ist,ien; 118 119 PetscCall(MatGetOwnershipRange(A,&ist,&ien)); 120 for (i = ist; i < ien; i++) { 121 for (j = 0; j < N; j++) { 122 PetscCall(MatSetValue(A,i,j,(*k)(dim,coords + i*dim,coords + j*dim,NULL),INSERT_VALUES)); 123 } 124 } 125 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 126 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 127 if (kernel) { 128 PetscCall(MatCreateH2OpusFromKernel(PETSC_COMM_WORLD,n,n,N,N,dim,coords + ist*dim,PETSC_TRUE,k,NULL,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,&B)); 129 } else { 130 PetscCall(MatCreateH2OpusFromMat(A,dim,coords + ist*dim,PETSC_TRUE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,&B)); 131 } 132 } else { 133 PetscInt ist; 134 135 PetscCall(MatGetOwnershipRange(A,&ist,NULL)); 136 PetscCall(MatSetRandom(A,NULL)); 137 if (Asymm) { 138 PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX,&B)); 139 PetscCall(MatAXPY(A,1.0,B,SAME_NONZERO_PATTERN)); 140 PetscCall(MatDestroy(&B)); 141 PetscCall(MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE)); 142 } 143 PetscCall(MatCreateH2OpusFromMat(A,dim,coords + ist*dim,PETSC_TRUE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,&B)); 144 } 145 PetscCall(PetscFree(coords)); 146 if (agpu) { 147 PetscCall(MatConvert(A,MATDENSECUDA,MAT_INPLACE_MATRIX,&A)); 148 } 149 PetscCall(MatViewFromOptions(A,NULL,"-A_view")); 150 151 PetscCall(MatSetOption(B,MAT_SYMMETRIC,symm)); 152 153 /* assemble the H-matrix */ 154 PetscCall(MatBindToCPU(B,(PetscBool)!bgpu)); 155 PetscCall(MatSetFromOptions(B)); 156 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 157 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 158 PetscCall(MatViewFromOptions(B,NULL,"-B_view")); 159 160 /* Test MatScale */ 161 PetscCall(MatScale(A,scale)); 162 PetscCall(MatScale(B,scale)); 163 164 /* Test MatMult */ 165 PetscCall(MatCreateVecs(A,&Ax,&Ay)); 166 PetscCall(MatCreateVecs(B,&Bx,&By)); 167 PetscCall(VecSetRandom(Ax,NULL)); 168 PetscCall(VecCopy(Ax,Bx)); 169 PetscCall(MatMult(A,Ax,Ay)); 170 PetscCall(MatMult(B,Bx,By)); 171 PetscCall(VecViewFromOptions(Ay,NULL,"-mult_vec_view")); 172 PetscCall(VecViewFromOptions(By,NULL,"-mult_vec_view")); 173 PetscCall(VecNorm(Ay,NORM_INFINITY,&nX)); 174 PetscCall(VecAXPY(Ay,-1.0,By)); 175 PetscCall(VecViewFromOptions(Ay,NULL,"-mult_vec_view")); 176 PetscCall(VecNorm(Ay,NORM_INFINITY,&err)); 177 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMult err %g\n",err/nX)); 178 PetscCall(VecScale(By,-1.0)); 179 PetscCall(MatMultAdd(B,Bx,By,By)); 180 PetscCall(VecNorm(By,NORM_INFINITY,&err)); 181 PetscCall(VecViewFromOptions(By,NULL,"-mult_vec_view")); 182 if (err > 10.*PETSC_SMALL) { 183 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMultAdd err %g\n",err)); 184 } 185 186 /* Test MatNorm */ 187 PetscCall(MatNorm(A,NORM_INFINITY,&norms[0])); 188 PetscCall(MatNorm(A,NORM_1,&norms[1])); 189 norms[2] = -1.; /* NORM_2 not supported */ 190 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])); 191 PetscCall(MatGetOperation(A,MATOP_NORM,&Anormfunc)); 192 PetscCall(MatGetOperation(B,MATOP_NORM,&approxnormfunc)); 193 PetscCall(MatSetOperation(A,MATOP_NORM,approxnormfunc)); 194 PetscCall(MatNorm(A,NORM_INFINITY,&norms[0])); 195 PetscCall(MatNorm(A,NORM_1,&norms[1])); 196 PetscCall(MatNorm(A,NORM_2,&norms[2])); 197 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])); 198 if (testnorm) { 199 PetscCall(MatNorm(B,NORM_INFINITY,&norms[0])); 200 PetscCall(MatNorm(B,NORM_1,&norms[1])); 201 PetscCall(MatNorm(B,NORM_2,&norms[2])); 202 } else { 203 norms[0] = -1.; 204 norms[1] = -1.; 205 norms[2] = -1.; 206 } 207 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])); 208 PetscCall(MatSetOperation(A,MATOP_NORM,Anormfunc)); 209 210 /* Test MatDuplicate */ 211 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&D)); 212 PetscCall(MatSetOption(D,MAT_SYMMETRIC,symm)); 213 PetscCall(MatMultEqual(B,D,10,&flg)); 214 if (!flg) { 215 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMult error after MatDuplicate\n")); 216 } 217 if (testtrans) { 218 PetscCall(MatMultTransposeEqual(B,D,10,&flg)); 219 if (!flg) { 220 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMultTranspose error after MatDuplicate\n")); 221 } 222 } 223 PetscCall(MatDestroy(&D)); 224 225 if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */ 226 PetscCall(VecSetRandom(Ay,NULL)); 227 PetscCall(VecCopy(Ay,By)); 228 PetscCall(MatMultTranspose(A,Ay,Ax)); 229 PetscCall(MatMultTranspose(B,By,Bx)); 230 PetscCall(VecViewFromOptions(Ax,NULL,"-multtrans_vec_view")); 231 PetscCall(VecViewFromOptions(Bx,NULL,"-multtrans_vec_view")); 232 PetscCall(VecNorm(Ax,NORM_INFINITY,&nX)); 233 PetscCall(VecAXPY(Ax,-1.0,Bx)); 234 PetscCall(VecViewFromOptions(Ax,NULL,"-multtrans_vec_view")); 235 PetscCall(VecNorm(Ax,NORM_INFINITY,&err)); 236 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMultTranspose err %g\n",err/nX)); 237 PetscCall(VecScale(Bx,-1.0)); 238 PetscCall(MatMultTransposeAdd(B,By,Bx,Bx)); 239 PetscCall(VecNorm(Bx,NORM_INFINITY,&err)); 240 if (err > 10.*PETSC_SMALL) { 241 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMultTransposeAdd err %g\n",err)); 242 } 243 } 244 PetscCall(VecDestroy(&Ax)); 245 PetscCall(VecDestroy(&Ay)); 246 PetscCall(VecDestroy(&Bx)); 247 PetscCall(VecDestroy(&By)); 248 249 /* Test MatMatMult */ 250 if (ldc) { 251 PetscCall(PetscMalloc1(nrhs*(n+ldc),&Cdata)); 252 } 253 PetscCall(MatCreateDense(PETSC_COMM_WORLD,n,PETSC_DECIDE,N,nrhs,Cdata,&C)); 254 PetscCall(MatDenseSetLDA(C,n+ldc)); 255 PetscCall(MatSetRandom(C,NULL)); 256 if (cgpu) { 257 PetscCall(MatConvert(C,MATDENSECUDA,MAT_INPLACE_MATRIX,&C)); 258 } 259 for (nt = 0; nt < ntrials; nt++) { 260 PetscCall(MatMatMult(B,C,nt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D)); 261 PetscCall(MatViewFromOptions(D,NULL,"-bc_view")); 262 PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)D,&flg,MATSEQDENSE,MATMPIDENSE,"")); 263 if (flg) { 264 PetscCall(MatCreateVecs(B,&x,&y)); 265 PetscCall(MatCreateVecs(D,NULL,&v)); 266 for (i = 0; i < nrhs; i++) { 267 PetscCall(MatGetColumnVector(D,v,i)); 268 PetscCall(MatGetColumnVector(C,x,i)); 269 PetscCall(MatMult(B,x,y)); 270 PetscCall(VecAXPY(y,-1.0,v)); 271 PetscCall(VecNorm(y,NORM_INFINITY,&err)); 272 if (err > 10.*PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMat err %" PetscInt_FMT " %g\n",i,err)); 273 } 274 PetscCall(VecDestroy(&y)); 275 PetscCall(VecDestroy(&x)); 276 PetscCall(VecDestroy(&v)); 277 } 278 } 279 PetscCall(MatDestroy(&D)); 280 281 /* Test MatTransposeMatMult */ 282 if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */ 283 for (nt = 0; nt < ntrials; nt++) { 284 PetscCall(MatTransposeMatMult(B,C,nt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D)); 285 PetscCall(MatViewFromOptions(D,NULL,"-btc_view")); 286 PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)D,&flg,MATSEQDENSE,MATMPIDENSE,"")); 287 if (flg) { 288 PetscCall(MatCreateVecs(B,&y,&x)); 289 PetscCall(MatCreateVecs(D,NULL,&v)); 290 for (i = 0; i < nrhs; i++) { 291 PetscCall(MatGetColumnVector(D,v,i)); 292 PetscCall(MatGetColumnVector(C,x,i)); 293 PetscCall(MatMultTranspose(B,x,y)); 294 PetscCall(VecAXPY(y,-1.0,v)); 295 PetscCall(VecNorm(y,NORM_INFINITY,&err)); 296 if (err > 10.*PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatTransMat err %" PetscInt_FMT " %g\n",i,err)); 297 } 298 PetscCall(VecDestroy(&y)); 299 PetscCall(VecDestroy(&x)); 300 PetscCall(VecDestroy(&v)); 301 } 302 } 303 PetscCall(MatDestroy(&D)); 304 } 305 306 /* Test basis orthogonalization */ 307 if (testorthog) { 308 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&D)); 309 PetscCall(MatSetOption(D,MAT_SYMMETRIC,symm)); 310 PetscCall(MatH2OpusOrthogonalize(D)); 311 PetscCall(MatMultEqual(B,D,10,&flg)); 312 if (!flg) { 313 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMult error after basis ortogonalization\n")); 314 } 315 PetscCall(MatDestroy(&D)); 316 } 317 318 /* Test matrix compression */ 319 if (testcompress) { 320 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&D)); 321 PetscCall(MatSetOption(D,MAT_SYMMETRIC,symm)); 322 PetscCall(MatH2OpusCompress(D,PETSC_SMALL)); 323 PetscCall(MatDestroy(&D)); 324 } 325 326 /* Test low-rank update */ 327 if (testhlru) { 328 Mat U, V; 329 PetscScalar *Udata = NULL, *Vdata = NULL; 330 331 if (ldu) { 332 PetscCall(PetscMalloc1(nlr*(n+ldu),&Udata)); 333 PetscCall(PetscMalloc1(nlr*(n+ldu+2),&Vdata)); 334 } 335 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&D)); 336 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)D),n,PETSC_DECIDE,N,nlr,Udata,&U)); 337 PetscCall(MatDenseSetLDA(U,n+ldu)); 338 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)D),n,PETSC_DECIDE,N,nlr,Vdata,&V)); 339 if (ldu) PetscCall(MatDenseSetLDA(V,n+ldu+2)); 340 PetscCall(MatSetRandom(U,NULL)); 341 PetscCall(MatSetRandom(V,NULL)); 342 PetscCall(MatH2OpusLowRankUpdate(D,U,V,0.5)); 343 PetscCall(MatH2OpusLowRankUpdate(D,U,V,-0.5)); 344 PetscCall(MatMultEqual(B,D,10,&flg)); 345 if (!flg) { 346 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMult error after low-rank update\n")); 347 } 348 PetscCall(MatDestroy(&D)); 349 PetscCall(MatDestroy(&U)); 350 PetscCall(PetscFree(Udata)); 351 PetscCall(MatDestroy(&V)); 352 PetscCall(PetscFree(Vdata)); 353 } 354 355 /* check explicit operator */ 356 if (checkexpl) { 357 Mat Be, Bet; 358 359 PetscCall(MatComputeOperator(B,MATDENSE,&D)); 360 PetscCall(MatDuplicate(D,MAT_COPY_VALUES,&Be)); 361 PetscCall(MatNorm(D,NORM_FROBENIUS,&nB)); 362 PetscCall(MatViewFromOptions(D,NULL,"-expl_view")); 363 PetscCall(MatAXPY(D,-1.0,A,SAME_NONZERO_PATTERN)); 364 PetscCall(MatViewFromOptions(D,NULL,"-diff_view")); 365 PetscCall(MatNorm(D,NORM_FROBENIUS,&nD)); 366 PetscCall(MatNorm(A,NORM_FROBENIUS,&nA)); 367 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Approximation error %g (%g / %g, %g)\n",nD/nA,nD,nA,nB)); 368 PetscCall(MatDestroy(&D)); 369 370 if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */ 371 PetscCall(MatTranspose(A,MAT_INPLACE_MATRIX,&A)); 372 PetscCall(MatComputeOperatorTranspose(B,MATDENSE,&D)); 373 PetscCall(MatDuplicate(D,MAT_COPY_VALUES,&Bet)); 374 PetscCall(MatNorm(D,NORM_FROBENIUS,&nB)); 375 PetscCall(MatViewFromOptions(D,NULL,"-expl_trans_view")); 376 PetscCall(MatAXPY(D,-1.0,A,SAME_NONZERO_PATTERN)); 377 PetscCall(MatViewFromOptions(D,NULL,"-diff_trans_view")); 378 PetscCall(MatNorm(D,NORM_FROBENIUS,&nD)); 379 PetscCall(MatNorm(A,NORM_FROBENIUS,&nA)); 380 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Approximation error transpose %g (%g / %g, %g)\n",nD/nA,nD,nA,nB)); 381 PetscCall(MatDestroy(&D)); 382 383 PetscCall(MatTranspose(Bet,MAT_INPLACE_MATRIX,&Bet)); 384 PetscCall(MatAXPY(Be,-1.0,Bet,SAME_NONZERO_PATTERN)); 385 PetscCall(MatViewFromOptions(Be,NULL,"-diff_expl_view")); 386 PetscCall(MatNorm(Be,NORM_FROBENIUS,&nB)); 387 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Approximation error B - (B^T)^T %g\n",nB)); 388 PetscCall(MatDestroy(&Be)); 389 PetscCall(MatDestroy(&Bet)); 390 } 391 } 392 PetscCall(MatDestroy(&A)); 393 PetscCall(MatDestroy(&B)); 394 PetscCall(MatDestroy(&C)); 395 PetscCall(PetscFree(Cdata)); 396 PetscCall(PetscFree(Adata)); 397 PetscCall(PetscFinalize()); 398 return 0; 399 } 400 401 /*TEST 402 403 build: 404 requires: h2opus 405 406 #tests from kernel 407 test: 408 requires: h2opus 409 nsize: 1 410 suffix: 1 411 args: -n {{17 33}} -kernel 1 -dim {{1 2 3}} -symm {{0 1}} -checkexpl -bgpu 0 412 413 test: 414 requires: h2opus 415 nsize: 1 416 suffix: 1_ld 417 output_file: output/ex66_1.out 418 args: -n 33 -kernel 1 -dim 1 -lda 13 -ldc 11 -ldu 17 -symm 0 -checkexpl -bgpu 0 419 420 test: 421 requires: h2opus cuda 422 nsize: 1 423 suffix: 1_cuda 424 output_file: output/ex66_1.out 425 args: -n {{17 33}} -kernel 1 -dim {{1 2 3}} -symm {{0 1}} -checkexpl -bgpu 1 426 427 test: 428 requires: h2opus cuda 429 nsize: 1 430 suffix: 1_cuda_ld 431 output_file: output/ex66_1.out 432 args: -n 33 -kernel 1 -dim 1 -lda 13 -ldc 11 -ldu 17 -symm 0 -checkexpl -bgpu 1 433 434 test: 435 requires: h2opus 436 nsize: {{2 3}} 437 suffix: 1_par 438 args: -n 64 -symm -kernel 1 -dim 1 -ldc 12 -testlayout {{0 1}} -bgpu 0 -cgpu 0 439 440 test: 441 requires: h2opus cuda 442 nsize: {{2 3}} 443 suffix: 1_par_cuda 444 args: -n 64 -symm -kernel 1 -dim 1 -ldc 12 -testlayout {{0 1}} -bgpu {{0 1}} -cgpu {{0 1}} 445 output_file: output/ex66_1_par.out 446 447 #tests from matrix sampling (parallel or unsymmetric not supported) 448 test: 449 requires: h2opus 450 nsize: 1 451 suffix: 2 452 args: -n {{17 33}} -kernel 0 -dim 2 -symm 1 -checkexpl -bgpu 0 453 454 test: 455 requires: h2opus cuda 456 nsize: 1 457 suffix: 2_cuda 458 output_file: output/ex66_2.out 459 args: -n {{17 33}} -kernel 0 -dim 2 -symm 1 -checkexpl -bgpu {{0 1}} -agpu {{0 1}} 460 461 #tests view operation 462 test: 463 requires: h2opus !cuda 464 filter: grep -v " MPI process" | grep -v "\[" | grep -v "\]" 465 nsize: {{1 2 3}} 466 suffix: view 467 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 468 469 test: 470 requires: h2opus cuda 471 filter: grep -v " MPI process" | grep -v "\[" | grep -v "\]" 472 nsize: {{1 2 3}} 473 suffix: view_cuda 474 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 475 476 TEST*/ 477