Problem with rectangular matrices

[Lisandro Dalcin]

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?


>>
>> 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