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