Re: Is ther any sense in *= and matrices?

Re: Is ther any sense in *= and matrices?

[tmoran]

To summarize some of the benchmarks of A=A+B vs A+=B, 3x3 matrices, 10**7 iterations
A=A+B  A+=B  ratio
55      19    2.9     Russ's special code on Sunblade, slightly faster than
                      an 866MHz Pentium III according to
                      http://www.phys.ocean.dal.ca/~kelley/linux/benchmark/benchmark_POM.html
4.4      3.3  1.3     MSVC++5.0 simple code on Tom's 866MHz Pentium III
1.4      0.8  1.6     Gnat 3.15p on Tom's 866MHz Pentium III
1.1      0.7  1.6     g++.5.8 (GCC) 3.2.2 using C templates on Hyman's "Sun workstation", speed unspecified
0.45     0.18 2.5     Ole-Hjalmar's machine ? gcc Ada or C the same?

Re: Is ther any sense in *= and matrices?

[Ole Kristensen]

tmoran@acm.org writes:

> To summarize some of the benchmarks of A=A+B vs A+=B, 3x3 matrices, 10**7 iterations
> A=A+B  A+=B  ratio
> 55      19    2.9     Russ's special code on Sunblade, slightly faster than
>                       an 866MHz Pentium III according to
>                       http://www.phys.ocean.dal.ca/~kelley/linux/benchmark/benchmark_POM.html
> 4.4      3.3  1.3     MSVC++5.0 simple code on Tom's 866MHz Pentium III
> 1.4      0.8  1.6     Gnat 3.15p on Tom's 866MHz Pentium III
> 1.1      0.7  1.6     g++.5.8 (GCC) 3.2.2 using C templates on Hyman's "Sun workstation", speed unspecified
> 0.45     0.18 2.5     Ole-Hjalmar's machine ? gcc Ada or C the same?

Your C++ code was compiled with g++ 3.2 under Cygwin, but I did not
set any flags to tell the compiler that my CPU is a Pentium4, so it
may not be optimal. The machine is a 1.8GHz Pentium 4.

Ole-Hj. Kristensen

Re: Is ther any sense in *= and matrices?

[Ole Kristensen]

tmoran@acm.org writes:

> To summarize some of the benchmarks of A=A+B vs A+=B, 3x3 matrices, 10**7 iterations
> A=A+B  A+=B  ratio
> 55      19    2.9     Russ's special code on Sunblade, slightly faster than
>                       an 866MHz Pentium III according to
>                       http://www.phys.ocean.dal.ca/~kelley/linux/benchmark/benchmark_POM.html
> 4.4      3.3  1.3     MSVC++5.0 simple code on Tom's 866MHz Pentium III
> 1.4      0.8  1.6     Gnat 3.15p on Tom's 866MHz Pentium III
> 1.1      0.7  1.6     g++.5.8 (GCC) 3.2.2 using C templates on Hyman's "Sun workstation", speed unspecified
> 0.45     0.18 2.5     Ole-Hjalmar's machine ? gcc Ada or C the same?

Your C++ code was compiled with g++ 3.2 under Cygwin, but I did not
set any flags to tell the compiler that my CPU is a Pentium4, so it
may not be optimal. The machine is a 1.8GHz Pentium 4.

Ole-Hj. Kristensen

Re: Is ther any sense in *= and matrices?

[Ole Kristensen]

tmoran@acm.org writes:

> To summarize some of the benchmarks of A=A+B vs A+=B, 3x3 matrices, 10**7 iterations
> A=A+B  A+=B  ratio
> 55      19    2.9     Russ's special code on Sunblade, slightly faster than
>                       an 866MHz Pentium III according to
>                       http://www.phys.ocean.dal.ca/~kelley/linux/benchmark/benchmark_POM.html
> 4.4      3.3  1.3     MSVC++5.0 simple code on Tom's 866MHz Pentium III
> 1.4      0.8  1.6     Gnat 3.15p on Tom's 866MHz Pentium III
> 1.1      0.7  1.6     g++.5.8 (GCC) 3.2.2 using C templates on Hyman's "Sun workstation", speed unspecified
> 0.45     0.18 2.5     Ole-Hjalmar's machine ? gcc Ada or C the same?

Another data point: My old home PC, a 266MHz pentium II, compiler g++
3.2, Tom's simple C++ code:

ole@linux:~> time ./tmoran
Go 1000000 4460000 6320000
real    0m13.626s
user    0m10.800s
sys     0m0.000s

If I read the results correctly, it actually beats the SunBlade 2000,
which it definitely should not do. But there is much less of a
difference between + and += on this machine than on the fast machine
at work.  Maybe I should check the g++ version on my machine at work
once more.

Turning optimization off produces even less difference:

ole@linux:~> g++ -o tmoran tmoran.cc
ole@linux:~> !time
time ./tmoran
Go 1000000 12850000 14120000
real    0m34.219s
user    0m26.950s
sys     0m0.040s

Ole-Hj. Kristensen

Is ther any sense in *= and matrices?

[Preben Randhol]

Given prior discussion I was thinking yesterday how do one make a:

   *= for matrices in Ada?

I mean: A *= B; or procedure Prod (A,B);

At least if you define your matrix type as:

 type Matrix_Type is array (Positive range <>, Positive range <>) of
       Integer;

I mean even if one define "*" for Matrix_Type one cannot do:

   A := A * B

without getting a Constraint_Error which is reasonable.

-- 
Preben Randhol                    http://www.pvv.org/~randhol/

Re: Is ther any sense in *= and matrices?

[Vinzent Hoefler]

Preben Randhol wrote:

>I mean even if one define "*" for Matrix_Type one cannot do:
>
>   A := A * B

You can, iff A and B are square matrices.


Vinzent.

Re: Is ther any sense in *= and matrices?

[Preben Randhol]

Vinzent Hoefler wrote:
> Preben Randhol wrote:
> 
>>I mean even if one define "*" for Matrix_Type one cannot do:
>>
>>   A := A * B
> 
> You can, iff A and B are square matrices.

Yes, but that is a special case.

-- 
Preben Randhol                    http://www.pvv.org/~randhol/

Re: Is ther any sense in *= and matrices?

[Vinzent Hoefler]

Preben Randhol wrote:

>Vinzent Hoefler wrote:
>> Preben Randhol wrote:
>> 
>>>I mean even if one define "*" for Matrix_Type one cannot do:
>>>
>>>   A := A * B
>> 
>> You can, iff A and B are square matrices.
>
>Yes, but that is a special case.

Of course it is. But because you have to write the Multiply function
yourself anyway and such can take care of that issue, I don't see why
this should be a problem. What is so unresonable about an instance of
a generic matrix multplying function that raises Contraint_Error in
case the index ranges do not match?


Vinzent.

Re: Is ther any sense in *= and matrices?

[Preben Randhol]

Vinzent Hoefler wrote:
> Of course it is. But because you have to write the Multiply function
> yourself anyway and such can take care of that issue, I don't see why
> this should be a problem. 

Because it doesn't make sense mathematically.

-- 
Preben Randhol                    http://www.pvv.org/~randhol/

Re: Is ther any sense in *= and matrices?

[Vinzent Hoefler]

Preben Randhol wrote:

>Vinzent Hoefler wrote:
>> Of course it is. But because you have to write the Multiply function
>> yourself anyway and such can take care of that issue, I don't see why
>> this should be a problem. 
>
>Because it doesn't make sense mathematically.

Dividing by zero doesn't make sense either. So IMO that's what the
Constraint_Error is for.


Vinzent.

Re: Is ther any sense in *= and matrices?

[Preben Randhol]

Vinzent Hoefler wrote:
> Preben Randhol wrote:
> 
>>Vinzent Hoefler wrote:
>>> Of course it is. But because you have to write the Multiply function
>>> yourself anyway and such can take care of that issue, I don't see why
>>> this should be a problem. 
>>
>>Because it doesn't make sense mathematically.
> 
> Dividing by zero doesn't make sense either. So IMO that's what the
> Constraint_Error is for.

The Contraint_Error here is for that *= would return a C (m x n) which
differes from the A (n x n). So *= is only usefule (/= too) in special
cases such as square matrix. Therefore you must make a function Prod
(A,B) at any rate and the *= will be of no use as Prod (A,B) is
general.

