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