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