lapackpp-devel — Discussion of Lapack++ development

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

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

With respect to void LaEigSolve, DGEEV/DGEEVX seem to return
eigenvalue/eigenvector pairs with the possibility that the eigenvalues might
be complex numbers.  

Maybe it would be better to implement it using dgesdd?  This may need a
transpose of the V matrix though, which is probably bad.

Alternatively, I could have it return a LaGenMatComplex for the eigenvector
component.

Any thoughts on this?

Jack


> -----Original Message-----
> From: Jacob (Jack) Gryn [mailto:jg...@cs...]
> Sent: March 30, 2005 12:37 PM
> To: 'Christian Stimming'
> Cc: 'lap...@li...'
> Subject: RE: [Lapackpp-devel] Pseudoinverse, Eigenvalues and Eigenvectors
> 
> 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

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

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

I guess my vision is to ultimately have something at a more abstracted
level, like matlab, where you could write, in C, something like =
C=3Dinv(A*B).


I guess this is another projected for another day :) =20

With respect to the inverse functions, I don't see why 2 lapack =
functions
can't be combined to make life simpler for people when they have the =
choice
not to use it.  In my particular application, it's not run very often =
(maybe
once in the entire program), so I would rather make my code easier to =
read
in such a case. =20

Ultimately, I'll leave the decision up to you, if you feel it's
inappropriate, I'll just stick the function in my own local library =
instead,
along with the pseudoinverse.

Jack

