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: }