← Back to branch summary

~ubuntu-branches/ubuntu/saucy/rheolef/saucy

« back to all changes in this revision

Viewing changes to skit/plib2/solver_pastix_seq.cc

  • Committer: Bazaar Package Importer
  • Author(s): Pierre Saramito
  • Date: 2011-03-23 11:14:26 UTC
  • mfrom: (1.1.4 upstream)
  • Revision ID: james.westby@ubuntu.com-20110323111426-cjvhey7lxt6077ty
Tags: 5.93-1
* New upstream release (minor changes):
  - some extra warning message deleted in heap_allocator
  - graphic output with mayavi2 fixed
  - add doc refman in .info and .pdf format

Show diffs side-by-side

added added

removed removed

Lines of Context:
skit/plib2/solver_pastix_seq.cc
 
1
///
 
2
/// This file is part of Rheolef.
 
3
///
 
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
5
///
 
6
/// Rheolef is free software; you can redistribute it and/or modify
 
7
/// it under the terms of the GNU General Public License as published by
 
8
/// the Free Software Foundation; either version 2 of the License, or
 
9
/// (at your option) any later version.
 
10
///
 
11
/// Rheolef is sequential in the hope that it will be useful,
 
12
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
/// GNU General Public License for more details.
 
15
///
 
16
/// You should have received a copy of the GNU General Public License
 
17
/// along with Rheolef; if not, write to the Free Software
 
18
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
19
///
 
20
/// =========================================================================
 
21
// c++ & boost::mpi adaptation of pastix "simple.c" sequential example
 
22
//
 
23
#include "rheolef/config.h"
 
24
#ifdef _RHEOLEF_HAVE_PASTIX
 
25
#include "solver_pastix.h"
 
26
#include "rheolef/pcg.h"
 
27
#include "rheolef/eye.h"
 
