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