00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #ifndef __UT_Matrix2_C__
00016 #define __UT_Matrix2_C__
00017
00018 #include <math.h>
00019 #include <stdio.h>
00020 #include <string.h>
00021 #include "UT_Defines.h"
00022 #include "UT_SysClone.h"
00023 #include "UT_Math.h"
00024 #include "UT_IStream.h"
00025 #include "UT_Matrix2.h"
00026 #include "UT_Matrix3.h"
00027 #include "UT_RootFinder.h"
00028
00029
00030
00031 template <typename T>
00032 UT_TMatrix2<T> &
00033 UT_TMatrix2<T>::operator=(const UT_Matrix3 &m)
00034 {
00035 matx[0][0]=m(0,0); matx[0][1]=m(0,1);
00036 matx[1][0]=m(1,0); matx[1][1]=m(1,1);
00037 return *this;
00038 }
00039
00040 template <typename T>
00041 UT_TMatrix2<T> &
00042 UT_TMatrix2<T>::operator*=(const UT_TMatrix2<T> &m)
00043 {
00044 T a, b;
00045 a = matx[0][0]; b = matx[0][1];
00046 matx[0][0] = a*m(0,0) + b*m(1,0); matx[0][1] = a*m(0,1) + b*m(1,1);
00047
00048 a = matx[1][0]; b = matx[1][1];
00049 matx[1][0] = a*m(0,0) + b*m(1,0); matx[1][1] = a*m(0,1) + b*m(1,1);
00050
00051 return *this;
00052 }
00053
00054 template <typename T>
00055 unsigned
00056 UT_TMatrix2<T>::operator==(const UT_TMatrix2<T> &m) const
00057 {
00058 return (&m == this) ? 1 : (
00059 matx[0][0]==m(0,0) && matx[0][1]==m(0,1) &&
00060 matx[1][0]==m(1,0) && matx[1][1]==m(1,1) );
00061 }
00062
00063 template <typename T>
00064 int
00065 UT_TMatrix2<T>::invert(void)
00066 {
00067 T det = determinant();
00068
00069 if (!UTequalZero(det, tolerance()))
00070 {
00071 T tmp = matx[1][1];
00072 matx[1][1] = matx[0][0]; matx[0][0] = tmp;
00073 matx[0][1] = -matx[0][1]; matx[1][0] = -matx[1][0];
00074
00075 (*this) *= (1.0f/det);
00076
00077 return 0;
00078 }
00079 else
00080 {
00081 return 1;
00082 }
00083 }
00084
00085 template <typename T>
00086 int
00087 UT_TMatrix2<T>::eigenvalues(UT_Vector2 &r, UT_Vector2 &i) const
00088 {
00089 int nroots;
00090
00091 nroots = UT_RootFinder::quadratic(1, -trace(), determinant(), r, i);
00092
00093 return nroots;
00094 }
00095
00096 template <typename T>
00097 int
00098 UT_TMatrix2<T>::invert(UT_TMatrix2<T> &m) const
00099 {
00100 T det = determinant();
00101
00102 if (!UTequalZero(det, tolerance()))
00103 {
00104 m = UT_TMatrix2<T>( matx[1][1],-matx[0][1],
00105 -matx[1][0], matx[0][0]);
00106 m *= (1.0f/det);
00107
00108 return 0;
00109 }
00110 else
00111 {
00112 return 1;
00113 }
00114 }
00115
00116 template <typename T>
00117 int
00118 UT_TMatrix2<T>::solve(const UT_Vector2 &b, UT_Vector2 &x) const
00119 {
00120 UT_TMatrix2<T> inv;
00121 if (invert(inv))
00122 {
00123 return 1;
00124 }
00125 else
00126 {
00127 x = colVecMult(inv,b);
00128 return 0;
00129 }
00130 }
00131
00132 template <typename T>
00133 UT_TMatrix2<T>
00134 UT_TMatrix2<T>::transpose() const
00135 {
00136 T r[2][2];
00137 r[0][1] = matx[1][0]; r[1][0] = matx[0][1];
00138 r[0][0] = matx[0][0]; r[1][1] = matx[1][1];
00139 return UT_TMatrix2<T>(r);
00140 }
00141
00142 template <typename T>
00143 T
00144 UT_TMatrix2<T>::getEuclideanNorm2() const
00145 {
00146 return
00147 matx[0][0]*matx[0][0] + matx[0][1]*matx[0][1] +
00148 matx[1][0]*matx[1][0] + matx[1][1]*matx[1][1];
00149 }
00150
00151 template <typename T>
00152 const char *
00153 UT_TMatrix2<T>::className(void)
00154 {
00155 return "UT_Matrix2";
00156 }
00157
00158 template <typename T>
00159 int
00160 UT_TMatrix2<T>::save(ostream &os, int binary) const
00161 {
00162 if ( binary )
00163 {
00164 os << className();
00165 if (!UTwrite(os, &matx[0][0], 4))
00166 return -1;
00167 }
00168 else
00169 {
00170 os << *this;
00171 os << endl;
00172 if ( !os.good() )
00173 return -1;
00174 }
00175 return 0;
00176 }
00177
00178 template <typename T>
00179 bool
00180 UT_TMatrix2<T>::load(UT_IStream &is)
00181 {
00182 if ( is.isBinary() )
00183 {
00184 char c[UT_SMALLBUF];
00185 size_t len = strlen(className());
00186
00187 if (is.bread(c, len) != len)
00188 return false;
00189 if (strncmp(c, className(), len))
00190 return false;
00191 }
00192 else
00193 {
00194 if (!is.checkToken(className()))
00195 return false;
00196 }
00197 if (is.read(&matx[0][0], 4) != 4)
00198 return false;
00199 return true;
00200 }
00201
00202 template <typename T>
00203 void
00204 UT_TMatrix2<T>::outAsciiNoName(ostream &os) const
00205 {
00206 int old_precision;
00207
00208
00209 old_precision = os.precision(SYS_DBL_DIG);
00210 for (unsigned i = 0; i < 2; i++)
00211 {
00212 for (unsigned j = 0; j < 2; j++)
00213 {
00214 if( i != 0 || j != 0 )
00215 os << ' ';
00216 os << matx[i][j];
00217 }
00218 }
00219 os.precision(old_precision);
00220 }
00221
00222
00223
00224
00225 template <typename T>
00226 UT_TMatrix2<T> operator+(const UT_TMatrix2<T> &m1, const UT_TMatrix2<T> &m2)
00227 {
00228 T m[2][2];
00229 m[0][0] = m1(0,0) + m2(0,0); m[0][1] = m1(0,1) + m2(0,1);
00230 m[1][0] = m1(1,0) + m2(1,0); m[1][1] = m1(1,1) + m2(1,1);
00231 return UT_TMatrix2<T>(m);
00232 }
00233 template <typename T>
00234 UT_TMatrix2<T> operator-(const UT_TMatrix2<T> &m1, const UT_TMatrix2<T> &m2)
00235 {
00236 T m[2][2];
00237 m[0][0] = m1(0,0) - m2(0,0); m[0][1] = m1(0,1) - m2(0,1);
00238 m[1][0] = m1(1,0) - m2(1,0); m[1][1] = m1(1,1) - m2(1,1);
00239 return UT_TMatrix2<T>(m);
00240 }
00241 template <typename T>
00242 UT_TMatrix2<T> operator*(const UT_TMatrix2<T> &m1, const UT_TMatrix2<T> &m2)
00243 {
00244 T m[2][2];
00245 m[0][0] = m1(0,0)*m2(0,0) + m1(0,1)*m2(1,0);
00246 m[0][1] = m1(0,0)*m2(0,1) + m1(0,1)*m2(1,1);
00247
00248 m[1][0] = m1(1,0)*m2(0,0) + m1(1,1)*m2(1,0);
00249 m[1][1] = m1(1,0)*m2(0,1) + m1(1,1)*m2(1,1);
00250 return UT_TMatrix2<T>(m);
00251 }
00252 template <typename T>
00253 UT_TMatrix2<T> operator+(const UT_TMatrix2<T> &m1, const UT_Vector2 &vec)
00254 {
00255 T m[2][2];
00256 m[0][0] = m1(0,0) + vec.x(); m[0][1] = m1(0,1) + vec.x();
00257 m[1][0] = m1(1,0) + vec.y(); m[1][1] = m1(1,1) + vec.y();
00258 return UT_TMatrix2<T>(m);
00259 }
00260 template <typename T>
00261 UT_TMatrix2<T> operator-(const UT_TMatrix2<T> &m1, const UT_Vector2 &vec)
00262 {
00263 T m[2][2];
00264 m[0][0] = m1(0,0) - vec.x(); m[0][1] = m1(0,1) - vec.x();
00265 m[1][0] = m1(1,0) - vec.y(); m[1][1] = m1(1,1) - vec.y();
00266 return UT_TMatrix2<T>(m);
00267 }
00268 template <typename T>
00269 UT_TMatrix2<T> operator-(const UT_Vector2 &vec, const UT_TMatrix2<T> &m1)
00270 {
00271 T m[2][2];
00272 m[0][0] = vec.x() - m1(0,0); m[0][1] = vec.x() - m1(0,1);
00273 m[1][0] = vec.y() - m1(1,0); m[1][1] = vec.y() - m1(1,1);
00274 return UT_TMatrix2<T>(m);
00275 }
00276 template <typename T>
00277 UT_TMatrix2<T> operator+(const UT_TMatrix2<T> &m1, T sc)
00278 {
00279 T m[2][2];
00280 m[0][0] = m1(0,0) + sc; m[0][1] = m1(0,1) + sc;
00281 m[1][0] = m1(1,0) + sc; m[1][1] = m1(1,1) + sc;
00282 return UT_TMatrix2<T>(m);
00283 }
00284 template <typename T>
00285 UT_TMatrix2<T> operator-(T sc, const UT_TMatrix2<T> &m1)
00286 {
00287 T m[2][2];
00288 m[0][0] = sc - m1(0,0); m[0][1] = sc - m1(0,1);
00289 m[1][0] = sc - m1(1,0); m[1][1] = sc - m1(1,1);
00290 return UT_TMatrix2<T>(m);
00291 }
00292 template <typename T>
00293 UT_TMatrix2<T> operator*(const UT_TMatrix2<T> &m1, T sc)
00294 {
00295 T m[2][2];
00296 m[0][0] = m1(0,0) * sc; m[0][1] = m1(0,1) * sc;
00297 m[1][0] = m1(1,0) * sc; m[1][1] = m1(1,1) * sc;
00298 return UT_TMatrix2<T>(m);
00299 }
00300 template <typename T>
00301 UT_TMatrix2<T> operator/(T sc, const UT_TMatrix2<T> &m1)
00302 {
00303 T m[2][2];
00304 m[0][0] = sc / m1(0,0); m[0][1] = sc / m1(0,1);
00305 m[1][0] = sc / m1(1,0); m[1][1] = sc / m1(1,1);
00306 return UT_TMatrix2<T>(m);
00307 }
00308
00309 #endif
