The VXL Project

Sure: since the two square roots of a complex number c are the two solutions of the equation x^2 - c = 0, you can use class vnl_cpoly_roots.

Example:
  To find the two complex square roots of 3+4i (which are 2+i and -2-i), compile&run the following:

#include <vnl/algo/vnl_cpoly_roots.h>
#include <vcl_iostream.h>
int main()
{
  vcl_complex<double> c(3.0,4.0);
  vnl_vector<vcl_complex<double> > equation(2);
  // Note: the coefficient in x^2 of the equation must be 1
  // (not to be placed in the equation vector);
  // equation[0] is the second highest coefficient, etc.
  equation[0] = 0; equation[1] = -c;
  vnl_cpoly_roots r(equation);
  vcl_cout << "First square root of 3+4i is " << r.solns[0]
           << "\nSecond square root is " << r.solns[1] << vcl_endl;
  return 0;
}

--	Peter.














--


--- Den ons 2010-12-22 skrev Martin Reuter <mreuter@...>:

> Från: Martin Reuter <mreuter@...>
> Ämne: [Vxl-users] complex squareroot?
> Till: vxl-users@...
> Datum: onsdag 22 december 2010 01:40
> Hi,
> 
> is there a way to compute a squareroot of a complex number
> in vnl?
> 
> Actually what I am trying to achieve is to use
> vnl_complex_eigensystem
> on a real matrix to compute eigenvectors and complex
> eigenvalues, then
> compute the squareroot of the eigenvalues and finally the
> squareroot of
> the original matrix. Of course it would be nice to have a
> sqrtm
> functions (see matlab) to do all this for me.
> 
> 
> Thanks for any hints, Martin
> 
> 
> ------------------------------------------------------------------------------
> Forrester recently released a report on the Return on
> Investment (ROI) of
> Google Apps. They found a 300% ROI, 38%-56% cost savings,
> and break-even
> within 7 months.  Over 3 million businesses have gone
> Google with Google Apps:
> an online email calendar, and document program that's
> accessible from your 
> browser. Read the Forrester report: http://p.sf.net/sfu/googleapps-sfnew
> _______________________________________________
> Vxl-users mailing list
> Vxl-users@...
> https://lists.sourceforge.net/lists/listinfo/vxl-users
>

Thanks,

is this (cpoly_roots) computationally efficient, or a lot of overhead
(have to solve this very often).

Another idea I had is to convert to std::complex and use the provided
sqrt function? Does that make sense?

This brings me to another question about vcl_complex. Is that a wrapper
around std::complex?

The vnl:numerics documentation says this:
"Compliance with the ANSI standard C++ library 
The current draft of the ANSI standard (as at May 1996) includes classes
for 1-dimensional vectors (valarray<T>) and complex numbers
(complex<T>). There is no standard for matrices. The current vnl classes
are not implemented in terms of valarray, as there is a potential
perfomance hit, but in the future they might be."

So what's with complex ?? (also performance is spelled with a 2nd 'r').

best, Martin

