00001
00002
00003
00004
00005
00006
00007
00008 #ifndef APFMATRIX_H
00009 #define APFMATRIX_H
00010
00014 #include "apfVector.h"
00015
00016 namespace apf {
00017
00030 template <std::size_t M, std::size_t N>
00031 class Matrix : public Array<Vector<N>,M>
00032 {
00033 public:
00035 Matrix() {}
00037 Matrix(double const (*array)[N])
00038 {
00039 for (std::size_t i=0; i < M; ++i)
00040 for (std::size_t j=0; j < N; ++j)
00041 (*this)[i][j] = array[i][j];
00042 }
00044 Matrix<M,N> operator+(Matrix<M,N> const& b) const
00045 {
00046 Matrix<M,N> r;
00047 for (std::size_t i=0; i < M; ++i)
00048 r[i] = (*this)[i] + b[i];
00049 return r;
00050 }
00052 Matrix<M,N> operator-(Matrix<M,N> const& b) const
00053 {
00054 Matrix<M,N> r;
00055 for (std::size_t i=0; i < M; ++i)
00056 r[i] = (*this)[i] - b[i];
00057 return r;
00058 }
00060 Matrix<M,N> operator*(double s) const
00061 {
00062 Matrix<M,N> r;
00063 for (std::size_t i=0; i < M; ++i)
00064 r[i] = (*this)[i] * s;
00065 return r;
00066 }
00068 Matrix<M,N> operator/(double s) const
00069 {
00070 Matrix<M,N> r;
00071 for (std::size_t i=0; i < M; ++i)
00072 r[i] = (*this)[i] / s;
00073 return r;
00074 }
00076 Vector<M> operator*(Vector<N> const& b) const
00077 {
00078 Vector<M> r;
00079 for (std::size_t i=0; i < M; ++i)
00080 r[i] = (*this)[i] * b;
00081 return r;
00082 }
00087 template <std::size_t O>
00088 Matrix<M,O> operator*(Matrix<N,O> const& b) const
00089 {
00090 Matrix<M,O> r;
00091 for (std::size_t i=0; i < M; ++i)
00092 for (std::size_t j=0; j < O; ++j)
00093 {
00094 r[i][j] = (*this)[i][0]*b[0][j];
00095 for (std::size_t k=1; k < N; ++k)
00096 r[i][j] += (*this)[i][k]*b[k][j];
00097 }
00098 return r;
00099 }
00100 };
00101
00103 template <std::size_t M, std::size_t N>
00104 Matrix<N,M> transpose(Matrix<M,N> const& m)
00105 {
00106 Matrix<N,M> r;
00107 for (std::size_t i=0; i < M; ++i)
00108 for (std::size_t j=0; j < N; ++j)
00109 r[j][i] = m[i][j];
00110 return r;
00111 }
00112
00114 template <std::size_t M, std::size_t N>
00115 Matrix<M,N> tensorProduct(Vector<M> const& a, Vector<N> const& b)
00116 {
00117 Matrix<M,N> r;
00118 for (std::size_t i=0; i < M; ++i)
00119 r[i] = b * a[i];
00120 return r;
00121 }
00122
00125 template <std::size_t M, std::size_t N>
00126 Matrix<M-1,N-1> getMinor(Matrix<M,N> const& A, std::size_t i, std::size_t j);
00127
00130 template <std::size_t M, std::size_t N>
00131 double getCofactor(Matrix<M,N> const& A, std::size_t i, std::size_t j);
00132
00135 template <std::size_t M, std::size_t N>
00136 double getDeterminant(Matrix<M,N> const& A);
00137
00139 inline Matrix<3,3> cofactor(Matrix<3,3> const &m)
00140 {
00141 Matrix<3,3> r;
00142 for (std::size_t i = 0; i < 3; ++i)
00143 for (std::size_t j = 0; j < 3; ++j)
00144 r[i][j] = getCofactor(m, i, j);
00145 return r;
00146 }
00147
00149 inline Matrix<2,2> invert(Matrix<2,2> const& m)
00150 {
00151 Matrix<2,2> a;
00152 a[0][0] = m[1][1]; a[0][1] = -m[0][1];
00153 a[1][0] = -m[1][0]; a[1][1] = m[0][0];
00154 return a / getDeterminant(m);
00155 }
00156
00158 inline Matrix<3,3> invert(Matrix<3,3> const& m)
00159 {
00160 Matrix<3,3> x = transpose(m);
00161 Matrix<3,3> r;
00162 r[0] = cross(x[1],x[2]);
00163 r[1] = cross(x[2],x[0]);
00164 r[2] = cross(x[0],x[1]);
00165 return r / getDeterminant(m);
00166 }
00167
00171 class Matrix3x3 : public Matrix<3,3>
00172 {
00173 public:
00175 Matrix3x3() {}
00178 Matrix3x3(double a11, double a12, double a13,
00179 double a21, double a22, double a23,
00180 double a31, double a32, double a33)
00181 {
00182 (*this)[0] = Vector3(a11,a12,a13);
00183 (*this)[1] = Vector3(a21,a22,a23);
00184 (*this)[2] = Vector3(a31,a32,a33);
00185 }
00187 Matrix3x3(Matrix<3,3> const& other):
00188 Matrix<3,3>(other)
00189 {}
00192 void toArray(double (*array)[3]) const
00193 {
00194 for (std::size_t i=0; i < 3; ++i)
00195 for (std::size_t j=0; j < 3; ++j)
00196 array[i][j] = (*this)[i][j];
00197 }
00198 };
00199
00201 template <std::size_t M, std::size_t N>
00202 double getInnerProduct(Matrix<M,N> const& a,
00203 Matrix<M,N> const& b)
00204 {
00205 double r = a[0]*b[0];
00206 for (std::size_t i=1; i < M; ++i)
00207 r += a[i]*b[i];
00208 return r;
00209 }
00210
00212 Matrix3x3 cross(Vector3 const& u);
00216 Matrix3x3 rotate(Vector3 const& u, double a);
00226 Matrix3x3 getFrame(Vector3 const& v);
00227
00229 int eigen(Matrix3x3 const& A,
00230 Vector<3>* eigenVectors,
00231 double* eigenValues);
00232
00233 }
00234
00235 template <std::size_t M, std::size_t N>
00236 std::ostream& operator<<(std::ostream& s, apf::Matrix<M,N> const& A)
00237 {
00238 for (std::size_t i = 0; i < M; ++i) {
00239 for (std::size_t j = 0; j < N; ++j)
00240 s << A[i][j] << ' ';
00241 s << '\n';
00242 }
00243 return s;
00244 }
00245
00246 #endif