> -----Original Message-----
> From: Christian Stimming [mailto:sti...@tu...]
> Sent: March 30, 2005 3:01 PM
> To: Jacob (Jack) Gryn
> Cc: lap...@li...
> Subject: Re: [Lapackpp-devel] Pseudoinverse, Eigenvalues and =
Eigenvectors
>=20
> Hi,
>=20
> Am Mittwoch, 30. M=E4rz 2005 19:37 schrieb Jacob (Jack) Gryn:
> > I'll work on LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV) =
and
> > LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV, =
LaVectorDouble
> > &work)
>=20
> Good.
>=20
> > How about adding
> > LaInverse(const LaGenMatDouble &A, LaGenMatDouble &invA)
> > LaInverseIP(LaGenMatDouble &A)
> > For those who want to do it all in one shot?
>=20
> Sigh... do you desperately need these? I would rather prefer not to =
offer
> these as extra functions, because it's just totally against the =
"spirit"
> of a
> numerically efficient library. In 95% of the cases people are better =
of by
> using the solver functions... and already the existence of the inverse
> functions will distract from that fact.
>=20
> > I'll also work on
> >
> > void LaEigSolve(const LaGenMatDouble &A, LaVectorDouble &eigvals,
> > LaGenMatDouble &eigvec)
>=20
> That's good. I'm looking forward to your contributions.
>=20
> Christian
>=20
> >
> > 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)=3Dinv(A'A)A', or to do the following (sorry =
for the
> > >
> > > matlab
> > >
> > > > style notation) [U,D,V']=3Dsvd(A).  pinv(A)=3DVpinv(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=3D2e-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=3D6595&alloc_id=3D14396&op=3Dclick
> > > _______________________________________________
> > > lapackpp-devel mailing list
> > > lap...@li...
> > > https://lists.sourceforge.net/lists/listinfo/lapackpp-devel
> >
> > -------------------------------------------------------
> > 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=3D6595&alloc_id=3D14396&op=3Dclick
> > _______________________________________________
> > lapackpp-devel mailing list
> > lap...@li...
> > https://lists.sourceforge.net/lists/listinfo/lapackpp-devel

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

[Christian S. <sti...@tu...>]

Hi,

Am Mittwoch, 30. M=E4rz 2005 19:37 schrieb Jacob (Jack) Gryn:
> I'll work on LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV) and
> LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV, LaVectorDouble
> &work)

Good.

> How about adding
> LaInverse(const LaGenMatDouble &A, LaGenMatDouble &invA)
> LaInverseIP(LaGenMatDouble &A)
> For those who want to do it all in one shot?

Sigh... do you desperately need these? I would rather prefer not to offer=20
these as extra functions, because it's just totally against the "spirit" of=
 a=20
numerically efficient library. In 95% of the cases people are better of by=
=20
using the solver functions... and already the existence of the inverse=20
functions will distract from that fact.

> I'll also work on
>
> void LaEigSolve(const LaGenMatDouble &A, LaVectorDouble &eigvals,
> LaGenMatDouble &eigvec)

That's good. I'm looking forward to your contributions.

Christian

>
> 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 Eigenvecto=
rs
> >
> > Hi,
> >
> > Jacob (Jack) Gryn schrieb:
> > > Apparently Lapack has a function named DGETRI for obtaining the inver=
se
> >
> > 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)=3Dinv(A'A)A', or to do the following (sorry for the
> >
> > matlab
> >
> > > style notation) [U,D,V']=3Dsvd(A).  pinv(A)=3DVpinv(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 lo=
op
> >
> > 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=3D2e-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=3D6595&alloc_id=3D14396&op=3Dclick
> > _______________________________________________
> > lapackpp-devel mailing list
> > lap...@li...
> > https://lists.sourceforge.net/lists/listinfo/lapackpp-devel
>
> -------------------------------------------------------
> 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=3D6595&alloc_id=3D14396&op=3Dclick
> _______________________________________________
> lapackpp-devel mailing list
> lap...@li...
> https://lists.sourceforge.net/lists/listinfo/lapackpp-devel

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

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

[Christian S. <sti...@tu...>]

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

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

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

Hi,

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?).

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.

Any thoughts on the above?  

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 ? 

Let me know

Jack

> -----Original Message-----
> From: lap...@li... [mailto:lapackpp-devel-
> ad...@li...] On Behalf Of Jacob (Jack) Gryn
> Sent: March 24, 2005 4:54 PM
> To: lap...@li...
> Subject: [Lapackpp-devel] Pseudoinverse, Eigenvalues and Eigenvectors
> 
> Just wondering if there's a fast way to do the following with lapack++ (if
> not, then in LAPACK) to do the following functions:
> 
> a) Pseudoinverse (maybe I should do inv(A'A)A', but inv() needs to be done
> via solving the linear equation AA'=I (is there a faster way to do this?)
> 
> b) Get eigenvectors and eigenvalues of a general matrix?  (or even of a
> symmetric matrix)?
> 
> Let me know
> 
> Thanks
> 
> Jack
> 
> 
> 
> 
> -------------------------------------------------------
> 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

[Lapackpp-devel] Re: question

[Christian S. <sti...@tu...>]

Dear Jorge,

sorry for repeating but I already said: Please use the mailing list=20
lapackpp-devel as noted on the project page=20
http://sourceforge.net/projects/lapackpp instead of writing to individual=20
developers.

As for documentation: http://lapackpp.sourceforge.net
Installing it under windows xp: Download the setup.exe and execute it.=20

I'm sorry, but you have to figure out these things for yourself. The projec=
t=20
is a voluntary effort and I and all the other developers are busy with othe=
r=20
things. Ask again if you have questions that are *not* answered in the=20
documentation already.

Christian Stimming

Am Sonntag, 27. M=E4rz 2005 05:12 schrieb Jorge Ivan Gaviria Orozco:
> Ok. Thanks.
> My question specifically is where i can find documentation about the
> library.
> i need to know  how install it under windows xp and what functions i can
> use thanks for your help
>
> ----- Original Message -----
> From: "Christian Stimming" <sti...@tu...>
> To: "Jorge Ivan Gaviria Orozco" <eji...@tu...>
> Cc: <lap...@li...>
> Sent: Saturday, March 26, 2005 12:10 PM
> Subject: Re: question
>
> Am Samstag, 26. M=E4rz 2005 17:29 schrieb Jorge Ivan Gaviria Orozco:
> > Hi, im a student of Antioquia University and in my class of numerical
> > methods we are working with the lapack++ library, i wish to know is you
>
> can
>
> > help us with sime documentation and information about the
> > library.______________________________
>
> Dear Sir,
>
> of course I can answer questions about Lapack++. However, I didn't get wh=
at
> you were trying to ask. Please say clearly what you need to know. Also,
> please use the mailing list lapackpp-devel as noted on the project page
> http://sourceforge.net/projects/lapackpp instead of writing to individual
> developers.
>
> Thanks,
>
> Christian Stimming
>
>
>
> ______________________________
> Visita http://www.tutopia.com y comienza a navegar m=E1s r=E1pido en Inte=
rnet.
> Tutopia es Internet para todos.

[Lapackpp-devel] Re: question

[Christian S. <sti...@tu...>]

Am Samstag, 26. M=E4rz 2005 17:29 schrieb Jorge Ivan Gaviria Orozco:
> Hi, im a student of Antioquia University and in my class of numerical
> methods we are working with the lapack++ library, i wish to know is you c=
an
> help us with sime documentation and information about the
> library.______________________________=20

Dear Sir,

of course I can answer questions about Lapack++. However, I didn't get what=
=20
you were trying to ask. Please say clearly what you need to know. Also,=20
please use the mailing list lapackpp-devel as noted on the project page=20
http://sourceforge.net/projects/lapackpp instead of writing to individual=20
developers.

Thanks,

Christian Stimming

Re: [Lapackpp-devel] Cholesky Decomposition A = D'D

[Christian S. <sti...@tu...>]

Am Donnerstag, 24. M=E4rz 2005 22:18 schrieb Jacob (Jack) Gryn:
> > -----Original Message-----
> > Err... No. The *Fact* classes are mentioned in the "Lapack++ 1.0 User
> > Guide", but that's about it.
>
> Strange, I did a google on
> LaSpdFactDouble and
> LaSymmBandFactDouble
>
> It turned up 0 results.

My point is that the only additional (and old) documentation is the one on=
=20
http://math.nist.gov/lapack++/ . Nothing else exists.

> > > As it stands, I don't really see anything to convert from a symmetric
> > > LaGenMatDouble to a LaSymmBandMatDouble or vice versa.
> >
> > A general conversion doesn't make sense, does it? I guess you mean if
> > and only if a general matrix happens to be a symmetric matrix, then you
> > would like to be able to convert it to the LaSymmBandMatDouble?=20
>
> Well, basically, I multiply a few LaGenMatDouble matrices together, and by
> construction the result is a Symmetric, Positive Definite matrix.
>
> Now, this matrix, I need to do a Cholesky factorization on.

Sure, I understand. Of course Lapack[pp] has no means of knowing that this=
=20
result happens to be symmetric, and it probably isn't even exactly symmetri=
c=20
because of rounding errors (on the order of +-1e-15). So in any case you ha=
ve=20
to construct a new LaSymmMatDouble class which you need to tell which part =
of=20
your general matrix should be used. By the way, maybe the constructor which=
=20
takes a (double*), i.e. a pointer to an array of doubles, is the right choi=
ce=20
for this, and A(0,0) will return the pointer to the full array of the=20
original matrix:

LaGenMatDouble A(10,20);
// do something with A; then=20
LaSymmMatDouble B( A(0,0), A.size(0), A.size(1) );

but I'm not sure which part of A is used in B and which one is ignored (upp=
er=20
vs. lower triangular matrix). You would have to look this up in the=20
documentation of the fortran lapack functions that are used.

