00001
00002 #ifndef vnl_diag_matrix_h_
00003 #define vnl_diag_matrix_h_
00004 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00005 #pragma interface
00006 #endif
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include <vcl_cassert.h>
00024 #include <vcl_iosfwd.h>
00025 #include <vnl/vnl_vector.h>
00026 #include <vnl/vnl_matrix.h>
00027
00028
00029 template <class T> class vnl_diag_matrix;
00030 template <class T> vnl_vector<T> operator*(vnl_diag_matrix<T> const&, vnl_vector<T> const&);
00031
00032
00033
00034
00035
00036
00037
00038 export
00039 template <class T>
00040 class vnl_diag_matrix
00041 {
00042 vnl_vector<T> diagonal_;
00043
00044 public:
00045 vnl_diag_matrix() : diagonal_() {}
00046
00047
00048 vnl_diag_matrix(unsigned nn) : diagonal_(nn) {}
00049
00050
00051 vnl_diag_matrix(unsigned nn, T const& value) : diagonal_(nn, value) {}
00052
00053
00054
00055 vnl_diag_matrix(vnl_vector<T> const& that): diagonal_(that) {}
00056 ~vnl_diag_matrix() {}
00057
00058 inline vnl_diag_matrix& operator=(vnl_diag_matrix<T> const& that) {
00059 this->diagonal_ = that.diagonal_;
00060 return *this;
00061 }
00062
00063
00064
00065
00066 inline vnl_diag_matrix<T>& operator*=(T v) { diagonal_ *= v; return *this; }
00067
00068 inline vnl_diag_matrix<T>& operator/=(T v) { diagonal_ /= v; return *this; }
00069
00070
00071
00072 vnl_diag_matrix& invert_in_place();
00073 T determinant() const;
00074 vnl_vector<T> solve(vnl_vector<T> const& b) const;
00075 void solve(vnl_vector<T> const& b, vnl_vector<T>* out) const;
00076
00077
00078
00079 inline T operator () (unsigned i, unsigned j) const {
00080 return (i != j) ? T(0) : diagonal_[i];
00081 }
00082
00083 inline T& operator () (unsigned i, unsigned j) {
00084 assert(i == j);
00085 return diagonal_[i];
00086 }
00087 inline T& operator() (unsigned i) { return diagonal_[i]; }
00088 inline T const& operator() (unsigned i) const { return diagonal_[i]; }
00089
00090 inline T& operator[] (unsigned i) { return diagonal_[i]; }
00091 inline T const& operator[] (unsigned i) const { return diagonal_[i]; }
00092
00093
00094 inline void put (unsigned r, unsigned c, T const& v) {
00095 assert(r == c); assert (r<size()); diagonal_[r] = v;
00096 }
00097
00098
00099 inline T get (unsigned r, unsigned c) const {
00100 assert(r == c); assert (r<size()); return diagonal_[r];
00101 }
00102
00103
00104 inline vnl_vector<T> get_diagonal() const { return diagonal_; }
00105
00106
00107 inline vnl_vector<T> const& diagonal() const { return diagonal_; }
00108
00109
00110 inline vnl_diag_matrix& fill_diagonal (T const& v) { diagonal_.fill(v); return *this; }
00111
00112
00113 inline vnl_diag_matrix& set_diagonal(vnl_vector<T> const& v) { diagonal_ = v; return *this; }
00114
00115
00116
00117 typedef typename vnl_vector<T>::iterator iterator;
00118 inline iterator begin() { return diagonal_.begin(); }
00119 inline iterator end() { return diagonal_.end(); }
00120 typedef typename vnl_vector<T>::const_iterator const_iterator;
00121 inline const_iterator begin() const { return diagonal_.begin(); }
00122 inline const_iterator end() const { return diagonal_.end(); }
00123
00124 inline unsigned size() const { return diagonal_.size(); }
00125 inline unsigned rows() const { return diagonal_.size(); }
00126 inline unsigned cols() const { return diagonal_.size(); }
00127 inline unsigned columns() const { return diagonal_.size(); }
00128
00129
00130 inline vnl_matrix<T> asMatrix() const;
00131
00132 inline vnl_matrix<T> as_ref() const { return asMatrix(); }
00133
00134
00135 inline operator vnl_matrix<T> () const { return asMatrix(); }
00136
00137 inline void set_size(int n) { diagonal_.set_size(n); }
00138
00139 inline void clear() { diagonal_.clear(); }
00140 inline vnl_diag_matrix& fill(T const &x) { diagonal_.fill(x); return *this; }
00141
00142
00143 inline T* data_block() { return diagonal_.data_block(); }
00144 inline T const* data_block() const { return diagonal_.data_block(); }
00145
00146
00147 inline vnl_diag_matrix& set(vnl_vector<T> const& v) { diagonal_=v; return *this; }
00148
00149 private:
00150 #if VCL_NEED_FRIEND_FOR_TEMPLATE_OVERLOAD
00151 friend vnl_vector<T> operator*(vnl_diag_matrix<T> const&,vnl_vector<T> const&);
00152 #endif
00153 };
00154
00155
00156
00157 template <class T>
00158 vcl_ostream& operator<< (vcl_ostream&, vnl_diag_matrix<T> const&);
00159
00160
00161 template <class T>
00162 inline vnl_matrix<T> vnl_diag_matrix<T>::asMatrix() const
00163 {
00164 unsigned len = diagonal_.size();
00165 vnl_matrix<T> ret(len, len);
00166 for (unsigned i = 0; i < len; ++i)
00167 {
00168 unsigned j;
00169 for (j = 0; j < i; ++j)
00170 ret(i,j) = T(0);
00171 for (j = i+1; j < len; ++j)
00172 ret(i,j) = T(0);
00173 ret(i,i) = diagonal_[i];
00174 }
00175 return ret;
00176 }
00177
00178
00179
00180 template <class T>
00181 inline vnl_diag_matrix<T>& vnl_diag_matrix<T>::invert_in_place()
00182 {
00183 unsigned len = diagonal_.size();
00184 T* d = data_block();
00185 T one = T(1);
00186 for (unsigned i = 0; i < len; ++i)
00187 d[i] = one / d[i];
00188 return *this;
00189 }
00190
00191
00192 template <class T>
00193 inline T vnl_diag_matrix<T>::determinant() const
00194 {
00195 T det = T(1);
00196 T const* d = data_block();
00197 unsigned len = diagonal_.size();
00198 for (unsigned i = 0; i < len; ++i)
00199 det *= d[i];
00200 return det;
00201 }
00202
00203
00204
00205 template <class T>
00206 inline vnl_diag_matrix<T> operator* (vnl_diag_matrix<T> const& A, vnl_diag_matrix<T> const& B)
00207 {
00208 assert(A.size() == B.size());
00209 vnl_diag_matrix<T> ret = A;
00210 for (unsigned i = 0; i < A.size(); ++i)
00211 ret(i,i) *= B(i,i);
00212 return ret;
00213 }
00214
00215
00216
00217
00218 template <class T>
00219 inline vnl_matrix<T> operator* (vnl_matrix<T> const& A, vnl_diag_matrix<T> const& D)
00220 {
00221 assert(A.columns() == D.size());
00222 vnl_matrix<T> ret(A.rows(), A.columns());
00223 for (unsigned i = 0; i < A.rows(); ++i)
00224 for (unsigned j = 0; j < A.columns(); ++j)
00225 ret(i,j) = A(i,j) * D(j,j);
00226 return ret;
00227 }
00228
00229
00230
00231
00232 template <class T>
00233 inline vnl_matrix<T> operator* (vnl_diag_matrix<T> const& D, vnl_matrix<T> const& A)
00234 {
00235 assert(A.rows() == D.size());
00236 vnl_matrix<T> ret(A.rows(), A.columns());
00237 T const* d = D.data_block();
00238 for (unsigned i = 0; i < A.rows(); ++i)
00239 for (unsigned j = 0; j < A.columns(); ++j)
00240 ret(i,j) = A(i,j) * d[i];
00241 return ret;
00242 }
00243
00244
00245
00246 template <class T>
00247 inline vnl_diag_matrix<T> operator+ (vnl_diag_matrix<T> const& A, vnl_diag_matrix<T> const& B)
00248 {
00249 assert(A.size() == B.size());
00250 vnl_diag_matrix<T> ret = A;
00251 for (unsigned i = 0; i < A.size(); ++i)
00252 ret(i,i) += B(i,i);
00253 return ret;
00254 }
00255
00256
00257
00258
00259 template <class T>
00260 inline vnl_matrix<T> operator+ (vnl_matrix<T> const& A, vnl_diag_matrix<T> const& D)
00261 {
00262 const unsigned n = D.size();
00263 assert(A.rows() == n); assert(A.columns() == n);
00264 vnl_matrix<T> ret(A);
00265 T const* d = D.data_block();
00266 for (unsigned j = 0; j < n; ++j)
00267 ret(j,j) += d[j];
00268 return ret;
00269 }
00270
00271
00272
00273
00274 template <class T>
00275 inline vnl_matrix<T> operator+ (vnl_diag_matrix<T> const& D, vnl_matrix<T> const& A)
00276 {
00277 return A + D;
00278 }
00279
00280
00281
00282 template <class T>
00283 inline vnl_diag_matrix<T> operator- (vnl_diag_matrix<T> const& A, vnl_diag_matrix<T> const& B)
00284 {
00285 assert(A.size() == B.size());
00286 vnl_diag_matrix<T> ret = A;
00287 for (unsigned i = 0; i < A.size(); ++i)
00288 ret(i,i) -= B(i,i);
00289 return ret;
00290 }
00291
00292
00293
00294
00295 template <class T>
00296 inline vnl_matrix<T> operator- (vnl_matrix<T> const& A, vnl_diag_matrix<T> const& D)
00297 {
00298 const unsigned n = D.size();
00299 assert(A.rows() == n); assert(A.columns() == n);
00300 vnl_matrix<T> ret(A);
00301 T const* d = D.data_block();
00302 for (unsigned j = 0; j < n; ++j)
00303 ret(j,j) -= d[j];
00304 return ret;
00305 }
00306
00307
00308
00309
00310 template <class T>
00311 inline vnl_matrix<T> operator- (vnl_diag_matrix<T> const& D, vnl_matrix<T> const& A)
00312 {
00313 const unsigned n = D.size();
00314 assert(A.rows() == n); assert(A.columns() == n);
00315 vnl_matrix<T> ret(n, n);
00316 T const* d = D.data_block();
00317 for (unsigned i = 0; i < n; ++i)
00318 {
00319 for (unsigned j = 0; j < i; ++j)
00320 ret(i,j) = -A(i,j);
00321 for (unsigned j = i+1; j < n; ++j)
00322 ret(i,j) = -A(i,j);
00323 ret(i,i) = d[i] - A(i,i);
00324 }
00325 return ret;
00326 }
00327
00328
00329
00330
00331 template <class T>
00332 inline vnl_vector<T> operator* (vnl_diag_matrix<T> const& D, vnl_vector<T> const& A)
00333 {
00334 assert(A.size() == D.size());
00335 return element_product(D.diagonal(), A);
00336 }
00337
00338
00339
00340
00341 template <class T>
00342 inline vnl_vector<T> operator* (vnl_vector<T> const& A, vnl_diag_matrix<T> const& D)
00343 {
00344 assert(A.size() == D.size());
00345 return element_product(D.diagonal(), A);
00346 }
00347
00348 #endif // vnl_diag_matrix_h_