- Hi all,

I've just joined this group. I am testing uBLAS and doing some speed

benchmark. In doing that, I found that inner_prod() works ok for

<float> but doesn't compile for std::complex<float>.

// testing inner_prod() for <float>

typedef vector<float>::value_type value_type;

typedef vector<float>::size_type size_type;

typedef type_traits<value_type>::real_type real_type;

vector<float> v1(1), v2(1), v3(1);

real_type f = inner_prod(v1, v2); // ok

// testing inner_prod() for <std::complex<float> >

typedef vector<std::complex<float> >::value_type value_type;

typedef vector<std::complex<float> >::size_type size_type;

typedef type_traits<value_type>::real_type real_type;

vector<std::complex<float> > v1(1), v2(1), v3(1);

real_type f = inner_prod(v1, v2); // type conversion error!

Am I missing anything here, or it is a bug?

Thank you,

Koji - pobdrah wrote:

> I've just joined this group. I am testing uBLAS and doing some speed

[...]

> benchmark. In doing that, I found that inner_prod() works ok for

> <float> but doesn't compile for std::complex<float>.

> // testing inner_prod() for <std::complex<float> >

Inner product is in uBLAS defined as

> typedef vector<std::complex<float> >::value_type value_type;

> typedef vector<std::complex<float> >::size_type size_type;

> typedef type_traits<value_type>::real_type real_type;

> vector<std::complex<float> > v1(1), v2(1), v3(1);

> real_type f = inner_prod(v1, v2); // type conversion error!

>

> Am I missing anything here, or it is a bug?

inner_prod (v1, v2) = sum (v1 [i] * v2 [i]),

that is, as

v1^T v2

both for real and complex vectors -- it is

the equivalent of `[cz]dotu' from BLAS.

The result type of the inner_prod() for two complex

vectors is a complex number, not real; e.g.:

v1 = { (1,-1), (1,-1), (1,-1), (1,-1), (1,-1), (1,-1) }

v2 = { (0,1), (0,1), (0,1), (0,1), (0,1), (0,1) }

inner_prod (v1, v2) = (-6, -6)

inner_prod (v1, v1) = (0, -12)

inner_prod (v2, v2) = (-6, 0)

As the last example shows, imaginary part can sometimes

be zero, but in general it is not.

If you need

v1^H v2

(i.e. `usual' definition of inner product for complex vectors,

i.e. equivalent of `[cz]dotc' from BLAS), you can write

inner_prod (conj (v1), v2)

e.g.:

inner_prod (conj (v1), v2) = (6,-6)

inner_prod (conj (v1), v1) = (12,0)

inner_prod (conj (v2), v2) = (6,0)

but the result is again complex number.

If you calculate v^H v, imaginary part is always zero

(second and third example), so you can say that the result

is real -- but this is special case.

Vector norms are real numbers both for real and complex

vectors; e.g. ||v||_2, that is, ublas::norm_2 (v), is defined as

sqrt (sum (v [i] * conj (v [i])))

i.e.

sqrt (inner_prod (v, conj (v))

Sincerely,

fres

PS. Joerg, there is a small and completely unimportant typo

in `vector_expression.hpp', line 1855: last `)' is missing in the

comment:

// inner_prod (v1, v2) = sum (v1 [i] * v2 [i] > pobdrah wrote:

[snip the excellent explanation]

>

> > I've just joined this group. I am testing uBLAS and doing some speed

> > benchmark. In doing that, I found that inner_prod() works ok for

> > <float> but doesn't compile for std::complex<float>.

> [...]

> > // testing inner_prod() for <std::complex<float> >

> > typedef vector<std::complex<float> >::value_type value_type;

> > typedef vector<std::complex<float> >::size_type size_type;

> > typedef type_traits<value_type>::real_type real_type;

> > vector<std::complex<float> > v1(1), v2(1), v3(1);

> > real_type f = inner_prod(v1, v2); // type conversion error!

> >

> > Am I missing anything here, or it is a bug?

> PS. Joerg, there is a small and completely unimportant typo

Fixed. How did COMO detect this one? ;-)

> in `vector_expression.hpp', line 1855: last `)' is missing in the

> comment:

>

> // inner_prod (v1, v2) = sum (v1 [i] * v2 [i]

Thanks,

Joerg

P.S.:I've just uploaded my current version containing mostly documentation

updates. Due to some bug reports regarding the examples in the documentation

I've now compiled every example and hopefully fixed the corresponding bugs

(BTW should I add the compilable examples to the distribution and to boost

CVS?). I've also changed the tool for editing the HTML documentation: I'm

using Mozilla Composer instead of MS Frontpage now.- Hi Joerg,

>>PS. Joerg, there is a small and completely unimportant typo

Compiler magick ;o) ... with a little help from a human friend.

>>in `vector_expression.hpp', line 1855: last `)' is missing in the

>>comment:

>>

>>// inner_prod (v1, v2) = sum (v1 [i] * v2 [i]

> Fixed. How did COMO detect this one? ;-)

> P.S.:I've just uploaded my current version containing mostly documentation

I think yes -- `examples' in `test?' and `bench?' directories are

> updates. Due to some bug reports regarding the examples in the documentation

> I've now compiled every example and hopefully fixed the corresponding bugs

> (BTW should I add the compilable examples to the distribution and to boost

> CVS?).

not very useful as introductory examples.

BTW, IMHO examples for hermitian matrix/adaptor are wrong

-- if A=A^H, then diagonal elements must have imaginary part

equal to zero [a + bi == a - bi]. Initialization loop can be:

for (int i = 0; i < ha.size1 (); ++ i)

for (int j = 0; j <= i; ++ j)

hal (i, j) = std::complex<double> (3 * i + j, 3 * i - 3 * j);

Regards,

fres - Hi Kresimir,

you wrote:

[snip]

> > P.S.:I've just uploaded my current version containing mostly

documentation

> > updates. Due to some bug reports regarding the examples in the

documentation

> > I've now compiled every example and hopefully fixed the corresponding

bugs

> > (BTW should I add the compilable examples to the distribution and to

boost

> > CVS?).

Done, new download available.

>

> I think yes -- `examples' in `test?' and `bench?' directories are

> not very useful as introductory examples.

> BTW, IMHO examples for hermitian matrix/adaptor are wrong

Corrected.

> -- if A=A^H, then diagonal elements must have imaginary part

> equal to zero [a + bi == a - bi]. Initialization loop can be:

>

> for (int i = 0; i < ha.size1 (); ++ i)

> for (int j = 0; j <= i; ++ j)

> hal (i, j) = std::complex<double> (3 * i + j, 3 * i - 3 * j);

Thanks,

Joerg