00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 #ifndef __SPMATRIX_H
00035 #define __SPMATRIX_H
00036
00037
00038 #include "Matrix.h"
00039 #include "CG.h"
00040 #include "BiCG.h"
00041 #include "CGS.h"
00042 #include "BiCG.h"
00043 #include "BiCGStab.h"
00044 #include "GMRes.h"
00045 #include "QMR.h"
00046
00047 namespace OFELI {
00048
00053 template<class T_> class Prec;
00054
00055 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00056 template<class T_>
00057 struct ErrorInSpMatrix {
00058 void Message(const char *file, int line, int code, int p1=0, int p2=0) {
00059 cerr << "\n*** Fatal Error in OFELI ***\nIn SpMatrix<>::";
00060 cerr << "File : " << file << ", line : " << line << "\nIn SkSMatrix<>::";
00061 switch (code) {
00062 case 1:
00063 cerr << "Set(i,j,x) : Index pair : (" << p1 << "," << p2 << ") is not compatible "
00064 << "with sparse storage." << endl;
00065 break;
00066 }
00067 cerr << "Program stops" << endl;
00068 exit(1);
00069 }
00070 };
00071 #endif
00072
00100 template<class T_> class SpMatrix : public Matrix<T_>
00101 {
00102 using Matrix<T_>::_rmin;
00103 using Matrix<T_>::_cmin;
00104 using Matrix<T_>::_rmax;
00105 using Matrix<T_>::_cmax;
00106 using Matrix<T_>::_nb_rows;
00107 using Matrix<T_>::_nb_cols;
00108 using Matrix<T_>::_size;
00109 using Matrix<T_>::_length;
00110 using Matrix<T_>::_zero;
00111 using Matrix<T_>::_temp;
00112 using Matrix<T_>::_a;
00113 using Matrix<T_>::_diag;
00114 using Matrix<T_>::_row_ptr;
00115 using Matrix<T_>::_col_ind;
00116 using Matrix<T_>::_ch;
00117 using Matrix<T_>::_dof_type;
00118
00119 public:
00120
00123 SpMatrix() { }
00124
00131 SpMatrix(size_t nb_rows, size_t nb_cols);
00132
00138 SpMatrix(size_t size);
00139
00152 SpMatrix(const class Mesh &mesh, size_t dof=0, int type=0);
00153
00154 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00155 SpMatrix(size_t dof, const class Mesh &mesh, int code=0);
00156 #endif
00157
00158 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00159 SpMatrix(size_t dof, size_t nb_eq, const class Mesh &mesh);
00160 #endif
00161
00173 SpMatrix(const Vect<size_t> &I, const Vect<size_t> &J, size_t nr=0, size_t nc=0, int opt=1);
00174
00188 SpMatrix(const Vect<size_t> &I, const Vect<size_t> &J, const Vect<T_> &a, size_t nr=0, size_t nc=0, int opt=1);
00189
00197 SpMatrix(size_t nr, size_t nc, const Vect<unsigned long> &row_ptr, const Vect<size_t> &col_ind);
00198
00206 SpMatrix(size_t nr, size_t nc, const Vect<unsigned long> &row_ptr, const Vect<size_t> &col_ind, const Vect<T_> &a);
00207
00213 SpMatrix(const Vect<unsigned long> &row_ptr, const Vect<size_t> &col_ind);
00214
00221 SpMatrix(const Vect<unsigned long> &row_ptr, const Vect<size_t> &col_ind, const Vect<T_> &a);
00222
00224 SpMatrix(const SpMatrix &m);
00225
00227 ~SpMatrix(void);
00228
00242 virtual void setMesh(const class Mesh &mesh, size_t dof=0, size_t type=0);
00243
00244 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00245 void setMesh(size_t dof, const class Mesh &mesh, int code=0);
00246 #endif
00247
00248 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00249 void setMesh(size_t dof, size_t nb_eq, const class Mesh &mesh);
00250 #endif
00251
00252 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00253
00258 void ExtendedGraph();
00259 #endif
00260
00262 void setOneDOF();
00263
00265 void setSides();
00266
00268 void setDiag();
00269
00271 Vect<T_> getColumn(size_t j) const;
00272
00283 void DiagPrescribe(const class Mesh &mesh, Vect<T_> &b, const Vect<T_> &u);
00284
00287 void Resize(size_t size);
00288
00292 void Resize(size_t nr, size_t nc);
00293
00295 Vect<T_> getRow(size_t i) const;
00296
00301 T_ & operator()(size_t i, size_t j);
00302
00307 T_ operator()(size_t i, size_t j) const;
00308
00314 const T_ operator()(size_t i) const;
00315
00321 const T_ operator[](size_t i) const;
00322
00325 void getMesh(const class Mesh &mesh);
00326
00331 void Mult(const Vect<T_> &v, Vect<T_> &w) const;
00332
00337 void MultAdd(const Vect<T_> &v, Vect<T_> &w) const;
00338
00344 void MultAdd(T_ a, const Vect<T_> &v, Vect<T_> &w) const;
00345
00350 void TMult(const Vect<T_> &v, Vect<T_> &w) const;
00351
00357 void Set(size_t i, size_t j, const T_ &x);
00358
00364 void Add(size_t i, size_t j, const T_ &x);
00365
00368 void operator=(const T_ &x);
00369
00372 size_t getColInd(size_t i) const;
00373
00375 size_t getRowPtr(size_t i) const;
00376
00377 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00379 int Factor();
00380 #endif
00381
00395 int Solve(Vect<T_> &b);
00396
00411 int Solve(const Vect<T_> &b, Vect<T_> &x);
00412
00435 void setSolver(Iteration solver=CG_SOLVER, Preconditioner prec=DIAG_PREC, int max_it=1000, double toler=1.e-8);
00436
00439 T_ *getArray() const;
00440
00442 T_ getEntry(size_t i, size_t j) const;
00443
00444 private:
00445
00446 #ifdef USE_MTL
00447 T_ *_A;
00448 #endif
00449 mutable ErrorInSpMatrix<T_> _e;
00450 int insert_(size_t i, size_t j);
00451 int _col_index(size_t i, size_t j) const
00452 {
00453 int k = 0;
00454 while (k<int(_row_ptr[i]-_row_ptr[i-1]))
00455 if (j==_col_ind[_row_ptr[i-1]+k-1])
00456 return k;
00457 else
00458 k++;
00459 return -1;
00460 }
00461
00462 void _set(size_t i, size_t j, const T_ &x) { _a[_row_ptr[i-1]+_col_index(i,j)-1] = x; }
00463 size_t _dof;
00464 int _type, _solver, _prec, _max_it;
00465 double _toler;
00466 bool _extended, _one_dof, _sides;
00467 };
00468
00469
00471
00473
00474
00475 template<class T_>
00476 SpMatrix<T_>::SpMatrix(size_t nb_rows, size_t nb_cols)
00477 {
00478 register size_t i, j;
00479 _nb_rows = nb_rows;
00480 _nb_cols = nb_cols;
00481 _size = 0;
00482 if (_nb_rows == _nb_cols)
00483 _size = _nb_rows;
00484 _length = nb_rows * nb_cols;
00485 _row_ptr.resize(_nb_rows+1);
00486 _col_ind.resize(_length);
00487 _row_ptr[0] = 1;
00488 for (i=1; i<=_nb_rows; i++)
00489 _row_ptr[i] = _row_ptr[i-1] + _nb_cols;
00490 size_t l = 0;
00491 for (i=0; i<_nb_rows; i++)
00492 for (j=0; j<_nb_cols; j++)
00493 _col_ind[l++] = j+1;
00494 _a.resize(_length,T_(0.));
00495 _rmin = 1;
00496 _rmax = _nb_rows;
00497 _cmin = 1;
00498 _cmax = _nb_cols;
00499 _solver = CG_SOLVER;
00500 _prec = DIAG_PREC;
00501 _max_it = 1000;
00502 _toler = 1.e-8;
00503 #ifdef USE_MTL
00504 _A = new MTL_Mat(_nb_rows, _nb_cols, _length, _a, _row_ptr, _col_ind);
00505 #endif
00506 #ifdef _OFELI_DEBUG
00507 std::clog << "An instance of class SpMatrix is constructed in\n";
00508 std::clog << "File :" << __FILE__ << ", Line: " << __LINE__ << endl;
00509 #endif
00510 }
00511
00512
00513 template<class T_>
00514 SpMatrix<T_>::SpMatrix(size_t size)
00515 {
00516 register size_t i;
00517 _size = _nb_rows = _nb_cols = size;
00518 _length = _size * _size;
00519 _row_ptr.resize(_nb_rows+1);
00520 _col_ind.resize(_length);
00521 _row_ptr[0] = 1;
00522 for (i=1; i<=_nb_rows; i++)
00523 _row_ptr[i] = _row_ptr[i-1] + _nb_cols;
00524 size_t l = 0;
00525 for (i=0; i<_nb_rows; i++)
00526 for (size_t j=0; j<_nb_cols; j++)
00527 _col_ind[l++] = j+1;
00528 _a.resize(_length,T_(0.));
00529 _rmin = 1;
00530 _rmax = _nb_rows;
00531 _cmin = 1;
00532 _cmax = _nb_cols;
00533 _solver = CG_SOLVER;
00534 _prec = DIAG_PREC;
00535 _max_it = 1000;
00536 _toler = 1.e-8;
00537 #ifdef USE_MTL
00538 _A = new MTL_Mat(_size, _size, _length, _a, _row_ptr, _col_ind);
00539 #endif
00540 #ifdef _OFELI_DEBUG
00541 std::clog << "An instance of class SpMatrix is constructed in\n";
00542 std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
00543 #endif
00544 }
00545
00546
00547 template<class T_>
00548 SpMatrix<T_>::SpMatrix(const class Mesh &mesh, size_t dof, int type)
00549 {
00550 setMesh(mesh,dof,type);
00551 #ifdef _OFELI_DEBUG
00552 std::clog << "An instance of class SpMatrix is constructed in\n";
00553 std::clog << "File : " << __FILE__<< ", Line : " << __LINE__ << endl;
00554 #endif
00555 #ifdef USE_MTL
00556 _A = new MTL_Mat(_size, _size, _length, _a, _row_ptr, _col_ind);
00557 #endif
00558 _rmin = 1;
00559 _rmax = _nb_rows;
00560 _cmin = 1;
00561 _cmax = _nb_cols;
00562 _solver = CG_SOLVER;
00563 _prec = DIAG_PREC;
00564 _max_it = 1000;
00565 _toler = 1.e-8;
00566 }
00567
00568
00569 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00570 template<class T_>
00571 SpMatrix<T_>::SpMatrix(size_t dof, const class Mesh &mesh, int code)
00572 {
00573 setMesh(dof,mesh,code);
00574 #ifdef USE_MTL
00575 _A = new MTL_Mat(_size, _size, _length, _a, _row_ptr, _col_ind);
00576 #endif
00577 #ifdef _OFELI_DEBUG
00578 std::clog << "An instance of class SpMatrix is constructed in\n";
00579 std::clog << "File : " << __FILE__<< ", Line : " << __LINE__ << endl;
00580 #endif
00581 _rmin = 1;
00582 _rmax = _nb_rows;
00583 _cmin = 1;
00584 _cmax = _nb_cols;
00585 _solver = CG_SOLVER;
00586 _prec = DIAG_PREC;
00587 _max_it = 1000;
00588 _toler = 1.e-8;
00589 }
00590 #endif
00591
00592
00593 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00594 template<class T_>
00595 SpMatrix<T_>::SpMatrix(size_t dof, size_t nb_eq, const class Mesh &mesh)
00596 {
00597 setMesh(dof,nb_eq,mesh);
00598 #ifdef USE_MTL
00599 _A = new MTL_Mat(_size, _size, _length, _a, _row_ptr, _col_ind);
00600 #endif
00601 #ifdef _OFELI_DEBUG
00602 std::clog << "An instance of class SpMatrix is constructed in\n";
00603 std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
00604 #endif
00605 _rmin = 1;
00606 _rmax = _nb_rows;
00607 _cmin = 1;
00608 _cmax = _nb_cols;
00609 _solver = CG_SOLVER;
00610 _prec = DIAG_PREC;
00611 _max_it = 1000;
00612 _toler = 1.e-8;
00613 }
00614 #endif
00615
00616
00617 template<class T_>
00618 SpMatrix<T_>::SpMatrix(const Vect<size_t> &I, const Vect<size_t> &J, size_t nr, size_t nc, int opt)
00619 {
00620 size_t i, n=I.size();
00621 std::vector<RC> pp(n);
00622 _size = 0;
00623 for (i=0; i<n; i++) {
00624 _size = max(_size,I[i]);
00625 pp[i] = RC(I[i]-1,J[i]-1);
00626 }
00627 _nb_rows = nr;
00628 _nb_cols = nc;
00629 if (nr==0)
00630 _nb_rows = _size;
00631 if (nc==0)
00632 _nb_rows = _size;
00633 _size = _nb_rows;
00634 if (opt==0) {
00635 std::vector<RC>::iterator it = pp.begin();
00636 sort(pp.begin(),pp.end());
00637 std::vector<RC>::iterator new_end = std::unique(pp.begin(),pp.end());
00638 pp.erase(new_end,pp.end());
00639 }
00640 _row_ptr.resize(_size+1);
00641 _col_ind.resize(pp.size());
00642 StoreGraph(_size,pp,_row_ptr,_col_ind);
00643 _length = pp.size();
00644 _a.resize(_length,0);
00645 for (i=0; i<_length; i++)
00646 _a[i] = 0;
00647 _rmin = 1;
00648 _rmax = _nb_rows;
00649 _cmin = 1;
00650 _cmax = _nb_cols;
00651 _solver = CG_SOLVER;
00652 _prec = DIAG_PREC;
00653 _max_it = 1000;
00654 _toler = 1.e-8;
00655 }
00656
00657
00658 template<class T_>
00659 SpMatrix<T_>::SpMatrix(const Vect<size_t> &I, const Vect<size_t> &J, const Vect<T_> &a, size_t nr, size_t nc, int opt)
00660 {
00661 size_t n = I.size();
00662 std::vector<RC> pp(n);
00663 _size = 0;
00664 for (size_t i=0; i<n; i++) {
00665 _size = max(_size,I[i]);
00666 pp[i] = RC(I[i]-1,J[i]-1);
00667 }
00668 _nb_rows = nr;
00669 _nb_cols = nc;
00670 if (nr==0)
00671 _nb_rows = _size;
00672 if (nc==0)
00673 _nb_rows = _size;
00674 _size = _nb_rows;
00675 if (opt==0) {
00676 sort(pp.begin(),pp.end());
00677 std::vector<RC>::iterator new_end = std::unique(pp.begin(),pp.end());
00678 pp.erase(new_end,pp.end());
00679 }
00680 _row_ptr.resize(_size+1);
00681 _col_ind.resize(n);
00682 StoreGraph(_size,pp,_row_ptr,_col_ind);
00683 _length = pp.size();
00684 _a.resize(_length);
00685 for (size_t j=0; j<n; j++)
00686 _a[_row_ptr[I[j]-1]+_col_index(I[j],J[j])-1] = a[j];
00687 _rmin = 1;
00688 _rmax = _nb_rows;
00689 _cmin = 1;
00690 _cmax = _nb_cols;
00691 _solver = CG_SOLVER;
00692 _prec = DIAG_PREC;
00693 _max_it = 1000;
00694 _toler = 1.e-8;
00695 }
00696
00697
00698 template<class T_>
00699 SpMatrix<T_>::SpMatrix(size_t nr, size_t nc, const Vect<unsigned long> &row_ptr, const Vect<size_t> &col_ind)
00700 : _type(0), _dof(0)
00701 {
00702 _nb_rows = nr;
00703 _nb_cols = nc;
00704 _size = 0;
00705 _length = col_ind.size();
00706 _row_ptr.resize(_nb_rows+1);
00707 _row_ptr = row_ptr;
00708 _col_ind.resize(_length);
00709 _col_ind = col_ind;
00710 _a.resize(_length,T_(0.));
00711 _diag.resize(_size);
00712 _rmin = 1;
00713 _rmax = _nb_rows;
00714 _cmin = 1;
00715 _cmax = _nb_cols;
00716 _solver = CG_SOLVER;
00717 _prec = DIAG_PREC;
00718 _max_it = 1000;
00719 _toler = 1.e-8;
00720 #ifdef USE_MTL
00721 _A = new MTL_Mat(_size, _size, _length, _a, _row_ptr, _col_ind);
00722 #endif
00723 #ifdef _OFELI_DEBUG
00724 std::clog << "An instance of class SpMatrix is constructed in\n";
00725 std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
00726 #endif
00727 }
00728
00729
00730 template<class T_>
00731 SpMatrix<T_>::SpMatrix(size_t nr, size_t nc, const Vect<unsigned long> &row_ptr, const Vect<size_t> &col_ind, const Vect<T_> &a)
00732 : _type(0), _dof(0)
00733 {
00734 _size = 0;
00735 _nb_rows = nr;
00736 _nb_cols = nc;
00737 _length = col_ind.size();
00738 _row_ptr.resize(_nb_rows+1);
00739 _row_ptr = row_ptr;
00740 _col_ind.resize(_length);
00741 _col_ind = col_ind;
00742 _a.resize(_length);
00743 _a = a;
00744 _diag.resize(_size);
00745 _rmin = 1;
00746 _rmax = _nb_rows;
00747 _cmin = 1;
00748 _cmax = _nb_cols;
00749 _solver = CG_SOLVER;
00750 _prec = DIAG_PREC;
00751 _max_it = 1000;
00752 _toler = 1.e-8;
00753 #ifdef USE_MTL
00754 _A = new MTL_Mat(_size, _size, _length, _a, _row_ptr, _col_ind);
00755 #endif
00756 #ifdef _OFELI_DEBUG
00757 std::clog << "An instance of class SpMatrix is constructed in\n";
00758 std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
00759 #endif
00760 }
00761
00762
00763 template<class T_>
00764 SpMatrix<T_>::SpMatrix(const Vect<unsigned long> &row_ptr, const Vect<size_t> &col_ind)
00765 : _type(0), _dof(0)
00766 {
00767 size_t i;
00768 _size = _nb_rows = _nb_cols = row_ptr.size()-1;
00769 _length = col_ind.size();
00770 _row_ptr.resize(_size+1);
00771 _row_ptr = row_ptr;
00772 _col_ind.resize(_length);
00773 _col_ind = col_ind;
00774 _a.resize(_length,T_(0.));
00775 _diag.resize(_size);
00776 _rmin = 1;
00777 _rmax = _nb_rows;
00778 _cmin = 1;
00779 _cmax = _nb_cols;
00780 _solver = CG_SOLVER;
00781 _prec = DIAG_PREC;
00782 _max_it = 1000;
00783 _toler = 1.e-8;
00784 #ifdef USE_MTL
00785 _A = new MTL_Mat(_size, _size, _length, _a, _row_ptr, _col_ind);
00786 #endif
00787 #ifdef _OFELI_DEBUG
00788 std::clog << "An instance of class SpMatrix is constructed in\n";
00789 std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00790 #endif
00791 }
00792
00793
00794 template<class T_>
00795 SpMatrix<T_>::SpMatrix(const Vect<unsigned long> &row_ptr, const Vect<size_t> &col_ind, const Vect<T_> &a)
00796 : _type(0), _dof(0)
00797 {
00798 size_t i;
00799 _size = _nb_rows = _nb_cols = row_ptr.size()-1;
00800 _length = col_ind.size();
00801 _row_ptr.resize(_size+1);
00802 _row_ptr = row_ptr;
00803 _col_ind.resize(_length);
00804 _col_ind = col_ind;
00805 _a.resize(_length);
00806 _a = a;
00807 _diag.resize(_size);
00808 _rmin = 1;
00809 _rmax = _nb_rows;
00810 _cmin = 1;
00811 _cmax = _nb_cols;
00812 _solver = CG_SOLVER;
00813 _prec = DIAG_PREC;
00814 _max_it = 1000;
00815 _toler = 1.e-8;
00816 #ifdef USE_MTL
00817 _A = new MTL_Mat(_size, _size, _length, _a, _row_ptr, _col_ind);
00818 #endif
00819 #ifdef _OFELI_DEBUG
00820 std::clog << "An instance of class SpMatrix is constructed in\n";
00821 std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
00822 #endif
00823 }
00824
00825
00826 template<class T_>
00827 SpMatrix<T_>::SpMatrix(const SpMatrix &m)
00828 {
00829 _dof = m._dof;
00830 _size = m._size;
00831 _length = m._length;
00832 _row_ptr.resize(_size+1);
00833 _col_ind.resize(_length);
00834 _col_ind = m._col_ind;
00835 _row_ptr = m._row_ptr;
00836 _a.resize(_length);
00837 _a = m._a;
00838 _diag.resize(_size);
00839 _diag = m._diag;
00840 _rmin = m._rmin;
00841 _rmax = m._rmax;
00842 _cmin = m._cmin;
00843 _cmax = m._cmax;
00844 _solver = CG_SOLVER;
00845 _prec = DIAG_PREC;
00846 _max_it = 1000;
00847 _toler = 1.e-8;
00848 _solver = m._solver;
00849 _prec = m._prec;
00850 _max_it = 1000;
00851 _toler = 1.e-8;
00852 #ifdef USE_MTL
00853 _A = new MTL_Mat(_size, _size, _length, _a, _row_ptr, _col_ind);
00854 #endif
00855 #ifdef _OFELI_DEBUG
00856 std::clog << "An instance of class SpMatrix is constructed in\n";
00857 std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
00858 #endif
00859 }
00860
00861
00862 template<class T_>
00863 SpMatrix<T_>::~SpMatrix(void)
00864 {
00865 #ifdef USE_MTL
00866 delete _A;
00867 _A = NULL;
00868 #endif
00869 }
00870
00871
00872
00873 template<class T_>
00874 void SpMatrix<T_>::Resize(size_t size)
00875 {
00876 _nb_rows = _nb_cols = _size = size;
00877 _length = lsize_t(_nb_rows*_nb_cols);
00878 _rmin = 1;
00879 _rmax = _nb_rows;
00880 _cmin = 1;
00881 _cmax = _nb_cols;
00882 _row_ptr.resize(_nb_rows+1);
00883 _col_ind.resize(_length);
00884 _row_ptr[0] = 1;
00885 size_t i;
00886 for (i=1; i<=_nb_rows; i++)
00887 _row_ptr[i] = _row_ptr[i-1] + _nb_cols;
00888 size_t l = 0;
00889 for (i=0; i<_nb_rows; i++)
00890 for (size_t j=0; j<_nb_cols; j++)
00891 _col_ind[l++] = j+1;
00892 _a.resize(_length,T_(0.));
00893 }
00894
00895
00896 template<class T_>
00897 void SpMatrix<T_>::Resize(size_t nr, size_t nc)
00898 {
00899 _nb_rows = nr;
00900 _nb_cols = nc;
00901 _size = 0;
00902 if (_nb_rows==_nb_cols)
00903 _size = _nb_rows;
00904 _length = lsize_t(_nb_rows*_nb_cols);
00905 _rmin = 1;
00906 _rmax = _nb_rows;
00907 _cmin = 1;
00908 _cmax = _nb_cols;
00909 _row_ptr.resize(_nb_rows+1);
00910 _col_ind.resize(_length);
00911 _row_ptr[0] = 1;
00912 size_t i;
00913 for (i=1; i<=_nb_rows; i++)
00914 _row_ptr[i] = _row_ptr[i-1] + _nb_cols;
00915 size_t l = 0;
00916 for (i=0; i<_nb_rows; i++)
00917 for (size_t j=0; j<_nb_cols; j++)
00918 _col_ind[l++] = j+1;
00919 _a.resize(_length,T_(0.));
00920 }
00921
00922
00923 template<class T_>
00924 void SpMatrix<T_>::setMesh(const class Mesh &mesh, size_t dof, size_t type)
00925 {
00926 Matrix<T_>::_init_set_mesh(mesh,dof,type);
00927 if (_dof_type==NODE_DOF) {
00928 if (type && dof)
00929 _length = XGraph(mesh, _row_ptr, _col_ind);
00930 else if (dof)
00931 _length = NodeGraphScal(mesh, _row_ptr, _col_ind);
00932 else if (type==0 && dof==0)
00933 _length = NodeGraph(mesh, _row_ptr, _col_ind);
00934 }
00935 else if (_dof_type==SIDE_DOF)
00936 _length = SideGraph(mesh, _row_ptr, _col_ind);
00937
00938 else if (_dof_type==ELEMENT_DOF)
00939 _length = ElementGraph(mesh, _row_ptr, _col_ind);
00940
00941 else;
00942
00943 _a.resize(_length,T_(0.));
00944 _diag.resize(_size);
00945 }
00946
00947
00948 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00949 template<class T_>
00950 void SpMatrix<T_>::setMesh(size_t dof, const class Mesh &mesh, int code)
00951 {
00952 _size = _nb_rows = _nb_cols = mesh.getNbEq();
00953 if (dof)
00954 _size = _nb_rows = _nb_cols = mesh.getNbNodes();
00955 if (code!=0)
00956 _length = XGraph(mesh, _row_ptr, _col_ind);
00957 else
00958 _length = NodeGraphScal(mesh, _row_ptr, _col_ind);
00959 _a.resize(_length,T_(0.));
00960 _diag.resize(_size);
00961 }
00962 #endif
00963
00964
00965 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00966 template<class T_>
00967 void SpMatrix<T_>::setMesh(size_t dof, size_t nb_eq, const class Mesh &mesh)
00968 {
00969 _type = 0;
00970 _dof = 0;
00971 _size = _nb_rows = _nb_cols = nb_eq;
00972 _length = NodeGraphScal(mesh, dof, nb_eq, _row_ptr, _col_ind);
00973 _a.resize(_length,T_(0.));
00974 _diag.resize(_size);
00975 }
00976 #endif
00977
00978
00979 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00980 template<class T_>
00981 void SpMatrix<T_>::ExtendedGraph()
00982 {
00983 _extended = true;
00984 }
00985 #endif
00986
00987
00988 template<class T_>
00989 void SpMatrix<T_>::setOneDOF()
00990 {
00991 _one_dof = true;
00992 }
00993
00994
00995 template<class T_>
00996 void SpMatrix<T_>::setSides()
00997 {
00998 _sides = true;
00999 }
01000
01001
01002 template<class T_>
01003 void SpMatrix<T_>::setDiag()
01004 {
01005 for (size_t i=0; i<_size; i++) {
01006 for (size_t j=0; j<_row_ptr[i+1]-_row_ptr[i]; j++) {
01007 if (i==j) {
01008 _diag[i] = _a[_row_ptr[i]+i-1];
01009 break;
01010 }
01011 }
01012 }
01013 }
01014
01015
01016 template<class T_>
01017 Vect<T_> SpMatrix<T_>::getColumn(size_t j) const
01018 {
01019 Vect<T_> v(_nb_rows);
01020 for (size_t i=1; i<=_nb_rows; i++)
01021 v(i) = getEntry(i,j);
01022 return v;
01023 }
01024
01025
01026 template<class T_>
01027 void SpMatrix<T_>::DiagPrescribe(const class Mesh &mesh, Vect<T_> &b, const Vect<T_> &u)
01028 {
01029 double p = 0;
01030 for (size_t j=1; j<=_nb_rows; j++)
01031 p = max(p,Abs(getEntry(j,j)));
01032
01033 const Node *nd;
01034 size_t k=0;
01035 for (mesh.topNode(); (nd=mesh.getNode());) {
01036 for (size_t i=1; i<=nd->getNbDOF(); ++i) {
01037 size_t ii = nd->getDOF(i)-1;
01038 for (size_t j=0; j<_row_ptr[ii+1]-_row_ptr[ii]; j++,k++)
01039 if (nd->getCode(i)>0) {
01040 b[ii] = p*u[ii];
01041 _a[k] = 0;
01042 if (ii+1==_col_ind[k])
01043 _a[k] = p;
01044 }
01045 }
01046 }
01047 }
01048
01049
01050 template<class T_>
01051 Vect<T_> SpMatrix<T_>::getRow(size_t i) const
01052 {
01053 return getColumn(i);
01054 }
01055
01056
01057 template<class T_>
01058 T_ & SpMatrix<T_>::operator()(size_t i, size_t j)
01059 {
01060 int k = _col_index(i,j);
01061 try {
01062 if (k<0)
01063 throw(1);
01064 else
01065 return _a[_row_ptr[i-1]+k-1];
01066 }
01067 catch(int n) { _e.Message(__FILE__,__LINE__,n,int(i),int(j)); }
01068 return _temp;
01069 }
01070
01071
01072 template<class T_>
01073 T_ SpMatrix<T_>::operator()(size_t i, size_t j) const
01074 {
01075 int k = _col_index(i,j);
01076 if (k<0)
01077 return _zero;
01078 else
01079 return _a[_row_ptr[i-1]+k-1];
01080 }
01081
01082
01083 template<class T_>
01084 const T_ SpMatrix<T_>::operator()(size_t i) const
01085 {
01086 return _a[i-1];
01087 }
01088
01089 template<class T_>
01090 const T_ SpMatrix<T_>::operator[](size_t i) const
01091 {
01092 return _a[i];
01093 }
01094
01095
01096 template<class T_>
01097 void SpMatrix<T_>::getMesh(const class Mesh &mesh)
01098 {
01099 if (_sides)
01100 SideGraph(mesh,_row_ptr,_col_ind);
01101 else {
01102 if (_dof)
01103 _size = _nb_rows = _nb_cols = mesh.getNbNodes();
01104 else
01105 _size = _nb_rows = _nb_cols = mesh.getNbEq();
01106 if (_type)
01107 XGraph(mesh,_row_ptr,_col_ind);
01108
01109 else
01110 NodeGraph(mesh,_row_ptr,_col_ind);
01111 }
01112 _a.resize(_length,T_(0.));
01113 _diag.resize(_size);
01114 }
01115
01116
01117 template<class T_>
01118 void SpMatrix<T_>::Mult(const Vect<T_> &v, Vect<T_> &w) const
01119 {
01120 #ifdef USE_MTL
01121 mtl::external_vec<double> dv(v.getArray(),v.getSize());
01122 mtl::external_vec<double> dw(w.getArray(),w.getSize());
01123 mtl::mult(*_A, dv, dw);
01124 #else
01125 Clear(w);
01126 MultAdd(v,w);
01127 #endif
01128 }
01129
01130
01131 template<class T_>
01132 void SpMatrix<T_>::MultAdd(const Vect<T_> &v, Vect<T_> &w) const
01133 {
01134 #ifdef USE_MTL
01135 mtl::external_vec<double> dv(v.Array(),v.Size());
01136 mtl::external_vec<double> dw(w.Array(),w.Size());
01137 mtl::mult_add(*_A, dv, dw) ;
01138 #else
01139 size_t l = 0;
01140 for (size_t i=0; i<_nb_rows; ++i)
01141 for (size_t j=0; j<_row_ptr[i+1]-_row_ptr[i]; ++j)
01142 w[i] += _a[_row_ptr[i]+j-1] * v(_col_ind[l++]);
01143 #endif
01144 }
01145
01146
01147 template<class T_>
01148 void SpMatrix<T_>::MultAdd(T_ a, const Vect<T_> &v, Vect<T_> &w) const
01149 {
01150 #ifdef USE_MTL
01151 mtl::external_vec<double> dv(v.Array(),v.Size());
01152 mtl::external_vec<double> dw(w.Array(),w.Size());
01153 cout << "Not implemented." << endl;
01154 mtl::mult_add(*_A, dv, dw) ;
01155 else
01156 size_t l = 0;
01157 for (size_t i=0; i<_nb_rows; i++)
01158 for (size_t j=0; j<_row_ptr[i+1]-_row_ptr[i]; j++)
01159 w[i] += a * _a[_row_ptr[i]+j-1] * v(_col_ind[l++]);
01160 #endif
01161 }
01162
01163
01164 template<class T_>
01165 void SpMatrix<T_>::TMult(const Vect<T_> &v, Vect<T_> &w) const
01166 {
01167 for (size_t i=0; i<_nb_rows; i++)
01168 for (size_t j=0; j<_row_ptr[i+1]-_row_ptr[i]; j++)
01169 w[_col_ind[j+_row_ptr[i]-1]-1] += _a[_row_ptr[i]+j-1] * v[i];
01170 }
01171
01172
01173 template<class T_>
01174 void SpMatrix<T_>::Set(size_t i, size_t j, const T_ &x)
01175 {
01176 int k = _col_index(i,j);
01177 if (k<0)
01178 _e.Message(__FILE__,__LINE__,1,int(i),int(j));
01179 else
01180 _a[_row_ptr[i-1]+_col_index(i,j)-1] = x;
01181 }
01182
01183
01184 template<class T_>
01185 void SpMatrix<T_>::Add(size_t i, size_t j, const T_ &x)
01186 {
01187 _a[_row_ptr[i-1]+_col_index(i,j)-1] += x;
01188 }
01189
01190
01191 template<class T_>
01192 void SpMatrix<T_>::operator=(const T_ &x)
01193 {
01194 for (size_t i=0; i<_length; i++)
01195 _a[i] = x;
01196 }
01197
01198
01199 template<class T_>
01200 size_t SpMatrix<T_>::getColInd(size_t i) const
01201 {
01202 #ifdef _BOUNDS_DEBUG
01203 assert(i<=_length && i>0);
01204 #endif
01205 return _col_ind[i-1];
01206 }
01207
01208
01209 template<class T_>
01210 size_t SpMatrix<T_>::getRowPtr(size_t i) const
01211 {
01212 #ifdef _BOUNDS_DEBUG
01213 assert(i<=_nb_rows+1 && i>0);
01214 #endif
01215 return _row_ptr[i-1];
01216 }
01217
01218
01219 #ifndef DOXYGEN_SHOULD_SKIP_THIS
01220 template<class T_>
01221 int SpMatrix<T_>::Factor()
01222 {
01223 return -1;
01224 }
01225 #endif
01226
01227
01228 template<class T_>
01229 int SpMatrix<T_>::Solve(Vect<T_> &b)
01230 {
01231 Vect<T_> x(b.getSize());
01232 int nb_it = Solve(b,x);
01233 b = x;
01234 return nb_it;
01235 }
01236
01237
01238 template<class T_>
01239 int SpMatrix<T_>::Solve(const Vect<T_> &b, Vect<T_> &x)
01240 {
01241 Prec<T_> *p = new Prec<T_>(*this,_prec);
01242 int nb_it = 0;
01243
01244 switch (_solver) {
01245
01246 case CG_SOLVER:
01247 nb_it = CG(*this,*p,b,x,_max_it,_toler,1);
01248 break;
01249
01250 case BICG_SOLVER:
01251 nb_it = BiCG(*this,*p,b,x,_max_it,_toler,1);
01252 break;
01253
01254 case CGS_SOLVER:
01255 nb_it = CGS(*this,*p,b,x,_max_it,_toler,1);
01256 break;
01257
01258 case BICG_STAB_SOLVER:
01259 nb_it = BiCGStab(*this,*p,b,x,_max_it,_toler,1);
01260 break;
01261
01262 case GMRES_SOLVER:
01263
01264 nb_it = GMRes(*this,*p,b,x,5,_max_it,_toler,1);
01265 break;
01266
01267 case QMR_SOLVER:
01268 nb_it = QMR(*this,*p,*p,b,x,_max_it,_toler,1);
01269 break;
01270 }
01271 delete p;
01272 return nb_it;
01273 }
01274
01275
01276 template<class T_>
01277 void SpMatrix<T_>::setSolver(Iteration solver, Preconditioner prec, int max_it, double toler)
01278 {
01279 _solver = solver;
01280 _prec = prec;
01281 _max_it = max_it;
01282 _toler = toler;
01283 }
01284
01285
01286 template<class T_>
01287 T_ * SpMatrix<T_>::getArray() const
01288 {
01289 return _a;
01290 }
01291
01292
01293 template<class T_>
01294 T_ SpMatrix<T_>::getEntry(size_t i, size_t j) const
01295 {
01296 int k = _col_index(i,j);
01297 if (k<0)
01298 return _zero;
01299 else
01300 return _a[_row_ptr[i-1]+k-1];
01301 }
01302
01303
01305
01307
01308
01313 template<class T_>
01314 ostream& operator<<(ostream& s, const SpMatrix<T_> &a)
01315 {
01316 #ifdef USE_MTL
01317 mtl::print_all_matrix(*(a._A));
01318 #else
01319 size_t rmin, rmax, cmin, cmax;
01320 a.getPrintView(rmin,rmax,cmin,cmax);
01321 s.setf(ios::right|ios::scientific);
01322 s << endl;
01323 for (size_t i=rmin; i<=rmax; i++) {
01324 s << "\nRow: " << setw(6) << i << endl;
01325 for (size_t j=cmin; j<=cmax; j++)
01326 s << " " << setprecision(8) << std::setfill(' ') << setw(18) << a(i,j);
01327 s << endl;
01328 }
01329 #endif
01330 return s;
01331 }
01332
01333 }
01334
01335 #endif
