Actual source code: ex43-aijcusparse.c

petsc-3.4.2 2013-07-02
  1: static char help[] = "Reads a PETSc matrix from a file and solves a linear system \n\
  2: using the aijcusparse class. This example also demonstrates how to set the storage \n\
  3: format on the GPU using the MatCUSPARSESetFormat method. Input parameters are:\n\
  4:   -f <input_file> : the file to load\n\n";

  6: /*
  7:   This code can be used to test PETSc interface to other packages.\n\
  8:   Examples of command line options:       \n\
  9:    ./ex43-aijcusparse -f DATAFILESPATH/matrices/cfd.2.10 -mat_cusparse_mult_storage_format ell  \n\
 10:   In a second example, one can read a symmetric matrix stored in upper triangular form.\n\
 11:   Then one can invoke the ICC preconditioner, however one has to indicate explicitly \n\
 12:   that the matrix is symmetric.
 13:    ./ex43-aijcusparse -f DATAFILESPATH/matrices/shallow_water1 -ksp_type cgs -pc_type icc -mat_symmetric -mat_cusparse_mult_storage_format ell  \n\
 14:    \n\n";
 15: */

 17: #include <petscksp.h>

 21: int main(int argc,char **argv)
 22: {
 23:   KSP                ksp;
 24:   Mat                A;
 25:   Vec                X,B;
 26:   PetscInt           m, its;
 27:   PetscReal          norm;
 28:   char               file[PETSC_MAX_PATH_LEN];
 29:   PetscBool          flg;
 30:   PetscViewer        fd;
 31:   PetscErrorCode     ierr;

 33:   PetscInitialize(&argc,&argv,0,help);

 35:   /* Load the data from a file */
 36:   PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
 37:   if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");
 38:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);

 40:   /* Build the matrix */
 41:   MatCreate(PETSC_COMM_WORLD,&A);
 42:   MatSetType(A,MATSEQAIJCUSPARSE);
 43:   MatCUSPARSESetFormat(A,MAT_CUSPARSE_MULT,MAT_CUSPARSE_ELL);

 45:   /* inform the matrix that it is symmetric */
 46:   PetscOptionsHasName(NULL, "-mat_symmetric", &flg);
 47:   if (flg) {
 48:     ierr=MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
 49:   }
 50:   MatSetFromOptions(A);
 51:   MatLoad(A,fd);

 53:   /* Build the vectors */
 54:   MatGetLocalSize(A,&m,PETSC_NULL);
 55:   VecCreateSeqCUSP(PETSC_COMM_WORLD,m,&B);
 56:   VecSet(B,1.0);
 57:   VecDuplicate(B,&X);

 59:   /* Build the KSP */
 60:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 61:   KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
 62:   KSPSetType(ksp,KSPGMRES);
 63:   KSPSetTolerances(ksp,1.0e-12,PETSC_DEFAULT,PETSC_DEFAULT,100);
 64:   KSPSetFromOptions(ksp);

 66:   /* Solve */
 67:   KSPSolve(ksp,B,X);

 69:   /* print out norm and the number of iterations */
 70:   KSPGetIterationNumber(ksp,&its);
 71:   KSPGetResidualNorm(ksp,&norm);
 72:   PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
 73:   PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);

 75:   /* Cleanup */
 76:   VecDestroy(&X);
 77:   VecDestroy(&B);
 78:   MatDestroy(&A);
 79:   KSPDestroy(&ksp);
 80:   PetscViewerDestroy(&fd);
 81:   PetscFinalize();
 82:   return 0;
 83: }