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