RE: [Lapackpp-devel] Pseudoinverse, Eigenvalues and Eigenvectors

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

I'll work on LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV) and
LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV, LaVectorDouble &work)

How about adding
LaInverse(const LaGenMatDouble &A, LaGenMatDouble &invA)
LaInverseIP(LaGenMatDouble &A)

For those who want to do it all in one shot?

I'll also work on

void LaEigSolve(const LaGenMatDouble &A, LaVectorDouble &eigvals,
LaGenMatDouble &eigvec)

Jack

> -----Original Message-----
> From: lap...@li... [mailto:lapackpp-devel-
> ad...@li...] On Behalf Of Christian Stimming
> Sent: March 30, 2005 3:55 AM
> To: Jacob (Jack) Gryn
> Cc: lap...@li...
> Subject: Re: [Lapackpp-devel] Pseudoinverse, Eigenvalues and Eigenvectors
> 
> Hi,
> 
> Jacob (Jack) Gryn schrieb:
> > Apparently Lapack has a function named DGETRI for obtaining the inverse
> of a
> > function.  Would it be reasonable to create a function called inv() in
> > LaGenMatDouble that makes use of DGETRI?  (If not somewhere else?).
> 
> Ok, I didn't know of DGETRI. Yes, that's reasonable to be implemented.
> However, the input to DGETRI is not any general matrix but the input is
> the result of the DGETRF operation, which in lapack++ is available e.g.
> in LUFactorizeIP in linslv.cc and laslv.h (huh, obviously I forgot to
> add the leading "La" to the name; sorry).
> 
> Therefore I would suggest to add a function like
> LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV) to
> laslv.h/linslv.cc, where the documentation clearly says that the input
> variables *must* come from LUFactorizeIP. I think that's better than
> adding an inv() method to LaGenMatDouble, because that one would always
> have to compute both dgetrf and dgetri but the fact that these are two
> separate steps should be visible to lapackpp's users, IMHO.
> 
> You can decide for yourself what to do with the WORK array of DGETRI --
> if this function is called several times, then you can gain a
> significant speedup if the application re-uses the same WORK array
> instead of allocating a new one each time (really!). Probably you can
> add two functions for the Inverse, one with the WORK array as argument
> and the other one without, where the latter internally allocates a new
> WORK array and then calls the former function.
> 
> > In addition, I would like to build a pseudoinverse function as well (not
> > implemented in lapack).  The two options would be to do it according to
> the
> > definition pinv(A)=inv(A'A)A', or to do the following (sorry for the
> matlab
> > style notation) [U,D,V']=svd(A).  pinv(A)=Vpinv(D)U', where to get
> pinv(D),
> > just find the inverse of each non-zero value in the diagonal matrix.
> >
> > The second way apparently is more computationally efficient, but a loop
> with
> > if statements to get the pseudoinverse of D may take longer than using
> > lapack to calculate the pseudoinverse by definition.
> 
> I'm still quite hesistant to add functions concerning "the inverse" to a
> numerical library (the exception being cases like above, when there
> already exists a LAPACK function for this). Really, all the numerics
> guy's faces turn red and they start sweating once you tell them you want
> to calculate the inverse :-))
> http://www.netlib.org/lapack/lug/node26.html#1232
> 
> But seriously, I found for myself that the definition of "the useful
> pseodo-inverse" depends much more on your actual application than you
> might think, which means the calculation of the pseudoinverse should be
> done in application code instead of lapackpp code.
> 
> The SVD way of calculating it is a good example, because when you want
> to calculate pinv(D) by inverting "each non-zero value", you have to
> think twice what "non-zero" actually means. Does it mean "exactly equal
> to the floating point number 0.0"? Probably not, but rather "smaller in
> absolute value than some value epsilon". If epsilon is the machine's
> precision, then it's approx. eps=2e-16 for doubles. But for your
> application it might turn out that a much better value is something
> smaller err larger, e.g. 1e-8 or 1e-4. But the whole point is that this
> really really really depends on your application. Also on
> http://www.netlib.org/lapack/lug/node53.html it says that there is no
> explicit LAPACK function for the pseudoinverse but "The effective rank,
> k, of A can be determined as the number of singular values which exceed
> a suitable threshold." So I would ask you not to implement such a
> function in lapackpp.
> 
> > Secondly, if I want to obtain the eigenvalues and eigenvectors of a
> general
> > matrix, is DGEEV the correct function to implement?  Should I create a
> > function to add to EigSlv.cc ?
> 
> Yes, and yes. (Or maybe dgeevx?) For dgeev/dgeevx (whichever seems
> easier for you) you can add a function right away. Just go ahead.
> 
> Thanks,
> 
> Christian
> 
> 
> 
> 
> -------------------------------------------------------
> SF email is sponsored by - The IT Product Guide
> Read honest & candid reviews on hundreds of IT Products from real users.
> Discover which products truly live up to the hype. Start reading now.
> http://ads.osdn.com/?ad_id=6595&alloc_id=14396&op=click
> _______________________________________________
> lapackpp-devel mailing list
> lap...@li...
> https://lists.sourceforge.net/lists/listinfo/lapackpp-devel