HeterogeneousMatrixCollection.cxx
1 // Copyright (C) 2010 INRIA
2 // Author(s): Marc Fragu
3 //
4 // This file is part of the linear-algebra library Seldon,
5 // http://seldon.sourceforge.net/.
6 //
7 // Seldon is free software; you can redistribute it and/or modify it under the
8 // terms of the GNU Lesser General Public License as published by the Free
9 // Software Foundation; either version 2.1 of the License, or (at your option)
10 // any later version.
11 //
12 // Seldon is distributed in the hope that it will be useful, but WITHOUT ANY
13 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15 // more details.
16 //
17 // You should have received a copy of the GNU Lesser General Public License
18 // along with Seldon. If not, see http://www.gnu.org/licenses/.
19 
20 
21 #ifndef SELDON_FILE_HETEROGENEOUS_MATRIX_COLLECTION_CXX
22 
23 #include "HeterogeneousMatrixCollection.hxx"
24 
25 namespace Seldon
26 {
27 
28 
30  // HETEROGENEOUSMATRIXCOLLECTION //
32 
33 
34  /****************
35  * CONSTRUCTORS *
36  ****************/
37 
38 
40 
43  template <class Prop0, class Storage0,
44  class Prop1, class Storage1,
45  template <class U> class Allocator>
46  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
47  ::HeterogeneousMatrixCollection():
48  Matrix_Base<double, Allocator<double> >(), Mlocal_(), Mlocal_sum_(1),
49  Nlocal_(), Nlocal_sum_(1), collection_(), float_dense_c_(),
50  float_sparse_c_(), double_dense_c_(), double_sparse_c_()
51  {
52  nz_ = 0;
53  Mmatrix_ = 0;
54  Nmatrix_ = 0;
55  Mlocal_sum_.Fill(0);
56  Nlocal_sum_.Fill(0);
57  collection_.Fill(-1);
58  }
59 
60 
62 
66  template <class Prop0, class Storage0,
67  class Prop1, class Storage1,
68  template <class U> class Allocator>
69  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
70  ::HeterogeneousMatrixCollection(int i, int j):
71  Matrix_Base<double, Allocator<double> >(i, j),
72  Mlocal_(i), Mlocal_sum_(i + 1),
73  Nlocal_(j), Nlocal_sum_(j + 1), collection_(i, j), float_dense_c_(i, j),
74  float_sparse_c_(i, j), double_dense_c_(i, j), double_sparse_c_(i, j)
75  {
76  nz_ = 0;
77  Mmatrix_ = i;
78  Nmatrix_ = j;
79  Mlocal_.Fill(0);
80  Nlocal_.Fill(0);
81  Mlocal_sum_.Fill(0);
82  Nlocal_sum_.Fill(0);
83  collection_.Fill(-1);
84  }
85 
86 
88  template <class Prop0, class Storage0,
89  class Prop1, class Storage1,
90  template <class U> class Allocator> void
91  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
92  ::Clear()
93  {
94  float_dense_c_.Nullify();
95  float_sparse_c_.Nullify();
96  double_dense_c_.Nullify();
97  double_sparse_c_.Nullify();
98 
99  nz_ = 0;
100  Mmatrix_ = 0;
101  Nmatrix_ = 0;
102  Mlocal_.Clear();
103  Nlocal_.Clear();
104  Mlocal_sum_.Clear();
105  Nlocal_sum_.Clear();
106  collection_.Clear();
107  }
108 
109 
111  template <class Prop0, class Storage0,
112  class Prop1, class Storage1,
113  template <class U> class Allocator> void
114  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
115  ::Nullify()
116  {
117  float_dense_c_.Nullify();
118  float_sparse_c_.Nullify();
119  double_dense_c_.Nullify();
120  double_sparse_c_.Nullify();
121 
122  nz_ = 0;
123  Mmatrix_ = 0;
124  Nmatrix_ = 0;
125  Mlocal_.Clear();
126  Nlocal_.Clear();
127  Mlocal_sum_.Clear();
128  Nlocal_sum_.Clear();
129  collection_.Clear();
130  }
131 
132 
134 
138  template <class Prop0, class Storage0,
139  class Prop1, class Storage1,
140  template <class U> class Allocator> void
141  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
142  ::Nullify(int i, int j)
143  {
144 #ifdef SELDON_CHECK_BOUNDS
145  if (i < 0 || i >= Mmatrix_)
146  throw WrongRow("HeterogeneousMatrixCollection::Nullify()",
147  string("Index should be in [0, ")
148  + to_str(Mmatrix_ - 1) + "], but is equal to "
149  + to_str(i) + ".");
150  if (j < 0 || j >= Nmatrix_)
151  throw WrongCol("HeterogeneousMatrixCollection::Nullify()",
152  string("Index should be in [0, ")
153  + to_str(Nmatrix_ - 1) + "], but is equal to "
154  + to_str(j) + ".");
155 #endif
156 
157  switch (collection_(i, j))
158  {
159  case 0:
160  nz_ -= float_dense_c_.GetMatrix(i, j).GetDataSize();
161  float_dense_c_.Nullify(i, j);
162  case 1:
163  nz_ -= float_sparse_c_.GetMatrix(i, j).GetDataSize();
164  float_sparse_c_.Nullify(i, j);
165  break;
166  case 2:
167  nz_ -= double_dense_c_.GetMatrix(i, j).GetDataSize();
168  double_dense_c_.Nullify(i, j);
169  break;
170  case 3:
171  nz_ -= double_sparse_c_.GetMatrix(i, j).GetDataSize();
172  double_sparse_c_.Nullify(i, j);
173  break;
174  }
175 
176  collection_(i, j) = -1;
177  }
178 
179 
181  template <class Prop0, class Storage0,
182  class Prop1, class Storage1,
183  template <class U> class Allocator> void
184  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
185  ::Deallocate()
186  {
187  float_dense_c_.Deallocate();
188  float_sparse_c_.Deallocate();
189  double_dense_c_.Deallocate();
190  double_sparse_c_.Deallocate();
191  this->Clear();
192  }
193 
194 
195  /*********************
196  * MEMORY MANAGEMENT *
197  *********************/
198 
199 
201 
207  template <class Prop0, class Storage0,
208  class Prop1, class Storage1,
209  template <class U> class Allocator> void
210  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
211  ::Reallocate(int i, int j)
212  {
213  nz_ = 0;
214  Mmatrix_ = i;
215  Nmatrix_ = j;
216  Mlocal_.Reallocate(i);
217  Nlocal_.Reallocate(j);
218  Mlocal_sum_.Reallocate(i + 1);
219  Nlocal_sum_.Reallocate(j + 1);
220  Mlocal_.Fill(0);
221  Nlocal_.Fill(0);
222  Mlocal_sum_.Fill(0);
223  Nlocal_sum_.Fill(0);
224 
225  collection_.Reallocate(i, j);
226  float_dense_c_.Reallocate(i, j);
227  float_sparse_c_.Reallocate(i, j);
228  double_dense_c_.Reallocate(i, j);
229  double_sparse_c_.Reallocate(i, j);
230  }
231 
232 
234 
239  template <class Prop0, class Storage0,
240  class Prop1, class Storage1,
241  template <class U> class Allocator> void
242  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
243  ::SetMatrix(int i, int j, const typename HeterogeneousMatrixCollection<
244  Prop0, Storage0, Prop1, Storage1, Allocator>::float_dense_m& A)
245  {
246 #ifdef SELDON_CHECK_BOUNDS
247  if (i < 0 || i >= Mmatrix_)
248  throw WrongRow("HeterogeneousMatrixCollection::"
249  "SetMatrix(float_dense_m)",
250  string("Index should be in [0, ")
251  + to_str(Mmatrix_ - 1) + "], but is equal to "
252  + to_str(i) + ".");
253  if (j < 0 || j >= Nmatrix_)
254  throw WrongCol("HeterogeneousMatrixCollection::"
255  "SetMatrix(float_dense_m)",
256  string("Index should be in [0, ")
257  + to_str(Nmatrix_ - 1) + "], but is equal to "
258  + to_str(j) + ".");
259  if ((Mlocal_(i) != 0) && (Mlocal_(i) != A.GetM()))
260  throw WrongDim("HeterogeneousMatrixCollection::"
261  "SetMatrix(float_dense_m)",
262  string("The matrix expected should have ")
263  + to_str(this->Mlocal_(i)) + " lines, but has "
264  + to_str(A.GetM()) + " lines.");
265  if ((Nlocal_(j) != 0) && (Nlocal_(j) != A.GetN()))
266  throw WrongDim("HeterogeneousMatrixCollection::"
267  "SetMatrix(float_dense_m)",
268  string("The matrix expected should have ")
269  + to_str(this->Nlocal_(j)) + " columns, but has "
270  + to_str(A.GetN()) + " columns.");
271 #endif
272 
273  Nullify(i, j);
274 
275  collection_(i, j) = 0;
276 
277  int Mdiff = A.GetM() - Mlocal_(i);
278  int Ndiff = A.GetN() - Nlocal_(j);
279 
280  Mlocal_(i) = A.GetM();
281  Nlocal_(j) = A.GetN();
282 
283  for (int k = i + 1; k < Mmatrix_ + 1; k++)
284  Mlocal_sum_(k) += Mdiff;
285 
286  for (int k = j + 1; k < Nmatrix_ + 1; k++)
287  Nlocal_sum_(k) += Ndiff;
288 
289  this->m_ = Mlocal_sum_(Mmatrix_);
290  this->n_ = Nlocal_sum_(Nmatrix_);
291 
292  float_dense_c_.SetMatrix(i, j, A);
293  }
294 
295 
297 
302  template <class Prop0, class Storage0,
303  class Prop1, class Storage1,
304  template <class U> class Allocator> void
305  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
306  ::SetMatrix(int i, int j, const typename HeterogeneousMatrixCollection<
307  Prop0, Storage0, Prop1, Storage1, Allocator>::float_sparse_m& A)
308  {
309 #ifdef SELDON_CHECK_BOUNDS
310  if (i < 0 || i >= Mmatrix_)
311  throw WrongRow("HeterogeneousMatrixCollection::"
312  "SetMatrix(float_sparse_m)",
313  string("Index should be in [0, ")
314  + to_str(Mmatrix_ - 1) + "], but is equal to "
315  + to_str(i) + ".");
316  if (j < 0 || j >= Nmatrix_)
317  throw WrongCol("HeterogeneousMatrixCollection::"
318  "SetMatrix(float_sparse_m)",
319  string("Index should be in [0, ")
320  + to_str(Nmatrix_ - 1) + "], but is equal to "
321  + to_str(j) + ".");
322  if ((Mlocal_(i) != 0) && (Mlocal_(i) != A.GetM()))
323  throw WrongDim("HeterogeneousMatrixCollection::"
324  "SetMatrix(float_sparse_m)",
325  string("The matrix expected should have ")
326  + to_str(this->Mlocal_(i)) + " lines, but has "
327  + to_str(A.GetM()) + " lines.");
328  if ((Nlocal_(j) != 0) && (Nlocal_(j) != A.GetN()))
329  throw WrongDim("HeterogeneousMatrixCollection::"
330  "SetMatrix(float_sparse_m)",
331  string("The matrix expected should have ")
332  + to_str(this->Nlocal_(j)) + " columns, but has "
333  + to_str(A.GetN()) + " columns.");
334 #endif
335 
336  Nullify(i, j);
337 
338  collection_(i, j) = 1;
339 
340  int Mdiff = A.GetM() - Mlocal_(i);
341  int Ndiff = A.GetN() - Nlocal_(j);
342 
343  Mlocal_(i) = A.GetM();
344  Nlocal_(j) = A.GetN();
345 
346  for (int k = i + 1; k < Mmatrix_ + 1; k++)
347  Mlocal_sum_(k) += Mdiff;
348 
349  for (int k = j + 1; k < Nmatrix_ + 1; k++)
350  Nlocal_sum_(k) += Ndiff;
351 
352  this->m_ = Mlocal_sum_(Mmatrix_);
353  this->n_ = Nlocal_sum_(Nmatrix_);
354 
355  float_sparse_c_.SetMatrix(i, j, A);
356  }
357 
358 
360 
365  template <class Prop0, class Storage0,
366  class Prop1, class Storage1,
367  template <class U> class Allocator> void
368  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
369  ::SetMatrix(int i, int j, const typename HeterogeneousMatrixCollection<
370  Prop0, Storage0, Prop1, Storage1, Allocator>::double_dense_m& A)
371  {
372 #ifdef SELDON_CHECK_BOUNDS
373  if (i < 0 || i >= Mmatrix_)
374  throw WrongRow("HeterogeneousMatrixCollection::"
375  "SetMatrix(double_dense_m)",
376  string("Index should be in [0, ")
377  + to_str(Mmatrix_ - 1) + "], but is equal to "
378  + to_str(i) + ".");
379  if (j < 0 || j >= Nmatrix_)
380  throw WrongCol("HeterogeneousMatrixCollection::"
381  "SetMatrix(double_dense_m)",
382  string("Index should be in [0, ")
383  + to_str(Nmatrix_ - 1) + "], but is equal to "
384  + to_str(j) + ".");
385  if ((Mlocal_(i) != 0) && (Mlocal_(i) != A.GetM()))
386  throw WrongDim("HeterogeneousMatrixCollection::"
387  "SetMatrix(double_dense_m)",
388  string("The matrix expected should have ")
389  + to_str(this->Mlocal_(i)) + " lines, but has "
390  + to_str(A.GetM()) + " lines.");
391  if ((Nlocal_(j) != 0) && (Nlocal_(j) != A.GetN()))
392  throw WrongDim("HeterogeneousMatrixCollection::"
393  "SetMatrix(double_dense_m)",
394  string("The matrix expected should have ")
395  + to_str(this->Nlocal_(j)) + " columns, but has "
396  + to_str(A.GetN()) + " columns.");
397 #endif
398 
399  Nullify(i, j);
400 
401  collection_(i, j) = 2;
402 
403  int Mdiff = A.GetM() - Mlocal_(i);
404  int Ndiff = A.GetN() - Nlocal_(j);
405 
406  Mlocal_(i) = A.GetM();
407  Nlocal_(j) = A.GetN();
408 
409  for (int k = i + 1; k < Mmatrix_ + 1; k++)
410  Mlocal_sum_(k) += Mdiff;
411 
412  for (int k = j + 1; k < Nmatrix_ + 1; k++)
413  Nlocal_sum_(k) += Ndiff;
414 
415  this->m_ = Mlocal_sum_(Mmatrix_);
416  this->n_ = Nlocal_sum_(Nmatrix_);
417 
418  double_dense_c_.SetMatrix(i, j, A);
419  }
420 
421 
423 
428  template <class Prop0, class Storage0,
429  class Prop1, class Storage1,
430  template <class U> class Allocator> void
431  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
432  ::SetMatrix(int i, int j,
433  const typename HeterogeneousMatrixCollection< Prop0, Storage0,
434  Prop1, Storage1, Allocator>::double_sparse_m& A)
435  {
436 #ifdef SELDON_CHECK_BOUNDS
437  if (i < 0 || i >= Mmatrix_)
438  throw WrongRow("HeterogeneousMatrixCollection::"
439  "SetMatrix(double_sparse_m)",
440  string("Index should be in [0, ")
441  + to_str(Mmatrix_ - 1) + "], but is equal to "
442  + to_str(i) + ".");
443  if (j < 0 || j >= Nmatrix_)
444  throw WrongCol("HeterogeneousMatrixCollection::"
445  "SetMatrix(double_sparse_m)",
446  string("Index should be in [0, ")
447  + to_str(Nmatrix_ - 1) + "], but is equal to "
448  + to_str(j) + ".");
449  if ((Mlocal_(i) != 0) && (Mlocal_(i) != A.GetM()))
450  throw WrongDim("HeterogeneousMatrixCollection::"
451  "SetMatrix(double_sparse_m)",
452  string("The matrix expected should have ")
453  + to_str(this->Mlocal_(i)) + " lines, but has "
454  + to_str(A.GetM()) + " lines.");
455  if ((Nlocal_(j) != 0) && (Nlocal_(j) != A.GetN()))
456  throw WrongDim("HeterogeneousMatrixCollection::"
457  "SetMatrix(double_sparse_m)",
458  string("The matrix expected should have ")
459  + to_str(this->Nlocal_(j)) + " columns, but has "
460  + to_str(A.GetN()) + " columns.");
461 #endif
462 
463  Nullify(i, j);
464 
465  collection_(i, j) = 3;
466 
467  int Mdiff = A.GetM() - Mlocal_(i);
468  int Ndiff = A.GetN() - Nlocal_(j);
469 
470  Mlocal_(i) = A.GetM();
471  Nlocal_(j) = A.GetN();
472 
473  for (int k = i + 1; k < Mmatrix_ + 1; k++)
474  Mlocal_sum_(k) += Mdiff;
475 
476  for (int k = j + 1; k < Nmatrix_ + 1; k++)
477  Nlocal_sum_(k) += Ndiff;
478 
479  this->m_ = Mlocal_sum_(Mmatrix_);
480  this->n_ = Nlocal_sum_(Nmatrix_);
481 
482  double_sparse_c_.SetMatrix(i, j, A);
483  }
484 
485 
486  /**********************************
487  * ELEMENT ACCESS AND AFFECTATION *
488  **********************************/
489 
490 
492 
498  template <class Prop0, class Storage0,
499  class Prop1, class Storage1,
500  template <class U> class Allocator> void
501  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
502  ::GetMatrix(int i, int j, typename HeterogeneousMatrixCollection<
503  Prop0, Storage0, Prop1, Storage1, Allocator>::float_dense_m& M)
504  const
505  {
506 #ifdef SELDON_CHECK_BOUNDS
507  if (i < 0 || i >= Mmatrix_)
508  throw WrongRow("HeterogeneousMatrixCollection::"
509  "GetMatrix(float_dense_m)",
510  string("Row index should be in [0, ")
511  + to_str(Mmatrix_ - 1) + "], but is equal to "
512  + to_str(i) + ".");
513  if (j < 0 || j >= Nmatrix_)
514  throw WrongCol("HeterogeneousMatrixCollection::"
515  "GetMatrix(float_dense_m)",
516  string("Column index should be in [0, ")
517  + to_str(Nmatrix_ - 1) + "], but is equal to "
518  + to_str(j) + ".");
519 #endif
520 
521  if (collection_(i, j) != 0)
522  {
523  string matrix_type;
524  switch(collection_(i, j))
525  {
526  case 1:
527  matrix_type = "float_sparse";
528  break;
529  case 2:
530  matrix_type = "double_dense";
531  break;
532  case 3:
533  matrix_type = "double_sparse";
534  break;
535  default:
536  throw
537  WrongArgument("HeterogeneousMatrixCollection::GetMatrix(i, j,"
538  "Matrix<Float, Dense> M)",
539  "Underlying matrix (" + to_str(i) + " ,"
540  + to_str(j) + " ) not defined.");
541  }
542 
543  throw WrongArgument("HeterogeneousMatrixCollection::GetMatrix(i, j, "
544  "Matrix<Float, Dense> M)",
545  string("Wrong type for matrix ")
546  + matrix_type + " M.");
547  }
548 
549  M.SetData(Mlocal_(i), Nlocal_(j),
550  float_dense_c_.GetMatrix(i, j).GetData());
551  }
552 
553 
555 
561  template <class Prop0, class Storage0,
562  class Prop1, class Storage1,
563  template <class U> class Allocator> void
564  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
565  ::GetMatrix(int i, int j, typename HeterogeneousMatrixCollection<
566  Prop0, Storage0, Prop1, Storage1, Allocator>::float_sparse_m& M)
567  const
568  {
569 #ifdef SELDON_CHECK_BOUNDS
570  if (i < 0 || i >= Mmatrix_)
571  throw WrongRow("HeterogeneousMatrixCollection::"
572  "GetMatrix(float_sparse_m)",
573  string("Row index should be in [0, ")
574  + to_str(Mmatrix_ - 1) + "], but is equal to "
575  + to_str(i) + ".");
576  if (j < 0 || j >= Nmatrix_)
577  throw WrongCol("HeterogeneousMatrixCollection::"
578  "GetMatrix(float_sparse_m)",
579  string("Column index should be in [0, ")
580  + to_str(Nmatrix_ - 1) + "], but is equal to "
581  + to_str(j) + ".");
582 #endif
583 
584  if (collection_(i, j) != 1)
585  {
586  string matrix_type;
587  switch(collection_(i, j))
588  {
589  case 0:
590  matrix_type = "float_dense";
591  break;
592  case 2:
593  matrix_type = "double_dense";
594  break;
595  case 3:
596  matrix_type = "double_sparse";
597  break;
598  default:
599  throw
600  WrongArgument("HeterogeneousMatrixCollection::GetMatrix(i, j, "
601  "Matrix<Float, Sparse> M)",
602  "Underlying matrix (" + to_str(i) + " ,"
603  + to_str(j) + " ) not defined.");
604  }
605 
606  throw WrongArgument("HeterogeneousMatrixCollection::GetMatrix(i, j, "
607  "Matrix<Float, Sparse> M)",
608  string("Wrong type for matrix ")
609  + matrix_type + " M.");
610  }
611 
612  M.SetData(float_sparse_c_.GetMatrix(i, j).GetM(),
613  float_sparse_c_.GetMatrix(i, j).GetN(),
614  float_sparse_c_.GetMatrix(i, j).GetNonZeros(),
615  float_sparse_c_.GetMatrix(i, j).GetData(),
616  float_sparse_c_.GetMatrix(i, j).GetPtr(),
617  float_sparse_c_.GetMatrix(i, j).GetInd());
618  }
619 
620 
622 
628  template <class Prop0, class Storage0,
629  class Prop1, class Storage1,
630  template <class U> class Allocator> void
631  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
632  ::GetMatrix(int i, int j, typename HeterogeneousMatrixCollection<
633  Prop0, Storage0, Prop1, Storage1, Allocator>::double_dense_m& M)
634  const
635  {
636 #ifdef SELDON_CHECK_BOUNDS
637  if (i < 0 || i >= Mmatrix_)
638  throw WrongRow("HeterogeneousMatrixCollection::"
639  "GetMatrix(double_dense_m)",
640  string("Row index should be in [0, ")
641  + to_str(Mmatrix_ - 1) + "], but is equal to "
642  + to_str(i) + ".");
643  if (j < 0 || j >= Nmatrix_)
644  throw WrongCol("HeterogeneousMatrixCollection::"
645  "GetMatrix(double_dense_m)",
646  string("Column index should be in [0, ")
647  + to_str(Nmatrix_ - 1) + "], but is equal to "
648  + to_str(j) + ".");
649 #endif
650 
651  if (collection_(i, j) != 2)
652  {
653  string matrix_type;
654  switch(collection_(i, j))
655  {
656  case 0:
657  matrix_type = "float_dense";
658  break;
659  case 1:
660  matrix_type = "float_sparse";
661  break;
662  case 3:
663  matrix_type = "double_sparse";
664  break;
665  default:
666  throw
667  WrongArgument("HeterogeneousMatrixCollection::GetMatrix(i, j, "
668  "Matrix<Double, Dense> M)",
669  "Underlying matrix (" + to_str(i) + " ,"
670  + to_str(j) + " ) not defined.");
671  }
672 
673  throw WrongArgument("HeterogeneousMatrixCollection::GetMatrix(i, j, "
674  "Matrix<Double, Dense> M)",
675  string("Wrong type for matrix ")
676  + matrix_type + " M.");
677  }
678 
679  M.SetData(Mlocal_(i), Nlocal_(j),
680  double_dense_c_.GetMatrix(i, j).GetData());
681  }
682 
683 
685 
691  template <class Prop0, class Storage0,
692  class Prop1, class Storage1,
693  template <class U> class Allocator> void
694  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
695  ::GetMatrix(int i, int j, typename HeterogeneousMatrixCollection<Prop0,
696  Storage0, Prop1, Storage1, Allocator>::double_sparse_m& M)
697  const
698  {
699 #ifdef SELDON_CHECK_BOUNDS
700  if (i < 0 || i >= Mmatrix_)
701  throw WrongRow("HeterogeneousMatrixCollection::"
702  "GetMatrix(double_sparse_m)",
703  string("Row index should be in [0, ")
704  + to_str(Mmatrix_ - 1) + "], but is equal to "
705  + to_str(i) + ".");
706  if (j < 0 || j >= Nmatrix_)
707  throw WrongCol("HeterogeneousMatrixCollection::"
708  "GetMatrix(double_sparse_m)",
709  string("Column index should be in [0, ")
710  + to_str(Nmatrix_ - 1) + "], but is equal to "
711  + to_str(j) + ".");
712 #endif
713 
714  if (collection_(i, j) != 3)
715  {
716  string matrix_type;
717  switch(collection_(i, j))
718  {
719  case 0:
720  matrix_type = "float_dense";
721  break;
722  case 1:
723  matrix_type = "float_sparse";
724  break;
725  case 2:
726  matrix_type = "double_dense";
727  break;
728  default:
729  throw
730  WrongArgument("HeterogeneousMatrixCollection::GetMatrix(i, j,"
731  "Matrix<Double, Sparse> M)",
732  "Underlying matrix (" + to_str(i) + " ,"
733  + to_str(j) + " ) not defined.");
734  }
735 
736  throw WrongArgument("HeterogeneousMatrixCollection::GetMatrix(i, j,"
737  "Matrix<Double, Sparse> M)",
738  string("Wrong type for matrix ")
739  + matrix_type + " M.");
740  }
741 
742  M.Nullify();
743  M.SetData(double_sparse_c_.GetMatrix(i, j).GetM(),
744  double_sparse_c_.GetMatrix(i, j).GetN(),
745  double_sparse_c_.GetMatrix(i, j).GetNonZeros(),
746  double_sparse_c_.GetMatrix(i, j).GetData(),
747  double_sparse_c_.GetMatrix(i, j).GetPtr(),
748  double_sparse_c_.GetMatrix(i, j).GetInd());
749  }
750 
751 
753 
759  template <class Prop0, class Storage0,
760  class Prop1, class Storage1,
761  template <class U> class Allocator>
762  double HeterogeneousMatrixCollection<Prop0, Storage0,
763  Prop1, Storage1, Allocator>
764  ::operator() (int i, int j) const
765  {
766 
767 #ifdef SELDON_CHECK_BOUNDS
768  if (i < 0 || i >= this->Mlocal_sum_(Mmatrix_))
769  throw WrongRow("HeterogeneousMatrixCollection::operator()",
770  string("Index should be in [0, ")
771  + to_str(this->Mlocal_sum_(Mmatrix_) - 1)
772  + "], but is equal to "
773  + to_str(i) + ".");
774  if (j < 0 || j >= this->Nlocal_sum_(Nmatrix_))
775  throw WrongCol("HeterogeneousMatrixCollection::operator()",
776  string("Index should be in [0, ")
777  + to_str(this->Nlocal_sum_(Nmatrix_) - 1)
778  + "], but is equal to "
779  + to_str(j) + ".");
780 #endif
781 
782  int i_global = 0;
783  while (i >= Mlocal_sum_(i_global))
784  i_global++;
785  i_global--;
786 
787  int j_global = 0;
788  while (j >= Nlocal_sum_(j_global))
789  j_global++;
790  j_global--;
791 
792  double res = 0.;
793  switch(collection_(i_global, j_global))
794  {
795  case 0:
796  res = double(float_dense_c_.GetMatrix(i_global, j_global)
797  (i - Mlocal_sum_(i_global), j - Nlocal_sum_(j_global)));
798  break;
799  case 1:
800  res = double(float_sparse_c_.GetMatrix(i_global, j_global)
801  (i - Mlocal_sum_(i_global), j - Nlocal_sum_(j_global)));
802  break;
803  case 2:
804  res = double_dense_c_.GetMatrix(i_global, j_global)
805  (i - Mlocal_sum_(i_global), j - Nlocal_sum_(j_global));
806  break;
807  case 3:
808  res = double_sparse_c_.GetMatrix(i_global, j_global)
809  (i - Mlocal_sum_(i_global), j - Nlocal_sum_(j_global));
810  break;
811  default:
812  throw
813  WrongArgument("HeterogeneousMatrixCollection::operator(int, int)",
814  "Underlying matrix (" + to_str(i) + " ,"
815  + to_str(j) + " ) not defined.");
816  }
817  return res;
818  }
819 
820 
822 
827  template <class Prop0, class Storage0,
828  class Prop1, class Storage1,
829  template <class U> class Allocator>
830  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>&
831  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
832  ::operator= (const HeterogeneousMatrixCollection<Prop0, Storage0,
833  Prop1, Storage1, Allocator>& A)
834  {
835  this->Copy(A);
836  return *this;
837  }
838 
839 
841 
846  template <class Prop0, class Storage0,
847  class Prop1, class Storage1,
848  template <class U> class Allocator>
849  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
850  ::Copy(const HeterogeneousMatrixCollection<Prop0, Storage0, Prop1,
851  Storage1, Allocator>& A)
852  {
853  Clear();
854  this->nz_ = A.nz_;
855  Mmatrix_ = A.Mmatrix_;
856  Nmatrix_ = A.Nmatrix_;
857  this->m_ = A.GetM();
858  this->n_ = A.GetN();
859 
860  this->Mlocal_ = A.Mlocal_;
861  this->Mlocal_sum_ = A.Mlocal_sum_;
862  this->Nlocal_ = A.Nlocal_;
863  this->Nlocal_sum_ = A.Nlocal_sum_;
864 
865  collection_.Copy(A.collection_);
866 
867  float_dense_c_.Reallocate(Mmatrix_, Nmatrix_);
868  float_sparse_c_.Reallocate(Mmatrix_, Nmatrix_);
869  double_dense_c_.Reallocate(Mmatrix_, Nmatrix_);
870  double_sparse_c_.Reallocate(Mmatrix_, Nmatrix_);
871 
872  float_dense_m m0a;
873  float_sparse_m m1a;
874  double_dense_m m2a;
875  double_sparse_m m3a;
876 
877  for (int i = 0; i < Mmatrix_; i++ )
878  for (int j = 0; j < Nmatrix_; j++)
879  {
880  switch (A.GetType(i, j))
881  {
882  case 0:
883  A.GetMatrix(i, j, m0a);
884  SetMatrix(i, j, m0a);
885  m0a.Nullify();
886  break;
887  case 1:
888  A.GetMatrix(i, j, m1a);
889  SetMatrix(i, j, m1a);
890  m1a.Nullify();
891  break;
892  case 2:
893  A.GetMatrix(i, j, m2a);
894  SetMatrix(i, j, m2a);
895  m2a.Nullify();
896  break;
897  case 3:
898  A.GetMatrix(i, j, m3a);
899  SetMatrix(i, j, m3a);
900  m3a.Nullify();
901  break;
902  default:
903  throw WrongArgument("Matrix<FloatDouble, DenseSparseCollection>"
904  "::MltAdd(alpha, A, B, beta, C) ",
905  "Underlying matrix C (" + to_str(i) + " ,"
906  + to_str(j) + " ) not defined.");
907  }
908  }
909  }
910 
911 
912  /************************
913  * CONVENIENT FUNCTIONS *
914  ************************/
915 
916 
918 
921  template <class Prop0, class Storage0,
922  class Prop1, class Storage1,
923  template <class U> class Allocator> void
924  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
925  ::Print() const
926  {
927  for (int i = 0; i < Mlocal_sum_(Mmatrix_); i++)
928  {
929  for (int j = 0; j < Nlocal_sum_(Nmatrix_); j++)
930  cout << (*this)(i, j) << endl;
931  cout << endl;
932  }
933  }
934 
935 
937 
945  template <class Prop0, class Storage0,
946  class Prop1, class Storage1,
947  template <class U> class Allocator> void
948  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
949  ::Write(string FileName, bool with_size) const
950  {
951  ofstream FileStream;
952  FileStream.open(FileName.c_str(), ofstream::binary);
953 
954 #ifdef SELDON_CHECK_IO
955  // Checks if the file was opened.
956  if (!FileStream.is_open())
957  throw IOError("Matrix_Pointers::Write(string FileName)",
958  string("Unable to open file \"") + FileName + "\".");
959 #endif
960 
961  this->Write(FileStream, with_size);
962 
963  FileStream.close();
964  }
965 
966 
968 
976  template <class Prop0, class Storage0,
977  class Prop1, class Storage1,
978  template <class U> class Allocator> void
979  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
980  ::Write(ostream& FileStream, bool with_size = true) const
981  {
982 
983 #ifdef SELDON_CHECK_IO
984  // Checks if the stream is ready.
985  if (!FileStream.good())
986  throw IOError("HeterogeneousMatrixCollection"
987  "::Write(ostream& FileStream)",
988  "The stream is not ready.");
989 #endif
990 
991  if (with_size)
992  {
993  FileStream.write(reinterpret_cast<char*>(const_cast<int*>(&Mmatrix_)),
994  sizeof(int));
995  FileStream.write(reinterpret_cast<char*>(const_cast<int*>(&Nmatrix_)),
996  sizeof(int));
997  }
998 
999  collection_.Write(FileStream, with_size);
1000 
1001  float_dense_m m0a;
1002  float_sparse_m m1a;
1003  double_dense_m m2a;
1004  double_sparse_m m3a;
1005 
1006  int i, j;
1007  for (i = 0; i < Mmatrix_; i++)
1008  for (j = 0; j < Nmatrix_; j++)
1009  {
1010  switch (GetType(i, j))
1011  {
1012  case 0:
1013  GetMatrix(i, j, m0a);
1014  m0a.Write(FileStream, with_size);
1015  m0a.Nullify();
1016  break;
1017  case 1:
1018  throw Undefined("Matrix<FloatDouble, DenseSparseCollection>"
1019  "Storage0, Prop1, Storage1, Allocator>"
1020  "::Write(ostream& FileStream, bool "
1021  "with_size = true) ");
1022  case 2:
1023  GetMatrix(i, j, m2a);
1024  m2a.Write(FileStream, with_size);
1025  m2a.Nullify();
1026  break;
1027  case 3:
1028  throw Undefined("Matrix<FloatDouble, DenseSparseCollection>"
1029  "Storage0, Prop1, Storage1, Allocator>"
1030  "::Write(ostream& FileStream, bool "
1031  "with_size = true) ");
1032  default:
1033  throw WrongArgument("Matrix<FloatDouble, DenseSparseCollection>"
1034  "::Write(ostream& FileStream, "
1035  "bool with_size = true) ",
1036  "Underlying matrix A (" + to_str(i) + " ,"
1037  + to_str(j) + " ) not defined.");
1038  }
1039  }
1040 
1041 
1042 #ifdef SELDON_CHECK_IO
1043  // Checks if data was written.
1044  if (!FileStream.good())
1045  throw IOError("HeterogeneousMatrixCollection"
1046  "::Write(ostream& FileStream)",
1047  "Output operation failed.");
1048 #endif
1049 
1050  }
1051 
1052 
1054 
1059  template <class Prop0, class Storage0,
1060  class Prop1, class Storage1,
1061  template <class U> class Allocator> void
1062  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
1063  ::WriteText(string FileName) const
1064  {
1065  ofstream FileStream;
1066  FileStream.precision(cout.precision());
1067  FileStream.flags(cout.flags());
1068  FileStream.open(FileName.c_str());
1069 
1070 #ifdef SELDON_CHECK_IO
1071  // Checks if the file was opened.
1072  if (!FileStream.is_open())
1073  throw IOError("HeterogeneousMatrixCollection"
1074  "::WriteText(string FileName)",
1075  string("Unable to open file \"") + FileName + "\".");
1076 #endif
1077 
1078  this->WriteText(FileStream);
1079 
1080  FileStream.close();
1081  }
1082 
1083 
1085 
1091  template <class Prop0, class Storage0,
1092  class Prop1, class Storage1,
1093  template <class U> class Allocator> void
1094  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
1095  ::WriteText(ostream& FileStream) const
1096  {
1097 
1098 #ifdef SELDON_CHECK_IO
1099  // Checks if the file is ready.
1100  if (!FileStream.good())
1101  throw IOError("HeterogeneousMatrixCollection"
1102  "::WriteText(ostream& FileStream)",
1103  "The stream is not ready.");
1104 #endif
1105 
1106  float_dense_m m0a;
1107  float_sparse_m m1a;
1108  double_dense_m m2a;
1109  double_sparse_m m3a;
1110 
1111  int i, j;
1112  for (i = 0; i < Mmatrix_; i++)
1113  for (j = 0; j < Nmatrix_; j++)
1114  {
1115  switch (GetType(i, j))
1116  {
1117  case 0:
1118  GetMatrix(i, j, m0a);
1119  m0a.WriteText(FileStream);
1120  m0a.Nullify();
1121  break;
1122  case 1:
1123  GetMatrix(i, j, m1a);
1124  m1a.WriteText(FileStream);
1125  m1a.Nullify();
1126  break;
1127  case 2:
1128  GetMatrix(i, j, m2a);
1129  m2a.WriteText(FileStream);
1130  m2a.Nullify();
1131  break;
1132  case 3:
1133  GetMatrix(i, j, m3a);
1134  m3a.WriteText(FileStream);
1135  m3a.Nullify();
1136  break;
1137  default:
1138  throw WrongArgument("Matrix<FloatDouble, DenseSparseCollection>"
1139  "::Write(ostream& FileStream, "
1140  "bool with_size = true) ",
1141  "Underlying matrix A (" + to_str(i) + " ,"
1142  + to_str(j) + " ) not defined.");
1143  }
1144  FileStream << endl;
1145  }
1146 
1147 #ifdef SELDON_CHECK_IO
1148  // Checks if data was written.
1149  if (!FileStream.good())
1150  throw IOError("HeterogeneousMatrixCollection"
1151  "::WriteText(ostream& FileStream)",
1152  "Output operation failed.");
1153 #endif
1154 
1155  }
1156 
1157 
1159 
1165  template <class Prop0, class Storage0,
1166  class Prop1, class Storage1,
1167  template <class U> class Allocator> void
1168  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
1169  ::Read(string FileName)
1170  {
1171  ifstream FileStream;
1172  FileStream.open(FileName.c_str(), ifstream::binary);
1173 
1174 #ifdef SELDON_CHECK_IO
1175  // Checks if the file was opened.
1176  if (!FileStream.is_open())
1177  throw IOError("HeterogeneousMatrixCollection<Prop0, Storage0, Prop1,"
1178  " Storage1, Allocator>::Read(string FileName)",
1179  string("Unable to open file \"") + FileName + "\".");
1180 #endif
1181 
1182  this->Read(FileStream);
1183 
1184  FileStream.close();
1185  }
1186 
1187 
1189 
1195  template <class Prop0, class Storage0,
1196  class Prop1, class Storage1,
1197  template <class U> class Allocator> void
1198  HeterogeneousMatrixCollection<Prop0, Storage0, Prop1, Storage1, Allocator>
1199  ::Read(istream& FileStream)
1200  {
1201 
1202 #ifdef SELDON_CHECK_IO
1203  // Checks if the stream is ready.
1204  if (!FileStream.good())
1205  throw IOError("HeterogeneousMatrixCollection<Prop0, Storage0, Prop1,"
1206  " Storage1, Allocator>::Read(istream& FileStream)",
1207  "The stream is not ready.");
1208 #endif
1209 
1210  int *new_m, *new_n;
1211  new_m = new int;
1212  new_n = new int;
1213 
1214  FileStream.read(reinterpret_cast<char*>(new_m), sizeof(int));
1215  FileStream.read(reinterpret_cast<char*>(new_n), sizeof(int));
1216 
1217  this->Reallocate(*new_m, *new_n);
1218 
1219  collection_.Read(FileStream);
1220 
1221  float_dense_m m0a;
1222  float_sparse_m m1a;
1223  double_dense_m m2a;
1224  double_sparse_m m3a;
1225  int i, j;
1226  for (i = 0; i < Mmatrix_; i++)
1227  for (j = 0; j < Nmatrix_; j++)
1228  {
1229  switch (GetType(i, j))
1230  {
1231  case 0:
1232  m0a.Read(FileStream);
1233  SetMatrix(i, j, m0a);
1234  m0a.Nullify();
1235  break;
1236  case 1:
1237  throw Undefined("Matrix<FloatDouble, DenseSparseCollection>"
1238  "Storage0, Prop1, Storage1, Allocator>"
1239  "::Read(istream& FileStream)");
1240  case 2:
1241  m2a.Read(FileStream);
1242  SetMatrix(i, j, m2a);
1243  m2a.Nullify();
1244  break;
1245  case 3:
1246  throw Undefined("Matrix<FloatDouble, DenseSparseCollection>"
1247  "Storage0, Prop1, Storage1, Allocator>"
1248  "::Read(istream& FileStream)");
1249  break;
1250  default:
1251  throw WrongArgument("HeterogeneousMatrixCollection<Prop0, "
1252  "Storage0, Prop1, Storage1, Allocator>"
1253  "::Read(istream& FileStream) ",
1254  "Underlying matrix A (" + to_str(i) + " ,"
1255  + to_str(j) + " ) not defined.");
1256  }
1257  }
1258 
1259 
1260  delete new_n;
1261  delete new_m;
1262 
1263 #ifdef SELDON_CHECK_IO
1264  // Checks if data was read.
1265  if (!FileStream.good())
1266  throw IOError("HeterogeneousMatrixCollection"
1267  "::Read(istream& FileStream)",
1268  "Input operation failed.");
1269 #endif
1270 
1271  }
1272 
1273 } // namespace Seldon.
1274 
1275 #define SELDON_FILE_MATRIX_HETEROGENEOUS_COLLECTION_CXX
1276 #endif
Seldon::HeterogeneousMatrixCollection::WriteText
void WriteText(string FileName) const
Writes the matrix collection in a file.
Definition: HeterogeneousMatrixCollection.cxx:1063
Seldon::HeterogeneousMatrixCollection::Print
void Print() const
Displays the matrix collection on the standard output.
Definition: HeterogeneousMatrixCollection.cxx:925
Seldon::Matrix_Base
Base class for all matrices.
Definition: Matrix_Base.hxx:143
Seldon::to_str
std::string to_str(const T &input)
Converts most types to string.
Definition: CommonInline.cxx:137
Seldon::WrongCol
Definition: Errors.hxx:138
Seldon::HeterogeneousMatrixCollection::Read
void Read(string FileName)
Reads the matrix collection from a file.
Definition: HeterogeneousMatrixCollection.cxx:1169
Seldon::Matrix
Definition: SeldonHeader.hxx:226
Seldon::HeterogeneousMatrixCollection::Reallocate
void Reallocate(int i, int j)
Reallocates memory to resize the matrix collection.
Definition: HeterogeneousMatrixCollection.cxx:211
Seldon::HeterogeneousMatrixCollection::operator=
HeterogeneousMatrixCollection< Prop0, Storage0, Prop1, Storage1, Allocator > & operator=(const HeterogeneousMatrixCollection< Prop0, Storage0, Prop1, Storage1, Allocator > &A)
Duplicates a matrix collection (assignment operator).
Definition: HeterogeneousMatrixCollection.cxx:832
Seldon::WrongRow
Definition: Errors.hxx:126
Seldon::HeterogeneousMatrixCollection
Matrix class made of an heterogeneous collection of matrices.
Definition: HeterogeneousMatrixCollection.hxx:42
Seldon::HeterogeneousMatrixCollection::HeterogeneousMatrixCollection
HeterogeneousMatrixCollection()
Default constructor.
Definition: HeterogeneousMatrixCollection.cxx:47
Seldon::Undefined
Definition: Errors.hxx:62
Seldon::HeterogeneousMatrixCollection::Write
void Write(string FileName, bool with_size) const
Writes the matrix collection in a file.
Definition: HeterogeneousMatrixCollection.cxx:949
Seldon::HeterogeneousMatrixCollection::Copy
void Copy(const HeterogeneousMatrixCollection< Prop0, Storage0, Prop1, Storage1, Allocator > &A)
Duplicates a matrix collection (assignment operator).
Definition: HeterogeneousMatrixCollection.cxx:850
Seldon::HeterogeneousMatrixCollection::Clear
void Clear()
Clears the matrix collection without releasing memory.
Definition: HeterogeneousMatrixCollection.cxx:92
Seldon::WrongArgument
Definition: Errors.hxx:76
Seldon::WrongDim
Definition: Errors.hxx:102
Seldon::HeterogeneousMatrixCollection::Deallocate
void Deallocate()
Deallocates underlying the matrices.
Definition: HeterogeneousMatrixCollection.cxx:185
Seldon
Seldon namespace.
Definition: Array.cxx:24
Seldon::IOError
Definition: Errors.hxx:150
Seldon::HeterogeneousMatrixCollection::Nullify
void Nullify()
Clears the matrix collection without releasing memory.
Definition: HeterogeneousMatrixCollection.cxx:115