[petsc-users] Create DM matrix

Sang pham van pvsang002 at gmail.com
Thu Feb 11 02:16:32 CST 2016


I insert the lines below into my code, but it does not work:
  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,&ltog);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/f6cd0d79/attachment.html>


More information about the petsc-users mailing list