source: ntrip/trunk/BNS/newmat/newmatap.h

Last change on this file was 810, checked in by mervart, 17 years ago

* empty log message *

File size: 7.3 KB
Line 
1/// \ingroup newmat
2///@{
3
4/// \file newmatap.h
5/// Definition file for advanced matrix functions.
6
7// Copyright (C) 1991,2,3,4,8: R B Davies
8
9#ifndef NEWMATAP_LIB
10#define NEWMATAP_LIB 0
11
12#include "newmat.h"
13
14#ifdef use_namespace
15namespace NEWMAT {
16#endif
17
18
19// ************************** applications *****************************/
20
21
22void QRZT(Matrix&, LowerTriangularMatrix&);
23
24void QRZT(const Matrix&, Matrix&, Matrix&);
25
26void QRZ(Matrix&, UpperTriangularMatrix&);
27
28void QRZ(const Matrix&, Matrix&, Matrix&);
29
30inline void HHDecompose(Matrix& X, LowerTriangularMatrix& L)
31{ QRZT(X,L); }
32
33inline void HHDecompose(const Matrix& X, Matrix& Y, Matrix& M)
34{ QRZT(X, Y, M); }
35
36void updateQRZT(Matrix& X, LowerTriangularMatrix& L);
37
38void updateQRZ(Matrix& X, UpperTriangularMatrix& U);
39
40inline void UpdateQRZT(Matrix& X, LowerTriangularMatrix& L)
41 { updateQRZT(X, L); }
42
43inline void UpdateQRZ(Matrix& X, UpperTriangularMatrix& U)
44 { updateQRZ(X, U); }
45
46// Matrix A's first n columns are orthonormal
47// so A.Columns(1,n).t() * A.Columns(1,n) is the identity matrix.
48// Fill out the remaining columns of A to make them orthonormal
49// so A.t() * A is the identity matrix
50void extend_orthonormal(Matrix& A, int n);
51
52
53ReturnMatrix Cholesky(const SymmetricMatrix&);
54
55ReturnMatrix Cholesky(const SymmetricBandMatrix&);
56
57
58// produces the Cholesky decomposition of A + x.t() * x where A = chol.t() * chol
59// and x is a RowVector
60void update_Cholesky(UpperTriangularMatrix& chol, RowVector x);
61inline void UpdateCholesky(UpperTriangularMatrix& chol, const RowVector& x)
62 { update_Cholesky(chol, x); }
63
64// produces the Cholesky decomposition of A - x.t() * x where A = chol.t() * chol
65// and x is a RowVector
66void downdate_Cholesky(UpperTriangularMatrix &chol, RowVector x);
67inline void DowndateCholesky(UpperTriangularMatrix &chol, const RowVector& x)
68 { downdate_Cholesky(chol, x); }
69
70// a RIGHT circular shift of the rows and columns from
71// 1,...,k-1,k,k+1,...l,l+1,...,p to
72// 1,...,k-1,l,k,k+1,...l-1,l+1,...p
73void right_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l);
74inline void RightCircularUpdateCholesky(UpperTriangularMatrix &chol,
75 int k, int l) { right_circular_update_Cholesky(chol, k, l); }
76
77// a LEFT circular shift of the rows and columns from
78// 1,...,k-1,k,k+1,...l,l+1,...,p to
79// 1,...,k-1,k+1,...l,k,l+1,...,p to
80void left_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l);
81inline void LeftCircularUpdateCholesky(UpperTriangularMatrix &chol,
82 int k, int l) { left_circular_update_Cholesky(chol, k, l); }
83
84
85void SVD(const Matrix&, DiagonalMatrix&, Matrix&, Matrix&,
86 bool=true, bool=true);
87
88void SVD(const Matrix&, DiagonalMatrix&);
89
90inline void SVD(const Matrix& A, DiagonalMatrix& D, Matrix& U,
91 bool withU = true) { SVD(A, D, U, U, withU, false); }
92
93void SortSV(DiagonalMatrix& D, Matrix& U, bool ascending = false);
94
95void SortSV(DiagonalMatrix& D, Matrix& U, Matrix& V, bool ascending = false);
96
97void Jacobi(const SymmetricMatrix&, DiagonalMatrix&);
98
99void Jacobi(const SymmetricMatrix&, DiagonalMatrix&, SymmetricMatrix&);
100
101void Jacobi(const SymmetricMatrix&, DiagonalMatrix&, Matrix&);
102
103void Jacobi(const SymmetricMatrix&, DiagonalMatrix&, SymmetricMatrix&,
104 Matrix&, bool=true);
105
106void eigenvalues(const SymmetricMatrix&, DiagonalMatrix&);
107
108void eigenvalues(const SymmetricMatrix&, DiagonalMatrix&, SymmetricMatrix&);
109
110void eigenvalues(const SymmetricMatrix&, DiagonalMatrix&, Matrix&);
111
112inline void EigenValues(const SymmetricMatrix& A, DiagonalMatrix& D)
113 { eigenvalues(A, D); }
114
115inline void EigenValues(const SymmetricMatrix& A, DiagonalMatrix& D,
116 SymmetricMatrix& S) { eigenvalues(A, D, S); }
117
118inline void EigenValues(const SymmetricMatrix& A, DiagonalMatrix& D, Matrix& V)
119 { eigenvalues(A, D, V); }
120
121class SymmetricEigenAnalysis
122// not implemented yet
123{
124public:
125 SymmetricEigenAnalysis(const SymmetricMatrix&);
126private:
127 DiagonalMatrix diag;
128 DiagonalMatrix offdiag;
129 SymmetricMatrix backtransform;
130 FREE_CHECK(SymmetricEigenAnalysis)
131};
132
133void sort_ascending(GeneralMatrix&);
134
135void sort_descending(GeneralMatrix&);
136
137inline void SortAscending(GeneralMatrix& gm) { sort_ascending(gm); }
138
139inline void SortDescending(GeneralMatrix& gm) { sort_descending(gm); }
140
141/// Decide which fft method to use and carry out new fft function
142class FFT_Controller
143{
144public:
145 static bool OnlyOldFFT;
146 static bool ar_1d_ft (int PTS, Real* X, Real *Y);
147 static bool CanFactor(int PTS);
148};
149
150void FFT(const ColumnVector&, const ColumnVector&,
151 ColumnVector&, ColumnVector&);
152
153void FFTI(const ColumnVector&, const ColumnVector&,
154 ColumnVector&, ColumnVector&);
155
156void RealFFT(const ColumnVector&, ColumnVector&, ColumnVector&);
157
158void RealFFTI(const ColumnVector&, const ColumnVector&, ColumnVector&);
159
160void DCT_II(const ColumnVector&, ColumnVector&);
161
162void DCT_II_inverse(const ColumnVector&, ColumnVector&);
163
164void DST_II(const ColumnVector&, ColumnVector&);
165
166void DST_II_inverse(const ColumnVector&, ColumnVector&);
167
168void DCT(const ColumnVector&, ColumnVector&);
169
170void DCT_inverse(const ColumnVector&, ColumnVector&);
171
172void DST(const ColumnVector&, ColumnVector&);
173
174void DST_inverse(const ColumnVector&, ColumnVector&);
175
176void FFT2(const Matrix& U, const Matrix& V, Matrix& X, Matrix& Y);
177
178void FFT2I(const Matrix& U, const Matrix& V, Matrix& X, Matrix& Y);
179
180
181// This class is used by the new FFT program
182
183// Suppose an integer is expressed as a sequence of digits with each
184// digit having a different radix.
185// This class supposes we are counting with this multi-radix number
186// but also keeps track of the number with the digits (and radices)
187// reversed.
188// The integer starts at zero
189// operator++() increases it by 1
190// Counter gives the number of increments
191// Reverse() gives the value with the digits in reverse order
192// Swap is true if reverse is less than counter
193// Finish is true when we have done a complete cycle and are back at zero
194
195class MultiRadixCounter
196{
197 const SimpleIntArray& Radix;
198 // radix of each digit
199 // n-1 highest order, 0 lowest order
200 SimpleIntArray& Value; // value of each digit
201 const int n; // number of digits
202 int reverse; // value when order of digits is reversed
203 int product; // product of radices
204 int counter; // counter
205 bool finish; // true when we have gone over whole range
206public:
207 MultiRadixCounter(int nx, const SimpleIntArray& rx,
208 SimpleIntArray& vx);
209 void operator++(); // increment the multi-radix counter
210 bool Swap() const { return reverse < counter; }
211 bool Finish() const { return finish; }
212 int Reverse() const { return reverse; }
213 int Counter() const { return counter; }
214};
215
216// multiplication by Helmert matrix
217ReturnMatrix Helmert(int n, bool full=false);
218ReturnMatrix Helmert(const ColumnVector& X, bool full=false);
219ReturnMatrix Helmert(int n, int j, bool full=false);
220ReturnMatrix Helmert_transpose(const ColumnVector& Y, bool full=false);
221Real Helmert_transpose(const ColumnVector& Y, int j, bool full=false);
222ReturnMatrix Helmert(const Matrix& X, bool full=false);
223ReturnMatrix Helmert_transpose(const Matrix& Y, bool full=false);
224
225
226
227
228#ifdef use_namespace
229}
230#endif
231
232
233
234#endif
235
236// body file: cholesky.cpp
237// body file: evalue.cpp
238// body file: fft.cpp
239// body file: hholder.cpp
240// body file: jacobi.cpp
241// body file: newfft.cpp
242// body file: sort.cpp
243// body file: svd.cpp
244// body file: nm_misc.cpp
245
246
247
248///@}
249
250
Note: See TracBrowser for help on using the repository browser.