Actual source code: ex18.c
petsc-3.4.2 2013-07-02
2: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
3: Input arguments are:\n\
4: -f <input_file> : file to load. For example see $PETSC_DIR/share/petsc/datafiles/matrices\n\n";
6: #include <petscmat.h>
7: #include <petscksp.h>
8: #include <petsctime.h>
12: int main(int argc,char **args)
13: {
15: PetscInt its,m,n,mvec;
16: PetscLogDouble time1,time2,time;
17: PetscReal norm;
18: Vec x,b,u;
19: Mat A;
20: KSP ksp;
21: char file[PETSC_MAX_PATH_LEN];
22: PetscViewer fd;
23: PetscLogStage stage1;
25: PetscInitialize(&argc,&args,(char*)0,help);
27: /* Read matrix and RHS */
28: PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,NULL);
29: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
30: MatCreate(PETSC_COMM_WORLD,&A);
31: MatSetType(A,MATSEQAIJ);
32: MatLoad(A,fd);
33: VecCreate(PETSC_COMM_WORLD,&b);
34: VecLoad(b,fd);
35: PetscViewerDestroy(&fd);
37: /*
38: If the load matrix is larger then the vector, due to being padded
39: to match the blocksize then create a new padded vector
40: */
41: MatGetSize(A,&m,&n);
42: VecGetSize(b,&mvec);
43: if (m > mvec) {
44: Vec tmp;
45: PetscScalar *bold,*bnew;
46: /* create a new vector b by padding the old one */
47: VecCreate(PETSC_COMM_WORLD,&tmp);
48: VecSetSizes(tmp,PETSC_DECIDE,m);
49: VecSetFromOptions(tmp);
50: VecGetArray(tmp,&bnew);
51: VecGetArray(b,&bold);
52: PetscMemcpy(bnew,bold,mvec*sizeof(PetscScalar));
53: VecDestroy(&b);
54: b = tmp;
55: }
57: /* Set up solution */
58: VecDuplicate(b,&x);
59: VecDuplicate(b,&u);
60: VecSet(x,0.0);
62: /* Solve system */
63: PetscLogStageRegister("Stage 1",&stage1);
64: PetscLogStagePush(stage1);
65: KSPCreate(PETSC_COMM_WORLD,&ksp);
66: KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
67: KSPSetFromOptions(ksp);
68: PetscTime(&time1);
69: KSPSolve(ksp,b,x);
70: PetscTime(&time2);
71: time = time2 - time1;
72: PetscLogStagePop();
74: /* Show result */
75: MatMult(A,x,u);
76: VecAXPY(u,-1.0,b);
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 %G\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: }