-- 
Preben Randhol                    http://www.pvv.org/~randhol/

Re: Is ther any sense in *= and matrices?

[Vinzent Hoefler]

Preben Randhol wrote:

>Vinzent Hoefler wrote:
>> Preben Randhol wrote:
>> 
>>>Vinzent Hoefler wrote:
>>>> Of course it is. But because you have to write the Multiply function
>>>> yourself anyway and such can take care of that issue, I don't see why
>>>> this should be a problem. 
>>>
>>>Because it doesn't make sense mathematically.
>> 
>> Dividing by zero doesn't make sense either. So IMO that's what the
>> Constraint_Error is for.
>
>The Contraint_Error here is for that *= would return a C (m x n) which
>differes from the A (n x n).

Ok, understood.

Depends on how you look at it:

Let's take an extreme example and suppose I have an integer type
defined only for range 0 .. 1. So writing a += would only make sense
for the "special" case that at least one of the operands is zero, not
for the general case, because then the result doesn't fit into the
original operand. But I wouldn't say, writing such a function wouldn't
make sense at all... ;)

In the matrix case it's just that you can't write a *= that works for
the more general case. Simple as that. I really wouldn't see a problem
with that.

>So *= is only usefule (/= too) in special
>cases such as square matrix.

Yes.

>Therefore you must make a function Prod
>(A,B) at any rate and the *= will be of no use as Prod (A,B) is
>general.

No, maybe you want to multiply square matrices only? ;)


Vinzent.

-- 
Parents strongly cautioned  --  this  posting  is  intended for mature
audiences  over  18.  It  may  contain some material that many parents
would not find suitable for children and may include intense violence,
sexual situations, coarse language and suggestive dialogue.

Re: Is ther any sense in *= and matrices?

[Preben Randhol]

Vinzent Hoefler wrote:
> 
> No, maybe you want to multiply square matrices only? ;)

Yes, but as you can either make a procedure or a function of Prod there
is really no need for the *= :-)

-- 
Preben Randhol                    http://www.pvv.org/~randhol/

Re: Is ther any sense in *= and matrices?

[Preben Randhol]

Vinzent Hoefler wrote:
> 
> Let's take an extreme example and suppose I have an integer type
> defined only for range 0 .. 1. So writing a += would only make sense
> for the "special" case that at least one of the operands is zero, not
> for the general case, because then the result doesn't fit into the
> original operand. But I wouldn't say, writing such a function wouldn't
> make sense at all... ;)
> 
> In the matrix case it's just that you can't write a *= that works for
> the more general case. Simple as that. I really wouldn't see a problem
> with that.

Well I don't see the point the *= syntax sugar either ;-)

-- 
Preben Randhol                    http://www.pvv.org/~randhol/

Re: Is ther any sense in *= and matrices?

[Vinzent Hoefler]

Preben Randhol wrote:

>Vinzent Hoefler wrote:
>> 
>> In the matrix case it's just that you can't write a *= that works for
>> the more general case. Simple as that. I really wouldn't see a problem
>> with that.
>
>Well I don't see the point the *= syntax sugar either ;-)

I agree. It was always confusing me from the first time I looked at
some C-source. =)

And especially I had a hard time translating an LZSS-implementation
from C to Pascal that didn't know all these fscking post- and
preincrements in indexing some array structures in a single line.

That's why I like Ada so much. All this obfuscating stuff simply
doesn't exist.


Vinzent.

-- 
Parents strongly cautioned  --  this  posting  is  intended for mature
audiences  over  18.  It  may  contain some material that many parents
would not find suitable for children and may include intense violence,
sexual situations, coarse language and suggestive dialogue.

Re: Is ther any sense in *= and matrices?

[Russ]

Preben Randhol <randhol+abuse@pvv.org> wrote in message news:<slrnbdupg6.88r.randhol+abuse@kiuk0152.chembio.ntnu.no>...
> Vinzent Hoefler wrote:
> > Preben Randhol wrote:
> > 
> >>Vinzent Hoefler wrote:
> >>> Of course it is. But because you have to write the Multiply function
> >>> yourself anyway and such can take care of that issue, I don't see why
> >>> this should be a problem. 
> >>
> >>Because it doesn't make sense mathematically.
> > 
> > Dividing by zero doesn't make sense either. So IMO that's what the
> > Constraint_Error is for.
> 
> The Contraint_Error here is for that *= would return a C (m x n) which
> differes from the A (n x n). So *= is only usefule (/= too) in special
> cases such as square matrix. Therefore you must make a function Prod
> (A,B) at any rate and the *= will be of no use as Prod (A,B) is
> general.

"*=" is not so useful for multipying non-square matrices. However, it
is very useful for multiplying a matrix (or a vector) by a scalar. For
example,

    A *= 2.0

can be used to multiply matrix A by two in place. This is potentially
more efficient than A := A * 2 for all the same reasons discussed with
respect to "+=".

Re: Is ther any sense in *= and matrices?

[John R. Strohm]

"Russ" <18k11tm001@sneakemail.com> wrote in message
news:bebbba07.0306071138.6bf9784f@posting.google.com...
> Preben Randhol <randhol+abuse@pvv.org> wrote in message
news:<slrnbdupg6.88r.randhol+abuse@kiuk0152.chembio.ntnu.no>...
> > Vinzent Hoefler wrote:
> > > Preben Randhol wrote:
> > >
> > >>Vinzent Hoefler wrote:
> > >>> Of course it is. But because you have to write the Multiply function
> > >>> yourself anyway and such can take care of that issue, I don't see
why
> > >>> this should be a problem.
> > >>
> > >>Because it doesn't make sense mathematically.
> > >
> > > Dividing by zero doesn't make sense either. So IMO that's what the
> > > Constraint_Error is for.
> >
> > The Contraint_Error here is for that *= would return a C (m x n) which
> > differes from the A (n x n). So *= is only usefule (/= too) in special
> > cases such as square matrix. Therefore you must make a function Prod
> > (A,B) at any rate and the *= will be of no use as Prod (A,B) is
> > general.
>
> "*=" is not so useful for multipying non-square matrices. However, it
> is very useful for multiplying a matrix (or a vector) by a scalar. For
> example,
>
>     A *= 2.0
>
> can be used to multiply matrix A by two in place. This is potentially
> more efficient than A := A * 2 for all the same reasons discussed with
> respect to "+=".

With all due respect, ladies and gentlemen, it has been known for a very
long time that the difference in "efficiency" between A := A + B and A += B
is lost in the noise floor compared to the improvements that can be gotten
by improving the algorithms involved.

And I'd be REALLY interested to know what you are doing with matrix
multiplication such that the product of a matrix and a scalar is the
slightest bit interesting.  (Usually, when I'm multiplying matrices, I'm
playing with direction cosine matrices, and I don't recall ever hearing
about any use for multiplying a DCM by a scalar.)

Re: Is ther any sense in *= and matrices?

[Russ]

"John R. Strohm" <strohm@airmail.net> wrote in message news:<bbumef$fnq@library2.airnews.net>...

> "Russ" <18k11tm001@sneakemail.com> wrote in message
> news:bebbba07.0306071138.6bf9784f@posting.google.com...

> > "*=" is not so useful for multipying non-square matrices. However, it
> > is very useful for multiplying a matrix (or a vector) by a scalar. For
> > example,
> >
> >     A *= 2.0
> >
> > can be used to multiply matrix A by two in place. This is potentially
> > more efficient than A := A * 2 for all the same reasons discussed with
> > respect to "+=".
> 
> With all due respect, ladies and gentlemen, it has been known for a very
> long time that the difference in "efficiency" between A := A + B and A += B
> is lost in the noise floor compared to the improvements that can be gotten
> by improving the algorithms involved.

Oh, really? I just did a test in C++ with 3x3 matrices. I added them
together 10,000,000 times using "+", then "+=". The "+=" version took
about 19 seconds, and the "+" version took about 55 seconds. That's
just shy of a factor of 3, folks. If that's your "noise floor," I
can't help wonder what kind of "algorithms" you are dealing with!

> And I'd be REALLY interested to know what you are doing with matrix
> multiplication such that the product of a matrix and a scalar is the
> slightest bit interesting.  (Usually, when I'm multiplying matrices, I'm
> playing with direction cosine matrices, and I don't recall ever hearing
> about any use for multiplying a DCM by a scalar.)