Again, why are you also talking about the banded matrices?

> Just to test, I created a function in linslv.cc / laslv.h called
> void LaMatFactorize(const LaGenMatDouble &A,LaGenMatDouble& AF)
> based on the function in spdfd.h, but used it with a LaGenMatDouble as
> input.
>
> The result seems to be 'almost' correct.  The function will only look at
> the upper triangle portion of the input matrix (it assumes symmetry), and
> the result is to be an upper-triangular matrix.
>
> The upper triangular portion of the resulting matrix is correct and does
> compare to the results obtained in matlab.  The bottom triangle needs to =
be
> zeroed out.

That results sounds 'exactly' correct to me, because that's the way how=20
fortran lapack stores symmetric matrices internally. One half of them is=20
simply disregarded, so it doesn't even need to be zeroed out. I guess you=20
refer to what you retrieve when you print it to an std::ostream? In that=20
case, the respective operator would have to treat this similar as the lapac=
k=20
functions, i.e. ignoring one half and printing the full matrix solely based=
=20
on the other half.

> So, as a result, I have a working function, the problem is, if I manually
> zero out the bottom triangle, it may not be as optimal as lapack was
> intended to be.
>
> My questions now would be
> a) Is there a quick way of doing this zero-out?

Again, that part is simply ignored. It shouldn't even be zeroed out, simply=
=20
ignored.

> b) Should I submit this as a patch?

No. See above. If you need this special case for yourself, then that's clea=
ry=20
well placed in your own code.

Thanks for asking, anyway. It's nice to see that the library is actually be=
ing=20
used :-))

Christian

[Lapackpp-devel] Pseudoinverse, Eigenvalues and Eigenvectors

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

Just wondering if there's a fast way to do the following with lapack++ (if
not, then in LAPACK) to do the following functions:

a) Pseudoinverse (maybe I should do inv(A'A)A', but inv() needs to be done
via solving the linear equation AA'=I (is there a faster way to do this?)

b) Get eigenvectors and eigenvalues of a general matrix?  (or even of a
symmetric matrix)?

Let me know

Thanks

Jack

RE: [Lapackpp-devel] Cholesky Decomposition A = D'D

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

> -----Original Message-----
> Err... No. The *Fact* classes are mentioned in the "Lapack++ 1.0 User
> Guide", but that's about it.

Strange, I did a google on
LaSpdFactDouble and
LaSymmBandFactDouble

It turned up 0 results.

> 
> > As it stands, I don't really see anything to convert from a symmetric
> > LaGenMatDouble to a LaSymmBandMatDouble or vice versa.
> 
> A general conversion doesn't make sense, does it? I guess you mean if
> and only if a general matrix happens to be a symmetric matrix, then you
> would like to be able to convert it to the LaSymmBandMatDouble? But
> "Band" also refers to a banded matrix, doesn't it? In that case, the
> internal storage structure is probably totally different anyways, so you
> have to make an element-wise copy. In theory some function for this
> could be added to Lapackpp, yes.

Well, basically, I multiply a few LaGenMatDouble matrices together, and by
construction the result is a Symmetric, Positive Definite matrix.

Now, this matrix, I need to do a Cholesky factorization on.

Just to test, I created a function in linslv.cc / laslv.h called
void LaMatFactorize(const LaGenMatDouble &A,LaGenMatDouble& AF)
based on the function in spdfd.h, but used it with a LaGenMatDouble as
input.

The result seems to be 'almost' correct.  The function will only look at the
upper triangle portion of the input matrix (it assumes symmetry), and the
result is to be an upper-triangular matrix.

The upper triangular portion of the resulting matrix is correct and does
compare to the results obtained in matlab.  The bottom triangle needs to be
zeroed out.

So, as a result, I have a working function, the problem is, if I manually
zero out the bottom triangle, it may not be as optimal as lapack was
intended to be.

My questions now would be
a) Is there a quick way of doing this zero-out?
b) Should I submit this as a patch?
c) If (b), is it ok in linslv.cc/laslv.h, or better elsewhere?

Jack

Re: [Lapackpp-devel] Cholesky Decomposition A = D'D

[Christian S. <sti...@tu...>]

Jacob (Jack) Gryn schrieb:
> Is there any documentation on the following matrices:
> LaSymmBandMatDouble
> LaSymmBandFactDouble
> LaSpdMatDouble
> LaSpdFactDouble

Err... No. The *Fact* classes are mentioned in the "Lapack++ 1.0 User 
Guide", but that's about it.

> As it stands, I don't really see anything to convert from a symmetric
> LaGenMatDouble to a LaSymmBandMatDouble or vice versa.  

A general conversion doesn't make sense, does it? I guess you mean if 
and only if a general matrix happens to be a symmetric matrix, then you 
would like to be able to convert it to the LaSymmBandMatDouble? But 
"Band" also refers to a banded matrix, doesn't it? In that case, the 
internal storage structure is probably totally different anyways, so you 
have to make an element-wise copy. In theory some function for this 
could be added to Lapackpp, yes.

> Is it possible to re-implement DPOTRF and/or DPBTRF to work with
> LaGenMatDouble's?

There is dgetrf, and LUFactorizeIP from laslv.h is using it. Is this 
what you were looking for?

Christian

RE: [Lapackpp-devel] Cholesky Decomposition A = D'D

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

Is there any documentation on the following matrices:
LaSymmBandMatDouble
LaSymmBandFactDouble
LaSpdMatDouble
LaSpdFactDouble

As it stands, I don't really see anything to convert from a symmetric
LaGenMatDouble to a LaSymmBandMatDouble or vice versa. =20

