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