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