Is it possible to re-implement DPOTRF and/or DPBTRF to work with
LaGenMatDouble's?

Jack

> -----Original Message-----
> From: lap...@li... =
[mailto:lapackpp-devel-
> ad...@li...] On Behalf Of Christian Stimming
> Sent: March 16, 2005 3:12 AM
> To: Jacob (Jack) Gryn
> Cc: lap...@li...
> Subject: Re: [Lapackpp-devel] Cholesky Decomposition A =3D D'D
>=20
> Hi,
>=20
> Jacob (Jack) Gryn schrieb:
> > I was wondering if anyone had implemented Cholesky Decomposition to
> > obtain D or D' from A, in the equation A=3DD'D (with D' as the =
transpose
> > of D)?  I noticed a Cholesky factorization in linslv.cc, but I don't
> > think it's what I'm looking for.
>=20
> In linslv.cc the name only refers to the fact that the linear solver
> internally uses such a decomposition, so that doesn't help you.
>=20
> > The LAPACK functions I most likely need are either
> >
> > DPOTRF  or  DPBTRF
>=20
> Oh, wait, these two functions are already all that you need? "grep
> dpotrf include/*.h" tells you that in fact these are available -- =
either
> as inline functions in fmd.h or in the classes in spdfd.h or sybfd.h.
>=20
> However, I never tested any of these, so you should verify they really
> do what they are supposed to. If there is any problem, I will happily
> assist you in fixing them.
>=20
> Christian
>=20
>=20
>=20
> -------------------------------------------------------
> 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_ide95&alloc_id=14396&op=3Dick
> _______________________________________________
> lapackpp-devel mailing list
> lap...@li...
> https://lists.sourceforge.net/lists/listinfo/lapackpp-devel

Re: [Lapackpp-devel] Updated display for lapackpp to display matrices in MATLAB and MAPLE formats

[Christian S. <sti...@tu...>]

Hi Jacob,

Jacob (Jack) Gryn schrieb:
>>- Could you please move the enum pFormat into the LaPreference class
>>instead of having it as a global type? 
> 
> I didn't know I can do typedefs inside class definitions, but I'll play
> around with it.  I may need to change the num type to an int, and use a few
> definitions.  The problem with this is that you'd always need to check if
> the value is valid.    In any case, when using the enum type, aren't the
> possible values restricted only to the scope of those variables?

Errr... IIRC you can surely have typedefs inside classes. That's what 
the STL uses all the time. Anyway, I've been using enums in classes 
before, see e.g.
http://simthetic.sourceforge.net/simthetic/api-doc/html/propertylist_8h-source.html

class Property {
   public:
     enum p_type { bool_t, int_t, /*...*/ };
     //...
}

(The expected result is obtained with or without typedef; I'm not 
totally sure why one or the other was better.)

>>- Right now the code of the LaGenMatDouble class is copied into
>>LaGenMatFloat. Since printing of matrices is not at all time-critical, I
>>was wondering whether you can get rid of this code duplication by
>>calling the LaGenMatDouble::operator<< from the
>>LaGenMatFloar::operator<< with a locally converted copy of the Float
>>matrix? 
> 
> I know, this bothered me as well that I had to make a copy.  Can I easily
> convert the datatypes from LaGenMatFloat to LaGenMatDouble?  

No, you cannot easily copy these all at once. You can only create a new 
matrix and then copy all elements one-by-one, because only the 
conversion from one float to one double is well-defined.

> This is why I
> think using a class inheritance system for all the matrices and vectors
> would be better.

We have discussed this before, haven't we? In short, this would only be 
possible if the element access functions were virtual functions, and 
this adds an extra vtable lookup for each element access, which is a 
considerable overhead which is not tolerable. Therefore the element 
access functions *must* stay inline functions which means they cannot be 
virtual, and this in turn means that any class inheritance is rather 
pointless.

>>- Your explanation is quite nice. Can you add that to laprefs.h in
>>doxygen format? That way, other people would immediately find this, too.
> 
> Ok.  I'll try and do this.
> 
> Jack

Thanks.

Christian

RE: [Lapackpp-devel] Updated display for lapackpp to display matrices in MATLAB and MAPLE formats

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

> -----Original Message-----
> From: lap...@li... =
[mailto:lapackpp-devel-
> ad...@li...] On Behalf Of Christian Stimming
> Sent: Thursday, March 17, 2005 3:24 AM
> To: Jacob (Jack) Gryn
> Cc: lap...@li...
> Subject: Re: [Lapackpp-devel] Updated display for lapackpp to display
> matrices in MATLAB and MAPLE formats
>=20
> Dear Jacob,
>=20
> ok, thanks f=FCr adding this. It looks quite nice. Just some small
> nitpicking from my side:
>=20
> - Could you please move the enum pFormat into the LaPreference class
> instead of having it as a global type? Currently, the symbols NORMAL,
> MATLAB, MAPLE are global values of the enum and might very well =
collide
> with some defines or whatever in some application. If they are inside
> the LaPreference class, they can be used as LaPreference::NORMAL
> everywhere and will avoid any collisions.

I didn't know I can do typedefs inside class definitions, but I'll play
around with it.  I may need to change the num type to an int, and use a =
few
definitions.  The problem with this is that you'd always need to check =
if
the value is valid.    In any case, when using the enum type, aren't the
possible values restricted only to the scope of those variables?

> - Right now the code of the LaGenMatDouble class is copied into
> LaGenMatFloat. Since printing of matrices is not at all time-critical, =
I
> was wondering whether you can get rid of this code duplication by
> calling the LaGenMatDouble::operator<< from the
> LaGenMatFloar::operator<< with a locally converted copy of the Float
> matrix? That way, the actual implementation of the output code would
> only exist once for doubles, and once for complex, but not more than =
that.

