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 */
21: PetscInitialize(&argc, &args, (char *)0, help);
22: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
23: MPI_Comm_size(PETSC_COMM_WORLD, &size);
25: /* Determine file from which we read the matrix A */
26: PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg);
29: /* Load matrix A */
30: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd);
31: MatCreate(PETSC_COMM_WORLD, &A);
32: MatLoad(A, fd);
33: VecCreate(PETSC_COMM_WORLD, &b);
34: VecLoad(b, fd);
35: PetscViewerDestroy(&fd);
36: MatGetLocalSize(A, &m, &n);
38: MatNorm(A, NORM_INFINITY, &Anorm);
40: /* Create vectors */
41: VecDuplicate(b, &x);
42: VecDuplicate(x, &u); /* save the true solution */
44: /* Test Cholesky Factorization */
45: MatGetOrdering(A, MATORDERINGNATURAL, &perm, &iperm);
47: if (rank == 0) printf(" Clique Cholesky:\n");
48: MatGetFactor(A, MATSOLVERCLIQUE, MAT_FACTOR_CHOLESKY, &F);
50: info.fill = 5.0;
51: MatCholeskyFactorSymbolic(F, A, perm, &info);
53: for (nfact = 0; nfact < 1; nfact++) {
54: if (rank == 0) printf(" %d-the Cholesky numfactorization \n", nfact);
55: MatCholeskyFactorNumeric(F, A, &info);
57: /* Test MatSolve() */
58: if (testMatSolve && nfact == 2) {
59: MatSolve(F, b, x);
61: /* Check the residual */
62: MatMult(A, x, u);
63: VecAXPY(u, -1.0, b);
64: VecNorm(u, NORM_INFINITY, &norm);
65: /* if (norm > tol) { */
66: if (rank == 0) 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: MatDestroy(&A);
73: MatDestroy(&F);
74: ISDestroy(&perm);
75: ISDestroy(&iperm);
76: VecDestroy(&x);
77: VecDestroy(&b);
78: VecDestroy(&u);
79: PetscFinalize();
80: return 0;
81: }