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