I know, this bothered me as well that I had to make a copy.  Can I =
easily
convert the datatypes from LaGenMatFloat to LaGenMatDouble?  This is why =
I
think using a class inheritance system for all the matrices and vectors
would be better.

> - Your explanation is quite nice. Can you add that to laprefs.h in
> doxygen format? That way, other people would immediately find this, =
too.

Ok.  I'll try and do this.

Jack

> Thanks for this contribution. If you are satisfied with the result, I
> would prepare a new release (2.1.3 or 2.2.0).
>=20
> Christian
>=20
>=20
> Jacob (Jack) Gryn schrieb:
>=20
> > Ok.
> >
> >
> >
> > So I=92ve committed the addition of an LaPreferences class, along =
with
> > modified LaGenMatDouble, LaGenMatFloat, and LaGenMatComplex classes =
to
> > support the new output formats (only operator<< has changed, along =
with
> > the additional include of LA_PREFS_H).
> >
> >
> >
> > The following is how one would modify the output display format to =
be
> > compatible with their favourite math program:
> >
> >
> >
> > Place
> >
> > #include LA_PREFS_H
> >
> > In the include statements, somewhere after =93lafnames.h=94
> >
> >
> >
> > At the beginning of your code, call
> >
> >
> >
> > LaPreferences::setPrintFormat(pFormat p, bool newlines);
> >
> >
> >
> > Where pFormat is an enum type, that can be one of the following =
{NORMAL,
> > MATLAB, MAPLE},
> >
> > and bool newlines toggles multiline matrix output (true =3D places a
> > newline after each matrix row, false uses only the appropriate
> > MATLAB/MAPLE delimiter) =96 this option is ignored if pFormat =3D=3D =
NORMAL
> >
> >
> >
> > Enjoy! J
> >
> >
> >
> > BTW, I only did those 3 classes for now, it would take a while for =
me to
> > implement and test all the classes, many of which I=92ve never used
> > before; but if anyone wants, feel free to make appropriate changes =
based
> > on what=92s been done to gmc.cc, gmd.cc, gmf.cc.
> >
> >
> >
> > I really think it would be simpler to do if the matrix classes were
> > derived from a single base class, so only one version of the =
operator<<
> > would be needed (except for the ones with complex #=92s).
> >
> >
> >
> > Jack
> >
>=20
>=20
>=20
> -------------------------------------------------------
> 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_ide95&alloc_id=14396&op=3Dick
> _______________________________________________
> lapackpp-devel mailing list
> lap...@li...
> https://lists.sourceforge.net/lists/listinfo/lapackpp-devel

Re: [Lapackpp-devel] Updated display for lapackpp to display matrices in MATLAB and MAPLE formats

[Christian S. <sti...@tu...>]

Dear Jacob,

ok, thanks f=FCr adding this. It looks quite nice. Just some small=20
nitpicking from my side:

- Could you please move the enum pFormat into the LaPreference class=20
instead of having it as a global type? Currently, the symbols NORMAL,=20
MATLAB, MAPLE are global values of the enum and might very well collide=20
with some defines or whatever in some application. If they are inside=20
the LaPreference class, they can be used as LaPreference::NORMAL=20
everywhere and will avoid any collisions.

- Right now the code of the LaGenMatDouble class is copied into=20
LaGenMatFloat. Since printing of matrices is not at all time-critical, I=20
was wondering whether you can get rid of this code duplication by=20
calling the LaGenMatDouble::operator<< from the=20
LaGenMatFloar::operator<< with a locally converted copy of the Float=20
matrix? That way, the actual implementation of the output code would=20
only exist once for doubles, and once for complex, but not more than that.

- Your explanation is quite nice. Can you add that to laprefs.h in=20
doxygen format? That way, other people would immediately find this, too.

Thanks for this contribution. If you are satisfied with the result, I=20
would prepare a new release (2.1.3 or 2.2.0).

Christian


Jacob (Jack) Gryn schrieb:

> Ok.
>=20
> =20
>=20
> So I=92ve committed the addition of an LaPreferences class, along with=20
> modified LaGenMatDouble, LaGenMatFloat, and LaGenMatComplex classes to=20
> support the new output formats (only operator<< has changed, along with=
=20
> the additional include of LA_PREFS_H).=20
>=20
> =20
>=20
> The following is how one would modify the output display format to be=20
> compatible with their favourite math program:
>=20
> =20
>=20
> Place
>=20
> #include LA_PREFS_H=20
>=20
> In the include statements, somewhere after =93lafnames.h=94
>=20
> =20
>=20
> At the beginning of your code, call
>=20
> =20
>=20
> LaPreferences::setPrintFormat(pFormat p, bool newlines);
>=20
> =20
>=20
> Where pFormat is an enum type, that can be one of the following {NORMAL=
,=20
> MATLAB, MAPLE},
>=20
> and bool newlines toggles multiline matrix output (true =3D places a=20
> newline after each matrix row, false uses only the appropriate=20
> MATLAB/MAPLE delimiter) =96 this option is ignored if pFormat =3D=3D NO=
RMAL
>=20
> =20
>=20
> Enjoy! J
>=20
> =20
>=20
> BTW, I only did those 3 classes for now, it would take a while for me t=
o=20
> implement and test all the classes, many of which I=92ve never used=20
> before; but if anyone wants, feel free to make appropriate changes base=
d=20
> on what=92s been done to gmc.cc, gmd.cc, gmf.cc. =20
>=20
> =20
>=20
> I really think it would be simpler to do if the matrix classes were=20
> derived from a single base class, so only one version of the operator<<=
=20
> would be needed (except for the ones with complex #=92s).
>=20
> =20
>=20
> Jack
>=20

[Lapackpp-devel] Updated display for lapackpp to display matrices in MATLAB and MAPLE formats

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

Ok.

 

So I've committed the addition of an LaPreferences class, along with
modified LaGenMatDouble, LaGenMatFloat, and LaGenMatComplex classes to
support the new output formats (only operator<< has changed, along with the
additional include of LA_PREFS_H).  

 

The following is how one would modify the output display format to be
compatible with their favourite math program:

 

Place

#include LA_PREFS_H  

In the include statements, somewhere after "lafnames.h"

 

At the beginning of your code, call

 

LaPreferences::setPrintFormat(pFormat p, bool newlines);

 

Where pFormat is an enum type, that can be one of the following {NORMAL,
MATLAB, MAPLE}, 

and bool newlines toggles multiline matrix output (true = places a newline
after each matrix row, false uses only the appropriate MATLAB/MAPLE
delimiter) - this option is ignored if pFormat == NORMAL

 

Enjoy! :-)

 

BTW, I only did those 3 classes for now, it would take a while for me to
implement and test all the classes, many of which I've never used before;
but if anyone wants, feel free to make appropriate changes based on what's
been done to gmc.cc, gmd.cc, gmf.cc.  

 

I really think it would be simpler to do if the matrix classes were derived
from a single base class, so only one version of the operator<< would be
needed (except for the ones with complex #'s).

 

Jack

Re: [Lapackpp-devel] Updating operator<< for matlab, maple, etc.

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

I like the reasoning of putting it LaException since it's an "exceptional"
case.  :) 

I may try the LaPref class. 

BTW, there is a class variable for each matrix type called info_ or
something like that; this is also looked at when printing matrices, shall
this possibly be moved into LaPref?

I may just end up doing it the easy way with LaException, I guess I'll let
you know once it's done.

Thanks

Jack

-----Original Message-----
From: Christian Stimming [mailto:sti...@tu...] 
Sent: Wednesday, March 16, 2005 9:52 AM
To: Jacob (Jack) Gryn
Cc: lap...@li...
Subject: Re: [Lapackpp-devel] Updating operator<< for matlab, maple, etc.

Jacob (Jack) Gryn schrieb:

> I guess I can put it in LaException, but that would be kind of wrong,
since
> exceptions generally have nothing to do with the display of matrices.
Maybe
> it would be better to create a globals.h file or put it in lafnames.h ?

I know, I know. On the other hand so far we avoided any extra globals.h 
file or a corresponding concept. There simply isn't anything global for 
the lapackpp library so far (except for LaException::_print). I'd 
therefore stick with this place so far unless we have more than this 
single reason for introducing more "globality" :-)

