[petsc-users] Malloc error with 'correct' preallocation?

[Smith, Barry F.]

> On Feb 27, 2018, at 7:05 PM, Appel, Thibaut <t.appel17 at imperial.ac.uk> wrote:
> 
> Dear PETSc developers and users,
> 
> I am forming a sparse matrix in complex, double precision arithmetic and can’t understand why I have «PETSC ERROR: Argument out of range - New nonzero at (X,X) caused a malloc» during the assembly of the matrix.
> 
> This matrix discretizes a PDE in 2D using a finite difference method in both spatial directions and, in short, here is the ensemble of routines I call:
> 
>    CALL MatCreate(PETSC_COMM_WORLD,MatA,ierr)
>    CALL MatSetType(MatA,MATAIJ,ierr)
>    CALL MatSetSizes(MatA,PETSC_DECIDE,PETSC_DECIDE,leading_dimension,leading_dimension,ierr)
>    CALL MatSeqAIJSetPreallocation(MatA,0,PETSC_NULL_INTEGER,ierr)
>    CALL MatMPIAIJSetPreallocation(MatA,0,PETSC_NULL_INTEGER,0,PETSC_NULL_INTEGER,ierr)
>    CALL MatGetOwnershipRange(MatA,istart,iend,ierr)
> 
> Then I preallocate my matrix doing tests on the derivative coefficient arrays like
> 
>    ALLOCATE(nnz(istart:iend-1), dnz(istart:iend-1), onz(istart:iend-1))
>    nnz = 0
>    dnz = 0
>    onz = 0
>    DO row= istart, iend-1
> 
> 	*detect the nonzero elements*
>        IF (ABS(this_derivative_coef) > 0.D0) THEN 
>          nnz(row) = nnz(row) + corresponding_stencil_size
>          DO *elements in stencil*
>              *compute corresponding column*
>              IF ((col >= istart) .AND. (col <= (iend-1))) THEN
>                dnz(row) = dnz(row) + 1
>              ELSE
>                onz(row) = onz(row) + 1
>              END IF
>          END DO
> 	END IF
> 
>   END DO
>   CALL MatSeqAIJSetPreallocation(MatA,PETSC_DEFAULT_INTEGER,nnz,ierr)
>   CALL MatMPIAIJSetPreallocation(MatA,PETSC_DEFAULT_INTEGER,dnz,PETSC_DEFAULT_INTEGER,onz,ierr)
>   CALL MatSetOption(MatA,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE,ierr)
> 
> And assemble the matrix, at each row, with the different derivatives terms (pure x derivative, pure y derivative, cross xy derivative…)
> 
>    DO row = istart, iend-1
> 
>        cols(0) = *compute corresponding column*
>        vals(0) = no_derivative_coef
>        CALL MatSetValues(MatA,1,row,1,cols,vals,ADD_VALUES,ierr)
> 
>        DO m=0,x_order
>          cols(m) = *compute corresponding column*
>          vals(m) = x_derivative_coef
>        END DO
>        CALL MatSetValues(MatA,1,row,x_order+1,cols,vals,ADD_VALUES,ierr)
> 
>        DO m=0,y_order
>          cols(m) = *compute corresponding column*
>          vals(m) = y_derivative_coef
>        END DO
>        CALL MatSetValues(MatA,1,row,y_order+1,cols,vals,ADD_VALUES,ierr)
> 
>        DO n=0,y_order
>          DO m=0,x_order
>            cols(..) =  *compute corresponding column*
>            vals(..)= xy_derivative_coef
>          END DO
>        END DO
>        CALL MatSetValues(MatA,1,row,(x_order+1)*(y_order+1),cols,vals,ADD_VALUES,ierr)
> 
>    END DO
> 
>    CALL MatAssemblyBegin(MatA,MAT_FINAL_ASSEMBLY,ierr)
>    CALL MatAssemblyEnd(MatA,MAT_FINAL_ASSEMBLY,ierr)
> 
> I am using ADD_VALUES as the different loops here-above can contribute to the same column.
> 
> The approach I chose is therefore preallocating without overestimating the non-zero elements and hope that the MAT_IGNORE_ZERO_ENTRIES option discards the 'vals(…)' who are exactly zero during the assembly (I read that the criteria PETSc uses to do so is ‘== 0.0’) so with the test I used everything should work fine.
> 
> However, when testing with 1 MPI process, I have this malloc problem appearing at a certain row.
> 
> I print the corresponding nnz(row) I allocated for this row, say NZ.
> I print for each (cols, vals) computed in the loops above the test (vals /= zero) to number the non-zero elements and notice that the malloc error appears when I insert a block of non-zero elements in which the last one is the NZth!
> 
> In other words, why a malloc if the 'correct' number of elements is allocated?

   Don't know

> Is there something wrong with my understanding of ADD_VALUES?

  Unlikely

> 
> I read somewhere that it is good to call both MatSeqAIJSetPreallocation and MatMPIAIJSetPreallocation as PETSc will automatically call the right one whether it is a serial or parallel computation.

   Yes this is true

> Is it better to use MatXAIJSetPreallocation?

   In your case no

> 
> Thank you in advance for any advice that could put me on the trail to correctness, and I would appreciate any correction should I do something that looks silly. 

  Nothing looks silly. 

   The easiest approach to tracking down the problem is to run in the debugger and put break points at the key points and make sure it behaves as expected at those points.

   Could you send the code? Preferably small and easy to build. If I can reproduce the problem I can track down what is going on.

   Barry

> 
> Kind regards,
> 
> 
> Thibaut
>