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.