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