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