Jed,
You pose several issues here.
1) By default SBAIJ matrices generate an error if you set a value in the lower
triangular portion.
I think this is the wrong default. I've changed it to ignore lower
triangular values unless you explicitly set the option.
2) MatSetType(A,MATSBAIJ); and friends deletes the matrix and creates a new one
if the matrix is MATSEQSBAIJ or MATMPISBAIJ.
This is silly. Given our model of setting types however there is no
trivial way to fix this. We could try to handle setting a heirarchy of types
properly (how to know if something is a subclass? etc?) or do something just
for Mat to handle the
rank = 1 and rank > 1 cases. I'm inclined to add a special Mat register
for example MatRegisterBaseName("SBAIJ","SeqSBAIJ","MPISBAIJ"); and have
MatSetType() first check if the type is in the BaseName list and if so grab the
correct
classname based on rank, then it will realize the matrix is not having a
type change and will not destroy it. Note that this means routines like
MatCreateSBAIJ() will disappear.
3) The general wordyness of setting up the matrix fill from the DM.
My planned model for the future was to change DMGetMatrix() to
DMFillMatrix() and have code something like
MatCreate(comm,&Ap);
MatSetOptionsPrefix(Ap,"Ap_");
MatSetFromOptions(Ap);
MatSetOption(Ap,MAT_SYMMETRIC,PETSC_TRUE);
DMFillMatrix(fsu,Ap);
If is actually not any shorter but I think it is a better model.
If no one objects I will implement 2).
3) I think needs more discussion before implementation.
Barry
On Apr 30, 2011, at 6:22 PM, Jed Brown wrote:
> I was tired of writing
> MatSetOption(Ap,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE) into my code to
> find preallocation issues, so I added -mat_new_nonzero_allocation_err to
> MatSetFromOptions(). Now I actually call that function. In particular, I do
> things like
>
> PetscOptionsList("-Ap_mat_type","Matrix type for velocity
> operator","",MatList,mattype_Ap,mattype_Ap,sizeof(mattype_Ap),NULL);
> ...
> DMGetMatrix(fsu,mattype_Ap,&Ap);
> MatSetOptionsPrefix(Ap,"Ap_");
> MatSetOption(Ap,MAT_SYMMETRIC,PETSC_TRUE);
> MatSetFromOptions(Ap);
> // Check type and set MAT_IGNORE_LOWER_TRIANGULAR if necessary
>
>
> Three problems:
>
> 1. This is too damn much work.
>
> 2. Checking the type is an ugly hack. Should MAT_IGNORE_LOWER_TRIANGULAR be
> default, or should nonsymmetric matrices ignore it, or should there be a new
> option with those semantics?
>
> 3. It's wrong to use "-Ap_mat_type sbaij" because the actual type becomes
> seqsbaij, and then preallocation is thrown away [*] when MatSetFromOptions
> finds "sbaij", calls MatSetType which compares it to "seqsbaij", and resets
> the type. This means that I have to use a different prefix to choose which
> type to pass to DMGetMatrix() than the prefix I subsequently set on the
> matrix.
>
>
> Any ideas how to make this work in a sane way?
>
> [*] I was getting a seg-fault, but I set mat->preallocated to false in
> MatSetType. It doesn't crash any more, but preallocation is still thrown away.
