153022affSStefano Zampini static char help[] = "Tests MATH2OPUS\n\n"; 2c4fb69feSStefano Zampini 3c4fb69feSStefano Zampini #include <petscmat.h> 4c4fb69feSStefano Zampini #include <petscsf.h> 5c4fb69feSStefano Zampini 6c4fb69feSStefano Zampini static PetscScalar GenEntry_Symm(PetscInt sdim, PetscReal x[], PetscReal y[], void *ctx) 7c4fb69feSStefano Zampini { 8c4fb69feSStefano Zampini PetscInt d; 9c4fb69feSStefano Zampini PetscReal clength = sdim == 3 ? 0.2 : 0.1; 10c4fb69feSStefano Zampini PetscReal dist, diff = 0.0; 11c4fb69feSStefano Zampini 12c4fb69feSStefano Zampini for (d = 0; d < sdim; d++) { diff += (x[d] - y[d]) * (x[d] - y[d]); } 13c4fb69feSStefano Zampini dist = PetscSqrtReal(diff); 14c4fb69feSStefano Zampini return PetscExpReal(-dist / clength); 15c4fb69feSStefano Zampini } 16c4fb69feSStefano Zampini 17c4fb69feSStefano Zampini static PetscScalar GenEntry_Unsymm(PetscInt sdim, PetscReal x[], PetscReal y[], void *ctx) 18c4fb69feSStefano Zampini { 19c4fb69feSStefano Zampini PetscInt d; 20c4fb69feSStefano Zampini PetscReal clength = sdim == 3 ? 0.2 : 0.1; 21c4fb69feSStefano Zampini PetscReal dist, diff = 0.0, nx = 0.0, ny = 0.0; 22c4fb69feSStefano Zampini 23c4fb69feSStefano Zampini for (d = 0; d < sdim; d++) { nx += x[d]*x[d]; } 24c4fb69feSStefano Zampini for (d = 0; d < sdim; d++) { ny += y[d]*y[d]; } 25c4fb69feSStefano Zampini for (d = 0; d < sdim; d++) { diff += (x[d] - y[d]) * (x[d] - y[d]); } 26c4fb69feSStefano Zampini dist = PetscSqrtReal(diff); 27c4fb69feSStefano Zampini return nx > ny ? PetscExpReal(-dist / clength) : PetscExpReal(-dist / clength) + 1.; 28c4fb69feSStefano Zampini } 29c4fb69feSStefano Zampini 30c4fb69feSStefano Zampini int main(int argc,char **argv) 31c4fb69feSStefano Zampini { 32c4fb69feSStefano Zampini Mat A,B,C,D; 33c4fb69feSStefano Zampini Vec v,x,y,Ax,Ay,Bx,By; 34c4fb69feSStefano Zampini PetscRandom r; 35c4fb69feSStefano Zampini PetscLayout map; 36c4fb69feSStefano Zampini PetscScalar *Adata = NULL, *Cdata = NULL, scale = 1.0; 37c4fb69feSStefano Zampini PetscReal *coords,nA,nD,nB,err,nX,norms[3]; 38300d917bSStefano Zampini PetscInt N, n = 64, dim = 1, i, j, nrhs = 11, lda = 0, ldc = 0, ldu = 0, nlr = 7, nt, ntrials = 2; 39c4fb69feSStefano Zampini PetscMPIInt size,rank; 40c4fb69feSStefano Zampini PetscBool testlayout = PETSC_FALSE, flg, symm = PETSC_FALSE, Asymm = PETSC_TRUE, kernel = PETSC_TRUE; 4153022affSStefano Zampini PetscBool checkexpl = PETSC_FALSE, agpu = PETSC_FALSE, bgpu = PETSC_FALSE, cgpu = PETSC_FALSE, flgglob; 42300d917bSStefano Zampini PetscBool testtrans, testnorm, randommat = PETSC_TRUE, testorthog, testcompress, testhlru; 43c4fb69feSStefano Zampini void (*approxnormfunc)(void); 44c4fb69feSStefano Zampini void (*Anormfunc)(void); 45c4fb69feSStefano Zampini 4653022affSStefano Zampini #if defined(PETSC_HAVE_MPI_INIT_THREAD) 47c4fb69feSStefano Zampini PETSC_MPI_THREAD_REQUIRED = MPI_THREAD_MULTIPLE; 4853022affSStefano Zampini #endif 49*327415f7SBarry Smith PetscFunctionBeginUser; 509566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc,&argv,(char*) 0,help)); 519566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-ng",&N,&flgglob)); 529566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 539566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL)); 549566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL)); 559566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-lda",&lda,NULL)); 569566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-ldc",&ldc,NULL)); 579566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-nlr",&nlr,NULL)); 589566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-ldu",&ldu,NULL)); 599566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-matmattrials",&ntrials,NULL)); 609566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-randommat",&randommat,NULL)); 619566063dSJacob Faibussowitsch if (!flgglob) PetscCall(PetscOptionsGetBool(NULL,NULL,"-testlayout",&testlayout,NULL)); 629566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-Asymm",&Asymm,NULL)); 639566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-symm",&symm,NULL)); 649566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-kernel",&kernel,NULL)); 659566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-checkexpl",&checkexpl,NULL)); 669566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-agpu",&agpu,NULL)); 679566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-bgpu",&bgpu,NULL)); 689566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-cgpu",&cgpu,NULL)); 699566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetScalar(NULL,NULL,"-scale",&scale,NULL)); 70c4fb69feSStefano Zampini if (!Asymm) symm = PETSC_FALSE; 71c4fb69feSStefano Zampini 729566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 73300d917bSStefano Zampini 74300d917bSStefano Zampini /* Disable tests for unimplemented variants */ 75c4fb69feSStefano Zampini testtrans = (PetscBool)(size == 1 || symm); 76c4fb69feSStefano Zampini testnorm = (PetscBool)(size == 1 || symm); 7753022affSStefano Zampini testorthog = (PetscBool)(size == 1 || symm); 7853022affSStefano Zampini testcompress = (PetscBool)(size == 1 || symm); 79300d917bSStefano Zampini testhlru = (PetscBool)(size == 1); 80c4fb69feSStefano Zampini 819566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 829566063dSJacob Faibussowitsch PetscCall(PetscLayoutCreate(PETSC_COMM_WORLD,&map)); 83c4fb69feSStefano Zampini if (testlayout) { 84c4fb69feSStefano Zampini if (rank%2) n = PetscMax(2*n-5*rank,0); 85c4fb69feSStefano Zampini else n = 2*n+rank; 86c4fb69feSStefano Zampini } 8753022affSStefano Zampini if (!flgglob) { 889566063dSJacob Faibussowitsch PetscCall(PetscLayoutSetLocalSize(map,n)); 899566063dSJacob Faibussowitsch PetscCall(PetscLayoutSetUp(map)); 909566063dSJacob Faibussowitsch PetscCall(PetscLayoutGetSize(map,&N)); 9153022affSStefano Zampini } else { 929566063dSJacob Faibussowitsch PetscCall(PetscLayoutSetSize(map,N)); 939566063dSJacob Faibussowitsch PetscCall(PetscLayoutSetUp(map)); 949566063dSJacob Faibussowitsch PetscCall(PetscLayoutGetLocalSize(map,&n)); 9553022affSStefano Zampini } 969566063dSJacob Faibussowitsch PetscCall(PetscLayoutDestroy(&map)); 97c4fb69feSStefano Zampini 98c4fb69feSStefano Zampini if (lda) { 999566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(N*(n+lda),&Adata)); 100c4fb69feSStefano Zampini } 1019566063dSJacob Faibussowitsch PetscCall(MatCreateDense(PETSC_COMM_WORLD,n,n,N,N,Adata,&A)); 1029566063dSJacob Faibussowitsch PetscCall(MatDenseSetLDA(A,n+lda)); 103c4fb69feSStefano Zampini 10453022affSStefano Zampini /* Create random points; these are replicated in order to populate a dense matrix and to compare sequential and dense runs 10553022affSStefano Zampini The constructor for MATH2OPUS can take as input the distributed coordinates and replicates them internally in case 10653022affSStefano Zampini a kernel construction is requested */ 1079566063dSJacob Faibussowitsch PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&r)); 1089566063dSJacob Faibussowitsch PetscCall(PetscRandomSetFromOptions(r)); 1099566063dSJacob Faibussowitsch PetscCall(PetscRandomSetSeed(r,123456)); 1109566063dSJacob Faibussowitsch PetscCall(PetscRandomSeed(r)); 1119566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(N*dim,&coords)); 1129566063dSJacob Faibussowitsch PetscCall(PetscRandomGetValuesReal(r,N*dim,coords)); 1139566063dSJacob Faibussowitsch PetscCall(PetscRandomDestroy(&r)); 11453022affSStefano Zampini 115c4fb69feSStefano Zampini if (kernel || !randommat) { 11653022affSStefano Zampini MatH2OpusKernel k = Asymm ? GenEntry_Symm : GenEntry_Unsymm; 117c4fb69feSStefano Zampini PetscInt ist,ien; 118c4fb69feSStefano Zampini 1199566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRange(A,&ist,&ien)); 120c4fb69feSStefano Zampini for (i = ist; i < ien; i++) { 121c4fb69feSStefano Zampini for (j = 0; j < N; j++) { 1229566063dSJacob Faibussowitsch PetscCall(MatSetValue(A,i,j,(*k)(dim,coords + i*dim,coords + j*dim,NULL),INSERT_VALUES)); 123c4fb69feSStefano Zampini } 124c4fb69feSStefano Zampini } 1259566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1269566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 127c4fb69feSStefano Zampini if (kernel) { 1289566063dSJacob Faibussowitsch PetscCall(MatCreateH2OpusFromKernel(PETSC_COMM_WORLD,n,n,N,N,dim,coords + ist*dim,PETSC_TRUE,k,NULL,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,&B)); 129c4fb69feSStefano Zampini } else { 1309566063dSJacob Faibussowitsch PetscCall(MatCreateH2OpusFromMat(A,dim,coords + ist*dim,PETSC_TRUE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,&B)); 131c4fb69feSStefano Zampini } 132c4fb69feSStefano Zampini } else { 13353022affSStefano Zampini PetscInt ist; 13453022affSStefano Zampini 1359566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRange(A,&ist,NULL)); 1369566063dSJacob Faibussowitsch PetscCall(MatSetRandom(A,NULL)); 137c4fb69feSStefano Zampini if (Asymm) { 1389566063dSJacob Faibussowitsch PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX,&B)); 1399566063dSJacob Faibussowitsch PetscCall(MatAXPY(A,1.0,B,SAME_NONZERO_PATTERN)); 1409566063dSJacob Faibussowitsch PetscCall(MatDestroy(&B)); 1419566063dSJacob Faibussowitsch PetscCall(MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE)); 142c4fb69feSStefano Zampini } 1439566063dSJacob Faibussowitsch PetscCall(MatCreateH2OpusFromMat(A,dim,coords + ist*dim,PETSC_TRUE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,&B)); 144c4fb69feSStefano Zampini } 1459566063dSJacob Faibussowitsch PetscCall(PetscFree(coords)); 146c4fb69feSStefano Zampini if (agpu) { 1479566063dSJacob Faibussowitsch PetscCall(MatConvert(A,MATDENSECUDA,MAT_INPLACE_MATRIX,&A)); 148c4fb69feSStefano Zampini } 1499566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(A,NULL,"-A_view")); 150c4fb69feSStefano Zampini 1519566063dSJacob Faibussowitsch PetscCall(MatSetOption(B,MAT_SYMMETRIC,symm)); 152c4fb69feSStefano Zampini 153c4fb69feSStefano Zampini /* assemble the H-matrix */ 1549566063dSJacob Faibussowitsch PetscCall(MatBindToCPU(B,(PetscBool)!bgpu)); 1559566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(B)); 1569566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1579566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 1589566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(B,NULL,"-B_view")); 159c4fb69feSStefano Zampini 160c4fb69feSStefano Zampini /* Test MatScale */ 1619566063dSJacob Faibussowitsch PetscCall(MatScale(A,scale)); 1629566063dSJacob Faibussowitsch PetscCall(MatScale(B,scale)); 163c4fb69feSStefano Zampini 164c4fb69feSStefano Zampini /* Test MatMult */ 1659566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(A,&Ax,&Ay)); 1669566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(B,&Bx,&By)); 1679566063dSJacob Faibussowitsch PetscCall(VecSetRandom(Ax,NULL)); 1689566063dSJacob Faibussowitsch PetscCall(VecCopy(Ax,Bx)); 1699566063dSJacob Faibussowitsch PetscCall(MatMult(A,Ax,Ay)); 1709566063dSJacob Faibussowitsch PetscCall(MatMult(B,Bx,By)); 1719566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(Ay,NULL,"-mult_vec_view")); 1729566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(By,NULL,"-mult_vec_view")); 1739566063dSJacob Faibussowitsch PetscCall(VecNorm(Ay,NORM_INFINITY,&nX)); 1749566063dSJacob Faibussowitsch PetscCall(VecAXPY(Ay,-1.0,By)); 1759566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(Ay,NULL,"-mult_vec_view")); 1769566063dSJacob Faibussowitsch PetscCall(VecNorm(Ay,NORM_INFINITY,&err)); 1779566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMult err %g\n",err/nX)); 1789566063dSJacob Faibussowitsch PetscCall(VecScale(By,-1.0)); 1799566063dSJacob Faibussowitsch PetscCall(MatMultAdd(B,Bx,By,By)); 1809566063dSJacob Faibussowitsch PetscCall(VecNorm(By,NORM_INFINITY,&err)); 1819566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(By,NULL,"-mult_vec_view")); 1824c3a5f3cSStefano Zampini if (err > 10.*PETSC_SMALL) { 1839566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMultAdd err %g\n",err)); 184c4fb69feSStefano Zampini } 185c4fb69feSStefano Zampini 186c4fb69feSStefano Zampini /* Test MatNorm */ 1879566063dSJacob Faibussowitsch PetscCall(MatNorm(A,NORM_INFINITY,&norms[0])); 1889566063dSJacob Faibussowitsch PetscCall(MatNorm(A,NORM_1,&norms[1])); 189c4fb69feSStefano Zampini norms[2] = -1.; /* NORM_2 not supported */ 1909566063dSJacob Faibussowitsch 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])); 1919566063dSJacob Faibussowitsch PetscCall(MatGetOperation(A,MATOP_NORM,&Anormfunc)); 1929566063dSJacob Faibussowitsch PetscCall(MatGetOperation(B,MATOP_NORM,&approxnormfunc)); 1939566063dSJacob Faibussowitsch PetscCall(MatSetOperation(A,MATOP_NORM,approxnormfunc)); 1949566063dSJacob Faibussowitsch PetscCall(MatNorm(A,NORM_INFINITY,&norms[0])); 1959566063dSJacob Faibussowitsch PetscCall(MatNorm(A,NORM_1,&norms[1])); 1969566063dSJacob Faibussowitsch PetscCall(MatNorm(A,NORM_2,&norms[2])); 1979566063dSJacob Faibussowitsch 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])); 198c4fb69feSStefano Zampini if (testnorm) { 1999566063dSJacob Faibussowitsch PetscCall(MatNorm(B,NORM_INFINITY,&norms[0])); 2009566063dSJacob Faibussowitsch PetscCall(MatNorm(B,NORM_1,&norms[1])); 2019566063dSJacob Faibussowitsch PetscCall(MatNorm(B,NORM_2,&norms[2])); 202c4fb69feSStefano Zampini } else { 203c4fb69feSStefano Zampini norms[0] = -1.; 204c4fb69feSStefano Zampini norms[1] = -1.; 205c4fb69feSStefano Zampini norms[2] = -1.; 206c4fb69feSStefano Zampini } 2079566063dSJacob Faibussowitsch 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])); 2089566063dSJacob Faibussowitsch PetscCall(MatSetOperation(A,MATOP_NORM,Anormfunc)); 209c4fb69feSStefano Zampini 210c4fb69feSStefano Zampini /* Test MatDuplicate */ 2119566063dSJacob Faibussowitsch PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&D)); 2129566063dSJacob Faibussowitsch PetscCall(MatSetOption(D,MAT_SYMMETRIC,symm)); 2139566063dSJacob Faibussowitsch PetscCall(MatMultEqual(B,D,10,&flg)); 214c4fb69feSStefano Zampini if (!flg) { 2159566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMult error after MatDuplicate\n")); 216c4fb69feSStefano Zampini } 217c4fb69feSStefano Zampini if (testtrans) { 2189566063dSJacob Faibussowitsch PetscCall(MatMultTransposeEqual(B,D,10,&flg)); 219c4fb69feSStefano Zampini if (!flg) { 2209566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMultTranspose error after MatDuplicate\n")); 221c4fb69feSStefano Zampini } 222c4fb69feSStefano Zampini } 2239566063dSJacob Faibussowitsch PetscCall(MatDestroy(&D)); 224c4fb69feSStefano Zampini 225c4fb69feSStefano Zampini if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */ 2269566063dSJacob Faibussowitsch PetscCall(VecSetRandom(Ay,NULL)); 2279566063dSJacob Faibussowitsch PetscCall(VecCopy(Ay,By)); 2289566063dSJacob Faibussowitsch PetscCall(MatMultTranspose(A,Ay,Ax)); 2299566063dSJacob Faibussowitsch PetscCall(MatMultTranspose(B,By,Bx)); 2309566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(Ax,NULL,"-multtrans_vec_view")); 2319566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(Bx,NULL,"-multtrans_vec_view")); 2329566063dSJacob Faibussowitsch PetscCall(VecNorm(Ax,NORM_INFINITY,&nX)); 2339566063dSJacob Faibussowitsch PetscCall(VecAXPY(Ax,-1.0,Bx)); 2349566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(Ax,NULL,"-multtrans_vec_view")); 2359566063dSJacob Faibussowitsch PetscCall(VecNorm(Ax,NORM_INFINITY,&err)); 2369566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMultTranspose err %g\n",err/nX)); 2379566063dSJacob Faibussowitsch PetscCall(VecScale(Bx,-1.0)); 2389566063dSJacob Faibussowitsch PetscCall(MatMultTransposeAdd(B,By,Bx,Bx)); 2399566063dSJacob Faibussowitsch PetscCall(VecNorm(Bx,NORM_INFINITY,&err)); 2404c3a5f3cSStefano Zampini if (err > 10.*PETSC_SMALL) { 2419566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMultTransposeAdd err %g\n",err)); 242c4fb69feSStefano Zampini } 243c4fb69feSStefano Zampini } 2449566063dSJacob Faibussowitsch PetscCall(VecDestroy(&Ax)); 2459566063dSJacob Faibussowitsch PetscCall(VecDestroy(&Ay)); 2469566063dSJacob Faibussowitsch PetscCall(VecDestroy(&Bx)); 2479566063dSJacob Faibussowitsch PetscCall(VecDestroy(&By)); 248c4fb69feSStefano Zampini 249c4fb69feSStefano Zampini /* Test MatMatMult */ 250c4fb69feSStefano Zampini if (ldc) { 2519566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(nrhs*(n+ldc),&Cdata)); 252c4fb69feSStefano Zampini } 2539566063dSJacob Faibussowitsch PetscCall(MatCreateDense(PETSC_COMM_WORLD,n,PETSC_DECIDE,N,nrhs,Cdata,&C)); 2549566063dSJacob Faibussowitsch PetscCall(MatDenseSetLDA(C,n+ldc)); 2559566063dSJacob Faibussowitsch PetscCall(MatSetRandom(C,NULL)); 256c4fb69feSStefano Zampini if (cgpu) { 2579566063dSJacob Faibussowitsch PetscCall(MatConvert(C,MATDENSECUDA,MAT_INPLACE_MATRIX,&C)); 258c4fb69feSStefano Zampini } 259c4fb69feSStefano Zampini for (nt = 0; nt < ntrials; nt++) { 2609566063dSJacob Faibussowitsch PetscCall(MatMatMult(B,C,nt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D)); 2619566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(D,NULL,"-bc_view")); 2629566063dSJacob Faibussowitsch PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)D,&flg,MATSEQDENSE,MATMPIDENSE,"")); 263c4fb69feSStefano Zampini if (flg) { 2649566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(B,&x,&y)); 2659566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(D,NULL,&v)); 266c4fb69feSStefano Zampini for (i = 0; i < nrhs; i++) { 2679566063dSJacob Faibussowitsch PetscCall(MatGetColumnVector(D,v,i)); 2689566063dSJacob Faibussowitsch PetscCall(MatGetColumnVector(C,x,i)); 2699566063dSJacob Faibussowitsch PetscCall(MatMult(B,x,y)); 2709566063dSJacob Faibussowitsch PetscCall(VecAXPY(y,-1.0,v)); 2719566063dSJacob Faibussowitsch PetscCall(VecNorm(y,NORM_INFINITY,&err)); 2729566063dSJacob Faibussowitsch if (err > 10.*PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMat err %" PetscInt_FMT " %g\n",i,err)); 273c4fb69feSStefano Zampini } 2749566063dSJacob Faibussowitsch PetscCall(VecDestroy(&y)); 2759566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x)); 2769566063dSJacob Faibussowitsch PetscCall(VecDestroy(&v)); 277c4fb69feSStefano Zampini } 278c4fb69feSStefano Zampini } 2799566063dSJacob Faibussowitsch PetscCall(MatDestroy(&D)); 280c4fb69feSStefano Zampini 281c4fb69feSStefano Zampini /* Test MatTransposeMatMult */ 282c4fb69feSStefano Zampini if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */ 283c4fb69feSStefano Zampini for (nt = 0; nt < ntrials; nt++) { 2849566063dSJacob Faibussowitsch PetscCall(MatTransposeMatMult(B,C,nt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D)); 2859566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(D,NULL,"-btc_view")); 2869566063dSJacob Faibussowitsch PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)D,&flg,MATSEQDENSE,MATMPIDENSE,"")); 287c4fb69feSStefano Zampini if (flg) { 2889566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(B,&y,&x)); 2899566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(D,NULL,&v)); 290c4fb69feSStefano Zampini for (i = 0; i < nrhs; i++) { 2919566063dSJacob Faibussowitsch PetscCall(MatGetColumnVector(D,v,i)); 2929566063dSJacob Faibussowitsch PetscCall(MatGetColumnVector(C,x,i)); 2939566063dSJacob Faibussowitsch PetscCall(MatMultTranspose(B,x,y)); 2949566063dSJacob Faibussowitsch PetscCall(VecAXPY(y,-1.0,v)); 2959566063dSJacob Faibussowitsch PetscCall(VecNorm(y,NORM_INFINITY,&err)); 2969566063dSJacob Faibussowitsch if (err > 10.*PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatTransMat err %" PetscInt_FMT " %g\n",i,err)); 297c4fb69feSStefano Zampini } 2989566063dSJacob Faibussowitsch PetscCall(VecDestroy(&y)); 2999566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x)); 3009566063dSJacob Faibussowitsch PetscCall(VecDestroy(&v)); 301c4fb69feSStefano Zampini } 302c4fb69feSStefano Zampini } 3039566063dSJacob Faibussowitsch PetscCall(MatDestroy(&D)); 304c4fb69feSStefano Zampini } 305c4fb69feSStefano Zampini 306300d917bSStefano Zampini /* Test basis orthogonalization */ 30753022affSStefano Zampini if (testorthog) { 3089566063dSJacob Faibussowitsch PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&D)); 3099566063dSJacob Faibussowitsch PetscCall(MatSetOption(D,MAT_SYMMETRIC,symm)); 3109566063dSJacob Faibussowitsch PetscCall(MatH2OpusOrthogonalize(D)); 3119566063dSJacob Faibussowitsch PetscCall(MatMultEqual(B,D,10,&flg)); 31253022affSStefano Zampini if (!flg) { 3139566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMult error after basis ortogonalization\n")); 31453022affSStefano Zampini } 3159566063dSJacob Faibussowitsch PetscCall(MatDestroy(&D)); 31653022affSStefano Zampini } 31753022affSStefano Zampini 318300d917bSStefano Zampini /* Test matrix compression */ 31953022affSStefano Zampini if (testcompress) { 3209566063dSJacob Faibussowitsch PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&D)); 3219566063dSJacob Faibussowitsch PetscCall(MatSetOption(D,MAT_SYMMETRIC,symm)); 3229566063dSJacob Faibussowitsch PetscCall(MatH2OpusCompress(D,PETSC_SMALL)); 3239566063dSJacob Faibussowitsch PetscCall(MatDestroy(&D)); 32453022affSStefano Zampini } 32553022affSStefano Zampini 326300d917bSStefano Zampini /* Test low-rank update */ 327300d917bSStefano Zampini if (testhlru) { 328300d917bSStefano Zampini Mat U, V; 329300d917bSStefano Zampini PetscScalar *Udata = NULL, *Vdata = NULL; 330300d917bSStefano Zampini 331300d917bSStefano Zampini if (ldu) { 3329566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(nlr*(n+ldu),&Udata)); 3339566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(nlr*(n+ldu+2),&Vdata)); 334300d917bSStefano Zampini } 3359566063dSJacob Faibussowitsch PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&D)); 3369566063dSJacob Faibussowitsch PetscCall(MatCreateDense(PetscObjectComm((PetscObject)D),n,PETSC_DECIDE,N,nlr,Udata,&U)); 3379566063dSJacob Faibussowitsch PetscCall(MatDenseSetLDA(U,n+ldu)); 3389566063dSJacob Faibussowitsch PetscCall(MatCreateDense(PetscObjectComm((PetscObject)D),n,PETSC_DECIDE,N,nlr,Vdata,&V)); 3399566063dSJacob Faibussowitsch if (ldu) PetscCall(MatDenseSetLDA(V,n+ldu+2)); 3409566063dSJacob Faibussowitsch PetscCall(MatSetRandom(U,NULL)); 3419566063dSJacob Faibussowitsch PetscCall(MatSetRandom(V,NULL)); 3429566063dSJacob Faibussowitsch PetscCall(MatH2OpusLowRankUpdate(D,U,V,0.5)); 3439566063dSJacob Faibussowitsch PetscCall(MatH2OpusLowRankUpdate(D,U,V,-0.5)); 3449566063dSJacob Faibussowitsch PetscCall(MatMultEqual(B,D,10,&flg)); 345300d917bSStefano Zampini if (!flg) { 3469566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMult error after low-rank update\n")); 347300d917bSStefano Zampini } 3489566063dSJacob Faibussowitsch PetscCall(MatDestroy(&D)); 3499566063dSJacob Faibussowitsch PetscCall(MatDestroy(&U)); 3509566063dSJacob Faibussowitsch PetscCall(PetscFree(Udata)); 3519566063dSJacob Faibussowitsch PetscCall(MatDestroy(&V)); 3529566063dSJacob Faibussowitsch PetscCall(PetscFree(Vdata)); 353300d917bSStefano Zampini } 354300d917bSStefano Zampini 355c4fb69feSStefano Zampini /* check explicit operator */ 356c4fb69feSStefano Zampini if (checkexpl) { 357c4fb69feSStefano Zampini Mat Be, Bet; 358c4fb69feSStefano Zampini 3599566063dSJacob Faibussowitsch PetscCall(MatComputeOperator(B,MATDENSE,&D)); 3609566063dSJacob Faibussowitsch PetscCall(MatDuplicate(D,MAT_COPY_VALUES,&Be)); 3619566063dSJacob Faibussowitsch PetscCall(MatNorm(D,NORM_FROBENIUS,&nB)); 3629566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(D,NULL,"-expl_view")); 3639566063dSJacob Faibussowitsch PetscCall(MatAXPY(D,-1.0,A,SAME_NONZERO_PATTERN)); 3649566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(D,NULL,"-diff_view")); 3659566063dSJacob Faibussowitsch PetscCall(MatNorm(D,NORM_FROBENIUS,&nD)); 3669566063dSJacob Faibussowitsch PetscCall(MatNorm(A,NORM_FROBENIUS,&nA)); 3679566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Approximation error %g (%g / %g, %g)\n",nD/nA,nD,nA,nB)); 3689566063dSJacob Faibussowitsch PetscCall(MatDestroy(&D)); 369c4fb69feSStefano Zampini 370c4fb69feSStefano Zampini if (testtrans) { /* MatMultTranspose for nonsymmetric matrices not implemented */ 3719566063dSJacob Faibussowitsch PetscCall(MatTranspose(A,MAT_INPLACE_MATRIX,&A)); 3729566063dSJacob Faibussowitsch PetscCall(MatComputeOperatorTranspose(B,MATDENSE,&D)); 3739566063dSJacob Faibussowitsch PetscCall(MatDuplicate(D,MAT_COPY_VALUES,&Bet)); 3749566063dSJacob Faibussowitsch PetscCall(MatNorm(D,NORM_FROBENIUS,&nB)); 3759566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(D,NULL,"-expl_trans_view")); 3769566063dSJacob Faibussowitsch PetscCall(MatAXPY(D,-1.0,A,SAME_NONZERO_PATTERN)); 3779566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(D,NULL,"-diff_trans_view")); 3789566063dSJacob Faibussowitsch PetscCall(MatNorm(D,NORM_FROBENIUS,&nD)); 3799566063dSJacob Faibussowitsch PetscCall(MatNorm(A,NORM_FROBENIUS,&nA)); 3809566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Approximation error transpose %g (%g / %g, %g)\n",nD/nA,nD,nA,nB)); 3819566063dSJacob Faibussowitsch PetscCall(MatDestroy(&D)); 382c4fb69feSStefano Zampini 3839566063dSJacob Faibussowitsch PetscCall(MatTranspose(Bet,MAT_INPLACE_MATRIX,&Bet)); 3849566063dSJacob Faibussowitsch PetscCall(MatAXPY(Be,-1.0,Bet,SAME_NONZERO_PATTERN)); 3859566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(Be,NULL,"-diff_expl_view")); 3869566063dSJacob Faibussowitsch PetscCall(MatNorm(Be,NORM_FROBENIUS,&nB)); 3879566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Approximation error B - (B^T)^T %g\n",nB)); 3889566063dSJacob Faibussowitsch PetscCall(MatDestroy(&Be)); 3899566063dSJacob Faibussowitsch PetscCall(MatDestroy(&Bet)); 390c4fb69feSStefano Zampini } 391c4fb69feSStefano Zampini } 3929566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A)); 3939566063dSJacob Faibussowitsch PetscCall(MatDestroy(&B)); 3949566063dSJacob Faibussowitsch PetscCall(MatDestroy(&C)); 3959566063dSJacob Faibussowitsch PetscCall(PetscFree(Cdata)); 3969566063dSJacob Faibussowitsch PetscCall(PetscFree(Adata)); 3979566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 398b122ec5aSJacob Faibussowitsch return 0; 399c4fb69feSStefano Zampini } 400c4fb69feSStefano Zampini 401c4fb69feSStefano Zampini /*TEST 402c4fb69feSStefano Zampini 403c4fb69feSStefano Zampini build: 40453022affSStefano Zampini requires: h2opus 405c4fb69feSStefano Zampini 406c4fb69feSStefano Zampini #tests from kernel 407c4fb69feSStefano Zampini test: 40853022affSStefano Zampini requires: h2opus 409c4fb69feSStefano Zampini nsize: 1 410c4fb69feSStefano Zampini suffix: 1 411c4fb69feSStefano Zampini args: -n {{17 33}} -kernel 1 -dim {{1 2 3}} -symm {{0 1}} -checkexpl -bgpu 0 412c4fb69feSStefano Zampini 413c4fb69feSStefano Zampini test: 41453022affSStefano Zampini requires: h2opus 415c4fb69feSStefano Zampini nsize: 1 416c4fb69feSStefano Zampini suffix: 1_ld 417c4fb69feSStefano Zampini output_file: output/ex66_1.out 418300d917bSStefano Zampini args: -n 33 -kernel 1 -dim 1 -lda 13 -ldc 11 -ldu 17 -symm 0 -checkexpl -bgpu 0 419c4fb69feSStefano Zampini 420c4fb69feSStefano Zampini test: 42153022affSStefano Zampini requires: h2opus cuda 422c4fb69feSStefano Zampini nsize: 1 423c4fb69feSStefano Zampini suffix: 1_cuda 424c4fb69feSStefano Zampini output_file: output/ex66_1.out 425c4fb69feSStefano Zampini args: -n {{17 33}} -kernel 1 -dim {{1 2 3}} -symm {{0 1}} -checkexpl -bgpu 1 426c4fb69feSStefano Zampini 427c4fb69feSStefano Zampini test: 42853022affSStefano Zampini requires: h2opus cuda 429c4fb69feSStefano Zampini nsize: 1 430c4fb69feSStefano Zampini suffix: 1_cuda_ld 431c4fb69feSStefano Zampini output_file: output/ex66_1.out 432300d917bSStefano Zampini args: -n 33 -kernel 1 -dim 1 -lda 13 -ldc 11 -ldu 17 -symm 0 -checkexpl -bgpu 1 433c4fb69feSStefano Zampini 434c4fb69feSStefano Zampini test: 43553022affSStefano Zampini requires: h2opus 43653022affSStefano Zampini nsize: {{2 3}} 437c4fb69feSStefano Zampini suffix: 1_par 43853022affSStefano Zampini args: -n 64 -symm -kernel 1 -dim 1 -ldc 12 -testlayout {{0 1}} -bgpu 0 -cgpu 0 439c4fb69feSStefano Zampini 440c4fb69feSStefano Zampini test: 44153022affSStefano Zampini requires: h2opus cuda 44253022affSStefano Zampini nsize: {{2 3}} 443c4fb69feSStefano Zampini suffix: 1_par_cuda 44453022affSStefano Zampini args: -n 64 -symm -kernel 1 -dim 1 -ldc 12 -testlayout {{0 1}} -bgpu {{0 1}} -cgpu {{0 1}} 445c4fb69feSStefano Zampini output_file: output/ex66_1_par.out 446c4fb69feSStefano Zampini 447c4fb69feSStefano Zampini #tests from matrix sampling (parallel or unsymmetric not supported) 448c4fb69feSStefano Zampini test: 44953022affSStefano Zampini requires: h2opus 450c4fb69feSStefano Zampini nsize: 1 451c4fb69feSStefano Zampini suffix: 2 452c4fb69feSStefano Zampini args: -n {{17 33}} -kernel 0 -dim 2 -symm 1 -checkexpl -bgpu 0 453c4fb69feSStefano Zampini 454c4fb69feSStefano Zampini test: 45553022affSStefano Zampini requires: h2opus cuda 456c4fb69feSStefano Zampini nsize: 1 457c4fb69feSStefano Zampini suffix: 2_cuda 458c4fb69feSStefano Zampini output_file: output/ex66_2.out 459c4fb69feSStefano Zampini args: -n {{17 33}} -kernel 0 -dim 2 -symm 1 -checkexpl -bgpu {{0 1}} -agpu {{0 1}} 460c4fb69feSStefano Zampini 46153022affSStefano Zampini #tests view operation 46253022affSStefano Zampini test: 46353022affSStefano Zampini requires: h2opus !cuda 4648cc725e6SPierre Jolivet filter: grep -v " MPI process" | grep -v "\[" | grep -v "\]" 46553022affSStefano Zampini nsize: {{1 2 3}} 46653022affSStefano Zampini suffix: view 46753022affSStefano Zampini 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 46853022affSStefano Zampini 46953022affSStefano Zampini test: 47053022affSStefano Zampini requires: h2opus cuda 4718cc725e6SPierre Jolivet filter: grep -v " MPI process" | grep -v "\[" | grep -v "\]" 47253022affSStefano Zampini nsize: {{1 2 3}} 47353022affSStefano Zampini suffix: view_cuda 47453022affSStefano Zampini 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 47553022affSStefano Zampini 476c4fb69feSStefano Zampini TEST*/ 477