On Wed, 2010-12-22 at 17:22 +0000, Peter Vanroose wrote:
> Sure: since the two square roots of a complex number c are the two solutions of the equation x^2 - c = 0, you can use class vnl_cpoly_roots.
> 
> Example:
>   To find the two complex square roots of 3+4i (which are 2+i and -2-i), compile&run the following:
> 
> #include <vnl/algo/vnl_cpoly_roots.h>
> #include <vcl_iostream.h>
> int main()
> {
>   vcl_complex<double> c(3.0,4.0);
>   vnl_vector<vcl_complex<double> > equation(2);
>   // Note: the coefficient in x^2 of the equation must be 1
>   // (not to be placed in the equation vector);
>   // equation[0] is the second highest coefficient, etc.
>   equation[0] = 0; equation[1] = -c;
>   vnl_cpoly_roots r(equation);
>   vcl_cout << "First square root of 3+4i is " << r.solns[0]
>            << "\nSecond square root is " << r.solns[1] << vcl_endl;
>   return 0;
> }
> 
> --	Peter.
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> --
> 
> 
> --- Den ons 2010-12-22 skrev Martin Reuter <mreuter@...>:
> 
> > Från: Martin Reuter <mreuter@...>
> > Ämne: [Vxl-users] complex squareroot?
> > Till: vxl-users@...
> > Datum: onsdag 22 december 2010 01:40
> > Hi,
> > 
> > is there a way to compute a squareroot of a complex number
> > in vnl?
> > 
> > Actually what I am trying to achieve is to use
> > vnl_complex_eigensystem
> > on a real matrix to compute eigenvectors and complex
> > eigenvalues, then
> > compute the squareroot of the eigenvalues and finally the
> > squareroot of
> > the original matrix. Of course it would be nice to have a
> > sqrtm
> > functions (see matlab) to do all this for me.
> > 
> > 
> > Thanks for any hints, Martin
> > 
> > 
> > ------------------------------------------------------------------------------
> > Forrester recently released a report on the Return on
> > Investment (ROI) of
> > Google Apps. They found a 300% ROI, 38%-56% cost savings,
> > and break-even
> > within 7 months.  Over 3 million businesses have gone
> > Google with Google Apps:
> > an online email calendar, and document program that's
> > accessible from your 
> > browser. Read the Forrester report: http://p.sf.net/sfu/googleapps-sfnew
> > _______________________________________________
> > Vxl-users mailing list
> > Vxl-users@...
> > https://lists.sourceforge.net/lists/listinfo/vxl-users
> > 
> 
>

> is this (cpoly_roots) computationally efficient, or a lot
> of overhead (have to solve this very often).

I guess it's indeed a lot of overhead, certainly for something as simple as square root. Didn't measure this, though.


> Another idea I had is to convert to std::complex and use
> the provided sqrt function? Does that make sense?

Of course. That's probably a lot better. But I'm not familiar enough with std::sqrt(std::complex<T>) to be sure about this; e.g.: which of the two square roots is returned by this function?


> This brings me to another question about vcl_complex. Is
> that a wrapper around std::complex?

Yes. Or better said: it is (in most cases) just a #define for std::complex.
So the two are (or should be) identical.


> The vnl:numerics documentation says this:
> "Compliance with the ANSI standard C++ library 
> The current draft of the ANSI standard (as at May 1996)
> includes classes for 1-dimensional vectors (valarray<T>)
> and complex numbers (complex<T>).
> There is no standard for matrices.
> The current vnl classes are not implemented in terms of
> valarray, as there is a potential performance hit,
> but in the future they might be."
>
> So what's with complex ??

There is no vnl_complex, just vcl_complex (which is std::complex).


--      Peter.

Yes, std::sqrt seems to work, only need one sqrt anyway.

Now I just need a fast way of computing the inverse of a complex 4x4
matrix. using vnl_inverse gives an error:

 In function ‘vnl_matrix_fixed<T, 4u, 4u> vnl_inverse(const
vnl_matrix_fixed<T, 4u, 4u>&) [with T = std::complex<double>]’:

vnl_inverse.h:114: error: no match for ‘operator==’ in ‘det == 0’

looks like it cannot deal with complex fixed matrices?

Martin

On Wed, 2010-12-22 at 20:50 +0000, Peter Vanroose wrote:
> > is this (cpoly_roots) computationally efficient, or a lot
> > of overhead (have to solve this very often).
> 
> I guess it's indeed a lot of overhead, certainly for something as simple as square root. Didn't measure this, though.
> 
> 
> > Another idea I had is to convert to std::complex and use
> > the provided sqrt function? Does that make sense?
> 
> Of course. That's probably a lot better. But I'm not familiar enough with std::sqrt(std::complex<T>) to be sure about this; e.g.: which of the two square roots is returned by this function?
> 
> 
> > This brings me to another question about vcl_complex. Is
> > that a wrapper around std::complex?
> 
> Yes. Or better said: it is (in most cases) just a #define for std::complex.
> So the two are (or should be) identical.
> 
> 
> > The vnl:numerics documentation says this:
> > "Compliance with the ANSI standard C++ library 
> > The current draft of the ANSI standard (as at May 1996)
> > includes classes for 1-dimensional vectors (valarray<T>)
> > and complex numbers (complex<T>).
> > There is no standard for matrices.
> > The current vnl classes are not implemented in terms of
> > valarray, as there is a potential performance hit,
> > but in the future they might be."
> >
> > So what's with complex ??
> 
> There is no vnl_complex, just vcl_complex (which is std::complex).
> 
> 
> --      Peter.
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
>

