Actual source code: ex18.c
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>
9: int main(int argc,char **args)
10: {
11: PetscInt its,m,n,mvec;
12: PetscReal norm;
13: Vec x,b,u;
14: Mat A;
15: KSP ksp;
16: char file[PETSC_MAX_PATH_LEN];
17: PetscViewer fd;
18: #if defined(PETSC_USE_LOG)
19: PetscLogStage stage1;
20: #endif
22: PetscInitialize(&argc,&args,(char*)0,help);
24: /* Read matrix and RHS */
25: PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),NULL);
26: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
27: MatCreate(PETSC_COMM_WORLD,&A);
28: MatSetType(A,MATSEQAIJ);
29: MatLoad(A,fd);
30: VecCreate(PETSC_COMM_WORLD,&b);
31: VecLoad(b,fd);
32: PetscViewerDestroy(&fd);
34: /*
35: If the load matrix is larger then the vector, due to being padded
36: to match the blocksize then create a new padded vector
37: */
38: MatGetSize(A,&m,&n);
39: VecGetSize(b,&mvec);
40: if (m > mvec) {
41: Vec tmp;
42: PetscScalar *bold,*bnew;
43: /* create a new vector b by padding the old one */
44: VecCreate(PETSC_COMM_WORLD,&tmp);
45: VecSetSizes(tmp,PETSC_DECIDE,m);
46: VecSetFromOptions(tmp);
47: VecGetArray(tmp,&bnew);
48: VecGetArray(b,&bold);
49: PetscArraycpy(bnew,bold,mvec);
50: VecDestroy(&b);
51: b = tmp;
52: }
54: /* Set up solution */
55: VecDuplicate(b,&x);
56: VecDuplicate(b,&u);
57: VecSet(x,0.0);
59: /* Solve system */
60: PetscLogStageRegister("Stage 1",&stage1);
61: PetscLogStagePush(stage1);
62: KSPCreate(PETSC_COMM_WORLD,&ksp);
63: KSPSetOperators(ksp,A,A);
64: KSPSetFromOptions(ksp);
65: KSPSolve(ksp,b,x);
66: PetscLogStagePop();
68: /* Show result */
69: MatMult(A,x,u);
70: VecAXPY(u,-1.0,b);
71: VecNorm(u,NORM_2,&norm);
72: KSPGetIterationNumber(ksp,&its);
73: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
74: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
76: /* Cleanup */
77: KSPDestroy(&ksp);
78: VecDestroy(&x);
79: VecDestroy(&b);
80: VecDestroy(&u);
81: MatDestroy(&A);
83: PetscFinalize();
84: return 0;
85: }
87: /*TEST
89: test:
90: args: -ksp_gmres_cgs_refinement_type refine_always -f ${DATAFILESPATH}/matrices/arco1 -ksp_monitor_short
91: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
93: TEST*/