28
 
 
29
namespace rheolef {
 
30
using namespace std;
 
31
 
 
32
template<class T, class M>
 
33
void
 
34
solver_pastix_base_rep<T,M>::resize (pastix_int_t n, pastix_int_t nnz)
 
35
{
 
36
   _n = n;
 
37
   _nnz  = nnz;
 
38
   _ptr.resize(n+1);
 
39
   _idx.resize(nnz);
 
40
   _val.resize(nnz);
 
41
} 
 
42
template<class T, class M>
 
43
void
 
44
solver_pastix_base_rep<T,M>::check () const
 
45
{
 
46
  warning_macro ("check...");
 
47
  if (_n == 0) return;
 
48
  /**
 
49
   * Pastix matrix needs :
 
50
   *    - to be in fortran numbering
 
51
   *    - to have only the lower triangular part in symmetric case
 
52
   *    - to have a graph with a symmetric structure in unsymmetric case
 
53
   */
 
54
  pastix_int_t  symmetry = (is_symmetric() ? API_SYM_YES : API_SYM_NO);
 
55
  pastix_int_t *ptr_begin = (pastix_int_t*) _ptr.begin().operator->();
 
56
  pastix_int_t *idx_begin = (pastix_int_t*) _idx.begin().operator->();
 
57
  T            *val_begin = (T*)            _val.begin().operator->(); 
 
58
  pastix_int_t status = d_pastix_checkMatrix (0, _opt.verbose_level, symmetry,  API_YES, 
 
59
               _n, &ptr_begin, &idx_begin, &val_begin, NULL, 1); 
 
60
  check_macro (status == 0, "pastix check returns error status = " << status);
 
61
  warning_macro ("pastix matrix check: ok");
 
62
}
 
63
template<class T, class M>
 
64
void
 
65
solver_pastix_base_rep<T,M>::load_both_continued (const csr<T,M>& a) 
 
66
{
 
67
  size_t first_dis_i = a.row_ownership().first_index();
 
68
  size_t first_dis_j = a.col_ownership().first_index();
 
69
  typename csr<T,M>::const_iterator aptr = a.begin();
 
70
  pastix_int_t bp = 0;
 
71
  _ptr [0] = bp + _base;
 
72
  for (size_t i = 0; i < a.nrow(); i++) {
 
73
    size_t dis_i = first_dis_i + i;
 
74
    for (typename csr<T,M>::const_data_iterator ap = aptr[i]; ap < aptr[i+1]; ap++) {
 
75
      size_t        j   = (*ap).first;
 
76
      const T& val = (*ap).second;
 
77
      size_t dis_j = first_dis_j + j;
 
78
      if (_is_sym && dis_i > dis_j) continue;
 
79
      _val      [bp] = val;
 
80
      _idx      [bp] = dis_j + _base;
 
81
      bp++;
 
82
    }
 
83
    _ptr [i+1] = bp + _base;
 
84
  }
 
85
  check_macro (bp == _nnz, "factorization: invalid nnz count");
 
86
}
 
87
template<class T, class M>
 
88
void
 
89
solver_pastix_base_rep<T,M>::load_unsymmetric (const csr<T,M>& a) 
 
90
{
 
91
  size_t n   = a.nrow();
 
92
  size_t nnz = a.nnz();
 
93
  resize (n, nnz);
 
94
  load_both_continued (a); 
 
95
}
 
96
template<class T, class M>
 
97
void
 
98
solver_pastix_base_rep<T,M>::load_symmetric (const csr<T,M>& a) 
 
99
{
 
100
  check_macro ((a.nnz()-a.ncol()) % 2 == 0, "factorization: symmetric csr has unsymmetric structure");
 
101
  // conserve only the lower part of the csc(pastix) = the upper past of the csr(rheolef)
 
102
  size_t n   = a.nrow();
 
103
  size_t nnz = a.nnz() - (a.nnz() - a.nrow())/2; // TODO: bug when matrix has some zero on the diag: fix it !
 
104
  resize (n, nnz);
 
105
  load_both_continued (a);
 
106
}
 
107
template<class T, class M>
 
108
void
 
109
solver_pastix_base_rep<T,M>::symbolic_factorization ()
 
110
{
 
111
  if (_n == 0) return;
 
112
  if (_have_pastix_bug_small_matrix) {
 
113
    warning_macro ("sym_fact: bug, circumvent !");
 
114
    return;
 
115
  }
 
116
  warning_macro ("sym_fact [1]...");
 
117
  const pastix_int_t nbthread = 1; // threads are not yet very well supported with openmpi & scotch
 
118
  const pastix_int_t ordering = 0; // scotch
 
119
 
 
120
  // tasks :
 
121
  //  0. set params to default values
 
122
  _iparm[IPARM_START_TASK]       = API_TASK_INIT;
 
123
  _iparm[IPARM_END_TASK]         = API_TASK_INIT;
 
124
  _iparm[IPARM_MODIFY_PARAMETER] = API_NO;
 
125
  d_pastix (&_pastix_data,
 
126
          0, 
 
127
          _n,
 
128
          _ptr.begin().operator->(),
 
129
          _idx.begin().operator->(),
 
130
          NULL, // _val.begin().operator->(),
 
131
          NULL,
 
132
          NULL,
 
133
          NULL,
 
134
          0,
 
135
          _iparm,
 
136
          _dparm);
 
137
 
 
138
  // Customize some parameters
 
139
  _iparm[IPARM_THREAD_NBR] = nbthread;
 
140
  if (is_symmetric()) {
 
141
      _iparm[IPARM_SYM]           = API_SYM_YES;
 
142
      _iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
 
143
  } else {
 
144
      _iparm[IPARM_SYM]           = API_SYM_NO;
 
145
      _iparm[IPARM_FACTORIZATION] = API_FACT_LU;
 
146
  }
 
147
  _iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
 
148
  _iparm[IPARM_VERBOSE]             = _opt.verbose_level;
 
149
  _iparm[IPARM_ORDERING]            = ordering;
 
150
  bool do_incomplete;
 
151
  if (_opt.incomplete != solver_option_type::decide) {
 
152
    do_incomplete = (_opt.incomplete != 0);
 
153
  } else {
 
154
    do_incomplete = (_pattern_dimension > 2); // 3d-pattern => iterative & IC(k) precond
 
155
  }
 
156
  warning_macro ("do_incomplete = " << do_incomplete);
 
157
  _iparm[IPARM_INCOMPLETE]          = (do_incomplete ? 1 : 0);
 
158
  _iparm[IPARM_OOC_LIMIT]           = _opt.ooc;
 
159
  if (_opt.incomplete == 1) {
 
160
      _dparm[DPARM_EPSILON_REFINEMENT] = _opt.tol;
 
161
  }
 
162
  _iparm[IPARM_LEVEL_OF_FILL]       = _opt.level_of_fill;
 
163
  _iparm[IPARM_AMALGAMATION_LEVEL]  = _opt.amalgamation;
 
164
  _iparm[IPARM_RHS_MAKING]          = API_RHS_B;
 
165
  // 1) get the i2new_dis_i array: its indexes are in the [0:dis_n[ range
 
166
  //     i2new_dis_i[i] = i2new_dis_i [i]
 
167
  _i2new_dis_i.resize (_n);
 
168
 
 
169
  warning_macro ("sym_fact [2] i2new_dis_i...");
 
170
  warning_macro ("sym_fact [2] n="<<_n);
 
171
  warning_macro ("sym_fact [2] _ptr.begin="<< _ptr.begin().operator->());
 
172
  warning_macro ("sym_fact [2] _idx.begin="<< _idx.begin().operator->());
 
173
  warning_macro ("sym_fact [2] _i2new_dis_i.begin="<< _i2new_dis_i.begin().operator->());
 
174
 
 
175
  pastix_int_t itmp = 0; // pastix does not manage NULL pointers: when n=0, vector::begin() retuns a NULL...
 
176
  T            vtmp = 0; // pastix does not manage NULL pointers: when n=0, vector::begin() retuns a NULL...
 
177
 
 
178
  // tasks :
 
179
  //  1. ordering
 
180
  //  2. symbolic factorization
 
181
  //  3. tasks mapping and scheduling
 
182
  _iparm[IPARM_START_TASK]          = API_TASK_ORDERING;
 
183
  _iparm[IPARM_END_TASK]            = API_TASK_ANALYSE;
 
184
 
 
185
  std::vector<pastix_int_t> new_i2i (_n);
 
186
  std::vector<double> dummy_rhs (_n);
 
187
  warning_macro ("sym_fact [3] call pastix... _n="<<_n);
 
188
  warning_macro ("sym_fact [3] call pastix... _ptr="<<_ptr.begin().operator->());
 
189
  warning_macro ("sym_fact [3] call pastix... _idx="<<_idx.begin().operator->());
 
190
  warning_macro ("sym_fact [3] call pastix... _i2new_i="<<_i2new_dis_i.begin().operator->());
 
191
  warning_macro ("sym_fact [3] call pastix... new_i2i="<<new_i2i.begin().operator->());
 
192
  d_pastix (&_pastix_data, 
 
193
          0, 
 
194
          _n, 
 
195
          _ptr.begin().operator->(),
 
196
          (_idx.begin().operator->() != NULL) ? _idx.begin().operator->() : &itmp,
 
197
          0,
 
198
          (_i2new_dis_i.begin().operator->()  != NULL) ? _i2new_dis_i.begin().operator->()  : &itmp,
 
199
          (new_i2i.begin().operator->()  != NULL) ? new_i2i.begin().operator->()  : &itmp,
 
200
          0,
 
201
          0,
 
202
          _iparm,
 
203
          _dparm);
 
204
  warning_macro ("sym_fact [3] call pastix done");
 
205
 
 
206
  warning_macro ("sym_fact done");
 
207
}
 
208
template<class T, class M>
 
209
void
 
210
solver_pastix_base_rep<T,M>::numeric_factorization ()
 
211
{
 
212
  if (_n == 0) return;
 
213
  if (_have_pastix_bug_small_matrix) {
 
214
    warning_macro ("num_fact: bug, circumvent !");
 
215
    return;
 
216
  }
 
217
  // pastix tasks:
 
218
  //    4. numerical factorization
 
219
  warning_macro ("num_fact [2] factorise...");
 
220
  pastix_int_t itmp = 0; // pastix does not manage NULL pointers: when n=0, vector::begin() retuns a NULL...
 
221
  T vtmp = 0; // pastix does not manage NULL rhs pointer: when new_n=0, but vector::begin() retuns a NULL...
 
222
  _iparm[IPARM_START_TASK]          = API_TASK_NUMFACT;
 
223
  _iparm[IPARM_END_TASK]            = API_TASK_NUMFACT;
 
224
  d_pastix (&_pastix_data, 
 
225
          0, 
 
226
          _n, 
 
227
          _ptr.begin().operator->(),
 
228
          (_idx.begin().operator->() != NULL) ? _idx.begin().operator->() : &itmp,
 
229
          (_val.begin().operator->() != NULL) ? _val.begin().operator->() : &vtmp,
 
230
          (_i2new_dis_i.begin().operator->()  != NULL) ? _i2new_dis_i.begin().operator->()  : &itmp,
 
231
          NULL,
 
232
          NULL,
 
233
          0,
 
234
          _iparm,
 
235
          _dparm);
 
236
 
 
237
  warning_macro ("num fact done");
 
238
}
 
239
template<class T, class M>
 
240
void
 
241
solver_pastix_base_rep<T,M>::load (const csr<T,M>& a, const solver_option_type& opt)
 
242
{
 
243
   _is_sym = a.is_symmetric(); 
 
244
   _pattern_dimension = a.pattern_dimension(); 
 
245
   warning_macro ("is_sym = " << _is_sym);
 
246
   warning_macro ("pattern_dim = " << _pattern_dimension);
 
247
   _csr_row_ownership = a.row_ownership();
 
248
   _csr_col_ownership = a.col_ownership();
 
249
   _opt = opt;
 
250
 
 
251
  check_macro (a.nrow() == a.ncol(), "factorization: only square matrix are supported");
 
252
 
 
253
  if (_is_sym) {
 
254
    load_symmetric(a);
 
255
  } else {
 
256
    load_unsymmetric(a);
 
257
  }
 
258
  if (_opt.do_check) {
 
259
    check ();
 
260
  }
 
261
  symbolic_factorization ();
 
262
  numeric_factorization();
 
263
}
 
264
template<class T, class M>
 
265
solver_pastix_base_rep<T,M>::solver_pastix_base_rep ()
 
266
 : _n(),
 
267
   _nnz(),
 
268
   _ptr(),
 
269
   _idx(),
 
270
   _val(),
 
271
   _is_sym(false),
 
272
   _pattern_dimension(0),
 
273
   _pastix_data(0),
 
274
   _iparm(),
 
275
   _dparm(),
 
276
   _csr_row_ownership(),
 
277
   _csr_col_ownership(),
 
278
   _opt(),
 
279
   _new_rhs(),
 
280
   _new_i2dis_i_base(),
 
281
   _i2new_dis_i(),
 
282
   _have_pastix_bug_small_matrix(false),
 
283
   _a_when_bug()
 
284
{
 
285
}
 
286
template<class T, class M>
 
287
solver_pastix_base_rep<T,M>::solver_pastix_base_rep (const csr<T,M>& a, const solver_option_type& opt)
 
288
 : _n(),
 
289
   _nnz(),
 
290
   _ptr(),
 
291
   _idx(),
 
292
   _val(),
 
293
   _is_sym(false),
 
294
   _pattern_dimension(0),
 
295
   _pastix_data(0),
 
296
   _iparm(),
 
297
   _dparm(),
 
298
   _csr_row_ownership(),
 
299
   _csr_col_ownership(),
 
300
   _opt(),
 
301
   _new_rhs(),
 
302
   _new_i2dis_i_base(),
 
303
   _i2new_dis_i(),
 
304
   _have_pastix_bug_small_matrix(false),
 
305
   _a_when_bug()
 
306
{
 
307
  load (a, opt);
 
308
}
 
309
template<class T, class M>
 
310
void
 
311
solver_pastix_base_rep<T,M>::update_values (const csr<T,M>& a) 
 
312
{
 
313
  if (_n == 0) return;
 
314
  check_macro (size_t(_n) == a.nrow() && size_t(_n) == a.ncol(),
 
315
    "update: local input matrix size distribution mismatch: ("<<a.nrow()<<","<<a.ncol()<<"), expect ("
 
316
    << _n << "," << _n << ")");
 
317
  size_t nnz_a;
 
318
  if (!_is_sym) {
 
319
    nnz_a = a.nnz();
 
320
  } else {
 
321
    // conserve only the lower part of the csc(pastix) = the upper past of the csr(rheolef)
 
322
    nnz_a = a.nnz() - (a.nnz() - a.nrow())/2; // TODO: fix when zero on diag
 
323
  }
 
324
  check_macro (size_t(_nnz) == nnz_a,
 
325
    "update: local input matrix nnz distribution mismatch: nnz="<<nnz_a<<", expect nnz="<<_nnz);
 
326
 
 
327
  size_t first_dis_i = a.row_ownership().first_index();
 
328
  size_t first_dis_j = a.col_ownership().first_index();
 
329
  pastix_int_t bp = 0;
 
330
  typename csr<T,M>::const_iterator aptr = a.begin();
 
331
  for (size_t i = 0; i < a.nrow(); i++) {
 
332
    size_t dis_i = first_dis_i + i;
 
333
    for (typename csr<T,M>::const_data_iterator ap = aptr[i]; ap < aptr[i+1]; ap++) {
 
334
      size_t     j = (*ap).first;
 
335
      size_t dis_j = first_dis_j + j;
 
336
      if (_is_sym && dis_i > dis_j) continue;
 
337
      _val [bp] = (*ap).second;
 
338
      bp++;
 
339
    }
 
340
  }
 
341
  numeric_factorization();
 
342
}
 
343
template<class T, class M>
 
344
vec<T,M>
 
345
solver_pastix_base_rep<T,M>::trans_solve (const vec<T,M>& rhs)
 
346
{
 
347
  if (_n == 0) return rhs;
 
348
  // ================================================================
 
349
  // solve
 
350
  // ================================================================
 
351
  warning_macro ("trans_solve [1]... n="<<_n);
 
352
  check_macro (rhs.size() == size_t(_n), "invalid rhs size="<<rhs.size()<<": expect size="<<_n);
 
353
 
 
354
  _new_rhs.resize (_n);
 
355
  for (pastix_int_t i = 0; i < _n; i++) {
 
356
    _new_rhs [i] = rhs [i];
 
357
  }
 
358
  // tasks:
 
359
  //    5. numerical solve
 
360
  //    6. numerical refinement
 
361
  //    7. clean
 
362
  warning_macro ("trans_solve [3] solve...new_rhs.ptr="<< _new_rhs.begin().operator->());
 
363
  pastix_int_t itmp = 0; // pastix does not manage NULL pointers: when n=0, vector::begin() retuns a NULL...
 
364
  T vtmp = 0; // pastix does not manage NULL rhs pointer: when new_n=0, but vector::begin() retuns a NULL...
 
365
  _iparm[IPARM_START_TASK]          = API_TASK_SOLVE;
 
366
  _iparm[IPARM_END_TASK]            = API_TASK_REFINE;
 
367
  d_pastix (&_pastix_data, 
 
368
          0, 
 
369
          _n, 
 
370
          _ptr.begin().operator->(),
 
371
          (_idx.begin().operator->() != NULL) ? _idx.begin().operator->() : &itmp,
 
372
          (_val.begin().operator->() != NULL) ? _val.begin().operator->() : &vtmp,
 
373
          (_i2new_dis_i.begin().operator->()  != NULL) ? _i2new_dis_i.begin().operator->()  : &itmp,
 
374
          NULL,
 
375
          (_n != 0) ? _new_rhs.begin().operator->() : &vtmp,
 
376
          1,
 
377
          _iparm,
 
378
          _dparm);
 
379
 
 
380
  // new_i2dis_i_base [new_i] - base = new_i2dis_i [new_i]
 
381
  warning_macro ("trans_solve [4] perm(x)...");
 
382
  vec<T,M> x (_csr_row_ownership);
 
383
  for (pastix_int_t i = 0; i < _n; i++) {
 
384
    x [i] = _new_rhs [i];
 
385
  }
 
386
  warning_macro ("trans_solve done");
 
387
  return x;
 
388
}
 
389
template<class T, class M>
 
390
vec<T,M>
 
391
solver_pastix_base_rep<T,M>::solve (const vec<T,M>& rhs)
 
392
{
 
393
  warning_macro ("solve...");
 
394
  // TODO: make a csc<T> wrapper arround csr<T> and use csc<T> in form & form_assembly
 
395
  // => avoid the transposition
 
396
  if (! _is_sym) {
 
397
    warning_macro ("solve: not yet supported for unsymmetric matrix; continue with trans_solve....");
 
398
  }
 
399
  return trans_solve (rhs);
 
400
}
 
401
template<class T, class M>
 
402
solver_pastix_base_rep<T,M>::~solver_pastix_base_rep()
 
403
{
 
404
  warning_macro ("dstor...");
 
405
  if (_pastix_data == 0) return;
 
406
  // tasks:
 
407
  //    7. clean
 
408
  _iparm[IPARM_START_TASK]          = API_TASK_CLEAN;
 
409
  _iparm[IPARM_END_TASK]            = API_TASK_CLEAN;
 
410
  d_pastix (&_pastix_data, 
 
411
          0, 
 
412
          _n, 
 
413
          0,
 
414
          0, 
 
415
          0, 
 
416
          NULL, 
 
417
          NULL,
 
418
          _new_rhs.begin().operator->(),
 
419
          1,
 
420
          _iparm,
 
421
          _dparm);
 
422
 
 
423
  // was allocated by d_cscd_redispatch()
 
424
  warning_macro ("dstor call free...");
 
425
  _pastix_data == 0;
 
426
  warning_macro ("dstor done");
 
427
}
 
428
// ----------------------------------------------------------------------------
 
429
// instanciation in library
 
430
// ----------------------------------------------------------------------------
 
431
// TODO: code is only valid here for T=double (d_pastix etc)
 
432
 
 
433
template class solver_pastix_base_rep<double,sequential>;
 
434
 
 
435
#ifdef _RHEOLEF_HAVE_MPI
 
436
template class solver_pastix_base_rep<double,distributed>;
 
437
#endif // _RHEOLEF_HAVE_MPI
 
438
 
 
439
} // namespace rheolef
 
440
#endif // _RHEOLEF_HAVE_PASTIX