[petsc-users] Creating SBAIJ matrix using 1-based csr arrays ia, ja and acsr

[Matthew Knepley]

On Tue, Nov 25, 2014 at 12:54 PM, Kirill Voronin <kvoronin at labchem.sscc.ru>
wrote:

>
> Thank you for the prompt answer!
>
> I shifted the indices but when calling MatSeqSBAIJSetPreallocationCSR got
> the error - nnz cannot be greater than block row length: local row 0 value
> 8 rowlength 1.
>

It sounds like the input you are giving is for actual rows, but this call
expects the
CSR to reference blocks. Does that make sense?


> Can anyone please make a bit more clear some details of block AIJ format
> in PETSC?
>
> My code looks like:
> -----------------------------------------------------
>   ierr = MatCreate(comm,&Afromcsr);CHKERRQ(ierr);
>   ierr = MatSetType(Afromcsr,"sbaij"); CHKERRQ(ierr);
>   ierr = MatSetSizes(Afromcsr,dim,dim,dim,dim);CHKERRQ(ierr);
>   ierr = MatSetFromOptions(Afromcsr);CHKERRQ(ierr);
>
>   printf("fill in for values of acsr \n");
>
>   for (i = 0; i <= dim; i++)
>         ia[i] -= 1;
>   for (i = 0; i < ia[dim]; i++)
>         ja[i] -= 1;
>   ierr = MatSeqSBAIJSetPreallocationCSR(Afromcsr, dim, ia, ja, acsr);
> CHKERRQ(ierr);
> -----------------------------------------------------------
> The questions:
>
> 1) Why no MatMPISBAIJSetPreaalocationCSR? Or should
> MatCreateMPISBAIJWithArrays be used instead?
>

Not written yet.


> 2) What is block size and why it is assumed to be 1 (according to the
> error I got).
>

Block size is the size of the dense blocks out of which the matrix is
constructed.
This can achieve better performance since indexing is only done on blocks,
not
individual elements, and dense algebra can be done on block-block
computation.


> 3) What is a local row when I call MatSeqSBAIJSetPreallocationCSR?
>

Its a row of blocks.


> 4) Can I consider my matrix as a single block when calling
> MatSeqSBAIJSetPreallocationCSR ?
>

The maximum block size is 13 I think.


> 5) Can I put my acsr array (with values of nonzero elements of the matrix)
> as the argument v[] in the MatSeqSBAIJSetPreallocationCSR? The manual
> pages told me that v[] should be of the form v[nnz][bs][bs] and I don't
> get why not only v[nnz]?
>

There is no v[] argument:


http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSeqSBAIJSetPreallocation.html

  Thanks,

     Matt


> Sorry for that set of questions but I missed an example of usage for a
> symmetric CSR matrix.
>
> Maybe a few lines of code for getting a symmetric (with upper part only
> stored) matrix from ia, ja and a arrays (usual CSR format) would help me
> to get the whole idea of PETSC matrix creating.
>
> Thanks in advance!
>
>
> Best regards,
>
> Kirill Voronin
>
>
> > On Tue, Nov 25, 2014 at 9:07 AM, Kirill Voronin
> > <kvoronin at labchem.sscc.ru>
> > wrote:
> >
> >
> >>
> >> Hello!
> >>
> >>
> >> I am trying to create a matrix for solving Ax = b with PETSC with given
> >>  1-based(!) indexed ia, ja and acsr arrays for a symmetric (only upper
> >> part is stored) complex matrix.
> >>
> >> Without memory preallocation it takes too long to fill in all the
> >> values of acsr into my Mat A object.
> >>
> >> But I failed to call
> >> ierr = MatSeqSBAIJSetPreallocationCSR(Afromcsr, dim, iaHSS, jaHSS,
> >> acHSS); because of the 1 based indexing. And I didn't find an option how
> >> to handle this issue.
> >>
> >> From the manual:
> >> "By default the internal data representation for the AIJ formats employs
> >>  zero-based indexing. For compatibility with standard Fortran storage,
> >> thus enabling use of external Fortran software packages such as
> SPARSKIT,
> >> the option -mat_aij_oneindex enables one-based indexing, where the
> >> stored row and column indices begin at one, not zero. All user calls to
> >> PETSc routines,
> >> regardless of this option, use zero-based indexing."
> >>
> >> So, what can you recommend to do? I obviously don't want to double the
> >> memory for storing ia and ja arrays in 0-based indexing. I don't want to
> >>  change the indexing everywhere to 0-based since I use a piece of
> >> external code which works only for 1-based indexing as a preconditioner
> >> in PETSC KSP solver.
> >>
> >>
> >
> > Just decrement ja, make the call, and increment ja.
> >
> >
> > Thanks,
> >
> >
> > Matt
> >
> >
> >
> >> Thanks in advance!
> >>
> >>
> >> With best regards,
> >>
> >>
> >> Kirill Voronin
> >>
> >>
> >>
> >>
> >>
> >
> >
> > --
> > 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
> >
> >
>
>
> --
>
>
>


-- 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20141125/fdbf2228/attachment.html>