Matrix multiplication

Pages: 12
fea12rs (13)
Hi people. New to this site, new to c++ so excuse me if i lack the normal decorum.

Im writing a program thats supposed to multiply 2x2 matrices together, multiple times, while updating the values of the elements in the matrix each loop. At the moment its just 2 alternating matrices. I then calculate some values at the end using the 'final matrix' elements.

I should output data which fits a curve with a maximum... a somewhat gaussian shape. The graphs produced are completely random. I think the problem may be in the implementation of the operator *= overload .

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
 
 #include "Matrix.h"
#include<complex>
#include<iostream>
#include<cmath>
MyMatrix::MyMatrix(double k, double n, double d)  // initialize values
{
	
	m11 = cos(k*d);
	m12= (1/n)*sin(k*d);
	m21= -n*sin(k*d);
	m22= cos(k*d);
}

void MyMatrix::operator*=(MyMatrix fred)  //overload the *= to carry out the following -premultiply?
{
	            
	MyMatrix newMat = MyMatrix(0,0,0);
    newMat.m11 = m11*fred.m11 + m12*fred.m21;
	newMat.m12 = m11*fred.m12 + m12*fred.m22;
	newMat.m21 = m21*fred.m11 + m22*fred.m21;
	newMat.m22 = m21*fred.m12 + m22*fred.m22;
 
	m11 = newMat.m11;
	m12 = newMat.m12;
	m21 = newMat.m21;
	m22 = newMat.m22;
 
}

double MyMatrix::calcR()
{
	double a = sqrt((m22*m22 + m12*m12)/(m11*m11 + m21*m21));
//	cout << "a  " << a << "asquared  " << a*a << endl;
	return  ((a-1)*(a-1))/((a+1)*(a+1));
	
	//std::cout<<"KKKKKK"<< (m22/c -m11)/(m11 + m22/c) << std::endl;
}

double MyMatrix::calcT()
{

	double R = calcR();
	return (1.0 - R);
}
double MyMatrix::transpose(MyMatrix jon)
{
	using namespace std;
	MyMatrix transp = MyMatrix(0,1.0,0);
    transp.m11 = m11*jon.m11 + m12*jon.m12;
	transp.m12 = m11*jon.m21 + m12*jon.m22;
	transp.m21 = m21*jon.m11 + m22*jon.m12;
	transp.m22 = m21*jon.m21 + m22*jon.m22;
    cout << " main diagonal sum " <<transp.m11 + transp.m22 <<endl <<" secondary diagonal sum " << transp.m12 +transp.m21 << endl << endl;
     
}


and the main 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
 
#include <iostream>
#include<fstream>
using namespace std;
#include "Matrix.h"
MyMatrix Resultant(double);
double Transpose();
int main()
{
	ofstream sela;
    sela.open("kProject.txt");
	for(double k=0.2; k<2; k+=0.01)
	{
		MyMatrix res = Resultant(k);
		cout << "wavenumber  "<<k<< "\t "<< " R=" << res.calcR() << "\t" <<" T=" << res.calcT()<<endl << endl;
		sela << k << "\t " << res.calcR() << "\t " << res.calcT()<<endl;
		
	
	}
	
}

MyMatrix Resultant(double k)
{	//k in meters
	MyMatrix m1 = MyMatrix(k, 1.0, 125.0);   // divide by 100nm
	MyMatrix m2 = MyMatrix(k, 2.0, 62.5);
	
	MyMatrix mT = m1;
     	mT*=m2;
	for(int i=0; i<10; i++)
	{
		mT*=m1;
		mT*=m2;
	}
     //  mT.transpose(mT);
     cout<< mT.m11 << " <m11  "<< mT.m12 << " <m12  "<< mT.m21 << " <m21  "<< mT.m22 << " <m22" <<endl;
   double det = (mT.m11*mT.m22 - mT.m12*mT.m21);
    cout << "det = " << det << endl << endl;
	return mT;
}


as I say i think the probllem is in the *= overload where I keep defining MyMatrix newMat each time I use the *= in the main. Would much appreciate any advice, please talk to me like a novice im new to this.

Also the det calculation gives 1, when I do it manually from the data outputted by the program i get 1.0004 or 0.999934 etc..
Duoas (8826)
lines 20 & 21 are in error.
fea12rs (13)
I cant see how.. i checked the calculation, im worried about the redeclaration of MyMatrix newMat = MyMatrix(0,0,0) every time I used *= ? Surely this will screw it up ..
Duoas (8826)
Check it again. (what column vector are you using?)
fea12rs (13)
I just did again, im sure its the correct expression for the matrix element?. where m21 is the bottom left and m22 is the bottom right element.

Duoas (8826)
fea12rs (13)
OK so your just taking the piss... cheers
Duoas (8826)
Nah, you're just refusing to accept the answer to your question with any grace. Your math is wrong.
fea12rs (13)
Ok other than the math, which im sorry, isnt wrong. Im using M_ij where i is the row and j the column maybe thats the confusion? Can you see where errors may be
James2250 (573)
Normally I wouldn't question Duoas but the math actually looks correct to me -- at least on lines 20/21. Or am I crazy?

