jiri/LA_library
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

vecmat3 transposed multiplications

Browse Source
This commit is contained in:
Jiri Pittner 2023-11-16 16:19:54 +01:00
parent 00b5163e31
commit c5a2865639
3 changed files with 91 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
t.cc
View File

@ -426,11 +426,12 @@ cout <<b;
cout <<d; cout <<d;
} }
if(0) if(1)
{ {
NRMat<double> a; NRMat<double> a;
cin >>a ; cin >>a ;
NRMat<double> abak=a; NRMat<double> abak=a;
NRMat<double> abak2=a;
NRMat<double> u(a.nrows(),a.nrows()),v(a.ncols(),a.ncols()); NRMat<double> u(a.nrows(),a.nrows()),v(a.ncols(),a.ncols());
NRVec<double>s(a.ncols()<a.nrows()?a.ncols():a.nrows()); NRVec<double>s(a.ncols()<a.nrows()?a.ncols():a.nrows());
singular_decomposition(a,&u,s,&v,0); singular_decomposition(a,&u,s,&v,0);
@ -440,6 +441,14 @@ sdiag.diagonalset(s);
cout <<sdiag; cout <<sdiag;
cout <<v; cout <<v;
cout << "Error "<<(u*sdiag*v-abak).norm()<<endl; cout << "Error "<<(u*sdiag*v-abak).norm()<<endl;
NRMat<double> ai=calcinverse(abak2);
cout <<"regular inverse "<<ai;
NRVec<double>ss(s);ss.copyonwrite();
for(int i=0; i<ss.size(); ++i) if(ss[i]>1e-12) ss[i]=1./ss[i]; else ss[i]=0;
sdiag.diagonalset(ss);
NRMat<double> aii = v.transpose() *sdiag*u.transpose();
cout <<"svd inverse "<<aii;
cout << "Error in inverse "<< (aii-ai).norm()<<endl;
} }
if(0) if(0)
@ -2087,7 +2096,7 @@ for(i=1; i<=n; ++i) v[i]=10.*i;
cout <<v.permuted(p); cout <<v.permuted(p);
} }
if(1) if(0)
{ {
int n; int n;
CyclePerm<int> c; CyclePerm<int> c;

75
vecmat3.cc
View File

@ -97,6 +97,18 @@ const Vec3<T> Vec3<T>::operator*(const Mat3<T> &rhs) const
return r; return r;
}; };
template<typename T>
const Vec3<T> Vec3<T>::timesT(const Mat3<T> &rhs) const
{
Vec3<T> r;
r[0] = q[0]*rhs.q[0][0] + q[1]*rhs.q[0][1] + q[2]*rhs.q[0][2];
r[1] = q[0]*rhs.q[1][0] + q[1]*rhs.q[1][1] + q[2]*rhs.q[1][2];
r[2] = q[0]*rhs.q[2][0] + q[1]*rhs.q[2][1] + q[2]*rhs.q[2][2];
return r;
};
template<typename T> template<typename T>
Mat3<T>& Mat3<T>::operator+=(const Mat3 &rhs) Mat3<T>& Mat3<T>::operator+=(const Mat3 &rhs)
@ -129,6 +141,58 @@ const Mat3<T> Mat3<T>::operator*(const Mat3 &rhs) const
return r; return r;
} }
template<typename T>
const Mat3<T> Mat3<T>::Ttimes(const Mat3 &rhs) const
{
Mat3<T> r;
r[0][0]= q[0][0]*rhs.q[0][0] + q[1][0]*rhs.q[1][0] + q[2][0]*rhs.q[2][0];
r[0][1]= q[0][0]*rhs.q[0][1] + q[1][0]*rhs.q[1][1] + q[2][0]*rhs.q[2][1];
r[0][2]= q[0][0]*rhs.q[0][2] + q[1][0]*rhs.q[1][2] + q[2][0]*rhs.q[2][2];
r[1][0]= q[0][1]*rhs.q[0][0] + q[1][1]*rhs.q[1][0] + q[2][1]*rhs.q[2][0];
r[1][1]= q[0][1]*rhs.q[0][1] + q[1][1]*rhs.q[1][1] + q[2][1]*rhs.q[2][1];
r[1][2]= q[0][1]*rhs.q[0][2] + q[1][1]*rhs.q[1][2] + q[2][1]*rhs.q[2][2];
r[2][0]= q[0][2]*rhs.q[0][0] + q[1][2]*rhs.q[1][0] + q[2][2]*rhs.q[2][0];
r[2][1]= q[0][2]*rhs.q[0][1] + q[1][2]*rhs.q[1][1] + q[2][2]*rhs.q[2][1];
r[2][2]= q[0][2]*rhs.q[0][2] + q[1][2]*rhs.q[1][2] + q[2][2]*rhs.q[2][2];
return r;
}
template<typename T>
const Mat3<T> Mat3<T>::timesT(const Mat3 &rhs) const
{
Mat3<T> r;
r[0][0]= q[0][0]*rhs.q[0][0] + q[0][1]*rhs.q[0][1] + q[0][2]*rhs.q[0][2];
r[0][1]= q[0][0]*rhs.q[1][0] + q[0][1]*rhs.q[1][1] + q[0][2]*rhs.q[1][2];
r[0][2]= q[0][0]*rhs.q[2][0] + q[0][1]*rhs.q[2][1] + q[0][2]*rhs.q[2][2];
r[1][0]= q[1][0]*rhs.q[0][0] + q[1][1]*rhs.q[0][1] + q[1][2]*rhs.q[0][2];
r[1][1]= q[1][0]*rhs.q[1][0] + q[1][1]*rhs.q[1][1] + q[1][2]*rhs.q[1][2];
r[1][2]= q[1][0]*rhs.q[2][0] + q[1][1]*rhs.q[2][1] + q[1][2]*rhs.q[2][2];
r[2][0]= q[2][0]*rhs.q[0][0] + q[2][1]*rhs.q[0][1] + q[2][2]*rhs.q[0][2];
r[2][1]= q[2][0]*rhs.q[1][0] + q[2][1]*rhs.q[1][1] + q[2][2]*rhs.q[1][2];
r[2][2]= q[2][0]*rhs.q[2][0] + q[2][1]*rhs.q[2][1] + q[2][2]*rhs.q[2][2];
return r;
}
template<typename T>
const Mat3<T> Mat3<T>::TtimesT(const Mat3 &rhs) const
{
Mat3<T> r;
r[0][0]= q[0][0]*rhs.q[0][0] + q[1][0]*rhs.q[0][1] + q[2][0]*rhs.q[0][2];
r[0][1]= q[0][0]*rhs.q[1][0] + q[1][0]*rhs.q[1][1] + q[2][0]*rhs.q[1][2];
r[0][2]= q[0][0]*rhs.q[2][0] + q[1][0]*rhs.q[2][1] + q[2][0]*rhs.q[2][2];
r[1][0]= q[0][1]*rhs.q[0][0] + q[1][1]*rhs.q[0][1] + q[2][1]*rhs.q[0][2];
r[1][1]= q[0][1]*rhs.q[1][0] + q[1][1]*rhs.q[1][1] + q[2][1]*rhs.q[1][2];
r[1][2]= q[0][1]*rhs.q[2][0] + q[1][1]*rhs.q[2][1] + q[2][1]*rhs.q[2][2];
r[2][0]= q[0][2]*rhs.q[0][0] + q[1][2]*rhs.q[0][1] + q[2][2]*rhs.q[0][2];
r[2][1]= q[0][2]*rhs.q[1][0] + q[1][2]*rhs.q[1][1] + q[2][2]*rhs.q[1][2];
r[2][2]= q[0][2]*rhs.q[2][0] + q[1][2]*rhs.q[2][1] + q[2][2]*rhs.q[2][2];
return r;
}
template<typename T> template<typename T>
T Mat3<T>::determinant() const T Mat3<T>::determinant() const
{ {
@ -165,6 +229,17 @@ const Vec3<T> Mat3<T>::operator*(const Vec3<T> &rhs) const
return r; return r;
} }
template<typename T>
const Vec3<T> Mat3<T>::Ttimes(const Vec3<T> &rhs) const
{
Vec3<T> r;
r[0] = q[0][0]*rhs.q[0] + q[1][0]*rhs.q[1] + q[2][0]*rhs.q[2];
r[1] = q[0][1]*rhs.q[0] + q[1][1]*rhs.q[1] + q[2][1]*rhs.q[2];
r[2] = q[0][2]*rhs.q[0] + q[1][2]*rhs.q[1] + q[2][2]*rhs.q[2];
return r;
}

