Problem with rectangular matrices

[Matthew Knepley]

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

> On Tue, May 26, 2009 at 2:15 PM, Matthew Knepley <knepley at gmail.com>
> wrote:
> > 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?
> >
>
> Could you elaborate a bit more what exactly are you going to try? Are
> you proposing that rectangular matrices should always have a square
> diagonal block?
>

Yes, I want to go into the code and every place we query the diagonal block,
replace
cMap by rMap. This is the only thing that makes sense to me. By my
definitions, diagonal
blocks are disjoint, so they must be square. Since we normally own rows, it
makes sense
to use row blocking to define them.

   Matt


> >>
> >> 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
> >
>
>
>
> --
> 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/3e430e64/attachment.html>