On Sep 14, 2019, at 7:19 AM, "Smith, Barry F." <bsm...@mcs.anl.gov> wrote:
>
>
>> On Sep 13, 2019, at 11:19 PM, Pierre Jolivet <pierre.joli...@enseeiht.fr>
>> wrote:
>>
>>
>>
>>> On 11 Sep 2019, at 3:46 PM, Smith, Barry F. <bsm...@mcs.anl.gov> wrote:
>>>
>>>
>>> The phrase Hermitian symmetric was originally an attempt to make the text
>>> clearer, but in fact it is confusing and wrong. The phrase Hermitian
>>> symmetric should be removed and just Hermitian used.
>>>
>>> Real matrices are either either both symmetric and Hermitian or not
>>> symmetric and not Hermitian
>>>
>>> Complex matrices can be set to symmetric or Hermitian. For SBAIJ this
>>> changes the meaning of the values in the matrix; for matrices with full
>>> storage this means the user is stating that the matrix values are actually
>>> symmetric or Hermitian in the storage.
>>>
>>> Currently for complex matrices one can state the matrix is both symmetric
>>> and Hermitian which would mean one is stating the imaginary parts are all
>>> zero; if this is the case then it doesn't matter if the
>>> MatMult_SeqSBAIJ_1_Hermitian() or the MatMult_SeqSBAIJ_1() is used for
>>> SeqBAIJ. Note that calling both of these doesn't mean the imaginary parts
>>> will be zeroed out, it is statement by the user that the imaginary parts
>>> are zero.
>>>
>>> For parallel matrices, even in debug node, checking the user is telling
>>> the truth when they set symmetric or Hermitian is considered too expensive
>>> so we just trust the user.
>>>
>>> With the replacement of Hermitian symmetric by Hermitian will the text
>>> now be clear or should more clarification be given in the manual pages (as
>>> above)?
>>
>> I think it makes the text more clear.
>> However, it’s not very clear to me:
>> 1) how to setup a SeqSBAIJ Mat as Hermitian, i.e., is MatSetOption(A,
>> MAT_SYMMETRIC, PETSC_FALSE); MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE); OK?
>
> That would work or I would just say the next overwrites the previous, so
> just call MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE); and it flips the
> symmetric flag also.
> Similar for setting symmetric.
So, how do I let PETSc and SLEPc know that my matrix is both symmetric _and_
Hermitian if you flip the flag internally?
> There should presumably really be only one flag inside the Mat object for
> tracking this something like enum {MAT_FORM_HERMITIAN, MAT_FORM_SYMMETRIC,
> MAT_FORM_GENERAL} MatSymmetricForm; and for SBAIJ it could only have the
> values MAT_FORM_HERMITIAN, MAT_FORM_SYMMETRIC
>
>> 2) why isn’t there more checks in MatSetOption_SeqSBAIJ.
>
> Oh, because hardly anyone ever uses it for complex numbers and people
> probably never have some symmetric matrices in the simulation and some
> Hermitian. My guess is 95% of all cases are one or the other only.
>
> Yes, definitely there should be a bunch more error checking.
>
>> As we can all agree one, if(MAT_SYMMETRIC and MAT_HERMITIAN), then there
>> should be no error for bs > 1, but on the contrary, isn’t (!MAT_SYMMETRIC
>> and MAT_HERMITIAN and bs > 1) problematic right now?
>
> Why does it matter what bs is?
I don't know, the people who put the SETERRQ in MatSetOption_SeqSBAIJ may be
able to answer that question.
Presumably because there are some BLAS/LAPACK calls with "T" instead of "H"?
Thanks,
Pierre
> Are you thinking of the diagonal blocks?
>
>
>> 3) why any of those checks should be performed in MatSetOption_SeqSBAIJ.
>> Calling MatAXPY with two matrices of the same type with the same flags is
>> feasible no matter the value of bs.
>
> True, but you would need to check the flags when it is called.
>>
>> Thanks,
>> Pierre
>>
>>> Barry
>>>
>>>
>>> case MAT_HERMITIAN:
>>> #if defined(PETSC_USE_COMPLEX) /* MAT_HERMITIAN is a synonym for
>>> MAT_SYMMETRIC with reals */
>>> if (!A->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must
>>> call MatAssemblyEnd() first");
>>> if (A->cmap->n < 65536 && A->cmap->bs == 1) {
>>> A->ops->mult = MatMult_SeqSBAIJ_1_Hermitian_ushort;
>>> } else if (A->cmap->bs == 1) {
>>> A->ops->mult = MatMult_SeqSBAIJ_1_Hermitian;
>>> } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for Hermitian
>>> with block size greater than 1");
>>> #endif
>>>
>>>
>>>
>>> case MAT_SYMMETRIC:
>>> mat->symmetric = flg;
>>> if (flg) mat->structurally_symmetric = PETSC_TRUE;
>>> mat->symmetric_set = PETSC_TRUE;
>>> mat->structurally_symmetric_set = flg;
>>> #if !defined(PETSC_USE_COMPLEX)
>>> mat->hermitian = flg;
>>> mat->hermitian_set = PETSC_TRUE;
>>> #endif
>>> break;
>>> case MAT_HERMITIAN:
>>> mat->hermitian = flg;
>>> if (flg) mat->structurally_symmetric = PETSC_TRUE;
>>> mat->hermitian_set = PETSC_TRUE;
>>> mat->structurally_symmetric_set = flg;
>>> #if !defined(PETSC_USE_COMPLEX)
>>> mat->symmetric = flg;
>>> mat->symmetric_set = PETSC_TRUE;
>>> #endif
>>>
>>>
>>>
>>>
>>>> On Sep 11, 2019, at 5:55 AM, Pierre Jolivet via petsc-dev
>>>> <petsc-dev@mcs.anl.gov> wrote:
>>>>
>>>> (Putting back petsc-dev on c/c)
>>>> I don’t think you are the one that needs to do the changes.
>>>> There should be clarifications from the PETSc side of things, first.
>>>> From the man page your are quoting
>>>> (https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MATSEQSBAIJ.html):
>>>> 1) how to make a SBAIJ Mat Hermitian but not symmetric? Is this handled by
>>>> PETSc? Is this legal to do MatSetOption(A, MAT_SYMMETRIC, PETSC_FALSE)
>>>> followed by MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE)? I think this is
>>>> the case where bs>1 is not handled internally by PETSc, so for this case
>>>> it should be OK to error out.
>>>> 2) “Hermitian symmetric” is a fancy way of saying that the matrix has no
>>>> imaginary part. Or is just “Hermitian" what is supposed to be written on
>>>> the page? In which case, how do you tell to PETSc that you are indeed
>>>> dealing with a matrix with no imaginary part.
>>>>
>>>> Thanks,
>>>> Pierre
>>>>
>>>>> On 11 Sep 2019, at 11:51 AM, Jose E. Roman <jro...@dsic.upv.es> wrote:
>>>>>
>>>>> Now I understand better. You are right. Maybe I should do some changes in
>>>>> SLEPc. I will think about this next week.
>>>>>
>>>>>
>>>>>> El 11 sept 2019, a las 11:38, Pierre Jolivet
>>>>>> <pierre.joli...@enseeiht.fr> escribió:
>>>>>>
>>>>>>
>>>>>>
>>>>>>> On 11 Sep 2019, at 11:21 AM, Jose E. Roman <jro...@dsic.upv.es> wrote:
>>>>>>>
>>>>>>> In
>>>>>>> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MATSEQSBAIJ.html
>>>>>>> it says "To make it Hermitian symmetric you can call MatSetOption(Mat,
>>>>>>> MAT_HERMITIAN);" So my understanding is that in SBAIJ assumes
>>>>>>> MAT_SYMMETRIC and the user must set MAT_HERMITIAN if desired.
>>>>>>
>>>>>> Right, that is exactly what you do in STMatSetHermitian. You get from
>>>>>> the user a SBAIJ matrix (so SYMMETRIC = true), and you set yourself
>>>>>> HERMITIAN = true.
>>>>>>
>>>>>>> I don't think it makes sense to set both MAT_HERMITIAN and
>>>>>>> MAT_SYMMETRIC.
>>>>>>
>>>>>> Why? A symmetric (or Hermitian) matrix with no imaginary part is both
>>>>>> Hermitian and symmetric, so at least this set of options defines
>>>>>> correctly the structure of the matrix (which is not the case if you only
>>>>>> assume it is _either_ symmetric or Hermitian).
>>>>>>
>>>>>>> But I may be wrong because I find this part of PETSc very confusing.
>>>>>>>
>>>>>>>
>>>>>>>> El 11 sept 2019, a las 11:08, Jose E. Roman <jro...@dsic.upv.es>
>>>>>>>> escribió:
>>>>>>>>
>>>>>>>> Do you mean setting both flags MAT_HERMITIAN and MAT_SYMMETRIC? Is
>>>>>>>> this possible?
>>>>>>
>>>>>> It was possible before my commit for bs == 1, but not for bs > 1.
>>>>>>
>>>>>> At least now my tests are not failing.
>>>>>>
>>>>>>>>> El 11 sept 2019, a las 11:04, Pierre Jolivet
>>>>>>>>> <pierre.joli...@enseeiht.fr> escribió:
>>>>>>>>>
>>>>>>>>> Symmetric <=> a_ij = a_ji
>>>>>>>>> Hermitian <=> a_ij = conj(a_ji)
>>>>>>>>>
>>>>>>>>> Symmetric + Hermitian <=> a_ij = conj(a_ji) = a_ji <=> imag(a_ji) = 0
>>>>>>>>> \forall ij
>>>>>>>>>
>>>>>>>>> Or am I stupid?
>>>>>>>>>
>>>>>>>>>> On 11 Sep 2019, at 10:59 AM, Jose E. Roman <jro...@dsic.upv.es>
>>>>>>>>>> wrote:
>>>>>>>>>>
>>>>>>>>>> Not sure if I understand you. Do you mean that a complex SBAIJ Mat
>>>>>>>>>> with MAT_HERMITIAN flag can be assumed to have zero imaginary part?
>>>>>>>>>> I don't think so. This matrix should have real diagonal entries, but
>>>>>>>>>> off-diagonal entries should be allowed to have nonzero imaginary
>>>>>>>>>> part. This is what is done in MatMult_SeqSBAIJ_1_Hermitian(), where
>>>>>>>>>> off-diagonal entries are conjugated when used for the strict lower
>>>>>>>>>> triangular part. So I guess the right fix is to implement
>>>>>>>>>> MatMult_SeqSBAIJ_2_Hermitian(), MatMult_SeqSBAIJ_3_Hermitian() and
>>>>>>>>>> so on with appropriate use of PetscConj().
>>>>>>>>>>
>>>>>>>>>> Jose
>>>>>>>>>>
>>>>>>>>>>> El 11 sept 2019, a las 10:36, Pierre Jolivet
>>>>>>>>>>> <pierre.joli...@enseeiht.fr> escribió:
>>>>>>>>>>>
>>>>>>>>>>> Nevermind, this is the wrong fix.
>>>>>>>>>>> The proper fix is in PETSc. It should not error out if the matrix
>>>>>>>>>>> is also symmetric.
>>>>>>>>>>> Indeed, complex symmetric Hermitian => complex with no imaginary
>>>>>>>>>>> part.
>>>>>>>>>>> Thus all operations like MatMult, MatMultHermitianTranspose,
>>>>>>>>>>> Cholesky… will work for bs > 1, since all is filled with zeroes.
>>>>>>>>>>> I will take care of this, I’m c/c’ing petsc-dev so that they don’t
>>>>>>>>>>> have to “reverse engineer” the trivial change to
>>>>>>>>>>> MatSetOption_SeqSBAIJ.
>>>>>>>>>>>
>>>>>>>>>>> Sorry about the noise.
>>>>>>>>>>>
>>>>>>>>>>> Thank you,
>>>>>>>>>>> Pierre
>>>>>>>>>>>
>>>>>>>>>>>> On 10 Sep 2019, at 8:37 AM, Pierre Jolivet
>>>>>>>>>>>> <pierre.joli...@enseeiht.fr> wrote:
>>>>>>>>>>>>
>>>>>>>>>>>> Hello,
>>>>>>>>>>>> Could you consider not setting MAT_HERMITIAN here
>>>>>>>>>>>> http://slepc.upv.es/documentation/current/src/sys/classes/st/interface/stsles.c.html#line276
>>>>>>>>>>>> when using SBAIJ matrices with bs > 1?
>>>>>>>>>>>> This makes PETSc error out with
>>>>>>>>>>>> # [1]PETSC ERROR: No support for this operation for this object
>>>>>>>>>>>> type
>>>>>>>>>>>> # [1]PETSC ERROR: No support for Hermitian with block size
>>>>>>>>>>>> greater than 1
>>>>>>>>>>>>
>>>>>>>>>>>> The change does not bring any regression, since PETSc is always
>>>>>>>>>>>> giving an error without it, but on the contrary, it improves the
>>>>>>>>>>>> range of applicability of SLEPc, e.g., for complex Hermitian
>>>>>>>>>>>> problems with SBAIJ matrices and bs > 1 that _don’t_ require the
>>>>>>>>>>>> flag MAT_HERMITIAN set to true.
>>>>>>>>>>>>
>>>>>>>>>>>> Thanks,
>>>>>>>>>>>> Pierre
>