On the other hand, Lapackpp is meant for high-performance numerics, and 
"printing" a matrix is therefore quite "exceptional" because it is no 
longer concerned with high-performance and/or numerics in general. 
Anyway, for now I'd ask to stick with the LaException class, with many 
explanations that if there is more than this single global flag then 
this might be moved to a separate LaPreference class, which should be 
added in a separate lapref.h header file (source in src/lapref.cc). 
Well, if you are up to adding this, you can just as well go ahead with 
that right now. I'm too lazy, but feel free to do this.

Christian



> 
> Jack
> 
> -----Original Message-----
> From: Christian Stimming [mailto:sti...@tu...] 
> Sent: Wednesday, March 16, 2005 3:28 AM
> To: Jacob (Jack) Gryn
> Cc: lap...@li...
> Subject: Re: [Lapackpp-devel] Updating operator<< for matlab, maple, etc.
> 
> Hi,
> 
> Jacob (Jack) Gryn schrieb:
> 
>>I've been pretty busy with projects lately, so I haven't had a chance to 
>>properly implement this into lapackpp; 
> 
> 
> No problem. Take your time :-)) as we all do.
> 
> 
>>but now that I have an 
>>opportunity, as I mentioned a few months back, I would like to change 
>>the operator<< output to give a display that can easily be copied/pasted 
>>into matlab or maple (or any other math software people may wish to use).
> 
> 
> Ok.
> 
> 
>>My question is this:
>>
>>I would like to make some sort of a global variable called something 
>>like DISP_FORMAT, that can be set to one of the following: {NORMAL, 
>>MATLAB, MAPLE, .}
> 
>  > Where should I put this variable that it can be accessible to the
matrix
> 
>>classes and most programs?
> 
> 
> So far we don't have any global variables *except* for the 
> LaException::_print variable.
> 
> Each of the matrix classes also has a static flag whether debug mode is 
> turned on or not, but a matrix class would probably be the wrong place 
> for the disp_format flag because that one should affect all matrix 
> classes at once.
> 
> 
>>Maybe it should be a member or class variable of the matrix classes 
>>instead, calling something like myGMD.setDisplay(MATLAB), or 
>>LaGenMatDouble::setDisplay(MATLAB). ?
> 
> 
> Probably not quite the best idea, because if you're using several 
> different matrix classes you would have to call
> LaGenMatDouble::setDisplay(bla);
> LaGenMatComplex::setDisplay(bla);
> LaSymmMatDouble::setDisplay(bla);
> 
> and so on. I would rather want to have only one single library-wide flag 
> for this and this means there is also only one single library-wide 
> function to set this. Hm... that one should either go into 
> LaGenMatDouble because that's the central matrix class, or we abuse the 
> LaException class for this because that's the only class in lapackpp 
> that is (potentially) in use by everything. I might actually prefer 
> LaException::setMatrixDisplay()... do you think that would work?
> 
> Christian
> 
> 
> 
> 

