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