Sorry, but I'm having some trouble with your "logic". DCM matrices
don't get multiplied by scalars, so you doubt that any vector or
matrix ever gets multiplied by a scalar?

Actually, it is more common to scale a vector than a matrix, but the
same principles regarding the efficiency of "+"and "+=" still apply.
Vectors often need to be scaled for conversion of units, for example.
Another very common example would be determining a change in position
by multiplying a velocity vector by a scalar time increment. Ditto for
determining a change in a velocity vector by multiplying an
acceleration vector by a time increment.

Since you know about DCMs, you must also be familiar with quaternions,
otherwise known as Euler parameters. Like DCMs, Euler parameter
"vectors" represents the attitude of a rigid body. Well, an Euler
parameter vector consists of 4 scalar components or "Euler
parameters", and the norm (roor sum square) of the 4 parameters is
unity by definition. However, due to numerical roundoff error, the
norm can drift away from one. However, you can renormalize by dividing
the vector by its own norm. This can be written very concisely and
efficiently in C++ as

    eulpar /= eulpar.norm()

Incidentally, this procedure is much simpler than the corresponding
procedure for orthogonalizing a DCM.

I could come up with other examples, but this should suffice.
Actually, I'm a bit dissappointed in myself for wasting as much time
as I did in replying to your meritless post. So is my wife.

Re: Is ther any sense in *= and matrices?

[tmoran]

>Oh, really? I just did a test in C++ with 3x3 matrices. I added them
>together 10,000,000 times using "+", then "+=". The "+=" version took
>about 19 seconds, and the "+" version took about 55 seconds. That's
  Would you be so kind as to post your code, and what C++ compiler
and what hardware you used?  Your results seem quite different from
other people's.  My old 900MHz Windows 2K machine using MSVC++ 5.0
took 4.38 and 3.28 seconds, a factor of 1.33, and using Gnat 3.15p
on the same machine took 1.38 and 0.85 seconds, a ratio of 1.6
Clearly, there's something substantially different between our
compilers/hardware/code.

Re: Is ther any sense in *= and matrices?

[Russ]

tmoran@acm.org wrote in message news:<AmNEa.67177$d51.133670@sccrnsc01>...
> >Oh, really? I just did a test in C++ with 3x3 matrices. I added them
> >together 10,000,000 times using "+", then "+=". The "+=" version took
> >about 19 seconds, and the "+" version took about 55 seconds. That's
>   Would you be so kind as to post your code, and what C++ compiler
> and what hardware you used?  Your results seem quite different from
> other people's.  My old 900MHz Windows 2K machine using MSVC++ 5.0
> took 4.38 and 3.28 seconds, a factor of 1.33, and using Gnat 3.15p
> on the same machine took 1.38 and 0.85 seconds, a ratio of 1.6
> Clearly, there's something substantially different between our
> compilers/hardware/code.

I'm using gcc 2.95.2 on a Sunblade 2000. I can't post the code, but it
is a pretty standard vector/matrix implementation in C++. Actually, it
is designed for very efficient indexing, perhaps at the expense of
slightly less efficient construction (it has a pointer for each row of
the matrix). That might explain part of the difference you are seeing,
but certainly not all. Perhaps your choice of C++ compiler is a factor
too.

Re: Is ther any sense in *= and matrices?

[John R. Strohm]

"Russ" <18k11tm001@sneakemail.com> wrote in message
news:bebbba07.0306082024.7cebb5df@posting.google.com...
> tmoran@acm.org wrote in message news:<AmNEa.67177$d51.133670@sccrnsc01>...
> > >Oh, really? I just did a test in C++ with 3x3 matrices. I added them
> > >together 10,000,000 times using "+", then "+=". The "+=" version took
> > >about 19 seconds, and the "+" version took about 55 seconds. That's
> >   Would you be so kind as to post your code, and what C++ compiler
> > and what hardware you used?  Your results seem quite different from
> > other people's.  My old 900MHz Windows 2K machine using MSVC++ 5.0
> > took 4.38 and 3.28 seconds, a factor of 1.33, and using Gnat 3.15p
> > on the same machine took 1.38 and 0.85 seconds, a ratio of 1.6
> > Clearly, there's something substantially different between our
> > compilers/hardware/code.
>
> I'm using gcc 2.95.2 on a Sunblade 2000. I can't post the code, but it
> is a pretty standard vector/matrix implementation in C++. Actually, it
> is designed for very efficient indexing, perhaps at the expense of
> slightly less efficient construction (it has a pointer for each row of
> the matrix). That might explain part of the difference you are seeing,
> but certainly not all. Perhaps your choice of C++ compiler is a factor
> too.

I think the point that Tom is trying to make is that your results are so far
out of line with expected reality that he suspects that your optimizations
for indexing may in fact be pessimizations.

Without seeing your code, we aren't going to be able to tell very much.

There is also this.  You're on a Sunblade.  I'm guessing that is a PowerPC.
PowerPC cache effects can do horrendous things to your runtime, and very
minor changes in code structure, that affect the cache hit pattern, can have
VERY disproportionate effects.  (It is not at all hard to see an order of
magnitude slowdown.)  Mercury Computer went to a great deal of trouble to
show this to us, and teach us how to exploit it properly, in the RACEway
programming class.  (And I'm the guy who took one of their examples, that
they believed was as good as it was possible to get it, and got another 5%
speedup out of it.  They were rather surprised.)

Re: Is ther any sense in *= and matrices?

[Ole-Hjalmar Kristensen]

"John R. Strohm" <strohm@airmail.net> writes:

> "Russ" <18k11tm001@sneakemail.com> wrote in message
> news:bebbba07.0306082024.7cebb5df@posting.google.com...
> > tmoran@acm.org wrote in message news:<AmNEa.67177$d51.133670@sccrnsc01>...
> > > >Oh, really? I just did a test in C++ with 3x3 matrices. I added them
> > > >together 10,000,000 times using "+", then "+=". The "+=" version took
> > > >about 19 seconds, and the "+" version took about 55 seconds. That's
> > >   Would you be so kind as to post your code, and what C++ compiler
> > > and what hardware you used?  Your results seem quite different from
> > > other people's.  My old 900MHz Windows 2K machine using MSVC++ 5.0
> > > took 4.38 and 3.28 seconds, a factor of 1.33, and using Gnat 3.15p
> > > on the same machine took 1.38 and 0.85 seconds, a ratio of 1.6
> > > Clearly, there's something substantially different between our
> > > compilers/hardware/code.
> >
> > I'm using gcc 2.95.2 on a Sunblade 2000. I can't post the code, but it
> > is a pretty standard vector/matrix implementation in C++. Actually, it
> > is designed for very efficient indexing, perhaps at the expense of
> > slightly less efficient construction (it has a pointer for each row of
> > the matrix). That might explain part of the difference you are seeing,
> > but certainly not all. Perhaps your choice of C++ compiler is a factor
> > too.
> 
> I think the point that Tom is trying to make is that your results are so far
> out of line with expected reality that he suspects that your optimizations
> for indexing may in fact be pessimizations.
> 

Yes, on modern hardware, having a pointer for each row of the matrix
usually leads to slower code. The reason is that the indexing
operation by multilply/add is actually faster than fetching a pointer
from an array and then adding, thereby incurring an extra memory fetch
and cache pollution. The above holds for large matrices, I haven't
investigated the effects on small matrices.

> Without seeing your code, we aren't going to be able to tell very much.
> 
> There is also this.  You're on a Sunblade.  I'm guessing that is a PowerPC.

It's a SPARC.

 <snip>

-- 
Ole-Hj. Kristensen

******************************************************************************
* This page intentionally left blank.
******************************************************************************

Re: Is ther any sense in *= and matrices?

[Wesley Groleau]

> Yes, on modern hardware, having a pointer for each row of the matrix
> usually leads to slower code. The reason is that the indexing
> operation by multilply/add is actually faster than fetching a pointer
> from an array and then adding, thereby incurring an extra memory fetch
> and cache pollution. The above holds for large matrices, I haven't

So Java's handling of arrays is inefficient
for matrices, even if bypassing the JVM by
compiling to native code?

Re: Is ther any sense in *= and matrices?

[Ole Kristensen]

Wesley Groleau <wesgroleau@despammed.com> writes:

