Actual source code: ex268.c

  1: static char help[] = "Tests MATFACTORHTOOL\n\n";

  3: #include <petscmat.h>

  5: static PetscErrorCode GenEntries(PetscInt sdim, PetscInt M, PetscInt N, const PetscInt *J, const PetscInt *K, PetscScalar *ptr, void *ctx)
  6: {
  7:   PetscInt  d, j, k;
  8:   PetscReal diff = 0.0, *coords = (PetscReal *)(ctx);

 10:   PetscFunctionBeginUser;
 11:   for (j = 0; j < M; j++) {
 12:     for (k = 0; k < N; k++) {
 13:       diff = 0.0;
 14:       for (d = 0; d < sdim; d++) diff += (coords[J[j] * sdim + d] - coords[K[k] * sdim + d]) * (coords[J[j] * sdim + d] - coords[K[k] * sdim + d]);
 15:       ptr[j + M * k] = 1.0 / (1.0e-1 + PetscSqrtReal(diff));
 16:     }
 17:   }
 18:   PetscFunctionReturn(PETSC_SUCCESS);
 19: }

 21: int main(int argc, char **argv)
 22: {
 23:   Mat               A, Ad, F, Fd, X, Xd, B;
 24:   Vec               x, xd, b;
 25:   PetscInt          m = 100, dim = 3, M, K = 10, begin, n = 0;
 26:   PetscMPIInt       size;
 27:   PetscReal        *coords, *gcoords, norm, epsilon;
 28:   MatHtoolKernelFn *kernel = GenEntries;
 29:   PetscBool         flg, sym = PETSC_FALSE;
 30:   PetscRandom       rdm;
 31:   MatSolverType     type;

 33:   PetscFunctionBeginUser;
 34:   PetscCall(PetscInitialize(&argc, &argv, (char *)NULL, help));
 35:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m_local", &m, NULL));
 36:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n_local", &n, NULL));
 37:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL));
 38:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-K", &K, NULL));
 39:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-symmetric", &sym, NULL));
 40:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-mat_htool_epsilon", &epsilon, NULL));
 41:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 42:   M = size * m;
 43:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
 44:   PetscCall(PetscMalloc1(m * dim, &coords));
 45:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
 46:   PetscCall(PetscRandomGetValuesReal(rdm, m * dim, coords));
 47:   PetscCall(PetscCalloc1(M * dim, &gcoords));
 48:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, M, K, NULL, &B));
 49:   PetscCall(MatSetRandom(B, rdm));
 50:   PetscCall(MatGetOwnershipRange(B, &begin, NULL));
 51:   PetscCall(PetscArraycpy(gcoords + begin * dim, coords, m * dim));
 52:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, gcoords, M * dim, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
 53:   PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, m, M, M, dim, coords, coords, kernel, gcoords, &A));
 54:   PetscCall(MatSetOption(A, MAT_SYMMETRIC, sym));
 55:   PetscCall(MatSetFromOptions(A));
 56:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 57:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 58:   PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &Ad));
 59:   PetscCall(MatMultEqual(A, Ad, 10, &flg));
 60:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Ax != Adx");
 61:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, M, K, NULL, &X));
 62:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, M, K, NULL, &Xd));
 63:   PetscCall(MatViewFromOptions(A, NULL, "-A"));
 64:   PetscCall(MatViewFromOptions(Ad, NULL, "-Ad"));
 65:   PetscCall(MatViewFromOptions(B, NULL, "-B"));
 66:   for (PetscInt i = 0; i < 2; ++i) {
 67:     PetscCall(MatGetFactor(A, MATSOLVERHTOOL, i == 0 ? MAT_FACTOR_LU : MAT_FACTOR_CHOLESKY, &F));
 68:     PetscCall(MatGetFactor(Ad, MATSOLVERPETSC, i == 0 ? MAT_FACTOR_LU : MAT_FACTOR_CHOLESKY, &Fd));
 69:     PetscCall(MatFactorGetSolverType(F, &type));
 70:     PetscCall(PetscStrncmp(type, MATSOLVERHTOOL, 5, &flg));
 71:     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "MATSOLVERHTOOL != htool");
 72:     if (i == 0) {
 73:       PetscCall(MatLUFactorSymbolic(F, A, NULL, NULL, NULL));
 74:       PetscCall(MatLUFactorNumeric(F, A, NULL));
 75:       PetscCall(MatLUFactorSymbolic(Fd, Ad, NULL, NULL, NULL));
 76:       PetscCall(MatLUFactorNumeric(Fd, Ad, NULL));
 77:     } else {
 78:       PetscCall(MatCholeskyFactorSymbolic(F, A, NULL, NULL));
 79:       PetscCall(MatCholeskyFactorNumeric(F, A, NULL));
 80:       PetscCall(MatCholeskyFactorSymbolic(Fd, Ad, NULL, NULL));
 81:       PetscCall(MatCholeskyFactorNumeric(Fd, Ad, NULL));
 82:     }
 83:     PetscCall(MatMatSolve(F, B, X));
 84:     PetscCall(MatMatSolve(Fd, B, Xd));
 85:     PetscCall(MatViewFromOptions(X, NULL, "-X"));
 86:     PetscCall(MatViewFromOptions(Xd, NULL, "-Xd"));
 87:     PetscCall(MatAXPY(Xd, -1.0, X, SAME_NONZERO_PATTERN));
 88:     PetscCall(MatNorm(Xd, NORM_INFINITY, &norm));
 89:     PetscCall(MatViewFromOptions(Xd, NULL, "-MatMatSolve"));
 90:     if (norm > 0.01) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: norm of residual for MatMatSolve %g\n", (double)norm));
 91:     if (!PetscDefined(USE_COMPLEX) || i == 0) {
 92:       PetscCall(MatMatSolveTranspose(F, B, X));
 93:       PetscCall(MatMatSolveTranspose(Fd, B, Xd));
 94:       PetscCall(MatAXPY(Xd, -1.0, X, SAME_NONZERO_PATTERN));
 95:       PetscCall(MatNorm(Xd, NORM_INFINITY, &norm));
 96:       PetscCall(MatViewFromOptions(Xd, NULL, "-MatMatSolveTranspose"));
 97:       if (norm > 0.01) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: norm of residual for MatMatSolveTranspose %g\n", (double)norm));
 98:     }
 99:     PetscCall(MatDenseGetColumnVecRead(B, 0, &b));
