Actual source code: ex18.c
petsc-3.3-p6 2013-02-11
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 example see $PETSC_DIR/share/petsc/datafiles/matrices\n\n";
8: #include <petscmat.h>
9: #include <petscksp.h>
13: int main(int argc,char **args)
14: {
16: PetscInt its,m,n,mvec;
17: PetscLogDouble time1,time2,time;
18: PetscReal norm;
19: Vec x,b,u;
20: Mat A;
21: KSP ksp;
22: char file[PETSC_MAX_PATH_LEN];
23: PetscViewer fd;
24: PetscLogStage stage1;
25:
26: PetscInitialize(&argc,&args,(char *)0,help);
28: /* Read matrix and RHS */
29: PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN,PETSC_NULL);
30: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
31: MatCreate(PETSC_COMM_WORLD,&A);
32: MatSetType(A,MATSEQAIJ);
33: MatLoad(A,fd);
34: VecCreate(PETSC_COMM_WORLD,&b);
35: VecLoad(b,fd);
36: PetscViewerDestroy(&fd);
38: /*
39: If the load matrix is larger then the vector, due to being padded
40: to match the blocksize then create a new padded vector
41: */
42: MatGetSize(A,&m,&n);
43: VecGetSize(b,&mvec);
44: if (m > mvec) {
45: Vec tmp;
46: PetscScalar *bold,*bnew;
47: /* create a new vector b by padding the old one */
48: VecCreate(PETSC_COMM_WORLD,&tmp);
49: VecSetSizes(tmp,PETSC_DECIDE,m);
50: VecSetFromOptions(tmp);
51: VecGetArray(tmp,&bnew);
52: VecGetArray(b,&bold);
53: PetscMemcpy(bnew,bold,mvec*sizeof(PetscScalar));
54: VecDestroy(&b);
55: b = tmp;
56: }
58: /* Set up solution */
59: VecDuplicate(b,&x);
60: VecDuplicate(b,&u);
61: VecSet(x,0.0);
63: /* Solve system */
64: PetscLogStageRegister("Stage 1",&stage1);
65: PetscLogStagePush(stage1);
66: KSPCreate(PETSC_COMM_WORLD,&ksp);
67: KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
68: KSPSetFromOptions(ksp);
69: PetscGetTime(&time1);
70: KSPSolve(ksp,b,x);
71: PetscGetTime(&time2);
72: time = time2 - time1;
73: PetscLogStagePop();
75: /* Show result */
76: MatMult(A,x,u);
77: VecAXPY(u,-1.0,b);
78: VecNorm(u,NORM_2,&norm);
79: KSPGetIterationNumber(ksp,&its);
80: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
81: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);
82: PetscPrintf(PETSC_COMM_WORLD,"Time for solve = %5.2f seconds\n",time);
84: /* Cleanup */
85: KSPDestroy(&ksp);
86: VecDestroy(&x);
87: VecDestroy(&b);
88: VecDestroy(&u);
89: MatDestroy(&A);
91: PetscFinalize();
92: return 0;
93: }
95: #else
96: #include <stdio.h>
97: int main(int argc,char **args)
98: {
99: fprintf(stdout,"This example does not work for complex numbers.\n");
100: return 0;
101: }
102: #endif