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

[Appel, Thibaut]

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? Is there something wrong with my understanding of ADD_VALUES?

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. Is it better to use MatXAIJSetPreallocation?

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. 

Kind regards,


Thibaut