PDSYNTRD not available for UPLO = 'U'

[naromero]

Hi,

This is a feature request. It would be nice to see:
PDSYNTRD

work with the fast algorithm for the UPLO = 'U' case.

Additionally, it would be good if all the eigensolver made use of the
faster House holder algorithm, i.e.
PDSYEV, PDSYEVD, PDSYEVX,

Nichols A. Romero
Argonne Leadership Computing Facility

[Julien Langou]

Hi Nichols,

Thanks for pointing this out. I believe you are referring to the master thesis of Jack Poulson for PDSYNTRD, and the paper "Accumulating Householder Transformations, Revisited" from Joffrain, Low, Quintana-Ortí, van de Geijn, and Van Zee for the "faster Householder algorithm".

Well, yes, we might definitely want to consider adding these. I am not too familiar with the work personally.

Thanks for pointing these, we will let you know the status of these items if it gets changed.
See: http://www.netlib.org/lapack/WishList/

Cheers,
Julien.

[jackpoulson]

Julien,

My master's thesis included work on an algorithm similar to PDSYNGST, not PDSYNTRD, but, at the time, I was unfortunately not aware of the existence of the work of Sears/Stanley/Henry and the resulting PDSYNGST/PZHENGST algorithm. Since ScaLAPACK already contains the lower-triangular routines, it is simply a matter of transcription to write the upper triangular routines. Because the upper triangle is the conjugate-transpose of the lower-triangle, simply conjugate-transposing each lower-triangular update yields the upper-triangular update, for both PDSYNGST/PZHENGST and PDSYNTRD/PZHENTRD.

For instance, suppose the lower-triangular update was:

CALL PZTRSM( 'Right', 'Lower', 'Conjugate transpose', 'Non-Unit',
$ K, JB, ONE, A, I, J, DESCA, A, I+JB,
$ J, DESCA )

This is simply A21 := A21 tril(A11)^{-H}, and so the conjugate-transpose update is A21^H := tril(A11)^{-1} A21^H, which is clearly equivalent to A12 := triu(A11)^{-H} A12. That update is translated back into PBLAS as

CALL PZTRSM( 'Left', 'Upper', 'No transpose', 'Non-Unit',
$ JB, K, ONE, A, I, J, DESCA, A, I,
$ J+JB, DESCA )

This is obviously a simple case (and forgive me if I translated back into PBLAS incorrectly), but the idea is pretty simple: conjugate-transpose the updates and you can easily translate a lower-triangular storage symmetric routine into an upper-triangular storage version.

This approach would require barely any work for PDSYNGST/PZHENGST, but PDSYNTRD/PZHENTRD is quite detailed...

As for the UT transform blocked Householder application algorithm, it relies on having computed 'true' (generalized) Householder reflectors. In particular, since identity is not a Householder transform, the identity shortcuts would have to be removed in favor of a legal reflection of the form:
H = I - 2 [1; zeros(n-1,1)] [1 zeros(1,n-1)], so that H [alpha11; 0] = [-alpha11; 0].

Jack

[Sarah2648]

Hi,

This is a feature request. It would be nice to see:
PDSYNTRD

work with the fast algorithm for the UPLO = 'U' case.

Additionally, it would be good if all the eigensolver made use of the
faster House holder algorithm, i.e.
PDSYEV, PDSYEVD, PDSYEVX,

Nichols A. Romero
Argonne Leadership Computing Facility

Hey guys,

Also looking for this feature if it has been added? I am also preparing for my masters thesis and i've checked the link that was posted here, but can't seem to find anything on the subject?

"Thanks for pointing these, we will let you know the status of these items if it gets changed.
See: http://www.netlib.org/lapack/WishList/"

Any help would be greatly appreciated

Sarah

[ldm001]

I would like to add my vote to having this included.