Actual source code: ex12.c
2: static char help[] = "Reads a PETSc matrix and vector from a file; expands the matrix with the vector\n\n";
4: /*
5: Include "petscmat.h" so that we can use matrices.
6: automatically includes:
7: petscsys.h - base PETSc routines petscvec.h - vectors
8: petscmat.h - matrices
9: petscis.h - index sets petscviewer.h - viewers
10: */
11: #include <petscmat.h>
13: /*
15: Adds a new column and row to the vector (the last) containing the vector
16: */
17: PetscErrorCode PadMatrix(Mat A, Vec v, PetscScalar c, Mat *B)
18: {
19: PetscInt n, i, *cnt, *indices, nc;
20: const PetscInt *aj;
21: const PetscScalar *vv, *aa;
23: PetscFunctionBegin;
24: PetscCall(MatGetSize(A, &n, NULL));
25: PetscCall(VecGetArrayRead(v, &vv));
26: PetscCall(PetscMalloc1(n, &indices));
27: for (i = 0; i < n; i++) indices[i] = i;
29: /* determine number of nonzeros per row in the new matrix */
30: PetscCall(PetscMalloc1(n + 1, &cnt));
31: for (i = 0; i < n; i++) {
32: PetscCall(MatGetRow(A, i, &nc, NULL, NULL));
33: cnt[i] = nc + (vv[i] != 0.0);
34: PetscCall(MatRestoreRow(A, i, &nc, NULL, NULL));
35: }
36: cnt[n] = 1;
37: for (i = 0; i < n; i++) cnt[n] += (vv[i] != 0.0);
38: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n + 1, n + 1, 0, cnt, B));
39: PetscCall(MatSetOption(*B, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
41: /* copy over the matrix entries from the matrix and then the vector */
42: for (i = 0; i < n; i++) {
43: PetscCall(MatGetRow(A, i, &nc, &aj, &aa));
44: PetscCall(MatSetValues(*B, 1, &i, nc, aj, aa, INSERT_VALUES));
45: PetscCall(MatRestoreRow(A, i, &nc, &aj, &aa));
46: }
47: PetscCall(MatSetValues(*B, 1, &n, n, indices, vv, INSERT_VALUES));
48: PetscCall(MatSetValues(*B, n, indices, 1, &n, vv, INSERT_VALUES));
49: PetscCall(MatSetValues(*B, 1, &n, 1, &n, &c, INSERT_VALUES));
51: PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
52: PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
53: PetscCall(VecRestoreArrayRead(v, &vv));
54: PetscCall(PetscFree(cnt));
55: PetscCall(PetscFree(indices));
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: int main(int argc, char **args)
60: {
61: Mat A, B;
62: PetscViewer fd; /* viewer */
63: char file[PETSC_MAX_PATH_LEN]; /* input file name */
64: PetscBool flg;
65: Vec v;
67: PetscFunctionBeginUser;
68: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
69: /*
70: Determine files from which we read the two linear systems
71: (matrix and right-hand-side vector).
72: */
73: PetscCall(PetscOptionsGetString(NULL, NULL, "-f0", file, sizeof(file), &flg));
74: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f0 option");
76: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
78: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
79: PetscCall(MatSetType(A, MATSEQAIJ));
80: PetscCall(MatLoad(A, fd));
81: PetscCall(VecCreate(PETSC_COMM_WORLD, &v));
82: PetscCall(VecLoad(v, fd));
83: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF));
84: PetscCall(PadMatrix(A, v, 3.0, &B));
85: PetscCall(MatView(B, PETSC_VIEWER_STDOUT_SELF));
86: PetscCall(MatDestroy(&B));
87: PetscCall(MatDestroy(&A));
88: PetscCall(VecDestroy(&v));
89: PetscCall(PetscViewerDestroy(&fd));
91: PetscCall(PetscFinalize());
92: return 0;
93: }
95: /*TEST
97: test:
98: args: -f0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/ns-real-int32-float64
99: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
101: TEST*/