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