xc
matdispZ.h
1 // -*-c++-*-
2 //----------------------------------------------------------------------------
3 // xc utils library; general purpose classes and functions.
4 //
5 // Copyright (C) Luis C. Pérez Tato
6 //
7 // XC utils is free software: you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation, either version 3 of the License, or
10 // (at your option) any later version.
11 //
12 // This software is distributed in the hope that it will be useful, but
13 // WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
16 //
17 // You should have received a copy of the GNU General Public License
18 // along with this program.
19 // If not, see <http://www.gnu.org/licenses/>.
20 //----------------------------------------------------------------------------
21 //matdispZ.h
22 //Matrices dispersas.
23 
24 #ifndef MATDISPZ_H
25 #define MATDISPZ_H
26 
27 #include "ZMatrix.h"
28 #include <map>
29 #include <deque>
30 #include <algorithm>
31 
32 template <class numero>
33 class matdispZ : public ZMatrix<numero>
34  {
35  public:
36  typedef ZMatrix<numero> ZMatrix_number;
37  typedef typename ZMatrix_number::size_type size_type;
38  typedef std::map<size_type,numero> map_elem;
39 
40  class sp_vector: public map_elem
41  {
42  public:
43  typedef typename map_elem::iterator iterator;
44  typedef typename map_elem::const_iterator const_iterator;
45  typedef typename map_elem::const_reverse_iterator const_reverse_iterator;
46  typedef typename map_elem::key_type key_type;
47  private:
48  inline int hasRow(const typename sp_vector::size_type f) const
49  { return (find(f)!=this->end()); }
50  public:
51  sp_vector &operator+=(const sp_vector &v);
52  sp_vector &operator-=(const sp_vector &v);
53  void QuitaElem(const numero &n);
54  void PutCol(const typename sp_vector::size_type c,ZMatrix_number &m) const;
55  sp_vector getNumberOfRows(const typename sp_vector::size_type f1,const typename sp_vector::size_type f2) const;
56  size_t ndiagL(const size_t &icol) const;
57  size_t ndiagU(const size_t &icol) const;
58  void writeCpp(std::ostream &os,const size_t &icol) const;
59  void PutColBanda(const size_t &sz,const size_t &i,const size_t &ndiagu,numero *vptr) const;
60  };
61  typedef std::map<size_type, sp_vector> column_map;
62  column_map columns;
63 
64 
65  typedef typename sp_vector::const_iterator const_f_iterator;
66  typedef typename sp_vector::iterator f_iterator;
67  typedef typename column_map::const_iterator const_c_iterator;
68  typedef typename column_map::iterator c_iterator;
69 
70  private:
71  inline virtual size_type Tam(size_type ,size_type )
72  { return 0; }
73  inline virtual size_type Indice(const size_type &iRow,const size_type &col) const
74  { return 0; }
75  inline bool hasColumn(const size_type &c) const
76  { return (columns.find(c)!=columns.end()); }
77  inline const numero &PorDefecto(void) const
78  { return (ZMatrix_number::operator()(1,1)); }
79  bool EqualTo(const matdispZ<numero> &other) const;
80 
81  template<class M>
82  matdispZ<numero> Post(const M &b) const;
83 
84  matdispZ<numero> Post(const matdispZ<numero> &b) const;
85 
86  public:
87 
88  matdispZ(size_type n_rows=1,size_type n_columns= 1)
89  : ZMatrix_number(1,1,numero())
90  { this->PutDim(n_rows,n_columns); }
91  matdispZ(const matdispZ<numero> &other)
92  : ZMatrix_number(other), columns(other.columns)
93  { this->PutDim(other.n_rows,other.n_columns); }
94  matdispZ<numero>& operator =(const matdispZ<numero> &m)
95  {
96  ZMatrix_number::operator =(m);
97  columns= m.columns;
98  return *this;
99  }
100  matdispZ<numero>& operator+=(const matdispZ<numero> &m);
101  matdispZ<numero>& operator-=(const matdispZ<numero> &m);
102 
103 
104  size_t ndiagL(void) const;
105  size_t ndiagU(void) const;
106  void FillVectorBanda(numero *vptr) const;
107  //Access to matrix elements.
108  numero &operator()(size_t iRow=1,size_t col=1)
109  {
110  assert(col>0 && col<=this->n_columns);
111  assert(iRow>0 && iRow<=this->n_rows);
112  sp_vector &column= columns[col];
113  numero &retval= column[iRow];
114  return retval;
115  }
116  const numero &operator()(size_t iRow=1,size_t col=1) const
117  {
118  const_c_iterator k1= columns.find(col);
119  if(k1!=columns.end())
120  {
121  const_f_iterator k2= (k1->second).find(iRow);
122  if(k2!= (k1->second).end())
123  return k2->second;
124  else
125  return PorDefecto();
126  }
127  else
128  return PorDefecto();
129  }
130  const_c_iterator columns_begin(void) const
131  { return columns.begin(); }
132  const_c_iterator columns_end(void) const
133  { return columns.end(); }
134  c_iterator columns_begin(void)
135  { return columns.begin(); }
136  c_iterator columns_end(void)
137  { return columns.end(); }
138  const_f_iterator rows_begin(const const_c_iterator &ci) const
139  { return ci->second.begin(); }
140  const_f_iterator find_row(const const_c_iterator &ci,size_t f) const
141  { return ci->second.find(f); }
142  const_f_iterator rows_end(const const_c_iterator &ci) const
143  { return ci->second.end(); }
144  f_iterator find_row(const c_iterator &ci,size_t f)
145  { return ci->second.find(f); }
146  sp_vector &getColumn(const c_iterator &ci)
147  { return ci->second; }
148 
149  size_type size(void) const;
150  void Compacta(void);
151  void Anula(void);
152  void Identity(void);
153  inline bool Cuadrada(void) const
154  { return (this->n_rows == this->n_columns); }
155  matdispZ<numero> &Trn(void);
156  matdispZ<numero> GetTrn(void)
157  {
158  matdispZ<numero> retval(*this);
159  retval.Trn();
160  return retval;
161  }
162  matdispZ<numero> getBox(size_t f1, size_t c1, size_t f2, size_t c2) const;
163  matdispZ<numero> getRow(size_t iRow) const
164  { return getBox(iRow,1,iRow,this->n_columns); }
165  matdispZ<numero> getColumn(size_t col) const
166  { return getBox(1,col,this->n_rows,col); }
167  void putBox(size_t f,size_t c,const matdispZ<numero> &box);
168  numero Traza(void) const;
169  ZMatrix_number GetCompleta(void) const;
170 
171  void writeCpp(std::ostream &os) const;
172 
173  template <class V>
174  numero dot(const V &v2) const
175  //Producto escalar de este por v2.
176  //v2: column vector.
177  {
178  check_dot(*this,v2);
179  return (this->Post(v2))(1,1);
180  }
181  numero dot(const matdispZ<numero> &v2) const
182  //Producto escalar de este por v2.
183  //v2: column vector.
184  {
185  check_dot(*this,v2);
186  return (this->Post(v2))(1,1);
187  }
188  friend matdispZ<numero> operator+(const matdispZ<numero> &m1,const matdispZ<numero> &m2)
189  {
190  check_sum(m1,m2);
191  matdispZ<numero> suma(m1);
192  suma+= m2;
193  suma.Compacta();
194  return suma;
195  }
196  friend matdispZ<numero> operator-(const matdispZ<numero> &m1,const matdispZ<numero> &m2)
197  {
198  check_dif(m1,m2);
199  matdispZ<numero> resta(m1);
200  resta-=m2;
201  resta.Compacta();
202  return resta;
203  }
204  inline friend bool operator==(const matdispZ<numero> &m1,const matdispZ<numero> &m2)
205  { return m1.EqualTo(m2); }
206  //Producto escalar de dos vectores.
207  //v1: row vector.
208  //v2: column vector.
209  friend numero dot(const matdispZ<numero> &v1,const matdispZ<numero> &v2)
210  { return v1.dot(v2); }
211  friend numero dot(const matdispZ<numero> &v1,const ZMatrix<numero> &v2)
212  { return v1.dot(v2); }
213  friend matdispZ<numero> operator*(const matdispZ<numero> &m1,const matdispZ<numero> &m2)
214  {
215  check_prod(m1,m2);
216  matdispZ<numero> producto(m1.Post(m2));
217  return producto;
218  }
219  friend matdispZ<numero> operator*(const matdispZ<numero> &m1,const ZMatrix<numero> &m2)
220  {
221  check_prod(m1,m2);
222  matdispZ<numero> producto(m1.Post(m2));
223  return producto;
224  }
225  /*
226  friend matdispZ<numero> operator*(const matdispZ<numero> &m1,const ZMatrix<numero> &m2)
227  {
228  check_prod(m1,m2);
229  matdispZ<numero> producto(m1.n_rows,m2.getNumberOfColumns());
230  matdispZ<numero>::size_type i=1,j=1;
231  for (i= 1;i<=m1.n_rows;i++)
232  for (j= 1;j<=m2.getNumberOfColumns();j++)
233  producto(i,j) = dot(m1.getRow(i),m2.getColumn(j));
234  return producto;
235  }
236  */
237  /*
238  friend matdispZ<numero> operator *(const matdispZ<numero> &m,const numero &p)
239  {
240  matdispZ<numero> producto(m);
241  producto.Prod(p);
242  return producto;
243  }
244  friend matdispZ<numero> operator *(const numero &p,const matdispZ<numero> &m)
245  { return m*p; }
246 
247  friend std::istream &operator >> (std::istream &stream,matdispZ<numero> &m)
248  { return ::operator >>(stream,(ZMatrix_number &) m); }
249  */
250  };
251 
252 //sp_vector
253 template<class numero>
254 typename matdispZ<numero>::sp_vector &matdispZ<numero>::sp_vector::operator+=(const matdispZ<numero>::sp_vector &v)
255  {
256  const_iterator f;
257  for(f= v.begin();f!=v.end();f++)
258  if(hasRow(f->first))
259  (*this)[f->first]+= f->second;
260  else
261  (*this)[f->first]= f->second;
262  return *this;
263  }
264 template<class numero>
265 typename matdispZ<numero>::sp_vector &matdispZ<numero>::sp_vector::operator-=(const matdispZ<numero>::sp_vector &v)
266  {
267  const_iterator f;
268  for(f= v.begin();f!=v.end();f++)
269  if(hasRow(f->first))
270  (*this)[f->first]-= f->second;
271  else
272  (*this)[f->first]= f->second;
273  return *this;
274  }
275 template<class numero>
276 void matdispZ<numero>::sp_vector::QuitaElem(const numero &n)
277  {
278  std::deque<key_type> markedForDeath;
279  for(const_iterator f= this->begin();f!=this->end();f++)
280  if(f->second == n)
281  markedForDeath.push_back(f->first);
282  for(auto i= markedForDeath.begin();i!=markedForDeath.end();i++)
283  this->erase(*i);
284  }
285 template<class numero>
286 void matdispZ<numero>::sp_vector::PutCol(const typename sp_vector::size_type c,ZMatrix_number &m) const
287  {
288  const_iterator f;
289  for(f= this->begin();f!=this->end();f++)
290  m(f->first,c)= f->second;
291  }
292 template<class numero>
293 typename matdispZ<numero>::sp_vector matdispZ<numero>::sp_vector::getNumberOfRows(const typename sp_vector::size_type f1,const typename sp_vector::size_type f2) const
294  {
295  sp_vector retval;
296  if(f2<f1) return retval;
297  const_iterator f;
298  for(f= this->begin();f!=this->end();f++)
299  if((f->first>=f1) && (f->first<=f2))
300  retval[f->first-f1+1]= f->second;
301  return retval;
302  }
303 
306 template<class numero>
307 size_t matdispZ<numero>::sp_vector::ndiagL(const size_t &icol) const
308  {
309  size_t retval= 0;
310  if(!this->empty())
311  {
312  const numero zero= numero();
313  for(const_reverse_iterator f= this->rbegin();f!=this->rend();f++)
314  {
315  const size_t iRow= f->first;
316  if(iRow>icol) //El elemento esta por debajo de la diagonal principal.
317  {
318  if(f->second != zero) //no nulo.
319  {
320  retval= iRow-icol;
321  break;
322  }
323  }
324  else //El elemento esta en la diagonal principal o por encima.
325  break;
326  }
327  }
328  return retval;
329  }
330 
333 template<class numero>
334 size_t matdispZ<numero>::sp_vector::ndiagU(const size_t &icol) const
335  {
336  size_t retval= 0;
337  if(!this->empty())
338  {
339  const numero zero= numero();
340  for(const_iterator f= this->begin();f!=this->end();f++)
341  {
342  const size_t iRow= f->first;
343  if(iRow<icol) //El elemento esta por encima de la diagonal principal.
344  {
345  if(f->second != zero) //no nulo.
346  {
347  retval= icol-iRow;
348  break;
349  }
350  }
351  else //El elemento esta en la diagonal principal o por debajo.
352  break;
353  }
354  }
355  return retval;
356  }
357 
359 template<class numero>
360 void matdispZ<numero>::sp_vector::writeCpp(std::ostream &os,const size_t &icol) const
361  {
362  if(!this->empty())
363  {
364  const numero zero= numero();
365  for(const_iterator f= this->begin();f!=this->end();f++)
366  if(f->second != zero)
367  os << "mtx(" << f->first << ',' << icol << ")= " << f->second << "; ";
368  }
369  }
370 
377 template<class numero>
378 void matdispZ<numero>::sp_vector::PutColBanda(const size_t &sz,const size_t &icol,const size_t &ndiagu,numero *vptr) const
379  {
380  if(!this->empty())
381  {
382  const numero zero= numero();
383  const int offset= (icol-1)*sz-icol+ndiagu;
384  for(const_iterator f= this->begin();f!=this->end();f++)
385  {
386  const size_t iRow_banda= f->first+offset;
387  if(f->second != zero)
388  vptr[iRow_banda]= f->second;
389  }
390  }
391  }
392 
393 //matdispZ
394 template <class numero>
395 matdispZ<numero>& matdispZ<numero>::operator+=(const matdispZ<numero> &m)
396  {
397  const_c_iterator c;
398  for(c= m.columns.begin();c!=m.columns.end();c++)
399  if(hasColumn(c->first))
400  columns[c->first]+= c->second;
401  else
402  columns[c->first]= c->second;
403  return *this;
404  }
405 
406 template <class numero>
407 matdispZ<numero>& matdispZ<numero>::operator-=(const matdispZ<numero> &m)
408  {
409  const_c_iterator c;
410  for(c= m.columns.begin();c!=m.columns.end();c++)
411  if(hasColumn(c->first))
412  columns[c->first]-= c->second;
413  else
414  columns[c->first]= c->second;
415  return *this;
416  }
417 
418 template <class numero>
419 typename matdispZ<numero>::size_type matdispZ<numero>::size(void) const
420  {
421  size_type retval= 0;
422  const_c_iterator c;
423  for(c= columns.begin();c!=columns.end();c++)
424  retval+=c->second.size();
425  return retval;
426  }
427 
428 template <class numero>
429 void matdispZ<numero>::Compacta(void)
430  {
431  const numero &n= PorDefecto();
432  typename column_map::iterator c;
433  for(c= columns.begin();c!=columns.end();c++)
434  c->second.QuitaElem(n);
435  }
436 
437 template <class numero>
438 void matdispZ<numero>::Anula(void)
439  { columns.erase(columns.begin(),columns.end()); }
440 
441 template <class numero>
442 void matdispZ<numero>::Identity(void)
443  {
444  this->Anula();
445  if(this->n_rows!=this->n_columns)
446  std::cerr << "matdispZ::" << __FUNCTION__
447  << "not a square matrix: "
448  << this->n_rows << " x " << this->n_columns << std::endl;
449  const size_t sz= std::min(this->n_rows,this->n_columns);
450  for(size_t i=1;i<=sz;i++)
451  (*this)(i,i)= 1.0;
452  Compacta();
453  }
454 
455 template <class numero>
456 matdispZ<numero> &matdispZ<numero>::Trn(void)
457  {
458  ZMatrix_number::Trn();
459  const_c_iterator c;
460  const_f_iterator f;
461  for(c= columns.begin();c!=columns.end();c++)
462  for(f= c->second.find(c->first+1);f!=c->second.end();f++)
463  std::swap((*this)(f->first,c->first),(*this)(c->first,f->first));
464  Compacta();
465  return *this;
466  }
467 
468 template <class numero>
469 size_t matdispZ<numero>::ndiagL(void) const
470  {
471  const_c_iterator c;
472  size_t retval= 0;
473  for(c= columns.begin();c!=columns.end();c++)
474  retval= std::max(retval,c->second.ndiagL(c->first));
475  return retval;
476  }
477 template <class numero>
478 size_t matdispZ<numero>::ndiagU(void) const
479  {
480  const_c_iterator c;
481  size_t retval= 0;
482  for(c= columns.begin();c!=columns.end();c++)
483  retval= std::max(retval,c->second.ndiagU(c->first));
484  return retval;
485  }
486 
488 template<class numero>
489 void matdispZ<numero>::FillVectorBanda(numero *vptr) const
490  {
491  const size_t ndiagu= ndiagU();
492  const size_t ndiagl= ndiagL();
493  const size_t ancho_banda= ndiagl+ndiagu+1;
494  for(const_c_iterator c= columns.begin();c!=columns.end();c++)
495  c->second.PutColBanda(ancho_banda,c->first,ndiagu,vptr);
496  }
497 
499 template<class numero>
500 void matdispZ<numero>::writeCpp(std::ostream &os) const
501  {
502  for(const_c_iterator c= columns.begin();c!=columns.end();c++)
503  c->second.writeCpp(os,c->first);
504  }
505 
506 
507 template <class numero>
508 matdispZ<numero> matdispZ<numero>::getBox(size_t f1, size_t c1, size_t f2, size_t c2) const
509  {
510  // XXX "mejorar esta rutina"
511  this->check_get_box(f1,c1,f2,c2);
512  matdispZ<numero> box(f2-f1+1,c2-c1+1);
513  const_c_iterator c;
514  const_f_iterator f;
515  for(c= columns.begin();c!=columns.end();c++)
516  for(f= c->second.begin();f!=c->second.end();f++)
517  if( (c->first>=c1) && (c->first<=c2) &&
518  (f->first>=f1) && (f->first<=f2))
519  box(f->first-f1+1,c->first-c1+1)= f->second;
520  return box;
521  }
522 
523 template <class numero>
524 void matdispZ<numero>::putBox(size_t f,size_t c,const matdispZ<numero> &box)
525  {
526  check_put_box(f,c,box);
527  size_t i,j;
528  for (i=1;i<=box.n_rows;i++)
529  for (j=1;j<=box.n_columns;j++)
530  (*this)(i+f-1,j+c-1)= box(i,j);
531  Compacta();
532  }
533 
534 template<class numero>
535 numero matdispZ<numero>::Traza(void) const
536  {
537  this->check_traza();
538  numero n= numero();
539  const_c_iterator c;
540  const_f_iterator f;
541  for(c= columns.begin();c!=columns.end();c++)
542  if((f= c->second.find(c->first)) != c->second.end())
543  n+=f->second;
544  return n;
545  }
546 
547 template<class numero>
548 typename matdispZ<numero>::ZMatrix_number matdispZ<numero>::GetCompleta(void) const
549  {
550  ZMatrix_number retval(this->n_rows,this->n_columns,PorDefecto());
551  const_c_iterator c;
552  for(c= columns.begin();c!=columns.end();c++)
553  c->second.PutCol(c->first,retval);
554  return retval;
555  }
556 
557 template<class numero>
558 bool matdispZ<numero>::EqualTo(const matdispZ<numero> &other) const
559  {
560  if(!CompDim(*this,other)) return false;
561  typename matdispZ<numero>::const_c_iterator c;
562  for(c= this->columns.begin();c!=this->columns.end();c++)
563  {
564  typename matdispZ<numero>::const_f_iterator f;
565  for(f= c->second.begin();f!=c->second.end();f++)
566  if(f->second != other(f->first,c->first)) return false;
567  }
568  return true;
569  }
570 
571 template<class numero>
572 //Return the product of this matrix by the matrix argument.
573 matdispZ<numero> matdispZ<numero>::Post(const matdispZ<numero> &b) const
574  {
575  check_prod(*this,b);
576  matdispZ<numero> ret(this->getNumberOfRows(),b.getNumberOfColumns());
577  const_c_iterator bc;
578  const_c_iterator c;
579  const_f_iterator f;
580  for(bc= b.columns.begin();bc!=b.columns.end();bc++)
581  for(c= columns.begin();c!=columns.end();c++)
582  for(f= c->second.begin();f!=c->second.end();f++)
583  ret(f->first,bc->first)+= f->second*b(c->first,bc->first);
584  return(ret);
585  }
586 
587 template<class numero>
588 template<class M>
589 //Return the product of this matrix by the matrix argument.
590 matdispZ<numero> matdispZ<numero>::Post(const M &b) const
591  {
592  check_prod(*this,b);
593  matdispZ<numero> ret(this->getNumberOfRows(),b.getNumberOfColumns());
594  size_type bc;
595  const_c_iterator c;
596  const_f_iterator f;
597  for(bc= 1;bc<=b.getNumberOfColumns();bc++)
598  for(c= columns.begin();c!=columns.end();c++)
599  for(f= c->second.begin();f!=c->second.end();f++)
600  ret(f->first,bc)+= f->second*b(c->first,bc);
601  return(ret);
602  }
603 
604 
605 #endif
matdispZ
Definition: matdispZ.h:33
ZMatrix
Matrix which element type has estructura de anillo respecto a las operaciones + y *...
Definition: ZMatrix.h:37
matdispZ::writeCpp
void writeCpp(std::ostream &os) const
Writes the matrix in C++ format (only non-zero components).
Definition: matdispZ.h:500
matdispZ::sp_vector::PutColBanda
void PutColBanda(const size_t &sz, const size_t &i, const size_t &ndiagu, numero *vptr) const
Coloca los elementos de the column que forman parte de la banda en el vector que is being passed as p...
Definition: matdispZ.h:378
matdispZ::sp_vector::writeCpp
void writeCpp(std::ostream &os, const size_t &icol) const
Escribe los elementos no nulos de the column en formato de C++.
Definition: matdispZ.h:360
matdispZ::sp_vector::ndiagL
size_t ndiagL(const size_t &icol) const
Return el número de diagonales, con algún elemento no nulo, que tiene the column por debajo de la dia...
Definition: matdispZ.h:307
matdispZ::FillVectorBanda
void FillVectorBanda(numero *vptr) const
Rellena el vector en banda que is being passed as parameter para su empleo en Arpack++.
Definition: matdispZ.h:489
matdispZ::sp_vector::ndiagU
size_t ndiagU(const size_t &icol) const
Return el número de diagonales, con algún elemento no nulo, que tiene the column por encima de la dia...
Definition: matdispZ.h:334
matdispZ::sp_vector
Definition: matdispZ.h:40
Generated by   doxygen 1.8.13