5
vecmat3.h
View File

@ -86,6 +86,7 @@ public:
Vec3& normalize(void) {*this /= norm(); return *this;}; Vec3& normalize(void) {*this /= norm(); return *this;};
Vec3& fast_normalize(void); Vec3& fast_normalize(void);
const Vec3 operator*(const Mat3<T> &rhs) const; const Vec3 operator*(const Mat3<T> &rhs) const;
const Vec3 timesT(const Mat3<T> &rhs) const; //with transpose
Mat3<T> outer(const Vec3 &rhs) const; //tensor product Mat3<T> outer(const Vec3 &rhs) const; //tensor product
void inertia(Mat3<T> &itensor, const T weight) const; //contribution to inertia tensor void inertia(Mat3<T> &itensor, const T weight) const; //contribution to inertia tensor
@ -151,7 +152,11 @@ public:
const Mat3 operator+(const Mat3 &rhs) const {return Mat3(*this) += rhs;}; const Mat3 operator+(const Mat3 &rhs) const {return Mat3(*this) += rhs;};
const Mat3 operator-(const Mat3 &rhs) const {return Mat3(*this) -= rhs;}; const Mat3 operator-(const Mat3 &rhs) const {return Mat3(*this) -= rhs;};
const Mat3 operator*(const Mat3 &rhs) const; //matrix product const Mat3 operator*(const Mat3 &rhs) const; //matrix product
const Mat3 timesT(const Mat3 &rhs) const; //matrix product with transpose
const Mat3 Ttimes(const Mat3 &rhs) const; //matrix product with transpose
const Mat3 TtimesT(const Mat3 &rhs) const; //matrix product with transpose
const Vec3<T> operator*(const Vec3<T> &rhs) const; //matrix times vector const Vec3<T> operator*(const Vec3<T> &rhs) const; //matrix times vector
const Vec3<T> Ttimes(const Vec3<T> &rhs) const; //matrix times vector with transpose
T trace() const {return q[0][0]+q[1][1]+q[2][2];}; T trace() const {return q[0][0]+q[1][1]+q[2][2];};
T determinant() const; T determinant() const;
void transposeme(); void transposeme();