> > Yes, on modern hardware, having a pointer for each row of the matrix
> > usually leads to slower code. The reason is that the indexing
> > operation by multilply/add is actually faster than fetching a pointer
> > from an array and then adding, thereby incurring an extra memory fetch
> > and cache pollution. The above holds for large matrices, I haven't
> 
> So Java's handling of arrays is inefficient
> for matrices, even if bypassing the JVM by
> compiling to native code?

Can't speak for Java, but when I tried a 128**3 matrix in Ada and C++,
the plain indexing scheme beat the clever pointer array by a
significant amount. This was on a 1.8GHz Pentium 4. I cannot say
anything definite about other systems, but in general, it pays to
think a lot about memory access patterns on fast CPU's.

Re: Is ther any sense in *= and matrices?

[Russ]

Ole-Hjalmar Kristensen <oleh@vlinux.voxelvision.no> wrote in message news:<7visre2rwv.fsf@vlinux.voxelvision.no>...
> "John R. Strohm" <strohm@airmail.net> writes:
> 
> > "Russ" <18k11tm001@sneakemail.com> wrote in message
> > news:bebbba07.0306082024.7cebb5df@posting.google.com...
> > > tmoran@acm.org wrote in message news:<AmNEa.67177$d51.133670@sccrnsc01>...
> > > > >Oh, really? I just did a test in C++ with 3x3 matrices. I added them
> > > > >together 10,000,000 times using "+", then "+=". The "+=" version took
> > > > >about 19 seconds, and the "+" version took about 55 seconds. That's
> > > >   Would you be so kind as to post your code, and what C++ compiler
> > > > and what hardware you used?  Your results seem quite different from
> > > > other people's.  My old 900MHz Windows 2K machine using MSVC++ 5.0
> > > > took 4.38 and 3.28 seconds, a factor of 1.33, and using Gnat 3.15p
> > > > on the same machine took 1.38 and 0.85 seconds, a ratio of 1.6
> > > > Clearly, there's something substantially different between our
> > > > compilers/hardware/code.
> > >
> > > I'm using gcc 2.95.2 on a Sunblade 2000. I can't post the code, but it
> > > is a pretty standard vector/matrix implementation in C++. Actually, it
> > > is designed for very efficient indexing, perhaps at the expense of
> > > slightly less efficient construction (it has a pointer for each row of
> > > the matrix). That might explain part of the difference you are seeing,
> > > but certainly not all. Perhaps your choice of C++ compiler is a factor
> > > too.
> > 
> > I think the point that Tom is trying to make is that your results are so far
> > out of line with expected reality that he suspects that your optimizations
> > for indexing may in fact be pessimizations.
> > 
> 
> Yes, on modern hardware, having a pointer for each row of the matrix
> usually leads to slower code. The reason is that the indexing
> operation by multilply/add is actually faster than fetching a pointer
> from an array and then adding, thereby incurring an extra memory fetch
> and cache pollution. The above holds for large matrices, I haven't
> investigated the effects on small matrices.

Thanks for clarifying that. I wrote that code about 10 years ago when
I was just starting to learn C++. It has other inefficiencies too. For
example, I used offset pointers to achieve indexing that starts with 1
rather than 0. This saved the "minus one" operation on each indexing,
but it put the "minus one" operations into the constructor. It makes
sense to do the "minus one" only once in the constructor, but it does
slow down the constructor (which increases the cost of a temporary).

If I had it to do over again, I'd probably just settle for indexing
that starts with zero. Ada allows offset indexing much more naturally
than C or C++, of course, but that's another topic.

I am certainly no efficiency expert, but I do believe that an basic
step in achieving reasonable efficiency is to eliminate excessive
generation of temporary objects.

Re: Is ther any sense in *= and matrices?

[John R. Strohm]

"Russ" <18k11tm001@sneakemail.com> wrote in message
news:bebbba07.0306101033.66054a83@posting.google.com...
> Ole-Hjalmar Kristensen <oleh@vlinux.voxelvision.no> wrote in message
news:<7visre2rwv.fsf@vlinux.voxelvision.no>...
> > "John R. Strohm" <strohm@airmail.net> writes:
> >
> > > "Russ" <18k11tm001@sneakemail.com> wrote in message
> > > news:bebbba07.0306082024.7cebb5df@posting.google.com...
> > > > tmoran@acm.org wrote in message
news:<AmNEa.67177$d51.133670@sccrnsc01>...
> > > > > >Oh, really? I just did a test in C++ with 3x3 matrices. I added
them
> > > > > >together 10,000,000 times using "+", then "+=". The "+=" version
took
> > > > > >about 19 seconds, and the "+" version took about 55 seconds.
That's
> > > > >   Would you be so kind as to post your code, and what C++ compiler
> > > > > and what hardware you used?  Your results seem quite different
from
> > > > > other people's.  My old 900MHz Windows 2K machine using MSVC++ 5.0
> > > > > took 4.38 and 3.28 seconds, a factor of 1.33, and using Gnat 3.15p
> > > > > on the same machine took 1.38 and 0.85 seconds, a ratio of 1.6
> > > > > Clearly, there's something substantially different between our
> > > > > compilers/hardware/code.
> > > >
> > > > I'm using gcc 2.95.2 on a Sunblade 2000. I can't post the code, but
it
> > > > is a pretty standard vector/matrix implementation in C++. Actually,
it
> > > > is designed for very efficient indexing, perhaps at the expense of
> > > > slightly less efficient construction (it has a pointer for each row
of
> > > > the matrix). That might explain part of the difference you are
seeing,
> > > > but certainly not all. Perhaps your choice of C++ compiler is a
factor
> > > > too.
> > >
> > > I think the point that Tom is trying to make is that your results are
so far
> > > out of line with expected reality that he suspects that your
optimizations
> > > for indexing may in fact be pessimizations.
> > >
> >
> > Yes, on modern hardware, having a pointer for each row of the matrix
> > usually leads to slower code. The reason is that the indexing
> > operation by multilply/add is actually faster than fetching a pointer
> > from an array and then adding, thereby incurring an extra memory fetch
> > and cache pollution. The above holds for large matrices, I haven't
> > investigated the effects on small matrices.
>
> Thanks for clarifying that. I wrote that code about 10 years ago when
> I was just starting to learn C++. It has other inefficiencies too. For
> example, I used offset pointers to achieve indexing that starts with 1
> rather than 0. This saved the "minus one" operation on each indexing,
> but it put the "minus one" operations into the constructor. It makes
> sense to do the "minus one" only once in the constructor, but it does
> slow down the constructor (which increases the cost of a temporary).

Sounds like you really ought to rewrite your code and see what happens to
your benchmark results.

> If I had it to do over again, I'd probably just settle for indexing
> that starts with zero. Ada allows offset indexing much more naturally
> than C or C++, of course, but that's another topic.
>
> I am certainly no efficiency expert, but I do believe that an basic
> step in achieving reasonable efficiency is to eliminate excessive
> generation of temporary objects.

It really depends on the temporary object in question.  For a 3x3 matrix
multiply (or a 3x3 matrix add), on a machine with LOTS of general-purpose
registers, the optimum temporary object is nine CPU registers.  ZERO
allocation cost, ZERO deallocation cost, ZERO impact.

Re: Is ther any sense in *= and matrices?

[tmoran]

> > >about 19 seconds, and the "+" version took about 55 seconds. That's
> > took 4.38 and 3.28 seconds, a factor of 1.33, and using Gnat 3.15p
> I'm using gcc 2.95.2 on a Sunblade 2000. I can't post the code,
  Here's my code.  Please try it on the Sunblade 2000 with gcc 2.95.2
and let us know the results.
#include <stdio.h>
#include <time.h>
#define N 3
        // 4.386/3.275= 1.33              10M            N=3

struct matrices{int x[N][N];};

struct matrices f(struct matrices left, struct matrices right){
        struct matrices result;
        int i,j;
        for(i=0;i<N;i++)
                for(j=0;j<N;j++)
                        result.x[i][j]=left.x[i][j]+right.x[i][j];
    return result;
}

void add_to(struct matrices left, struct matrices right){
        int i,j;
        for(i=0;i<N;i++)
                for(j=0;j<N;j++)
                        left.x[i][j]+=right.x[i][j];
}