Sorry to bother you guys again, but I found out that matlabs sqrtm is
based on complex schur decomposition (taken from LAPACK routines). That
is more general than the eigenvalue decomposition and can deal with
defective cases.

Is Schur implemented in vnl? Could not really find it, only as part of
solving generalized eigensystems.

Thanks, Martin

On Wed, 2010-12-22 at 17:33 -0500, Martin Reuter wrote:
> Yes, std::sqrt seems to work, only need one sqrt anyway.
> 
> Now I just need a fast way of computing the inverse of a complex 4x4
> matrix. using vnl_inverse gives an error:
> 
>  In function ‘vnl_matrix_fixed<T, 4u, 4u> vnl_inverse(const
> vnl_matrix_fixed<T, 4u, 4u>&) [with T = std::complex<double>]’:
> 
> vnl_inverse.h:114: error: no match for ‘operator==’ in ‘det == 0’
> 
> looks like it cannot deal with complex fixed matrices?
> 
> Martin
> 
> On Wed, 2010-12-22 at 20:50 +0000, Peter Vanroose wrote:
> > > is this (cpoly_roots) computationally efficient, or a lot
> > > of overhead (have to solve this very often).
> > 
> > I guess it's indeed a lot of overhead, certainly for something as simple as square root. Didn't measure this, though.
> > 
> > 
> > > Another idea I had is to convert to std::complex and use
> > > the provided sqrt function? Does that make sense?
> > 
> > Of course. That's probably a lot better. But I'm not familiar enough with std::sqrt(std::complex<T>) to be sure about this; e.g.: which of the two square roots is returned by this function?
> > 
> > 
> > > This brings me to another question about vcl_complex. Is
> > > that a wrapper around std::complex?
> > 
> > Yes. Or better said: it is (in most cases) just a #define for std::complex.
> > So the two are (or should be) identical.
> > 
> > 
> > > The vnl:numerics documentation says this:
> > > "Compliance with the ANSI standard C++ library 
> > > The current draft of the ANSI standard (as at May 1996)
> > > includes classes for 1-dimensional vectors (valarray<T>)
> > > and complex numbers (complex<T>).
> > > There is no standard for matrices.
> > > The current vnl classes are not implemented in terms of
> > > valarray, as there is a potential performance hit,
> > > but in the future they might be."
> > >
> > > So what's with complex ??
> > 
> > There is no vnl_complex, just vcl_complex (which is std::complex).
> > 
> > 
> > --      Peter.
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> 
> 
> ------------------------------------------------------------------------------
> Learn how Oracle Real Application Clusters (RAC) One Node allows customers
> to consolidate database storage, standardize their database environment, and, 
> should the need arise, upgrade to a full multi-node Oracle RAC database 
> without downtime or disruption
> http://p.sf.net/sfu/oracle-sfdevnl
> _______________________________________________
> Vxl-users mailing list
> Vxl-users@...
> https://lists.sourceforge.net/lists/listinfo/vxl-users

Martin Reuter wrote:
>  [...] I found out that matlabs sqrtm is based on
> complex schur decomposition (taken from LAPACK routines).
> Is Schur implemented in vnl?

Yes, but not the complex Schur.
See core/vnl/algo/vnl_generalized_schur.h
According to its documentation, it
   "solves the generalized eigenproblem det(t A - s B) = 0".
Its implementation is based on v3p/netlib/lapack/double/dgges.c
It should not be too difficult to add
   v3p/netlib/lapack/complex16/zgges.c
and rewrite core/vnl/algo/vnl_generalized_schur into
core/vnl/algo/vnl_generalized_complex_schur
to also support the "complex" case.

Contribution to work this out and test it is very much welcomed!

--      Peter.

> Now I just need a fast way of computing the inverse of a
> complex 4x4 matrix. using vnl_inverse gives an error;
> looks like it cannot deal with complex fixed matrices?

That's right.
We have been careful to not #include <complex> all over vnl;
ideally, a variant vnl_complex_inverse.h should be written out, which implements (just) the complex matrix case. Should not be too difficult, I assume.

--      Peter.