[petsc-users] Create DM matrix
Dave May
dave.mayhem23 at gmail.com
Thu Feb 11 02:29:34 CST 2016
On 11 February 2016 at 09:16, Sang pham van <pvsang002 at 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.mayhem23 at 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 at 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 <knepley at gmail.com> wrote:
>>
>>> On Wed, Feb 10, 2016 at 8:59 PM, Sang pham van <pvsang002 at 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" <knepley at gmail.com> wrote:
>>>>
>>>>> On Wed, Feb 10, 2016 at 8:44 PM, Sang pham van <pvsang002 at 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" <knepley at gmail.com> wrote:
>>>>>>
>>>>>>> On Wed, Feb 10, 2016 at 6:14 PM, Sang pham van <pvsang002 at 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
>>>
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160211/d33e6163/attachment-0001.html>
More information about the petsc-users
mailing list