Re: [Lapackpp-devel] Updating operator<< for matlab, maple, etc.

[Christian S. <sti...@tu...>]

Jacob (Jack) Gryn schrieb:

> I guess I can put it in LaException, but that would be kind of wrong, since
> exceptions generally have nothing to do with the display of matrices.  Maybe
> it would be better to create a globals.h file or put it in lafnames.h ?

I know, I know. On the other hand so far we avoided any extra globals.h 
file or a corresponding concept. There simply isn't anything global for 
the lapackpp library so far (except for LaException::_print). I'd 
therefore stick with this place so far unless we have more than this 
single reason for introducing more "globality" :-)

On the other hand, Lapackpp is meant for high-performance numerics, and 
"printing" a matrix is therefore quite "exceptional" because it is no 
longer concerned with high-performance and/or numerics in general. 
Anyway, for now I'd ask to stick with the LaException class, with many 
explanations that if there is more than this single global flag then 
this might be moved to a separate LaPreference class, which should be 
added in a separate lapref.h header file (source in src/lapref.cc). 
Well, if you are up to adding this, you can just as well go ahead with 
that right now. I'm too lazy, but feel free to do this.

Christian



> 
> Jack
> 
> -----Original Message-----
> From: Christian Stimming [mailto:sti...@tu...] 
> Sent: Wednesday, March 16, 2005 3:28 AM
> To: Jacob (Jack) Gryn
> Cc: lap...@li...
> Subject: Re: [Lapackpp-devel] Updating operator<< for matlab, maple, etc.
> 
> Hi,
> 
> Jacob (Jack) Gryn schrieb:
> 
>>I've been pretty busy with projects lately, so I haven't had a chance to 
>>properly implement this into lapackpp; 
> 
> 
> No problem. Take your time :-)) as we all do.
> 
> 
>>but now that I have an 
>>opportunity, as I mentioned a few months back, I would like to change 
>>the operator<< output to give a display that can easily be copied/pasted 
>>into matlab or maple (or any other math software people may wish to use).
> 
> 
> Ok.
> 
> 
>>My question is this:
>>
>>I would like to make some sort of a global variable called something 
>>like DISP_FORMAT, that can be set to one of the following: {NORMAL, 
>>MATLAB, MAPLE, .}
> 
>  > Where should I put this variable that it can be accessible to the matrix
> 
>>classes and most programs?
> 
> 
> So far we don't have any global variables *except* for the 
> LaException::_print variable.
> 
> Each of the matrix classes also has a static flag whether debug mode is 
> turned on or not, but a matrix class would probably be the wrong place 
> for the disp_format flag because that one should affect all matrix 
> classes at once.
> 
> 
>>Maybe it should be a member or class variable of the matrix classes 
>>instead, calling something like myGMD.setDisplay(MATLAB), or 
>>LaGenMatDouble::setDisplay(MATLAB). ?
> 
> 
> Probably not quite the best idea, because if you're using several 
> different matrix classes you would have to call
> LaGenMatDouble::setDisplay(bla);
> LaGenMatComplex::setDisplay(bla);
> LaSymmMatDouble::setDisplay(bla);
> 
> and so on. I would rather want to have only one single library-wide flag 
> for this and this means there is also only one single library-wide 
> function to set this. Hm... that one should either go into 
> LaGenMatDouble because that's the central matrix class, or we abuse the 
> LaException class for this because that's the only class in lapackpp 
> that is (potentially) in use by everything. I might actually prefer 
> LaException::setMatrixDisplay()... do you think that would work?
> 
> Christian
> 
> 
> 
> 

RE: [Lapackpp-devel] Updating operator<< for matlab, maple, etc.

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

I guess I can put it in LaException, but that would be kind of wrong, since
exceptions generally have nothing to do with the display of matrices.  Maybe
it would be better to create a globals.h file or put it in lafnames.h ?

Jack

-----Original Message-----
From: Christian Stimming [mailto:sti...@tu...] 
Sent: Wednesday, March 16, 2005 3:28 AM
To: Jacob (Jack) Gryn
Cc: lap...@li...
Subject: Re: [Lapackpp-devel] Updating operator<< for matlab, maple, etc.

Hi,

Jacob (Jack) Gryn schrieb:
> I've been pretty busy with projects lately, so I haven't had a chance to 
> properly implement this into lapackpp; 

No problem. Take your time :-)) as we all do.

> but now that I have an 
> opportunity, as I mentioned a few months back, I would like to change 
> the operator<< output to give a display that can easily be copied/pasted 
> into matlab or maple (or any other math software people may wish to use).

Ok.

> My question is this:
> 
> I would like to make some sort of a global variable called something 
> like DISP_FORMAT, that can be set to one of the following: {NORMAL, 
> MATLAB, MAPLE, .}
 > Where should I put this variable that it can be accessible to the matrix
> classes and most programs?

So far we don't have any global variables *except* for the 
LaException::_print variable.

Each of the matrix classes also has a static flag whether debug mode is 
turned on or not, but a matrix class would probably be the wrong place 
for the disp_format flag because that one should affect all matrix 
classes at once.

> Maybe it should be a member or class variable of the matrix classes 
> instead, calling something like myGMD.setDisplay(MATLAB), or 
> LaGenMatDouble::setDisplay(MATLAB). ?

Probably not quite the best idea, because if you're using several 
different matrix classes you would have to call
LaGenMatDouble::setDisplay(bla);
LaGenMatComplex::setDisplay(bla);
LaSymmMatDouble::setDisplay(bla);

