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: PetscLogStage stage1;
26:
27: PetscInitialize(&argc,&args,(char *)0,help);
29: /* Read matrix and RHS */
30: PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN,PETSC_NULL);
31: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
32: MatCreate(PETSC_COMM_WORLD,&A);
33: MatSetType(A,MATSEQAIJ);
34: MatLoad(A,fd);
35: VecCreate(PETSC_COMM_WORLD,&b);
36: VecLoad(b,fd);
37: PetscViewerDestroy(&fd);
39: /*
40: If the load matrix is larger then the vector, due to being padded
41: to match the blocksize then create a new padded vector
42: */
43: MatGetSize(A,&m,&n);
44: VecGetSize(b,&mvec);
45: if (m > mvec) {
46: Vec tmp;
47: PetscScalar *bold,*bnew;
48: /* create a new vector b by padding the old one */
49: VecCreate(PETSC_COMM_WORLD,&tmp);
50: VecSetSizes(tmp,PETSC_DECIDE,m);
51: VecSetFromOptions(tmp);
52: VecGetArray(tmp,&bnew);
53: VecGetArray(b,&bold);
54: PetscMemcpy(bnew,bold,mvec*sizeof(PetscScalar));
55: VecDestroy(&b);
56: b = tmp;
57: }
59: /* Set up solution */
60: VecDuplicate(b,&x);
61: VecDuplicate(b,&u);
62: VecSet(x,0.0);
64: /* Solve system */
65: PetscLogStageRegister("Stage 1",&stage1);
66: PetscLogStagePush(stage1);
67: KSPCreate(PETSC_COMM_WORLD,&ksp);
68: KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
69: KSPSetFromOptions(ksp);
70: PetscGetTime(&time1);
71: KSPSolve(ksp,b,x);
72: PetscGetTime(&time2);
73: time = time2 - time1;
74: PetscLogStagePop();
76: /* Show result */
77: MatMult(A,x,u);
78: VecAXPY(u,-1.0,b);
79: VecNorm(u,NORM_2,&norm);
80: KSPGetIterationNumber(ksp,&its);
81: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
82: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %A\n",norm);
83: PetscPrintf(PETSC_COMM_WORLD,"Time for solve = %5.2f seconds\n",time);
85: /* Cleanup */
86: KSPDestroy(&ksp);
87: VecDestroy(&x);
88: VecDestroy(&b);
89: VecDestroy(&u);
90: MatDestroy(&A);
92: PetscFinalize();
93: return 0;
94: }
96: #else
97: #include <stdio.h>
98: int main(int argc,char **args)
99: {
100: fprintf(stdout,"This example does not work for complex numbers.\n");
101: return 0;
102: }
103: #endif