Matrix_ComplexSparse.cxx
1 // Copyright (C) 2001-2011 Vivien Mallet
2 // Copyright (C) 2003-2011 Marc Duruflé
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_MATRIX_COMPLEXSPARSE_CXX
22 
23 #include "Matrix_ComplexSparse.hxx"
24 
25 namespace Seldon
26 {
27 
28 
29  /****************
30  * CONSTRUCTORS *
31  ****************/
32 
33 
35 
38  template <class T, class Prop, class Storage, class Allocator>
39  Matrix_ComplexSparse<T, Prop, Storage, Allocator>
40  ::Matrix_ComplexSparse(): Matrix_Base<T, Allocator>()
41  {
42  real_nz_ = 0;
43  imag_nz_ = 0;
44  real_ptr_ = NULL;
45  imag_ptr_ = NULL;
46  real_ind_ = NULL;
47  imag_ind_ = NULL;
48  real_data_ = NULL;
49  imag_data_ = NULL;
50  }
51 
52 
54 
59  template <class T, class Prop, class Storage, class Allocator>
60  Matrix_ComplexSparse<T, Prop, Storage, Allocator>
61  ::Matrix_ComplexSparse(int i, int j): Matrix_Base<T, Allocator>()
62  {
63  real_nz_ = 0;
64  imag_nz_ = 0;
65  real_ptr_ = NULL;
66  imag_ptr_ = NULL;
67  real_ind_ = NULL;
68  imag_ind_ = NULL;
69  real_data_ = NULL;
70  imag_data_ = NULL;
71 
72  Reallocate(i, j);
73  }
74 
75 
77 
87  template <class T, class Prop, class Storage, class Allocator>
88  Matrix_ComplexSparse<T, Prop, Storage, Allocator>::
89  Matrix_ComplexSparse(int i, int j, long real_nz, long imag_nz):
90  Matrix_Base<T, Allocator>()
91  {
92  real_nz_ = 0;
93  imag_nz_ = 0;
94  real_ptr_ = NULL;
95  imag_ptr_ = NULL;
96  real_ind_ = NULL;
97  imag_ind_ = NULL;
98  real_data_ = NULL;
99  imag_data_ = NULL;
100 
101  Reallocate(i, j, real_nz, imag_nz);
102  }
103 
104 
106 
124  template <class T, class Prop, class Storage, class Allocator>
125  template <class Storage0, class Allocator0,
126  class Storage1, class Allocator1,
127  class Storage2, class Allocator2>
128  Matrix_ComplexSparse<T, Prop, Storage, Allocator>::
129  Matrix_ComplexSparse(int i, int j,
130  Vector<value_type, Storage0, Allocator0>& real_values,
131  Vector<long, Storage1, Allocator1>& real_ptr,
132  Vector<int, Storage2, Allocator2>& real_ind,
133  Vector<value_type, Storage0, Allocator0>& imag_values,
134  Vector<long, Storage1, Allocator1>& imag_ptr,
135  Vector<int, Storage2, Allocator2>& imag_ind):
136  Matrix_Base<T, Allocator>(i, j)
137  {
138 
139  real_nz_ = real_values.GetM();
140  imag_nz_ = imag_values.GetM();
141 
142 #ifdef SELDON_CHECK_DIMENSIONS
143  // Checks whether vector sizes are acceptable.
144 
145  if (real_ind.GetM() != real_nz_)
146  {
147  this->m_ = 0;
148  this->n_ = 0;
149  real_nz_ = 0;
150  imag_nz_ = 0;
151  real_ptr_ = NULL;
152  imag_ptr_ = NULL;
153  real_ind_ = NULL;
154  imag_ind_ = NULL;
155  this->real_data_ = NULL;
156  this->imag_data_ = NULL;
157  throw WrongDim(string("Matrix_ComplexSparse::")
158  + string("Matrix_ComplexSparse(int, int, ")
159  + string("const Vector&, const Vector&, const Vector&")
160  + ", const Vector&, const Vector&, const Vector&)",
161  string("There are ") + to_str(real_nz_)
162  + " values (real part) but "
163  + to_str(real_ind.GetM())
164  + " row or column indices.");
165  }
166 
167  if (imag_ind.GetM() != imag_nz_)
168  {
169  this->m_ = 0;
170  this->n_ = 0;
171  real_nz_ = 0;
172  imag_nz_ = 0;
173  real_ptr_ = NULL;
174  imag_ptr_ = NULL;
175  real_ind_ = NULL;
176  imag_ind_ = NULL;
177  this->real_data_ = NULL;
178  this->imag_data_ = NULL;
179  throw WrongDim(string("Matrix_ComplexSparse::")
180  + string("Matrix_ComplexSparse(int, int, ")
181  + string("const Vector&, const Vector&, const Vector&")
182  + ", const Vector&, const Vector&, const Vector&)",
183  string("There are ") + to_str(imag_nz_)
184  + " values (imaginary part) but "
185  + to_str(imag_ind.GetM())
186  + " row or column indices.");
187  }
188 
189  if (real_ptr.GetM()-1 != Storage::GetFirst(i, j))
190  {
191  this->m_ = 0;
192  this->n_ = 0;
193  real_nz_ = 0;
194  imag_nz_ = 0;
195  real_ptr_ = NULL;
196  imag_ptr_ = NULL;
197  real_ind_ = NULL;
198  imag_ind_ = NULL;
199  this->real_data_ = NULL;
200  this->imag_data_ = NULL;
201  throw WrongDim(string("Matrix_ComplexSparse::")
202  + string("Matrix_ComplexSparse(int, int, ")
203  + string("const Vector&, const Vector&, const Vector&")
204  + ", const Vector&, const Vector&, const Vector&)",
205  string("The vector of start indices (real part)")
206  + " contains " + to_str(real_ptr.GetM()-1)
207  + string(" row or column start indices (plus the")
208  + " number of non-zero entries) but there are "
209  + to_str(Storage::GetFirst(i, j))
210  + " rows or columns ("
211  + to_str(i) + " by " + to_str(j) + " matrix).");
212  }
213 
214  if (imag_ptr.GetM()-1 != Storage::GetFirst(i, j))
215  {
216  this->m_ = 0;
217  this->n_ = 0;
218  real_nz_ = 0;
219  imag_nz_ = 0;
220  real_ptr_ = NULL;
221  imag_ptr_ = NULL;
222  real_ind_ = NULL;
223  imag_ind_ = NULL;
224  this->real_data_ = NULL;
225  this->imag_data_ = NULL;
226  throw WrongDim(string("Matrix_ComplexSparse::")
227  + string("Matrix_ComplexSparse(int, int, ")
228  + string("const Vector&, const Vector&, const Vector&")
229  + ", const Vector&, const Vector&, const Vector&)",
230  string("The vector of start indices (imaginary part)")
231  + " contains " + to_str(imag_ptr.GetM()-1)
232  + string(" row or column start indices (plus the")
233  + " number of non-zero entries) but there are "
234  + to_str(Storage::GetFirst(i, j))
235  + " rows or columns ("
236  + to_str(i) + " by " + to_str(j) + " matrix).");
237  }
238 
239  if ((real_nz_ > 0
240  && (j == 0
241  || static_cast<long int>(real_nz_-1) / static_cast<long int>(j)
242  >= static_cast<long int>(i)))
243  ||
244  (imag_nz_ > 0
245  && (j == 0
246  || static_cast<long int>(imag_nz_-1) / static_cast<long int>(j)
247  >= static_cast<long int>(i))))
248  {
249  this->m_ = 0;
250  this->n_ = 0;
251  real_nz_ = 0;
252  imag_nz_ = 0;
253  real_ptr_ = NULL;
254  imag_ptr_ = NULL;
255  real_ind_ = NULL;
256  imag_ind_ = NULL;
257  this->real_data_ = NULL;
258  this->imag_data_ = NULL;
259  throw WrongDim(string("Matrix_ComplexSparse::")
260  + string("Matrix_ComplexSparse(int, int, ")
261  + string("const Vector&, const Vector&, const Vector&")
262  + ", const Vector&, const Vector&, const Vector&)",
263  string("There are more values (")
264  + to_str(real_values.GetM())
265  + " values for the real part and "
266  + to_str(real_values.GetM()) + string(" values")
267  + string(" for the imaginary part) than elements")
268  + " in the matrix ("
269  + to_str(i) + " by " + to_str(j) + ").");
270  }
271 #endif
272 
273  this->real_ptr_ = real_ptr.GetData();
274  this->imag_ptr_ = imag_ptr.GetData();
275  this->real_ind_ = real_ind.GetData();
276  this->imag_ind_ = imag_ind.GetData();
277  this->real_data_ = real_values.GetData();
278  this->imag_data_ = imag_values.GetData();
279 
280  real_ptr.Nullify();
281  imag_ptr.Nullify();
282  real_ind.Nullify();
283  imag_ind.Nullify();
284  real_values.Nullify();
285  imag_values.Nullify();
286  }
287 
288 
290  template <class T, class Prop, class Storage, class Allocator>
291  Matrix_ComplexSparse<T, Prop, Storage, Allocator>
292  ::Matrix_ComplexSparse(const Matrix_ComplexSparse<T, Prop,
293  Storage, Allocator>& A)
294  : Matrix_Base<T, Allocator>()
295  {
296  this->m_ = 0;
297  this->n_ = 0;
298  real_nz_ = 0;
299  imag_nz_ = 0;
300  real_ptr_ = NULL;
301  imag_ptr_ = NULL;
302  real_ind_ = NULL;
303  imag_ind_ = NULL;
304  real_data_ = NULL;
305  imag_data_ = NULL;
306 
307  this->Copy(A);
308  }
309 
310 
312 
315  template <class T, class Prop, class Storage, class Allocator>
316  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::Clear()
317  {
318 
319 #ifdef SELDON_CHECK_MEMORY
320  try
321  {
322 #endif
323 
324  if (real_ptr_ != NULL)
325  {
326  AllocatorLong::deallocate(real_ptr_, this->m_+1);
327  real_ptr_ = NULL;
328  }
329 
330 #ifdef SELDON_CHECK_MEMORY
331  }
332  catch (...)
333  {
334  real_ptr_ = NULL;
335  }
336 #endif
337 
338 #ifdef SELDON_CHECK_MEMORY
339  try
340  {
341 #endif
342 
343  if (imag_ptr_ != NULL)
344  {
345  AllocatorLong::deallocate(imag_ptr_, this->m_+1);
346  imag_ptr_ = NULL;
347  }
348 
349 #ifdef SELDON_CHECK_MEMORY
350  }
351  catch (...)
352  {
353  imag_ptr_ = NULL;
354  }
355 #endif
356 
357 #ifdef SELDON_CHECK_MEMORY
358  try
359  {
360 #endif
361 
362  if (real_ind_ != NULL)
363  {
364  AllocatorInt::deallocate(real_ind_, this->real_nz_);
365  real_ind_ = NULL;
366  }
367 
368 #ifdef SELDON_CHECK_MEMORY
369  }
370  catch (...)
371  {
372  real_ind_ = NULL;
373  }
374 #endif
375 
376 #ifdef SELDON_CHECK_MEMORY
377  try
378  {
379 #endif
380 
381  if (imag_ind_ != NULL)
382  {
383  AllocatorInt::deallocate(imag_ind_, this->imag_nz_);
384  imag_ind_ = NULL;
385  }
386 
387 #ifdef SELDON_CHECK_MEMORY
388  }
389  catch (...)
390  {
391  imag_ind_ = NULL;
392  }
393 #endif
394 
395 #ifdef SELDON_CHECK_MEMORY
396  try
397  {
398 #endif
399 
400  if (this->real_data_ != NULL)
401  {
402  Allocator::deallocate(this->real_data_, real_nz_);
403  this->real_data_ = NULL;
404  }
405 
406 #ifdef SELDON_CHECK_MEMORY
407  }
408  catch (...)
409  {
410  this->real_nz_ = 0;
411  this->real_data_ = NULL;
412  }
413 #endif
414 
415 #ifdef SELDON_CHECK_MEMORY
416  try
417  {
418 #endif
419 
420  if (this->imag_data_ != NULL)
421  {
422  Allocator::deallocate(this->imag_data_, imag_nz_);
423  this->imag_data_ = NULL;
424  }
425 
426 #ifdef SELDON_CHECK_MEMORY
427  }
428  catch (...)
429  {
430  this->imag_nz_ = 0;
431  this->imag_data_ = NULL;
432  }
433 #endif
434 
435  this->m_ = 0;
436  this->n_ = 0;
437  this->real_nz_ = 0;
438  this->imag_nz_ = 0;
439  }
440 
441 
442  /*******************
443  * BASIC FUNCTIONS *
444  *******************/
445 
446 
448  template<class T, class Prop, class Storage, class Allocator>
449  size_t Matrix_ComplexSparse<T, Prop, Storage, Allocator>
450  ::GetMemorySize() const
451  {
452  size_t taille = sizeof(*this) + 2*this->GetRealPtrSize()*sizeof(long);
453  int coef = sizeof(value_type) + sizeof(int); // for each non-zero entry
454  taille += coef*size_t(this->real_nz_ + this->imag_nz_);
455  return taille;
456  }
457 
458 
459  /*********************
460  * MEMORY MANAGEMENT *
461  *********************/
462 
463 
465 
482  template <class T, class Prop, class Storage, class Allocator>
483  template <class Storage0, class Allocator0,
484  class Storage1, class Allocator1,
485  class Storage2, class Allocator2>
486  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::
487  SetData(int i, int j,
488  Vector<value_type, Storage0, Allocator0>& real_values,
489  Vector<long, Storage1, Allocator1>& real_ptr,
490  Vector<int, Storage2, Allocator2>& real_ind,
491  Vector<value_type, Storage0, Allocator0>& imag_values,
492  Vector<long, Storage1, Allocator1>& imag_ptr,
493  Vector<int, Storage2, Allocator2>& imag_ind)
494  {
495  this->Clear();
496  this->m_ = i;
497  this->n_ = j;
498  real_nz_ = real_values.GetM();
499  imag_nz_ = imag_values.GetM();
500 
501 #ifdef SELDON_CHECK_DIMENSIONS
502  // Checks whether vector sizes are acceptable.
503 
504  if (real_ind.GetM() != real_nz_)
505  {
506  this->m_ = 0;
507  this->n_ = 0;
508  real_nz_ = 0;
509  imag_nz_ = 0;
510  real_ptr_ = NULL;
511  imag_ptr_ = NULL;
512  real_ind_ = NULL;
513  imag_ind_ = NULL;
514  this->real_data_ = NULL;
515  this->imag_data_ = NULL;
516  throw WrongDim(string("Matrix_ComplexSparse::SetData(int, int, ")
517  + string("const Vector&, const Vector&, const Vector&")
518  + ", const Vector&, const Vector&, const Vector&)",
519  string("There are ") + to_str(real_nz_)
520  + " values (real part) but "
521  + to_str(real_ind.GetM())
522  + " row or column indices.");
523  }
524 
525  if (imag_ind.GetM() != imag_nz_)
526  {
527  this->m_ = 0;
528  this->n_ = 0;
529  real_nz_ = 0;
530  imag_nz_ = 0;
531  real_ptr_ = NULL;
532  imag_ptr_ = NULL;
533  real_ind_ = NULL;
534  imag_ind_ = NULL;
535  this->real_data_ = NULL;
536  this->imag_data_ = NULL;
537  throw WrongDim(string("Matrix_ComplexSparse::SetData(int, int, ")
538  + string("const Vector&, const Vector&, const Vector&")
539  + ", const Vector&, const Vector&, const Vector&)",
540  string("There are ") + to_str(imag_nz_)
541  + " values (imaginary part) but "
542  + to_str(imag_ind.GetM())
543  + " row or column indices.");
544  }
545 
546  if (real_ptr.GetM()-1 != Storage::GetFirst(i, j))
547  {
548  this->m_ = 0;
549  this->n_ = 0;
550  real_nz_ = 0;
551  imag_nz_ = 0;
552  real_ptr_ = NULL;
553  imag_ptr_ = NULL;
554  real_ind_ = NULL;
555  imag_ind_ = NULL;
556  this->real_data_ = NULL;
557  this->imag_data_ = NULL;
558  throw WrongDim(string("Matrix_ComplexSparse::SetData(int, int, ")
559  + string("const Vector&, const Vector&, const Vector&")
560  + ", const Vector&, const Vector&, const Vector&)",
561  string("The vector of start indices (real part)")
562  + " contains " + to_str(real_ptr.GetM()-1)
563  + string(" row or column start indices (plus the")
564  + " number of non-zero entries) but there are "
565  + to_str(Storage::GetFirst(i, j))
566  + " rows or columns ("
567  + to_str(i) + " by " + to_str(j) + " matrix).");
568  }
569 
570  if (imag_ptr.GetM()-1 != Storage::GetFirst(i, j))
571  {
572  this->m_ = 0;
573  this->n_ = 0;
574  real_nz_ = 0;
575  imag_nz_ = 0;
576  real_ptr_ = NULL;
577  imag_ptr_ = NULL;
578  real_ind_ = NULL;
579  imag_ind_ = NULL;
580  this->real_data_ = NULL;
581  this->imag_data_ = NULL;
582  throw WrongDim(string("Matrix_ComplexSparse::SetData(int, int, ")
583  + string("const Vector&, const Vector&, const Vector&")
584  + ", const Vector&, const Vector&, const Vector&)",
585  string("The vector of start indices (imaginary part)")
586  + " contains " + to_str(imag_ptr.GetM()-1)
587  + string(" row or column start indices (plus the")
588  + " number of non-zero entries) but there are "
589  + to_str(Storage::GetFirst(i, j))
590  + " rows or columns ("
591  + to_str(i) + " by " + to_str(j) + " matrix).");
592  }
593 
594  if ((real_nz_ > 0
595  && (j == 0
596  || static_cast<long int>(real_nz_-1) / static_cast<long int>(j)
597  >= static_cast<long int>(i)))
598  ||
599  (imag_nz_ > 0
600  && (j == 0
601  || static_cast<long int>(imag_nz_-1) / static_cast<long int>(j)
602  >= static_cast<long int>(i))))
603  {
604  this->m_ = 0;
605  this->n_ = 0;
606  real_nz_ = 0;
607  imag_nz_ = 0;
608  real_ptr_ = NULL;
609  imag_ptr_ = NULL;
610  real_ind_ = NULL;
611  imag_ind_ = NULL;
612  this->real_data_ = NULL;
613  this->imag_data_ = NULL;
614  throw WrongDim(string("Matrix_ComplexSparse::SetData(int, int, ")
615  + string("const Vector&, const Vector&, const Vector&")
616  + ", const Vector&, const Vector&, const Vector&)",
617  string("There are more values (")
618  + to_str(real_values.GetM())
619  + " values for the real part and "
620  + to_str(real_values.GetM()) + string(" values")
621  + string(" for the imaginary part) than elements")
622  + " in the matrix ("
623  + to_str(i) + " by " + to_str(j) + ").");
624  }
625 #endif
626 
627  this->real_ptr_ = real_ptr.GetData();
628  this->imag_ptr_ = imag_ptr.GetData();
629  this->real_ind_ = real_ind.GetData();
630  this->imag_ind_ = imag_ind.GetData();
631  this->real_data_ = real_values.GetData();
632  this->imag_data_ = imag_values.GetData();
633 
634  real_ptr.Nullify();
635  imag_ptr.Nullify();
636  real_ind.Nullify();
637  imag_ind.Nullify();
638  real_values.Nullify();
639  imag_values.Nullify();
640  }
641 
642 
644 
665  template <class T, class Prop, class Storage, class Allocator>
666  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::
667  SetData(int i, int j, long real_nz,
668  typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
669  ::pointer real_values,
670  long* real_ptr, int* real_ind, long imag_nz,
671  typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
672  ::pointer imag_values,
673  long* imag_ptr, int* imag_ind)
674  {
675  this->Clear();
676 
677  this->m_ = i;
678  this->n_ = j;
679 
680  this->real_nz_ = real_nz;
681  this->imag_nz_ = imag_nz;
682 
683  real_data_ = real_values;
684  imag_data_ = imag_values;
685  real_ind_ = real_ind;
686  imag_ind_ = imag_ind;
687  real_ptr_ = real_ptr;
688  imag_ptr_ = imag_ptr;
689  }
690 
691 
693 
697  template <class T, class Prop, class Storage, class Allocator>
698  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::Nullify()
699  {
700  this->data_ = NULL;
701  this->m_ = 0;
702  this->n_ = 0;
703  real_nz_ = 0;
704  real_ptr_ = NULL;
705  real_ind_ = NULL;
706  imag_nz_ = 0;
707  imag_ptr_ = NULL;
708  imag_ind_ = NULL;
709  real_data_ = NULL;
710  imag_data_ = NULL;
711  }
712 
713 
715  template <class T, class Prop, class Storage, class Allocator>
716  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
717  ::Reallocate(int i, int j)
718  {
719  // previous entries are removed
720  Clear();
721 
722  this->m_ = i;
723  this->n_ = j;
724 
725  // we try to allocate real_ptr_ and imag_ptr_
726 #ifdef SELDON_CHECK_MEMORY
727  try
728  {
729 #endif
730 
731  real_ptr_
732  = reinterpret_cast<long*>( AllocatorLong::
733  allocate(Storage::GetFirst(i, j)+1) );
734 
735  imag_ptr_
736  = reinterpret_cast<long*>( AllocatorLong::
737  allocate(Storage::GetFirst(i, j)+1) );
738 
739 #ifdef SELDON_CHECK_MEMORY
740  }
741  catch (...)
742  {
743  this->m_ = 0;
744  this->n_ = 0;
745  real_nz_ = 0;
746  real_ptr_ = NULL;
747  real_ind_ = NULL;
748  imag_nz_ = 0;
749  imag_ptr_ = NULL;
750  imag_ind_ = NULL;
751  real_data_ = NULL;
752  imag_data_ = NULL;
753  }
754  if ((real_ptr_ == NULL) || (imag_ptr_ == NULL))
755  {
756  this->m_ = 0;
757  this->n_ = 0;
758  real_nz_ = 0;
759  real_ptr_ = NULL;
760  real_ind_ = NULL;
761  imag_nz_ = 0;
762  imag_ptr_ = NULL;
763  imag_ind_ = NULL;
764  real_data_ = NULL;
765  imag_data_ = NULL;
766  }
767  if (((real_ptr_ == NULL) || (imag_ptr_ == NULL)) && i != 0 && j != 0)
768  throw NoMemory("Matrix_ComplexSparse::Reallocate(int, int)",
769  string("Unable to allocate ")
770  + to_str(sizeof(int) * (Storage::GetFirst(i, j)+1) )
771  + " bytes to store " + to_str(Storage::GetFirst(i, j)+1)
772  + " row or column start indices, for a "
773  + to_str(i) + " by " + to_str(j) + " matrix.");
774 #endif
775 
776  // then filing real_ptr_ with 0
777  for (int k = 0; k <= Storage::GetFirst(i, j); k++)
778  {
779  real_ptr_[k] = 0;
780  imag_ptr_[k] = 0;
781  }
782  }
783 
784 
786  template <class T, class Prop, class Storage, class Allocator>
787  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
788  ::Reallocate(int i, int j, long real_nz, long imag_nz)
789  {
790  // previous entries are removed
791  Clear();
792 
793  this->m_ = i;
794  this->n_ = j;
795  this->real_nz_ = real_nz;
796  this->imag_nz_ = imag_nz;
797 
798 #ifdef SELDON_CHECK_DIMENSIONS
799  if (real_nz_ < 0 || imag_nz_ < 0)
800  {
801  this->m_ = 0;
802  this->n_ = 0;
803  real_nz_ = 0;
804  imag_nz_ = 0;
805  real_ptr_ = NULL;
806  imag_ptr_ = NULL;
807  real_ind_ = NULL;
808  imag_ind_ = NULL;
809  this->real_data_ = NULL;
810  this->imag_data_ = NULL;
811  throw WrongDim(string("Matrix_ComplexSparse::")
812  + "Matrix_ComplexSparse(int, int, int, int)",
813  "Invalid number of non-zero elements: "
814  + to_str(real_nz) + " in the real part and "
815  + to_str(imag_nz) + " in the imaginary part.");
816  }
817  if ((real_nz_ > 0
818  && (j == 0
819  || static_cast<long int>(real_nz_-1) / static_cast<long int>(j)
820  >= static_cast<long int>(i)))
821  ||
822  (imag_nz_ > 0
823  && (j == 0
824  || static_cast<long int>(imag_nz_-1) / static_cast<long int>(j)
825  >= static_cast<long int>(i))))
826  {
827  this->m_ = 0;
828  this->n_ = 0;
829  real_nz_ = 0;
830  imag_nz_ = 0;
831  real_ptr_ = NULL;
832  imag_ptr_ = NULL;
833  real_ind_ = NULL;
834  imag_ind_ = NULL;
835  this->real_data_ = NULL;
836  this->imag_data_ = NULL;
837  throw WrongDim(string("Matrix_ComplexSparse::")
838  + "Matrix_ComplexSparse(int, int, int, int)",
839  string("There are more values (") + to_str(real_nz)
840  + " values for the real part and " + to_str(imag_nz)
841  + string(" values for the imaginary part) than")
842  + " elements in the matrix (" + to_str(i) + " by "
843  + to_str(j) + ").");
844  }
845 #endif
846 
847 #ifdef SELDON_CHECK_MEMORY
848  try
849  {
850 #endif
851 
852  real_ptr_
853  = reinterpret_cast<long*>( AllocatorLong::
854  allocate(Storage::GetFirst(i, j)+1) );
855 
856 #ifdef SELDON_CHECK_MEMORY
857  }
858  catch (...)
859  {
860  this->m_ = 0;
861  this->n_ = 0;
862  real_nz_ = 0;
863  imag_nz_ = 0;
864  real_ptr_ = NULL;
865  imag_ptr_ = NULL;
866  real_ind_ = NULL;
867  imag_ind_ = NULL;
868  this->real_data_ = NULL;
869  this->imag_data_ = NULL;
870  }
871  if (real_ptr_ == NULL)
872  {
873  this->m_ = 0;
874  this->n_ = 0;
875  real_nz_ = 0;
876  imag_nz_ = 0;
877  imag_ptr_ = 0;
878  real_ind_ = NULL;
879  imag_ind_ = NULL;
880  this->real_data_ = NULL;
881  this->imag_data_ = NULL;
882  }
883  if (real_ptr_ == NULL && i != 0 && j != 0)
884  throw NoMemory(string("Matrix_ComplexSparse::")
885  + "Matrix_ComplexSparse(int, int, int, int)",
886  string("Unable to allocate ")
887  + to_str(sizeof(int) * (Storage::GetFirst(i, j)+1))
888  + " bytes to store " + to_str(Storage::GetFirst(i, j)+1)
889  + string(" row or column start indices (for the real")
890  + " part), for a "
891  + to_str(i) + " by " + to_str(j) + " matrix.");
892 #endif
893 
894 #ifdef SELDON_CHECK_MEMORY
895  try
896  {
897 #endif
898 
899  imag_ptr_
900  = reinterpret_cast<long*>( AllocatorLong::
901  allocate(Storage::GetFirst(i, j)+1) );
902 
903 #ifdef SELDON_CHECK_MEMORY
904  }
905  catch (...)
906  {
907  this->m_ = 0;
908  this->n_ = 0;
909  real_nz_ = 0;
910  imag_nz_ = 0;
911  AllocatorLong::deallocate(real_ptr_,
912  Storage::GetFirst(i, j)+1);
913  real_ptr_ = NULL;
914  imag_ptr_ = NULL;
915  real_ind_ = NULL;
916  imag_ind_ = NULL;
917  this->real_data_ = NULL;
918  this->imag_data_ = NULL;
919  }
920  if (imag_ptr_ == NULL)
921  {
922  this->m_ = 0;
923  this->n_ = 0;
924  real_nz_ = 0;
925  imag_nz_ = 0;
926  AllocatorLong::
927  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
928  real_ptr_ = 0;
929  real_ind_ = NULL;
930  imag_ind_ = NULL;
931  this->real_data_ = NULL;
932  this->imag_data_ = NULL;
933  }
934  if (imag_ptr_ == NULL && i != 0 && j != 0)
935  throw NoMemory(string("Matrix_ComplexSparse::")
936  + "Matrix_ComplexSparse(int, int, int, int)",
937  string("Unable to allocate ")
938  + to_str(sizeof(int) * (Storage::GetFirst(i, j)+1))
939  + " bytes to store " + to_str(Storage::GetFirst(i, j)+1)
940  + string(" row or column start indices (for the")
941  + string(" imaginary part), for a ")
942  + to_str(i) + " by " + to_str(j) + " matrix.");
943 #endif
944 
945 #ifdef SELDON_CHECK_MEMORY
946  try
947  {
948 #endif
949 
950  real_ind_ =
951  reinterpret_cast<int*>( AllocatorInt::allocate(real_nz_, this) );
952 
953 #ifdef SELDON_CHECK_MEMORY
954  }
955  catch (...)
956  {
957  this->m_ = 0;
958  this->n_ = 0;
959  real_nz_ = 0;
960  imag_nz_ = 0;
961  AllocatorLong::
962  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
963  AllocatorLong::
964  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
965  real_ptr_ = NULL;
966  imag_ptr_ = NULL;
967  real_ind_ = NULL;
968  imag_ind_ = NULL;
969  this->real_data_ = NULL;
970  this->imag_data_ = NULL;
971  }
972  if (real_ind_ == NULL)
973  {
974  this->m_ = 0;
975  this->n_ = 0;
976  real_nz_ = 0;
977  imag_nz_ = 0;
978  AllocatorLong::
979  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
980  AllocatorLong::
981  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
982  real_ptr_ = NULL;
983  imag_ptr_ = NULL;
984  this->real_data_ = NULL;
985  this->imag_data_ = NULL;
986  }
987  if (real_ind_ == NULL && i != 0 && j != 0)
988  throw NoMemory(string("Matrix_ComplexSparse::")
989  + "Matrix_ComplexSparse(int, int, int, int)",
990  string("Unable to allocate ")
991  + to_str(sizeof(int) * real_nz)
992  + " bytes to store " + to_str(real_nz)
993  + " row or column indices (real part), for a "
994  + to_str(i) + " by " + to_str(j) + " matrix.");
995 #endif
996 
997 #ifdef SELDON_CHECK_MEMORY
998  try
999  {
1000 #endif
1001 
1002  imag_ind_ =
1003  reinterpret_cast<int*>( AllocatorInt::allocate(imag_nz_, this) );
1004 
1005 #ifdef SELDON_CHECK_MEMORY
1006  }
1007  catch (...)
1008  {
1009  this->m_ = 0;
1010  this->n_ = 0;
1011  real_nz_ = 0;
1012  imag_nz_ = 0;
1013  AllocatorLong::
1014  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1015  AllocatorLong::
1016  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1017  real_ptr_ = NULL;
1018  imag_ptr_ = NULL;
1019  AllocatorInt::deallocate(imag_ind_, imag_nz);
1020  real_ind_ = NULL;
1021  imag_ind_ = NULL;
1022  this->real_data_ = NULL;
1023  this->imag_data_ = NULL;
1024  }
1025  if (imag_ind_ == NULL)
1026  {
1027  this->m_ = 0;
1028  this->n_ = 0;
1029  real_nz_ = 0;
1030  imag_nz_ = 0;
1031  AllocatorLong::
1032  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1033  AllocatorLong::
1034  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1035  real_ptr_ = NULL;
1036  imag_ptr_ = NULL;
1037  AllocatorInt::deallocate(imag_ind_, imag_nz);
1038  imag_ind_ = NULL;
1039  this->real_data_ = NULL;
1040  this->imag_data_ = NULL;
1041  }
1042  if (imag_ind_ == NULL && i != 0 && j != 0)
1043  throw NoMemory(string("Matrix_ComplexSparse::")
1044  + "Matrix_ComplexSparse(int, int, int, int)",
1045  string("Unable to allocate ")
1046  + to_str(sizeof(int) * imag_nz)
1047  + " bytes to store " + to_str(imag_nz)
1048  + " row or column indices (imaginary part), for a "
1049  + to_str(i) + " by " + to_str(j) + " matrix.");
1050 #endif
1051 
1052 #ifdef SELDON_CHECK_MEMORY
1053  try
1054  {
1055 #endif
1056 
1057  this->real_data_ = Allocator::allocate(real_nz_, this);
1058 
1059 #ifdef SELDON_CHECK_MEMORY
1060  }
1061  catch (...)
1062  {
1063  this->m_ = 0;
1064  this->n_ = 0;
1065  AllocatorLong::
1066  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1067  AllocatorLong::
1068  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1069  real_ptr_ = NULL;
1070  imag_ptr_ = NULL;
1071  AllocatorInt::deallocate(real_ind_, real_nz);
1072  AllocatorInt::deallocate(imag_ind_, imag_nz);
1073  real_ind_ = NULL;
1074  imag_ind_ = NULL;
1075  this->real_data_ = NULL;
1076  this->imag_data_ = NULL;
1077  }
1078  if (real_data_ == NULL)
1079  {
1080  this->m_ = 0;
1081  this->n_ = 0;
1082  AllocatorLong::
1083  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1084  AllocatorLong::
1085  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1086  real_ptr_ = NULL;
1087  imag_ptr_ = NULL;
1088  AllocatorInt::deallocate(real_ind_, real_nz);
1089  AllocatorInt::deallocate(imag_ind_, imag_nz);
1090  real_ind_ = NULL;
1091  imag_ind_ = NULL;
1092  imag_data_ = NULL;
1093  }
1094  if (real_data_ == NULL && i != 0 && j != 0)
1095  throw NoMemory(string("Matrix_ComplexSparse::")
1096  + "Matrix_ComplexSparse(int, int, int, int)",
1097  string("Unable to allocate ")
1098  + to_str(sizeof(int) * real_nz)
1099  + " bytes to store " + to_str(real_nz)
1100  + " values (real part), for a "
1101  + to_str(i) + " by " + to_str(j) + " matrix.");
1102 #endif
1103 
1104 #ifdef SELDON_CHECK_MEMORY
1105  try
1106  {
1107 #endif
1108 
1109  this->imag_data_ = Allocator::allocate(imag_nz_, this);
1110 
1111 #ifdef SELDON_CHECK_MEMORY
1112  }
1113  catch (...)
1114  {
1115  this->m_ = 0;
1116  this->n_ = 0;
1117  AllocatorLong::
1118  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1119  AllocatorLong::
1120  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1121  real_ptr_ = NULL;
1122  imag_ptr_ = NULL;
1123  AllocatorInt::deallocate(real_ind_, real_nz);
1124  AllocatorInt::deallocate(imag_ind_, imag_nz);
1125  real_ind_ = NULL;
1126  imag_ind_ = NULL;
1127  Allocator::deallocate(this->real_data_, real_nz_);
1128  this->real_data_ = NULL;
1129  this->imag_data_ = NULL;
1130  }
1131  if (real_data_ == NULL)
1132  {
1133  this->m_ = 0;
1134  this->n_ = 0;
1135  AllocatorLong::
1136  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1137  AllocatorLong::
1138  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1139  real_ptr_ = NULL;
1140  imag_ptr_ = NULL;
1141  AllocatorInt::deallocate(real_ind_, real_nz);
1142  AllocatorInt::deallocate(imag_ind_, imag_nz);
1143  real_ind_ = NULL;
1144  imag_ind_ = NULL;
1145  Allocator::deallocate(this->real_data_, real_nz_);
1146  real_data_ = NULL;
1147  }
1148  if (imag_data_ == NULL && i != 0 && j != 0)
1149  throw NoMemory(string("Matrix_ComplexSparse::")
1150  + "Matrix_ComplexSparse(int, int, int, int)",
1151  string("Unable to allocate ")
1152  + to_str(sizeof(int) * imag_nz)
1153  + " bytes to store " + to_str(imag_nz)
1154  + " values (imaginary part), for a "
1155  + to_str(i) + " by " + to_str(j) + " matrix.");
1156 #endif
1157 
1158  // then filing real_ptr_ with 0
1159  for (int k = 0; k <= Storage::GetFirst(i, j); k++)
1160  {
1161  real_ptr_[k] = 0;
1162  imag_ptr_[k] = 0;
1163  }
1164  }
1165 
1166 
1168 
1172  template <class T, class Prop, class Storage, class Allocator>
1173  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
1174  ::Resize(int i, int j)
1175  {
1176  if (Storage::GetFirst(i, j) < Storage::GetFirst(this->m_, this->n_))
1177  Resize(i, j, real_ptr_[Storage::GetFirst(i, j)],
1178  imag_ptr_[Storage::GetFirst(i, j)]);
1179  else
1180  Resize(i, j, real_nz_, imag_nz_);
1181  }
1182 
1183 
1185 
1190  template <class T, class Prop, class Storage, class Allocator>
1191  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
1192  ::Resize(int i, int j, long real_nz, long imag_nz)
1193  {
1194 #ifdef SELDON_CHECK_DIMENSIONS
1195  if (real_nz < 0 || imag_nz < 0)
1196  {
1197  this->m_ = 0;
1198  this->n_ = 0;
1199  real_nz_ = 0;
1200  imag_nz_ = 0;
1201  real_ptr_ = NULL;
1202  imag_ptr_ = NULL;
1203  real_ind_ = NULL;
1204  imag_ind_ = NULL;
1205  this->real_data_ = NULL;
1206  this->imag_data_ = NULL;
1207  throw WrongDim(string("Matrix_ComplexSparse::")
1208  + "Resize(int, int, int, int)",
1209  "Invalid number of non-zero elements: "
1210  + to_str(real_nz) + " in the real part and "
1211  + to_str(imag_nz) + " in the imaginary part.");
1212  }
1213  if ((real_nz > 0
1214  && (j == 0
1215  || static_cast<long int>(real_nz-1) / static_cast<long int>(j)
1216  >= static_cast<long int>(i)))
1217  ||
1218  (imag_nz > 0
1219  && (j == 0
1220  || static_cast<long int>(imag_nz-1) / static_cast<long int>(j)
1221  >= static_cast<long int>(i))))
1222  {
1223  this->m_ = 0;
1224  this->n_ = 0;
1225  real_nz_ = 0;
1226  imag_nz_ = 0;
1227  real_ptr_ = NULL;
1228  imag_ptr_ = NULL;
1229  real_ind_ = NULL;
1230  imag_ind_ = NULL;
1231  this->real_data_ = NULL;
1232  this->imag_data_ = NULL;
1233  throw WrongDim(string("Matrix_ComplexSparse::")
1234  + "Resize(int, int, int, int)",
1235  string("There are more values (") + to_str(real_nz)
1236  + " values for the real part and " + to_str(imag_nz)
1237  + string(" values for the imaginary part) than")
1238  + " elements in the matrix (" + to_str(i) + " by "
1239  + to_str(j) + ").");
1240  }
1241 #endif
1242 
1243  if (Storage::GetFirst(this->m_, this->n_) != Storage::GetFirst(i, j))
1244  {
1245 #ifdef SELDON_CHECK_MEMORY
1246  try
1247  {
1248 #endif
1249 
1250  real_ptr_ =
1251  reinterpret_cast<long*>( AllocatorLong::
1252  reallocate(real_ptr_,
1253  (Storage::GetFirst(i, j)+1)) );
1254 
1255 #ifdef SELDON_CHECK_MEMORY
1256  }
1257  catch (...)
1258  {
1259  this->m_ = 0;
1260  this->n_ = 0;
1261  real_nz_ = 0;
1262  imag_nz_ = 0;
1263  real_ptr_ = NULL;
1264  imag_ptr_ = NULL;
1265  real_ind_ = NULL;
1266  imag_ind_ = NULL;
1267  this->real_data_ = NULL;
1268  this->imag_data_ = NULL;
1269  }
1270  if (real_ptr_ == NULL)
1271  {
1272  this->m_ = 0;
1273  this->n_ = 0;
1274  real_nz_ = 0;
1275  imag_nz_ = 0;
1276  imag_ptr_ = 0;
1277  real_ind_ = NULL;
1278  imag_ind_ = NULL;
1279  this->real_data_ = NULL;
1280  this->imag_data_ = NULL;
1281  }
1282  if (real_ptr_ == NULL && i != 0 && j != 0)
1283  throw NoMemory(string("Matrix_ComplexSparse::")
1284  + "Resize(int, int, int, int)",
1285  string("Unable to allocate ")
1286  + to_str(sizeof(int) * (Storage::GetFirst(i, j)+1))
1287  + " bytes to store " + to_str(Storage::GetFirst(i, j)+1)
1288  + string(" row or column start indices (for the real")
1289  + " part), for a "
1290  + to_str(i) + " by " + to_str(j) + " matrix.");
1291 #endif
1292 
1293 #ifdef SELDON_CHECK_MEMORY
1294  try
1295  {
1296 #endif
1297 
1298  imag_ptr_ =
1299  reinterpret_cast<long*>( AllocatorLong::
1300  reallocate(imag_ptr_,
1301  (Storage::GetFirst(i, j)+1)) );
1302 
1303 #ifdef SELDON_CHECK_MEMORY
1304  }
1305  catch (...)
1306  {
1307  this->m_ = 0;
1308  this->n_ = 0;
1309  real_nz_ = 0;
1310  imag_nz_ = 0;
1311  AllocatorLong::
1312  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1313  real_ptr_ = NULL;
1314  imag_ptr_ = NULL;
1315  real_ind_ = NULL;
1316  imag_ind_ = NULL;
1317  this->real_data_ = NULL;
1318  this->imag_data_ = NULL;
1319  }
1320  if (imag_ptr_ == NULL)
1321  {
1322  this->m_ = 0;
1323  this->n_ = 0;
1324  real_nz_ = 0;
1325  imag_nz_ = 0;
1326  AllocatorLong::
1327  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1328  real_ptr_ = 0;
1329  real_ind_ = NULL;
1330  imag_ind_ = NULL;
1331  this->real_data_ = NULL;
1332  this->imag_data_ = NULL;
1333  }
1334  if (imag_ptr_ == NULL && i != 0 && j != 0)
1335  throw NoMemory(string("Matrix_ComplexSparse::")
1336  + "Resize(int, int, int, int)",
1337  string("Unable to allocate ")
1338  + to_str(sizeof(int) * (Storage::GetFirst(i, j)+1))
1339  + " bytes to store " + to_str(Storage::GetFirst(i, j)+1)
1340  + string(" row or column start indices (for the")
1341  + string(" imaginary part), for a ")
1342  + to_str(i) + " by " + to_str(j) + " matrix.");
1343 #endif
1344  }
1345 
1346  if (real_nz != real_nz_)
1347  {
1348 #ifdef SELDON_CHECK_MEMORY
1349  try
1350  {
1351 #endif
1352 
1353  real_ind_ =
1354  reinterpret_cast<int*>( AllocatorInt::
1355  reallocate(real_ind_,
1356  real_nz) );
1357 
1358 #ifdef SELDON_CHECK_MEMORY
1359  }
1360  catch (...)
1361  {
1362  this->m_ = 0;
1363  this->n_ = 0;
1364  real_nz_ = 0;
1365  imag_nz_ = 0;
1366  AllocatorLong::
1367  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1368  AllocatorLong::
1369  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1370  real_ptr_ = NULL;
1371  imag_ptr_ = NULL;
1372  real_ind_ = NULL;
1373  imag_ind_ = NULL;
1374  this->real_data_ = NULL;
1375  this->imag_data_ = NULL;
1376  }
1377  if (real_ind_ == NULL)
1378  {
1379  this->m_ = 0;
1380  this->n_ = 0;
1381  real_nz_ = 0;
1382  imag_nz_ = 0;
1383  AllocatorLong::
1384  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1385  AllocatorLong::
1386  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1387  real_ptr_ = NULL;
1388  imag_ptr_ = NULL;
1389  this->real_data_ = NULL;
1390  this->imag_data_ = NULL;
1391  }
1392  if (real_ind_ == NULL && i != 0 && j != 0)
1393  throw NoMemory(string("Matrix_ComplexSparse::")
1394  + "Resize(int, int, int, int)",
1395  string("Unable to allocate ")
1396  + to_str(sizeof(int) * real_nz)
1397  + " bytes to store " + to_str(real_nz)
1398  + " row or column indices (real part), for a "
1399  + to_str(i) + " by " + to_str(j) + " matrix.");
1400 #endif
1401  }
1402 
1403  if (imag_nz != imag_nz_)
1404  {
1405 #ifdef SELDON_CHECK_MEMORY
1406  try
1407  {
1408 #endif
1409 
1410  imag_ind_ =
1411  reinterpret_cast<int*>( AllocatorInt::
1412  reallocate(imag_ind_,
1413  imag_nz) );
1414 #ifdef SELDON_CHECK_MEMORY
1415  }
1416  catch (...)
1417  {
1418  this->m_ = 0;
1419  this->n_ = 0;
1420  real_nz_ = 0;
1421  imag_nz_ = 0;
1422  AllocatorLong::
1423  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1424  AllocatorLong::
1425  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1426  real_ptr_ = NULL;
1427  imag_ptr_ = NULL;
1428  AllocatorInt::deallocate(imag_ind_, imag_nz);
1429  real_ind_ = NULL;
1430  imag_ind_ = NULL;
1431  this->real_data_ = NULL;
1432  this->imag_data_ = NULL;
1433  }
1434  if (imag_ind_ == NULL)
1435  {
1436  this->m_ = 0;
1437  this->n_ = 0;
1438  real_nz_ = 0;
1439  imag_nz_ = 0;
1440  AllocatorLong::
1441  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1442  AllocatorLong::
1443  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1444  real_ptr_ = NULL;
1445  imag_ptr_ = NULL;
1446  AllocatorInt::deallocate(imag_ind_, imag_nz);
1447  imag_ind_ = NULL;
1448  this->real_data_ = NULL;
1449  this->imag_data_ = NULL;
1450  }
1451  if (imag_ind_ == NULL && i != 0 && j != 0)
1452  throw NoMemory(string("Matrix_ComplexSparse::")
1453  + "Resize(int, int, int, int)",
1454  string("Unable to allocate ")
1455  + to_str(sizeof(int) * imag_nz)
1456  + " bytes to store " + to_str(imag_nz)
1457  + " row or column indices (imaginary part), for a "
1458  + to_str(i) + " by " + to_str(j) + " matrix.");
1459 #endif
1460  }
1461 
1462  if (real_nz != real_nz_)
1463  {
1464  Vector<value_type, VectFull, Allocator> val;
1465  val.SetData(real_nz_, real_data_);
1466  val.Resize(real_nz);
1467 
1468  real_data_ = val.GetData();
1469  val.Nullify();
1470  }
1471 
1472  if (imag_nz != imag_nz_)
1473  {
1474  Vector<value_type, VectFull, Allocator> val;
1475  val.SetData(imag_nz_, imag_data_);
1476  val.Resize(imag_nz);
1477 
1478  imag_data_ = val.GetData();
1479  val.Nullify();
1480  }
1481 
1482  // then filing last values of ptr_ with nz_
1483  for (int k = Storage::GetFirst(this->m_, this->n_);
1484  k <= Storage::GetFirst(i, j); k++)
1485  {
1486  real_ptr_[k] = real_nz_;
1487  imag_ptr_[k] = imag_nz_;
1488  }
1489 
1490  this->m_ = i;
1491  this->n_ = j;
1492  real_nz_ = real_nz;
1493  imag_nz_ = imag_nz;
1494  }
1495 
1496 
1498  template <class T, class Prop, class Storage, class Allocator>
1499  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::
1500  Copy(const Matrix_ComplexSparse<T, Prop, Storage, Allocator>& A)
1501  {
1502  this->Clear();
1503  int i = A.m_;
1504  int j = A.n_;
1505  real_nz_ = A.real_nz_;
1506  imag_nz_ = A.imag_nz_;
1507  this->m_ = i;
1508  this->n_ = j;
1509  if ((i == 0)||(j == 0))
1510  {
1511  this->m_ = 0;
1512  this->n_ = 0;
1513  this->real_nz_ = 0;
1514  this->imag_nz_ = 0;
1515  return;
1516  }
1517 
1518 #ifdef SELDON_CHECK_DIMENSIONS
1519  if ((real_nz_ > 0
1520  && (j == 0
1521  || static_cast<long int>(real_nz_-1) / static_cast<long int>(j)
1522  >= static_cast<long int>(i)))
1523  ||
1524  (imag_nz_ > 0
1525  && (j == 0
1526  || static_cast<long int>(imag_nz_-1) / static_cast<long int>(j)
1527  >= static_cast<long int>(i))))
1528  {
1529  this->m_ = 0;
1530  this->n_ = 0;
1531  real_nz_ = 0;
1532  imag_nz_ = 0;
1533  real_ptr_ = NULL;
1534  imag_ptr_ = NULL;
1535  real_ind_ = NULL;
1536  imag_ind_ = NULL;
1537  this->real_data_ = NULL;
1538  this->imag_data_ = NULL;
1539  throw WrongDim(string("Matrix_ComplexSparse::")
1540  + "Matrix_ComplexSparse(int, int, int, int)",
1541  string("There are more values (") + to_str(real_nz_)
1542  + " values for the real part and " + to_str(imag_nz_)
1543  + string(" values for the imaginary part) than")
1544  + " elements in the matrix (" + to_str(i) + " by "
1545  + to_str(j) + ").");
1546  }
1547 #endif
1548 
1549 #ifdef SELDON_CHECK_MEMORY
1550  try
1551  {
1552 #endif
1553 
1554  real_ptr_ =
1555  reinterpret_cast<long*>( AllocatorLong::
1556  allocate(Storage::GetFirst(i, j)+1) );
1557 
1558  AllocatorLong::
1559  memorycpy(this->real_ptr_, A.real_ptr_,
1560  (Storage::GetFirst(i, j) + 1) );
1561 
1562 #ifdef SELDON_CHECK_MEMORY
1563  }
1564  catch (...)
1565  {
1566  this->m_ = 0;
1567  this->n_ = 0;
1568  real_nz_ = 0;
1569  imag_nz_ = 0;
1570  real_ptr_ = NULL;
1571  imag_ptr_ = NULL;
1572  real_ind_ = NULL;
1573  imag_ind_ = NULL;
1574  this->real_data_ = NULL;
1575  this->imag_data_ = NULL;
1576  }
1577  if (real_ptr_ == NULL)
1578  {
1579  this->m_ = 0;
1580  this->n_ = 0;
1581  real_nz_ = 0;
1582  imag_nz_ = 0;
1583  imag_ptr_ = 0;
1584  real_ind_ = NULL;
1585  imag_ind_ = NULL;
1586  this->real_data_ = NULL;
1587  this->imag_data_ = NULL;
1588  }
1589  if (real_ptr_ == NULL && i != 0 && j != 0)
1590  throw NoMemory(string("Matrix_ComplexSparse::")
1591  + "Matrix_ComplexSparse(int, int, int, int)",
1592  string("Unable to allocate ")
1593  + to_str(sizeof(int) * (Storage::GetFirst(i, j)+1))
1594  + " bytes to store " + to_str(Storage::GetFirst(i, j)+1)
1595  + string(" row or column start indices (for the real")
1596  + " part), for a "
1597  + to_str(i) + " by " + to_str(j) + " matrix.");
1598 #endif
1599 
1600 #ifdef SELDON_CHECK_MEMORY
1601  try
1602  {
1603 #endif
1604 
1605  imag_ptr_ =
1606  reinterpret_cast<long*>( AllocatorLong::
1607  allocate(Storage::GetFirst(i, j)+1) );
1608 
1609  AllocatorLong::
1610  memorycpy(this->imag_ptr_, A.imag_ptr_,
1611  (Storage::GetFirst(i, j) + 1) );
1612 
1613 #ifdef SELDON_CHECK_MEMORY
1614  }
1615  catch (...)
1616  {
1617  this->m_ = 0;
1618  this->n_ = 0;
1619  real_nz_ = 0;
1620  imag_nz_ = 0;
1621  AllocatorLong::
1622  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1623  real_ptr_ = NULL;
1624  imag_ptr_ = NULL;
1625  real_ind_ = NULL;
1626  imag_ind_ = NULL;
1627  this->real_data_ = NULL;
1628  this->imag_data_ = NULL;
1629  }
1630  if (imag_ptr_ == NULL)
1631  {
1632  this->m_ = 0;
1633  this->n_ = 0;
1634  real_nz_ = 0;
1635  imag_nz_ = 0;
1636  AllocatorLong::
1637  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1638  real_ptr_ = 0;
1639  real_ind_ = NULL;
1640  imag_ind_ = NULL;
1641  this->real_data_ = NULL;
1642  this->imag_data_ = NULL;
1643  }
1644  if (imag_ptr_ == NULL && i != 0 && j != 0)
1645  throw NoMemory(string("Matrix_ComplexSparse::")
1646  + "Matrix_ComplexSparse(int, int, int, int)",
1647  string("Unable to allocate ")
1648  + to_str(sizeof(int) * (Storage::GetFirst(i, j)+1))
1649  + " bytes to store " + to_str(Storage::GetFirst(i, j)+1)
1650  + string(" row or column start indices (for the")
1651  + string(" imaginary part), for a ")
1652  + to_str(i) + " by " + to_str(j) + " matrix.");
1653 #endif
1654 
1655 #ifdef SELDON_CHECK_MEMORY
1656  try
1657  {
1658 #endif
1659 
1660  real_ind_ =
1661  reinterpret_cast<int*>( AllocatorInt::
1662  allocate(real_nz_, this) );
1663 
1664  AllocatorInt::
1665  memorycpy(this->real_ind_, A.real_ind_, real_nz_);
1666 
1667 #ifdef SELDON_CHECK_MEMORY
1668  }
1669  catch (...)
1670  {
1671  this->m_ = 0;
1672  this->n_ = 0;
1673  real_nz_ = 0;
1674  imag_nz_ = 0;
1675  AllocatorLong::
1676  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1677  AllocatorLong::
1678  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1679  real_ptr_ = NULL;
1680  imag_ptr_ = NULL;
1681  real_ind_ = NULL;
1682  imag_ind_ = NULL;
1683  this->real_data_ = NULL;
1684  this->imag_data_ = NULL;
1685  }
1686  if (real_ind_ == NULL)
1687  {
1688  this->m_ = 0;
1689  this->n_ = 0;
1690  real_nz_ = 0;
1691  imag_nz_ = 0;
1692  AllocatorLong::
1693  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1694  AllocatorLong::
1695  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1696  real_ptr_ = NULL;
1697  imag_ptr_ = NULL;
1698  this->real_data_ = NULL;
1699  this->imag_data_ = NULL;
1700  }
1701  if (real_ind_ == NULL && i != 0 && j != 0)
1702  throw NoMemory(string("Matrix_ComplexSparse::")
1703  + "Matrix_ComplexSparse(int, int, int, int)",
1704  string("Unable to allocate ")
1705  + to_str(sizeof(int) * real_nz_)
1706  + " bytes to store " + to_str(real_nz_)
1707  + " row or column indices (real part), for a "
1708  + to_str(i) + " by " + to_str(j) + " matrix.");
1709 #endif
1710 
1711 #ifdef SELDON_CHECK_MEMORY
1712  try
1713  {
1714 #endif
1715 
1716  imag_ind_ =
1717  reinterpret_cast<int*>( AllocatorInt::
1718  allocate(imag_nz_, this) );
1719 
1720  AllocatorInt::
1721  memorycpy(this->imag_ind_, A.imag_ind_, imag_nz_);
1722 
1723 #ifdef SELDON_CHECK_MEMORY
1724  }
1725  catch (...)
1726  {
1727  this->m_ = 0;
1728  this->n_ = 0;
1729  real_nz_ = 0;
1730  imag_nz_ = 0;
1731  AllocatorLong::
1732  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1733  AllocatorLong::
1734  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1735  real_ptr_ = NULL;
1736  imag_ptr_ = NULL;
1737  AllocatorInt::deallocate(imag_ind_, A.imag_nz_);
1738  real_ind_ = NULL;
1739  imag_ind_ = NULL;
1740  this->real_data_ = NULL;
1741  this->imag_data_ = NULL;
1742  }
1743  if (real_ind_ == NULL)
1744  {
1745  this->m_ = 0;
1746  this->n_ = 0;
1747  real_nz_ = 0;
1748  imag_nz_ = 0;
1749  AllocatorLong::
1750  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1751  AllocatorLong::
1752  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1753  real_ptr_ = NULL;
1754  imag_ptr_ = NULL;
1755  AllocatorInt::deallocate(imag_ind_, A.imag_nz_);
1756  imag_ind_ = NULL;
1757  this->real_data_ = NULL;
1758  this->imag_data_ = NULL;
1759  }
1760  if (imag_ind_ == NULL && i != 0 && j != 0)
1761  throw NoMemory(string("Matrix_ComplexSparse::")
1762  + "Matrix_ComplexSparse(int, int, int, int)",
1763  string("Unable to allocate ")
1764  + to_str(sizeof(int) * imag_nz_)
1765  + " bytes to store " + to_str(imag_nz_)
1766  + " row or column indices (imaginary part), for a "
1767  + to_str(i) + " by " + to_str(j) + " matrix.");
1768 #endif
1769 
1770 #ifdef SELDON_CHECK_MEMORY
1771  try
1772  {
1773 #endif
1774 
1775  this->real_data_ = Allocator::allocate(real_nz_, this);
1776  Allocator::memorycpy(this->real_data_, A.real_data_, real_nz_);
1777 
1778 #ifdef SELDON_CHECK_MEMORY
1779  }
1780  catch (...)
1781  {
1782  this->m_ = 0;
1783  this->n_ = 0;
1784  AllocatorLong::
1785  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1786  AllocatorLong::
1787  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1788  real_ptr_ = NULL;
1789  imag_ptr_ = NULL;
1790  AllocatorInt::deallocate(real_ind_, A.real_nz_);
1791  AllocatorInt::deallocate(imag_ind_, A.imag_nz_);
1792  real_ind_ = NULL;
1793  imag_ind_ = NULL;
1794  this->real_data_ = NULL;
1795  this->imag_data_ = NULL;
1796  }
1797  if (real_data_ == NULL)
1798  {
1799  this->m_ = 0;
1800  this->n_ = 0;
1801  AllocatorLong::
1802  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1803  AllocatorLong::
1804  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1805  real_ptr_ = NULL;
1806  imag_ptr_ = NULL;
1807  AllocatorInt::deallocate(real_ind_, A.real_nz_);
1808  AllocatorInt::deallocate(imag_ind_, A.imag_nz_);
1809  real_ind_ = NULL;
1810  imag_ind_ = NULL;
1811  imag_data_ = NULL;
1812  }
1813  if (real_data_ == NULL && i != 0 && j != 0)
1814  throw NoMemory(string("Matrix_ComplexSparse::")
1815  + "Matrix_ComplexSparse(int, int, int, int)",
1816  string("Unable to allocate ")
1817  + to_str(sizeof(int) * real_nz_)
1818  + " bytes to store " + to_str(real_nz_)
1819  + " values (real part), for a "
1820  + to_str(i) + " by " + to_str(j) + " matrix.");
1821 #endif
1822 
1823 #ifdef SELDON_CHECK_MEMORY
1824  try
1825  {
1826 #endif
1827 
1828  this->imag_data_ = Allocator::allocate(imag_nz_, this);
1829  Allocator::memorycpy(this->imag_data_, A.imag_data_, imag_nz_);
1830 
1831 #ifdef SELDON_CHECK_MEMORY
1832  }
1833  catch (...)
1834  {
1835  this->m_ = 0;
1836  this->n_ = 0;
1837  AllocatorLong::
1838  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1839  AllocatorLong::
1840  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1841  real_ptr_ = NULL;
1842  imag_ptr_ = NULL;
1843  AllocatorInt::deallocate(real_ind_, A.real_nz_);
1844  AllocatorInt::deallocate(imag_ind_, A.imag_nz_);
1845  real_ind_ = NULL;
1846  imag_ind_ = NULL;
1847  Allocator::deallocate(this->real_data_, real_nz_);
1848  this->real_data_ = NULL;
1849  this->imag_data_ = NULL;
1850  }
1851  if (real_data_ == NULL)
1852  {
1853  this->m_ = 0;
1854  this->n_ = 0;
1855  AllocatorLong::
1856  deallocate(real_ptr_, Storage::GetFirst(i, j)+1);
1857  AllocatorLong::
1858  deallocate(imag_ptr_, Storage::GetFirst(i, j)+1);
1859  real_ptr_ = NULL;
1860  imag_ptr_ = NULL;
1861  AllocatorInt::deallocate(real_ind_, A.real_nz_);
1862  AllocatorInt::deallocate(imag_ind_, A.imag_nz_);
1863  real_ind_ = NULL;
1864  imag_ind_ = NULL;
1865  Allocator::deallocate(this->real_data_, real_nz_);
1866  real_data_ = NULL;
1867  }
1868  if (imag_data_ == NULL && i != 0 && j != 0)
1869  throw NoMemory(string("Matrix_ComplexSparse::")
1870  + "Matrix_ComplexSparse(int, int, int, int)",
1871  string("Unable to allocate ")
1872  + to_str(sizeof(int) * imag_nz_)
1873  + " bytes to store " + to_str(imag_nz_)
1874  + " values (imaginary part), for a "
1875  + to_str(i) + " by " + to_str(j) + " matrix.");
1876 #endif
1877 
1878  }
1879 
1880 
1881  /**********************************
1882  * ELEMENT ACCESS AND AFFECTATION *
1883  **********************************/
1884 
1885 
1887 
1893  template <class T, class Prop, class Storage, class Allocator>
1894  const typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>::entry_type
1895  Matrix_ComplexSparse<T, Prop, Storage, Allocator>
1896  ::operator() (int i, int j) const
1897  {
1898 
1899 #ifdef SELDON_CHECK_BOUNDS
1900  if (i < 0 || i >= this->m_)
1901  throw WrongRow("Matrix_ComplexSparse::operator()",
1902  string("Index should be in [0, ") + to_str(this->m_-1)
1903  + "], but is equal to " + to_str(i) + ".");
1904  if (j < 0 || j >= this->n_)
1905  throw WrongCol("Matrix_ComplexSparse::operator()",
1906  string("Index should be in [0, ") + to_str(this->n_-1)
1907  + "], but is equal to " + to_str(j) + ".");
1908 #endif
1909 
1910  long real_k, imag_k; int l;
1911  long real_a, real_b;
1912  long imag_a, imag_b;
1913 
1914  real_a = real_ptr_[Storage::GetFirst(i, j)];
1915  real_b = real_ptr_[Storage::GetFirst(i, j) + 1];
1916 
1917  imag_a = imag_ptr_[Storage::GetFirst(i, j)];
1918  imag_b = imag_ptr_[Storage::GetFirst(i, j) + 1];
1919  if (real_a != real_b)
1920  {
1921  l = Storage::GetSecond(i, j);
1922  for (real_k = real_a;
1923  (real_k <real_b-1) && (real_ind_[real_k] < l);
1924  real_k++);
1925  if (imag_a != imag_b)
1926  {
1927  for (imag_k = imag_a;
1928  (imag_k < imag_b-1) && (imag_ind_[imag_k] < l);
1929  imag_k++);
1930  if (real_ind_[real_k] == l)
1931  {
1932  if (imag_ind_[imag_k] == l)
1933  return entry_type(real_data_[real_k], imag_data_[imag_k]);
1934  else
1935  return entry_type(real_data_[real_k], value_type(0));
1936  }
1937  else
1938  if (imag_ind_[imag_k] == l)
1939  return entry_type(value_type(0), imag_data_[imag_k]);
1940  else
1941  return entry_type(value_type(0), value_type(0));
1942  }
1943  else
1944  {
1945  if (real_ind_[real_k] == l)
1946  return entry_type(real_data_[real_k], value_type(0));
1947  else
1948  return entry_type(value_type(0), value_type(0));
1949  }
1950  }
1951  else
1952  {
1953  if (imag_a != imag_b)
1954  {
1955  l = Storage::GetSecond(i, j);
1956  for (imag_k = imag_a;
1957  (imag_k < imag_b-1) && (imag_ind_[imag_k] < l);
1958  imag_k++);
1959  if (imag_ind_[imag_k] == l)
1960  return entry_type(value_type(0), imag_data_[imag_k]);
1961  else
1962  return entry_type(value_type(0), value_type(0));
1963  }
1964  else
1965  return entry_type(value_type(0), value_type(0));
1966  }
1967 
1968  }
1969 
1970 
1972 
1980  template <class T, class Prop, class Storage, class Allocator>
1981  typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
1982  ::value_type&
1983  Matrix_ComplexSparse<T, Prop, Storage, Allocator>::ValReal(int i, int j)
1984  {
1985 
1986 #ifdef SELDON_CHECK_BOUNDS
1987  if (i < 0 || i >= this->m_)
1988  throw WrongRow("Matrix_ComplexSparse::ValReal(int, int)",
1989  string("Index should be in [0, ") + to_str(this->m_-1)
1990  + "], but is equal to " + to_str(i) + ".");
1991  if (j < 0 || j >= this->n_)
1992  throw WrongCol("Matrix_ComplexSparse::ValReal(int, int)",
1993  string("Index should be in [0, ") + to_str(this->n_-1)
1994  + "], but is equal to " + to_str(j) + ".");
1995 #endif
1996 
1997  long k; int l;
1998  long a, b;
1999 
2000  a = real_ptr_[Storage::GetFirst(i, j)];
2001  b = real_ptr_[Storage::GetFirst(i, j) + 1];
2002 
2003  if (a == b)
2004  throw WrongArgument("Matrix_ComplexSparse::ValReal(int, int)",
2005  "No reference to element (" + to_str(i) + ", "
2006  + to_str(j)
2007  + ") can be returned: it is a zero entry.");
2008 
2009  l = Storage::GetSecond(i, j);
2010 
2011  for (k = a; (k < b-1) && (real_ind_[k] < l); k++);
2012 
2013  if (real_ind_[k] == l)
2014  return this->real_data_[k];
2015  else
2016  throw WrongArgument("Matrix_ComplexSparse::ValReal(int, int)",
2017  "No reference to element (" + to_str(i) + ", "
2018  + to_str(j)
2019  + ") can be returned: it is a zero entry.");
2020  }
2021 
2022 
2024 
2032  template <class T, class Prop, class Storage, class Allocator>
2033  const typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
2034  ::value_type&
2035  Matrix_ComplexSparse<T, Prop, Storage, Allocator>::ValReal(int i, int j) const
2036  {
2037 
2038 #ifdef SELDON_CHECK_BOUNDS
2039  if (i < 0 || i >= this->m_)
2040  throw WrongRow("Matrix_ComplexSparse::ValReal(int, int)",
2041  string("Index should be in [0, ") + to_str(this->m_-1)
2042  + "], but is equal to " + to_str(i) + ".");
2043  if (j < 0 || j >= this->n_)
2044  throw WrongCol("Matrix_ComplexSparse::ValReal(int, int)",
2045  string("Index should be in [0, ") + to_str(this->n_-1)
2046  + "], but is equal to " + to_str(j) + ".");
2047 #endif
2048 
2049  long k; int l;
2050  long a, b;
2051 
2052  a = real_ptr_[Storage::GetFirst(i, j)];
2053  b = real_ptr_[Storage::GetFirst(i, j) + 1];
2054 
2055  if (a == b)
2056  throw WrongArgument("Matrix_ComplexSparse::ValReal(int, int)",
2057  "No reference to element (" + to_str(i) + ", "
2058  + to_str(j)
2059  + ") can be returned: it is a zero entry.");
2060 
2061  l = Storage::GetSecond(i, j);
2062 
2063  for (k = a; (k < b-1) && (real_ind_[k] < l); k++);
2064 
2065  if (real_ind_[k] == l)
2066  return this->real_data_[k];
2067  else
2068  throw WrongArgument("Matrix_ComplexSparse::ValReal(int, int)",
2069  "No reference to element (" + to_str(i) + ", "
2070  + to_str(j)
2071  + ") can be returned: it is a zero entry.");
2072  }
2073 
2074 
2076 
2084  template <class T, class Prop, class Storage, class Allocator>
2085  typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
2086  ::value_type&
2087  Matrix_ComplexSparse<T, Prop, Storage, Allocator>::ValImag(int i, int j)
2088  {
2089 
2090 #ifdef SELDON_CHECK_BOUNDS
2091  if (i < 0 || i >= this->m_)
2092  throw WrongRow("Matrix_ComplexSparse::ValImag(int, int)",
2093  string("Index should be in [0, ") + to_str(this->m_-1)
2094  + "], but is equal to " + to_str(i) + ".");
2095  if (j < 0 || j >= this->n_)
2096  throw WrongCol("Matrix_ComplexSparse::ValImag(int, int)",
2097  string("Index should be in [0, ") + to_str(this->n_-1)
2098  + "], but is equal to " + to_str(j) + ".");
2099 #endif
2100 
2101  long k; int l;
2102  long a, b;
2103 
2104  a = imag_ptr_[Storage::GetFirst(i, j)];
2105  b = imag_ptr_[Storage::GetFirst(i, j) + 1];
2106 
2107  if (a == b)
2108  throw WrongArgument("Matrix_ComplexSparse::ValImag(int, int)",
2109  "No reference to element (" + to_str(i) + ", "
2110  + to_str(j)
2111  + ") can be returned: it is a zero entry.");
2112 
2113  l = Storage::GetSecond(i, j);
2114 
2115  for (k = a; (k < b-1) && (imag_ind_[k] < l); k++);
2116 
2117  if (imag_ind_[k] == l)
2118  return this->imag_data_[k];
2119  else
2120  throw WrongArgument("Matrix_ComplexSparse::ValImag(int, int)",
2121  "No reference to element (" + to_str(i) + ", "
2122  + to_str(j)
2123  + ") can be returned: it is a zero entry.");
2124  }
2125 
2126 
2128 
2136  template <class T, class Prop, class Storage, class Allocator>
2137  const typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
2138  ::value_type&
2139  Matrix_ComplexSparse<T, Prop, Storage, Allocator>::ValImag(int i, int j) const
2140  {
2141 
2142 #ifdef SELDON_CHECK_BOUNDS
2143  if (i < 0 || i >= this->m_)
2144  throw WrongRow("Matrix_ComplexSparse::ValImag(int, int)",
2145  string("Index should be in [0, ") + to_str(this->m_-1)
2146  + "], but is equal to " + to_str(i) + ".");
2147  if (j < 0 || j >= this->n_)
2148  throw WrongCol("Matrix_ComplexSparse::ValImag(int, int)",
2149  string("Index should be in [0, ") + to_str(this->n_-1)
2150  + "], but is equal to " + to_str(j) + ".");
2151 #endif
2152 
2153  long k; int l;
2154  long a, b;
2155 
2156  a = imag_ptr_[Storage::GetFirst(i, j)];
2157  b = imag_ptr_[Storage::GetFirst(i, j) + 1];
2158 
2159  if (a == b)
2160  throw WrongArgument("Matrix_ComplexSparse::ValImag(int, int)",
2161  "No reference to element (" + to_str(i) + ", "
2162  + to_str(j)
2163  + ") can be returned: it is a zero entry.");
2164 
2165  l = Storage::GetSecond(i, j);
2166 
2167  for (k = a; (k < b-1) && (imag_ind_[k] < l); k++);
2168 
2169  if (imag_ind_[k] == l)
2170  return this->imag_data_[k];
2171  else
2172  throw WrongArgument("Matrix_ComplexSparse::ValImag(int, int)",
2173  "No reference to element (" + to_str(i) + ", "
2174  + to_str(j)
2175  + ") can be returned: it is a zero entry.");
2176  }
2177 
2178 
2180 
2187  template <class T, class Prop, class Storage, class Allocator>
2188  typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
2189  ::value_type&
2190  Matrix_ComplexSparse<T, Prop, Storage, Allocator>::GetReal(int i, int j)
2191  {
2192 
2193 #ifdef SELDON_CHECK_BOUNDS
2194  if (i < 0 || i >= this->m_)
2195  throw WrongRow("Matrix_ComplexSparse::GetReal(int, int)",
2196  string("Index should be in [0, ") + to_str(this->m_-1)
2197  + "], but is equal to " + to_str(i) + ".");
2198  if (j < 0 || j >= this->n_)
2199  throw WrongCol("Matrix_ComplexSparse::GetReal(int, int)",
2200  string("Index should be in [0, ") + to_str(this->n_-1)
2201  + "], but is equal to " + to_str(j) + ".");
2202 #endif
2203 
2204  long k; int l;
2205  long a, b;
2206 
2207  a = real_ptr_[Storage::GetFirst(i, j)];
2208  b = real_ptr_[Storage::GetFirst(i, j) + 1];
2209 
2210  if (a < b)
2211  {
2212  l = Storage::GetSecond(i, j);
2213 
2214  for (k = a; (k < b) && (real_ind_[k] < l); k++);
2215 
2216  if ( (k < b) && (real_ind_[k] == l))
2217  return this->real_data_[k];
2218  }
2219  else
2220  k = a;
2221 
2222  // adding a non-zero entry
2223  Resize(this->m_, this->n_, real_nz_+1, imag_nz_);
2224 
2225  for (int m = Storage::GetFirst(i, j)+1;
2226  m <= Storage::GetFirst(this->m_, this->n_); m++)
2227  real_ptr_[m]++;
2228 
2229  for (long m = real_nz_-1; m >= k+1; m--)
2230  {
2231  real_ind_[m] = real_ind_[m-1];
2232  this->real_data_[m] = this->real_data_[m-1];
2233  }
2234 
2235  real_ind_[k] = Storage::GetSecond(i, j);
2236 
2237  // value of new non-zero entry is set to 0
2238  SetComplexZero(this->real_data_[k]);
2239 
2240  return this->real_data_[k];
2241  }
2242 
2243 
2245 
2252  template <class T, class Prop, class Storage, class Allocator>
2253  typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
2254  ::value_type&
2255  Matrix_ComplexSparse<T, Prop, Storage, Allocator>::GetImag(int i, int j)
2256  {
2257 
2258 #ifdef SELDON_CHECK_BOUNDS
2259  if (i < 0 || i >= this->m_)
2260  throw WrongRow("Matrix_ComplexSparse::GetImag(int, int)",
2261  string("Index should be in [0, ") + to_str(this->m_-1)
2262  + "], but is equal to " + to_str(i) + ".");
2263  if (j < 0 || j >= this->n_)
2264  throw WrongCol("Matrix_ComplexSparse::GetImag(int, int)",
2265  string("Index should be in [0, ") + to_str(this->n_-1)
2266  + "], but is equal to " + to_str(j) + ".");
2267 #endif
2268 
2269  long k; int l;
2270  long a, b;
2271 
2272  a = imag_ptr_[Storage::GetFirst(i, j)];
2273  b = imag_ptr_[Storage::GetFirst(i, j) + 1];
2274 
2275  if (a < b)
2276  {
2277  l = Storage::GetSecond(i, j);
2278 
2279  for (k = a; (k < b) && (imag_ind_[k] < l); k++);
2280 
2281  if ( (k < b) && (imag_ind_[k] == l))
2282  return this->imag_data_[k];
2283  }
2284  else
2285  k = a;
2286 
2287  // adding a non-zero entry
2288  Resize(this->m_, this->n_, real_nz_, imag_nz_+1);
2289  for (int m = Storage::GetFirst(i, j)+1;
2290  m <= Storage::GetFirst(this->m_, this->n_); m++)
2291  imag_ptr_[m]++;
2292 
2293  for (long m = imag_nz_-1; m >= k+1; m--)
2294  {
2295  imag_ind_[m] = imag_ind_[m-1];
2296  this->imag_data_[m] = this->imag_data_[m-1];
2297  }
2298 
2299  imag_ind_[k] = Storage::GetSecond(i, j);
2300 
2301  // value of new non-zero entry is set to 0
2302  SetComplexZero(this->imag_data_[k]);
2303 
2304  return this->imag_data_[k];
2305  }
2306 
2307 
2308  /************************
2309  * CONVENIENT FUNCTIONS *
2310  ************************/
2311 
2312 
2314 
2315  template <class T, class Prop, class Storage, class Allocator>
2316  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::Zero()
2317  {
2318  Allocator::memoryset(this->real_data_, char(0),
2319  this->real_nz_ * sizeof(value_type));
2320 
2321  Allocator::memoryset(this->imag_data_, char(0),
2322  this->imag_nz_ * sizeof(value_type));
2323  }
2324 
2325 
2327 
2330  template <class T, class Prop, class Storage, class Allocator>
2331  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::SetIdentity()
2332  {
2333  int m = this->m_;
2334  int n = this->n_;
2335  int nz = min(m, n);
2336 
2337  if (nz == 0)
2338  return;
2339 
2340  Clear();
2341 
2342  Vector<value_type, VectFull, Allocator> real_values(nz), imag_values;
2343  Vector<long, VectFull, AllocatorLong>
2344  real_ptr(Storage::GetFirst(m, n) + 1);
2345  Vector<int, VectFull, AllocatorInt> real_ind(nz);
2346  Vector<long, VectFull, AllocatorLong> imag_ptr(real_ptr);
2347  Vector<int, VectFull, AllocatorInt> imag_ind;
2348 
2349  real_values.Fill(value_type(1));
2350  real_ind.Fill();
2351  imag_ind.Zero();
2352  imag_ptr.Zero();
2353  int i;
2354  for (i = 0; i < nz + 1; i++)
2355  real_ptr(i) = i;
2356 
2357  for (i = nz + 1; i < real_ptr.GetM(); i++)
2358  real_ptr(i) = nz;
2359 
2360  SetData(m, n, real_values, real_ptr, real_ind,
2361  imag_values, imag_ptr, imag_ind);
2362  }
2363 
2364 
2366 
2369  template <class T, class Prop, class Storage, class Allocator>
2370  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::Fill()
2371  {
2372  for (long i = 0; i < this->real_nz_; i++)
2373  SetComplexReal(i, this->real_data_[i]);
2374 
2375  for (long i = 0; i < this->imag_nz_; i++)
2376  this->imag_data_[i] = value_type(0);
2377  }
2378 
2379 
2381 
2384  template <class T, class Prop, class Storage, class Allocator>
2385  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
2386  ::Fill(const entry_type& x)
2387  {
2388  for (long i = 0; i < this->real_nz_; i++)
2389  this->real_data_[i] = real(x);
2390 
2391  for (long i = 0; i < this->imag_nz_; i++)
2392  this->imag_data_[i] = imag(x);
2393  }
2394 
2395 
2397 
2400  template <class T, class Prop, class Storage, class Allocator>
2401  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::FillRand()
2402  {
2403 #ifndef SELDON_WITHOUT_REINIT_RANDOM
2404  srand(time(NULL));
2405 #endif
2406  for (long i = 0; i < this->real_nz_; i++)
2407  SetComplexReal(rand(), this->real_data_[i]);
2408 
2409  for (long i = 0; i < this->imag_nz_; i++)
2410  SetComplexReal(rand(), this->imag_data_[i]);
2411  }
2412 
2413 
2415 
2420  template <class T, class Prop, class Storage, class Allocator>
2421  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::Print() const
2422  {
2423  for (int i = 0; i < this->m_; i++)
2424  {
2425  for (int j = 0; j < this->n_; j++)
2426  cout << (*this)(i, j) << "\t";
2427  cout << endl;
2428  }
2429  }
2430 
2431 
2433 
2437  template <class T, class Prop, class Storage, class Allocator>
2438  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
2439  ::Write(string FileName) const
2440  {
2441  ofstream FileStream;
2442  FileStream.open(FileName.c_str(), ofstream::binary);
2443 
2444 #ifdef SELDON_CHECK_IO
2445  // Checks if the file was opened.
2446  if (!FileStream.is_open())
2447  throw IOError("Matrix_ComplexSparse::Write(string FileName)",
2448  string("Unable to open file \"") + FileName + "\".");
2449 #endif
2450 
2451  this->Write(FileStream);
2452 
2453  FileStream.close();
2454  }
2455 
2456 
2458 
2462  template <class T, class Prop, class Storage, class Allocator>
2463  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
2464  ::Write(ostream& FileStream) const
2465  {
2466 #ifdef SELDON_CHECK_IO
2467  // Checks if the stream is ready.
2468  if (!FileStream.good())
2469  throw IOError("Matrix_ComplexSparse::Write(ofstream& FileStream)",
2470  "Stream is not ready.");
2471 #endif
2472 
2473  FileStream.write(reinterpret_cast<char*>(const_cast<int*>(&this->m_)),
2474  sizeof(int));
2475  FileStream.write(reinterpret_cast<char*>(const_cast<int*>(&this->n_)),
2476  sizeof(int));
2477  FileStream.write(reinterpret_cast<char*>(const_cast<long*>(&this->real_nz_)),
2478  sizeof(long));
2479  FileStream.write(reinterpret_cast<char*>(const_cast<long*>(&this->imag_nz_)),
2480  sizeof(long));
2481 
2482  FileStream.write(reinterpret_cast<char*>(this->real_ptr_),
2483  sizeof(long)*(Storage::GetFirst(this->m_, this->n_)+1));
2484  FileStream.write(reinterpret_cast<char*>(this->real_ind_),
2485  sizeof(int)*this->real_nz_);
2486  FileStream.write(reinterpret_cast<char*>(this->real_data_),
2487  sizeof(value_type)*this->real_nz_);
2488 
2489  FileStream.write(reinterpret_cast<char*>(this->imag_ptr_),
2490  sizeof(long)*(Storage::GetFirst(this->m_, this->n_)+1));
2491  FileStream.write(reinterpret_cast<char*>(this->imag_ind_),
2492  sizeof(int)*this->imag_nz_);
2493  FileStream.write(reinterpret_cast<char*>(this->imag_data_),
2494  sizeof(value_type)*this->imag_nz_);
2495  }
2496 
2497 
2499 
2505  template <class T, class Prop, class Storage, class Allocator>
2506  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::
2507  WriteText(string FileName, bool cplx) const
2508  {
2509  ofstream FileStream; FileStream.precision(14);
2510  FileStream.open(FileName.c_str());
2511 
2512  // changing precision
2513  FileStream.precision(cout.precision());
2514 
2515 #ifdef SELDON_CHECK_IO
2516  // Checks if the file was opened.
2517  if (!FileStream.is_open())
2518  throw IOError("Matrix_ComplexSparse::Write(string FileName)",
2519  string("Unable to open file \"") + FileName + "\".");
2520 #endif
2521 
2522  this->WriteText(FileStream, cplx);
2523 
2524  FileStream.close();
2525  }
2526 
2527 
2529 
2535  template <class T, class Prop, class Storage, class Allocator>
2536  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::
2537  WriteText(ostream& FileStream, bool cplx) const
2538  {
2539 
2540 #ifdef SELDON_CHECK_IO
2541  // Checks if the stream is ready.
2542  if (!FileStream.good())
2543  throw IOError("Matrix_ComplexSparse::Write(ofstream& FileStream)",
2544  "Stream is not ready.");
2545 #endif
2546 
2547  // conversion in coordinate format (1-index convention)
2548  const Matrix<T, Prop, Storage, Allocator>& leaf_class =
2549  static_cast<const Matrix<T, Prop, Storage, Allocator>& >(*this);
2550 
2551  entry_type zero; int index = 1;
2552  WriteCoordinateMatrix(leaf_class, FileStream, zero, index, cplx);
2553  }
2554 
2555 
2557 
2561  template <class T, class Prop, class Storage, class Allocator>
2562  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::
2563  Read(string FileName)
2564  {
2565  ifstream FileStream;
2566  FileStream.open(FileName.c_str(), ifstream::binary);
2567 
2568 #ifdef SELDON_CHECK_IO
2569  // Checks if the file was opened.
2570  if (!FileStream.is_open())
2571  throw IOError("Matrix_ComplexSparse::Read(string FileName)",
2572  string("Unable to open file \"") + FileName + "\".");
2573 #endif
2574 
2575  this->Read(FileStream);
2576 
2577  FileStream.close();
2578  }
2579 
2580 
2582 
2586  template <class T, class Prop, class Storage, class Allocator>
2587  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
2588  ::Read(istream& FileStream)
2589  {
2590 
2591 #ifdef SELDON_CHECK_IO
2592  // Checks if the stream is ready.
2593  if (!FileStream.good())
2594  throw IOError("Matrix_ComplexSparse::Read(istream& FileStream)",
2595  "Stream is not ready.");
2596 #endif
2597 
2598  int m, n, real_nz, imag_nz;
2599  FileStream.read(reinterpret_cast<char*>(&m), sizeof(int));
2600  FileStream.read(reinterpret_cast<char*>(&n), sizeof(int));
2601  FileStream.read(reinterpret_cast<char*>(&real_nz), sizeof(long));
2602  FileStream.read(reinterpret_cast<char*>(&imag_nz), sizeof(long));
2603 
2604  Reallocate(m, n, real_nz, imag_nz);
2605 
2606  FileStream.read(reinterpret_cast<char*>(real_ptr_),
2607  sizeof(long)*(Storage::GetFirst(m, n)+1));
2608  FileStream.read(reinterpret_cast<char*>(real_ind_), sizeof(int)*real_nz);
2609  FileStream.read(reinterpret_cast<char*>(this->real_data_),
2610  sizeof(value_type)*real_nz);
2611 
2612  FileStream.read(reinterpret_cast<char*>(imag_ptr_),
2613  sizeof(long)*(Storage::GetFirst(m, n)+1));
2614  FileStream.read(reinterpret_cast<char*>(imag_ind_), sizeof(int)*imag_nz);
2615  FileStream.read(reinterpret_cast<char*>(this->imag_data_),
2616  sizeof(value_type)*imag_nz);
2617 
2618 #ifdef SELDON_CHECK_IO
2619  // Checks if data was read.
2620  if (!FileStream.good())
2621  throw IOError("Matrix_ComplexSparse::Read(istream& FileStream)",
2622  string("Input operation failed.")
2623  + string(" The input file may have been removed")
2624  + " or may not contain enough data.");
2625 #endif
2626 
2627  }
2628 
2629 
2631 
2635  template <class T, class Prop, class Storage, class Allocator>
2636  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
2637  ::ReadText(string FileName, bool cplx)
2638  {
2639  ifstream FileStream;
2640  FileStream.open(FileName.c_str());
2641 
2642 #ifdef SELDON_CHECK_IO
2643  // Checks if the file was opened.
2644  if (!FileStream.is_open())
2645  throw IOError("Matrix_ComplexSparse::ReadText(string FileName)",
2646  string("Unable to open file \"") + FileName + "\".");
2647 #endif
2648 
2649  this->ReadText(FileStream, cplx);
2650 
2651  FileStream.close();
2652  }
2653 
2654 
2656 
2660  template <class T, class Prop, class Storage, class Allocator>
2661  void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
2662  ::ReadText(istream& FileStream, bool cplx)
2663  {
2664  Matrix<T, Prop, Storage, Allocator>& leaf_class =
2665  static_cast<Matrix<T, Prop, Storage, Allocator>& >(*this);
2666 
2667  entry_type zero; int index = 1;
2668  ReadCoordinateMatrix(leaf_class, FileStream, zero, index, -1, cplx);
2669  }
2670 
2671 
2672 } // namespace Seldon.
2673 
2674 #define SELDON_FILE_MATRIX_COMPLEXSPARSE_CXX
2675 #endif
Seldon::Matrix_ComplexSparse
Complex sparse-matrix class.
Definition: Matrix_ComplexSparse.hxx:112
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::Vector
Definition: SeldonHeader.hxx:207
Seldon::Matrix
Definition: SeldonHeader.hxx:226
Seldon::WrongRow
Definition: Errors.hxx:126
Seldon::Matrix_ComplexSparse::Matrix_ComplexSparse
Matrix_ComplexSparse()
Default constructor.
Definition: Matrix_ComplexSparse.cxx:40
Seldon::WrongArgument
Definition: Errors.hxx:76
Seldon::WrongDim
Definition: Errors.hxx:102
Seldon
Seldon namespace.
Definition: Array.cxx:24
Seldon::IOError
Definition: Errors.hxx:150
Seldon::NoMemory
Definition: Errors.hxx:90
Seldon::SetComplexReal
void SetComplexReal(int n, std::complex< T > &number)
Sets a complex number to (n, 0).
Definition: CommonInline.cxx:234