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:   MatGetSize(A, &n, NULL);
 24:   VecGetArrayRead(v, &vv);
 25:   PetscMalloc1(n, &indices);
 26:   for (i = 0; i < n; i++) indices[i] = i;

 28:   /* determine number of nonzeros per row in the new matrix */
 29:   PetscMalloc1(n + 1, &cnt);
 30:   for (i = 0; i < n; i++) {
 31:     MatGetRow(A, i, &nc, NULL, NULL);
 32:     cnt[i] = nc + (vv[i] != 0.0);
 33:     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:   MatCreateSeqAIJ(PETSC_COMM_SELF, n + 1, n + 1, 0, cnt, B);
 38:   MatSetOption(*B, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);

 40:   /* copy over the matrix entries from the matrix and then the vector */
 41:   for (i = 0; i < n; i++) {
 42:     MatGetRow(A, i, &nc, &aj, &aa);
 43:     MatSetValues(*B, 1, &i, nc, aj, aa, INSERT_VALUES);
 44:     MatRestoreRow(A, i, &nc, &aj, &aa);
 45:   }
 46:   MatSetValues(*B, 1, &n, n, indices, vv, INSERT_VALUES);
 47:   MatSetValues(*B, n, indices, 1, &n, vv, INSERT_VALUES);
 48:   MatSetValues(*B, 1, &n, 1, &n, &c, INSERT_VALUES);

 50:   MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
 51:   MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY);
 52:   VecRestoreArrayRead(v, &vv);
 53:   PetscFree(cnt);
 54:   PetscFree(indices);
 55:   return 0;
 56: }

 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;

 67:   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:   PetscOptionsGetString(NULL, NULL, "-f0", file, sizeof(file), &flg);

 75:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd);

 77:   MatCreate(PETSC_COMM_WORLD, &A);
 78:   MatSetType(A, MATSEQAIJ);
 79:   MatLoad(A, fd);
 80:   VecCreate(PETSC_COMM_WORLD, &v);
 81:   VecLoad(v, fd);
 82:   MatView(A, PETSC_VIEWER_STDOUT_SELF);
 83:   PadMatrix(A, v, 3.0, &B);
 84:   MatView(B, PETSC_VIEWER_STDOUT_SELF);
 85:   MatDestroy(&B);
 86:   MatDestroy(&A);
 87:   VecDestroy(&v);
 88:   PetscViewerDestroy(&fd);

 90:   PetscFinalize();
 91:   return 0;
 92: }

 94: /*TEST

 96:    test:
 97:       args: -f0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/ns-real-int32-float64
 98:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)

100: TEST*/