void set(struct matrices m, int v){
        int i,j;
        for(i=0;i<N;i++)
                for(j=0;j<N;j++)
                        m.x[i][j]=v;
}

void main() {
        int i;
        clock_t t0,t1,t2;
        matrices a;
        matrices b;
        printf("Go ");
        set(a,0);
        set(b,0);
        add_to(a,b);
        a = f(a,b);
        t0 = clock();
        for(i=0;i<10000000;i++)
          add_to(a,b);
        t1 = clock();
        for(i=0;i<10000000;i++)
          a = f(a,b);
        t2=clock();
        printf("%i %i %i",(int)CLOCKS_PER_SEC,(int)(t1-t0),(int)(t2-t1));
}

Re: Is ther any sense in *= and matrices?

[John R. Strohm]

"Russ" <18k11tm001@sneakemail.com> wrote in message
news:bebbba07.0306081051.5f3ac24b@posting.google.com...
> "John R. Strohm" <strohm@airmail.net> wrote in message
news:<bbumef$fnq@library2.airnews.net>...
>
> > "Russ" <18k11tm001@sneakemail.com> wrote in message
> > news:bebbba07.0306071138.6bf9784f@posting.google.com...
>
> > > "*=" is not so useful for multipying non-square matrices. However, it
> > > is very useful for multiplying a matrix (or a vector) by a scalar. For
> > > example,
> > >
> > >     A *= 2.0
> > >
> > > can be used to multiply matrix A by two in place. This is potentially
> > > more efficient than A := A * 2 for all the same reasons discussed with
> > > respect to "+=".
> >
> > With all due respect, ladies and gentlemen, it has been known for a very
> > long time that the difference in "efficiency" between A := A + B and A
+= B
> > is lost in the noise floor compared to the improvements that can be
gotten
> > by improving the algorithms involved.
>
> Oh, really? I just did a test in C++ with 3x3 matrices. I added them
> together 10,000,000 times using "+", then "+=". The "+=" version took
> about 19 seconds, and the "+" version took about 55 seconds. That's
> just shy of a factor of 3, folks. If that's your "noise floor," I
> can't help wonder what kind of "algorithms" you are dealing with!

19 seconds divided by 10,000,000 is 1.9 microseconds per matrix addition.
55 seconds divided by 10,000,000 is 5.5 microseconds per matrix addition.  I
find it VERY difficult to believe that your total processing chain is going
to be such that 3 microseconds per matrix addition is a DOMINANT part of the
timeline.

In other words, if your TOTAL timeline is, say, 50 microseconds, 3
microseconds isn't going to buy you that much.  I'm still wondering what you
are doing that makes that 3 usec interesting.

> > And I'd be REALLY interested to know what you are doing with matrix
> > multiplication such that the product of a matrix and a scalar is the
> > slightest bit interesting.  (Usually, when I'm multiplying matrices, I'm
> > playing with direction cosine matrices, and I don't recall ever hearing
> > about any use for multiplying a DCM by a scalar.)
>
> Sorry, but I'm having some trouble with your "logic". DCM matrices
> don't get multiplied by scalars, so you doubt that any vector or
> matrix ever gets multiplied by a scalar?

Since you want to get nasty personal about it...

I am fully aware that VECTOR * scalar happens, frequently.  I was asking for
an application that did MATRIX * scalar.

> Actually, it is more common to scale a vector than a matrix, but the
> same principles regarding the efficiency of "+"and "+=" still apply.
> Vectors often need to be scaled for conversion of units, for example.
> Another very common example would be determining a change in position
> by multiplying a velocity vector by a scalar time increment. Ditto for
> determining a change in a velocity vector by multiplying an
> acceleration vector by a time increment.

Determining a change in a position vector by multiplying a velocity vector
by a scalar time increment, or determining a change in a velocity vector by
multiplying an acceleration vector by a time increment, just like that, is
Euler's Method for numerical integration.  It is about the worst method
known to man, BOTH in terms of processing time AND error buildup (and
demonstrating this fact is a standard assignment in the undergraduate
numerical methods survey course).  If you are using Euler's Method for a
high-performance sim, you are out of your mind.  Runge-Kutta methods offer
orders of magnitude better performance, in time AND error control.

And don't change the subject.  We are talking about MATRIX * scalar, not
vector * scalar.  Quaternions are vectors, not matrices.

> Since you know about DCMs, you must also be familiar with quaternions,
> otherwise known as Euler parameters. Like DCMs, Euler parameter
> "vectors" represents the attitude of a rigid body. Well, an Euler
> parameter vector consists of 4 scalar components or "Euler
> parameters", and the norm (roor sum square) of the 4 parameters is
> unity by definition. However, due to numerical roundoff error, the
> norm can drift away from one. However, you can renormalize by dividing
> the vector by its own norm. This can be written very concisely and
> efficiently in C++ as
>
>     eulpar /= eulpar.norm()

Concisely, yes.  Efficiently, not necessarily.  You still have to compute
the norm, which involves taking a square root, which, off the top of my
head, probably costs quite a bit more than four divisions.  If you are
really into microefficiency, you would code this as renorm(eulpar), or
eulpar.renorm(), and do it in-place, computing 1/sqrt(norm) directly and
multplying, rather than computing sqrt() and dividing.

And all that is irrelevant to the original question, because quaternion
renormalization is VECTOR * scalar, not MATRIX * scalar.

> Incidentally, this procedure is much simpler than the corresponding
> procedure for orthogonalizing a DCM.

No argument on the microissue.  I haven't looked at quaternion methods in a
few years, so I don't recall offhand what all is involved in building a sim
around quaternions as opposed to DCMs.  My recollection is that you still
need the DCMs, and you still have to compute between quaternions and DCMs.
You use quaternions because (a) you can interpolate/integrate them easily
(good) and (b) they don't suffer from gimbal lock (much better in some
applications).  You pay a price for that flexibility.

> I could come up with other examples, but this should suffice.
> Actually, I'm a bit dissappointed in myself for wasting as much time
> as I did in replying to your meritless post. So is my wife.

Re: Is ther any sense in *= and matrices?

[Russ]

"John R. Strohm" <strohm@airmail.net> wrote in message news:<bc0dql$qnh@library1.airnews.net>...
> "Russ" <18k11tm001@sneakemail.com> wrote in message
> news:bebbba07.0306081051.5f3ac24b@posting.google.com...
> > "John R. Strohm" <strohm@airmail.net> wrote in message
>  news:<bbumef$fnq@library2.airnews.net>...
> >
> > > "Russ" <18k11tm001@sneakemail.com> wrote in message
> > > news:bebbba07.0306071138.6bf9784f@posting.google.com...
>  
> > > > "*=" is not so useful for multipying non-square matrices. However, it
> > > > is very useful for multiplying a matrix (or a vector) by a scalar. For
> > > > example,
> > > >
> > > >     A *= 2.0
> > > >
> > > > can be used to multiply matrix A by two in place. This is potentially
> > > > more efficient than A := A * 2 for all the same reasons discussed with
> > > > respect to "+=".
> > >
> > > With all due respect, ladies and gentlemen, it has been known for a very
> > > long time that the difference in "efficiency" between A := A + B and A
>  += B
> > > is lost in the noise floor compared to the improvements that can be
>  gotten
> > > by improving the algorithms involved.
> >
> > Oh, really? I just did a test in C++ with 3x3 matrices. I added them
> > together 10,000,000 times using "+", then "+=". The "+=" version took
> > about 19 seconds, and the "+" version took about 55 seconds. That's
> > just shy of a factor of 3, folks. If that's your "noise floor," I
> > can't help wonder what kind of "algorithms" you are dealing with!
> 
> 19 seconds divided by 10,000,000 is 1.9 microseconds per matrix addition.
> 55 seconds divided by 10,000,000 is 5.5 microseconds per matrix addition.  I
> find it VERY difficult to believe that your total processing chain is going
> to be such that 3 microseconds per matrix addition is a DOMINANT part of the
> timeline.
> 
> In other words, if your TOTAL timeline is, say, 50 microseconds, 3
> microseconds isn't going to buy you that much.  I'm still wondering what you
> are doing that makes that 3 usec interesting.

