On Thu, Oct 30, 2014 at 7:25 AM, Barry Smith <bsm...@mcs.anl.gov> wrote:
>
> > On Oct 30, 2014, at 3:36 AM, Florian Lindner <mailingli...@xgm.de>
> wrote:
> >
> > 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?
>
> No, flush assembly does not clean out any unused preallocated space.
> final assembly does
>
Something else must be going on in your full code (as shown by the small
example). Are you running it
in serial or parallel?
Thanks,
Matt
> Barry
>
> > 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 <mailingli...@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
> >>
>
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
