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 __SKSMATRIX_H
00035 #define __SKSMATRIX_H
00036
00037
00038 #include "Matrix.h"
00039 #include "Vect.h"
00040
00041 namespace OFELI {
00042
00071 template<class T_> class SkSMatrix;
00072
00073 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00074 template<class T_>
00075 struct ErrorInSkSMatrix {
00076 void Message(const char *file, size_t line, int code, int p1=0, int p2=0, T_ p3=0) {
00077 cerr << "\n*** Fatal Error in OFELI ***\n";
00078 cerr << "File : " << file << ", line : " << line << "\nIn SkSMatrix<>::";
00079 switch (code) {
00080
00081 case 51:
00082 cerr << "Set(i,j,x) : Index pair : (" << p1 << "," << p2 << ") is not compatible "
00083 << "with skyline symmeric storage." << endl;
00084 break;
00085
00086 case 61:
00087 cerr << "Factor() : The " << p1 << "th pivot = " << p3 << " is too small\n";
00088 break;
00089
00090 case 71:
00091 cerr << "Solve(b) : Warning : Matrix has not been factorized\n";
00092 return;
00093 }
00094 cerr << "Program stops" << endl;
00095 exit(1);
00096 }
00097 };
00098 #endif
00099
00100 template<class T_>
00101 class SkSMatrix : public Matrix<T_>
00102 {
00103
00104 using Matrix<T_>::_rmin;
00105 using Matrix<T_>::_cmin;
00106 using Matrix<T_>::_rmax;
00107 using Matrix<T_>::_cmax;
00108 using Matrix<T_>::_nb_rows;
00109 using Matrix<T_>::_nb_cols;
00110 using Matrix<T_>::_size;
00111 using Matrix<T_>::_length;
00112 using Matrix<T_>::_zero;
00113 using Matrix<T_>::_dof_type;
00114 using Matrix<T_>::_temp;
00115 using Matrix<T_>::_fact;
00116 using Matrix<T_>::_a;
00117 using Matrix<T_>::_diag;
00118 using Matrix<T_>::_ch;
00119 using Matrix<T_>::_is_diagonal;
00120
00121 public:
00122
00124 SkSMatrix();
00125
00132 SkSMatrix(size_t size, bool is_diagonal=false);
00133
00143 SkSMatrix(const class Mesh &mesh, size_t dof=0, bool is_diagonal=false);
00144
00149 SkSMatrix(const Vect<size_t> &ColHt);
00150
00158 SkSMatrix(const Vect<size_t> &I, const Vect<size_t> &J, int opt=1);
00159
00169 SkSMatrix(const Vect<size_t> &I, const Vect<size_t> &J, const Vect<T_> &a, int opt=1);
00170
00172 SkSMatrix(const SkSMatrix<T_> &m);
00173
00175 ~SkSMatrix();
00176
00187 virtual void setMesh(const class Mesh &mesh, size_t dof=0, size_t type=0);
00188
00189 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00190 void setMesh(size_t dof, const class Mesh &mesh, int code=0);
00191 #endif
00192
00193 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00194 void setMesh(size_t dof, size_t nb_eq, const class Mesh &mesh);
00195 #endif
00196
00199 void setSkyline(const class Mesh &mesh);
00200
00202 void setDiag();
00203
00209 void Set(size_t i, size_t j, const T_ &x);
00210
00211 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00212 void SSet(size_t i, size_t j, const T_ &x);
00213 #endif
00214
00219 void MultAdd(const Vect<T_> &x, Vect<T_> &y) const;
00220
00226 void MultAdd(T_ a, const Vect<T_> &x, Vect<T_> &y) const;
00227
00232 void Mult(const Vect<T_> &x, Vect<T_> &y) const;
00233
00238 void TMult(const Vect<T_> &x, Vect<T_> &y) const;
00239
00245 void Add(size_t i, size_t j, const T_ &x);
00246
00249 size_t getColHeight(size_t i) const;
00250
00252 Vect<T_> getColumn(size_t j) const;
00253
00255 Vect<T_> getRow(size_t i) const;
00256
00262 const T_ operator()(size_t i) const { return _a[i-1]; }
00263
00268 T_ & operator()(size_t i, size_t j);
00269
00274 T_ operator()(size_t i, size_t j) const;
00275
00278 SkSMatrix<T_> & operator=(const SkSMatrix<T_> &m);
00279
00282 SkSMatrix<T_> & operator=(const T_ &x);
00283
00286 SkSMatrix<T_> & operator+=(const SkSMatrix<T_> &m);
00287
00290 SkSMatrix<T_> & operator*=(const T_ &x);
00291
00296 int Factor();
00297
00303 int Solve(Vect<T_> &b);
00304
00311 int Solve(const Vect<T_> &b, Vect<T_> &x);
00312
00315 T_ *getArray() const { return _a; }
00316
00318 void setArray(size_t i, T_ x) { _a[i] = x; }
00319
00321 T_ getEntry(size_t i, size_t j) const;
00322
00323 private:
00324
00325 mutable ErrorInSkSMatrix<T_> _e;
00326 int _dof;
00327 void _set(size_t i, size_t j, const T_ &x) { _a[_ch[i-1]+j-i] = x; }
00328 };
00329
00331
00333
00334
00335 template<class T_>
00336 SkSMatrix<T_>::SkSMatrix() : _dof(0)
00337 {
00338 _fact = false;
00339 _is_diagonal = false;
00340 _dof_type = NODE_DOF;
00341 #ifdef _OFELI_DEBUG
00342 std::clog << "An instance of class SkSMatrix is constructed.\n";
00343 std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00344 #endif
00345 }
00346
00347
00348 template<class T_>
00349 SkSMatrix<T_>::SkSMatrix(size_t size, bool is_diagonal) : _dof(0)
00350 {
00351 _zero = 0;
00352 _fact = false;
00353 _is_diagonal = false;
00354 _dof_type = NODE_DOF;
00355 _size = size;
00356 _ch.resize(size);
00357 _diag.resize(_size);
00358 _ch[0] = 0;
00359 for (size_t i=1; i<_size; i++)
00360 _ch[i] = _ch[i-1] + i + 1;
00361 _length = _ch[_size-1] + 1;
00362 _a.resize(_length);
00363 Clear(_a);
00364 #ifdef _OFELI_DEBUG
00365 std::clog << "An instance of class SkSMatrix is constructed.\n";
00366 std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00367 #endif
00368 }
00369
00370
00371 template<class T_>
00372 SkSMatrix<T_>::SkSMatrix(const class Mesh &mesh, size_t dof, bool is_diagonal)
00373 {
00374 _is_diagonal = is_diagonal;
00375 setMesh(mesh,dof);
00376 }
00377
00378
00379 template<class T_>
00380 SkSMatrix<T_>::SkSMatrix(const Vect<size_t> &ColHt) : _dof(0)
00381 {
00382 _fact = false;
00383 _is_diagonal = false;
00384 _dof_type = NODE_DOF;
00385 _size = ColHt.getSize();
00386 _zero = 0;
00387 _ch.resize(_size);
00388 _ch[0] = 0;
00389 for (size_t i=1; i<_size; ++i)
00390 _ch[i] = _ch[i-1] + ColHt[i];
00391 _length = _ch[_size-1]+1;
00392 _a.resize(_length,0);
00393 _diag.resize(_size);
00394 #ifdef _OFELI_DEBUG
00395 std::clog << "An instance of class SkSMatrix is constructed.\n";
00396 std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00397 #endif
00398 }
00399
00400
00401 template<class T_>
00402 SkSMatrix<T_>::SkSMatrix(const Vect<size_t> &I, const Vect<size_t> &J, int opt)
00403 : _dof(0)
00404 {
00405 _size = 0;
00406 _fact = false;
00407 _is_diagonal = false;
00408 _dof_type = NODE_DOF;
00409 size_t i;
00410 size_t n = I.size();
00411 std::vector<RC> pp(n);
00412 for (i=0; i<n; i++) {
00413 _size = max(_size,I[i]);
00414 pp[i] = RC(I[i]-1,J[i]-1);
00415 }
00416 _ch.resize(_size,0);
00417 _nb_rows = _nb_cols = _size;
00418 if (opt==0) {
00419 std::vector<RC>::iterator it = pp.begin();
00420 sort(pp.begin(),pp.end());
00421 vector<RC>::iterator new_end = std::unique(pp.begin(),pp.end());
00422 pp.erase(new_end,pp.end());
00423 }
00424 for (i=0; i<n; i++) {
00425 if (I[i]>J[i])
00426 _ch[I[i]-1] = max(static_cast<unsigned>(abs(int(I[i])-int(J[i]))),_ch[I[i]-1]);
00427 }
00428 _ch[0] = 0;
00429 for (i=1; i<_size; ++i)
00430 _ch[i] += _ch[i-1] + 1;
00431 _length = _ch[_size-1]+1;
00432 _a.resize(_length,0);
00433 _diag.resize(_size);
00434 }
00435
00436
00437 template<class T_>
00438 SkSMatrix<T_>::SkSMatrix(const Vect<size_t> &I, const Vect<size_t> &J, const Vect<T_> &a, int opt)
00439 : _dof(0)
00440 {
00441 _fact = false;
00442 _dof_type = NODE_DOF;
00443 size_t i;
00444 size_t n = I.size();
00445 std::vector<RC> pp(n);
00446 _size = 0;
00447 for (i=0; i<n; i++) {
00448 _size = max(_size,I[i]);
00449 pp[i] = RC(I[i]-1,J[i]-1);
00450 }
00451 _ch.resize(_size,0);
00452 _nb_rows = _nb_cols = _size;
00453 if (opt==0) {
00454 std::vector<RC>::iterator it = pp.begin();
00455 sort(pp.begin(),pp.end());
00456 vector<RC>::iterator new_end = std::unique(pp.begin(),pp.end());
00457 pp.erase(new_end,pp.end());
00458 }
00459 for (i=0; i<n; i++)
00460 if (I[i]>J[i])
00461 _ch[I[i]-1] = max(static_cast<unsigned>(abs(int(I[i])-int(J[i]))),_ch[I[i]-1]);
00462 _ch[0] = 0;
00463 for (i=1; i<_size; ++i)
00464 _ch[i] += _ch[i-1] + 1;
00465 _length = _ch[_size-1]+1;
00466 _a.resize(_length,0);
00467 size_t k=0;
00468 for (i=0; i<n; i++)
00469 Set(I[i],J[i],a[k++]);
00470 _diag.resize(_size);
00471 }
00472
00473
00474 template<class T_>
00475 SkSMatrix<T_>::SkSMatrix(const SkSMatrix<T_> &m)
00476 {
00477 _length = m._length;
00478 _size = m._size;
00479 _ch.resize(_size);
00480 _ch = m._ch;
00481 _a.resize(_length);
00482 _a = m._a;
00483 _diag.resize(_size);
00484 _diag = m._diag;
00485 _fact = m._fact;
00486 _dof = m._dof;
00487 _zero = T_(0);
00488 _is_diagonal = m._is_diagonal;
00489 #ifdef _OFELI_DEBUG
00490 std::clog << "An instance of class SkSMatrix is constructed.\n";
00491 std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00492 #endif
00493 }
00494
00495
00496 template<class T_>
00497 SkSMatrix<T_>::~SkSMatrix()
00498 {
00499 #ifdef _OFELI_DEBUG
00500 std::clog << "An instance of class SkSMatrix is constructed.\n";
00501 std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00502 #endif
00503 }
00504
00505
00506 template<class T_>
00507 void SkSMatrix<T_>::setMesh(const class Mesh &mesh, size_t dof, size_t type)
00508 {
00509 Matrix<T_>::_init_set_mesh(mesh,dof,type);
00510 if (mesh.NodesAreDOF())
00511 if (dof)
00512 _length = NodeSkyline(mesh,_ch,dof);
00513 else
00514 _length = NodeSkyline(mesh,_ch);
00515 else if (mesh.SidesAreDOF())
00516 if (dof)
00517 _length = SideSkyline(mesh,_ch,dof);
00518 else
00519 _length = SideSkyline(mesh,_ch);
00520 else if (mesh.ElementsAreDOF())
00521 if (dof)
00522 _length = ElementSkyline(mesh,_ch,dof);
00523 else
00524 _length = ElementSkyline(mesh,_ch);
00525 else;
00526 _diag.resize(_size);
00527 _a.resize(_length);
00528 _ch[0] = 0;
00529 for (size_t i=1; i<_size; i++)
00530 _ch[i] += _ch[i-1];
00531 Clear(_a);
00532 _fact = false;
00533 _dof = 0;
00534 _nb_rows = _nb_cols = _size;
00535 }
00536
00537
00538 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00539 template<class T_>
00540 void SkSMatrix<T_>::setMesh(size_t dof, const class Mesh &mesh, int code)
00541 {
00542
00543 dof = 0;
00544 if (mesh.getDim()==0) { }
00545 code = 0;
00546 }
00547 #endif
00548
00549
00550 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00551 template<class T_>
00552 void SkSMatrix<T_>::setMesh(size_t dof, size_t nb_eq, const class Mesh &mesh)
00553 {
00554
00555 dof = 0;
00556 nb_eq = 0;
00557 if (mesh.getDim()==0) { }
00558 }
00559 #endif
00560
00561
00562 template<class T_>
00563 void SkSMatrix<T_>::setSkyline(const class Mesh &mesh)
00564 {
00565 _zero = 0;
00566 bool set_nodes = mesh.NodesAreDOF();
00567 bool set_sides = mesh.SidesAreDOF();
00568 bool set_elements = mesh.ElementsAreDOF();
00569 _size = mesh.getNbEq();
00570 if (_dof)
00571 _size = mesh.getNbNodes();
00572 if (set_sides) {
00573 if (_dof) {
00574 _size = mesh.getNbSides();
00575 _length = SideSkyline(mesh,_ch,_dof);
00576 }
00577 _length = SideSkyline(mesh,_ch);
00578 }
00579 else {
00580 if (_dof) {
00581 _size = mesh.getNbNodes();
00582 _length = NodeSkyline(mesh,_ch,_dof);
00583 }
00584 _length = NodeSkyline(mesh,_ch);
00585 }
00586 _diag.resize(_size);
00587 _ch[0] = 0;
00588 for (size_t i=1; i<_size; i++)
00589 _ch[i] += _ch[i-1];
00590 _a.resize(_length);
00591 Clear(_a);
00592 _fact = false;
00593 }
00594
00595
00596 template<class T_>
00597 void SkSMatrix<T_>::setDiag()
00598 {
00599 for (size_t i=0; i<_size; i++)
00600 _diag[i] = _a[_ch[i]];
00601 }
00602
00603
00604 template<class T_>
00605 void SkSMatrix<T_>::Set(size_t i, size_t j, const T_ &x)
00606 {
00607 try {
00608 if (i>=j)
00609 _a[_ch[i-1]+j-i] = x;
00610 else
00611 throw(51);
00612 }
00613 catch(int n) { _e.Message(__FILE__,__LINE__,n,int(i),int(j)); }
00614 }
00615
00616
00617 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00618 template<class T_>
00619 void SkSMatrix<T_>::SSet(size_t i, size_t j, const T_ &x)
00620 {
00621 _a[_ch[i-1]+j-i] = x;
00622 }
00623 #endif
00624
00625
00626 template<class T_>
00627 void SkSMatrix<T_>::MultAdd(const Vect<T_> &x, Vect<T_> &y) const
00628 {
00629 for (size_t i=0; i<_size; i++) {
00630 for (size_t j=0; j<_size; j++)
00631 y[i] += operator()(i+1,j+1)*x[j];
00632 }
00633 }
00634
00635
00636 template<class T_>
00637 void SkSMatrix<T_>::MultAdd(T_ a, const Vect<T_> &x, Vect<T_> &y) const
00638 {
00639 for (size_t i=0; i<_size; i++)
00640 for (size_t j=0; j<_size; j++)
00641 y[i] += a * operator()(i+1,j+1)*x[j];
00642 }
00643
00644
00645 template<class T_>
00646 void SkSMatrix<T_>::Mult(const Vect<T_> &x, Vect<T_> &y) const
00647 {
00648 Clear(y);
00649 MultAdd(x,y);
00650 }
00651
00652
00653 template<class T_>
00654 void SkSMatrix<T_>::TMult(const Vect<T_> &x, Vect<T_> &y) const
00655 {
00656 Mult(x,y);
00657 }
00658
00659
00660 template<class T_>
00661 void SkSMatrix<T_>::Add(size_t i, size_t j, const T_ &x)
00662 {
00663 if (i>=j)
00664 _a[_ch[i-1]+j-i] += x;
00665 }
00666
00667
00668 template<class T_>
00669 size_t SkSMatrix<T_>::getColHeight(size_t i) const
00670 {
00671 #ifdef _BOUNDS_DEBUG
00672 assert(i>0);
00673 assert(i<=_size);
00674 #endif
00675 if (i==1)
00676 return 1;
00677 else
00678 return _ch[i-1]-_ch[i-2];
00679 }
00680
00681
00682 template<class T_>
00683 Vect<T_> SkSMatrix<T_>::getColumn(size_t j) const
00684 {
00685 Vect<T_> v(_nb_rows);
00686 for (size_t i=1; i<=_nb_rows; i++)
00687 v(i) = (*this)(i,j);
00688 return v;
00689 }
00690
00691
00692 template<class T_>
00693 Vect<T_> SkSMatrix<T_>::getRow(size_t i) const
00694 {
00695 return getColumn(i);
00696 }
00697
00698
00699 template<class T_>
00700 T_ & SkSMatrix<T_>::operator()(size_t i, size_t j)
00701 {
00702 #ifdef _BOUNDS_DEBUG
00703 assert(i>0);
00704 assert(i<=_size);
00705 assert(j>0);
00706 assert(j<=_size);
00707 #endif
00708 int k=0, l=0;
00709 if (i>1)
00710 k = int(j-i+_ch[i-1]-_ch[i-2]-1);
00711 if (j>1)
00712 l = int(i-j+_ch[j-1]-_ch[j-2]-1);
00713 if (k>=0 && i>=j)
00714 return _a[_ch[i-1]+j-i];
00715 else if (l>=0 && i<j)
00716 return _a[_ch[j-1]+i-j];
00717 else
00718 return _zero;
00719 }
00720
00721
00722 template<class T_>
00723 T_ SkSMatrix<T_>::operator()(size_t i, size_t j) const
00724 {
00725 #ifdef _BOUNDS_DEBUG
00726 assert(i>0);
00727 assert(i<=_size);
00728 assert(j>0);
00729 assert(j<=_size);
00730 #endif
00731 int k=0, l=0;
00732 if (i>1)
00733 k = int(j-i+_ch[i-1]-_ch[i-2]-1);
00734 if (j>1)
00735 l = int(i-j+_ch[j-1]-_ch[j-2]-1);
00736 if (k>=0 && i>=j)
00737 return _a[_ch[i-1]+j-i];
00738 else if (l>=0 && i<j)
00739 return _a[_ch[j-1]+i-j];
00740 else
00741 return _zero;
00742 }
00743
00744
00745 template<class T_>
00746 SkSMatrix<T_> & SkSMatrix<T_>::operator=(const SkSMatrix<T_> &m)
00747 {
00748 _a = m._a;
00749 return *this;
00750 }
00751
00752
00753 template<class T_>
00754 SkSMatrix<T_> & SkSMatrix<T_>::operator=(const T_ &x)
00755 {
00756 for (size_t i=0; i<_length; i++)
00757 _a[i] = x;
00758 return *this;
00759 }
00760
00761
00762 template<class T_>
00763 SkSMatrix<T_> & SkSMatrix<T_>::operator+=(const SkSMatrix<T_> &m)
00764 {
00765 #ifdef _BOUNDS_DEBUG
00766 assert(m.Dim() >= _size);
00767 #endif
00768 for (size_t i=0; i<_length; i++)
00769 _a[i] += m._a[i];
00770 return *this;
00771 }
00772
00773
00774 template<class T_>
00775 SkSMatrix<T_> & SkSMatrix<T_>::operator*=(const T_ &x)
00776 {
00777 for (size_t i=0; i<_length; i++)
00778 _a[i] *= x;
00779 return *this;
00780 }
00781
00782
00783 template<class T_>
00784 int SkSMatrix<T_>::Factor()
00785 {
00786 size_t di=0, dij;
00787 T_ s;
00788 T_ pivot = _a[_ch[0]];
00789 if (Abs(pivot) < OFELI_TOLERANCE)
00790 _e.Message(__FILE__,__LINE__,61,1,0,pivot);
00791 _a[_ch[0]] = T_(1.)/pivot;
00792
00793 for (size_t i=1; i<_size; i++) {
00794 size_t dj = 0;
00795 for (size_t j=di=i+1+_ch[i-1]-_ch[i]; j<i; j++) {
00796 if (j>0)
00797 dj = j+1+_ch[j-1]-_ch[j];
00798 dij = max(di,dj);
00799 for (size_t k=0; k<j-dij; k++)
00800 _a[_ch[i]+j-i] -= _a[_ch[i]+dij+k-i]*_a[_ch[j]+dij+k-j];
00801 }
00802
00803 pivot = _a[_ch[i]];
00804 for (size_t k=di; k<i; k++) {
00805 s = _a[_ch[i]+k-i]*_a[_ch[k]];
00806 pivot -= s*_a[_ch[i]+k-i];
00807 _a[_ch[i]+k-i] = s;
00808 }
00809 if (Abs(pivot) < OFELI_TOLERANCE)
00810 _e.Message(__FILE__,__LINE__,61,int(i+1),0,pivot);
00811 _a[_ch[i]] = T_(1.)/pivot;
00812 }
00813 _fact = true;
00814 return 0;
00815 }
00816
00817
00818 template<class T_>
00819 int SkSMatrix<T_>::Solve(Vect<T_> &b)
00820 {
00821 if (!_fact)
00822 _e.Message(__FILE__,__LINE__,71);
00823 size_t i, j, di;
00824 for (i=1; i<_size; i++) {
00825 di = i+1+_ch[i-1]-_ch[i];
00826 T_ s = 0;
00827 for (j=0; j<i-di; j++)
00828 s += _a[_ch[i]+di+j-i] * b[di+j];
00829 b[i] -= s;
00830 }
00831 for (i=0; i<_size; i++)
00832 b[i] *= _a[_ch[i]];
00833 for (int k=int(_size-1); k>0; k--) {
00834 di = k+1+_ch[k-1]-_ch[k];
00835 for (j=0; j<k-di; j++)
00836 b[j+di] -= b[k] * _a[_ch[k]+di+j-k];
00837 }
00838 return 0;
00839 }
00840
00841
00842 template<class T_>
00843 int SkSMatrix<T_>::Solve(const Vect<T_> &b, Vect<T_> &x)
00844 {
00845 x = b;
00846 return Solve(x);
00847 }
00848
00849
00850 template<class T_>
00851 T_ SkSMatrix<T_>::getEntry(size_t i, size_t j) const
00852 {
00853 #ifdef _BOUNDS_DEBUG
00854 assert(i>0);
00855 assert(i<=_size);
00856 assert(j>0);
00857 assert(j<=_size);
00858 #endif
00859 if (i>=j)
00860 return _a[_ch[i-1]+j-i];
00861 else
00862 return _a[_ch[j-1]+i-j];
00863 }
00864
00866
00868
00872 template<class T_>
00873 ostream& operator<<(ostream& s, const SkSMatrix<T_> & a)
00874 {
00875 size_t rmin, rmax, cmin, cmax;
00876 a.getPrintView(rmin,rmax,cmin,cmax);
00877 s.setf(ios::right|ios::scientific);
00878 s << endl;
00879 for (size_t i=rmin; i<=rmax; i++) {
00880 s << "\nRow: " << setw(6) << i << endl;
00881 for (size_t j=cmin; j<=cmax; j++)
00882 s << " " << setprecision(8) << std::setfill(' ') << setw(18) << a(i,j);
00883 s << endl;
00884 }
00885 return s;
00886 }
00887
00888 }
00889
00890 #endif
