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 petscsles.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: SLES sles;
23: char file[128];
24: PetscViewer fd;
25: PetscTruth table,flg;
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_BINARY_RDONLY,&fd);
41: MatLoad(fd,MATSEQAIJ,&A);
42: VecLoad(fd,&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: SLESCreate(PETSC_COMM_WORLD,&sles);
81: SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTERN);
82: SLESSetFromOptions(sles);
83: SLESSetUp(sles,b,x);
84: SLESSetUpOnBlocks(sles);
85: PetscGetTime(&tsetup2);
86: tsetup = tsetup2 -tsetup1;
87: PetscLogStagePop();
88: PetscBarrier((PetscObject)A);
90: PetscLogStageRegister(&stage2,"mystage 2");
91: PetscLogStagePush(stage2);
92: PetscGetTime(&tsolve1);
93: SLESSolve(sles,b,x,&its);
94: PetscGetTime(&tsolve2);
95: tsolve = tsolve2 - tsolve1;
96: PetscLogStagePop();
98: /* Show result */
99: MatMult(A,x,u);
100: VecAXPY(&none,b,u);
101: VecNorm(u,NORM_2,&norm);
102: /* matrix PC KSP Options its residual setuptime solvetime */
103: if (table) {
104: char *matrixname,slesinfo[120];
105: PetscViewer viewer;
106: PetscViewerStringOpen(PETSC_COMM_WORLD,slesinfo,120,&viewer);
107: SLESView(sles,viewer);
108: PetscStrrchr(file,'/',&matrixname);
109: PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3d %2.0e %2.1e %2.1e %2.1e %s \n",
110: matrixname,its,norm,tsetup+tsolve,tsetup,tsolve,slesinfo);
111: PetscViewerDestroy(viewer);
112: } else {
113: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3d\n",its);
114: PetscPrintf(PETSC_COMM_WORLD,"Residual norm = %A\n",norm);
115: }
117: /* Cleanup */
118: SLESDestroy(sles);
119: VecDestroy(x);
120: VecDestroy(b);
121: VecDestroy(u);
122: MatDestroy(A);
124: PetscFinalize();
125: #endif
126: return 0;
127: }