On 11 February 2016 at 09:16, Sang pham van <pvsang...@gmail.com> wrote:
> I insert the lines below into my code, but it does not work: > Look more carefully at the calling pattern of the standard implementation DMCreateMatrix_DA() calls DMCreateMatrix_DA_3d_MPIAIJ() DMCreateMatrix_DA() allocates the matrix. If in your user code you simply call DMCreateMatrix_DA_3d_MPIAIJ_pvs() from your user code without first calling performing a similar setup and allocation as that occurring DMCreateMatrix_DA(), it clearly won't work. The error tells you everything. It indicates the first argument is not allocated, e.g. pMat has not been allocated. > PetscErrorCode ierr; > PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; > PetscInt m,n,dim,s,*cols = NULL,k,nc,*rows = > NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL; > PetscInt > istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; > MPI_Comm comm; > PetscScalar *values; > DMDABoundaryType bx,by,bz; > ISLocalToGlobalMapping ltog,ltogb; > DMDAStencilType st; > > PetscFunctionBegin; > ierr = > DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); > col = 2*s + 1; > > ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); > ierr = > DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); > ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); > > ierr = > PetscMalloc2(nc,PetscInt,&rows,col*col*col*nc*nc,PetscInt,&cols);CHKERRQ(ierr); > ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); > > ierr = MatSetBlockSize(pMat,nc);CHKERRQ(ierr); /// error happens ! > /* determine the matrix preallocation information */ > ierr = > MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);CHKERRQ(ierr); > for (i=xs; i<xs+nx; i++) { > istart = (bx == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); > iend = (bx == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); > for (j=ys; j<ys+ny; j++) { > jstart = (by == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); > jend = (by == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); > for (k=zs; k<zs+nz; k++) { > kstart = (bz == DMDA_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); > kend = (bz == DMDA_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); > > slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); > > cnt = 0; > for (l=0; l<nc; l++) { > for (ii=istart; ii<iend+1; ii++) { > for (jj=jstart; jj<jend+1; jj++) { > for (kk=kstart; kk<kend+1; kk++) { > if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && > !kk) || (!ii && !kk))) {/* entries on star*/ > cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk); > } > } > } > } > rows[l] = l + nc*(slot); > } > ierr = > MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); > } > } > } > ierr = MatSetBlockSize(pMat,nc);CHKERRQ(ierr); > ierr = MatSeqAIJSetPreallocation(pMat,0,dnz);CHKERRQ(ierr); > ierr = MatMPIAIJSetPreallocation(pMat,0,dnz,0,onz);CHKERRQ(ierr); > ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); > ierr = MatSetLocalToGlobalMapping(pMat,ltog,ltog);CHKERRQ(ierr); > ierr = PetscFree2(rows,cols);CHKERRQ(ierr); > > An error reported when it runs into ierr = > MatSetBlockSize(pMat,nc);CHKERRQ(ierr); > > "unknowndirectory/"src/mat_vec/sasMatVecPetsc.cpp:154: > __FUNCT__="DMCreateMatrix_DA_3d_MPIAIJ_pvs" does not agree with > __PRETTY_FUNCTION__="PetscErrorCode > sasMatVecPetsc::DMCreateMatrix_DA_3d_MPIAIJ_pvs(DM, sasSmesh*, > sasVector<int>&, sasVector<int>&, sasVector<int>&, sasVector<int>&)" > [0]PETSC ERROR: --------------------- Error Message > ------------------------------------ > [0]PETSC ERROR: Null argument, when expecting valid pointer! > [0]PETSC ERROR: Null Object: Parameter # 1! > [0]PETSC ERROR: > ------------------------------------------------------------------------ > [0]PETSC ERROR: Petsc Release Version 3.4.5, Jun, 29, 2014 > [0]PETSC ERROR: See docs/changes/index.html for recent updates. > [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting. > [0]PETSC ERROR: See docs/index.html for manual pages. > [0]PETSC ERROR: > ------------------------------------------------------------------------ > [0]PETSC ERROR: ./sasMAIN on a arch-cxx-linux-deb named sang by sang Thu > Feb 11 14:54:12 2016 > [0]PETSC ERROR: Libraries linked from > /home/sang/petsc/petsc-3.4.5/arch-cxx-linux-deb/lib > [0]PETSC ERROR: Configure run at Thu Feb 11 11:44:35 2016 > [0]PETSC ERROR: Configure options --download-mpich=1 > --with-scalar-type=real --with-clanguage=cxx --download-mumps=1 > --download-blacs=1 --download-parmetis=1 --download-scalapack=1 > --with-debugging=1 --download-hypre=1 --with-fc=gfortran --download-metis=1 > -download-cmake=1 --download-f-blas-lapack=1 > [0]PETSC ERROR: > ------------------------------------------------------------------------ > [0]PETSC ERROR: MatSetBlockSize() line 6686 in > /home/sang/petsc/petsc-3.4.5/src/mat/interface/matrix.c > [0]PETSC ERROR: PetscErrorCode > sasMatVecPetsc::DMCreateMatrix_DA_3d_MPIAIJ_pvs(DM, sasSmesh*, > sasVector<int>&, sasVector<int>&, sasVector<int>&, sasVector<int>&)() line > 165 in "unknowndirectory/"src/mat_vec/sasMatVecPetsc.cpp > > where is the mistake? > > Many thanks. > > Pham > > On Thu, Feb 11, 2016 at 1:18 PM, Dave May <dave.mayhe...@gmail.com> wrote: > >> I think he wants the source location so that he can copy and >> implementation and "tweak" it slightly >> >> The location is here >> ${PETSC_DIR}/src/dm/impls/da/fdda.c >> >> /Users/dmay/software/petsc-3.6.0/src >> dmay@nikkan:~/software/petsc-3.6.0/src $ grep -r >> DMCreateMatrix_DA_3d_MPIAIJ * >> dm/impls/da/fdda.c:extern PetscErrorCode >> DMCreateMatrix_DA_3d_MPIAIJ(DM,Mat); >> dm/impls/da/fdda.c:extern PetscErrorCode >> DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM,Mat); >> dm/impls/da/fdda.c: ierr = >> DMCreateMatrix_DA_3d_MPIAIJ_Fill(da,A);CHKERRQ(ierr); >> dm/impls/da/fdda.c: ierr = >> DMCreateMatrix_DA_3d_MPIAIJ(da,A);CHKERRQ(ierr); >> >> >> *dm/impls/da/fdda.c:#define __FUNCT__ >> "DMCreateMatrix_DA_3d_MPIAIJ"dm/impls/da/fdda.c:PetscErrorCode >> DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J)*dm/impls/da/fdda.c:#define >> __FUNCT__ "DMCreateMatrix_DA_3d_MPIAIJ_Fill" >> dm/impls/da/fdda.c:PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM >> da,Mat J) >> >> >> On 11 February 2016 at 04:08, Matthew Knepley <knep...@gmail.com> wrote: >> >>> On Wed, Feb 10, 2016 at 8:59 PM, Sang pham van <pvsang...@gmail.com> >>> wrote: >>> >>>> The irregular rows is quite many. The matrix really needs to be >>>> preallocated. >>>> Could you show me how to use DMCreateMatrix_DA_3d_MPIAIJ() directly? >>>> >>>> Just put the declaration right into your source. >>> >>> Matt >>> >>>> Pham >>>> On Feb 11, 2016 9:52 AM, "Matthew Knepley" <knep...@gmail.com> wrote: >>>> >>>>> On Wed, Feb 10, 2016 at 8:44 PM, Sang pham van <pvsang...@gmail.com> >>>>> wrote: >>>>> >>>>>> That is because my matrix has some rows which need more entries than >>>>>> usual. >>>>>> >>>>> >>>>> If its only a few, you could just turn off the allocation error. >>>>> >>>>> >>>>>> Where can i find source code of DMCreateMatrix()? >>>>>> >>>>>> >>>>> >>>>> https://bitbucket.org/petsc/petsc/src/827b69d6bb12709ff9b9a0dede31640477fc2b74/src/dm/impls/da/fdda.c?at=master&fileviewer=file-view-default#fdda.c-1024 >>>>> >>>>> Matt >>>>> >>>>>> Pham. >>>>>> On Feb 11, 2016 8:35 AM, "Matthew Knepley" <knep...@gmail.com> wrote: >>>>>> >>>>>>> On Wed, Feb 10, 2016 at 6:14 PM, Sang pham van <pvsang...@gmail.com> >>>>>>> wrote: >>>>>>> >>>>>>>> Hi, >>>>>>>> >>>>>>>> I am trying to create a DM matrix with DMCreateMatrix_DA_3d_MPIAIJ() >>>>>>>> instead of using DMCreateMatrix(). >>>>>>>> >>>>>>> >>>>>>> Why, that should be called automatically by DMCreateMatrix()? >>>>>>> >>>>>>> Matt >>>>>>> >>>>>>> >>>>>>>> Which header file should I include to use that routine? also, what >>>>>>>> is the source file containing the DMCreateMatrix() routine? >>>>>>>> >>>>>>>> Many thanks in advance. >>>>>>>> Pham >>>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> -- >>>>>>> What most experimenters take for granted before they begin their >>>>>>> experiments is infinitely more interesting than any results to which >>>>>>> their >>>>>>> experiments lead. >>>>>>> -- Norbert Wiener >>>>>>> >>>>>> >>>>> >>>>> >>>>> -- >>>>> What most experimenters take for granted before they begin their >>>>> experiments is infinitely more interesting than any results to which their >>>>> experiments lead. >>>>> -- Norbert Wiener >>>>> >>>> >>> >>> >>> -- >>> What most experimenters take for granted before they begin their >>> experiments is infinitely more interesting than any results to which their >>> experiments lead. >>> -- Norbert Wiener >>> >> >> >