> On Sep 14, 2019, at 8:56 AM, Pierre Jolivet <pierre.joli...@enseeiht.fr> > wrote: > > > > 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?
Hmm, Is that a reasonable case to worry about? As was noted before doesn't this mean that the matrix is real off the diagonal? Yes, you are right if that is a useful case to indicate in this way then we need to keep all the flags and not simplify as I suggested. > >> 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"? Looked at the code, yes we never wrote and tested versions of the matrix methods for Hermitian matrices for bs > 1. They would need to be implemented. The simple way would be cut and paste and put the appropriate conjugates in the correct places for small bs and then change the BLAS call for the large bs. > > 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