00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef SELDON_FILE_MATRIX_COMPLEXSPARSE_CXX
00022
00023 #include "Matrix_ComplexSparse.hxx"
00024
00025 namespace Seldon
00026 {
00027
00028
00029
00030
00031
00032
00033
00035
00038 template <class T, class Prop, class Storage, class Allocator>
00039 inline Matrix_ComplexSparse<T, Prop, Storage, Allocator>
00040 ::Matrix_ComplexSparse(): Matrix_Base<T, Allocator>()
00041 {
00042 real_nz_ = 0;
00043 imag_nz_ = 0;
00044 real_ptr_ = NULL;
00045 imag_ptr_ = NULL;
00046 real_ind_ = NULL;
00047 imag_ind_ = NULL;
00048 real_data_ = NULL;
00049 imag_data_ = NULL;
00050 }
00051
00052
00054
00059 template <class T, class Prop, class Storage, class Allocator>
00060 inline Matrix_ComplexSparse<T, Prop, Storage, Allocator>
00061 ::Matrix_ComplexSparse(int i, int j): Matrix_Base<T, Allocator>()
00062 {
00063 real_nz_ = 0;
00064 imag_nz_ = 0;
00065 real_ptr_ = NULL;
00066 imag_ptr_ = NULL;
00067 real_ind_ = NULL;
00068 imag_ind_ = NULL;
00069 real_data_ = NULL;
00070 imag_data_ = NULL;
00071
00072 Reallocate(i, j);
00073 }
00074
00075
00077
00087 template <class T, class Prop, class Storage, class Allocator>
00088 inline Matrix_ComplexSparse<T, Prop, Storage, Allocator>::
00089 Matrix_ComplexSparse(int i, int j, int real_nz, int imag_nz):
00090 Matrix_Base<T, Allocator>()
00091 {
00092 real_nz_ = 0;
00093 imag_nz_ = 0;
00094 real_ptr_ = NULL;
00095 imag_ptr_ = NULL;
00096 real_ind_ = NULL;
00097 imag_ind_ = NULL;
00098 real_data_ = NULL;
00099 imag_data_ = NULL;
00100
00101 Reallocate(i, j, real_nz, imag_nz);
00102 }
00103
00104
00106
00124 template <class T, class Prop, class Storage, class Allocator>
00125 template <class Storage0, class Allocator0,
00126 class Storage1, class Allocator1,
00127 class Storage2, class Allocator2>
00128 inline Matrix_ComplexSparse<T, Prop, Storage, Allocator>::
00129 Matrix_ComplexSparse(int i, int j,
00130 Vector<T, Storage0, Allocator0>& real_values,
00131 Vector<int, Storage1, Allocator1>& real_ptr,
00132 Vector<int, Storage2, Allocator2>& real_ind,
00133 Vector<T, Storage0, Allocator0>& imag_values,
00134 Vector<int, Storage1, Allocator1>& imag_ptr,
00135 Vector<int, Storage2, Allocator2>& imag_ind):
00136 Matrix_Base<T, Allocator>(i, j)
00137 {
00138
00139 real_nz_ = real_values.GetLength();
00140 imag_nz_ = imag_values.GetLength();
00141
00142 #ifdef SELDON_CHECK_DIMENSIONS
00143
00144
00145 if (real_ind.GetLength() != real_nz_)
00146 {
00147 this->m_ = 0;
00148 this->n_ = 0;
00149 real_nz_ = 0;
00150 imag_nz_ = 0;
00151 real_ptr_ = NULL;
00152 imag_ptr_ = NULL;
00153 real_ind_ = NULL;
00154 imag_ind_ = NULL;
00155 this->real_data_ = NULL;
00156 this->imag_data_ = NULL;
00157 throw WrongDim(string("Matrix_ComplexSparse::")
00158 + string("Matrix_ComplexSparse(int, int, ")
00159 + string("const Vector&, const Vector&, const Vector&")
00160 + ", const Vector&, const Vector&, const Vector&)",
00161 string("There are ") + to_str(real_nz_)
00162 + " values (real part) but "
00163 + to_str(real_ind.GetLength())
00164 + " row or column indices.");
00165 }
00166
00167 if (imag_ind.GetLength() != imag_nz_)
00168 {
00169 this->m_ = 0;
00170 this->n_ = 0;
00171 real_nz_ = 0;
00172 imag_nz_ = 0;
00173 real_ptr_ = NULL;
00174 imag_ptr_ = NULL;
00175 real_ind_ = NULL;
00176 imag_ind_ = NULL;
00177 this->real_data_ = NULL;
00178 this->imag_data_ = NULL;
00179 throw WrongDim(string("Matrix_ComplexSparse::")
00180 + string("Matrix_ComplexSparse(int, int, ")
00181 + string("const Vector&, const Vector&, const Vector&")
00182 + ", const Vector&, const Vector&, const Vector&)",
00183 string("There are ") + to_str(imag_nz_)
00184 + " values (imaginary part) but "
00185 + to_str(imag_ind.GetLength())
00186 + " row or column indices.");
00187 }
00188
00189 if (real_ptr.GetLength()-1 != Storage::GetFirst(i, j))
00190 {
00191 this->m_ = 0;
00192 this->n_ = 0;
00193 real_nz_ = 0;
00194 imag_nz_ = 0;
00195 real_ptr_ = NULL;
00196 imag_ptr_ = NULL;
00197 real_ind_ = NULL;
00198 imag_ind_ = NULL;
00199 this->real_data_ = NULL;
00200 this->imag_data_ = NULL;
00201 throw WrongDim(string("Matrix_ComplexSparse::")
00202 + string("Matrix_ComplexSparse(int, int, ")
00203 + string("const Vector&, const Vector&, const Vector&")
00204 + ", const Vector&, const Vector&, const Vector&)",
00205 string("The vector of start indices (real part)")
00206 + " contains " + to_str(real_ptr.GetLength()-1)
00207 + string(" row or column start indices (plus the")
00208 + " number of non-zero entries) but there are "
00209 + to_str(Storage::GetFirst(i, j))
00210 + " rows or columns ("
00211 + to_str(i) + " by " + to_str(j) + " matrix).");
00212 }
00213
00214 if (imag_ptr.GetLength()-1 != Storage::GetFirst(i, j))
00215 {
00216 this->m_ = 0;
00217 this->n_ = 0;
00218 real_nz_ = 0;
00219 imag_nz_ = 0;
00220 real_ptr_ = NULL;
00221 imag_ptr_ = NULL;
00222 real_ind_ = NULL;
00223 imag_ind_ = NULL;
00224 this->real_data_ = NULL;
00225 this->imag_data_ = NULL;
00226 throw WrongDim(string("Matrix_ComplexSparse::")
00227 + string("Matrix_ComplexSparse(int, int, ")
00228 + string("const Vector&, const Vector&, const Vector&")
00229 + ", const Vector&, const Vector&, const Vector&)",
00230 string("The vector of start indices (imaginary part)")
00231 + " contains " + to_str(imag_ptr.GetLength()-1)
00232 + string(" row or column start indices (plus the")
00233 + " number of non-zero entries) but there are "
00234 + to_str(Storage::GetFirst(i, j))
00235 + " rows or columns ("
00236 + to_str(i) + " by " + to_str(j) + " matrix).");
00237 }
00238
00239 if ((real_nz_ > 0
00240 && (j == 0
00241 || static_cast<long int>(real_nz_-1) / static_cast<long int>(j)
00242 >= static_cast<long int>(i)))
00243 ||
00244 (imag_nz_ > 0
00245 && (j == 0
00246 || static_cast<long int>(imag_nz_-1) / static_cast<long int>(j)
00247 >= static_cast<long int>(i))))
00248 {
00249 this->m_ = 0;
00250 this->n_ = 0;
00251 real_nz_ = 0;
00252 imag_nz_ = 0;
00253 real_ptr_ = NULL;
00254 imag_ptr_ = NULL;
00255 real_ind_ = NULL;
00256 imag_ind_ = NULL;
00257 this->real_data_ = NULL;
00258 this->imag_data_ = NULL;
00259 throw WrongDim(string("Matrix_ComplexSparse::")
00260 + string("Matrix_ComplexSparse(int, int, ")
00261 + string("const Vector&, const Vector&, const Vector&")
00262 + ", const Vector&, const Vector&, const Vector&)",
00263 string("There are more values (")
00264 + to_str(real_values.GetLength())
00265 + " values for the real part and "
00266 + to_str(real_values.GetLength()) + string(" values")
00267 + string(" for the imaginary part) than elements")
00268 + " in the matrix ("
00269 + to_str(i) + " by " + to_str(j) + ").");
00270 }
00271 #endif
00272
00273 this->real_ptr_ = real_ptr.GetData();
00274 this->imag_ptr_ = imag_ptr.GetData();
00275 this->real_ind_ = real_ind.GetData();
00276 this->imag_ind_ = imag_ind.GetData();
00277 this->real_data_ = real_values.GetData();
00278 this->imag_data_ = imag_values.GetData();
00279
00280 real_ptr.Nullify();
00281 imag_ptr.Nullify();
00282 real_ind.Nullify();
00283 imag_ind.Nullify();
00284 real_values.Nullify();
00285 imag_values.Nullify();
00286 }
00287
00288
00290 template <class T, class Prop, class Storage, class Allocator>
00291 Matrix_ComplexSparse<T, Prop, Storage, Allocator>
00292 ::Matrix_ComplexSparse(const Matrix_ComplexSparse<T, Prop,
00293 Storage, Allocator>& A)
00294 {
00295 this->m_ = 0;
00296 this->n_ = 0;
00297 real_nz_ = 0;
00298 imag_nz_ = 0;
00299 real_ptr_ = NULL;
00300 imag_ptr_ = NULL;
00301 real_ind_ = NULL;
00302 imag_ind_ = NULL;
00303 real_data_ = NULL;
00304 imag_data_ = NULL;
00305
00306 this->Copy(A);
00307 }
00308
00309
00310
00311
00312
00313
00314
00316 template <class T, class Prop, class Storage, class Allocator>
00317 inline Matrix_ComplexSparse<T, Prop, Storage, Allocator>
00318 ::~Matrix_ComplexSparse()
00319 {
00320 this->m_ = 0;
00321 this->n_ = 0;
00322
00323 #ifdef SELDON_CHECK_MEMORY
00324 try
00325 {
00326 #endif
00327
00328 if (real_ptr_ != NULL)
00329 {
00330 free(real_ptr_);
00331 real_ptr_ = NULL;
00332 }
00333
00334 #ifdef SELDON_CHECK_MEMORY
00335 }
00336 catch (...)
00337 {
00338 real_ptr_ = NULL;
00339 }
00340 #endif
00341
00342 #ifdef SELDON_CHECK_MEMORY
00343 try
00344 {
00345 #endif
00346
00347 if (imag_ptr_ != NULL)
00348 {
00349 free(imag_ptr_);
00350 imag_ptr_ = NULL;
00351 }
00352
00353 #ifdef SELDON_CHECK_MEMORY
00354 }
00355 catch (...)
00356 {
00357 imag_ptr_ = NULL;
00358 }
00359 #endif
00360
00361 #ifdef SELDON_CHECK_MEMORY
00362 try
00363 {
00364 #endif
00365
00366 if (real_ind_ != NULL)
00367 {
00368 free(real_ind_);
00369 real_ind_ = NULL;
00370 }
00371
00372 #ifdef SELDON_CHECK_MEMORY
00373 }
00374 catch (...)
00375 {
00376 real_ind_ = NULL;
00377 }
00378 #endif
00379
00380 #ifdef SELDON_CHECK_MEMORY
00381 try
00382 {
00383 #endif
00384
00385 if (imag_ind_ != NULL)
00386 {
00387 free(imag_ind_);
00388 imag_ind_ = NULL;
00389 }
00390
00391 #ifdef SELDON_CHECK_MEMORY
00392 }
00393 catch (...)
00394 {
00395 imag_ind_ = NULL;
00396 }
00397 #endif
00398
00399 #ifdef SELDON_CHECK_MEMORY
00400 try
00401 {
00402 #endif
00403
00404 if (this->real_data_ != NULL)
00405 {
00406 this->allocator_.deallocate(this->real_data_, real_nz_);
00407 this->real_data_ = NULL;
00408 }
00409
00410 #ifdef SELDON_CHECK_MEMORY
00411 }
00412 catch (...)
00413 {
00414 this->real_nz_ = 0;
00415 this->real_data_ = NULL;
00416 }
00417 #endif
00418
00419 #ifdef SELDON_CHECK_MEMORY
00420 try
00421 {
00422 #endif
00423
00424 if (this->imag_data_ != NULL)
00425 {
00426 this->allocator_.deallocate(this->imag_data_, imag_nz_);
00427 this->imag_data_ = NULL;
00428 }
00429
00430 #ifdef SELDON_CHECK_MEMORY
00431 }
00432 catch (...)
00433 {
00434 this->imag_nz_ = 0;
00435 this->imag_data_ = NULL;
00436 }
00437 #endif
00438
00439 this->real_nz_ = 0;
00440 this->imag_nz_ = 0;
00441 }
00442
00443
00445
00448 template <class T, class Prop, class Storage, class Allocator>
00449 inline void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::Clear()
00450 {
00451 this->~Matrix_ComplexSparse();
00452 }
00453
00454
00455
00456
00457
00458
00459
00461
00478 template <class T, class Prop, class Storage, class Allocator>
00479 template <class Storage0, class Allocator0,
00480 class Storage1, class Allocator1,
00481 class Storage2, class Allocator2>
00482 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::
00483 SetData(int i, int j,
00484 Vector<T, Storage0, Allocator0>& real_values,
00485 Vector<int, Storage1, Allocator1>& real_ptr,
00486 Vector<int, Storage2, Allocator2>& real_ind,
00487 Vector<T, Storage0, Allocator0>& imag_values,
00488 Vector<int, Storage1, Allocator1>& imag_ptr,
00489 Vector<int, Storage2, Allocator2>& imag_ind)
00490 {
00491 this->Clear();
00492 this->m_ = i;
00493 this->n_ = j;
00494 real_nz_ = real_values.GetLength();
00495 imag_nz_ = imag_values.GetLength();
00496
00497 #ifdef SELDON_CHECK_DIMENSIONS
00498
00499
00500 if (real_ind.GetLength() != real_nz_)
00501 {
00502 this->m_ = 0;
00503 this->n_ = 0;
00504 real_nz_ = 0;
00505 imag_nz_ = 0;
00506 real_ptr_ = NULL;
00507 imag_ptr_ = NULL;
00508 real_ind_ = NULL;
00509 imag_ind_ = NULL;
00510 this->real_data_ = NULL;
00511 this->imag_data_ = NULL;
00512 throw WrongDim(string("Matrix_ComplexSparse::SetData(int, int, ")
00513 + string("const Vector&, const Vector&, const Vector&")
00514 + ", const Vector&, const Vector&, const Vector&)",
00515 string("There are ") + to_str(real_nz_)
00516 + " values (real part) but "
00517 + to_str(real_ind.GetLength())
00518 + " row or column indices.");
00519 }
00520
00521 if (imag_ind.GetLength() != imag_nz_)
00522 {
00523 this->m_ = 0;
00524 this->n_ = 0;
00525 real_nz_ = 0;
00526 imag_nz_ = 0;
00527 real_ptr_ = NULL;
00528 imag_ptr_ = NULL;
00529 real_ind_ = NULL;
00530 imag_ind_ = NULL;
00531 this->real_data_ = NULL;
00532 this->imag_data_ = NULL;
00533 throw WrongDim(string("Matrix_ComplexSparse::SetData(int, int, ")
00534 + string("const Vector&, const Vector&, const Vector&")
00535 + ", const Vector&, const Vector&, const Vector&)",
00536 string("There are ") + to_str(imag_nz_)
00537 + " values (imaginary part) but "
00538 + to_str(imag_ind.GetLength())
00539 + " row or column indices.");
00540 }
00541
00542 if (real_ptr.GetLength()-1 != Storage::GetFirst(i, j))
00543 {
00544 this->m_ = 0;
00545 this->n_ = 0;
00546 real_nz_ = 0;
00547 imag_nz_ = 0;
00548 real_ptr_ = NULL;
00549 imag_ptr_ = NULL;
00550 real_ind_ = NULL;
00551 imag_ind_ = NULL;
00552 this->real_data_ = NULL;
00553 this->imag_data_ = NULL;
00554 throw WrongDim(string("Matrix_ComplexSparse::SetData(int, int, ")
00555 + string("const Vector&, const Vector&, const Vector&")
00556 + ", const Vector&, const Vector&, const Vector&)",
00557 string("The vector of start indices (real part)")
00558 + " contains " + to_str(real_ptr.GetLength()-1)
00559 + string(" row or column start indices (plus the")
00560 + " number of non-zero entries) but there are "
00561 + to_str(Storage::GetFirst(i, j))
00562 + " rows or columns ("
00563 + to_str(i) + " by " + to_str(j) + " matrix).");
00564 }
00565
00566 if (imag_ptr.GetLength()-1 != Storage::GetFirst(i, j))
00567 {
00568 this->m_ = 0;
00569 this->n_ = 0;
00570 real_nz_ = 0;
00571 imag_nz_ = 0;
00572 real_ptr_ = NULL;
00573 imag_ptr_ = NULL;
00574 real_ind_ = NULL;
00575 imag_ind_ = NULL;
00576 this->real_data_ = NULL;
00577 this->imag_data_ = NULL;
00578 throw WrongDim(string("Matrix_ComplexSparse::SetData(int, int, ")
00579 + string("const Vector&, const Vector&, const Vector&")
00580 + ", const Vector&, const Vector&, const Vector&)",
00581 string("The vector of start indices (imaginary part)")
00582 + " contains " + to_str(imag_ptr.GetLength()-1)
00583 + string(" row or column start indices (plus the")
00584 + " number of non-zero entries) but there are "
00585 + to_str(Storage::GetFirst(i, j))
00586 + " rows or columns ("
00587 + to_str(i) + " by " + to_str(j) + " matrix).");
00588 }
00589
00590 if ((real_nz_ > 0
00591 && (j == 0
00592 || static_cast<long int>(real_nz_-1) / static_cast<long int>(j)
00593 >= static_cast<long int>(i)))
00594 ||
00595 (imag_nz_ > 0
00596 && (j == 0
00597 || static_cast<long int>(imag_nz_-1) / static_cast<long int>(j)
00598 >= static_cast<long int>(i))))
00599 {
00600 this->m_ = 0;
00601 this->n_ = 0;
00602 real_nz_ = 0;
00603 imag_nz_ = 0;
00604 real_ptr_ = NULL;
00605 imag_ptr_ = NULL;
00606 real_ind_ = NULL;
00607 imag_ind_ = NULL;
00608 this->real_data_ = NULL;
00609 this->imag_data_ = NULL;
00610 throw WrongDim(string("Matrix_ComplexSparse::SetData(int, int, ")
00611 + string("const Vector&, const Vector&, const Vector&")
00612 + ", const Vector&, const Vector&, const Vector&)",
00613 string("There are more values (")
00614 + to_str(real_values.GetLength())
00615 + " values for the real part and "
00616 + to_str(real_values.GetLength()) + string(" values")
00617 + string(" for the imaginary part) than elements")
00618 + " in the matrix ("
00619 + to_str(i) + " by " + to_str(j) + ").");
00620 }
00621 #endif
00622
00623 this->real_ptr_ = real_ptr.GetData();
00624 this->imag_ptr_ = imag_ptr.GetData();
00625 this->real_ind_ = real_ind.GetData();
00626 this->imag_ind_ = imag_ind.GetData();
00627 this->real_data_ = real_values.GetData();
00628 this->imag_data_ = imag_values.GetData();
00629
00630 real_ptr.Nullify();
00631 imag_ptr.Nullify();
00632 real_ind.Nullify();
00633 imag_ind.Nullify();
00634 real_values.Nullify();
00635 imag_values.Nullify();
00636 }
00637
00638
00640
00661 template <class T, class Prop, class Storage, class Allocator>
00662 inline void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::
00663 SetData(int i, int j, int real_nz,
00664 typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
00665 ::pointer real_values,
00666 int* real_ptr, int* real_ind, int imag_nz,
00667 typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
00668 ::pointer imag_values,
00669 int* imag_ptr, int* imag_ind)
00670 {
00671 this->Clear();
00672
00673 this->m_ = i;
00674 this->n_ = j;
00675
00676 this->real_nz_ = real_nz;
00677 this->imag_nz_ = imag_nz;
00678
00679 real_data_ = real_values;
00680 imag_data_ = imag_values;
00681 real_ind_ = real_ind;
00682 imag_ind_ = imag_ind;
00683 real_ptr_ = real_ptr;
00684 imag_ptr_ = imag_ptr;
00685 }
00686
00687
00689
00693 template <class T, class Prop, class Storage, class Allocator>
00694 inline void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::Nullify()
00695 {
00696 this->data_ = NULL;
00697 this->m_ = 0;
00698 this->n_ = 0;
00699 real_nz_ = 0;
00700 real_ptr_ = NULL;
00701 real_ind_ = NULL;
00702 imag_nz_ = 0;
00703 imag_ptr_ = NULL;
00704 imag_ind_ = NULL;
00705 real_data_ = NULL;
00706 imag_data_ = NULL;
00707 }
00708
00709
00711 template <class T, class Prop, class Storage, class Allocator>
00712 inline void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
00713 ::Reallocate(int i, int j)
00714 {
00715
00716 Clear();
00717
00718 this->m_ = i;
00719 this->n_ = j;
00720
00721
00722 #ifdef SELDON_CHECK_MEMORY
00723 try
00724 {
00725 #endif
00726
00727 real_ptr_ = reinterpret_cast<int*>( calloc(Storage::GetFirst(i, j)+1,
00728 sizeof(int)) );
00729
00730 imag_ptr_ = reinterpret_cast<int*>( calloc(Storage::GetFirst(i, j)+1,
00731 sizeof(int)) );
00732
00733 #ifdef SELDON_CHECK_MEMORY
00734 }
00735 catch (...)
00736 {
00737 this->m_ = 0;
00738 this->n_ = 0;
00739 real_nz_ = 0;
00740 real_ptr_ = NULL;
00741 real_ind_ = NULL;
00742 imag_nz_ = 0;
00743 imag_ptr_ = NULL;
00744 imag_ind_ = NULL;
00745 real_data_ = NULL;
00746 imag_data_ = NULL;
00747 }
00748 if ((real_ptr_ == NULL) || (imag_ptr_ == NULL))
00749 {
00750 this->m_ = 0;
00751 this->n_ = 0;
00752 real_nz_ = 0;
00753 real_ptr_ = NULL;
00754 real_ind_ = NULL;
00755 imag_nz_ = 0;
00756 imag_ptr_ = NULL;
00757 imag_ind_ = NULL;
00758 real_data_ = NULL;
00759 imag_data_ = NULL;
00760 }
00761 if (((real_ptr_ == NULL) || (imag_ptr_ == NULL)) && i != 0 && j != 0)
00762 throw NoMemory("Matrix_ComplexSparse::Reallocate(int, int)",
00763 string("Unable to allocate ")
00764 + to_str(sizeof(int) * (Storage::GetFirst(i, j)+1) )
00765 + " bytes to store " + to_str(Storage::GetFirst(i, j)+1)
00766 + " row or column start indices, for a "
00767 + to_str(i) + " by " + to_str(j) + " matrix.");
00768 #endif
00769
00770
00771 for (int k = 0; k <= Storage::GetFirst(i, j); k++)
00772 {
00773 real_ptr_[k] = 0;
00774 imag_ptr_[k] = 0;
00775 }
00776 }
00777
00778
00780 template <class T, class Prop, class Storage, class Allocator>
00781 inline void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
00782 ::Reallocate(int i, int j, int real_nz, int imag_nz)
00783 {
00784
00785 Clear();
00786
00787 this->m_ = i;
00788 this->n_ = j;
00789 this->real_nz_ = real_nz;
00790 this->imag_nz_ = imag_nz;
00791
00792 #ifdef SELDON_CHECK_DIMENSIONS
00793 if (real_nz_ < 0 || imag_nz_ < 0)
00794 {
00795 this->m_ = 0;
00796 this->n_ = 0;
00797 real_nz_ = 0;
00798 imag_nz_ = 0;
00799 real_ptr_ = NULL;
00800 imag_ptr_ = NULL;
00801 real_ind_ = NULL;
00802 imag_ind_ = NULL;
00803 this->real_data_ = NULL;
00804 this->imag_data_ = NULL;
00805 throw WrongDim(string("Matrix_ComplexSparse::")
00806 + "Matrix_ComplexSparse(int, int, int, int)",
00807 "Invalid number of non-zero elements: "
00808 + to_str(real_nz) + " in the real part and "
00809 + to_str(imag_nz) + " in the imaginary part.");
00810 }
00811 if ((real_nz_ > 0
00812 && (j == 0
00813 || static_cast<long int>(real_nz_-1) / static_cast<long int>(j)
00814 >= static_cast<long int>(i)))
00815 ||
00816 (imag_nz_ > 0
00817 && (j == 0
00818 || static_cast<long int>(imag_nz_-1) / static_cast<long int>(j)
00819 >= static_cast<long int>(i))))
00820 {
00821 this->m_ = 0;
00822 this->n_ = 0;
00823 real_nz_ = 0;
00824 imag_nz_ = 0;
00825 real_ptr_ = NULL;
00826 imag_ptr_ = NULL;
00827 real_ind_ = NULL;
00828 imag_ind_ = NULL;
00829 this->real_data_ = NULL;
00830 this->imag_data_ = NULL;
00831 throw WrongDim(string("Matrix_ComplexSparse::")
00832 + "Matrix_ComplexSparse(int, int, int, int)",
00833 string("There are more values (") + to_str(real_nz)
00834 + " values for the real part and " + to_str(imag_nz)
00835 + string(" values for the imaginary part) than")
00836 + " elements in the matrix (" + to_str(i) + " by "
00837 + to_str(j) + ").");
00838 }
00839 #endif
00840
00841 #ifdef SELDON_CHECK_MEMORY
00842 try
00843 {
00844 #endif
00845
00846 real_ptr_ = reinterpret_cast<int*>( calloc(Storage::GetFirst(i, j)+1,
00847 sizeof(int)) );
00848
00849 #ifdef SELDON_CHECK_MEMORY
00850 }
00851 catch (...)
00852 {
00853 this->m_ = 0;
00854 this->n_ = 0;
00855 real_nz_ = 0;
00856 imag_nz_ = 0;
00857 real_ptr_ = NULL;
00858 imag_ptr_ = NULL;
00859 real_ind_ = NULL;
00860 imag_ind_ = NULL;
00861 this->real_data_ = NULL;
00862 this->imag_data_ = NULL;
00863 }
00864 if (real_ptr_ == NULL)
00865 {
00866 this->m_ = 0;
00867 this->n_ = 0;
00868 real_nz_ = 0;
00869 imag_nz_ = 0;
00870 imag_ptr_ = 0;
00871 real_ind_ = NULL;
00872 imag_ind_ = NULL;
00873 this->real_data_ = NULL;
00874 this->imag_data_ = NULL;
00875 }
00876 if (real_ptr_ == NULL && i != 0 && j != 0)
00877 throw NoMemory(string("Matrix_ComplexSparse::")
00878 + "Matrix_ComplexSparse(int, int, int, int)",
00879 string("Unable to allocate ")
00880 + to_str(sizeof(int) * (Storage::GetFirst(i, j)+1))
00881 + " bytes to store " + to_str(Storage::GetFirst(i, j)+1)
00882 + string(" row or column start indices (for the real")
00883 + " part), for a "
00884 + to_str(i) + " by " + to_str(j) + " matrix.");
00885 #endif
00886
00887 #ifdef SELDON_CHECK_MEMORY
00888 try
00889 {
00890 #endif
00891
00892 imag_ptr_ = reinterpret_cast<int*>( calloc(Storage::GetFirst(i, j)+1,
00893 sizeof(int)) );
00894
00895 #ifdef SELDON_CHECK_MEMORY
00896 }
00897 catch (...)
00898 {
00899 this->m_ = 0;
00900 this->n_ = 0;
00901 real_nz_ = 0;
00902 imag_nz_ = 0;
00903 free(real_ptr_);
00904 real_ptr_ = NULL;
00905 imag_ptr_ = NULL;
00906 real_ind_ = NULL;
00907 imag_ind_ = NULL;
00908 this->real_data_ = NULL;
00909 this->imag_data_ = NULL;
00910 }
00911 if (imag_ptr_ == NULL)
00912 {
00913 this->m_ = 0;
00914 this->n_ = 0;
00915 real_nz_ = 0;
00916 imag_nz_ = 0;
00917 free(real_ptr_);
00918 real_ptr_ = 0;
00919 real_ind_ = NULL;
00920 imag_ind_ = NULL;
00921 this->real_data_ = NULL;
00922 this->imag_data_ = NULL;
00923 }
00924 if (imag_ptr_ == NULL && i != 0 && j != 0)
00925 throw NoMemory(string("Matrix_ComplexSparse::")
00926 + "Matrix_ComplexSparse(int, int, int, int)",
00927 string("Unable to allocate ")
00928 + to_str(sizeof(int) * (Storage::GetFirst(i, j)+1))
00929 + " bytes to store " + to_str(Storage::GetFirst(i, j)+1)
00930 + string(" row or column start indices (for the")
00931 + string(" imaginary part), for a ")
00932 + to_str(i) + " by " + to_str(j) + " matrix.");
00933 #endif
00934
00935 #ifdef SELDON_CHECK_MEMORY
00936 try
00937 {
00938 #endif
00939
00940 real_ind_ = reinterpret_cast<int*>( calloc(real_nz_, sizeof(int)) );
00941
00942 #ifdef SELDON_CHECK_MEMORY
00943 }
00944 catch (...)
00945 {
00946 this->m_ = 0;
00947 this->n_ = 0;
00948 real_nz_ = 0;
00949 imag_nz_ = 0;
00950 free(real_ptr_);
00951 free(imag_ptr_);
00952 real_ptr_ = NULL;
00953 imag_ptr_ = NULL;
00954 real_ind_ = NULL;
00955 imag_ind_ = NULL;
00956 this->real_data_ = NULL;
00957 this->imag_data_ = NULL;
00958 }
00959 if (real_ind_ == NULL)
00960 {
00961 this->m_ = 0;
00962 this->n_ = 0;
00963 real_nz_ = 0;
00964 imag_nz_ = 0;
00965 free(real_ptr_);
00966 free(imag_ptr_);
00967 real_ptr_ = NULL;
00968 imag_ptr_ = NULL;
00969 this->real_data_ = NULL;
00970 this->imag_data_ = NULL;
00971 }
00972 if (real_ind_ == NULL && i != 0 && j != 0)
00973 throw NoMemory(string("Matrix_ComplexSparse::")
00974 + "Matrix_ComplexSparse(int, int, int, int)",
00975 string("Unable to allocate ")
00976 + to_str(sizeof(int) * real_nz)
00977 + " bytes to store " + to_str(real_nz)
00978 + " row or column indices (real part), for a "
00979 + to_str(i) + " by " + to_str(j) + " matrix.");
00980 #endif
00981
00982 #ifdef SELDON_CHECK_MEMORY
00983 try
00984 {
00985 #endif
00986
00987 imag_ind_ = reinterpret_cast<int*>( calloc(imag_nz_, sizeof(int)) );
00988
00989 #ifdef SELDON_CHECK_MEMORY
00990 }
00991 catch (...)
00992 {
00993 this->m_ = 0;
00994 this->n_ = 0;
00995 real_nz_ = 0;
00996 imag_nz_ = 0;
00997 free(real_ptr_);
00998 free(imag_ptr_);
00999 real_ptr_ = NULL;
01000 imag_ptr_ = NULL;
01001 free(imag_ind_);
01002 real_ind_ = NULL;
01003 imag_ind_ = NULL;
01004 this->real_data_ = NULL;
01005 this->imag_data_ = NULL;
01006 }
01007 if (real_ind_ == NULL)
01008 {
01009 this->m_ = 0;
01010 this->n_ = 0;
01011 real_nz_ = 0;
01012 imag_nz_ = 0;
01013 free(real_ptr_);
01014 free(imag_ptr_);
01015 real_ptr_ = NULL;
01016 imag_ptr_ = NULL;
01017 free(imag_ind_);
01018 imag_ind_ = NULL;
01019 this->real_data_ = NULL;
01020 this->imag_data_ = NULL;
01021 }
01022 if (imag_ind_ == NULL && i != 0 && j != 0)
01023 throw NoMemory(string("Matrix_ComplexSparse::")
01024 + "Matrix_ComplexSparse(int, int, int, int)",
01025 string("Unable to allocate ")
01026 + to_str(sizeof(int) * imag_nz)
01027 + " bytes to store " + to_str(imag_nz)
01028 + " row or column indices (imaginary part), for a "
01029 + to_str(i) + " by " + to_str(j) + " matrix.");
01030 #endif
01031
01032 #ifdef SELDON_CHECK_MEMORY
01033 try
01034 {
01035 #endif
01036
01037 this->real_data_ = this->allocator_.allocate(real_nz_, this);
01038
01039 #ifdef SELDON_CHECK_MEMORY
01040 }
01041 catch (...)
01042 {
01043 this->m_ = 0;
01044 this->n_ = 0;
01045 free(real_ptr_);
01046 free(imag_ptr_);
01047 real_ptr_ = NULL;
01048 imag_ptr_ = NULL;
01049 free(real_ind_);
01050 free(imag_ind_);
01051 real_ind_ = NULL;
01052 imag_ind_ = NULL;
01053 this->real_data_ = NULL;
01054 this->imag_data_ = NULL;
01055 }
01056 if (real_data_ == NULL)
01057 {
01058 this->m_ = 0;
01059 this->n_ = 0;
01060 free(real_ptr_);
01061 free(imag_ptr_);
01062 real_ptr_ = NULL;
01063 imag_ptr_ = NULL;
01064 free(real_ind_);
01065 free(imag_ind_);
01066 real_ind_ = NULL;
01067 imag_ind_ = NULL;
01068 imag_data_ = NULL;
01069 }
01070 if (real_data_ == NULL && i != 0 && j != 0)
01071 throw NoMemory(string("Matrix_ComplexSparse::")
01072 + "Matrix_ComplexSparse(int, int, int, int)",
01073 string("Unable to allocate ")
01074 + to_str(sizeof(int) * real_nz)
01075 + " bytes to store " + to_str(real_nz)
01076 + " values (real part), for a "
01077 + to_str(i) + " by " + to_str(j) + " matrix.");
01078 #endif
01079
01080 #ifdef SELDON_CHECK_MEMORY
01081 try
01082 {
01083 #endif
01084
01085 this->imag_data_ = this->allocator_.allocate(imag_nz_, this);
01086
01087 #ifdef SELDON_CHECK_MEMORY
01088 }
01089 catch (...)
01090 {
01091 this->m_ = 0;
01092 this->n_ = 0;
01093 free(real_ptr_);
01094 free(imag_ptr_);
01095 real_ptr_ = NULL;
01096 imag_ptr_ = NULL;
01097 free(real_ind_);
01098 free(imag_ind_);
01099 real_ind_ = NULL;
01100 imag_ind_ = NULL;
01101 this->allocator_.deallocate(this->real_data_, real_nz_);
01102 this->real_data_ = NULL;
01103 this->imag_data_ = NULL;
01104 }
01105 if (real_data_ == NULL)
01106 {
01107 this->m_ = 0;
01108 this->n_ = 0;
01109 free(real_ptr_);
01110 free(imag_ptr_);
01111 real_ptr_ = NULL;
01112 imag_ptr_ = NULL;
01113 free(real_ind_);
01114 free(imag_ind_);
01115 real_ind_ = NULL;
01116 imag_ind_ = NULL;
01117 this->allocator_.deallocate(this->real_data_, real_nz_);
01118 real_data_ = NULL;
01119 }
01120 if (imag_data_ == NULL && i != 0 && j != 0)
01121 throw NoMemory(string("Matrix_ComplexSparse::")
01122 + "Matrix_ComplexSparse(int, int, int, int)",
01123 string("Unable to allocate ")
01124 + to_str(sizeof(int) * imag_nz)
01125 + " bytes to store " + to_str(imag_nz)
01126 + " values (imaginary part), for a "
01127 + to_str(i) + " by " + to_str(j) + " matrix.");
01128 #endif
01129 }
01130
01131
01133
01137 template <class T, class Prop, class Storage, class Allocator>
01138 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
01139 ::Resize(int i, int j)
01140 {
01141 if (Storage::GetFirst(i, j) < Storage::GetFirst(this->m_, this->n_))
01142 Resize(i, j, real_ptr_[Storage::GetFirst(i, j)],
01143 imag_ptr_[Storage::GetFirst(i, j)]);
01144 else
01145 Resize(i, j, real_nz_, imag_nz_);
01146 }
01147
01148
01150
01155 template <class T, class Prop, class Storage, class Allocator>
01156 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
01157 ::Resize(int i, int j, int real_nz, int imag_nz)
01158 {
01159 #ifdef SELDON_CHECK_DIMENSIONS
01160 if (real_nz < 0 || imag_nz < 0)
01161 {
01162 this->m_ = 0;
01163 this->n_ = 0;
01164 real_nz_ = 0;
01165 imag_nz_ = 0;
01166 real_ptr_ = NULL;
01167 imag_ptr_ = NULL;
01168 real_ind_ = NULL;
01169 imag_ind_ = NULL;
01170 this->real_data_ = NULL;
01171 this->imag_data_ = NULL;
01172 throw WrongDim(string("Matrix_ComplexSparse::")
01173 + "Resize(int, int, int, int)",
01174 "Invalid number of non-zero elements: "
01175 + to_str(real_nz) + " in the real part and "
01176 + to_str(imag_nz) + " in the imaginary part.");
01177 }
01178 if ((real_nz > 0
01179 && (j == 0
01180 || static_cast<long int>(real_nz-1) / static_cast<long int>(j)
01181 >= static_cast<long int>(i)))
01182 ||
01183 (imag_nz > 0
01184 && (j == 0
01185 || static_cast<long int>(imag_nz-1) / static_cast<long int>(j)
01186 >= static_cast<long int>(i))))
01187 {
01188 this->m_ = 0;
01189 this->n_ = 0;
01190 real_nz_ = 0;
01191 imag_nz_ = 0;
01192 real_ptr_ = NULL;
01193 imag_ptr_ = NULL;
01194 real_ind_ = NULL;
01195 imag_ind_ = NULL;
01196 this->real_data_ = NULL;
01197 this->imag_data_ = NULL;
01198 throw WrongDim(string("Matrix_ComplexSparse::")
01199 + "Resize(int, int, int, int)",
01200 string("There are more values (") + to_str(real_nz)
01201 + " values for the real part and " + to_str(imag_nz)
01202 + string(" values for the imaginary part) than")
01203 + " elements in the matrix (" + to_str(i) + " by "
01204 + to_str(j) + ").");
01205 }
01206 #endif
01207
01208 if (Storage::GetFirst(this->m_, this->n_) != Storage::GetFirst(i, j))
01209 {
01210 #ifdef SELDON_CHECK_MEMORY
01211 try
01212 {
01213 #endif
01214
01215 real_ptr_
01216 = reinterpret_cast<int*>( realloc(real_ptr_,
01217 (Storage::GetFirst(i, j)+1)
01218 *sizeof(int)) );
01219
01220 #ifdef SELDON_CHECK_MEMORY
01221 }
01222 catch (...)
01223 {
01224 this->m_ = 0;
01225 this->n_ = 0;
01226 real_nz_ = 0;
01227 imag_nz_ = 0;
01228 real_ptr_ = NULL;
01229 imag_ptr_ = NULL;
01230 real_ind_ = NULL;
01231 imag_ind_ = NULL;
01232 this->real_data_ = NULL;
01233 this->imag_data_ = NULL;
01234 }
01235 if (real_ptr_ == NULL)
01236 {
01237 this->m_ = 0;
01238 this->n_ = 0;
01239 real_nz_ = 0;
01240 imag_nz_ = 0;
01241 imag_ptr_ = 0;
01242 real_ind_ = NULL;
01243 imag_ind_ = NULL;
01244 this->real_data_ = NULL;
01245 this->imag_data_ = NULL;
01246 }
01247 if (real_ptr_ == NULL && i != 0 && j != 0)
01248 throw NoMemory(string("Matrix_ComplexSparse::")
01249 + "Resize(int, int, int, int)",
01250 string("Unable to allocate ")
01251 + to_str(sizeof(int) * (Storage::GetFirst(i, j)+1))
01252 + " bytes to store " + to_str(Storage::GetFirst(i, j)+1)
01253 + string(" row or column start indices (for the real")
01254 + " part), for a "
01255 + to_str(i) + " by " + to_str(j) + " matrix.");
01256 #endif
01257
01258 #ifdef SELDON_CHECK_MEMORY
01259 try
01260 {
01261 #endif
01262
01263 imag_ptr_
01264 = reinterpret_cast<int*>( realloc(imag_ptr_,
01265 (Storage::GetFirst(i, j)+1)
01266 *sizeof(int)) );
01267
01268 #ifdef SELDON_CHECK_MEMORY
01269 }
01270 catch (...)
01271 {
01272 this->m_ = 0;
01273 this->n_ = 0;
01274 real_nz_ = 0;
01275 imag_nz_ = 0;
01276 free(real_ptr_);
01277 real_ptr_ = NULL;
01278 imag_ptr_ = NULL;
01279 real_ind_ = NULL;
01280 imag_ind_ = NULL;
01281 this->real_data_ = NULL;
01282 this->imag_data_ = NULL;
01283 }
01284 if (imag_ptr_ == NULL)
01285 {
01286 this->m_ = 0;
01287 this->n_ = 0;
01288 real_nz_ = 0;
01289 imag_nz_ = 0;
01290 free(real_ptr_);
01291 real_ptr_ = 0;
01292 real_ind_ = NULL;
01293 imag_ind_ = NULL;
01294 this->real_data_ = NULL;
01295 this->imag_data_ = NULL;
01296 }
01297 if (imag_ptr_ == NULL && i != 0 && j != 0)
01298 throw NoMemory(string("Matrix_ComplexSparse::")
01299 + "Resize(int, int, int, int)",
01300 string("Unable to allocate ")
01301 + to_str(sizeof(int) * (Storage::GetFirst(i, j)+1))
01302 + " bytes to store " + to_str(Storage::GetFirst(i, j)+1)
01303 + string(" row or column start indices (for the")
01304 + string(" imaginary part), for a ")
01305 + to_str(i) + " by " + to_str(j) + " matrix.");
01306 #endif
01307 }
01308
01309 if (real_nz != real_nz_)
01310 {
01311 #ifdef SELDON_CHECK_MEMORY
01312 try
01313 {
01314 #endif
01315
01316 real_ind_
01317 = reinterpret_cast<int*>( realloc(real_ind_,
01318 real_nz*sizeof(int)) );
01319
01320 #ifdef SELDON_CHECK_MEMORY
01321 }
01322 catch (...)
01323 {
01324 this->m_ = 0;
01325 this->n_ = 0;
01326 real_nz_ = 0;
01327 imag_nz_ = 0;
01328 free(real_ptr_);
01329 free(imag_ptr_);
01330 real_ptr_ = NULL;
01331 imag_ptr_ = NULL;
01332 real_ind_ = NULL;
01333 imag_ind_ = NULL;
01334 this->real_data_ = NULL;
01335 this->imag_data_ = NULL;
01336 }
01337 if (real_ind_ == NULL)
01338 {
01339 this->m_ = 0;
01340 this->n_ = 0;
01341 real_nz_ = 0;
01342 imag_nz_ = 0;
01343 free(real_ptr_);
01344 free(imag_ptr_);
01345 real_ptr_ = NULL;
01346 imag_ptr_ = NULL;
01347 this->real_data_ = NULL;
01348 this->imag_data_ = NULL;
01349 }
01350 if (real_ind_ == NULL && i != 0 && j != 0)
01351 throw NoMemory(string("Matrix_ComplexSparse::")
01352 + "Resize(int, int, int, int)",
01353 string("Unable to allocate ")
01354 + to_str(sizeof(int) * real_nz)
01355 + " bytes to store " + to_str(real_nz)
01356 + " row or column indices (real part), for a "
01357 + to_str(i) + " by " + to_str(j) + " matrix.");
01358 #endif
01359 }
01360
01361 if (imag_nz != imag_nz_)
01362 {
01363 #ifdef SELDON_CHECK_MEMORY
01364 try
01365 {
01366 #endif
01367
01368 imag_ind_
01369 = reinterpret_cast<int*>( realloc(imag_ind_,
01370 imag_nz*sizeof(int)) );
01371
01372 #ifdef SELDON_CHECK_MEMORY
01373 }
01374 catch (...)
01375 {
01376 this->m_ = 0;
01377 this->n_ = 0;
01378 real_nz_ = 0;
01379 imag_nz_ = 0;
01380 free(real_ptr_);
01381 free(imag_ptr_);
01382 real_ptr_ = NULL;
01383 imag_ptr_ = NULL;
01384 free(imag_ind_);
01385 real_ind_ = NULL;
01386 imag_ind_ = NULL;
01387 this->real_data_ = NULL;
01388 this->imag_data_ = NULL;
01389 }
01390 if (real_ind_ == NULL)
01391 {
01392 this->m_ = 0;
01393 this->n_ = 0;
01394 real_nz_ = 0;
01395 imag_nz_ = 0;
01396 free(real_ptr_);
01397 free(imag_ptr_);
01398 real_ptr_ = NULL;
01399 imag_ptr_ = NULL;
01400 free(imag_ind_);
01401 imag_ind_ = NULL;
01402 this->real_data_ = NULL;
01403 this->imag_data_ = NULL;
01404 }
01405 if (imag_ind_ == NULL && i != 0 && j != 0)
01406 throw NoMemory(string("Matrix_ComplexSparse::")
01407 + "Resize(int, int, int, int)",
01408 string("Unable to allocate ")
01409 + to_str(sizeof(int) * imag_nz)
01410 + " bytes to store " + to_str(imag_nz)
01411 + " row or column indices (imaginary part), for a "
01412 + to_str(i) + " by " + to_str(j) + " matrix.");
01413 #endif
01414 }
01415
01416 if (real_nz != real_nz_)
01417 {
01418 Vector<T, VectFull, Allocator> val;
01419 val.SetData(real_nz_, real_data_);
01420 val.Resize(real_nz);
01421
01422 real_data_ = val.GetData();
01423 val.Nullify();
01424 }
01425
01426 if (imag_nz != imag_nz_)
01427 {
01428 Vector<T, VectFull, Allocator> val;
01429 val.SetData(imag_nz_, imag_data_);
01430 val.Resize(imag_nz);
01431
01432 imag_data_ = val.GetData();
01433 val.Nullify();
01434 }
01435
01436
01437 for (int k = Storage::GetFirst(this->m_, this->n_);
01438 k <= Storage::GetFirst(i, j); k++)
01439 {
01440 real_ptr_[k] = real_nz_;
01441 imag_ptr_[k] = imag_nz_;
01442 }
01443
01444 this->m_ = i;
01445 this->n_ = j;
01446 real_nz_ = real_nz;
01447 imag_nz_ = imag_nz;
01448 }
01449
01450
01452 template <class T, class Prop, class Storage, class Allocator>
01453 inline void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::
01454 Copy(const Matrix_ComplexSparse<T, Prop, Storage, Allocator>& A)
01455 {
01456 this->Clear();
01457 int i = A.m_;
01458 int j = A.n_;
01459 real_nz_ = A.real_nz_;
01460 imag_nz_ = A.imag_nz_;
01461 this->m_ = i;
01462 this->n_ = j;
01463 if ((i == 0)||(j == 0))
01464 {
01465 this->m_ = 0;
01466 this->n_ = 0;
01467 this->real_nz_ = 0;
01468 this->imag_nz_ = 0;
01469 return;
01470 }
01471
01472 #ifdef SELDON_CHECK_DIMENSIONS
01473 if ((real_nz_ > 0
01474 && (j == 0
01475 || static_cast<long int>(real_nz_-1) / static_cast<long int>(j)
01476 >= static_cast<long int>(i)))
01477 ||
01478 (imag_nz_ > 0
01479 && (j == 0
01480 || static_cast<long int>(imag_nz_-1) / static_cast<long int>(j)
01481 >= static_cast<long int>(i))))
01482 {
01483 this->m_ = 0;
01484 this->n_ = 0;
01485 real_nz_ = 0;
01486 imag_nz_ = 0;
01487 real_ptr_ = NULL;
01488 imag_ptr_ = NULL;
01489 real_ind_ = NULL;
01490 imag_ind_ = NULL;
01491 this->real_data_ = NULL;
01492 this->imag_data_ = NULL;
01493 throw WrongDim(string("Matrix_ComplexSparse::")
01494 + "Matrix_ComplexSparse(int, int, int, int)",
01495 string("There are more values (") + to_str(real_nz_)
01496 + " values for the real part and " + to_str(imag_nz_)
01497 + string(" values for the imaginary part) than")
01498 + " elements in the matrix (" + to_str(i) + " by "
01499 + to_str(j) + ").");
01500 }
01501 #endif
01502
01503 #ifdef SELDON_CHECK_MEMORY
01504 try
01505 {
01506 #endif
01507
01508 real_ptr_ = reinterpret_cast<int*>( calloc(Storage::GetFirst(i, j)+1,
01509 sizeof(int)) );
01510 memcpy(this->real_ptr_, A.real_ptr_,
01511 (Storage::GetFirst(i, j) + 1) * sizeof(int));
01512
01513 #ifdef SELDON_CHECK_MEMORY
01514 }
01515 catch (...)
01516 {
01517 this->m_ = 0;
01518 this->n_ = 0;
01519 real_nz_ = 0;
01520 imag_nz_ = 0;
01521 real_ptr_ = NULL;
01522 imag_ptr_ = NULL;
01523 real_ind_ = NULL;
01524 imag_ind_ = NULL;
01525 this->real_data_ = NULL;
01526 this->imag_data_ = NULL;
01527 }
01528 if (real_ptr_ == NULL)
01529 {
01530 this->m_ = 0;
01531 this->n_ = 0;
01532 real_nz_ = 0;
01533 imag_nz_ = 0;
01534 imag_ptr_ = 0;
01535 real_ind_ = NULL;
01536 imag_ind_ = NULL;
01537 this->real_data_ = NULL;
01538 this->imag_data_ = NULL;
01539 }
01540 if (real_ptr_ == NULL && i != 0 && j != 0)
01541 throw NoMemory(string("Matrix_ComplexSparse::")
01542 + "Matrix_ComplexSparse(int, int, int, int)",
01543 string("Unable to allocate ")
01544 + to_str(sizeof(int) * (Storage::GetFirst(i, j)+1))
01545 + " bytes to store " + to_str(Storage::GetFirst(i, j)+1)
01546 + string(" row or column start indices (for the real")
01547 + " part), for a "
01548 + to_str(i) + " by " + to_str(j) + " matrix.");
01549 #endif
01550
01551 #ifdef SELDON_CHECK_MEMORY
01552 try
01553 {
01554 #endif
01555
01556 imag_ptr_ = reinterpret_cast<int*>( calloc(Storage::GetFirst(i, j)+1,
01557 sizeof(int)) );
01558 memcpy(this->imag_ptr_, A.imag_ptr_,
01559 (Storage::GetFirst(i, j) + 1) * sizeof(int));
01560
01561 #ifdef SELDON_CHECK_MEMORY
01562 }
01563 catch (...)
01564 {
01565 this->m_ = 0;
01566 this->n_ = 0;
01567 real_nz_ = 0;
01568 imag_nz_ = 0;
01569 free(real_ptr_);
01570 real_ptr_ = NULL;
01571 imag_ptr_ = NULL;
01572 real_ind_ = NULL;
01573 imag_ind_ = NULL;
01574 this->real_data_ = NULL;
01575 this->imag_data_ = NULL;
01576 }
01577 if (imag_ptr_ == NULL)
01578 {
01579 this->m_ = 0;
01580 this->n_ = 0;
01581 real_nz_ = 0;
01582 imag_nz_ = 0;
01583 free(real_ptr_);
01584 real_ptr_ = 0;
01585 real_ind_ = NULL;
01586 imag_ind_ = NULL;
01587 this->real_data_ = NULL;
01588 this->imag_data_ = NULL;
01589 }
01590 if (imag_ptr_ == NULL && i != 0 && j != 0)
01591 throw NoMemory(string("Matrix_ComplexSparse::")
01592 + "Matrix_ComplexSparse(int, int, int, int)",
01593 string("Unable to allocate ")
01594 + to_str(sizeof(int) * (Storage::GetFirst(i, j)+1))
01595 + " bytes to store " + to_str(Storage::GetFirst(i, j)+1)
01596 + string(" row or column start indices (for the")
01597 + string(" imaginary part), for a ")
01598 + to_str(i) + " by " + to_str(j) + " matrix.");
01599 #endif
01600
01601 #ifdef SELDON_CHECK_MEMORY
01602 try
01603 {
01604 #endif
01605
01606 real_ind_ = reinterpret_cast<int*>( calloc(real_nz_, sizeof(int)) );
01607 memcpy(this->real_ind_, A.real_ind_, real_nz_ * sizeof(int));
01608
01609 #ifdef SELDON_CHECK_MEMORY
01610 }
01611 catch (...)
01612 {
01613 this->m_ = 0;
01614 this->n_ = 0;
01615 real_nz_ = 0;
01616 imag_nz_ = 0;
01617 free(real_ptr_);
01618 free(imag_ptr_);
01619 real_ptr_ = NULL;
01620 imag_ptr_ = NULL;
01621 real_ind_ = NULL;
01622 imag_ind_ = NULL;
01623 this->real_data_ = NULL;
01624 this->imag_data_ = NULL;
01625 }
01626 if (real_ind_ == NULL)
01627 {
01628 this->m_ = 0;
01629 this->n_ = 0;
01630 real_nz_ = 0;
01631 imag_nz_ = 0;
01632 free(real_ptr_);
01633 free(imag_ptr_);
01634 real_ptr_ = NULL;
01635 imag_ptr_ = NULL;
01636 this->real_data_ = NULL;
01637 this->imag_data_ = NULL;
01638 }
01639 if (real_ind_ == NULL && i != 0 && j != 0)
01640 throw NoMemory(string("Matrix_ComplexSparse::")
01641 + "Matrix_ComplexSparse(int, int, int, int)",
01642 string("Unable to allocate ")
01643 + to_str(sizeof(int) * real_nz_)
01644 + " bytes to store " + to_str(real_nz_)
01645 + " row or column indices (real part), for a "
01646 + to_str(i) + " by " + to_str(j) + " matrix.");
01647 #endif
01648
01649 #ifdef SELDON_CHECK_MEMORY
01650 try
01651 {
01652 #endif
01653
01654 imag_ind_ = reinterpret_cast<int*>( calloc(imag_nz_, sizeof(int)) );
01655 memcpy(this->imag_ind_, A.imag_ind_, imag_nz_ * sizeof(int));
01656
01657 #ifdef SELDON_CHECK_MEMORY
01658 }
01659 catch (...)
01660 {
01661 this->m_ = 0;
01662 this->n_ = 0;
01663 real_nz_ = 0;
01664 imag_nz_ = 0;
01665 free(real_ptr_);
01666 free(imag_ptr_);
01667 real_ptr_ = NULL;
01668 imag_ptr_ = NULL;
01669 free(imag_ind_);
01670 real_ind_ = NULL;
01671 imag_ind_ = NULL;
01672 this->real_data_ = NULL;
01673 this->imag_data_ = NULL;
01674 }
01675 if (real_ind_ == NULL)
01676 {
01677 this->m_ = 0;
01678 this->n_ = 0;
01679 real_nz_ = 0;
01680 imag_nz_ = 0;
01681 free(real_ptr_);
01682 free(imag_ptr_);
01683 real_ptr_ = NULL;
01684 imag_ptr_ = NULL;
01685 free(imag_ind_);
01686 imag_ind_ = NULL;
01687 this->real_data_ = NULL;
01688 this->imag_data_ = NULL;
01689 }
01690 if (imag_ind_ == NULL && i != 0 && j != 0)
01691 throw NoMemory(string("Matrix_ComplexSparse::")
01692 + "Matrix_ComplexSparse(int, int, int, int)",
01693 string("Unable to allocate ")
01694 + to_str(sizeof(int) * imag_nz_)
01695 + " bytes to store " + to_str(imag_nz_)
01696 + " row or column indices (imaginary part), for a "
01697 + to_str(i) + " by " + to_str(j) + " matrix.");
01698 #endif
01699
01700 #ifdef SELDON_CHECK_MEMORY
01701 try
01702 {
01703 #endif
01704
01705 this->real_data_ = this->allocator_.allocate(real_nz_, this);
01706 this->allocator_.memorycpy(this->real_data_, A.real_data_, real_nz_);
01707
01708 #ifdef SELDON_CHECK_MEMORY
01709 }
01710 catch (...)
01711 {
01712 this->m_ = 0;
01713 this->n_ = 0;
01714 free(real_ptr_);
01715 free(imag_ptr_);
01716 real_ptr_ = NULL;
01717 imag_ptr_ = NULL;
01718 free(real_ind_);
01719 free(imag_ind_);
01720 real_ind_ = NULL;
01721 imag_ind_ = NULL;
01722 this->real_data_ = NULL;
01723 this->imag_data_ = NULL;
01724 }
01725 if (real_data_ == NULL)
01726 {
01727 this->m_ = 0;
01728 this->n_ = 0;
01729 free(real_ptr_);
01730 free(imag_ptr_);
01731 real_ptr_ = NULL;
01732 imag_ptr_ = NULL;
01733 free(real_ind_);
01734 free(imag_ind_);
01735 real_ind_ = NULL;
01736 imag_ind_ = NULL;
01737 imag_data_ = NULL;
01738 }
01739 if (real_data_ == NULL && i != 0 && j != 0)
01740 throw NoMemory(string("Matrix_ComplexSparse::")
01741 + "Matrix_ComplexSparse(int, int, int, int)",
01742 string("Unable to allocate ")
01743 + to_str(sizeof(int) * real_nz_)
01744 + " bytes to store " + to_str(real_nz_)
01745 + " values (real part), for a "
01746 + to_str(i) + " by " + to_str(j) + " matrix.");
01747 #endif
01748
01749 #ifdef SELDON_CHECK_MEMORY
01750 try
01751 {
01752 #endif
01753
01754 this->imag_data_ = this->allocator_.allocate(imag_nz_, this);
01755 this->allocator_.memorycpy(this->imag_data_, A.imag_data_, imag_nz_);
01756
01757 #ifdef SELDON_CHECK_MEMORY
01758 }
01759 catch (...)
01760 {
01761 this->m_ = 0;
01762 this->n_ = 0;
01763 free(real_ptr_);
01764 free(imag_ptr_);
01765 real_ptr_ = NULL;
01766 imag_ptr_ = NULL;
01767 free(real_ind_);
01768 free(imag_ind_);
01769 real_ind_ = NULL;
01770 imag_ind_ = NULL;
01771 this->allocator_.deallocate(this->real_data_, real_nz_);
01772 this->real_data_ = NULL;
01773 this->imag_data_ = NULL;
01774 }
01775 if (real_data_ == NULL)
01776 {
01777 this->m_ = 0;
01778 this->n_ = 0;
01779 free(real_ptr_);
01780 free(imag_ptr_);
01781 real_ptr_ = NULL;
01782 imag_ptr_ = NULL;
01783 free(real_ind_);
01784 free(imag_ind_);
01785 real_ind_ = NULL;
01786 imag_ind_ = NULL;
01787 this->allocator_.deallocate(this->real_data_, real_nz_);
01788 real_data_ = NULL;
01789 }
01790 if (imag_data_ == NULL && i != 0 && j != 0)
01791 throw NoMemory(string("Matrix_ComplexSparse::")
01792 + "Matrix_ComplexSparse(int, int, int, int)",
01793 string("Unable to allocate ")
01794 + to_str(sizeof(int) * imag_nz_)
01795 + " bytes to store " + to_str(imag_nz_)
01796 + " values (imaginary part), for a "
01797 + to_str(i) + " by " + to_str(j) + " matrix.");
01798 #endif
01799
01800 }
01801
01802
01803
01804
01805
01806
01807
01809
01815 template <class T, class Prop, class Storage, class Allocator>
01816 int Matrix_ComplexSparse<T, Prop, Storage, Allocator>::GetDataSize() const
01817 {
01818 return real_nz_ + imag_nz_;
01819 }
01820
01821
01823
01827 template <class T, class Prop, class Storage, class Allocator>
01828 int* Matrix_ComplexSparse<T, Prop, Storage, Allocator>::GetRealPtr() const
01829 {
01830 return real_ptr_;
01831 }
01832
01833
01835
01839 template <class T, class Prop, class Storage, class Allocator>
01840 int* Matrix_ComplexSparse<T, Prop, Storage, Allocator>::GetImagPtr() const
01841 {
01842 return imag_ptr_;
01843 }
01844
01845
01847
01854 template <class T, class Prop, class Storage, class Allocator>
01855 int* Matrix_ComplexSparse<T, Prop, Storage, Allocator>::GetRealInd() const
01856 {
01857 return real_ind_;
01858 }
01859
01860
01863
01870 template <class T, class Prop, class Storage, class Allocator>
01871 int* Matrix_ComplexSparse<T, Prop, Storage, Allocator>::GetImagInd() const
01872 {
01873 return imag_ind_;
01874 }
01875
01876
01878
01881 template <class T, class Prop, class Storage, class Allocator>
01882 int Matrix_ComplexSparse<T, Prop, Storage, Allocator>
01883 ::GetRealPtrSize() const
01884 {
01885 return (Storage::GetFirst(this->m_, this->n_) + 1);
01886 }
01887
01888
01890
01893 template <class T, class Prop, class Storage, class Allocator>
01894 int Matrix_ComplexSparse<T, Prop, Storage, Allocator>
01895 ::GetImagPtrSize() const
01896 {
01897 return (Storage::GetFirst(this->m_, this->n_) + 1);
01898 }
01899
01900
01903
01912 template <class T, class Prop, class Storage, class Allocator>
01913 int Matrix_ComplexSparse<T, Prop, Storage, Allocator>
01914 ::GetRealIndSize() const
01915 {
01916 return real_nz_;
01917 }
01918
01919
01922
01931 template <class T, class Prop, class Storage, class Allocator>
01932 int Matrix_ComplexSparse<T, Prop, Storage, Allocator>
01933 ::GetImagIndSize() const
01934 {
01935 return imag_nz_;
01936 }
01937
01938
01941
01950 template <class T, class Prop, class Storage, class Allocator>
01951 int Matrix_ComplexSparse<T, Prop, Storage, Allocator>
01952 ::GetRealDataSize() const
01953 {
01954 return real_nz_;
01955 }
01956
01957
01960
01969 template <class T, class Prop, class Storage, class Allocator>
01970 int Matrix_ComplexSparse<T, Prop, Storage, Allocator>
01971 ::GetImagDataSize() const
01972 {
01973 return imag_nz_;
01974 }
01975
01976
01978
01981 template <class T, class Prop, class Storage, class Allocator>
01982 T* Matrix_ComplexSparse<T, Prop, Storage, Allocator>::GetRealData() const
01983 {
01984 return real_data_;
01985 }
01986
01987
01989
01992 template <class T, class Prop, class Storage, class Allocator>
01993 T* Matrix_ComplexSparse<T, Prop, Storage, Allocator>::GetImagData() const
01994 {
01995 return imag_data_;
01996 }
01997
01998
01999
02000
02001
02002
02003
02005
02011 template <class T, class Prop, class Storage, class Allocator>
02012 inline complex<typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02013 ::value_type>
02014 Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02015 ::operator() (int i, int j) const
02016 {
02017
02018 #ifdef SELDON_CHECK_BOUNDS
02019 if (i < 0 || i >= this->m_)
02020 throw WrongRow("Matrix_ComplexSparse::operator()",
02021 string("Index should be in [0, ") + to_str(this->m_-1)
02022 + "], but is equal to " + to_str(i) + ".");
02023 if (j < 0 || j >= this->n_)
02024 throw WrongCol("Matrix_ComplexSparse::operator()",
02025 string("Index should be in [0, ") + to_str(this->n_-1)
02026 + "], but is equal to " + to_str(j) + ".");
02027 #endif
02028
02029 int real_k, imag_k, l;
02030 int real_a, real_b;
02031 int imag_a, imag_b;
02032
02033 real_a = real_ptr_[Storage::GetFirst(i, j)];
02034 real_b = real_ptr_[Storage::GetFirst(i, j) + 1];
02035
02036 imag_a = imag_ptr_[Storage::GetFirst(i, j)];
02037 imag_b = imag_ptr_[Storage::GetFirst(i, j) + 1];
02038
02039 if (real_a != real_b)
02040 {
02041 l = Storage::GetSecond(i, j);
02042 for (real_k = real_a;
02043 (real_k <real_b-1) && (real_ind_[real_k] < l);
02044 real_k++);
02045 if (imag_a != imag_b)
02046 {
02047 for (imag_k = imag_a;
02048 (imag_k < imag_b-1) && (imag_ind_[imag_k] < l);
02049 imag_k++);
02050 if (real_ind_[real_k] == l)
02051 {
02052 if (imag_ind_[imag_k] == l)
02053 return complex<T>(real_data_[real_k], imag_data_[imag_k]);
02054 else
02055 return complex<T>(real_data_[real_k], T(0));
02056 }
02057 else
02058 if (imag_ind_[imag_k] == l)
02059 return complex<T>(T(0), imag_data_[imag_k]);
02060 else
02061 return complex<T>(T(0), T(0));
02062 }
02063 else
02064 {
02065 if (real_ind_[real_k] == l)
02066 return complex<T>(real_data_[real_k], T(0));
02067 else
02068 return complex<T>(T(0), T(0));
02069 }
02070 }
02071 else
02072 {
02073 if (imag_a != imag_b)
02074 {
02075 l = Storage::GetSecond(i, j);
02076 for (imag_k = imag_a;
02077 (imag_k < imag_b-1) && (imag_ind_[imag_k] < l);
02078 imag_k++);
02079 if (imag_ind_[imag_k] == l)
02080 return complex<T>(T(0), imag_data_[imag_k]);
02081 else
02082 return complex<T>(T(0), T(0));
02083 }
02084 else
02085 return complex<T>(T(0), T(0));
02086 }
02087
02088 }
02089
02090
02092
02100 template <class T, class Prop, class Storage, class Allocator>
02101 inline typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02102 ::value_type&
02103 Matrix_ComplexSparse<T, Prop, Storage, Allocator>::ValReal(int i, int j)
02104 {
02105
02106 #ifdef SELDON_CHECK_BOUNDS
02107 if (i < 0 || i >= this->m_)
02108 throw WrongRow("Matrix_ComplexSparse::ValReal(int, int)",
02109 string("Index should be in [0, ") + to_str(this->m_-1)
02110 + "], but is equal to " + to_str(i) + ".");
02111 if (j < 0 || j >= this->n_)
02112 throw WrongCol("Matrix_ComplexSparse::ValReal(int, int)",
02113 string("Index should be in [0, ") + to_str(this->n_-1)
02114 + "], but is equal to " + to_str(j) + ".");
02115 #endif
02116
02117 int k, l;
02118 int a, b;
02119
02120 a = real_ptr_[Storage::GetFirst(i, j)];
02121 b = real_ptr_[Storage::GetFirst(i, j) + 1];
02122
02123 if (a == b)
02124 throw WrongArgument("Matrix_ComplexSparse::ValReal(int, int)",
02125 "No reference to element (" + to_str(i) + ", "
02126 + to_str(j)
02127 + ") can be returned: it is a zero entry.");
02128
02129 l = Storage::GetSecond(i, j);
02130
02131 for (k = a; (k < b-1) && (real_ind_[k] < l); k++);
02132
02133 if (real_ind_[k] == l)
02134 return this->real_data_[k];
02135 else
02136 throw WrongArgument("Matrix_ComplexSparse::ValReal(int, int)",
02137 "No reference to element (" + to_str(i) + ", "
02138 + to_str(j)
02139 + ") can be returned: it is a zero entry.");
02140 }
02141
02142
02144
02152 template <class T, class Prop, class Storage, class Allocator>
02153 inline const typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02154 ::value_type&
02155 Matrix_ComplexSparse<T, Prop, Storage, Allocator>::ValReal(int i, int j) const
02156 {
02157
02158 #ifdef SELDON_CHECK_BOUNDS
02159 if (i < 0 || i >= this->m_)
02160 throw WrongRow("Matrix_ComplexSparse::ValReal(int, int)",
02161 string("Index should be in [0, ") + to_str(this->m_-1)
02162 + "], but is equal to " + to_str(i) + ".");
02163 if (j < 0 || j >= this->n_)
02164 throw WrongCol("Matrix_ComplexSparse::ValReal(int, int)",
02165 string("Index should be in [0, ") + to_str(this->n_-1)
02166 + "], but is equal to " + to_str(j) + ".");
02167 #endif
02168
02169 int k, l;
02170 int a, b;
02171
02172 a = real_ptr_[Storage::GetFirst(i, j)];
02173 b = real_ptr_[Storage::GetFirst(i, j) + 1];
02174
02175 if (a == b)
02176 throw WrongArgument("Matrix_ComplexSparse::ValReal(int, int)",
02177 "No reference to element (" + to_str(i) + ", "
02178 + to_str(j)
02179 + ") can be returned: it is a zero entry.");
02180
02181 l = Storage::GetSecond(i, j);
02182
02183 for (k = a; (k < b-1) && (real_ind_[k] < l); k++);
02184
02185 if (real_ind_[k] == l)
02186 return this->real_data_[k];
02187 else
02188 throw WrongArgument("Matrix_ComplexSparse::ValReal(int, int)",
02189 "No reference to element (" + to_str(i) + ", "
02190 + to_str(j)
02191 + ") can be returned: it is a zero entry.");
02192 }
02193
02194
02196
02204 template <class T, class Prop, class Storage, class Allocator>
02205 inline typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02206 ::value_type&
02207 Matrix_ComplexSparse<T, Prop, Storage, Allocator>::ValImag(int i, int j)
02208 {
02209
02210 #ifdef SELDON_CHECK_BOUNDS
02211 if (i < 0 || i >= this->m_)
02212 throw WrongRow("Matrix_ComplexSparse::ValImag(int, int)",
02213 string("Index should be in [0, ") + to_str(this->m_-1)
02214 + "], but is equal to " + to_str(i) + ".");
02215 if (j < 0 || j >= this->n_)
02216 throw WrongCol("Matrix_ComplexSparse::ValImag(int, int)",
02217 string("Index should be in [0, ") + to_str(this->n_-1)
02218 + "], but is equal to " + to_str(j) + ".");
02219 #endif
02220
02221 int k, l;
02222 int a, b;
02223
02224 a = imag_ptr_[Storage::GetFirst(i, j)];
02225 b = imag_ptr_[Storage::GetFirst(i, j) + 1];
02226
02227 if (a == b)
02228 throw WrongArgument("Matrix_ComplexSparse::ValImag(int, int)",
02229 "No reference to element (" + to_str(i) + ", "
02230 + to_str(j)
02231 + ") can be returned: it is a zero entry.");
02232
02233 l = Storage::GetSecond(i, j);
02234
02235 for (k = a; (k < b-1) && (imag_ind_[k] < l); k++);
02236
02237 if (imag_ind_[k] == l)
02238 return this->imag_data_[k];
02239 else
02240 throw WrongArgument("Matrix_ComplexSparse::ValImag(int, int)",
02241 "No reference to element (" + to_str(i) + ", "
02242 + to_str(j)
02243 + ") can be returned: it is a zero entry.");
02244 }
02245
02246
02248
02256 template <class T, class Prop, class Storage, class Allocator>
02257 inline const typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02258 ::value_type&
02259 Matrix_ComplexSparse<T, Prop, Storage, Allocator>::ValImag(int i, int j) const
02260 {
02261
02262 #ifdef SELDON_CHECK_BOUNDS
02263 if (i < 0 || i >= this->m_)
02264 throw WrongRow("Matrix_ComplexSparse::ValImag(int, int)",
02265 string("Index should be in [0, ") + to_str(this->m_-1)
02266 + "], but is equal to " + to_str(i) + ".");
02267 if (j < 0 || j >= this->n_)
02268 throw WrongCol("Matrix_ComplexSparse::ValImag(int, int)",
02269 string("Index should be in [0, ") + to_str(this->n_-1)
02270 + "], but is equal to " + to_str(j) + ".");
02271 #endif
02272
02273 int k, l;
02274 int a, b;
02275
02276 a = imag_ptr_[Storage::GetFirst(i, j)];
02277 b = imag_ptr_[Storage::GetFirst(i, j) + 1];
02278
02279 if (a == b)
02280 throw WrongArgument("Matrix_ComplexSparse::ValImag(int, int)",
02281 "No reference to element (" + to_str(i) + ", "
02282 + to_str(j)
02283 + ") can be returned: it is a zero entry.");
02284
02285 l = Storage::GetSecond(i, j);
02286
02287 for (k = a; (k < b-1) && (imag_ind_[k] < l); k++);
02288
02289 if (imag_ind_[k] == l)
02290 return this->imag_data_[k];
02291 else
02292 throw WrongArgument("Matrix_ComplexSparse::ValImag(int, int)",
02293 "No reference to element (" + to_str(i) + ", "
02294 + to_str(j)
02295 + ") can be returned: it is a zero entry.");
02296 }
02297
02298
02300
02307 template <class T, class Prop, class Storage, class Allocator>
02308 inline typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02309 ::value_type&
02310 Matrix_ComplexSparse<T, Prop, Storage, Allocator>::GetReal(int i, int j)
02311 {
02312
02313 #ifdef SELDON_CHECK_BOUNDS
02314 if (i < 0 || i >= this->m_)
02315 throw WrongRow("Matrix_ComplexSparse::GetReal(int, int)",
02316 string("Index should be in [0, ") + to_str(this->m_-1)
02317 + "], but is equal to " + to_str(i) + ".");
02318 if (j < 0 || j >= this->n_)
02319 throw WrongCol("Matrix_ComplexSparse::GetReal(int, int)",
02320 string("Index should be in [0, ") + to_str(this->n_-1)
02321 + "], but is equal to " + to_str(j) + ".");
02322 #endif
02323
02324 int k, l;
02325 int a, b;
02326
02327 a = real_ptr_[Storage::GetFirst(i, j)];
02328 b = real_ptr_[Storage::GetFirst(i, j) + 1];
02329
02330 if (a < b)
02331 {
02332 l = Storage::GetSecond(i, j);
02333
02334 for (k = a; (k < b) && (real_ind_[k] < l); k++);
02335
02336 if ( (k < b) && (real_ind_[k] == l))
02337 return this->real_data_[k];
02338 }
02339 else
02340 k = a;
02341
02342
02343 Resize(this->m_, this->n_, real_nz_+1, imag_nz_);
02344
02345 for (int m = Storage::GetFirst(i, j)+1;
02346 m <= Storage::GetFirst(this->m_, this->n_); m++)
02347 real_ptr_[m]++;
02348
02349 for (int m = real_nz_-1; m >= k+1; m--)
02350 {
02351 real_ind_[m] = real_ind_[m-1];
02352 this->real_data_[m] = this->real_data_[m-1];
02353 }
02354
02355 real_ind_[k] = Storage::GetSecond(i, j);
02356
02357
02358 SetComplexZero(this->real_data_[k]);
02359
02360 return this->real_data_[k];
02361 }
02362
02363
02365
02372 template <class T, class Prop, class Storage, class Allocator>
02373 inline const typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02374 ::value_type&
02375 Matrix_ComplexSparse<T, Prop, Storage, Allocator>::GetReal(int i, int j) const
02376 {
02377 return ValReal(i, j);
02378 }
02379
02380
02382
02389 template <class T, class Prop, class Storage, class Allocator>
02390 inline typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02391 ::value_type&
02392 Matrix_ComplexSparse<T, Prop, Storage, Allocator>::GetImag(int i, int j)
02393 {
02394
02395 #ifdef SELDON_CHECK_BOUNDS
02396 if (i < 0 || i >= this->m_)
02397 throw WrongRow("Matrix_ComplexSparse::GetImag(int, int)",
02398 string("Index should be in [0, ") + to_str(this->m_-1)
02399 + "], but is equal to " + to_str(i) + ".");
02400 if (j < 0 || j >= this->n_)
02401 throw WrongCol("Matrix_ComplexSparse::GetImag(int, int)",
02402 string("Index should be in [0, ") + to_str(this->n_-1)
02403 + "], but is equal to " + to_str(j) + ".");
02404 #endif
02405
02406 int k, l;
02407 int a, b;
02408
02409 a = imag_ptr_[Storage::GetFirst(i, j)];
02410 b = imag_ptr_[Storage::GetFirst(i, j) + 1];
02411
02412 if (a < b)
02413 {
02414 l = Storage::GetSecond(i, j);
02415
02416 for (k = a; (k < b) && (imag_ind_[k] < l); k++);
02417
02418 if ( (k < b) && (imag_ind_[k] == l))
02419 return this->imag_data_[k];
02420 }
02421 else
02422 k = a;
02423
02424
02425 Resize(this->m_, this->n_, real_nz_, imag_nz_+1);
02426
02427 for (int m = Storage::GetFirst(i, j)+1;
02428 m <= Storage::GetFirst(this->m_, this->n_); m++)
02429 imag_ptr_[m]++;
02430
02431 for (int m = imag_nz_-1; m >= k+1; m--)
02432 {
02433 imag_ind_[m] = imag_ind_[m-1];
02434 this->imag_data_[m] = this->imag_data_[m-1];
02435 }
02436
02437 imag_ind_[k] = Storage::GetSecond(i, j);
02438
02439
02440 SetComplexZero(this->imag_data_[k]);
02441
02442 return this->imag_data_[k];
02443 }
02444
02445
02447
02454 template <class T, class Prop, class Storage, class Allocator>
02455 inline const typename Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02456 ::value_type&
02457 Matrix_ComplexSparse<T, Prop, Storage, Allocator>::GetImag(int i, int j) const
02458 {
02459 return ValImag(i, j);
02460 }
02461
02462
02464
02471 template <class T, class Prop, class Storage, class Allocator>
02472 inline void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02473 ::AddInteraction(int i, int j, const complex<T>& val)
02474 {
02475 if (real(val) != T(0))
02476 GetReal(i, j) += real(val);
02477
02478 if (imag(val) != T(0))
02479 GetImag(i, j) += imag(val);
02480 }
02481
02482
02484
02489 template <class T, class Prop, class Storage, class Allocator>
02490 inline void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02491 ::Set(int i, int j, const complex<T>& val)
02492 {
02493 GetReal(i, j) = real(val);
02494 GetImag(i, j) = imag(val);
02495 }
02496
02497
02499
02504 template <class T, class Prop, class Storage, class Allocator>
02505 inline Matrix_ComplexSparse<T, Prop, Storage, Allocator>&
02506 Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02507 ::operator= (const Matrix_ComplexSparse<T, Prop, Storage, Allocator>& A)
02508 {
02509 this->Copy(A);
02510
02511 return *this;
02512 }
02513
02514
02515
02516
02517
02518
02519
02521
02522 template <class T, class Prop, class Storage, class Allocator>
02523 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::Zero()
02524 {
02525 this->allocator_.memoryset(this->real_data_, char(0),
02526 this->real_nz_ * sizeof(value_type));
02527
02528 this->allocator_.memoryset(this->imag_data_, char(0),
02529 this->imag_nz_ * sizeof(value_type));
02530 }
02531
02532
02534
02537 template <class T, class Prop, class Storage, class Allocator>
02538 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::SetIdentity()
02539 {
02540 int m = this->m_;
02541 int n = this->n_;
02542 int nz = min(m, n);
02543
02544 if (nz == 0)
02545 return;
02546
02547 Clear();
02548
02549 Vector<T, VectFull, Allocator> real_values(nz), imag_values;
02550 Vector<int, VectFull, CallocAlloc<int> >
02551 real_ptr(Storage::GetFirst(m, n) + 1);
02552 Vector<int, VectFull, CallocAlloc<int> > real_ind(nz);
02553 Vector<int, VectFull, CallocAlloc<int> > imag_ptr(real_ptr);
02554 Vector<int, VectFull, CallocAlloc<int> > imag_ind;
02555
02556 real_values.Fill(T(1));
02557 real_ind.Fill();
02558 imag_ind.Zero();
02559 int i;
02560 for (i = 0; i < nz + 1; i++)
02561 real_ptr(i) = i;
02562 for (i = nz + 1; i < real_ptr.GetLength(); i++)
02563 real_ptr(i) = nz;
02564
02565 SetData(m, n, real_values, real_ptr, real_ind,
02566 imag_values, imag_ptr, imag_ind);
02567 }
02568
02569
02571
02574 template <class T, class Prop, class Storage, class Allocator>
02575 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::Fill()
02576 {
02577 for (int i = 0; i < this->real_nz_; i++)
02578 this->real_data_[i] = i;
02579
02580 for (int i = 0; i < this->imag_nz_; i++)
02581 this->imag_data_[i] = T(0);
02582 }
02583
02584
02586
02589 template <class T, class Prop, class Storage, class Allocator>
02590 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02591 ::Fill(const complex<T>& x)
02592 {
02593 for (int i = 0; i < this->real_nz_; i++)
02594 this->real_data_[i] = real(x);
02595
02596 for (int i = 0; i < this->imag_nz_; i++)
02597 this->imag_data_[i] = imag(x);
02598 }
02599
02600
02602
02605 template <class T, class Prop, class Storage, class Allocator>
02606 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::FillRand()
02607 {
02608 srand(time(NULL));
02609 for (int i = 0; i < this->real_nz_; i++)
02610 this->real_data_[i] = rand();
02611
02612 for (int i = 0; i < this->imag_nz_; i++)
02613 this->imag_data_[i] = rand();
02614 }
02615
02616
02618
02623 template <class T, class Prop, class Storage, class Allocator>
02624 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::Print() const
02625 {
02626 for (int i = 0; i < this->m_; i++)
02627 {
02628 for (int j = 0; j < this->n_; j++)
02629 cout << (*this)(i, j) << "\t";
02630 cout << endl;
02631 }
02632 }
02633
02634
02636
02640 template <class T, class Prop, class Storage, class Allocator>
02641 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02642 ::Write(string FileName) const
02643 {
02644 ofstream FileStream;
02645 FileStream.open(FileName.c_str());
02646
02647 #ifdef SELDON_CHECK_IO
02648
02649 if (!FileStream.is_open())
02650 throw IOError("Matrix_ComplexSparse::Write(string FileName)",
02651 string("Unable to open file \"") + FileName + "\".");
02652 #endif
02653
02654 this->Write(FileStream);
02655
02656 FileStream.close();
02657 }
02658
02659
02661
02665 template <class T, class Prop, class Storage, class Allocator>
02666 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02667 ::Write(ostream& FileStream) const
02668 {
02669 #ifdef SELDON_CHECK_IO
02670
02671 if (!FileStream.good())
02672 throw IOError("Matrix_ComplexSparse::Write(ofstream& FileStream)",
02673 "Stream is not ready.");
02674 #endif
02675
02676 FileStream.write(reinterpret_cast<char*>(const_cast<int*>(&this->m_)),
02677 sizeof(int));
02678 FileStream.write(reinterpret_cast<char*>(const_cast<int*>(&this->n_)),
02679 sizeof(int));
02680 FileStream.write(reinterpret_cast<char*>(const_cast<int*>(&this->real_nz_)),
02681 sizeof(int));
02682 FileStream.write(reinterpret_cast<char*>(const_cast<int*>(&this->imag_nz_)),
02683 sizeof(int));
02684
02685 FileStream.write(reinterpret_cast<char*>(this->real_ptr_),
02686 sizeof(int)*(Storage::GetFirst(this->m_, this->n_)+1));
02687 FileStream.write(reinterpret_cast<char*>(this->real_ind_),
02688 sizeof(int)*this->real_nz_);
02689 FileStream.write(reinterpret_cast<char*>(this->real_data_),
02690 sizeof(T)*this->real_nz_);
02691
02692 FileStream.write(reinterpret_cast<char*>(this->imag_ptr_),
02693 sizeof(int)*(Storage::GetFirst(this->m_, this->n_)+1));
02694 FileStream.write(reinterpret_cast<char*>(this->imag_ind_),
02695 sizeof(int)*this->imag_nz_);
02696 FileStream.write(reinterpret_cast<char*>(this->imag_data_),
02697 sizeof(T)*this->imag_nz_);
02698 }
02699
02700
02702
02708 template <class T, class Prop, class Storage, class Allocator>
02709 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::
02710 WriteText(string FileName) const
02711 {
02712 ofstream FileStream; FileStream.precision(14);
02713 FileStream.open(FileName.c_str());
02714
02715 #ifdef SELDON_CHECK_IO
02716
02717 if (!FileStream.is_open())
02718 throw IOError("Matrix_ComplexSparse::Write(string FileName)",
02719 string("Unable to open file \"") + FileName + "\".");
02720 #endif
02721
02722 this->WriteText(FileStream);
02723
02724 FileStream.close();
02725 }
02726
02727
02729
02735 template <class T, class Prop, class Storage, class Allocator>
02736 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::
02737 WriteText(ostream& FileStream) const
02738 {
02739
02740 #ifdef SELDON_CHECK_IO
02741
02742 if (!FileStream.good())
02743 throw IOError("Matrix_ComplexSparse::Write(ofstream& FileStream)",
02744 "Stream is not ready.");
02745 #endif
02746
02747
02748 IVect IndRow, IndCol; Vector<complex<T> > Value;
02749 const Matrix<T, Prop, Storage, Allocator>& leaf_class =
02750 static_cast<const Matrix<T, Prop, Storage, Allocator>& >(*this);
02751
02752 ConvertMatrix_to_Coordinates(leaf_class, IndRow, IndCol,
02753 Value, 1);
02754
02755 for (int i = 0; i < IndRow.GetM(); i++)
02756 FileStream << IndRow(i) << " " << IndCol(i) << " " << Value(i) << '\n';
02757
02758 }
02759
02760
02762
02766 template <class T, class Prop, class Storage, class Allocator>
02767 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>::
02768 Read(string FileName)
02769 {
02770 ifstream FileStream;
02771 FileStream.open(FileName.c_str());
02772
02773 #ifdef SELDON_CHECK_IO
02774
02775 if (!FileStream.is_open())
02776 throw IOError("Matrix_ComplexSparse::Read(string FileName)",
02777 string("Unable to open file \"") + FileName + "\".");
02778 #endif
02779
02780 this->Read(FileStream);
02781
02782 FileStream.close();
02783 }
02784
02785
02787
02791 template <class T, class Prop, class Storage, class Allocator>
02792 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02793 ::Read(istream& FileStream)
02794 {
02795
02796 #ifdef SELDON_CHECK_IO
02797
02798 if (!FileStream.good())
02799 throw IOError("Matrix_ComplexSparse::Read(istream& FileStream)",
02800 "Stream is not ready.");
02801 #endif
02802
02803 int m, n, real_nz, imag_nz;
02804 FileStream.read(reinterpret_cast<char*>(&m), sizeof(int));
02805 FileStream.read(reinterpret_cast<char*>(&n), sizeof(int));
02806 FileStream.read(reinterpret_cast<char*>(&real_nz), sizeof(int));
02807 FileStream.read(reinterpret_cast<char*>(&imag_nz), sizeof(int));
02808
02809 Reallocate(m, n, real_nz, imag_nz);
02810
02811 FileStream.read(reinterpret_cast<char*>(real_ptr_),
02812 sizeof(int)*(Storage::GetFirst(m, n)+1));
02813 FileStream.read(reinterpret_cast<char*>(real_ind_), sizeof(int)*real_nz);
02814 FileStream.read(reinterpret_cast<char*>(this->real_data_), sizeof(T)*real_nz);
02815
02816 FileStream.read(reinterpret_cast<char*>(imag_ptr_),
02817 sizeof(int)*(Storage::GetFirst(m, n)+1));
02818 FileStream.read(reinterpret_cast<char*>(imag_ind_), sizeof(int)*imag_nz);
02819 FileStream.read(reinterpret_cast<char*>(this->imag_data_), sizeof(T)*imag_nz);
02820
02821 #ifdef SELDON_CHECK_IO
02822
02823 if (!FileStream.good())
02824 throw IOError("Matrix_ComplexSparse::Read(istream& FileStream)",
02825 string("Input operation failed.")
02826 + string(" The input file may have been removed")
02827 + " or may not contain enough data.");
02828 #endif
02829
02830 }
02831
02832
02834
02838 template <class T, class Prop, class Storage, class Allocator>
02839 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02840 ::ReadText(string FileName)
02841 {
02842 ifstream FileStream;
02843 FileStream.open(FileName.c_str());
02844
02845 #ifdef SELDON_CHECK_IO
02846
02847 if (!FileStream.is_open())
02848 throw IOError("Matrix_ComplexSparse::ReadText(string FileName)",
02849 string("Unable to open file \"") + FileName + "\".");
02850 #endif
02851
02852 this->ReadText(FileStream);
02853
02854 FileStream.close();
02855 }
02856
02857
02859
02863 template <class T, class Prop, class Storage, class Allocator>
02864 void Matrix_ComplexSparse<T, Prop, Storage, Allocator>
02865 ::ReadText(istream& FileStream)
02866 {
02867 Matrix<T, Prop, Storage, Allocator>& leaf_class =
02868 static_cast<Matrix<T, Prop, Storage, Allocator>& >(*this);
02869
02870 complex<T> zero; int index = 1;
02871 ReadCoordinateMatrix(leaf_class, FileStream, zero, index);
02872 }
02873
02874
02876
02878
02879
02880
02881
02882
02883
02885
02888 template <class T, class Prop, class Allocator>
02889 Matrix<T, Prop, ColComplexSparse, Allocator>::Matrix():
02890 Matrix_ComplexSparse<T, Prop, ColComplexSparse, Allocator>()
02891 {
02892 }
02893
02894
02896
02900 template <class T, class Prop, class Allocator>
02901 Matrix<T, Prop, ColComplexSparse, Allocator>::Matrix(int i, int j):
02902 Matrix_ComplexSparse<T, Prop, ColComplexSparse, Allocator>(i, j, 0, 0)
02903 {
02904 }
02905
02906
02908
02916 template <class T, class Prop, class Allocator>
02917 Matrix<T, Prop, ColComplexSparse, Allocator>::Matrix(int i, int j,
02918 int real_nz,
02919 int imag_nz):
02920 Matrix_ComplexSparse<T, Prop, ColComplexSparse, Allocator>(i, j,
02921 real_nz,
02922 imag_nz)
02923 {
02924 }
02925
02926
02928
02946 template <class T, class Prop, class Allocator>
02947 template <class Storage0, class Allocator0,
02948 class Storage1, class Allocator1,
02949 class Storage2, class Allocator2>
02950 Matrix<T, Prop, ColComplexSparse, Allocator>::
02951 Matrix(int i, int j,
02952 Vector<T, Storage0, Allocator0>& real_values,
02953 Vector<int, Storage1, Allocator1>& real_ptr,
02954 Vector<int, Storage2, Allocator2>& real_ind,
02955 Vector<T, Storage0, Allocator0>& imag_values,
02956 Vector<int, Storage1, Allocator1>& imag_ptr,
02957 Vector<int, Storage2, Allocator2>& imag_ind):
02958 Matrix_ComplexSparse<T, Prop, ColComplexSparse, Allocator>(i, j,
02959 real_values,
02960 real_ptr,
02961 real_ind,
02962 imag_values,
02963 imag_ptr,
02964 imag_ind)
02965 {
02966 }
02967
02968
02969
02971
02973
02974
02975
02976
02977
02978
02980
02983 template <class T, class Prop, class Allocator>
02984 Matrix<T, Prop, RowComplexSparse, Allocator>::Matrix():
02985 Matrix_ComplexSparse<T, Prop, RowComplexSparse, Allocator>()
02986 {
02987 }
02988
02989
02991
02995 template <class T, class Prop, class Allocator>
02996 Matrix<T, Prop, RowComplexSparse, Allocator>::Matrix(int i, int j):
02997 Matrix_ComplexSparse<T, Prop, RowComplexSparse, Allocator>(i, j, 0, 0)
02998 {
02999 }
03000
03001
03010 template <class T, class Prop, class Allocator>
03011 Matrix<T, Prop, RowComplexSparse, Allocator>::Matrix(int i, int j,
03012 int real_nz,
03013 int imag_nz):
03014 Matrix_ComplexSparse<T, Prop, RowComplexSparse, Allocator>(i, j,
03015 real_nz,
03016 imag_nz)
03017 {
03018 }
03019
03020
03022
03040 template <class T, class Prop, class Allocator>
03041 template <class Storage0, class Allocator0,
03042 class Storage1, class Allocator1,
03043 class Storage2, class Allocator2>
03044 Matrix<T, Prop, RowComplexSparse, Allocator>::
03045 Matrix(int i, int j,
03046 Vector<T, Storage0, Allocator0>& real_values,
03047 Vector<int, Storage1, Allocator1>& real_ptr,
03048 Vector<int, Storage2, Allocator2>& real_ind,
03049 Vector<T, Storage0, Allocator0>& imag_values,
03050 Vector<int, Storage1, Allocator1>& imag_ptr,
03051 Vector<int, Storage2, Allocator2>& imag_ind):
03052 Matrix_ComplexSparse<T, Prop, RowComplexSparse, Allocator>(i, j,
03053 real_values,
03054 real_ptr,
03055 real_ind,
03056 imag_values,
03057 imag_ptr,
03058 imag_ind)
03059 {
03060 }
03061
03062
03063 }
03064
03065 #define SELDON_FILE_MATRIX_COMPLEXSPARSE_CXX
03066 #endif
