$search
00001 00002 00003 00006 00007 // Copyright (C) 1991,2,3,4,8: R B Davies 00008 00009 #ifndef NEWMATAP_LIB 00010 #define NEWMATAP_LIB 0 00011 00012 #include "newmat.h" 00013 00014 #ifdef use_namespace 00015 namespace NEWMAT { 00016 #endif 00017 00018 00019 // ************************** applications *****************************/ 00020 00021 00022 void QRZT(Matrix&, LowerTriangularMatrix&); 00023 00024 void QRZT(const Matrix&, Matrix&, Matrix&); 00025 00026 void QRZ(Matrix&, UpperTriangularMatrix&); 00027 00028 void QRZ(const Matrix&, Matrix&, Matrix&); 00029 00030 inline void HHDecompose(Matrix& X, LowerTriangularMatrix& L) 00031 { QRZT(X,L); } 00032 00033 inline void HHDecompose(const Matrix& X, Matrix& Y, Matrix& M) 00034 { QRZT(X, Y, M); } 00035 00036 void updateQRZT(Matrix& X, LowerTriangularMatrix& L); 00037 00038 void updateQRZ(Matrix& X, UpperTriangularMatrix& U); 00039 00040 inline void UpdateQRZT(Matrix& X, LowerTriangularMatrix& L) 00041 { updateQRZT(X, L); } 00042 00043 inline void UpdateQRZ(Matrix& X, UpperTriangularMatrix& U) 00044 { updateQRZ(X, U); } 00045 00046 // Matrix A's first n columns are orthonormal 00047 // so A.Columns(1,n).t() * A.Columns(1,n) is the identity matrix. 00048 // Fill out the remaining columns of A to make them orthonormal 00049 // so A.t() * A is the identity matrix 00050 void extend_orthonormal(Matrix& A, int n); 00051 00052 00053 ReturnMatrix Cholesky(const SymmetricMatrix&); 00054 00055 ReturnMatrix Cholesky(const SymmetricBandMatrix&); 00056 00057 00058 // produces the Cholesky decomposition of A + x.t() * x where A = chol.t() * chol 00059 // and x is a RowVector 00060 void update_Cholesky(UpperTriangularMatrix& chol, RowVector x); 00061 inline void UpdateCholesky(UpperTriangularMatrix& chol, const RowVector& x) 00062 { update_Cholesky(chol, x); } 00063 00064 // produces the Cholesky decomposition of A - x.t() * x where A = chol.t() * chol 00065 // and x is a RowVector 00066 void downdate_Cholesky(UpperTriangularMatrix &chol, RowVector x); 00067 inline void DowndateCholesky(UpperTriangularMatrix &chol, const RowVector& x) 00068 { downdate_Cholesky(chol, x); } 00069 00070 // a RIGHT circular shift of the rows and columns from 00071 // 1,...,k-1,k,k+1,...l,l+1,...,p to 00072 // 1,...,k-1,l,k,k+1,...l-1,l+1,...p 00073 void right_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l); 00074 inline void RightCircularUpdateCholesky(UpperTriangularMatrix &chol, 00075 int k, int l) { right_circular_update_Cholesky(chol, k, l); } 00076 00077 // a LEFT circular shift of the rows and columns from 00078 // 1,...,k-1,k,k+1,...l,l+1,...,p to 00079 // 1,...,k-1,k+1,...l,k,l+1,...,p to 00080 void left_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l); 00081 inline void LeftCircularUpdateCholesky(UpperTriangularMatrix &chol, 00082 int k, int l) { left_circular_update_Cholesky(chol, k, l); } 00083 00084 00085 void SVD(const Matrix&, DiagonalMatrix&, Matrix&, Matrix&, 00086 bool=true, bool=true); 00087 00088 void SVD(const Matrix&, DiagonalMatrix&); 00089 00090 inline void SVD(const Matrix& A, DiagonalMatrix& D, Matrix& U, 00091 bool withU = true) { SVD(A, D, U, U, withU, false); } 00092 00093 void SortSV(DiagonalMatrix& D, Matrix& U, bool ascending = false); 00094 00095 void SortSV(DiagonalMatrix& D, Matrix& U, Matrix& V, bool ascending = false); 00096 00097 void Jacobi(const SymmetricMatrix&, DiagonalMatrix&); 00098 00099 void Jacobi(const SymmetricMatrix&, DiagonalMatrix&, SymmetricMatrix&); 00100 00101 void Jacobi(const SymmetricMatrix&, DiagonalMatrix&, Matrix&); 00102 00103 void Jacobi(const SymmetricMatrix&, DiagonalMatrix&, SymmetricMatrix&, 00104 Matrix&, bool=true); 00105 00106 void eigenvalues(const SymmetricMatrix&, DiagonalMatrix&); 00107 00108 void eigenvalues(const SymmetricMatrix&, DiagonalMatrix&, SymmetricMatrix&); 00109 00110 void eigenvalues(const SymmetricMatrix&, DiagonalMatrix&, Matrix&); 00111 00112 inline void EigenValues(const SymmetricMatrix& A, DiagonalMatrix& D) 00113 { eigenvalues(A, D); } 00114 00115 inline void EigenValues(const SymmetricMatrix& A, DiagonalMatrix& D, 00116 SymmetricMatrix& S) { eigenvalues(A, D, S); } 00117 00118 inline void EigenValues(const SymmetricMatrix& A, DiagonalMatrix& D, Matrix& V) 00119 { eigenvalues(A, D, V); } 00120 00121 class SymmetricEigenAnalysis 00122 // not implemented yet 00123 { 00124 public: 00125 SymmetricEigenAnalysis(const SymmetricMatrix&); 00126 private: 00127 DiagonalMatrix diag; 00128 DiagonalMatrix offdiag; 00129 SymmetricMatrix backtransform; 00130 FREE_CHECK(SymmetricEigenAnalysis) 00131 }; 00132 00133 void sort_ascending(GeneralMatrix&); 00134 00135 void sort_descending(GeneralMatrix&); 00136 00137 inline void SortAscending(GeneralMatrix& gm) { sort_ascending(gm); } 00138 00139 inline void SortDescending(GeneralMatrix& gm) { sort_descending(gm); } 00140 00142 class FFT_Controller 00143 { 00144 public: 00145 static bool OnlyOldFFT; 00146 static bool ar_1d_ft (int PTS, Real* X, Real *Y); 00147 static bool CanFactor(int PTS); 00148 }; 00149 00150 void FFT(const ColumnVector&, const ColumnVector&, 00151 ColumnVector&, ColumnVector&); 00152 00153 void FFTI(const ColumnVector&, const ColumnVector&, 00154 ColumnVector&, ColumnVector&); 00155 00156 void RealFFT(const ColumnVector&, ColumnVector&, ColumnVector&); 00157 00158 void RealFFTI(const ColumnVector&, const ColumnVector&, ColumnVector&); 00159 00160 void DCT_II(const ColumnVector&, ColumnVector&); 00161 00162 void DCT_II_inverse(const ColumnVector&, ColumnVector&); 00163 00164 void DST_II(const ColumnVector&, ColumnVector&); 00165 00166 void DST_II_inverse(const ColumnVector&, ColumnVector&); 00167 00168 void DCT(const ColumnVector&, ColumnVector&); 00169 00170 void DCT_inverse(const ColumnVector&, ColumnVector&); 00171 00172 void DST(const ColumnVector&, ColumnVector&); 00173 00174 void DST_inverse(const ColumnVector&, ColumnVector&); 00175 00176 void FFT2(const Matrix& U, const Matrix& V, Matrix& X, Matrix& Y); 00177 00178 void FFT2I(const Matrix& U, const Matrix& V, Matrix& X, Matrix& Y); 00179 00180 00181 // This class is used by the new FFT program 00182 00183 // Suppose an integer is expressed as a sequence of digits with each 00184 // digit having a different radix. 00185 // This class supposes we are counting with this multi-radix number 00186 // but also keeps track of the number with the digits (and radices) 00187 // reversed. 00188 // The integer starts at zero 00189 // operator++() increases it by 1 00190 // Counter gives the number of increments 00191 // Reverse() gives the value with the digits in reverse order 00192 // Swap is true if reverse is less than counter 00193 // Finish is true when we have done a complete cycle and are back at zero 00194 00195 class MultiRadixCounter 00196 { 00197 const SimpleIntArray& Radix; 00198 // radix of each digit 00199 // n-1 highest order, 0 lowest order 00200 SimpleIntArray& Value; // value of each digit 00201 const int n; // number of digits 00202 int reverse; // value when order of digits is reversed 00203 int product; // product of radices 00204 int counter; // counter 00205 bool finish; // true when we have gone over whole range 00206 public: 00207 MultiRadixCounter(int nx, const SimpleIntArray& rx, 00208 SimpleIntArray& vx); 00209 void operator++(); // increment the multi-radix counter 00210 bool Swap() const { return reverse < counter; } 00211 bool Finish() const { return finish; } 00212 int Reverse() const { return reverse; } 00213 int Counter() const { return counter; } 00214 }; 00215 00216 // multiplication by Helmert matrix 00217 ReturnMatrix Helmert(int n, bool full=false); 00218 ReturnMatrix Helmert(const ColumnVector& X, bool full=false); 00219 ReturnMatrix Helmert(int n, int j, bool full=false); 00220 ReturnMatrix Helmert_transpose(const ColumnVector& Y, bool full=false); 00221 Real Helmert_transpose(const ColumnVector& Y, int j, bool full=false); 00222 ReturnMatrix Helmert(const Matrix& X, bool full=false); 00223 ReturnMatrix Helmert_transpose(const Matrix& Y, bool full=false); 00224 00225 00226 00227 00228 #ifdef use_namespace 00229 } 00230 #endif 00231 00232 00233 00234 #endif 00235 00236 // body file: cholesky.cpp 00237 // body file: evalue.cpp 00238 // body file: fft.cpp 00239 // body file: hholder.cpp 00240 // body file: jacobi.cpp 00241 // body file: newfft.cpp 00242 // body file: sort.cpp 00243 // body file: svd.cpp 00244 // body file: nm_misc.cpp 00245 00246 00247 00249 00250