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 */ 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 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, (char *)0, 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