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(Vector<N> const* vectors)
00045 {
00046 for (std::size_t i=0; i < M; ++i)
00047 for (std::size_t j=0; j < N; ++j)
00048 (*this)[i][j] = vectors[i][j];
00049 }
00051 Matrix<M,N> operator+(Matrix<M,N> const& b) const
00052 {
00053 Matrix<M,N> r;
00054 for (std::size_t i=0; i < M; ++i)
00055 r[i] = (*this)[i] + b[i];
00056 return r;
00057 }
00059 Matrix<M,N> operator-(Matrix<M,N> const& b) const
00060 {
00061 Matrix<M,N> r;
00062 for (std::size_t i=0; i < M; ++i)
00063 r[i] = (*this)[i] - b[i];
00064 return r;
00065 }
00067 Matrix<M,N> operator*(double s) const
00068 {
00069 Matrix<M,N> r;
00070 for (std::size_t i=0; i < M; ++i)
00071 r[i] = (*this)[i] * s;
00072 return r;
00073 }
00075 Matrix<M,N> operator/(double s) const
00076 {
00077 Matrix<M,N> r;
00078 for (std::size_t i=0; i < M; ++i)
00079 r[i] = (*this)[i] / s;
00080 return r;
00081 }
00083 Vector<M> operator*(Vector<N> const& b) const
00084 {
00085 Vector<M> r;
00086 for (std::size_t i=0; i < M; ++i)
00087 r[i] = (*this)[i] * b;
00088 return r;
00089 }
00094 template <std::size_t O>
00095 Matrix<M,O> operator*(Matrix<N,O> const& b) const
00096 {
00097 Matrix<M,O> r;
00098 for (std::size_t i=0; i < M; ++i)
00099 for (std::size_t j=0; j < O; ++j)
00100 {
00101 r[i][j] = (*this)[i][0]*b[0][j];
00102 for (std::size_t k=1; k < N; ++k)
00103 r[i][j] += (*this)[i][k]*b[k][j];
00104 }
00105 return r;
00106 }
00107 };
00108
00110 template <std::size_t M, std::size_t N>
00111 Matrix<N,M> transpose(Matrix<M,N> const& m)
00112 {
00113 Matrix<N,M> r;
00114 for (std::size_t i=0; i < M; ++i)
00115 for (std::size_t j=0; j < N; ++j)
00116 r[j][i] = m[i][j];
00117 return r;
00118 }
00119
00121 template <std::size_t M, std::size_t N>
00122 Matrix<M,N> tensorProduct(Vector<M> const& a, Vector<N> const& b)
00123 {
00124 Matrix<M,N> r;
00125 for (std::size_t i=0; i < M; ++i)
00126 r[i] = b * a[i];
00127 return r;
00128 }
00129
00132 template <std::size_t M, std::size_t N>
00133 Matrix<M-1,N-1> getMinor(Matrix<M,N> const& A, std::size_t i, std::size_t j);
00134
00137 template <std::size_t M, std::size_t N>
00138 double getCofactor(Matrix<M,N> const& A, std::size_t i, std::size_t j);
00139
00142 template <std::size_t M, std::size_t N>
00143 double getDeterminant(Matrix<M,N> const& A);
00144
00146 inline Matrix<3,3> cofactor(Matrix<3,3> const &m)
00147 {
00148 Matrix<3,3> r;
00149 for (std::size_t i = 0; i < 3; ++i)
00150 for (std::size_t j = 0; j < 3; ++j)
00151 r[i][j] = getCofactor(m, i, j);
00152 return r;
00153 }
00154
00156 inline Matrix<2,2> invert(Matrix<2,2> const& m)
00157 {
00158 Matrix<2,2> a;
00159 a[0][0] = m[1][1]; a[0][1] = -m[0][1];
00160 a[1][0] = -m[1][0]; a[1][1] = m[0][0];
00161 return a / getDeterminant(m);
00162 }
00163
00165 inline Matrix<3,3> invert(Matrix<3,3> const& m)
00166 {
00167 Matrix<3,3> x = transpose(m);
00168 Matrix<3,3> r;
00169 r[0] = cross(x[1],x[2]);
00170 r[1] = cross(x[2],x[0]);
00171 r[2] = cross(x[0],x[1]);
00172 return r / getDeterminant(m);
00173 }
00174
00178 class Matrix3x3 : public Matrix<3,3>
00179 {
00180 public:
00182 Matrix3x3() {}
00185 Matrix3x3(double a11, double a12, double a13,
00186 double a21, double a22, double a23,
00187 double a31, double a32, double a33)
00188 {
00189 (*this)[0] = Vector3(a11,a12,a13);
00190 (*this)[1] = Vector3(a21,a22,a23);
00191 (*this)[2] = Vector3(a31,a32,a33);
00192 }
00194 Matrix3x3(Matrix<3,3> const& other):
00195 Matrix<3,3>(other)
00196 {}
00199 void toArray(double (*array)[3]) const
00200 {
00201 for (std::size_t i=0; i < 3; ++i)
00202 for (std::size_t j=0; j < 3; ++j)
00203 array[i][j] = (*this)[i][j];
00204 }
00205 };
00206
00208 template <std::size_t M, std::size_t N>
00209 double getInnerProduct(Matrix<M,N> const& a,
00210 Matrix<M,N> const& b)
00211 {
00212 double r = a[0]*b[0];
00213 for (std::size_t i=1; i < M; ++i)
00214 r += a[i]*b[i];
00215 return r;
00216 }
00217
00219 Matrix3x3 cross(Vector3 const& u);
00223 Matrix3x3 rotate(Vector3 const& u, double a);
00233 Matrix3x3 getFrame(Vector3 const& v);
00234
00236 int eigen(Matrix3x3 const& A,
00237 Vector<3>* eigenVectors,
00238 double* eigenValues);
00239
00240 }
00241
00242 template <std::size_t M, std::size_t N>
00243 std::ostream& operator<<(std::ostream& s, apf::Matrix<M,N> const& A)
00244 {
00245 for (std::size_t i = 0; i < M; ++i) {
00246 for (std::size_t j = 0; j < N; ++j)
00247 s << A[i][j] << ' ';
00248 s << '\n';
00249 }
00250 return s;
00251 }
00252
00253 #endif