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