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
23: PetscInitialize(&argc, &args, (char *)0, help);
25: /* Read matrix and RHS */
26: PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), NULL);
27: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd);
28: MatCreate(PETSC_COMM_WORLD, &A);
29: MatSetType(A, MATSEQAIJ);
30: MatLoad(A, fd);
31: VecCreate(PETSC_COMM_WORLD, &b);
32: VecLoad(b, fd);
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: PetscArraycpy(bnew, bold, mvec);
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: PetscLogStageRegister("Stage 1", &stage1);
62: PetscLogStagePush(stage1);
63: KSPCreate(PETSC_COMM_WORLD, &ksp);
64: KSPSetOperators(ksp, A, A);
65: KSPSetFromOptions(ksp);
66: KSPSolve(ksp, b, x);
67: PetscLogStagePop();
69: /* Show result */
70: MatMult(A, x, u);
71: VecAXPY(u, -1.0, b);
72: VecNorm(u, NORM_2, &norm);
73: KSPGetIterationNumber(ksp, &its);
74: PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its);
75: PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g\n", (double)norm);
77: /* Cleanup */
78: KSPDestroy(&ksp);
79: VecDestroy(&x);
80: VecDestroy(&b);
81: VecDestroy(&u);
82: MatDestroy(&A);
84: PetscFinalize();
85: return 0;
86: }
88: /*TEST
90: test:
91: args: -ksp_gmres_cgs_refinement_type refine_always -f ${DATAFILESPATH}/matrices/arco1 -ksp_monitor_short
92: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
94: TEST*/