Re: A plea for gsvd

[Pascal Dupuis]

On Sat, Apr 22, 2006 at 06:33:07PM +0200, David Bateman wrote:

> There were a couple of very simple bugs in your code that I fixed (size
> of arg to dggsvd and number of F77 args incorrect), fixing that the code
> now appears to work correctly and not seg-fault. Can you check. Also the
> help of the gsvd function needs to be updated as it references svd and
> in particular suggests using "svd(hilb(3))" as an example. Could you
> update this.
> 

Hello,

revised code included. This one works.

Some rought points:
- the doc
- according to dggsvd.f doc, the purpose is to get
   U���*A*Q = D1*( 0 R ),    V���*B*Q = D2*( 0 R )
  As the function also computes R, I added a 6th output argument to
  return R as well. This only works if  M-K-L >= 0 (see the doc),
  otherwise part of R is stored in B
- there are two formulations, this one and another as explained at the
  end of the doc, U���*A*X = diag(alpha_i),   V���*B*X = diag(beta_i), but
  the 'X' in the latter case requires to compute Q*[ I   0; 0 R]. X is a
  factor which simultaneously diagonalise A and B, and is not
  orthogonal. Which one to choose ?
- in the case of 1 output arg, returns the extended singular values
  alpha_i/beta_i, infinite values not being considered as
  problematic. If B is invertible, they are the singular values of A*inv(B)

Best regards

Pascal

-- 
Dr. ir. Pascal Dupuis, conseiller technologique/technological advisor
U. C. Louvain AC/ADRE :  http://www.adre.ucl.ac.be/
Place de L'Universit��, 1; B-1348 Louvain-la-Neuve Belgium
T??l. +32-10-47 92 33; Fax +32-10-47 48 30

gsvd.tar.bz2
Description: Binary data