Actual source code: ex18.c
petsc-3.5.4 2015-05-23
2: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
3: Input arguments are:\n\
4: -f <input_file> : file to load. For example see $PETSC_DIR/share/petsc/datafiles/matrices\n\n";
6: #include <petscmat.h>
7: #include <petscksp.h>
11: int main(int argc,char **args)
12: {
14: PetscInt its,m,n,mvec;
15: PetscReal norm;
16: Vec x,b,u;
17: Mat A;
18: KSP ksp;
19: char file[PETSC_MAX_PATH_LEN];
20: PetscViewer fd;
21: PetscLogStage stage1;
23: PetscInitialize(&argc,&args,(char*)0,help);
25: /* Read matrix and RHS */
26: PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,NULL);
27: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
28: MatCreate(PETSC_COMM_WORLD,&A);
29: MatSetType(A,MATSEQAIJ);
30: MatLoad(A,fd);
31: VecCreate(PETSC_COMM_WORLD,&b);
32: VecLoad(b,fd);
33: PetscViewerDestroy(&fd);
35: /*
36: If the load matrix is larger then the vector, due to being padded
37: to match the blocksize then create a new padded vector
38: */
39: MatGetSize(A,&m,&n);
40: VecGetSize(b,&mvec);
41: if (m > mvec) {
42: Vec tmp;
43: PetscScalar *bold,*bnew;
44: /* create a new vector b by padding the old one */
45: VecCreate(PETSC_COMM_WORLD,&tmp);
46: VecSetSizes(tmp,PETSC_DECIDE,m);
47: VecSetFromOptions(tmp);
48: VecGetArray(tmp,&bnew);
49: VecGetArray(b,&bold);
50: PetscMemcpy(bnew,bold,mvec*sizeof(PetscScalar));
51: VecDestroy(&b);
52: b = tmp;
53: }
55: /* Set up solution */
56: VecDuplicate(b,&x);
57: VecDuplicate(b,&u);
58: VecSet(x,0.0);
60: /* Solve system */
61: PetscLogStageRegister("Stage 1",&stage1);
62: PetscLogStagePush(stage1);
63: KSPCreate(PETSC_COMM_WORLD,&ksp);
64: KSPSetOperators(ksp,A,A);
65: KSPSetFromOptions(ksp);
66: KSPSolve(ksp,b,x);
67: PetscLogStagePop();
69: /* Show result */
70: MatMult(A,x,u);
71: VecAXPY(u,-1.0,b);
72: VecNorm(u,NORM_2,&norm);
73: KSPGetIterationNumber(ksp,&its);
74: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
75: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
77: /* Cleanup */
78: KSPDestroy(&ksp);
79: VecDestroy(&x);
80: VecDestroy(&b);
81: VecDestroy(&u);
82: MatDestroy(&A);
84: PetscFinalize();
85: return 0;
86: }