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