Apart from your transpose function not doing anything (should be returning a value..but you aren't using it anyway) I don't see any major C++ errors.

Your *= function looks ok...you could do it without MyMatrix newMat = MyMatrix(0,0,0) if you wanted to but I don't see that hurting anything. Are you sure all of your calculations are correct?
norm b (284)
The *= function looks okay to me too. I did one iteration with the armadillo lib. Does the output agree with your results? 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
 
#include <armadillo>
#include <cmath>

using namespace arma;

int main() {
    double k = 0.2;
    double n = 1.0;
    double d = 125.0;
    mat m1 { cos(k*d), -n*sin(k*d), (1/n)*sin(k*d), cos(k*d) };
    n = 2.0;
    d = 62.5;
    mat m2 { cos(k*d),  -n*sin(k*d),(1/n)*sin(k*d), cos(k*d) };
    m1.reshape(2,2);
    m1.print("\nm1");
    m2.reshape(2, 2);
    m2.print("\nm2");

    mat mT = m1;
    mT *= m2;
    for(int i=0; i<10; ++i){
        mT *= m1;
        mT *= m2;
    }
    mT.print("\nmT");
}


m1
0.9912 -0.1324
0.1324 0.9912

m2
0.9978 -0.0332
0.1326 0.9978

mT
-0.6964 -0.5854
0.9354 -0.6497
Last edited on
fea12rs (13)
Im pretty sure but checking again cant hurt i guess. I used transpose to chek whether the final matrix is correct as itz supposed to symmetric. Its wierd because using the m_ij values outputted by the function, when i calculate det manually its almost 1 i.e 0.9997 ,wen the det function does it its exactly 1, using the same values??? Wtf. Also i get the ocasional m_ij = 234 or 2e+360 .. and ur nt crazy.. unless im crazy and its circular

The graph looks like a sinusoid on meth dno whats happening
Thanks
fea12rs (13)
Thanks norm il chek tomorrow as its 4am my time nw.
fea12rs (13)
couldnt resist.. naa i get different values..no idea wat the problem is. ill try the armadillo library see what the graph looks like.
norm b (284)
Tried your code and I get the same result. Shouldn't make any difference but the following is the header file that I used. The rest of the code was c/p from above.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
 
#ifndef MYMATRIX_H
#define MYMATRIX_H

class MyMatrix
{
public:
    MyMatrix();
    MyMatrix(double k, double n, double d);
    void operator *=(MyMatrix fred);
    double calcR();
    double calcT();
    double transpose(MyMatrix jon);
    double m11, m12, m21, m22;
};

#endif // MYMATRIX_H


Output:
-0.696393 <m11 -0.5854 <m12 0.935398 <m21 -0.649659 <m22
det = 1

wavenumber 0.2 R=0.0204275 T=0.979572
fea12rs (13)
Yeah sorry you're right. It does give those results. It stil has the det error, where when i do it by calculator its 0.99999 but the goddamm machine gives out 1 exactly. thats odd?? Also every nnow and again m12 or m22 is like 500.. what could cause these random errors any ideas?

Thanks
norm b (284)
Thought you were going to bed.

fea12rs wrote:
when i do it by calculator its 0.99999
Probably need more significant digits. Try putting the following before the loop in main():
cout << setprecision(12);
You'll need #include <iomanip> also.
http://www.cplusplus.com/reference/iomanip/setprecision/?kw=setprecision


Also every nnow and again m12 or m22 is like 500.. what could cause these random errors any ideas?
Dunno. Just the way the math works out, I guess. For k=0.37 armadillo yields (which agrees with your output):
m1
-0.6418 0.7668
-0.7668 -0.6418

m2
-0.4232 -0.4530
1.8121 -0.4232

mT
308.669819683394 -9.742399962588
-242.216884689234 7.648217019585



fea12rs (13)
its the m1 or m2 elements which are bllown up sometimes, but for the most part say 95% of the time are in the normal range for sines and cosines.

The setprecision() worked for the matrix elements, but for det it still gives 1. not even 1.00000 just 1 ?
Is there something wrong with the way ive written it to output det?

Im at a loss, is there any way i can upload armadillo library to my crappy bloodshed? Might be the only solution
norm b (284)
See std::fixed: http://www.cplusplus.com/reference/ios/fixed/?kw=fixed

cout << fixed << det << endl;

What OS are you running? My Linux distro has an Armadillo package.
You can download and compile it: http://arma.sourceforge.net/download.html

Edit: there's also the Eigen lib: http://eigen.tuxfamily.org/index.php?title=Main_Page
Last edited on
fea12rs (13)
windows 8. May aswell be using a gameboy. Yeah ill try just get armadillo running. thanks though. still random errors .
Pages: 12

C++

Forum