← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty

« back to all changes in this revision

Viewing changes to skit/plib2/csr_seq.cc

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito
  • Date: 2012-04-06 09:12:21 UTC
  • mfrom: (1.1.5)
  • Revision ID: package-import@ubuntu.com-20120406091221-m58me99p1nxqui49
Tags: 6.0-1
http://bugs.debian.org/667356
http://bugs.debian.org/661689
* New upstream release 6.0 (major changes):
  - massively distributed and parallel support
  - full FEM characteristic method (Lagrange-Gakerkin method) support
  - enhanced users documentation 
  - source code supports g++-4.7 (closes: #667356)
* debian/control: dependencies for MPI distributed solvers added
* debian/rules: build commands simplified
* debian/librheolef-dev.install: man1/* to man9/* added
* debian/changelog: package description rewritted (closes: #661689)

Show diffs side-by-side

added added

removed removed

Lines of Context:
skit/plib2/csr_seq.cc
26
26
#include "rheolef/asr_to_csr.h"
27
27
#include "rheolef/csr_to_asr.h"
28
28
#include "rheolef/csr_amux.h"
 
29
#include "rheolef/csr_cumul_trans_mult.h"
29
30
using namespace std;
30
31
namespace rheolef {
31
32
// ----------------------------------------------------------------------------
50
51
   _row_ownership = row_ownership;
51
52
   _col_ownership = col_ownership;
52
53
   _data.resize (nnz1);
 
54
   // first pointer points to the beginning of the data:
 
55
   vector_of_iterator<pair_type>::operator[](0) = _data.begin().operator->();
53
56
}
54
57
template<class T>
55
58
csr_seq_rep<T>::csr_seq_rep(size_type loc_nrow1, size_type loc_ncol1, size_type loc_nnz1) 
69
72
   _row_ownership = distributor (distributor::decide, communicator(), loc_nrow1);
70
73
   _col_ownership = distributor (distributor::decide, communicator(), loc_ncol1);
71
74
   _data.resize (loc_nnz1);
 
75
   // first pointer points to the beginning of the data:
 
76
   vector_of_iterator<pair_type>::operator[](0) = _data.begin().operator->();
72
77
}
73
78
template<class T>
74
 
csr_seq_rep<T>::csr_seq_rep(const csr_seq_rep<T>& a)
75
 
 : vector_of_iterator<pair_type>(a.nrow()+1),
76
 
   _row_ownership (a.row_ownership()),
77
 
   _col_ownership (a.col_ownership()),
78
 
   _data(a._data),
79
 
   _is_symmetric (a._is_symmetric),
80
 
   _pattern_dimension (a._pattern_dimension)
 
79
csr_seq_rep<T>::csr_seq_rep(const csr_seq_rep<T>& b)
 
80
 : vector_of_iterator<pair_type>(b.nrow()+1),
 
81
   _row_ownership (b.row_ownership()),
 
82
   _col_ownership (b.col_ownership()),
 
83
   _data(b._data),
 
84
   _is_symmetric (b._is_symmetric),
 
85
   _pattern_dimension (b._pattern_dimension)
81
86
{
82
 
   // TODO: performs copy of arrays
83
 
   fatal_macro ("physical copy of csr");
84
 
}     
 
87
  // physical copy of csr
 
88
  typedef typename csr_seq_rep<T>::size_type size_type;
 
89
  typename csr_seq_rep<T>::const_iterator ib = b.begin();
 
90
  typename csr_seq_rep<T>::iterator       ia = begin();
 
91
  ia[0] = _data.begin().operator->();
 
92
  for (size_type i = 0, n = b.nrow(); i < n; i++) {
 
93
    ia [i+1] = ia[0] + (ib[i+1] - ib[0]);
 
94
  }
 
95
}
85
96
template<class T>
86
97
csr_seq_rep<T>::csr_seq_rep(const asr_seq_rep<T>& a)
87
98
  : vector_of_iterator<pair_type>(a.nrow()+1),
137
148
csr_seq_rep<T>::put (odiststream& ops, size_type istart) const
138
149
{
139
150
    std::ostream& os = ops.os();
140
 
    os << "%%MatrixMarket matrix coordinate real general" << std::endl
 
151
    os << setprecision (std::numeric_limits<T>::digits10)
 
152
       << "%%MatrixMarket matrix coordinate real general" << std::endl
141
153
       << nrow() << " " << ncol() << " " << nnz() << endl;
142
154
    const_iterator ia = begin();
143
155
    const size_type base = 1;
179
191
        set_op<T,T>(),
180
192
        y.begin());
181
193
}
 
194
template<class T>
 
195
void
 
196
csr_seq_rep<T>::trans_mult(
 
197
    const vec<T,sequential>& x,
 
198
    vec<T,sequential>&       y)
 
199
    const
 
200
{
 
201
    // reset y and then cumul
 
202
    std::fill (y.begin(), y.end(), T(0));
 
203
    csr_cumul_trans_mult (
 
204
        vector_of_iterator<pair_type>::begin(), 
 
205
        vector_of_iterator<pair_type>::end(), 
 
206
        x.begin(), 
 
207
        set_add_op<T,T>(),
 
208
        y.begin());
 
209
}
 
210
template<class T>
 
211
csr_seq_rep<T>&
 
212
csr_seq_rep<T>::operator*= (const T& lambda)
 
213
{
 
214
  iterator ia = begin();
 
215
  for (data_iterator p = ia[0], last_p = ia[nrow()]; p != last_p; p++) {
 
216
    (*p).second *= lambda;
 
217
  }
 
218
  return *this;
 
219
}
182
220
// ----------------------------------------------------------------------------
183
221
// expression c=a+b and c=a-b with a temporary c=*this
184
222
// ----------------------------------------------------------------------------
248
286
    set_pattern_dimension (std::max(a.pattern_dimension(), b.pattern_dimension()));
249
287
}
250
288
// ----------------------------------------------------------------------------
 
289
// trans(a)
 
290
// ----------------------------------------------------------------------------
 
291
/*@! 
 
292
 \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
 
293
  \caption{{\tt sort}: sort rows by increasing column order}
 
294
  \begin{algorithmic}
 
295
    \INPUT {the matrix in CSR format}
 
296
      ia(0:nrow), ja(0:nnz-1), a(0:nnz-1)
 
297
    \ENDINPUT
 
298
    \OUTPUT {the sorted CSR matrix}
 
299
      iao(0:nrow), jao(0:nnzl-1), ao(0:nnzl-1)
 
300
    \ENDOUTPUT
 
301
    \BEGIN 
 
302
      {\em first: reset iao} \\
 
303
      \FORTO {i := 0} {nrow}
 
304
          iao(i) := 0;
 
305
      \ENDFOR
 
306
        
 
307
      {\em second: compute lengths from pointers} \\
 
308
      \FORTO {i := 0} {nrow-1}
 
309
        \FORTO {p := ia(i)} {ia(i+1)-1}
 
310
            iao (ja(p)+1)++;
 
311
        \ENDFOR
 
312
      \ENDFOR
 
313
 
 
314
      {\em third: compute pointers from lengths} \\
 
315
      \FORTO {i := 0} {nrow-1}
 
316
          iao(i+1) += iao(i)
 
317
      \ENDFOR
 
318
 
 
319
      {\em fourth: copy values} \\
 
320
      \FORTO {i := 0} {nrow-1}
 
321
        \FORTO {p := ia(i)} {ia(i+1)-1}
 
322
          j := ja(p) \\
 
323
          q := iao(j) \\
 
324
          jao(q) := i \\
 
325
          ao(q) := a(p) \\
 
326
          iao (j)++
 
327
        \ENDFOR
 
328
      \ENDFOR
 
329
 
 
330
      {\em fiveth: reshift pointers} \\
 
331
      \FORTOSTEP {i := nrow-1} {0} {-1}
 
332
        iao(i+1) := iao(i);
 
333
      \ENDFOR
 
334
      iao(0) := 0
 
335
    \END
 
336
 \end{algorithmic} \end{algorithm}
 
337
 \vfill \pagebreak \mbox{} \vfill
 
338
*/
 
339
 
 
340
// NOTE: transposed matrix has always rows sorted by increasing column indexes
 
341
//       even if original matrix has not
 
342
template<class T>
 
343
void
 
344
csr_seq_rep<T>::build_transpose (csr_seq_rep<T>& b) const
 
345
{
 
346
  b.resize (col_ownership(), row_ownership(), nnz()); 
 
347
 
 
348
  // first pass: set ib(*) to ib(0)
 
349
  iterator ib = b.begin();
 
350
  for (size_type j = 0, m = b.nrow(); j < m; j++) {
 
351
    ib[j+1] = ib[0];
 
352
  }
 
353
  // second pass: compute lengths of row(i) of a^T in ib(i+1)-ib(0)
 
354
  const_iterator ia = begin();
 
355
  for (size_type i = 0, n = nrow(); i < n; i++) {
 
356
    for (const_data_iterator p = ia[i], last_p = ia[i+1]; p != last_p; p++) {
 
357
      size_type j = (*p).first;
 
358
      ib [j+1]++;
 
359
    }
 
360
  }
 
361
  // third pass: compute pointers from lengths
 
362
  for (size_type j = 0, m = b.nrow(); j < m; j++) {
 
363
    ib [j+1] += (ib[j]-ib[0]);
 
364
  }
 
365
  // fourth pass: store values
 
366
  data_iterator q0 = ib[0];
 
367
  for (size_type i = 0, n = nrow(); i < n; i++) {
 
368
    for (const_data_iterator p = ia[i], last_p = ia[i+1]; p != last_p; p++) {
 
369
      size_type j = (*p).first;
 
370
      data_iterator q = ib[j];
 
371
      (*q).first  = i;
 
372
      (*q).second = (*p).second;
 
373
      ib[j]++;
 
374
    }
 
375
  }
 
376
  // fiveth: shift pointers
 
377
  for (long int j = b.nrow()-1; j >= 0; j--) {
 
378
    ib[j+1] = ib[j];
 
379
  }
 
380
  ib[0] = q0;
 
381
}
 
382
// ----------------------------------------------------------------------------
251
383
// instanciation in library
252
384
// ----------------------------------------------------------------------------
253
385
template class csr_seq_rep<Float>;