Actual source code: ex6.c
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 a 5X5 example of the 5-pt. stencil,\n\
5: use the file petsc/src/mat/examples/matbinary.ex\n\n";
7: #include <petscksp.h>
8: #include <petsclog.h>
12: int main(int argc,char **args)
13: {
14: #if !defined(PETSC_USE_COMPLEX)
16: PetscInt its;
17: PetscLogStage stage1,stage2;
18: PetscReal norm;
19: PetscLogDouble tsetup1,tsetup2,tsetup,tsolve1,tsolve2,tsolve;
20: Vec x,b,u;
21: Mat A;
22: char file[PETSC_MAX_PATH_LEN];
23: PetscViewer fd;
24: PetscBool table = PETSC_FALSE,flg;
25: KSP ksp;
26: #endif
28: PetscInitialize(&argc,&args,(char *)0,help);
30: #if defined(PETSC_USE_COMPLEX)
31: SETERRQ(PETSC_COMM_WORLD,1,"This example does not work with complex numbers");
32: #else
33: PetscOptionsGetBool(PETSC_NULL,"-table",&table,PETSC_NULL);
36: /* Read matrix and RHS */
37: PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
38: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");
39: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
40: MatCreate(PETSC_COMM_WORLD,&A);
41: MatLoad(A,fd);
42: VecCreate(PETSC_COMM_WORLD,&b);
43: VecLoad(b,fd);
44: PetscViewerDestroy(&fd);
46: /*
47: If the load matrix is larger then the vector, due to being padded
48: to match the blocksize then create a new padded vector
49: */
50: {
51: PetscInt m,n,j,mvec,start,end,indx;
52: Vec tmp;
53: PetscScalar *bold;
55: MatGetLocalSize(A,&m,&n);
56: VecCreate(PETSC_COMM_WORLD,&tmp);
57: VecSetSizes(tmp,m,PETSC_DECIDE);
58: VecSetFromOptions(tmp);
59: VecGetOwnershipRange(b,&start,&end);
60: VecGetLocalSize(b,&mvec);
61: VecGetArray(b,&bold);
62: for (j=0; j<mvec; j++) {
63: indx = start+j;
64: VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
65: }
66: VecRestoreArray(b,&bold);
67: VecDestroy(&b);
68: VecAssemblyBegin(tmp);
69: VecAssemblyEnd(tmp);
70: b = tmp;
71: }
72: VecDuplicate(b,&x);
73: VecDuplicate(b,&u);
75: VecSet(x,0.0);
76: PetscBarrier((PetscObject)A);
78: PetscLogStageRegister("mystage 1",&stage1);
79: PetscLogStagePush(stage1);
80: PetscGetTime(&tsetup1);
81: KSPCreate(PETSC_COMM_WORLD,&ksp);
82: KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
83: KSPSetFromOptions(ksp);
84: KSPSetUp(ksp);
85: KSPSetUpOnBlocks(ksp);
86: PetscGetTime(&tsetup2);
87: tsetup = tsetup2 -tsetup1;
88: PetscLogStagePop();
89: PetscBarrier((PetscObject)A);
91: PetscLogStageRegister("mystage 2",&stage2);
92: PetscLogStagePush(stage2);
93: PetscGetTime(&tsolve1);
94: KSPSolve(ksp,b,x);
95: PetscGetTime(&tsolve2);
96: tsolve = tsolve2 - tsolve1;
97: PetscLogStagePop();
99: /* Show result */
100: MatMult(A,x,u);
101: VecAXPY(u,-1.0,b);
102: VecNorm(u,NORM_2,&norm);
103: KSPGetIterationNumber(ksp,&its);
104: /* matrix PC KSP Options its residual setuptime solvetime */
105: if (table) {
106: char *matrixname,kspinfo[120];
107: PetscViewer viewer;
108: PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,120,&viewer);
109: KSPView(ksp,viewer);
110: PetscStrrchr(file,'/',&matrixname);
111: PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3D %2.0e %2.1e %2.1e %2.1e %s \n",
112: matrixname,its,norm,tsetup+tsolve,tsetup,tsolve,kspinfo);
113: PetscViewerDestroy(&viewer);
114: } else {
115: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
116: PetscPrintf(PETSC_COMM_WORLD,"Residual norm = %A\n",norm);
117: }
119: /* Cleanup */
120: KSPDestroy(&ksp);
121: VecDestroy(&x);
122: VecDestroy(&b);
123: VecDestroy(&u);
124: MatDestroy(&A);
126: PetscFinalize();
127: #endif
128: return 0;
129: }