RE: [Lapackpp-devel] SVD that ignores U, Sigma or VT

[Jacob \(Jack\) G. <jg...@cs...>]

It might get a little confusing to implement, or to understand, but as long
as things don't change for those who want the FULL SVD, it should be ok.

The only other probs I can think of is that the LWORK variable may be
different if JOBZ is different; also, if a number is too big (hence the
reason to ignore part of the SVD, LWORK may also get much larger.

Would you mind implementing it?

Thanks

Jack

-----Original Message-----
From: Christian Stimming [mailto:sti...@tu...] 
Sent: September 3, 2004 4:23 AM
To: Jacob (Jack) Gryn
Cc: lap...@li...
Subject: Re: [Lapackpp-devel] SVD that ignores U, Sigma or VT

Jacob (Jack) Gryn schrieb:
> I'm writing some code that only uses data from VT of an SVD, and don't 
> need the data in U and Sigma.  Although I can ignore it; the resulting U 
> is significantly large (512^2x512^2), which is too large a matrix for 
> the system to handle; and is clearly a waste.

Absolutely agreed. I guess you mean the JOBZ argument of the LAPACK 
function dgesdd http://docs.sun.com/source/817-6700/dgesdd.html, where 
currently only the settings 'A' and 'N' are used, and instead you would 
like to use either the 'O' or the 'S' setting.

> I've noticed that there is a function that ignores both U and VT, by 
> essentially setting ldv and ldu to 0, that seemingly tells LAPACK to 
> ignore the data.  

The point here is that the argument JOBZ was set to 'N'. Only this is 
the reason for LAPACK not to use the U and VT arguments. Look into the 
man page mentioned above.

> Maybe it would be better just to use a single LaSVD_IP 
> function, and allow uninitialized LaGenMat[Double/Complex] as input. 
>  If, for example U is declared with the default constructor, since 
> U.gdim(0) == 0, I would guess that LAPACK will ignore writing output to U.

For people who only want the singular values but not the vectors, I 
would rather want to leave that extra function in there, since it is 
totally unambiguous what this function will do. But for your case where 
you want only some of the singular vectors, there are several 
possibilities. As you say, one possibility would be to do different 
things according to what the input arguments U and VT have been set to. 
Let's assume A to be M-by-N and M != N. Right now we do the following:

If (U.size == M-by-M) && (VT.size == N-by-N)
   calculate all singular vectors (JOBZ='A')
else
   error.

Instead, we could also do the following:

MNmin = min(M,N);
If (U.size == M-by-M) && (VT.size == N-by-N)
   calculate all singular vectors (JOBZ='A')
else if (U.size == M-by-MNmin) && (VT.size == MNmin-by-N)
   calculate the MNmin first singular vectors (JOBZ='S')
else if (M >= N) && (U.size = 0) && (VT.size == N-by-N)
   calculate full VT and first N vectors of U (JOBZ='O')
else if (M < N) && (U.size = M-by-M) && (VT.size = 0)
   calculate full U and first M vectors of VT (JOBZ='O')
else
   error

Does this meet your requirements? On the other hand, this increases the 
risk of the programmer to be totally confused about what the arguments 
are supposed to be. But nevertheless I think it would be possible. Do 
you want me to implement this? Or do you want to give it a try for yourself?

Christian