1 2 static char help[] = "Reads a PETSc matrix and vector from a file; expands the matrix with the vector\n\n"; 3 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> 12 13 /* 14 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; 22 23 PetscFunctionBegin; 24 PetscCall(MatGetSize(A,&n,NULL)); 25 PetscCall(VecGetArrayRead(v,&vv)); 26 PetscCall(PetscMalloc1(n,&indices)); 27 for (i=0; i<n; i++) indices[i] = i; 28 29 /* determine number of nonzeros per row in the new matrix */ 30 PetscCall(PetscMalloc1(n+1,&cnt)); 31 for (i=0; i<n; i++) { 32 PetscCall(MatGetRow(A,i,&nc,NULL,NULL)); 33 cnt[i] = nc + (vv[i] != 0.0); 34 PetscCall(MatRestoreRow(A,i,&nc,NULL,NULL)); 35 } 36 cnt[n] = 1; 37 for (i=0; i<n; i++) { 38 cnt[n] += (vv[i] != 0.0); 39 } 40 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF,n+1,n+1,0,cnt,B)); 41 PetscCall(MatSetOption(*B,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE)); 42 43 /* copy over the matrix entries from the matrix and then the vector */ 44 for (i=0; i<n; i++) { 45 PetscCall(MatGetRow(A,i,&nc,&aj,&aa)); 46 PetscCall(MatSetValues(*B,1,&i,nc,aj,aa,INSERT_VALUES)); 47 PetscCall(MatRestoreRow(A,i,&nc,&aj,&aa)); 48 } 49 PetscCall(MatSetValues(*B,1,&n,n,indices,vv,INSERT_VALUES)); 50 PetscCall(MatSetValues(*B,n,indices,1,&n,vv,INSERT_VALUES)); 51 PetscCall(MatSetValues(*B,1,&n,1,&n,&c,INSERT_VALUES)); 52 53 PetscCall(MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY)); 54 PetscCall(MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY)); 55 PetscCall(VecRestoreArrayRead(v,&vv)); 56 PetscCall(PetscFree(cnt)); 57 PetscCall(PetscFree(indices)); 58 PetscFunctionReturn(0); 59 } 60 61 int main(int argc,char **args) 62 { 63 Mat A,B; 64 PetscViewer fd; /* viewer */ 65 char file[PETSC_MAX_PATH_LEN]; /* input file name */ 66 PetscBool flg; 67 Vec v; 68 69 PetscFunctionBeginUser; 70 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 71 /* 72 Determine files from which we read the two linear systems 73 (matrix and right-hand-side vector). 74 */ 75 PetscCall(PetscOptionsGetString(NULL,NULL,"-f0",file,sizeof(file),&flg)); 76 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f0 option"); 77 78 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd)); 79 80 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 81 PetscCall(MatSetType(A,MATSEQAIJ)); 82 PetscCall(MatLoad(A,fd)); 83 PetscCall(VecCreate(PETSC_COMM_WORLD,&v)); 84 PetscCall(VecLoad(v,fd)); 85 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF)); 86 PetscCall(PadMatrix(A,v,3.0,&B)); 87 PetscCall(MatView(B,PETSC_VIEWER_STDOUT_SELF)); 88 PetscCall(MatDestroy(&B)); 89 PetscCall(MatDestroy(&A)); 90 PetscCall(VecDestroy(&v)); 91 PetscCall(PetscViewerDestroy(&fd)); 92 93 PetscCall(PetscFinalize()); 94 return 0; 95 } 96 97 /*TEST 98 99 test: 100 args: -f0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/ns-real-int32-float64 101 requires: double !complex !defined(PETSC_USE_64BIT_INDICES) 102 103 TEST*/ 104