1c4762a1bSJed Brown static char help[] = "Tests external Clique direct solvers. Simplified from ex130.c\n\
2c4762a1bSJed Brown Example: mpiexec -n <np> ./ex168 -f <matrix binary file> \n\n";
3c4762a1bSJed Brown
4c4762a1bSJed Brown #include <petscmat.h>
5c4762a1bSJed Brown
main(int argc,char ** args)6d71ae5a4SJacob Faibussowitsch int main(int argc, char **args)
7d71ae5a4SJacob Faibussowitsch {
8c4762a1bSJed Brown Mat A, F;
9c4762a1bSJed Brown Vec u, x, b;
10c4762a1bSJed Brown PetscMPIInt rank, size;
11c4762a1bSJed Brown PetscInt m, n, nfact;
12c4762a1bSJed Brown PetscReal norm, tol = 1.e-12, Anorm;
13c4762a1bSJed Brown IS perm, iperm;
14c4762a1bSJed Brown MatFactorInfo info;
15c4762a1bSJed Brown PetscBool flg, testMatSolve = PETSC_TRUE;
16c4762a1bSJed Brown PetscViewer fd; /* viewer */
17c4762a1bSJed Brown char file[PETSC_MAX_PATH_LEN]; /* input file name */
18c4762a1bSJed Brown
19327415f7SBarry Smith PetscFunctionBeginUser;
20*c8025a54SPierre Jolivet PetscCall(PetscInitialize(&argc, &args, NULL, help));
219566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
229566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
23c4762a1bSJed Brown
24c4762a1bSJed Brown /* Determine file from which we read the matrix A */
259566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
2628b400f6SJacob Faibussowitsch PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
27c4762a1bSJed Brown
28c4762a1bSJed Brown /* Load matrix A */
299566063dSJacob Faibussowitsch PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
309566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
319566063dSJacob Faibussowitsch PetscCall(MatLoad(A, fd));
329566063dSJacob Faibussowitsch PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
339566063dSJacob Faibussowitsch PetscCall(VecLoad(b, fd));
349566063dSJacob Faibussowitsch PetscCall(PetscViewerDestroy(&fd));
359566063dSJacob Faibussowitsch PetscCall(MatGetLocalSize(A, &m, &n));
3608401ef6SPierre Jolivet PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
379566063dSJacob Faibussowitsch PetscCall(MatNorm(A, NORM_INFINITY, &Anorm));
38c4762a1bSJed Brown
39c4762a1bSJed Brown /* Create vectors */
409566063dSJacob Faibussowitsch PetscCall(VecDuplicate(b, &x));
419566063dSJacob Faibussowitsch PetscCall(VecDuplicate(x, &u)); /* save the true solution */
42c4762a1bSJed Brown
43c4762a1bSJed Brown /* Test Cholesky Factorization */
449566063dSJacob Faibussowitsch PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &perm, &iperm));
45c4762a1bSJed Brown
46dd400576SPatrick Sanan if (rank == 0) printf(" Clique Cholesky:\n");
479566063dSJacob Faibussowitsch PetscCall(MatGetFactor(A, MATSOLVERCLIQUE, MAT_FACTOR_CHOLESKY, &F));
48c4762a1bSJed Brown
49c4762a1bSJed Brown info.fill = 5.0;
509566063dSJacob Faibussowitsch PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
51c4762a1bSJed Brown
52c4762a1bSJed Brown for (nfact = 0; nfact < 1; nfact++) {
53dd400576SPatrick Sanan if (rank == 0) printf(" %d-the Cholesky numfactorization \n", nfact);
549566063dSJacob Faibussowitsch PetscCall(MatCholeskyFactorNumeric(F, A, &info));
55c4762a1bSJed Brown
56c4762a1bSJed Brown /* Test MatSolve() */
57c4762a1bSJed Brown if (testMatSolve && nfact == 2) {
589566063dSJacob Faibussowitsch PetscCall(MatSolve(F, b, x));
59c4762a1bSJed Brown
60c4762a1bSJed Brown /* Check the residual */
619566063dSJacob Faibussowitsch PetscCall(MatMult(A, x, u));
629566063dSJacob Faibussowitsch PetscCall(VecAXPY(u, -1.0, b));
639566063dSJacob Faibussowitsch PetscCall(VecNorm(u, NORM_INFINITY, &norm));
64c4762a1bSJed Brown /* if (norm > tol) { */
6548a46eb9SPierre Jolivet if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSolve: rel residual %g/%g = %g, LU numfact %d\n", norm, Anorm, norm / Anorm, nfact));
66c4762a1bSJed Brown /*} */
67c4762a1bSJed Brown }
68c4762a1bSJed Brown }
69c4762a1bSJed Brown
70c4762a1bSJed Brown /* Free data structures */
719566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A));
729566063dSJacob Faibussowitsch PetscCall(MatDestroy(&F));
739566063dSJacob Faibussowitsch PetscCall(ISDestroy(&perm));
749566063dSJacob Faibussowitsch PetscCall(ISDestroy(&iperm));
759566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x));
769566063dSJacob Faibussowitsch PetscCall(VecDestroy(&b));
779566063dSJacob Faibussowitsch PetscCall(VecDestroy(&u));
789566063dSJacob Faibussowitsch PetscCall(PetscFinalize());
79b122ec5aSJacob Faibussowitsch return 0;
80c4762a1bSJed Brown }
81