All Classes Namespaces Functions Variables 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 
66  typedef AT AtomicType;
67 
72  Matrix() : memEle(), memI(), memJ(), ele(0), I(0), J(0), m(0), n(0), k(0) { }
73 
74 // template<class Ini=All<AT> >
75 // 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_) {
76 // init(ini);
77 // }
78 // template<class Ini=All<AT> >
79 // 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_) {
80 // init(ini);
81 // }
82 
83  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_) { }
84  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); }
85 
86  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_) { }
87  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); }
88 
89  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_) { }
90  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); }
91 
92 
100  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) {
101  }
102 
106  }
107 
108  Matrix<Sparse,Ref,Ref,AT>& resize() {
109  m = n = k = 0;
110  memEle.resize(0);
111  memI.resize(0);
112  memJ.resize(0);
113  ele = 0;
114  I = J = 0;
115  return *this;
116  }
117 
118  Matrix<Sparse,Ref,Ref,AT>& resize(int n_, int k_, Noinit) {
119  m = n_; n = n_; k = k_;
120  memEle.resize(k);
121  memI.resize(m+1);
122  memJ.resize(k);
123  ele = (AT*)memEle.get();
124  I = (int*)memI.get();
125  J = (int*)memJ.get();
126  return *this;
127  }
128 
129  Matrix<Sparse,Ref,Ref,AT>& resize(int n, int k, Init ini=INIT, const AT &a=0) { return resize(m,k,Noinit()).init(a); }
130 
137  inline Matrix<Sparse,Ref,Ref,AT>& operator<<(const Matrix<Sparse,Ref,Ref,AT> &A);
138 
145  inline Matrix<Sparse,Ref,Ref,AT>& operator>>(const Matrix<Sparse,Ref,Ref,AT> &A);
146 
151  inline Matrix<Sparse,Ref,Ref,AT>& operator<<(const SquareMatrix<Ref,AT> &A);
152 
157  inline Matrix<Sparse,Ref,Ref,AT>& operator<<(const Matrix<Symmetric, Var, Var ,AT> &A);
158 
167  inline Matrix<Sparse,Ref,Ref,AT>& operator=(const Matrix<Sparse,Ref,Ref,AT> &A);
168 
173  const AT* operator()() const {
174  return ele;
175  };
176 
182  AT* operator()() {
183  return ele;
184  };
185 
186  AT& operator()(int i, int j) {
187  throw std::runtime_error("Matrix<Sparse, Ref, Ref, AT>::operator(int i, int j) is not implemented.");
188  }
189 
190  bool transposed() {
191  throw std::runtime_error("Matrix<Sparse, Ref, Ref, AT>::transposed() cannot be called.");
192  }
193 
194  int ldim() {
195  throw std::runtime_error("Matrix<Sparse, Ref, Ref, AT>::ldim() cannot be called.");
196  }
197 
202  const int* Ip() const {
203  return I;
204  };
205 
210  int* Ip() {
211  return I;
212  };
213 
218  const int* Jp() const {
219  return J;
220  };
221 
226  int* Jp() {
227  return J;
228  };
229 
234  int rows() const {return m;};
235 
240  int cols() const {return n;};
241 
247  Matrix<Sparse,Ref,Ref,AT> copy() const;
248 
256  Matrix<Sparse,Ref,Ref,AT>& init(const AT &a);
257  inline Matrix<Sparse,Ref,Ref,AT>& init(Init, const AT &a=0) { return init(a); }
258  inline Matrix<Sparse,Ref,Ref,AT>& init(Noinit, const AT &a=0) { return *this; }
259 
260  };
261  // ------------------------- Constructors -------------------------------------
262  // ----------------------------------------------------------------------------
263 
264  template <class AT>
266 
267  m=A.m;
268  n=A.n;
269  k=A.k;
270  memEle = A.memEle;
271  memI = A.memI;
272  memJ = A.memJ;
273  ele = A.ele;
274  I = A.I;
275  J = A.J;
276  k = A.k;
277 
278  return *this;
279  }
280 
281  template <class AT>
283 
284  if(!ele) {
285  m = A.m;
286  n = A.n;
287  k = A.k;
288  memEle.resize(k);
289  memI.resize(m+1);
290  memJ.resize(k);
291  ele = (AT*)memEle.get();
292  I = (int*)memI.get();
293  J = (int*)memJ.get();
294  }
295  else {
296 #ifndef FMATVEC_NO_SIZE_CHECK
297  assert(m == A.m);
298  assert(n == A.n);
299  assert(k == A.k);
300 #endif
301  }
302 
303  deepCopy(A);
304 
305  return *this;
306  }
307 
308  template <class AT>
310 
311  if(m!=A.m || n!=A.n || k!=A.k) {
312  m = A.m;
313  n = A.n;
314  k = A.k;
315  memEle.resize(k);
316  memI.resize(m+1);
317  memJ.resize(k);
318  ele = (AT*)memEle.get();
319  I = (int*)memI.get();
320  J = (int*)memJ.get();
321  }
322 
323  deepCopy(A);
324 
325  return *this;
326  }
327 
328  template <class AT>
330 
331  if(m!=A.rows() || n!=A.cols()) {
332  m = A.rows();
333  n = A.cols();
334  k = countElements(A);
335  memEle.resize(k);
336  memI.resize(m+1);
337  memJ.resize(k);
338  ele = (AT*)memEle.get();
339  I = (int*)memI.get();
340  J = (int*)memJ.get();
341  } else {
342  int k_ = countElements(A);
343  if(k != k_) {
344  k = k_;
345  memEle.resize(k);
346  memJ.resize(k);
347  ele = (AT*)memEle.get();
348  J = (int*)memJ.get();
349  }
350  }
351 
352  deepCopy(A);
353 
354  return *this;
355  }
356 
357  template <class AT>
359 
360  if (m != A.rows() || n != A.cols()) {
361  m = A.rows();
362  n = A.cols();
363  k = countElementsLT(A) * 2 - A.rows();
364  memEle.resize(k);
365  memI.resize(m + 1);
366  memJ.resize(k);
367  ele = (AT*) memEle.get();
368  I = (int*) memI.get();
369  J = (int*) memJ.get();
370  }
371  else {
372  int k_ = countElementsLT(A) * 2 - A.rows();
373  if (k != k_) {
374  k = k_;
375  memEle.resize(k);
376  memJ.resize(k);
377  ele = (AT*) memEle.get();
378  J = (int*) memJ.get();
379  }
380  }
381 
382  deepCopy(A);
383 
384  return *this;
385  }
386 
387  template <class AT>
389 
390  Matrix<Sparse,Ref,Ref,AT> A(m,NONINIT);
391 
392  A.deepCopy(*this);
393 
394  return A;
395  }
396 
397  template <class AT>
399  for(int i=0; i<=m; i++) {
400  I[i] = A.I[i];
401  }
402  for(int i=0; i<k; i++) {
403  ele[i] = A.ele[i];
404  J[i] = A.J[i];
405  }
406  }
407 
408  template <class AT> void Matrix<Sparse,Ref,Ref,AT>::deepCopy(const SquareMatrix<Ref,AT> &A) {
409  int k=0;
410  int i;
411  for(i=0; i<A.size(); i++) {
412  ele[k]=A(i,i);
413  J[k]=i;
414  I[i]=k++;
415  for(int j=0; j<i; j++) {
416  // TODO eps
417  if(fabs(A(i,j))>1e-16) {
418  ele[k]=A(i,j);
419  J[k++]=j;
420  }
421  }
422  for(int j=i+1; j<A.size(); j++) {
423  // TODO eps
424  if(fabs(A(i,j))>1e-16) {
425  ele[k]=A(i,j);
426  J[k++]=j;
427  }
428  }
429  }
430  if(n) I[i]=k;
431  }
432 
433  template <class AT> void Matrix<Sparse,Ref,Ref,AT>::deepCopy(const Matrix<Symmetric, Var, Var,AT> &A) {
434  int k=0;
435  int i;
436  for(i=0; i<A.size(); i++) {
437  ele[k]=A(i,i);
438  J[k]=i;
439  I[i]=k++;
440  for(int j=0; j<i; j++) {
441  // TODO eps
442  if(fabs(A(i,j))>1e-16) {
443  ele[k]=A(i,j);
444  J[k++]=j;
445  }
446  }
447  for(int j=i+1; j<A.size(); j++) {
448  // TODO eps
449  if(fabs(A(i,j))>1e-16) {
450  ele[k]=A(i,j);
451  J[k++]=j;
452  }
453  }
454  }
455  if(n) I[i]=k;
456  }
457 
458  template <class AT>
460  for(int i=0; i<k; i++) {
461  ele[i] = val;
462  }
463  return *this;
464  }
465 
466 /* template <class AT> template <class Type> void Matrix<Sparse,Ref,Ref,AT>::deepCopy(const Matrix<Type, AT> &A) {
467  k=0;
468  AT *eleBuf = new AT[m*n];
469  int *JBuf = new int[m*n];
470  int i;
471  for(i=0; i<A.rows(); i++) {
472  eleBuf[k]=A(i,i);
473  JBuf[k]=i;
474  I[i]=k++;
475  for(int j=0; j<i; j++) {
476  // TODO eps
477  if(fabs(A(i,j))>1e-16) {
478  eleBuf[k]=A(i,j);
479  JBuf[k++]=j;
480  }
481  }
482  for(int j=i+1; j<A.cols(); j++) {
483  // TODO eps
484  if(fabs(A(i,j))>1e-16) {
485  eleBuf[k]=A(i,j);
486  JBuf[k++]=j;
487  }
488  }
489  }
490  I[i]=k;
491  memEle.resize(k);
492  memJ.resize(k);
493  ele = (AT*)memEle.get();
494  J = (int*)memJ.get();
495  for(int i=0; i<k; i++) {
496  ele[i] = eleBuf[i];
497  J[i] = JBuf[i];
498  }
499 
500  delete [] eleBuf;
501  delete [] JBuf;
502  }
503 
504  */
505 
506  //template <> extern void Matrix<Sparse,Ref,Ref,double>::deepCopy(const Matrix<Sparse,Ref,Ref,double> &A);
507 
508 }
509 
510 #endif
Matrix< Sparse, Ref, Ref, AT > & init(const AT &a)
Initialization.
Definition: sparse_matrix.h:459
Matrix()
Standard constructor.
Definition: sparse_matrix.h:72
This is a matrix class for symmetric matrices.
Definition: var_symmetric_matrix.h:37
This is the basic matrix class for arbitrary matrices.
Definition: fmatvec.h:41
int rows() const
Number of rows.
Definition: var_symmetric_matrix.h:245
AT & operator()(int i, int j)
Standard constructor.
Definition: matrix.h:86
int rows() const
Number of rows.
Definition: general_matrix.h:272
Definition: matrix.h:39
AT * operator()()
Pointer operator.
Definition: sparse_matrix.h:182
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:210
int rows() const
Number of rows.
Definition: sparse_matrix.h:234
const int * Jp() const
Pointer operator.
Definition: sparse_matrix.h:218
int * Jp()
Pointer operator.
Definition: sparse_matrix.h:226
std::istream & operator>>(std::istream &is, Matrix< Type, Row, Col, AT > &A)
Matrix input.
Definition: matrix.h:171
const int * Ip() const
Pointer operator.
Definition: sparse_matrix.h:202
Shape class for sparse matrices.
Definition: types.h:140
Matrix(const Matrix< Sparse, Ref, Ref, AT > &A)
Copy Constructor.
Definition: sparse_matrix.h:100
Definition: matrix.h:38
~Matrix()
Destructor.
Definition: sparse_matrix.h:105
This is a matrix class of general quadratic matrices.
Definition: square_matrix.h:39
int cols() const
Number of columns.
Definition: var_symmetric_matrix.h:251
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:173
int cols() const
Number of columns.
Definition: sparse_matrix.h:240
int cols() const
Number of columns.
Definition: general_matrix.h:278

Impressum / Disclaimer / Datenschutz Generated by doxygen 1.8.5 Valid HTML