[petsc-users] MatSetValues() for MATMPICUSP

[recrusader]

Dear Jed,

When I removed 'if (NONEW == -2)
SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"New nonzero at
(%D,%D) caused a malloc",ROW,COL); \' from the following function in
src/mat/impls/aij/seq/aij.h. It works in GPU mode. Do you have any
comments? Thanks a lot.

#define MatSeqXAIJReallocateAIJ(Amat,AM,BS2,NROW,ROW,COL,RMAX,AA,AI,AJ,RP,AP,AIMAX,NONEW,datatype)
\
  if (NROW >= RMAX) {\
        Mat_SeqAIJ *Ain = (Mat_SeqAIJ*)Amat->data;\
        /* there is no extra room in row, therefore enlarge */ \
        PetscInt   CHUNKSIZE = 15,new_nz = AI[AM] +
CHUNKSIZE,len,*new_i=0,*new_j=0; \
        datatype   *new_a; \
 \
        /*if (NONEW == -2)
SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"New nonzero at
(%D,%D) caused a malloc",ROW,COL); \  */ \
        /* malloc new storage space */ \
        ierr = PetscMalloc3(BS2*new_nz,datatype,&new_a,new_nz,PetscInt,&new_j,AM+1,PetscInt,&new_i);CHKERRQ(ierr);\
 \
        /* copy over old data into new slots */ \
        for (ii=0; ii<ROW+1; ii++) {new_i[ii] = AI[ii];} \
        for (ii=ROW+1; ii<AM+1; ii++) {new_i[ii] = AI[ii]+CHUNKSIZE;} \
        ierr = PetscMemcpy(new_j,AJ,(AI[ROW]+NROW)*sizeof(PetscInt));CHKERRQ(ierr);
\
        len = (new_nz - CHUNKSIZE - AI[ROW] - NROW); \
        ierr = PetscMemcpy(new_j+AI[ROW]+NROW+CHUNKSIZE,AJ+AI[ROW]+NROW,len*sizeof(PetscInt));CHKERRQ(ierr);
\
        ierr = PetscMemcpy(new_a,AA,BS2*(AI[ROW]+NROW)*sizeof(datatype));CHKERRQ(ierr);
\
        ierr = PetscMemzero(new_a+BS2*(AI[ROW]+NROW),BS2*CHUNKSIZE*sizeof(datatype));CHKERRQ(ierr);\
        ierr = PetscMemcpy(new_a+BS2*(AI[ROW]+NROW+CHUNKSIZE),AA+BS2*(AI[ROW]+NROW),BS2*len*sizeof(datatype));CHKERRQ(ierr);
 \
        /* free up old matrix storage */ \
        ierr = MatSeqXAIJFreeAIJ(A,&Ain->a,&Ain->j,&Ain->i);CHKERRQ(ierr);\
        AA = new_a; \
        Ain->a = (MatScalar*) new_a;               \
        AI = Ain->i = new_i; AJ = Ain->j = new_j;  \
        Ain->singlemalloc = PETSC_TRUE; \
 \
        RP          = AJ + AI[ROW]; AP = AA + BS2*AI[ROW]; \
        RMAX        = AIMAX[ROW] = AIMAX[ROW] + CHUNKSIZE; \
        Ain->maxnz += BS2*CHUNKSIZE; \
        Ain->reallocs++; \
      } \

Best,
Yujie


On 2/11/12, Jed Brown <jedbrown at mcs.anl.gov> wrote:
> On Sat, Feb 11, 2012 at 00:00, recrusader <recrusader at gmail.com> wrote:
>
>> in
>> http://www.mcs.anl.gov/petsc/petsc-dev/src/ksp/ksp/examples/tutorials/ex43.c.html
>>
>> Why are two matrices set to MATAIJ? if I set MATAIJCUSP, is them changed?
>>
>> 1468:
>> DMCreateMatrix<http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCreateMatrix.html#DMCreateMatrix>
>> (da_Stokes,MATAIJ<http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/Mat/MATAIJ.html#MATAIJ>
>> ,&A);
>>
>> 1469:   DMCreateMatrix
>> <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCreateMatrix.html#DMCreateMatrix>(da_Stokes,MATAIJ
>> <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/Mat/MATAIJ.html#MATAIJ>,&B);
>>
>>
> You have to follow the code down to where those matrices are used.
>
> 1478:   AssembleA_Stokes(A,da_Stokes,da_prop,properties);
> 1479:   AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
>
> This uses a "block preconditioner", the B matrix is intentionally different
> from the A matrix in order to handle the saddle point.
>