|
|
Blitz Support : |
From: Julian Cummings (cummings_at_[hidden])
Date: 2005-05-06 13:49:33
Hello Marc,
> -----Original Message-----
> From: blitz-support-bounces_at_[hidden]
> [mailto:blitz-support-bounces_at_[hidden]] On Behalf Of Marc Vinyes
> Sent: Wednesday, May 04, 2005 8:08 AM
> To: blitz-support_at_[hidden]
> Subject: [Blitz-support] matrix product (matrix dimensions
> posted againwithin the list)
>
> Apologies for cross-posting. I already sent a message to the
> list without subscribing but now I have read more things and
> I have some other questions.
>
> http://www.oonumerics.org/MailArchives/blitz-support/2004/01/0993.php
>
> I got the same problem:
>
> Array<float,2> B(2,3),C(3,1);
> B=1,2,3,4,5,6;
> cout<<B;
> C=1,2,3;
> cout<<C;
> Array<float,2> A(B*C);
> cout<<A;
>
> //Has as a result:
> // 2 x 3
> // [ 1 2 3
> // 4 5 6 ]
> // 3 x 1
> // [ 1
> // 2
> // 3 ]
> // 1 x 1
> // [ 1 ]
This one has the wrong output Array shape because the Array A was never
allocated.
>
> And,also, for example:
>
> Array<int,2> B(2,3),C(3,1);
> B=1,2,3,4,5,6;
> cout<<B;
> C=1,2,3;
> cout<<C;
> Array<int,2> A(2,1);
> A=B*C;
> cout<<A;
> //obtaining:
> //2 x 3
> //[ 1 2 3
> // 4 5 6 ]
> //3 x 1
> //[ 1
> // 2
> // 3 ]
> //2 x 1
> //[ 1
> // 4 ]
>
> when "the result sould be intuitively":
> 2x1
> [ 14
> 32]
Blitz Arrays are intended to perform elementwise math operations and that is
what the operator*() does. You get a 2x1 result Array because that is what
you allocated. The result Array A is filled by marching through Arrays B
and C in elementwise order and performing multiplication. These are
row-major Arrays, so we get 1x1=1 and 2x2=4 for the results that are stored
into A. Of course, this is not the semantics you are looking for. What you
want is to loop over the rows of B and the columns of C and perform inner
products to get the resulting element values of A. You can do this with the
proper use of Array slicing and summation.
void matmatProduct(Array<int,2> A, Array<int,2> B, Array<int,2> C) {
int rows = B.rows();
int columns = C.columns();
assert(B.columns() == C.rows()); // sanity check
assert(A.rows()==rows && A.columns()==columns); // could reshape A here
instead
for (int i=0; i<rows; ++i)
for (int j=0; j<columns; ++j)
A(i,j) = sum(B(i,Range::all())*C(Range::all(),j));
}
Sorry there is no simple way to do this by overloading the standard *
notation, but this function will give you the result you expected.
Regards, Julian C.
Dr. Julian C. Cummings Office: PB-111
Caltech/CACR, MC 158-79 Phone: 626-395-2543
1200 E. California Blvd. Fax: 626-584-5917
Pasadena, CA 91125