and so on. I would rather want to have only one single library-wide flag 
for this and this means there is also only one single library-wide 
function to set this. Hm... that one should either go into 
LaGenMatDouble because that's the central matrix class, or we abuse the 
LaException class for this because that's the only class in lapackpp 
that is (potentially) in use by everything. I might actually prefer 
LaException::setMatrixDisplay()... do you think that would work?

Christian

Re: [Lapackpp-devel] Updating operator<< for matlab, maple, etc.

[Christian S. <sti...@tu...>]

Hi,

Jacob (Jack) Gryn schrieb:
> I=92ve been pretty busy with projects lately, so I haven=92t had a chan=
ce to=20
> properly implement this into lapackpp;=20

No problem. Take your time :-)) as we all do.

> but now that I have an=20
> opportunity, as I mentioned a few months back, I would like to change=20
> the operator<< output to give a display that can easily be copied/paste=
d=20
> into matlab or maple (or any other math software people may wish to use=
).

Ok.

> My question is this:
>=20
> I would like to make some sort of a global variable called something=20
> like DISP_FORMAT, that can be set to one of the following: {NORMAL,=20
> MATLAB, MAPLE, =85}
 > Where should I put this variable that it can be accessible to the matr=
ix
> classes and most programs?

So far we don't have any global variables *except* for the=20
LaException::_print variable.

Each of the matrix classes also has a static flag whether debug mode is=20
turned on or not, but a matrix class would probably be the wrong place=20
for the disp_format flag because that one should affect all matrix=20
classes at once.

> Maybe it should be a member or class variable of the matrix classes=20
> instead, calling something like myGMD.setDisplay(MATLAB), or=20
> LaGenMatDouble::setDisplay(MATLAB). ?

Probably not quite the best idea, because if you're using several=20
different matrix classes you would have to call
LaGenMatDouble::setDisplay(bla);
LaGenMatComplex::setDisplay(bla);
LaSymmMatDouble::setDisplay(bla);

and so on. I would rather want to have only one single library-wide flag=20
for this and this means there is also only one single library-wide=20
function to set this. Hm... that one should either go into=20
LaGenMatDouble because that's the central matrix class, or we abuse the=20
LaException class for this because that's the only class in lapackpp=20
that is (potentially) in use by everything. I might actually prefer=20
LaException::setMatrixDisplay()... do you think that would work?

Christian

Re: [Lapackpp-devel] Cholesky Decomposition A = D'D

[Christian S. <sti...@tu...>]

Hi,

Jacob (Jack) Gryn schrieb:
> I was wondering if anyone had implemented Cholesky Decomposition to=20
> obtain D or D=92 from A, in the equation A=3DD=92D (with D=92 as the tr=
anspose=20
> of D)?  I noticed a Cholesky factorization in linslv.cc, but I don=92t=20
> think it=92s what I=92m looking for.

In linslv.cc the name only refers to the fact that the linear solver=20
internally uses such a decomposition, so that doesn't help you.

> The LAPACK functions I most likely need are either
>=20
> DPOTRF  or  DPBTRF

Oh, wait, these two functions are already all that you need? "grep=20
dpotrf include/*.h" tells you that in fact these are available -- either=20
as inline functions in fmd.h or in the classes in spdfd.h or sybfd.h.

However, I never tested any of these, so you should verify they really=20
do what they are supposed to. If there is any problem, I will happily=20
assist you in fixing them.

Christian

[Lapackpp-devel] Cholesky Decomposition A = D'D

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

I was wondering if anyone had implemented Cholesky Decomposition to obtain D
or D' from A, in the equation A=D'D (with D' as the transpose of D)?  I
noticed a Cholesky factorization in linslv.cc, but I don't think it's what
I'm looking for.

 

The LAPACK functions I most likely need are either 

DPOTRF  or  DPBTRF

 

Let me know if the existing functions can do what I want, if not, let me
know if you have any tips as to how I can implement one of these functions
into lapackpp.

 

Thanks

 

Jack

[Lapackpp-devel] Updating operator<< for matlab, maple, etc.

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

Hi,

 

I've been pretty busy with projects lately, so I haven't had a chance to
properly implement this into lapackpp; but now that I have an opportunity,
as I mentioned a few months back, I would like to change the operator<<
output to give a display that can easily be copied/pasted into matlab or
maple (or any other math software people may wish to use).

 

My question is this:

 

I would like to make some sort of a global variable called something like
DISP_FORMAT, that can be set to one of the following: {NORMAL, MATLAB,
MAPLE, .}

 

Where should I put this variable that it can be accessible to the matrix
classes and most programs?

 

Maybe it should be a member or class variable of the matrix classes instead,
calling something like myGMD.setDisplay(MATLAB), or
LaGenMatDouble::setDisplay(MATLAB). ?

 

Any thoughts before I begin coding anything?

 

Thanks

 

Jack

[Lapackpp-devel] Re: lapackpp question

[Christian S. <sti...@tu...>]

Dear Neil,

NEIL ANTHONY KLINGENSMITH schrieb:
> I understand that you are an active participant in the development of
> the lapackpp package, and I'd just like to know if it provides
> functionality for inverting square matricies.  I was unable to find any
> information about that in the documentation, but it seems such an
> essential operation that I can't believe it would be left out.

The functionality of inverting a matrix is included. However, with some 
background in numerics you will surely have heard that inverting a 
matrix should not be calculated explicitly; rather, at the point where 
an inverted matrix is used, a linear equation system has to be solved 
instead. Also, if you solve a linear equation system with an identity 
matrix as the right hand side your result will be the inverted matrix.

The interesting question rather is: are you looking for general 
real-valued matrices or anything more special (symmetric, sparse, 
complex-valued)?

Anyway, inverting a matrix means solving a system of linear equations, 
and this is supported for all matrices in Lapack++.

Christian Stimming