All Classes Namespaces Functions Typedefs Enumerations Pages
sparse_matrix.h
1 /* Copyright (C) 2003-2005 Martin Förg
2 
3  * This library is free software; you can redistribute it and/or
4  * modify it under the terms of the GNU Lesser General Public
5  * License as published by the Free Software Foundation; either
6  * version 2.1 of the License, or (at your option) any later version.
7  *
8  * This library is distributed in the hope that it will be useful,
9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11  * Lesser General Public License for more details.
12  *
13  * You should have received a copy of the GNU Lesser General Public
14  * License along with this library; if not, write to the Free Software
15  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
16  *
17  * Contact:
18  * martin.o.foerg@googlemail.com
19  *
20  */
21 
22 #ifndef sparse_matrix_h
23 #define sparse_matrix_h
24 
25 #include "matrix.h"
26 #include "square_matrix.h"
27 #include "types.h"
28 #include "_memory.h"
29 
30 namespace fmatvec {
31 
43  template <class AT>
44  class Matrix<Sparse,Ref,Ref,AT> {
45 
46  protected:
47 
49 
50  Memory<AT> memEle;
51  Memory<int> memI, memJ;
52  AT *ele;
53  int *I, *J;
54  int m, n, k;
55 
56  void deepCopy(const Matrix<Sparse,Ref,Ref,AT> &x);
57 
58  void deepCopy(const SquareMatrix<Ref,AT> &A);
59  void deepCopy(const Matrix<Symmetric, Var, Var,AT> &A);
60 
62 
63 
64  public:
65 
70  Matrix() : memEle(), memI(), memJ(), ele(0), I(0), J(0), m(0), n(0), k(0) { }
71 
72 // template<class Ini=All<AT> >
73 // Matrix(int n_, Ini ini=All<AT>()) : memEle(n_*n_), memI(n_+1), memJ(n_*n_), ele((AT*)memEle.get()), I((int*)memI.get()), J((int*)memJ.get()), m(n_), n(n_), k(n_*n_) {
74 // init(ini);
75 // }
76 // template<class Ini=All<AT> >
77 // Matrix(int m_, int n_, Ini ini) : memEle(n_*n_), memI(n_+1), memJ(n_*n_), ele((AT*)memEle.get()), I((int*)memI.get()), J((int*)memJ.get()), m(n_), n(n_), k(n_*n_) {
78 // init(ini);
79 // }
80 
81  Matrix(int m_, int n_, int k_, Noinit) : memEle(k_), memI(n_+1), memJ(k_), ele((AT*)memEle.get()), I((int*)memI.get()), J((int*)memJ.get()), m(n_), n(n_), k(k_) { }
82  Matrix(int m_, int n_, int k_, Init ini=INIT, const AT &a=0) : memEle(k_), memI(n_+1), memJ(k_), ele((AT*)memEle.get()), I((int*)memI.get()), J((int*)memJ.get()), m(n_), n(n_), k(k_) { init(a); }
83 
84  Matrix(int n_, Noinit) : memEle(n_*n_), memI(n_+1), memJ(n_*n_), ele((AT*)memEle.get()), I((int*)memI.get()), J((int*)memJ.get()), m(n_), n(n_), k(n_*n_) { }
85  Matrix(int n_, Init ini=INIT, const AT &a=0) : memEle(n_*n_), memI(n_+1), memJ(n_*n_), ele((AT*)memEle.get()), I((int*)memI.get()), J((int*)memJ.get()), m(n_), n(n_), k(n_*n_) { init(a); }
86 
87  Matrix(int m_, int n_, Noinit) : memEle(n_*n_), memI(n_+1), memJ(n_*n_), ele((AT*)memEle.get()), I((int*)memI.get()), J((int*)memJ.get()), m(n_), n(n_), k(n_*n_) { }
88  Matrix(int m_, int n_, Init ini=INIT, const AT &a=0) : memEle(n_*n_), memI(n_+1), memJ(n_*n_), ele((AT*)memEle.get()), I((int*)memI.get()), J((int*)memJ.get()), m(n_), n(n_), k(n_*n_) { init(a); }
89 
90 
98  Matrix(const Matrix<Sparse,Ref,Ref,AT> &A) : memEle(A.memEle), memI(A.memI), memJ(A.memJ), ele(A.ele), I(A.I), J(A.J), m(A.m), n(A.n), k(A.k) {
99  }
100 
104  }
105 
106  Matrix<Sparse,Ref,Ref,AT>& resize() {
107  m = n = k = 0;
108  memEle.resize(0);
109  memI.resize(0);
110  memJ.resize(0);
111  ele = 0;
112  I = J = 0;
113  return *this;
114  }
115 
116  Matrix<Sparse,Ref,Ref,AT>& resize(int n_, int k_, Noinit) {
117  m = n_; n = n_; k = k_;
118  memEle.resize(k);
119  memI.resize(m+1);
120  memJ.resize(k);
121  ele = (AT*)memEle.get();
122  I = (int*)memI.get();
123  J = (int*)memJ.get();
124  return *this;
125  }
126 
127  Matrix<Sparse,Ref,Ref,AT>& resize(int n, int k, Init ini=INIT, const AT &a=0) { return resize(m,k,Noinit()).init(a); }
128 
135  inline Matrix<Sparse,Ref,Ref,AT>& operator<<(const Matrix<Sparse,Ref,Ref,AT> &A);
136 
143  inline Matrix<Sparse,Ref,Ref,AT>& operator>>(const Matrix<Sparse,Ref,Ref,AT> &A);
144 
149  inline Matrix<Sparse,Ref,Ref,AT>& operator<<(const SquareMatrix<Ref,AT> &A);
150 
155  inline Matrix<Sparse,Ref,Ref,AT>& operator<<(const Matrix<Symmetric, Var, Var ,AT> &A);
156 
165  inline Matrix<Sparse,Ref,Ref,AT>& operator=(const Matrix<Sparse,Ref,Ref,AT> &A);
166 
171  const AT* operator()() const {
172  return ele;
173  };
174 
180  AT* operator()() {
181  return ele;
182  };
183 
188  const int* Ip() const {
189  return I;
190  };
191 
196  int* Ip() {
197  return I;
198  };
199 
204  const int* Jp() const {
205  return J;
206  };
207 
212  int* Jp() {
213  return J;
214  };
215 
220  int rows() const {return m;};
221 
226  int cols() const {return n;};
227 
233  Matrix<Sparse,Ref,Ref,AT> copy() const;
234 
242  Matrix<Sparse,Ref,Ref,AT>& init(const AT &a);
243  inline Matrix<Sparse,Ref,Ref,AT>& init(Init, const AT &a=0) { return init(a); }
244  inline Matrix<Sparse,Ref,Ref,AT>& init(Noinit, const AT &a=0) { return *this; }
245 
246  };
247  // ------------------------- Constructors -------------------------------------
248  // ----------------------------------------------------------------------------
249 
250  template <class AT>
252 
253  m=A.m;
254  n=A.n;
255  k=A.k;
256  memEle = A.memEle;
257  memI = A.memI;
258  memJ = A.memJ;
259  ele = A.ele;
260  I = A.I;
261  J = A.J;
262  k = A.k;
263 
264  return *this;
265  }
266 
267  template <class AT>
269 
270  if(!ele) {
271  m = A.m;
272  n = A.n;
273  k = A.k;
274  memEle.resize(k);
275  memI.resize(m+1);
276  memJ.resize(k);
277  ele = (AT*)memEle.get();
278  I = (int*)memI.get();
279  J = (int*)memJ.get();
280  }
281  else {
282 #ifndef FMATVEC_NO_SIZE_CHECK
283  assert(m == A.m);
284  assert(n == A.n);
285  assert(k == A.k);
286 #endif
287  }
288 
289  deepCopy(A);
290 
291  return *this;
292  }
293 
294  template <class AT>
296 
297  if(m!=A.m || n!=A.n || k!=A.k) {
298  m = A.m;
299  n = A.n;
300  k = A.k;
301  memEle.resize(k);
302  memI.resize(m+1);
303  memJ.resize(k);
304  ele = (AT*)memEle.get();
305  I = (int*)memI.get();
306  J = (int*)memJ.get();
307  }
308 
309  deepCopy(A);
310 
311  return *this;
312  }
313 
314  template <class AT>
316 
317  if(m!=A.rows() || n!=A.cols()) {
318  m = A.rows();
319  n = A.cols();
320  k = countElements(A);
321  memEle.resize(k);
322  memI.resize(m+1);
323  memJ.resize(k);
324  ele = (AT*)memEle.get();
325  I = (int*)memI.get();
326  J = (int*)memJ.get();
327  } else {
328  int k_ = countElements(A);
329  if(k != k_) {
330  k = k_;
331  memEle.resize(k);
332  memJ.resize(k);
333  ele = (AT*)memEle.get();
334  J = (int*)memJ.get();
335  }
336  }
337 
338  deepCopy(A);
339 
340  return *this;
341  }
342 
343  template <class AT>
345 
346  if (m != A.rows() || n != A.cols()) {
347  m = A.rows();
348  n = A.cols();
349  k = countElementsLT(A) * 2 - A.rows();
350  memEle.resize(k);
351  memI.resize(m + 1);
352  memJ.resize(k);
353  ele = (AT*) memEle.get();
354  I = (int*) memI.get();
355  J = (int*) memJ.get();
356  }
357  else {
358  int k_ = countElementsLT(A) * 2 - A.rows();
359  if (k != k_) {
360  k = k_;
361  memEle.resize(k);
362  memJ.resize(k);
363  ele = (AT*) memEle.get();
364  J = (int*) memJ.get();
365  }
366  }
367 
368  deepCopy(A);
369 
370  return *this;
371  }
372 
373  template <class AT>
375 
376  Matrix<Sparse,Ref,Ref,AT> A(m,NONINIT);
377 
378  A.deepCopy(*this);
379 
380  return A;
381  }
382 
383  template <class AT>
385  for(int i=0; i<=m; i++) {
386  I[i] = A.I[i];
387  }
388  for(int i=0; i<k; i++) {
389  ele[i] = A.ele[i];
390  J[i] = A.J[i];
391  }
392  }
393 
394  template <class AT> void Matrix<Sparse,Ref,Ref,AT>::deepCopy(const SquareMatrix<Ref,AT> &A) {
395  int k=0;
396  int i;
397  for(i=0; i<A.size(); i++) {
398  ele[k]=A(i,i);
399  J[k]=i;
400  I[i]=k++;
401  for(int j=0; j<i; j++) {
402  // TODO eps
403  if(fabs(A(i,j))>1e-16) {
404  ele[k]=A(i,j);
405  J[k++]=j;
406  }
407  }
408  for(int j=i+1; j<A.size(); j++) {
409  // TODO eps
410  if(fabs(A(i,j))>1e-16) {
411  ele[k]=A(i,j);
412  J[k++]=j;
413  }
414  }
415  }
416  if(n) I[i]=k;
417  }
418 
419  template <class AT> void Matrix<Sparse,Ref,Ref,AT>::deepCopy(const Matrix<Symmetric, Var, Var,AT> &A) {
420  int k=0;
421  int i;
422  for(i=0; i<A.size(); i++) {
423  ele[k]=A(i,i);
424  J[k]=i;
425  I[i]=k++;
426  for(int j=0; j<i; j++) {
427  // TODO eps
428  if(fabs(A(i,j))>1e-16) {
429  ele[k]=A(i,j);
430  J[k++]=j;
431  }
432  }
433  for(int j=i+1; j<A.size(); j++) {
434  // TODO eps
435  if(fabs(A(i,j))>1e-16) {
436  ele[k]=A(i,j);
437  J[k++]=j;
438  }
439  }
440  }
441  if(n) I[i]=k;
442  }
443 
444  template <class AT>
446  for(int i=0; i<k; i++) {
447  ele[i] = val;
448  }
449  return *this;
450  }
451 
452 /* template <class AT> template <class Type> void Matrix<Sparse,Ref,Ref,AT>::deepCopy(const Matrix<Type, AT> &A) {
453  k=0;
454  AT *eleBuf = new AT[m*n];
455  int *JBuf = new int[m*n];
456  int i;
457  for(i=0; i<A.rows(); i++) {
458  eleBuf[k]=A(i,i);
459  JBuf[k]=i;
460  I[i]=k++;
461  for(int j=0; j<i; j++) {
462  // TODO eps
463  if(fabs(A(i,j))>1e-16) {
464  eleBuf[k]=A(i,j);
465  JBuf[k++]=j;
466  }
467  }
468  for(int j=i+1; j<A.cols(); j++) {
469  // TODO eps
470  if(fabs(A(i,j))>1e-16) {
471  eleBuf[k]=A(i,j);
472  JBuf[k++]=j;
473  }
474  }
475  }
476  I[i]=k;
477  memEle.resize(k);
478  memJ.resize(k);
479  ele = (AT*)memEle.get();
480  J = (int*)memJ.get();
481  for(int i=0; i<k; i++) {
482  ele[i] = eleBuf[i];
483  J[i] = JBuf[i];
484  }
485 
486  delete [] eleBuf;
487  delete [] JBuf;
488  }
489 
490  */
491 
492  //template <> extern void Matrix<Sparse,Ref,Ref,double>::deepCopy(const Matrix<Sparse,Ref,Ref,double> &A);
493 
494 }
495 
496 #endif
Matrix< Sparse, Ref, Ref, AT > & init(const AT &a)
Initialization.
Definition: sparse_matrix.h:445
Matrix()
Standard constructor.
Definition: sparse_matrix.h:70
This is a matrix class for symmetric matrices.
Definition: var_symmetric_matrix.h:37
This is the basic matrix class for arbitrary matrices.
Definition: matrix.h:56
int rows() const
Number of rows.
Definition: var_symmetric_matrix.h:231
int rows() const
Number of rows.
Definition: general_matrix.h:270
Definition: matrix.h:39
AT * operator()()
Pointer operator.
Definition: sparse_matrix.h:180
int countElements(const SquareMatrix< Ref, AT > &A)
Count nonzero elements.
Definition: linear_algebra.h:1839
This is a matrix class for sparse quadratic matrices.
Definition: sparse_matrix.h:44
Definition: types.h:62
int * Ip()
Pointer operator.
Definition: sparse_matrix.h:196
int rows() const
Number of rows.
Definition: sparse_matrix.h:220
const int * Jp() const
Pointer operator.
Definition: sparse_matrix.h:204
int * Jp()
Pointer operator.
Definition: sparse_matrix.h:212
std::istream & operator>>(std::istream &is, Matrix< Type, Row, Col, AT > &A)
Matrix input.
Definition: matrix.h:170
const int * Ip() const
Pointer operator.
Definition: sparse_matrix.h:188
Shape class for sparse matrices.
Definition: types.h:140
Matrix(const Matrix< Sparse, Ref, Ref, AT > &A)
Copy Constructor.
Definition: sparse_matrix.h:98
Definition: matrix.h:38
~Matrix()
Destructor.
Definition: sparse_matrix.h:103
This is a matrix class of general quadratic matrices.
Definition: square_matrix.h:38
int cols() const
Number of columns.
Definition: var_symmetric_matrix.h:237
int countElementsLT(const Matrix< Symmetric, Ref, Ref, AT > &A)
Count nonzero elements of the low triangular part of a symmetric matrix.
Definition: linear_algebra.h:1857
const AT * operator()() const
Pointer operator.
Definition: sparse_matrix.h:171
int cols() const
Number of columns.
Definition: sparse_matrix.h:226
int cols() const
Number of columns.
Definition: general_matrix.h:276

Impressum / Disclaimer / Datenschutz Generated by doxygen 1.8.5 Valid HTML