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++) {
 37:     cnt[n] += (vv[i] != 0.0);
 38:   }
 39:   MatCreateSeqAIJ(PETSC_COMM_SELF,n+1,n+1,0,cnt,B);
 40:   MatSetOption(*B,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);

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

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

 60: int main(int argc,char **args)
 61: {
 62:   Mat            A,B;
 63:   PetscViewer    fd;                        /* viewer */
 64:   char           file[PETSC_MAX_PATH_LEN];  /* input file name */
 65:   PetscBool      flg;
 66:   Vec            v;

 68:   PetscInitialize(&argc,&args,(char*)0,help);
 69:   /*
 70:      Determine files from which we read the two linear systems
 71:      (matrix and right-hand-side vector).
 72:   */
 73:   PetscOptionsGetString(NULL,NULL,"-f0",file,sizeof(file),&flg);

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

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

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

 95: /*TEST

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

101: TEST*/