Problem with rectangular matrices

[Matthew Knepley]

On Tue, May 26, 2009 at 12:07 PM, Lisandro Dalcin <dalcinl at gmail.com> wrote:

> The only real issue I see here is that the term "diagonal block" could
> lead to misinterpretation. But the behaviour you commented is the one
> I would expect.


Okay, anyone else not want me to try and do this?


>
> The rest of this mail is perhaps off-topic, but it is somewhat related.
>
> The other issues is that some Mat preallocation routines are lazy
> about properly setting the cmap. I had to implement some hackery in
> petsc4py to make sure that created Mat's (at least these created from
> petsc4py calls) have a cmap consistent with the local/global
> row/column sizes.
>
> Moreover, supose you create a global square matrix with the following
> row/col distribution (I know, it is perhaps a nonsense, but IMHO is a
> valid distribution)
> from petsc4py import PETSc
>
> size = PETSc.COMM_WORLD.getSize()
> assert size == 2
> rank = PETSc.COMM_WORLD.getRank()
> Ms = [2, 6]
> Ns = [6, 2]
>
> m = Ms[rank]
> n = Ns[rank]
>
> A = PETSc.Mat().createAIJ([(m, None), (n, None)]) # None is
> intrepreted as PETSc.DECIDE
> for i in range(8):
>    A[i,i] = 1.0/(i+1)
> A.assemble()
>
> x, b = A.getVecs()
> x.set(1)
> A.mult(x,b)
>
> Up to now, all works as expected, all rmap/cmap are consistent, the
> 'x' and 'b' vecs have the same global size, though different local
> sizes, and the MatMult succeed. All fine... But now, suppose you do:
>
> ksp = PETSc.KSP().create()
> ksp.setOperators(A)
> ksp.pc.type = 'none'
> ksp.solve(b,x)
>
> You are lost, this does not work (even if you try with nonzero initial
> guess)... But it should, right?
>

No, I  do not think that this should work, because the Krylov methods impose
a more stringent condition on the matrix A, namely that the layout of Ax
should
be the same as the layout of x. This necessitates a symmetric partition, and
is
probably at the root of our problems with rectangular matrices.

  Matt


> On Tue, May 26, 2009 at 10:35 AM, Matthew Knepley <knepley at gmail.com>
> wrote:
> > I think there is an inconsistency in our definitions. If you create a
> matrix
> > which
> > is 4x8 on 2 procs, then you end up with a row and column map like this:
> >
> >     rowMap           colMap
> >   [0..2) [2..4)      [0..4) [4..8)
> >
> > However, we use the column map for queries which then give nonsense. The
> > most
> > damaging is in MatSetValues(). If I give a value for (2,2) on proc 1,
> then
> > since
> > 2 < cstart = 4, this value is considered "off-diagonal" and can cause a
> > malloc(). So
> > I think we have a choice. We can define "diagonal block" to have this
> weird
> > shape,
> > or we can change all queries for the diagonal block to use only the
> rowMap.
> > I vote
> > for the latter, but does anyone have an objection?
> >
> >    Matt
> >
> > --
> > 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
> >
>
>
>
> --
> Lisandro Dalcín
> ---------------
> Centro Internacional de Métodos Computacionales en Ingeniería (CIMEC)
> Instituto de Desarrollo Tecnológico para la Industria Química (INTEC)
> Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)
> PTLC - Güemes 3450, (3000) Santa Fe, Argentina
> Tel/Fax: +54-(0)342-451.1594
>



-- 
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-dev/attachments/20090526/259f047b/attachment.html>