by jackpoulson » Tue Feb 15, 2011 12:07 pm
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