Functions.h
Go to the documentation of this file.
00001 #ifndef __NEW_MAT_FUNCTIONS_H__
00002 #define __NEW_MAT_FUNCTIONS_H__
00003 
00004 #ifdef use_namespace
00005 namespace NEWMAT {
00006 #endif
00007 
00008 
00009 // ************************ functions ************************************ //
00010 
00011 bool operator==(const GeneralMatrix& A, const GeneralMatrix& B);
00012 bool operator==(const BaseMatrix& A, const BaseMatrix& B);
00013 inline bool operator!=(const GeneralMatrix& A, const GeneralMatrix& B)
00014    { return ! (A==B); }
00015 inline bool operator!=(const BaseMatrix& A, const BaseMatrix& B)
00016    { return ! (A==B); }
00017 
00018    // inequality operators are dummies included for compatibility
00019    // with STL. They throw an exception if actually called.
00020 inline bool operator<=(const BaseMatrix& A, const BaseMatrix&)
00021    { A.IEQND(); return true; }
00022 inline bool operator>=(const BaseMatrix& A, const BaseMatrix&)
00023    { A.IEQND(); return true; }
00024 inline bool operator<(const BaseMatrix& A, const BaseMatrix&)
00025    { A.IEQND(); return true; }
00026 inline bool operator>(const BaseMatrix& A, const BaseMatrix&)
00027    { A.IEQND(); return true; }
00028 
00029 bool is_zero(const BaseMatrix& A);
00030 inline bool IsZero(const BaseMatrix& A) { return is_zero(A); }
00031 
00032 Real dotproduct(const Matrix& A, const Matrix& B);
00033 Matrix crossproduct(const Matrix& A, const Matrix& B);
00034 ReturnMatrix crossproduct_rows(const Matrix& A, const Matrix& B);
00035 ReturnMatrix crossproduct_columns(const Matrix& A, const Matrix& B);
00036 
00037 inline Real DotProduct(const Matrix& A, const Matrix& B)
00038    { return dotproduct(A, B); }
00039 inline Matrix CrossProduct(const Matrix& A, const Matrix& B)
00040    { return crossproduct(A, B); }
00041 inline ReturnMatrix CrossProductRows(const Matrix& A, const Matrix& B)
00042    { return crossproduct_rows(A, B); }
00043 inline ReturnMatrix CrossProductColumns(const Matrix& A, const Matrix& B)
00044    { return crossproduct_columns(A, B); }
00045 
00046 void newmat_block_copy(int n, Real* from, Real* to);
00047 
00048 // ********************* friend functions ******************************** //
00049 
00050 // Functions declared as friends - G++ wants them declared externally as well
00051 
00052 bool Rectangular(MatrixType a, MatrixType b, MatrixType c);
00053 bool Compare(const MatrixType&, MatrixType&);
00054 Real dotproduct(const Matrix& A, const Matrix& B);
00055 SPMatrix SP(const BaseMatrix&, const BaseMatrix&);
00056 KPMatrix KP(const BaseMatrix&, const BaseMatrix&);
00057 ShiftedMatrix operator+(Real f, const BaseMatrix& BM);
00058 NegShiftedMatrix operator-(Real, const BaseMatrix&);
00059 ScaledMatrix operator*(Real f, const BaseMatrix& BM);
00060 
00061 // ********************* inline functions ******************************** //
00062 
00063 inline LogAndSign log_determinant(const BaseMatrix& B)
00064    { return B.log_determinant(); }
00065 inline LogAndSign LogDeterminant(const BaseMatrix& B)
00066    { return B.log_determinant(); }
00067 inline Real determinant(const BaseMatrix& B)
00068    { return B.determinant(); }
00069 inline Real Determinant(const BaseMatrix& B)
00070    { return B.determinant(); }
00071 inline Real sum_square(const BaseMatrix& B) { return B.sum_square(); }
00072 inline Real SumSquare(const BaseMatrix& B) { return B.sum_square(); }
00073 inline Real norm_Frobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
00074 inline Real norm_frobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
00075 inline Real NormFrobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
00076 inline Real trace(const BaseMatrix& B) { return B.trace(); }
00077 inline Real Trace(const BaseMatrix& B) { return B.trace(); }
00078 inline Real sum_absolute_value(const BaseMatrix& B)
00079    { return B.sum_absolute_value(); }
00080 inline Real SumAbsoluteValue(const BaseMatrix& B)
00081    { return B.sum_absolute_value(); }
00082 inline Real sum(const BaseMatrix& B)
00083    { return B.sum(); }
00084 inline Real Sum(const BaseMatrix& B)
00085    { return B.sum(); }
00086 inline Real maximum_absolute_value(const BaseMatrix& B)
00087    { return B.maximum_absolute_value(); }
00088 inline Real MaximumAbsoluteValue(const BaseMatrix& B)
00089    { return B.maximum_absolute_value(); }
00090 inline Real minimum_absolute_value(const BaseMatrix& B)
00091    { return B.minimum_absolute_value(); }
00092 inline Real MinimumAbsoluteValue(const BaseMatrix& B)
00093    { return B.minimum_absolute_value(); }
00094 inline Real maximum(const BaseMatrix& B) { return B.maximum(); }
00095 inline Real Maximum(const BaseMatrix& B) { return B.maximum(); }
00096 inline Real minimum(const BaseMatrix& B) { return B.minimum(); }
00097 inline Real Minimum(const BaseMatrix& B) { return B.minimum(); }
00098 inline Real norm1(const BaseMatrix& B) { return B.norm1(); }
00099 inline Real Norm1(const BaseMatrix& B) { return B.norm1(); }
00100 inline Real norm1(RowVector& RV) { return RV.maximum_absolute_value(); }
00101 inline Real Norm1(RowVector& RV) { return RV.maximum_absolute_value(); }
00102 inline Real norm_infinity(const BaseMatrix& B) { return B.norm_infinity(); }
00103 inline Real NormInfinity(const BaseMatrix& B) { return B.norm_infinity(); }
00104 inline Real norm_infinity(ColumnVector& CV)
00105    { return CV.maximum_absolute_value(); }
00106 inline Real NormInfinity(ColumnVector& CV)
00107    { return CV.maximum_absolute_value(); }
00108 inline bool IsZero(const GeneralMatrix& A) { return A.IsZero(); }
00109 inline bool is_zero(const GeneralMatrix& A) { return A.is_zero(); }
00110 
00111 
00112 inline MatrixInput MatrixInput::operator<<(int f) { return *this << (Real)f; }
00113 inline MatrixInput GeneralMatrix::operator<<(int f) { return *this << (Real)f; }
00114 inline MatrixInput BandMatrix::operator<<(int f) { return *this << (Real)f; }
00115 inline MatrixInput GetSubMatrix::operator<<(int f) { return *this << (Real)f; }
00116 
00117 inline ReversedMatrix BaseMatrix::Reverse() const { return reverse(); }
00118 inline RowedMatrix BaseMatrix::AsRow() const { return as_row(); }
00119 inline ColedMatrix BaseMatrix::AsColumn() const { return as_column(); }
00120 inline DiagedMatrix BaseMatrix::AsDiagonal() const { return as_diagonal(); }
00121 inline MatedMatrix BaseMatrix::AsMatrix(int m, int n) const
00122    { return as_matrix(m, n); }
00123 inline GetSubMatrix BaseMatrix::SubMatrix(int fr, int lr, int fc, int lc) const
00124    { return submatrix(fr, lr, fc, lc); }
00125 inline GetSubMatrix BaseMatrix::SymSubMatrix(int f, int l) const
00126    { return sym_submatrix(f, l); }
00127 inline GetSubMatrix BaseMatrix::Row(int f) const { return row(f); }
00128 inline GetSubMatrix BaseMatrix::Rows(int f, int l) const { return rows(f, l); }
00129 inline GetSubMatrix BaseMatrix::Column(int f) const { return column(f); }
00130 inline GetSubMatrix BaseMatrix::Columns(int f, int l) const
00131    { return columns(f, l); }
00132 inline Real BaseMatrix::AsScalar() const { return as_scalar(); }
00133 
00134 inline ReturnMatrix GeneralMatrix::ForReturn() const { return for_return(); }
00135 
00136 inline void swap(Matrix& A, Matrix& B) { A.swap(B); }
00137 inline void swap(SquareMatrix& A, SquareMatrix& B) { A.swap(B); }
00138 inline void swap(nricMatrix& A, nricMatrix& B) { A.swap(B); }
00139 inline void swap(UpperTriangularMatrix& A, UpperTriangularMatrix& B)
00140    { A.swap(B); }
00141 inline void swap(LowerTriangularMatrix& A, LowerTriangularMatrix& B)
00142    { A.swap(B); }
00143 inline void swap(SymmetricMatrix& A, SymmetricMatrix& B) { A.swap(B); }
00144 inline void swap(DiagonalMatrix& A, DiagonalMatrix& B) { A.swap(B); }
00145 inline void swap(RowVector& A, RowVector& B) { A.swap(B); }
00146 inline void swap(ColumnVector& A, ColumnVector& B) { A.swap(B); }
00147 inline void swap(CroutMatrix& A, CroutMatrix& B) { A.swap(B); }
00148 inline void swap(BandMatrix& A, BandMatrix& B) { A.swap(B); }
00149 inline void swap(UpperBandMatrix& A, UpperBandMatrix& B) { A.swap(B); }
00150 inline void swap(LowerBandMatrix& A, LowerBandMatrix& B) { A.swap(B); }
00151 inline void swap(SymmetricBandMatrix& A, SymmetricBandMatrix& B) { A.swap(B); }
00152 inline void swap(BandLUMatrix& A, BandLUMatrix& B) { A.swap(B); }
00153 inline void swap(IdentityMatrix& A, IdentityMatrix& B) { A.swap(B); }
00154 inline void swap(GenericMatrix& A, GenericMatrix& B) { A.swap(B); }
00155 
00156 #ifdef OPT_COMPATIBLE                    // for compatibility with opt++
00157 
00158 inline Real Norm2(const ColumnVector& CV) { return CV.norm_Frobenius(); }
00159 inline Real Dot(ColumnVector& CV1, ColumnVector& CV2)
00160    { return dotproduct(CV1, CV2); }
00161 
00162 #endif
00163 
00164 // ************************** applications *****************************/
00165 
00166 
00167 void QRZT(Matrix&, LowerTriangularMatrix&);
00168 
00169 void QRZT(const Matrix&, Matrix&, Matrix&);
00170 
00171 void QRZ(Matrix&, UpperTriangularMatrix&);
00172 
00173 void QRZ(const Matrix&, Matrix&, Matrix&);
00174 
00175 inline void HHDecompose(Matrix& X, LowerTriangularMatrix& L)
00176 { QRZT(X,L); }
00177 
00178 inline void HHDecompose(const Matrix& X, Matrix& Y, Matrix& M)
00179 { QRZT(X, Y, M); }
00180 
00181 void updateQRZT(Matrix& X, LowerTriangularMatrix& L);
00182 
00183 void updateQRZ(Matrix& X, UpperTriangularMatrix& U);
00184 
00185 inline void UpdateQRZT(Matrix& X, LowerTriangularMatrix& L)
00186    { updateQRZT(X, L); }
00187 
00188 inline void UpdateQRZ(Matrix& X, UpperTriangularMatrix& U)
00189    { updateQRZ(X, U); }
00190 
00191 // Matrix A's first n columns are orthonormal
00192 // so A.Columns(1,n).t() * A.Columns(1,n) is the identity matrix.
00193 // Fill out the remaining columns of A to make them orthonormal
00194 // so A.t() * A is the identity matrix
00195 void extend_orthonormal(Matrix& A, int n);
00196 
00197 
00198 ReturnMatrix Cholesky(const SymmetricMatrix&);
00199 
00200 ReturnMatrix Cholesky(const SymmetricBandMatrix&);
00201 
00202 
00203 // produces the Cholesky decomposition of A + x.t() * x where A = chol.t() * chol
00204 // and x is a RowVector
00205 void update_Cholesky(UpperTriangularMatrix& chol, RowVector x);
00206 inline void UpdateCholesky(UpperTriangularMatrix& chol, const RowVector& x)
00207    { update_Cholesky(chol, x); }
00208 
00209 // produces the Cholesky decomposition of A - x.t() * x where A = chol.t() * chol
00210 // and x is a RowVector
00211 void downdate_Cholesky(UpperTriangularMatrix &chol, RowVector x);
00212 inline void DowndateCholesky(UpperTriangularMatrix &chol, const RowVector& x)
00213    { downdate_Cholesky(chol, x); }
00214 
00215 // a RIGHT circular shift of the rows and columns from
00216 // 1,...,k-1,k,k+1,...l,l+1,...,p to
00217 // 1,...,k-1,l,k,k+1,...l-1,l+1,...p
00218 void right_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l);
00219 inline void RightCircularUpdateCholesky(UpperTriangularMatrix &chol,
00220   int k, int l) { right_circular_update_Cholesky(chol, k, l); }
00221 
00222 // a LEFT circular shift of the rows and columns from
00223 // 1,...,k-1,k,k+1,...l,l+1,...,p to
00224 // 1,...,k-1,k+1,...l,k,l+1,...,p to
00225 void left_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l);
00226 inline void LeftCircularUpdateCholesky(UpperTriangularMatrix &chol,
00227    int k, int l) { left_circular_update_Cholesky(chol, k, l); }
00228 
00229 
00230 void SVD(const Matrix&, DiagonalMatrix&, Matrix&, Matrix&,
00231     bool=true, bool=true);
00232 
00233 void SVD(const Matrix&, DiagonalMatrix&);
00234 
00235 inline void SVD(const Matrix& A, DiagonalMatrix& D, Matrix& U,
00236    bool withU = true) { SVD(A, D, U, U, withU, false); }
00237 
00238 void SortSV(DiagonalMatrix& D, Matrix& U, bool ascending = false);
00239 
00240 void SortSV(DiagonalMatrix& D, Matrix& U, Matrix& V, bool ascending = false);
00241 
00242 void Jacobi(const SymmetricMatrix&, DiagonalMatrix&);
00243 
00244 void Jacobi(const SymmetricMatrix&, DiagonalMatrix&, SymmetricMatrix&);
00245 
00246 void Jacobi(const SymmetricMatrix&, DiagonalMatrix&, Matrix&);
00247 
00248 void Jacobi(const SymmetricMatrix&, DiagonalMatrix&, SymmetricMatrix&,
00249    Matrix&, bool=true);
00250 
00251 void eigenvalues(const SymmetricMatrix&, DiagonalMatrix&);
00252 
00253 void eigenvalues(const SymmetricMatrix&, DiagonalMatrix&, SymmetricMatrix&);
00254 
00255 void eigenvalues(const SymmetricMatrix&, DiagonalMatrix&, Matrix&);
00256 
00257 inline void EigenValues(const SymmetricMatrix& A, DiagonalMatrix& D)
00258    { eigenvalues(A, D); }
00259 
00260 inline void EigenValues(const SymmetricMatrix& A, DiagonalMatrix& D,
00261    SymmetricMatrix& S) { eigenvalues(A, D, S); }
00262 
00263 inline void EigenValues(const SymmetricMatrix& A, DiagonalMatrix& D, Matrix& V)
00264    { eigenvalues(A, D, V); }
00265 
00266 class SymmetricEigenAnalysis
00267 // not implemented yet
00268 {
00269 public:
00270    SymmetricEigenAnalysis(const SymmetricMatrix&);
00271 private:
00272    DiagonalMatrix diag;
00273    DiagonalMatrix offdiag;
00274    SymmetricMatrix backtransform;
00275    FREE_CHECK(SymmetricEigenAnalysis)
00276 };
00277 
00278 void sort_ascending(GeneralMatrix&);
00279 
00280 void sort_descending(GeneralMatrix&);
00281 
00282 inline void SortAscending(GeneralMatrix& gm) { sort_ascending(gm); }
00283 
00284 inline void SortDescending(GeneralMatrix& gm) { sort_descending(gm); }
00285 
00287 class FFT_Controller
00288 {
00289 public:
00290    static bool OnlyOldFFT;
00291    static bool ar_1d_ft (int PTS, Real* X, Real *Y);
00292    static bool CanFactor(int PTS);
00293 };
00294 
00295 void FFT(const ColumnVector&, const ColumnVector&,
00296    ColumnVector&, ColumnVector&);
00297 
00298 void FFTI(const ColumnVector&, const ColumnVector&,
00299    ColumnVector&, ColumnVector&);
00300 
00301 void RealFFT(const ColumnVector&, ColumnVector&, ColumnVector&);
00302 
00303 void RealFFTI(const ColumnVector&, const ColumnVector&, ColumnVector&);
00304 
00305 void DCT_II(const ColumnVector&, ColumnVector&);
00306 
00307 void DCT_II_inverse(const ColumnVector&, ColumnVector&);
00308 
00309 void DST_II(const ColumnVector&, ColumnVector&);
00310 
00311 void DST_II_inverse(const ColumnVector&, ColumnVector&);
00312 
00313 void DCT(const ColumnVector&, ColumnVector&);
00314 
00315 void DCT_inverse(const ColumnVector&, ColumnVector&);
00316 
00317 void DST(const ColumnVector&, ColumnVector&);
00318 
00319 void DST_inverse(const ColumnVector&, ColumnVector&);
00320 
00321 void FFT2(const Matrix& U, const Matrix& V, Matrix& X, Matrix& Y);
00322 
00323 void FFT2I(const Matrix& U, const Matrix& V, Matrix& X, Matrix& Y);
00324 
00325 
00326 // This class is used by the new FFT program
00327 
00328 // Suppose an integer is expressed as a sequence of digits with each
00329 // digit having a different radix.
00330 // This class supposes we are counting with this multi-radix number
00331 // but also keeps track of the number with the digits (and radices)
00332 // reversed.
00333 // The integer starts at zero
00334 // operator++() increases it by 1
00335 // Counter gives the number of increments
00336 // Reverse() gives the value with the digits in reverse order
00337 // Swap is true if reverse is less than counter
00338 // Finish is true when we have done a complete cycle and are back at zero
00339 
00340 class MultiRadixCounter
00341 {
00342    const SimpleIntArray& Radix;
00343                               // radix of each digit
00344                               // n-1 highest order, 0 lowest order
00345    SimpleIntArray& Value;     // value of each digit
00346    const int n;               // number of digits
00347    int reverse;               // value when order of digits is reversed
00348    int product;               // product of radices
00349    int counter;               // counter
00350    bool finish;               // true when we have gone over whole range
00351 public:
00352    MultiRadixCounter(int nx, const SimpleIntArray& rx,
00353       SimpleIntArray& vx);
00354    void operator++();         // increment the multi-radix counter
00355    bool Swap() const { return reverse < counter; }
00356    bool Finish() const { return finish; }
00357    int Reverse() const { return reverse; }
00358    int Counter() const { return counter; }
00359 };
00360 
00361 // multiplication by Helmert matrix
00362 ReturnMatrix Helmert(int n, bool full=false);
00363 ReturnMatrix Helmert(const ColumnVector& X, bool full=false);
00364 ReturnMatrix Helmert(int n, int j, bool full=false);
00365 ReturnMatrix Helmert_transpose(const ColumnVector& Y, bool full=false);
00366 Real Helmert_transpose(const ColumnVector& Y, int j, bool full=false);
00367 ReturnMatrix Helmert(const Matrix& X, bool full=false);
00368 ReturnMatrix Helmert_transpose(const Matrix& Y, bool full=false);
00369 
00370 
00371 
00372 // Copyright (C) 1991,2,3,4: R B Davies
00373 
00374 #ifndef WANT_STREAM
00375 #define WANT_STREAM
00376 #endif
00377 
00378 // **************************** input/output *****************************/
00379 
00380 ostream& operator<<(ostream&, const BaseMatrix&);
00381 
00382 ostream& operator<<(ostream&, const GeneralMatrix&);
00383 
00384 
00385 /*  Use in some old versions of G++ without complete iomanipulators
00386 
00387 class Omanip_precision
00388 {
00389    int x;
00390 public:
00391    Omanip_precision(int i) : x(i) {}
00392    friend ostream& operator<<(ostream& os, Omanip_precision i);
00393 };
00394 
00395 
00396 Omanip_precision setprecision(int i);
00397 
00398 class Omanip_width
00399 {
00400    int x;
00401 public:
00402    Omanip_width(int i) : x(i) {}
00403    friend ostream& operator<<(ostream& os, Omanip_width i);
00404 };
00405 
00406 Omanip_width setw(int i);
00407 
00408 */
00409 
00410 #ifdef use_namespace
00411 }
00412 #endif
00413 
00414 #ifdef use_namespace
00415 namespace ROBOOP {
00416   using namespace NEWMAT;
00417 #endif
00418 
00419 #ifndef M_PI
00420 #define M_PI 3.14159265358979
00421 #endif
00422 
00423 #define GRAVITY 9.81
00424 
00425 // global variables
00426 extern Real fourbyfourident[];
00427 extern Real threebythreeident[];
00428 
00429 // angle conversion
00430 inline double deg2rad(const double angle_deg){ return angle_deg*M_PI/180; }
00431 inline double rad2deg(const double angle_rad){ return angle_rad*180/M_PI; }
00432 
00433 // vector operation
00434 
00435 ReturnMatrix x_prod_matrix(const ColumnVector & x);
00436 
00437 ReturnMatrix pinv(const Matrix & M);
00438 
00439 // numerical analysis tools
00440 
00441 ReturnMatrix Integ_Trap(const ColumnVector & present, ColumnVector & past, const Real dt);
00442 
00443 void Runge_Kutta4(ReturnMatrix (*xdot)(Real time, const Matrix & xin),
00444                   const Matrix & xo, Real to, Real tf, int nsteps,
00445                   RowVector & tout, Matrix & xout);
00446 
00447 void Runge_Kutta4_Real_time(ReturnMatrix (*xdot)(Real time, const Matrix & xin),
00448                             const Matrix & xo, Real to, Real tf, int nsteps);
00449 
00450 void Runge_Kutta4_Real_time(ReturnMatrix (*xdot)(Real time, const Matrix & xin,
00451                             bool & exit, bool & init),
00452                             const Matrix & xo, Real to, Real tf, int nsteps);
00453 
00454 void odeint(ReturnMatrix (*xdot)(Real time, const Matrix & xin),
00455             Matrix & xo, Real to, Real tf, Real eps, Real h1, Real hmin,
00456             int & nok, int & nbad,
00457             RowVector & tout, Matrix & xout, Real dtsav);
00458 
00459 ReturnMatrix sign(const Matrix & x);
00460 
00461 short sign(const Real x);
00462 
00463 const double epsilon = 0.0000001;
00464 
00465 inline bool isZero(const double x)
00466 {
00467     if ( fabs(x) < epsilon)
00468     {
00469         return true;
00470     }
00471     return false;
00472 }
00473 
00474 
00475 // translation
00476 ReturnMatrix trans(const ColumnVector & a);
00477 
00478 // rotation matrices
00479 ReturnMatrix rotx(const Real alpha);
00480 ReturnMatrix roty(const Real beta);
00481 ReturnMatrix rotz(const Real gamma);
00482 ReturnMatrix rotk(const Real theta, const ColumnVector & k);
00483 
00484 ReturnMatrix rpy(const ColumnVector & a);
00485 ReturnMatrix eulzxz(const ColumnVector & a);
00486 ReturnMatrix rotd(const Real theta, const ColumnVector & k1, const ColumnVector & k2);
00487 
00488 
00489 // inverse on rotation matrices
00490 ReturnMatrix irotk(const Matrix & R);
00491 ReturnMatrix irpy(const Matrix & R);
00492 ReturnMatrix ieulzxz(const Matrix & R);
00493 
00494 #ifdef use_namespace
00495 }
00496 #endif
00497 
00498 #endif
00499 


lo
Author(s): U. Klank
autogenerated on Mon Oct 6 2014 10:44:13