Ten years ago I was using my C++ vector/matrix code for an
experimental, real-time GPS/INS precision landing and guidance system,
including a 16-state Kalman filter (if I recall correctly). We had a
dozen or so tasks, and the guidance update rate was 20 Hz. It was a
fairly ambitious project at the time, and it worked out well. It might
NOT have worked had we been generating temporary vectors and matrices
at a high rate.

But that's all beside the point. Who are you to question whether a
certain level of efficiency is enough? Just because *you* don't need
high efficiency for your work, that doesn't mean that *nobody* needs
it. And by the way, I'm not talking about ultra-efficiency here. I'm
just talking about the basics. You could write the slickest algorithm
ever written, but if you generate temporary vectors and matrices at a
high rate, it's like running a 100-yard dash carrying a freakin'
bowling ball.

> > > And I'd be REALLY interested to know what you are doing with matrix
> > > multiplication such that the product of a matrix and a scalar is the
> > > slightest bit interesting.  (Usually, when I'm multiplying matrices, I'm
> > > playing with direction cosine matrices, and I don't recall ever hearing
> > > about any use for multiplying a DCM by a scalar.)
> >
> > Sorry, but I'm having some trouble with your "logic". DCM matrices
> > don't get multiplied by scalars, so you doubt that any vector or
> > matrix ever gets multiplied by a scalar?
> 
> Since you want to get nasty personal about it...
> 
> I am fully aware that VECTOR * scalar happens, frequently.  I was asking for
> an application that did MATRIX * scalar.

Well, first of all, a vector is really just a matrix with one column.
More importantly, as I wrote earlier (and which you apparently
missed), all the efficiency concerns regarding "=" and "+=" are
EXACTLY the same for both vectors *and* matrices. Not one iota of
difference. So it really doesn't matter whether we are talking about
vectors or matrices. The fact that I referred to matrices in my
original post is irrrelevant.

> > Actually, it is more common to scale a vector than a matrix, but the
> > same principles regarding the efficiency of "+"and "+=" still apply.
> > Vectors often need to be scaled for conversion of units, for example.
> > Another very common example would be determining a change in position
> > by multiplying a velocity vector by a scalar time increment. Ditto for
> > determining a change in a velocity vector by multiplying an
> > acceleration vector by a time increment.
> 
> Determining a change in a position vector by multiplying a velocity vector
> by a scalar time increment, or determining a change in a velocity vector by
> multiplying an acceleration vector by a time increment, just like that, is
> Euler's Method for numerical integration.  It is about the worst method
> known to man, BOTH in terms of processing time AND error buildup (and
> demonstrating this fact is a standard assignment in the undergraduate
> numerical methods survey course).  If you are using Euler's Method for a
> high-performance sim, you are out of your mind.  Runge-Kutta methods offer
> orders of magnitude better performance, in time AND error control.

Oh, good lord. You still need to multiply velocity by time, regardless
of which method of numerical integration you are using. In fact, you
need to do it more with the higher-order methods.

> And don't change the subject.  We are talking about MATRIX * scalar, not
> vector * scalar.  Quaternions are vectors, not matrices.

No, *you* are talking only about multiplying matrices by scalars, but
I can't understand why, since the principle in question applies
equally to both vectors *and* matrices.

> > Since you know about DCMs, you must also be familiar with quaternions,
> > otherwise known as Euler parameters. Like DCMs, Euler parameter
> > "vectors" represents the attitude of a rigid body. Well, an Euler
> > parameter vector consists of 4 scalar components or "Euler
> > parameters", and the norm (roor sum square) of the 4 parameters is
> > unity by definition. However, due to numerical roundoff error, the
> > norm can drift away from one. However, you can renormalize by dividing
> > the vector by its own norm. This can be written very concisely and
> > efficiently in C++ as
> >
> >     eulpar /= eulpar.norm()
> 
> Concisely, yes.  Efficiently, not necessarily.  You still have to compute
> the norm, which involves taking a square root, which, off the top of my
> head, probably costs quite a bit more than four divisions.  If you are
> really into microefficiency, you would code this as renorm(eulpar), or
> eulpar.renorm(), and do it in-place, computing 1/sqrt(norm) directly and
> multplying, rather than computing sqrt() and dividing.

Taking a square root is ultra-fast on modern processors, but that is
irrelevant because there is no way around it for normalizing a
quaternion. And, no, I'm not into micro-efficiency at all. However,
you may have a point: if multiplicatioin is slightly more efficient
than division (is it?), then it *would* make sense to do the division
once and then do multiplication after that. But *that's*
micro-efficiency.

Re: Is ther any sense in *= and matrices?

[Mango Jones]

That reminds me of a time at www.spacemag.co.uk

18k11tm001@sneakemail.com (Russ) wrote in message news:<bebbba07.0306082206.69603145@posting.google.com>...
> "John R. Strohm" <strohm@airmail.net> wrote in message news:<bc0dql$qnh@library1.airnews.net>...
> > "Russ" <18k11tm001@sneakemail.com> wrote in message
> > news:bebbba07.0306081051.5f3ac24b@posting.google.com...
> > > "John R. Strohm" <strohm@airmail.net> wrote in message
>  news:<bbumef$fnq@library2.airnews.net>...
> > >
> > > > "Russ" <18k11tm001@sneakemail.com> wrote in message
> > > > news:bebbba07.0306071138.6bf9784f@posting.google.com...
>  
> > > > > "*=" is not so useful for multipying non-square matrices. However, it
> > > > > is very useful for multiplying a matrix (or a vector) by a scalar. For
> > > > > example,
> > > > >
> > > > >     A *= 2.0
> > > > >
> > > > > can be used to multiply matrix A by two in place. This is potentially
> > > > > more efficient than A := A * 2 for all the same reasons discussed with
> > > > > respect to "+=".
> > > >
> > > > With all due respect, ladies and gentlemen, it has been known for a very
> > > > long time that the difference in "efficiency" between A := A + B and A
>  += B
> > > > is lost in the noise floor compared to the improvements that can be
>  gotten
> > > > by improving the algorithms involved.
> > >
> > > Oh, really? I just did a test in C++ with 3x3 matrices. I added them
> > > together 10,000,000 times using "+", then "+=". The "+=" version took
> > > about 19 seconds, and the "+" version took about 55 seconds. That's
> > > just shy of a factor of 3, folks. If that's your "noise floor," I
> > > can't help wonder what kind of "algorithms" you are dealing with!
> > 
> > 19 seconds divided by 10,000,000 is 1.9 microseconds per matrix addition.
> > 55 seconds divided by 10,000,000 is 5.5 microseconds per matrix addition.  I
> > find it VERY difficult to believe that your total processing chain is going
> > to be such that 3 microseconds per matrix addition is a DOMINANT part of the
> > timeline.
> > 
> > In other words, if your TOTAL timeline is, say, 50 microseconds, 3
> > microseconds isn't going to buy you that much.  I'm still wondering what you
> > are doing that makes that 3 usec interesting.
> 
> Ten years ago I was using my C++ vector/matrix code for an
> experimental, real-time GPS/INS precision landing and guidance system,
> including a 16-state Kalman filter (if I recall correctly). We had a
> dozen or so tasks, and the guidance update rate was 20 Hz. It was a
> fairly ambitious project at the time, and it worked out well. It might
> NOT have worked had we been generating temporary vectors and matrices
> at a high rate.
> 
> But that's all beside the point. Who are you to question whether a
> certain level of efficiency is enough? Just because *you* don't need
> high efficiency for your work, that doesn't mean that *nobody* needs
> it. And by the way, I'm not talking about ultra-efficiency here. I'm
> just talking about the basics. You could write the slickest algorithm
> ever written, but if you generate temporary vectors and matrices at a
> high rate, it's like running a 100-yard dash carrying a freakin'
> bowling ball.
> 
> > > > And I'd be REALLY interested to know what you are doing with matrix
> > > > multiplication such that the product of a matrix and a scalar is the
> > > > slightest bit interesting.  (Usually, when I'm multiplying matrices, I'm
> > > > playing with direction cosine matrices, and I don't recall ever hearing
> > > > about any use for multiplying a DCM by a scalar.)
> > >
> > > Sorry, but I'm having some trouble with your "logic". DCM matrices
> > > don't get multiplied by scalars, so you doubt that any vector or
> > > matrix ever gets multiplied by a scalar?
> > 
> > Since you want to get nasty personal about it...
> > 
> > I am fully aware that VECTOR * scalar happens, frequently.  I was asking for
> > an application that did MATRIX * scalar.
> 
> Well, first of all, a vector is really just a matrix with one column.
> More importantly, as I wrote earlier (and which you apparently
> missed), all the efficiency concerns regarding "=" and "+=" are
> EXACTLY the same for both vectors *and* matrices. Not one iota of
> difference. So it really doesn't matter whether we are talking about
> vectors or matrices. The fact that I referred to matrices in my
> original post is irrrelevant.
> 
> > > Actually, it is more common to scale a vector than a matrix, but the
> > > same principles regarding the efficiency of "+"and "+=" still apply.
> > > Vectors often need to be scaled for conversion of units, for example.
> > > Another very common example would be determining a change in position
> > > by multiplying a velocity vector by a scalar time increment. Ditto for
> > > determining a change in a velocity vector by multiplying an
> > > acceleration vector by a time increment.
> > 
> > Determining a change in a position vector by multiplying a velocity vector
> > by a scalar time increment, or determining a change in a velocity vector by
> > multiplying an acceleration vector by a time increment, just like that, is
> > Euler's Method for numerical integration.  It is about the worst method
> > known to man, BOTH in terms of processing time AND error buildup (and
> > demonstrating this fact is a standard assignment in the undergraduate
> > numerical methods survey course).  If you are using Euler's Method for a
> > high-performance sim, you are out of your mind.  Runge-Kutta methods offer
> > orders of magnitude better performance, in time AND error control.
> 
> Oh, good lord. You still need to multiply velocity by time, regardless
> of which method of numerical integration you are using. In fact, you
> need to do it more with the higher-order methods.
> 
> > And don't change the subject.  We are talking about MATRIX * scalar, not
> > vector * scalar.  Quaternions are vectors, not matrices.
> 
> No, *you* are talking only about multiplying matrices by scalars, but
> I can't understand why, since the principle in question applies
> equally to both vectors *and* matrices.
> 
> > > Since you know about DCMs, you must also be familiar with quaternions,
> > > otherwise known as Euler parameters. Like DCMs, Euler parameter
> > > "vectors" represents the attitude of a rigid body. Well, an Euler
> > > parameter vector consists of 4 scalar components or "Euler
> > > parameters", and the norm (roor sum square) of the 4 parameters is
> > > unity by definition. However, due to numerical roundoff error, the
> > > norm can drift away from one. However, you can renormalize by dividing
> > > the vector by its own norm. This can be written very concisely and
> > > efficiently in C++ as
> > >
> > >     eulpar /= eulpar.norm()
> > 
> > Concisely, yes.  Efficiently, not necessarily.  You still have to compute
> > the norm, which involves taking a square root, which, off the top of my
> > head, probably costs quite a bit more than four divisions.  If you are
> > really into microefficiency, you would code this as renorm(eulpar), or
> > eulpar.renorm(), and do it in-place, computing 1/sqrt(norm) directly and
> > multplying, rather than computing sqrt() and dividing.
> 
> Taking a square root is ultra-fast on modern processors, but that is
> irrelevant because there is no way around it for normalizing a
> quaternion. And, no, I'm not into micro-efficiency at all. However,
> you may have a point: if multiplicatioin is slightly more efficient
> than division (is it?), then it *would* make sense to do the division
> once and then do multiplication after that. But *that's*
> micro-efficiency.

