Actual source code: ex6.c
1: /*$Id: ex6.c,v 1.74 2001/09/11 16:33:24 bsmith Exp $*/
3: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
4: Input arguments are:\n\
5: -f <input_file> : file to load. For a 5X5 example of the 5-pt. stencil,\n\
6: use the file petsc/src/mat/examples/matbinary.ex\n\n";
8: #include petscksp.h
9: #include petsclog.h
13: int main(int argc,char **args)
14: {
15: #if !defined(PETSC_USE_COMPLEX)
16: int ierr,its,stage1,stage2;
17: PetscReal norm;
18: PetscLogDouble tsetup1,tsetup2,tsetup,tsolve1,tsolve2,tsolve;
19: PetscScalar zero = 0.0,none = -1.0;
20: Vec x,b,u;
21: Mat A;
22: char file[128];
23: PetscViewer fd;
24: PetscTruth table,flg;
25: KSP ksp;
26: #endif
28: PetscInitialize(&argc,&args,(char *)0,help);
30: #if defined(PETSC_USE_COMPLEX)
31: SETERRQ(1,"This example does not work with complex numbers");
32: #else
33: PetscOptionsHasName(PETSC_NULL,"-table",&table);
36: /* Read matrix and RHS */
37: PetscOptionsGetString(PETSC_NULL,"-f",file,127,&flg);
38: if (!flg) SETERRQ(1,"Must indicate binary file with the -f option");
39: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,PETSC_FILE_RDONLY,&fd);
41: MatLoad(fd,MATSEQAIJ,&A);
42: VecLoad(fd,PETSC_NULL,&b);
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: int 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(&zero,x);
75: PetscBarrier((PetscObject)A);
77: PetscLogStageRegister(&stage1,"mystage 1");
78: PetscLogStagePush(stage1);
79: PetscGetTime(&tsetup1);
80: KSPCreate(PETSC_COMM_WORLD,&ksp);
81: KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
82: KSPSetFromOptions(ksp);
83: KSPSetRhs(ksp,b);
84: KSPSetSolution(ksp,x);
85: KSPSetUp(ksp);
86: KSPSetUpOnBlocks(ksp);
87: PetscGetTime(&tsetup2);
88: tsetup = tsetup2 -tsetup1;
89: PetscLogStagePop();
90: PetscBarrier((PetscObject)A);
92: PetscLogStageRegister(&stage2,"mystage 2");
93: PetscLogStagePush(stage2);
94: PetscGetTime(&tsolve1);
95: KSPSolve(ksp);
96: PetscGetTime(&tsolve2);
97: tsolve = tsolve2 - tsolve1;
98: PetscLogStagePop();
100: /* Show result */
101: MatMult(A,x,u);
102: VecAXPY(&none,b,u);
103: VecNorm(u,NORM_2,&norm);
104: KSPGetIterationNumber(ksp,&its);
105: /* matrix PC KSP Options its residual setuptime solvetime */
106: if (table) {
107: char *matrixname,kspinfo[120];
108: PetscViewer viewer;
109: PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,120,&viewer);
110: KSPView(ksp,viewer);
111: PetscStrrchr(file,'/',&matrixname);
112: PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3d %2.0e %2.1e %2.1e %2.1e %s \n",
113: matrixname,its,norm,tsetup+tsolve,tsetup,tsolve,kspinfo);
114: PetscViewerDestroy(viewer);
115: } else {
116: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3d\n",its);
117: PetscPrintf(PETSC_COMM_WORLD,"Residual norm = %A\n",norm);
118: }
120: /* Cleanup */
121: KSPDestroy(ksp);
122: VecDestroy(x);
123: VecDestroy(b);
124: VecDestroy(u);
125: MatDestroy(A);
127: PetscFinalize();
128: #endif
129: return 0;
130: }