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 PetscBool flg; 75 Vec v; 76 77 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd)); 86 87 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 88 CHKERRQ(MatSetType(A,MATSEQAIJ)); 89 CHKERRQ(MatLoad(A,fd)); 90 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&v)); 91 CHKERRQ(VecLoad(v,fd)); 92 CHKERRQ(MatView(A,PETSC_VIEWER_STDOUT_SELF)); 93 CHKERRQ(PadMatrix(A,v,3.0,&B)); 94 CHKERRQ(MatView(B,PETSC_VIEWER_STDOUT_SELF)); 95 CHKERRQ(MatDestroy(&B)); 96 CHKERRQ(MatDestroy(&A)); 97 CHKERRQ(VecDestroy(&v)); 98 CHKERRQ(PetscViewerDestroy(&fd)); 99 100 CHKERRQ(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