[petsc-users] Problem with Preallocating

[Florian Lindner]

Am Mittwoch, 29. Oktober 2014, 13:00:17 schrieb Barry Smith:
> 
>    MatAssembly flushes out any preallocated space that was not used. Hence when you try to put a diagonal entry in later it is as if you did not preallocate. You should set zeros in the diagonal as part of your initial setting values in the matrix before calling MatAssembly.

Even with FLUSH_ASSEMBLY? I had this idea too, but was not able to reproduce the error I have in my application in a small example. 

Why does this program works? Shouldn't the MatDiagonalSet produce an error/warning that there were mallocs involved? Since the preallocation was lost?

Thanks,
Florian

#include <iostream>
#include "petscmat.h"   
#include "petscviewer.h"
#include "petscksp.h"

using namespace std;

// Compiling with: mpic++ -g3 -Wall -I ~/software/petsc/include -I ~/software/petsc/arch-linux2-c-debug/include -L ~/software/petsc/arch-linux2-c-debug/lib -lpetsc test.cpp

int main(int argc, char **args)
{
  PetscInitialize(&argc, &args, "", NULL);

  PetscErrorCode ierr = 0;
  int N = 9;
  Mat matrix;

  ierr = MatCreate(PETSC_COMM_SELF, &matrix);
  ierr = MatSetType(matrix, MATSBAIJ); CHKERRQ(ierr);
  ierr = MatSetSizes(matrix, PETSC_DECIDE, PETSC_DECIDE, N, N); CHKERRQ(ierr); 
  ierr = MatSetOption(matrix, MAT_SYMMETRY_ETERNAL, PETSC_TRUE); CHKERRQ(ierr);
  
  PetscInt nnz[N];
  nnz[0] = 2; nnz[1] = 1; nnz[2] = 1;
  nnz[3] = 1; nnz[4] = 1; nnz[5] = 1;
  nnz[6] = 1; nnz[7] = 1; nnz[8] = 1;
  
  for (int r = 0; r < N; r++) {
    cout << "Preallocated row " << r << " with " << nnz[r]  << " elements." << endl;
  }
  MatSeqSBAIJSetPreallocation(matrix, 1, PETSC_DEFAULT, nnz);

  MatSetValue(matrix, 0, 5, 10, INSERT_VALUES);

  ierr = MatAssemblyBegin(matrix, MAT_FLUSH_ASSEMBLY); CHKERRQ(ierr);
  ierr = MatAssemblyEnd(matrix, MAT_FLUSH_ASSEMBLY); CHKERRQ(ierr); 
  Vec zeros;
  MatGetVecs(matrix, NULL, &zeros);

  MatDiagonalSet(matrix, zeros, ADD_VALUES);
  ierr = MatAssemblyBegin(matrix, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ierr = MatAssemblyEnd(matrix, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 
  
  ierr = MatView(matrix, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

  MatDestroy(&matrix);
  PetscFinalize();  
  return 0;
}




> 
>   Barry
> 
> > On Oct 29, 2014, at 8:00 AM, Florian Lindner <mailinglists at xgm.de> wrote:
> > 
> > Hello,
> > 
> > I try to preallocate a sparse matrix like it was recommended in another posting, but get an error which kind of surprises me. Somehow I think it might be related to the order of assembly calls...
> > 
> > My code creates the matrix:
> > 
> >  MatSetType(_matrixC.matrix, MATSBAIJ);
> >  MatSetSizes(_matrixC.matrix, PETSC_DECIDE, PETSC_DECIDE, n, n); 
> >  MatSetOption(_matrixC.matrix, MAT_SYMMETRY_ETERNAL, PETSC_TRUE); 
> > 
> > Than I enter a loop to get the number of elements per row.  The entire code is a bit too long, but right before the SetPreallocaiton I do:
> > 
> >  for (int r = 0; r < n; r++) {
> >    cout << "Preallocated row " << r << " with " << nnz[r]  << " elements." << endl;
> >  }
> > 
> > which results in the output:
> > 
> > Matrix Size = 9
> > Preallocated row 0 with 9 elements.
> > Preallocated row 1 with 8 elements.
> > Preallocated row 2 with 7 elements.
> > Preallocated row 3 with 6 elements.
> > Preallocated row 4 with 5 elements.
> > Preallocated row 5 with 4 elements.
> > Preallocated row 6 with 1 elements.
> > Preallocated row 7 with 1 elements.
> > Preallocated row 8 with 1 elements.
> > 
> > at least in rows 0 to 5 it's a dense upper triangular matrix (incl. main diagonal).
> > 
> >  ierr = MatSeqSBAIJSetPreallocation(_matrixC.matrix, 1, PETSC_DEFAULT, nnz);
> > 
> > Now more or less the same loops as above with MatSetValues call for each row. 
> > 
> > Since the KSPSolver requires that the diagonal is set, even if it's merely with zeros I do:
> > 
> >  _matrixC.assemble(MAT_FLUSH_ASSEMBLY)  // my helper function, but I'm sure you get the point ;-)
> >  petsc::Vector zeros(_matrixC);
> >  MatDiagonalSet(_matrixC.matrix, zeros.vector, ADD_VALUES);
> >  ierr = MatAssemblyBegin(_matrixC.matrix, MAT_FINAL_ASSEMBLY); CHKERRV(ierr); 
> >  [ some more code ]
> >  ierr = MatAssemblyEnd(_matrixC.matrix, MAT_FINAL_ASSEMBLY); CHKERRV(ierr); 
> > 
> > But when I run the code I get an error message:
> > 
> > [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
> > [0]PETSC ERROR: Argument out of range
> > [0]PETSC ERROR: New nonzero at (0,0) caused a malloc
> > Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off this check
> > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
> > [0]PETSC ERROR: Petsc Release Version 3.5.2, unknown 
> > [0]PETSC ERROR: ./petBench on a arch-linux2-c-debug named helium by lindnefn Wed Oct 29 13:48:49 2014
> > [0]PETSC ERROR: Configure options --with-debugging=1
> > [0]PETSC ERROR: #1 MatSetValues_SeqSBAIJ() line 977 in /data2/scratch/lindner/petsc/src/mat/impls/sbaij/seq/sbaij.c
> > [0]PETSC ERROR: #2 MatSetValues() line 1135 in /data2/scratch/lindner/petsc/src/mat/interface/matrix.c
> > [0]PETSC ERROR: #3 MatDiagonalSet_Default() line 203 in /data2/scratch/lindner/petsc/src/mat/utils/axpy.c
> > [0]PETSC ERROR: #4 MatDiagonalSet() line 241 in /data2/scratch/lindner/petsc/src/mat/utils/axpy.c
> > 
> > I'm surprised that I get this error message since I do have preallocated 9 elements and I do not understand why inserting an element at (0, 0) could be a problem / needing a malloc. Has the preallocation got lost somewhere?
> > 
> > Thanks once again....
> > 
> > Florian
>