Re: Is ther any sense in *= and matrices?

[Bobby D. Bryant]

On Sun, 08 Jun 2003 11:51:31 -0700, Russ wrote:

> "John R. Strohm" <strohm@airmail.net> wrote in message
> news:<bbumef$fnq@library2.airnews.net>...

>> With all due respect, ladies and gentlemen, it has been known for a
>> very long time that the difference in "efficiency" between A := A + B
>> and A += B is lost in the noise floor compared to the improvements that
>> can be gotten by improving the algorithms involved.
> 
> Oh, really? I just did a test in C++ with 3x3 matrices. I added them
> together 10,000,000 times using "+", then "+=". The "+=" version took
> about 19 seconds, and the "+" version took about 55 seconds. That's just
> shy of a factor of 3, folks. If that's your "noise floor," I can't help
> wonder what kind of "algorithms" you are dealing with!

I'm just curious why the compiler didn't generate the same code for both
versions.

-- 
Bobby Bryant
Austin, Texas

Re: Is ther any sense in *= and matrices?

[Russ]

"Bobby D. Bryant" <bdbryant@mail.utexas.edu> wrote in message news:<pan.2003.06.08.22.56.37.297066@mail.utexas.edu>...
> On Sun, 08 Jun 2003 11:51:31 -0700, Russ wrote:
> 
> > "John R. Strohm" <strohm@airmail.net> wrote in message
> > news:<bbumef$fnq@library2.airnews.net>...
>  
> >> With all due respect, ladies and gentlemen, it has been known for a
> >> very long time that the difference in "efficiency" between A := A + B
> >> and A += B is lost in the noise floor compared to the improvements that
> >> can be gotten by improving the algorithms involved.
> > 
> > Oh, really? I just did a test in C++ with 3x3 matrices. I added them
> > together 10,000,000 times using "+", then "+=". The "+=" version took
> > about 19 seconds, and the "+" version took about 55 seconds. That's just
> > shy of a factor of 3, folks. If that's your "noise floor," I can't help
> > wonder what kind of "algorithms" you are dealing with!
> 
> I'm just curious why the compiler didn't generate the same code for both
> versions.

The compiler didn't generate the code. I did.

Re: Is ther any sense in *= and matrices?

[John R. Strohm]

