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 CHKERRQ(MatGetSize(A,&n,NULL)); 33 CHKERRQ(VecGetArrayRead(v,&vv)); 34 CHKERRQ(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 CHKERRQ(PetscMalloc1(n+1,&cnt)); 39 for (i=0; i<n; i++) { 40 CHKERRQ(MatGetRow(A,i,&nc,NULL,NULL)); 41 cnt[i] = nc + (vv[i] != 0.0); 42 CHKERRQ(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 CHKERRQ(MatCreateSeqAIJ(PETSC_COMM_SELF,n+1,n+1,0,cnt,B)); 49 CHKERRQ(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 CHKERRQ(MatGetRow(A,i,&nc,&aj,&aa)); 54 CHKERRQ(MatSetValues(*B,1,&i,nc,aj,aa,INSERT_VALUES)); 55 CHKERRQ(MatRestoreRow(A,i,&nc,&aj,&aa)); 56 } 57 CHKERRQ(MatSetValues(*B,1,&n,n,indices,vv,INSERT_VALUES)); 58 CHKERRQ(MatSetValues(*B,n,indices,1,&n,vv,INSERT_VALUES)); 59 CHKERRQ(MatSetValues(*B,1,&n,1,&n,&c,INSERT_VALUES)); 60 61 CHKERRQ(MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY)); 62 CHKERRQ(MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY)); 63 CHKERRQ(VecRestoreArrayRead(v,&vv)); 64 CHKERRQ(PetscFree(cnt)); 65 CHKERRQ(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 PetscErrorCode ierr; 75 PetscBool flg; 76 Vec v; 77 78 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 79 /* 80 Determine files from which we read the two linear systems 81 (matrix and right-hand-side vector). 82 */ 83 CHKERRQ(PetscOptionsGetString(NULL,NULL,"-f0",file,sizeof(file),&flg)); 84 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f0 option"); 85 86 CHKERRQ(PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd)); 87 88 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 89 CHKERRQ(MatSetType(A,MATSEQAIJ)); 90 CHKERRQ(MatLoad(A,fd)); 91 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&v)); 92 CHKERRQ(VecLoad(v,fd)); 93 CHKERRQ(MatView(A,PETSC_VIEWER_STDOUT_SELF)); 94 CHKERRQ(PadMatrix(A,v,3.0,&B)); 95 CHKERRQ(MatView(B,PETSC_VIEWER_STDOUT_SELF)); 96 CHKERRQ(MatDestroy(&B)); 97 CHKERRQ(MatDestroy(&A)); 98 CHKERRQ(VecDestroy(&v)); 99 CHKERRQ(PetscViewerDestroy(&fd)); 100 101 ierr = PetscFinalize(); 102 return ierr; 103 } 104 105 /*TEST 106 107 test: 108 args: -f0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/ns-real-int32-float64 109 requires: double !complex !defined(PETSC_USE_64BIT_INDICES) 110 111 TEST*/ 112