Actual source code: ex18.c
1: /*$Id: ex18.c,v 1.27 2001/08/07 21:30:50 bsmith Exp $*/
3: #if !defined(PETSC_USE_COMPLEX)
5: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
6: Input arguments are:\n\
7: -f <input_file> : file to load. For a 5X5 example of the 5-pt. stencil,\n\
8: use the file petsc/src/mat/examples/matbinary.ex\n\n";
10: #include petscmat.h
11: #include petscksp.h
15: int main(int argc,char **args)
16: {
17: int ierr,its,m,n,mvec;
18: PetscLogDouble time1,time2,time;
19: PetscReal norm;
20: PetscScalar zero = 0.0,none = -1.0;
21: Vec x,b,u;
22: Mat A;
23: KSP ksp;
24: char file[128];
25: PetscViewer fd;
26:
27: PetscInitialize(&argc,&args,(char *)0,help);
29: /* Read matrix and RHS */
30: PetscOptionsGetString(PETSC_NULL,"-f",file,127,PETSC_NULL);
31: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,PETSC_FILE_RDONLY,&fd);
32: MatLoad(fd,MATSEQAIJ,&A);
33: VecLoad(fd,PETSC_NULL,&b);
34: PetscViewerDestroy(fd);
36: /*
37: If the load matrix is larger then the vector, due to being padded
38: to match the blocksize then create a new padded vector
39: */
40: MatGetSize(A,&m,&n);
41: VecGetSize(b,&mvec);
42: if (m > mvec) {
43: Vec tmp;
44: PetscScalar *bold,*bnew;
45: /* create a new vector b by padding the old one */
46: VecCreate(PETSC_COMM_WORLD,&tmp);
47: VecSetSizes(tmp,PETSC_DECIDE,m);
48: VecSetFromOptions(tmp);
49: VecGetArray(tmp,&bnew);
50: VecGetArray(b,&bold);
51: PetscMemcpy(bnew,bold,mvec*sizeof(PetscScalar));
52: VecDestroy(b);
53: b = tmp;
54: }
56: /* Set up solution */
57: VecDuplicate(b,&x);
58: VecDuplicate(b,&u);
59: VecSet(&zero,x);
61: /* Solve system */
62: PetscLogStagePush(1);
63: KSPCreate(PETSC_COMM_WORLD,&ksp);
64: KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
65: KSPSetFromOptions(ksp);
66: PetscGetTime(&time1);
67: KSPSetRhs(ksp,b);
68: KSPSetSolution(ksp,x);
69: KSPSolve(ksp);
70: PetscGetTime(&time2);
71: time = time2 - time1;
72: PetscLogStagePop();
74: /* Show result */
75: MatMult(A,x,u);
76: VecAXPY(&none,b,u);
77: VecNorm(u,NORM_2,&norm);
78: KSPGetIterationNumber(ksp,&its);
79: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3d\n",its);
80: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %A\n",norm);
81: PetscPrintf(PETSC_COMM_WORLD,"Time for solve = %5.2f seconds\n",time);
83: /* Cleanup */
84: KSPDestroy(ksp);
85: VecDestroy(x);
86: VecDestroy(b);
87: VecDestroy(u);
88: MatDestroy(A);
90: PetscFinalize();
91: return 0;
92: }
94: #else
95: #include <stdio.h>
96: int main(int argc,char **args)
97: {
98: fprintf(stdout,"This example does not work for complex numbers.\n");
99: return 0;
100: }
101: #endif