Actual source code: ex12.c

      1: static char help[] = "Reads a PETSc matrix and vector from a file; expands the matrix with the vector\n\n";

      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>

     12: /*

     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;

     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;

     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));

     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));

     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: }

     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;

     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");

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

     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));

     90:   PetscCall(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*/