100:     PetscCall(MatDenseGetColumnVecWrite(X, 0, &x));
101:     PetscCall(MatDenseGetColumnVecWrite(Xd, 0, &xd));
102:     PetscCall(MatSolve(F, b, x));
103:     PetscCall(MatSolve(Fd, b, xd));
104:     PetscCall(VecAXPY(xd, -1.0, x));
105:     PetscCall(VecNorm(xd, NORM_INFINITY, &norm));
106:     PetscCall(MatViewFromOptions(Xd, NULL, "-MatSolve"));
107:     if (norm > 0.01) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: norm of residual for MatSolve %g\n", (double)norm));
108:     if (!PetscDefined(USE_COMPLEX) || i == 0) {
109:       PetscCall(MatSolveTranspose(F, b, x));
110:       PetscCall(MatSolveTranspose(Fd, b, xd));
111:       PetscCall(VecAXPY(xd, -1.0, x));
112:       PetscCall(VecNorm(xd, NORM_INFINITY, &norm));
113:       PetscCall(MatViewFromOptions(Xd, NULL, "-MatSolveTranspose"));
114:       if (norm > 0.01) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: norm of residual for MatSolveTranspose %g\n", (double)norm));
115:     }
116:     PetscCall(MatDenseRestoreColumnVecWrite(Xd, 0, &xd));
117:     PetscCall(MatDenseRestoreColumnVecWrite(X, 0, &x));
118:     PetscCall(MatDenseRestoreColumnVecRead(B, 0, &b));
119:     PetscCall(MatDestroy(&Fd));
120:     PetscCall(MatDestroy(&F));
121:   }
122:   PetscCall(MatDestroy(&Xd));
123:   PetscCall(MatDestroy(&X));
124:   PetscCall(PetscRandomDestroy(&rdm));
125:   PetscCall(MatDestroy(&Ad));
126:   PetscCall(MatDestroy(&A));
127:   PetscCall(MatDestroy(&B));
128:   PetscCall(PetscFree(gcoords));
129:   PetscCall(PetscFree(coords));
130:   PetscCall(PetscFinalize());
131:   return 0;
132: }

134: /*TEST

136:    build:
137:       requires: htool

139:    test:
140:       requires: htool
141:       suffix: 1
142:       nsize: 1
143:       args: -mat_htool_epsilon 1.0e-11
144:       output_file: output/empty.out

146: TEST*/