X-A-Notice: References line has been trimed due to 512 byte limitation
Abuse-Reports-To: abuse at airmail.net to report improper postings
NNTP-Proxy-Relay: library2.airnews.net
NNTP-Posting-Time: Mon, 09 Jun 2003 00:20:49 -0500 (CDT)
NNTP-Posting-Host: !\p8l1k-V\n`(V3 (Encoded at Airnews!)
X-Priority: 3
X-MSMail-Priority: Normal
X-Newsreader: Microsoft Outlook Express 6.00.2800.1106
X-MimeOLE: Produced By Microsoft MimeOLE V6.00.2800.1106

"Russ" <18k11tm001@sneakemail.com> wrote in message
news:bebbba07.0306082027.2c50eb6b@posting.google.com...
> "Bobby D. Bryant" <bdbryant@mail.utexas.edu> wrote in message
news:<pan.2003.06.08.22.56.37.297066@mail.utexas.edu>...
> > On Sun, 08 Jun 2003 11:51:31 -0700, Russ wrote:
> >
> > > "John R. Strohm" <strohm@airmail.net> wrote in message
> > > news:<bbumef$fnq@library2.airnews.net>...
> >
> > >> With all due respect, ladies and gentlemen, it has been known for a
> > >> very long time that the difference in "efficiency" between A := A + B
> > >> and A += B is lost in the noise floor compared to the improvements
that
> > >> can be gotten by improving the algorithms involved.
> > >
> > > Oh, really? I just did a test in C++ with 3x3 matrices. I added them
> > > together 10,000,000 times using "+", then "+=". The "+=" version took
> > > about 19 seconds, and the "+" version took about 55 seconds. That's
just
> > > shy of a factor of 3, folks. If that's your "noise floor," I can't
help
> > > wonder what kind of "algorithms" you are dealing with!
> >
> > I'm just curious why the compiler didn't generate the same code for both
> > versions.
>
> The compiler didn't generate the code. I did.

Why?

If you are hand-generating assembly code for a PowerPC processor, I would
bet a fair number of doughnuts that your cache utilization is shaky at best,
and THAT is what is generating your factor of 3 difference.

Re: Is ther any sense in *= and matrices?

[Bobby D. Bryant]

On Sun, 08 Jun 2003 21:27:13 -0700, Russ wrote:

> "Bobby D. Bryant" <bdbryant@mail.utexas.edu> wrote in message
> news:<pan.2003.06.08.22.56.37.297066@mail.utexas.edu>...
>> On Sun, 08 Jun 2003 11:51:31 -0700, Russ wrote:
>> 
>> > "John R. Strohm" <strohm@airmail.net> wrote in message
>> > news:<bbumef$fnq@library2.airnews.net>...
>>  
>> >> With all due respect, ladies and gentlemen, it has been known for a
>> >> very long time that the difference in "efficiency" between A := A +
>> >> B and A += B is lost in the noise floor compared to the improvements
>> >> that can be gotten by improving the algorithms involved.
>> > 
>> > Oh, really? I just did a test in C++ with 3x3 matrices. I added them
>> > together 10,000,000 times using "+", then "+=". The "+=" version took
>> > about 19 seconds, and the "+" version took about 55 seconds. That's
>> > just shy of a factor of 3, folks. If that's your "noise floor," I
>> > can't help wonder what kind of "algorithms" you are dealing with!
>> 
>> I'm just curious why the compiler didn't generate the same code for
>> both versions.
> 
> The compiler didn't generate the code. I did.

No, I'm talking about the compiler's machine-code output.  You wrote in
C++, right?

-- 
Bobby Bryant
Austin, Texas

Re: Is ther any sense in *= and matrices?

[Russ]

"Bobby D. Bryant" <bdbryant@mail.utexas.edu> wrote in message news:<pan.2003.06.09.14.53.15.912223@mail.utexas.edu>...
> On Sun, 08 Jun 2003 21:27:13 -0700, Russ wrote:
> 
> > "Bobby D. Bryant" <bdbryant@mail.utexas.edu> wrote in message
> > news:<pan.2003.06.08.22.56.37.297066@mail.utexas.edu>...
> >> On Sun, 08 Jun 2003 11:51:31 -0700, Russ wrote:
> >> 
> >> > "John R. Strohm" <strohm@airmail.net> wrote in message
> >> > news:<bbumef$fnq@library2.airnews.net>...
>  
> >> >> With all due respect, ladies and gentlemen, it has been known for a
> >> >> very long time that the difference in "efficiency" between A := A +
> >> >> B and A += B is lost in the noise floor compared to the improvements
> >> >> that can be gotten by improving the algorithms involved.
> >> > 
> >> > Oh, really? I just did a test in C++ with 3x3 matrices. I added them
> >> > together 10,000,000 times using "+", then "+=". The "+=" version took
> >> > about 19 seconds, and the "+" version took about 55 seconds. That's
> >> > just shy of a factor of 3, folks. If that's your "noise floor," I
> >> > can't help wonder what kind of "algorithms" you are dealing with!
> >> 
> >> I'm just curious why the compiler didn't generate the same code for
> >> both versions.
> > 
> > The compiler didn't generate the code. I did.
> 
> No, I'm talking about the compiler's machine-code output.  You wrote in
> C++, right?

Yes, I wrote it in C++. I would not have expected the compiler to
generate the same code for both versions. I did not use any
optimization flag, by the way, but even if I had, I still would not
have expected the same executable. If you think it should, your beef
is with the gcc folks.

Re: Is ther any sense in *= and matrices?

[Ole-Hjalmar Kristensen]

18k11tm001@sneakemail.com (Russ) writes:

> "Bobby D. Bryant" <bdbryant@mail.utexas.edu> wrote in message news:<pan.2003.06.09.14.53.15.912223@mail.utexas.edu>...
> > On Sun, 08 Jun 2003 21:27:13 -0700, Russ wrote:
> > 
> > > "Bobby D. Bryant" <bdbryant@mail.utexas.edu> wrote in message
> > > news:<pan.2003.06.08.22.56.37.297066@mail.utexas.edu>...
> > >> On Sun, 08 Jun 2003 11:51:31 -0700, Russ wrote:
> > >> 
> > >> > "John R. Strohm" <strohm@airmail.net> wrote in message
> > >> > news:<bbumef$fnq@library2.airnews.net>...
> >  
> > >> >> With all due respect, ladies and gentlemen, it has been known for a
> > >> >> very long time that the difference in "efficiency" between A := A +
> > >> >> B and A += B is lost in the noise floor compared to the improvements
> > >> >> that can be gotten by improving the algorithms involved.
> > >> > 
> > >> > Oh, really? I just did a test in C++ with 3x3 matrices. I added them
> > >> > together 10,000,000 times using "+", then "+=". The "+=" version took
> > >> > about 19 seconds, and the "+" version took about 55 seconds. That's
> > >> > just shy of a factor of 3, folks. If that's your "noise floor," I
> > >> > can't help wonder what kind of "algorithms" you are dealing with!
> > >> 
> > >> I'm just curious why the compiler didn't generate the same code for
> > >> both versions.
> > > 
> > > The compiler didn't generate the code. I did.
> > 
> > No, I'm talking about the compiler's machine-code output.  You wrote in
> > C++, right?
> 
> Yes, I wrote it in C++. I would not have expected the compiler to
> generate the same code for both versions. I did not use any
> optimization flag, by the way, but even if I had, I still would not
> have expected the same executable. If you think it should, your beef
> is with the gcc folks.

You cannot do source level optimizations based on runs with the
optimizer turned off.  For example, you usually need some level of
optimization to turn inlining on, and without inlining, the peephole
optimizer cannot do as much.

When that is said, my run of tmoran's code shows a relatively large
performance penalty for the function version, which supports your argument.
It just shows that you cannot expect the same results on different systems.

oleh@VOLGA /cygdrive/c/div
$ g++ -O3 -funroll-loops -o moran-add moran-add.cpp

oleh@VOLGA /cygdrive/c/div
$ time ./moran-add
Go 1000 180 450
real    0m0.655s
user    0m0.660s
sys     0m0.010s

Ole-Hj. Kristensen

******************************************************************************
* You cannot consistently believe this sentence.
******************************************************************************

Re: Is ther any sense in *= and matrices?

[John R. Strohm]

"Preben Randhol" <randhol+abuse@pvv.org> wrote in message
news:slrnbduaav.6dd.randhol+abuse@kiuk0152.chembio.ntnu.no...
> Vinzent Hoefler wrote:
> > Preben Randhol wrote:
> >
> >>I mean even if one define "*" for Matrix_Type one cannot do:
> >>
> >>   A := A * B
> >
> > You can, iff A and B are square matrices.
>
> Yes, but that is a special case.

Not quite.

A[m x n] = B[m x p] * C[p x n]

You would want *= to work for a general result, which is

A[m x n] = A[m x n] * B[n x n]

If the matrices are not compatible for multiplication, then their product is
not mathematically defined, at all.

Usually, of course, you are only interested in square matrices.

Re: Is ther any sense in *= and matrices?

[Wesley Groleau]

>>I mean even if one define "*" for Matrix_Type one cannot do:
>>
>>  A := A * B
> 
> You can, iff A and B are square matrices.

But you cannot do it efficiently "in-place"
because each element in the result is the sum
of multiply operators on several pairs of
elements in the two addends.  So if I compute

A(1,1) before A (1,2) then the latter has already
had one of its inputs changed.  To overcome this
without introducing a temporary "A" one would
require a rather complicated set of temporary
_parts_ of A.

Re: Is ther any sense in *= and matrices?

[David C. Hoos]

"Wesley Groleau" <wesgroleau@despammed.com> wrote in message
news:qgqdnUqY3I2EBEKjXTWcpg@gbronline.com...
>
> >>I mean even if one define "*" for Matrix_Type one cannot do:
> >>
> >>  A := A * B
> >
> > You can, iff A and B are square matrices.
>
> But you cannot do it efficiently "in-place"
> because each element in the result is the sum
> of multiply operators on several pairs of
> elements in the two addends.  So if I compute

The operands of * are called "factors", not "addends."

>
> A(1,1) before A (1,2) then the latter has already
> had one of its inputs changed.  To overcome this
> without introducing a temporary "A" one would
> require a rather complicated set of temporary
> _parts_ of A.
>
> _______________________________________________
> comp.lang.ada mailing list
> comp.lang.ada@ada.eu.org
> http://ada.eu.org/mailman/listinfo/comp.lang.ada
>

Re: Is ther any sense in *= and matrices?

[Wesley Groleau]

> The operands of * are called "factors", not "addends."

oops, sloppy terminology.  I haven't multiplied
matrices in ten years, but I think that (other
than jargon issues), my observation is true:

>> if A(1,1) is computed before A (1,2) 
(in A := A * B), then the latter has already
>>had one of its inputs changed.  To overcome this
>>without introducing a temporary "A" one would
>>require a rather complicated set of temporary
>>_parts_ of A.