00001 #ifndef TNT_SPARSE_MATRIX_H
00002 #define TNT_SPARSE_MATRIX_H
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include <vector>
00026 #include "tnt_vector.h"
00027 #include "tnt_sparse_vector.h"
00028
00029 namespace TNT
00030 {
00031
00032
00033
00034
00035
00036 #if 0
00037 template <class T, class Integer>
00038 class Sparse_Matrix_Coordinate_Element
00039 {
00040 private:
00041
00042 T val_;
00043 Integer row_index_;
00044 Integer col_index_;
00045
00046 public:
00047
00048 Sparse_Matrix_Coordinate_Element(
00049 const T& a, const Integer &i, const Integer &j) :
00050 val_(a), row_index_(i), col_index(i) {}
00051
00052
00053 const T& value() const { return val_; }
00054 Integer row_index() const { return row_index_; }
00055 Integer col_index() const { return col_index_; }
00056
00057 T& value() { return val_; }
00058 Integer& row_index() { return row_index_; }
00059 Integer& col_index() { return col_index_; }
00060
00061 };
00062 #endif
00063
00081 template <class T>
00082 class Sparse_Matrix
00083 {
00084 private:
00085
00086
00087
00088 std::vector< Sparse_Vector<T> > S_;
00089
00090 int num_rows_;
00091 int num_cols_;
00092 int num_nonzeros_;
00093
00094
00095
00096
00097
00098 int internal_state_;
00099
00100
00101
00102 public:
00103
00104
00105 Sparse_Matrix(Subscript M, Subscript N):
00106 S_(M),
00107 num_rows_(M), num_cols_(N), num_nonzeros_(0),
00108 internal_state_(1) {};
00109
00110
00111
00112 Sparse_Matrix(Subscript M, Subscript N, Subscript nz, const T* val,
00113 const Subscript *r, const Subscript *c):
00114 S_(M),
00115 num_rows_(M),
00116 num_cols_(N),
00117 num_nonzeros_(0),
00118 internal_state_(1)
00119 {
00120
00121 insert(nz, val, r, c);
00122 close();
00123 };
00124
00125
00126 int is_closed() { return internal_state_; }
00127
00128 void insert(const T& val, Subscript i, Subscript j)
00129 {
00130 if (internal_state_ == 0) return;
00131
00132
00133 S_[i].insert(val, j);
00134 num_nonzeros_++;
00135 }
00136
00137 void insert(Subscript nz, const T* val, const Subscript *i,
00138 const Subscript *j)
00139 {
00140 if (internal_state_ == 0) return;
00141
00142 for (int count=0; count<nz; count++)
00143 {
00144 insert(val[count], i[count], j[count]);
00145 }
00146 }
00147
00148 void insert_base_one(const T& val, Subscript i, Subscript j)
00149 {
00150 insert_one_base(val, i-1, j-1);
00151 }
00152
00153 void insert_base_one(Subscript nz, const T* val, const Subscript *i,
00154 const Subscript *j)
00155 {
00156 for (int count=0; count<nz; count++)
00157 {
00158 insert(val[count], i[count]-1, j[count]-1);
00159 }
00160 }
00161
00162
00163 void close()
00164 {
00165
00166
00167
00168
00169
00170
00171 internal_state_ = 0;
00172 }
00173
00174 inline int num_rows() const {return num_rows_;}
00175 inline int num_cols() const {return num_cols_;}
00176 inline int num_columns() const {return num_cols_;}
00177 int num_nonzeros() const {return num_nonzeros_;}
00178
00179
00180 Vector<T> diag() const
00181 {
00182 int minMN = num_rows() < num_columns() ? num_rows() : num_columns();
00183 Vector<T> diag_(minMN, T(0));
00184
00185 for (int i=0; i<minMN; i++)
00186 {
00187 for (typename Sparse_Vector<T>::const_iterator p = S_[i].begin();
00188 p < S_[i].end(); p++ )
00189 {
00190 if (p->index() == i)
00191 {
00192 diag_[i] += p->value();
00193 break;
00194 }
00195 }
00196 }
00197 return diag_;
00198 }
00199
00200 Vector<T> mult(const Vector<T> &x) const
00201 {
00202 int M = num_rows();
00203 Vector<T> y(M);
00204 for (int i=0; i<M; i++)
00205 {
00206 y[i] = dot_product(S_[i], x);
00207 }
00208 return y;
00209 }
00210
00211
00212 inline double norm() const
00213 {
00214 T sum(0.0);
00215 for (int i=0; i<num_rows_; i++)
00216 {
00217 for (typename Sparse_Vector<T>::const_iterator p = S_[i].begin();
00218 p < S_[i].end(); p++ )
00219 {
00220 sum += p->value() * p->value();
00221 }
00222 }
00223 return sqrt(sum);
00224 }
00225
00226 std::ostream & print(std::ostream &s) const
00227 {
00228 for (int i=0; i<num_rows_; i++)
00229 {
00230 for (typename Sparse_Vector<T>::const_iterator p = S_[i].begin();
00231 p < S_[i].end(); p++ )
00232 {
00233 s << "( " << p->value() << " , " << i << ", " << p->index() << " )\n";
00234 }
00235 }
00236
00237 return s;
00238 }
00239
00240 std::ostream & print_base_one(std::ostream &s) const
00241 {
00242 for (int i=0; i<num_rows_; i++)
00243 {
00244 for (typename Sparse_Vector<T>::const_iterator p = S_[i].begin();
00245 p < S_[i].end(); p++ )
00246 {
00247 s <<"( "<<p->value()<<" , "<<i+1<<", "<< p->index()+1 << " )\n";
00248 }
00249 }
00250
00251 return s;
00252 }
00253
00254
00255 };
00256
00257
00258 template <class T>
00259 inline Vector<T> operator*(const Sparse_Matrix<T> &S, const Vector<T> &x)
00260 {
00261 return S.mult(x);
00262 }
00263
00264 template <class T>
00265 inline double norm(const Sparse_Matrix<T> &S)
00266 {
00267 return S.norm();
00268 }
00269
00270 template <class T>
00271 inline std::ostream& operator<<(std::ostream &s, const Sparse_Matrix<T> &A)
00272 {
00273 return A.print(s);
00274 }
00275
00276
00277
00278
00279 }
00280
00281 #endif