Re: matrix multiplication question

This isn't an MSVC6 problem -- your example will also crash under g++.  
I think the symmetric multiply routine assumes that there will be an 
element on the diagonal.  We'll try to get a fix -- but you are welcome 
to also look -- routine mult_symm__ in mtl.h (starting around line 
1589).

So, if you add a statement A(1,1) = 0; -- to put an element 
(structurally) in place, the code will run correctly -- at least on 
g++.

Best Wishes,
Andrew Lumsdaine

On Dec 2, 2003, at 3:11 PM, Sven Olsen wrote:

> Hello MTL developers,
>
> I've just started trying to hack together a thin plate spline 
> interpolator using MTL and ITL.
> I've run into a number of problems, some of which seem to stem from 
> the multiplication operator for symmetric matrices.
>
> Here's my current attempt to isolate one of them:
>
> Given A=
> 1 0
> 0 0
>
> and x=
> 1
> 0
>
> attempts to calculate Ax seem to crash, if A is defined to be a 
> symmetric matrix.
>
> Is this a bug or a feature?
>
> -Sven
>
> (Here's the code I'm using.  It's worth noting that I am working with 
> MSVC6, which I know is not really MTL's most natural habitat.)
>
> /*
> 	The goal here is to perform the multiplication:
> 	1 0		x	1	=	1
> 	0 0			0		0
> */
>
>
> //A good line to have if you are using MSVC6
> #ifdef WIN32
> #pragma warning(disable: 4786)
> #endif
>
> #include "mtl/matrix.h"
> #include "mtl/mtl.h"
> #include "mtl/utils.h"
>
> using namespace mtl;
> typedef  double Type;
>
> //If we use Matrix1 as our matrix class, the program crashes
> typedef matrix< Type,
>                 symmetric<upper>,  compressed<>,
>                 row_major >::type Matrix1;
>
> //But if we use Matrix2, everything is fine.
> typedef matrix < Type,
>                   rectangle <>,
>                   compressed<>,
> 	row_major >::type Matrix2;
>
>
>
> int main ()
> {
>   using std::cout;
>   using std::endl;
>
>   //initilize the matrix and vectors
>   const int size=2;
>   Matrix1 A(size,size);
>   A(0,0)=1;
>
>   dense1D<Type> x(A.ncols());
>   x[0]=1;
>   x[1]=0;
>   dense1D<Type> b1(A.ncols());
>
>   //do the multiplication
>   mtl::mult(A, x, b1);
>
>   //output the results
>   for (int ii=0; ii<size; ii++) {
>     for (int j=0; j<size; j++) {
>       cout.width(6);
>       cout << A(ii, j) << " ";
>     }
>     cout << "  x  ";
>     cout.width(16);
>       cout << x[ii] << "  =  ";
>       cout.width(6);
>       cout << b1[ii] << endl;
>
>   }
>
>   return 0;
> }
>
> _______________________________________________
> This list is archived at http://www.osl.iu.edu/MailArchives/mtl-devel/

_______________________________________________
This list is archived at http://www.osl.iu.edu/MailArchives/mtl-devel/