Actual source code: ex168.c
2: static char help[] = "Tests external Clique direct solvers. Simplified from ex130.c\n\
3: Example: mpiexec -n <np> ./ex168 -f <matrix binary file> \n\n";
5: #include <petscmat.h>
7: int main(int argc, char **args)
8: {
9: Mat A, F;
10: Vec u, x, b;
11: PetscMPIInt rank, size;
12: PetscInt m, n, nfact;
13: PetscReal norm, tol = 1.e-12, Anorm;
14: IS perm, iperm;
15: MatFactorInfo info;
16: PetscBool flg, testMatSolve = PETSC_TRUE;
17: PetscViewer fd; /* viewer */
18: char file[PETSC_MAX_PATH_LEN]; /* input file name */
20: PetscFunctionBeginUser;
21: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
22: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
23: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
25: /* Determine file from which we read the matrix A */
26: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
27: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
29: /* Load matrix A */
30: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
31: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
32: PetscCall(MatLoad(A, fd));
33: PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
34: PetscCall(VecLoad(b, fd));
35: PetscCall(PetscViewerDestroy(&fd));
36: PetscCall(MatGetLocalSize(A, &m, &n));
37: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
38: PetscCall(MatNorm(A, NORM_INFINITY, &Anorm));
40: /* Create vectors */
41: PetscCall(VecDuplicate(b, &x));
42: PetscCall(VecDuplicate(x, &u)); /* save the true solution */
44: /* Test Cholesky Factorization */
45: PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &perm, &iperm));
47: if (rank == 0) printf(" Clique Cholesky:\n");
48: PetscCall(MatGetFactor(A, MATSOLVERCLIQUE, MAT_FACTOR_CHOLESKY, &F));
50: info.fill = 5.0;
51: PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
53: for (nfact = 0; nfact < 1; nfact++) {
54: if (rank == 0) printf(" %d-the Cholesky numfactorization \n", nfact);
55: PetscCall(MatCholeskyFactorNumeric(F, A, &info));
57: /* Test MatSolve() */
58: if (testMatSolve && nfact == 2) {
59: PetscCall(MatSolve(F, b, x));
61: /* Check the residual */
62: PetscCall(MatMult(A, x, u));
63: PetscCall(VecAXPY(u, -1.0, b));
64: PetscCall(VecNorm(u, NORM_INFINITY, &norm));
65: /* if (norm > tol) { */
66: if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSolve: rel residual %g/%g = %g, LU numfact %d\n", norm, Anorm, norm / Anorm, nfact));
67: /*} */
68: }
69: }
71: /* Free data structures */
72: PetscCall(MatDestroy(&A));
73: PetscCall(MatDestroy(&F));
74: PetscCall(ISDestroy(&perm));
75: PetscCall(ISDestroy(&iperm));
76: PetscCall(VecDestroy(&x));
77: PetscCall(VecDestroy(&b));
78: PetscCall(VecDestroy(&u));
79: PetscCall(PetscFinalize());
80: return 0;
81: }