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