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

  4: /*T
  5:    Concepts: Mat^ordering a matrix - loading a binary matrix and vector;
  6:    Concepts: Mat^loading a binary matrix and vector;
  7:    Concepts: Vectors^loading a binary vector;
  8:    Concepts: PetscLog^preloading executable
  9:    Processors: 1
 10: T*/



 14: /*
 15:   Include "petscmat.h" so that we can use matrices.
 16:   automatically includes:
 17:      petscsys.h       - base PETSc routines   petscvec.h    - vectors
 18:      petscmat.h    - matrices
 19:      petscis.h     - index sets            petscviewer.h - viewers
 20: */
 21:  #include <petscmat.h>

 23: /*

 25:    Adds a new column and row to the vector (the last) containing the vector
 26: */
 27: PetscErrorCode PadMatrix(Mat A,Vec v,PetscScalar c,Mat *B)
 28: {
 29:   PetscErrorCode    ierr;
 30:   PetscInt          n,i,*cnt,*indices,nc;
 31:   const PetscInt    *aj;
 32:   const PetscScalar *vv,*aa;

 35:   MatGetSize(A,&n,NULL);
 36:   VecGetArrayRead(v,&vv);
 37:   PetscMalloc1(n,&indices);
 38:   for (i=0; i<n; i++) indices[i] = i;

 40:   /* determine number of nonzeros per row in the new matrix */
 41:   PetscMalloc1(n+1,&cnt);
 42:   for (i=0; i<n; i++) {
 43:     MatGetRow(A,i,&nc,NULL,NULL);
 44:     cnt[i] = nc + (vv[i] != 0.0);
 45:     MatRestoreRow(A,i,&nc,NULL,NULL);
 46:   }
 47:   cnt[n] = 1;
 48:   for (i=0; i<n; i++) {
 49:     cnt[n] += (vv[i] != 0.0);
 50:   }
 51:   MatCreateSeqAIJ(PETSC_COMM_SELF,n+1,n+1,0,cnt,B);
 52:   MatSetOption(*B,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);

 54:   /* copy over the matrix entries from the matrix and then the vector */
 55:   for (i=0; i<n; i++) {
 56:     MatGetRow(A,i,&nc,&aj,&aa);
 57:     MatSetValues(*B,1,&i,nc,aj,aa,INSERT_VALUES);
 58:     MatRestoreRow(A,i,&nc,&aj,&aa);
 59:   }
 60:   MatSetValues(*B,1,&n,n,indices,vv,INSERT_VALUES);
 61:   MatSetValues(*B,n,indices,1,&n,vv,INSERT_VALUES);
 62:   MatSetValues(*B,1,&n,1,&n,&c,INSERT_VALUES);

 64:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
 65:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
 66:   VecRestoreArrayRead(v,&vv);
 67:   PetscFree(cnt);
 68:   PetscFree(indices);
 69:   return(0);
 70: }


 73: int main(int argc,char **args)
 74: {
 75:   Mat            A,B;
 76:   PetscViewer    fd;                        /* viewer */
 77:   char           file[PETSC_MAX_PATH_LEN];  /* input file name */
 79:   PetscBool      flg;
 80:   Vec            v;

 82:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 83:   /*
 84:      Determine files from which we read the two linear systems
 85:      (matrix and right-hand-side vector).
 86:   */
 87:   PetscOptionsGetString(NULL,NULL,"-f0",file,PETSC_MAX_PATH_LEN,&flg);
 88:   if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f0 option");

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

 92:   MatCreate(PETSC_COMM_WORLD,&A);
 93:   MatSetType(A,MATSEQAIJ);
 94:   MatLoad(A,fd);
 95:   VecCreate(PETSC_COMM_WORLD,&v);
 96:   VecLoad(v,fd);
 97:   MatView(A,PETSC_VIEWER_STDOUT_SELF);
 98:   PadMatrix(A,v,3.0,&B);
 99:   MatView(B,PETSC_VIEWER_STDOUT_SELF);
100:   MatDestroy(&B);
101:   MatDestroy(&A);
102:   VecDestroy(&v);
103:   PetscViewerDestroy(&fd);

105:   PetscFinalize();
106:   return ierr;
107: }



111: /*TEST

113:    test:
114:       args: -f0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/ns-real-int32-float64
115:       requires: double !complex !define(PETSC_USE_64BIT_INDICES)

117: TEST*/