1 /*
2  *
3  * Copyright (c) Kresimir Fresl 2003
4  *
5  * Distributed under the Boost Software License, Version 1.0.
6  * (See accompanying file LICENSE_1_0.txt or copy at
7  * http://www.boost.org/LICENSE_1_0.txt)
8  *
9  * Author acknowledges the support of the Faculty of Civil Engineering,
10  * University of Zagreb, Croatia.
11  *
12  */
13 
14 /* for UMFPACK Copyright, License and Availability see umfpack_inc.hpp */
15 
16 
17 #ifndef BOOST_NUMERIC_BINDINGS_UMFPACK_HPP
18 #define BOOST_NUMERIC_BINDINGS_UMFPACK_HPP
19 
20 
21 #include <boost/noncopyable.hpp>
22 #include <boost/numeric/bindings/umfpack/umfpack_overloads.hpp>
23 #include <boost/numeric/bindings/value_type.hpp>
24 #include <boost/numeric/bindings/begin.hpp>
25 #include <boost/numeric/bindings/end.hpp>
26 #include <boost/numeric/bindings/size.hpp>
27 #include <boost/numeric/bindings/data_order.hpp>
28 #include <boost/numeric/bindings/index_base.hpp>
29 
30 
31 namespace boost { namespace numeric { namespace bindings {  namespace umfpack {
32 
33   template <typename MatrA>
check_umfpack_structure()34   void check_umfpack_structure()
35   {
36     BOOST_STATIC_ASSERT((boost::is_same<
37       typename bindings::detail::property_at< MatrA, tag::matrix_type >::type,
38       tag::general
39     >::value));
40     BOOST_STATIC_ASSERT((boost::is_same<
41       typename bindings::result_of::data_order<MatrA>::type,
42       tag::column_major
43     >::value));
44     typedef typename bindings::result_of::index_base<MatrA>::type index_b;
45     BOOST_STATIC_ASSERT(index_b::value == 0);
46     typedef typename bindings::detail::property_at<
47       MatrA, tag::data_structure >::type storage_f;
48     BOOST_STATIC_ASSERT(
49       (boost::is_same<storage_f, tag::compressed_sparse>::value ||
50        boost::is_same<storage_f, tag::coordinate_sparse>::value ));
51   }
52 
53   template <typename T = double>
54   struct symbolic_type : private noncopyable {
55     void *ptr;
symbolic_typeboost::numeric::bindings::umfpack::symbolic_type56     symbolic_type():ptr(0){}
~symbolic_typeboost::numeric::bindings::umfpack::symbolic_type57     ~symbolic_type() {
58       if (ptr)
59         detail::free_symbolic (T(), 0, &ptr);
60     }
freeboost::numeric::bindings::umfpack::symbolic_type61     void free() {
62       if (ptr)
63         detail::free_symbolic (T(), 0, &ptr);
64       ptr = 0;
65     }
66   };
67 
68   template <typename T>
free(symbolic_type<T> & s)69   void free (symbolic_type<T>& s) { s.free(); }
70 
71   template <typename T = double>
72   struct numeric_type : private noncopyable {
73     void *ptr;
numeric_typeboost::numeric::bindings::umfpack::numeric_type74     numeric_type():ptr(0){}
~numeric_typeboost::numeric::bindings::umfpack::numeric_type75     ~numeric_type() {
76       if (ptr)
77         detail::free_numeric (T(), 0, &ptr);
78     }
freeboost::numeric::bindings::umfpack::numeric_type79     void free() {
80       if (ptr)
81         detail::free_numeric (T(), 0, &ptr);
82       ptr = 0;
83     }
84   };
85 
86   template <typename T>
free(numeric_type<T> & n)87   void free (numeric_type<T>& n) { n.free(); }
88 
89 
90   template <typename T = double>
91   struct control_type : private noncopyable {
92     double ptr[UMFPACK_CONTROL];
control_typeboost::numeric::bindings::umfpack::control_type93     control_type() { detail::defaults (T(), 0, ptr); }
operator []boost::numeric::bindings::umfpack::control_type94     double operator[] (int i) const { return ptr[i]; }
operator []boost::numeric::bindings::umfpack::control_type95     double& operator[] (int i) { return ptr[i]; }
defaultsboost::numeric::bindings::umfpack::control_type96     void defaults() { detail::defaults (T(), 0, ptr); }
97   };
98 
99   template <typename T>
defaults(control_type<T> & c)100   void defaults (control_type<T>& c) { c.defaults(); }
101 
102   template <typename T = double>
103   struct info_type : private noncopyable {
104     double ptr[UMFPACK_INFO];
operator []boost::numeric::bindings::umfpack::info_type105     double operator[] (int i) const { return ptr[i]; }
operator []boost::numeric::bindings::umfpack::info_type106     double& operator[] (int i) { return ptr[i]; }
107   };
108 
109 
110   /////////////////////////////////////
111   // solving system of linear equations
112   /////////////////////////////////////
113 
114 
115   // symbolic
116   /*
117    * Given nonzero pattern of a sparse matrix A in column-oriented form,
118    * umfpack_*_symbolic performs a column pre-ordering to reduce fill-in
119    * (using COLAMD or AMD) and a symbolic factorisation.  This is required
120    * before the matrix can be numerically factorised with umfpack_*_numeric.
121    */
122   namespace detail {
123 
124     template <typename MatrA>
125     inline
symbolic(tag::compressed_sparse,MatrA const & A,void ** Symbolic,double const * Control=0,double * Info=0)126     int symbolic (tag::compressed_sparse,
127                   MatrA const& A, void **Symbolic,
128                   double const* Control = 0, double* Info = 0)
129     {
130       return detail::symbolic (bindings::size_row (A),
131                                bindings::size_column (A),
132                                bindings::begin_compressed_index_major (A),
133                                bindings::begin_index_minor (A),
134                                bindings::begin_value (A),
135                                Symbolic, Control, Info);
136     }
137 
138     template <typename MatrA, typename QVec>
139     inline
symbolic(tag::compressed_sparse,MatrA const & A,QVec const & Qinit,void ** Symbolic,double const * Control=0,double * Info=0)140     int symbolic (tag::compressed_sparse,
141                   MatrA const& A, QVec const& Qinit, void **Symbolic,
142                   double const* Control = 0, double* Info = 0)
143     {
144 #ifdef CHECK_TEST_COVERAGE
145       typedef typename MatrA::not_yet_tested i_m_still_here;
146 #endif
147       return detail::qsymbolic (bindings::size_row (A),
148                                 bindings::size_column (A),
149                                 bindings::begin_compressed_index_major (A),
150                                 bindings::begin_index_minor (A),
151                                 bindings::begin_value (A),
152                                 bindings::begin_value (Qinit),
153                                 Symbolic, Control, Info);
154     }
155 
156     template <typename MatrA>
157     inline
symbolic(tag::coordinate_sparse,MatrA const & A,void ** Symbolic,double const * Control=0,double * Info=0)158     int symbolic (tag::coordinate_sparse,
159                   MatrA const& A, void **Symbolic,
160                   double const* Control = 0, double* Info = 0)
161     {
162       int n_row = bindings::size_row (A);
163       int n_col = bindings::size_column (A);
164       int nnz = bindings::end_value (A) - bindings::begin_value (A);
165 
166       typedef typename bindings::value_type<MatrA>::type val_t;
167 
168       int const* Ti = bindings::begin_index_minor (A);
169       int const* Tj = bindings::begin_index_major (A);
170       bindings::detail::array<int> Ap (n_col+1);
171       if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory;
172       bindings::detail::array<int> Ai (nnz);
173       if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory;
174 
175       int status = detail::triplet_to_col (n_row, n_col, nnz,
176                                            Ti, Tj, static_cast<val_t*> (0),
177                                            Ap.storage(), Ai.storage(),
178                                            static_cast<val_t*> (0), 0);
179       if (status != UMFPACK_OK) return status;
180 
181       return detail::symbolic (n_row, n_col,
182                                Ap.storage(), Ai.storage(),
183                                bindings::begin_value (A),
184                                Symbolic, Control, Info);
185     }
186 
187     template <typename MatrA, typename QVec>
188     inline
symbolic(tag::coordinate_sparse,MatrA const & A,QVec const & Qinit,void ** Symbolic,double const * Control=0,double * Info=0)189     int symbolic (tag::coordinate_sparse,
190                   MatrA const& A, QVec const& Qinit, void **Symbolic,
191                   double const* Control = 0, double* Info = 0)
192     {
193 #ifdef CHECK_TEST_COVERAGE
194       typedef typename MatrA::not_yet_tested i_m_still_here;
195 #endif
196       int n_row = bindings::size_row (A);
197       int n_col = bindings::size_column (A);
198       int nnz = bindings::end_value (A) - bindings::begin_value (A);
199 
200       typedef typename bindings::value_type<MatrA>::type val_t;
201 
202       int const* Ti = bindings::begin_index_minor (A);
203       int const* Tj = bindings::begin_index_major (A);
204       bindings::detail::array<int> Ap (n_col+1);
205       if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory;
206       bindings::detail::array<int> Ai (nnz);
207       if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory;
208 
209       int status = detail::triplet_to_col (n_row, n_col, nnz,
210                                            Ti, Tj, static_cast<val_t*> (0),
211                                            Ap.storage(), Ai.storage(),
212                                            static_cast<val_t*> (0), 0);
213       if (status != UMFPACK_OK) return status;
214 
215       return detail::qsymbolic (n_row, n_col,
216                                 Ap.storage(), Ai.storage(),
217                                 bindings::begin_value (A),
218                                 bindings::begin_value (Qinit),
219                                 Symbolic, Control, Info);
220     }
221 
222   } // detail
223 
224   template <typename MatrA>
225   inline
symbolic(MatrA const & A,symbolic_type<typename bindings::value_type<MatrA>::type> & Symbolic,double const * Control=0,double * Info=0)226   int symbolic (MatrA const& A,
227                 symbolic_type<
228                   typename bindings::value_type<MatrA>::type
229                 >& Symbolic,
230                 double const* Control = 0, double* Info = 0)
231   {
232 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
233     check_umfpack_structure<MatrA>();
234 #endif
235     typedef typename bindings::detail::property_at<
236       MatrA, tag::data_structure >::type storage_f;
237 
238     return detail::symbolic (storage_f(), A, &Symbolic.ptr, Control, Info);
239   }
240 
241   template <typename MatrA>
242   inline
symbolic(MatrA const & A,symbolic_type<typename bindings::value_type<MatrA>::type> & Symbolic,control_type<typename bindings::value_type<MatrA>::type> const & Control,info_type<typename bindings::value_type<MatrA>::type> & Info)243   int symbolic (MatrA const& A,
244                 symbolic_type<
245                   typename bindings::value_type<MatrA>::type
246                 >& Symbolic,
247                 control_type<
248                   typename bindings::value_type<MatrA>::type
249                 > const& Control,
250                 info_type<
251                   typename bindings::value_type<MatrA>::type
252                 >& Info)
253   {
254     return symbolic (A, Symbolic, Control.ptr, Info.ptr);
255   }
256 
257   template <typename MatrA>
258   inline
symbolic(MatrA const & A,symbolic_type<typename bindings::value_type<MatrA>::type> & Symbolic,control_type<typename bindings::value_type<MatrA>::type> const & Control)259   int symbolic (MatrA const& A,
260                 symbolic_type<
261                   typename bindings::value_type<MatrA>::type
262                 >& Symbolic,
263                 control_type<
264                   typename bindings::value_type<MatrA>::type
265                 > const& Control)
266   {
267     return symbolic (A, Symbolic, Control.ptr);
268   }
269 
270   template <typename MatrA, typename QVec>
271   inline
symbolic(MatrA const & A,QVec const & Qinit,symbolic_type<typename bindings::value_type<MatrA>::type> & Symbolic,double const * Control=0,double * Info=0)272   int symbolic (MatrA const& A, QVec const& Qinit,
273                 symbolic_type<
274                   typename bindings::value_type<MatrA>::type
275                 >& Symbolic,
276                 double const* Control = 0, double* Info = 0)
277   {
278 #ifdef CHECK_TEST_COVERAGE
279       typedef typename MatrA::not_yet_tested i_m_still_here;
280 #endif
281 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
282     check_umfpack_structure<MatrA>();
283 #endif
284     typedef typename bindings::detail::property_at<
285       MatrA, tag::data_structure >::type storage_f;
286 
287     assert (bindings::size_column (A) == bindings::size (Qinit));
288 
289     return detail::symbolic (storage_f(), A, Qinit,
290                              &Symbolic.ptr, Control, Info);
291   }
292 
293   template <typename MatrA, typename QVec>
294   inline
symbolic(MatrA const & A,QVec const & Qinit,symbolic_type<typename bindings::value_type<MatrA>::type> & Symbolic,control_type<typename bindings::value_type<MatrA>::type> const & Control,info_type<typename bindings::value_type<MatrA>::type> & Info)295   int symbolic (MatrA const& A, QVec const& Qinit,
296                 symbolic_type<
297                   typename bindings::value_type<MatrA>::type
298                 >& Symbolic,
299                 control_type<
300                   typename bindings::value_type<MatrA>::type
301                 > const& Control,
302                 info_type<
303                   typename bindings::value_type<MatrA>::type
304                 >& Info)
305   {
306     return symbolic (A, Qinit, Symbolic, Control.ptr, Info.ptr);
307   }
308 
309   template <typename MatrA, typename QVec>
310   inline
symbolic(MatrA const & A,QVec const & Qinit,symbolic_type<typename bindings::value_type<MatrA>::type> & Symbolic,control_type<typename bindings::value_type<MatrA>::type> const & Control)311   int symbolic (MatrA const& A, QVec const& Qinit,
312                 symbolic_type<
313                   typename bindings::value_type<MatrA>::type
314                 >& Symbolic,
315                 control_type<
316                   typename bindings::value_type<MatrA>::type
317                 > const& Control)
318   {
319     return symbolic (A, Qinit, Symbolic, Control.ptr);
320   }
321 
322 
323   // numeric
324   /*
325    * Given a sparse matrix A in column-oriented form, and a symbolic analysis
326    * computed by umfpack_*_*symbolic, the umfpack_*_numeric routine performs
327    * the numerical factorisation, PAQ=LU, PRAQ=LU, or P(R\A)Q=LU, where P
328    * and Q are permutation matrices (represented as permutation vectors),
329    * R is the row scaling, L is unit-lower triangular, and U is upper
330    * triangular.  This is required before the system Ax=b (or other related
331    * linear systems) can be solved.
332    */
333   namespace detail {
334 
335     template <typename MatrA>
336     inline
numeric(tag::compressed_sparse,MatrA const & A,void * Symbolic,void ** Numeric,double const * Control=0,double * Info=0)337     int numeric (tag::compressed_sparse, MatrA const& A,
338                  void *Symbolic, void** Numeric,
339                  double const* Control = 0, double* Info = 0)
340     {
341       return detail::numeric (bindings::size_row (A),
342                               bindings::size_column (A),
343                               bindings::begin_compressed_index_major (A),
344                               bindings::begin_index_minor (A),
345                               bindings::begin_value (A),
346                               Symbolic, Numeric, Control, Info);
347     }
348 
349     template <typename MatrA>
350     inline
numeric(tag::coordinate_sparse,MatrA const & A,void * Symbolic,void ** Numeric,double const * Control=0,double * Info=0)351     int numeric (tag::coordinate_sparse, MatrA const& A,
352                  void *Symbolic, void** Numeric,
353                  double const* Control = 0, double* Info = 0)
354     {
355       int n_row = bindings::size_row (A);
356       int n_col = bindings::size_column (A);
357       int nnz = bindings::end_value (A) - bindings::begin_value (A);
358 
359       typedef typename bindings::value_type<MatrA>::type val_t;
360 
361       int const* Ti = bindings::begin_index_minor (A);
362       int const* Tj = bindings::begin_index_major (A);
363       bindings::detail::array<int> Ap (n_col+1);
364       if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory;
365       bindings::detail::array<int> Ai (nnz);
366       if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory;
367 
368       int status = detail::triplet_to_col (n_row, n_col, nnz,
369                                            Ti, Tj, static_cast<val_t*> (0),
370                                            Ap.storage(), Ai.storage(),
371                                            static_cast<val_t*> (0), 0);
372       if (status != UMFPACK_OK) return status;
373 
374       return detail::numeric (n_row, n_col,
375                               Ap.storage(), Ai.storage(),
376                               bindings::begin_value (A),
377                               Symbolic, Numeric, Control, Info);
378     }
379 
380   } // detail
381 
382   template <typename MatrA>
383   inline
numeric(MatrA const & A,symbolic_type<typename bindings::value_type<MatrA>::type> const & Symbolic,numeric_type<typename bindings::value_type<MatrA>::type> & Numeric,double const * Control=0,double * Info=0)384   int numeric (MatrA const& A,
385                symbolic_type<
386                  typename bindings::value_type<MatrA>::type
387                > const& Symbolic,
388                numeric_type<
389                  typename bindings::value_type<MatrA>::type
390                >& Numeric,
391                double const* Control = 0, double* Info = 0)
392   {
393 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
394     check_umfpack_structure<MatrA>();
395 #endif
396     typedef typename bindings::detail::property_at<
397       MatrA, tag::data_structure >::type storage_f;
398 
399     return detail::numeric (storage_f(), A,
400                             Symbolic.ptr, &Numeric.ptr, Control, Info);
401   }
402 
403   template <typename MatrA>
404   inline
numeric(MatrA const & A,symbolic_type<typename bindings::value_type<MatrA>::type> const & Symbolic,numeric_type<typename bindings::value_type<MatrA>::type> & Numeric,control_type<typename bindings::value_type<MatrA>::type> const & Control,info_type<typename bindings::value_type<MatrA>::type> & Info)405   int numeric (MatrA const& A,
406                symbolic_type<
407                  typename bindings::value_type<MatrA>::type
408                > const& Symbolic,
409                numeric_type<
410                  typename bindings::value_type<MatrA>::type
411                >& Numeric,
412                control_type<
413                  typename bindings::value_type<MatrA>::type
414                > const& Control,
415                info_type<
416                  typename bindings::value_type<MatrA>::type
417                >& Info)
418 
419   {
420     // g++ (3.2) is unable to distinguish
421     //           function numeric() and namespace boost::numeric ;o)
422     return umfpack::numeric (A, Symbolic, Numeric, Control.ptr, Info.ptr);
423   }
424 
425   template <typename MatrA>
426   inline
numeric(MatrA const & A,symbolic_type<typename bindings::value_type<MatrA>::type> const & Symbolic,numeric_type<typename bindings::value_type<MatrA>::type> & Numeric,control_type<typename bindings::value_type<MatrA>::type> const & Control)427   int numeric (MatrA const& A,
428                symbolic_type<
429                  typename bindings::value_type<MatrA>::type
430                > const& Symbolic,
431                numeric_type<
432                  typename bindings::value_type<MatrA>::type
433                >& Numeric,
434                control_type<
435                  typename bindings::value_type<MatrA>::type
436                > const& Control)
437   {
438     return umfpack::numeric (A, Symbolic, Numeric, Control.ptr);
439   }
440 
441 
442   // factor
443   /*
444    * symbolic and numeric
445    */
446   namespace detail {
447 
448     template <typename MatrA>
449     inline
factor(tag::compressed_sparse,MatrA const & A,void ** Numeric,double const * Control=0,double * Info=0)450     int factor (tag::compressed_sparse, MatrA const& A,
451                 void** Numeric, double const* Control = 0, double* Info = 0)
452     {
453 #ifdef CHECK_TEST_COVERAGE
454       typedef typename MatrA::not_yet_tested i_m_still_here;
455 #endif
456       symbolic_type<typename bindings::value_type<MatrA>::type>
457         Symbolic;
458 
459       int status;
460       status = detail::symbolic (bindings::size_row (A),
461                                  bindings::size_column (A),
462                                  bindings::begin_compressed_index_major (A),
463                                  bindings::begin_index_minor (A),
464                                  bindings::begin_value (A),
465                                  &Symbolic.ptr, Control, Info);
466       if (status != UMFPACK_OK) return status;
467 
468       return detail::numeric (bindings::size_row (A),
469                               bindings::size_column (A),
470                               bindings::begin_compressed_index_major (A),
471                               bindings::begin_index_minor (A),
472                               bindings::begin_value (A),
473                               Symbolic.ptr, Numeric, Control, Info);
474     }
475 
476     template <typename MatrA>
477     inline
factor(tag::coordinate_sparse,MatrA const & A,void ** Numeric,double const * Control=0,double * Info=0)478     int factor (tag::coordinate_sparse, MatrA const& A,
479                 void** Numeric, double const* Control = 0, double* Info = 0)
480     {
481 #ifdef CHECK_TEST_COVERAGE
482       typedef typename MatrA::not_yet_tested i_m_still_here;
483 #endif
484       int n_row = bindings::size_row (A);
485       int n_col = bindings::size_column (A);
486       int nnz = bindings::end_value (A) - bindings::begin_value (A);
487 
488       typedef typename bindings::value_type<MatrA>::type val_t;
489 
490       int const* Ti = bindings::begin_index_minor (A);
491       int const* Tj = bindings::begin_index_major (A);
492       bindings::detail::array<int> Ap (n_col+1);
493       if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory;
494       bindings::detail::array<int> Ai (nnz);
495       if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory;
496 
497       int status = detail::triplet_to_col (n_row, n_col, nnz,
498                                            Ti, Tj, static_cast<val_t*> (0),
499                                            Ap.storage(), Ai.storage(),
500                                            static_cast<val_t*> (0), 0);
501       if (status != UMFPACK_OK) return status;
502 
503       symbolic_type<typename bindings::value_type<MatrA>::type>
504         Symbolic;
505 
506       status = detail::symbolic (n_row, n_col,
507                                  Ap.storage(), Ai.storage(),
508                                  bindings::begin_value (A),
509                                  &Symbolic.ptr, Control, Info);
510       if (status != UMFPACK_OK) return status;
511 
512       return detail::numeric (n_row, n_col,
513                               Ap.storage(), Ai.storage(),
514                               bindings::begin_value (A),
515                               Symbolic.ptr, Numeric, Control, Info);
516     }
517 
518   } // detail
519 
520   template <typename MatrA>
521   inline
factor(MatrA const & A,numeric_type<typename bindings::value_type<MatrA>::type> & Numeric,double const * Control=0,double * Info=0)522   int factor (MatrA const& A,
523               numeric_type<
524                 typename bindings::value_type<MatrA>::type
525               >& Numeric,
526               double const* Control = 0, double* Info = 0)
527   {
528 #ifdef CHECK_TEST_COVERAGE
529       typedef typename MatrA::not_yet_tested i_m_still_here;
530 #endif
531 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
532     check_umfpack_structure<MatrA>();
533 #endif
534     typedef typename bindings::detail::property_at<
535       MatrA, tag::data_structure >::type storage_f;
536 
537     return detail::factor (storage_f(), A, &Numeric.ptr, Control, Info);
538   }
539 
540   template <typename MatrA>
541   inline
factor(MatrA const & A,numeric_type<typename bindings::value_type<MatrA>::type> & Numeric,control_type<typename bindings::value_type<MatrA>::type> const & Control,info_type<typename bindings::value_type<MatrA>::type> & Info)542   int factor (MatrA const& A,
543               numeric_type<
544                 typename bindings::value_type<MatrA>::type
545               >& Numeric,
546               control_type<
547                 typename bindings::value_type<MatrA>::type
548               > const& Control,
549               info_type<
550                 typename bindings::value_type<MatrA>::type
551               >& Info)
552   {
553     return factor (A, Numeric, Control.ptr, Info.ptr);
554   }
555 
556   template <typename MatrA>
557   inline
factor(MatrA const & A,numeric_type<typename bindings::value_type<MatrA>::type> & Numeric,control_type<typename bindings::value_type<MatrA>::type> const & Control)558   int factor (MatrA const& A,
559               numeric_type<
560                 typename bindings::value_type<MatrA>::type
561               >& Numeric,
562               control_type<
563                 typename bindings::value_type<MatrA>::type
564               > const& Control)
565   {
566     return factor (A, Numeric, Control.ptr);
567   }
568 
569 
570   // solve
571   /*
572    * Given LU factors computed by umfpack_*_numeric and the right-hand-side,
573    * B, solve a linear system for the solution X.  Iterative refinement is
574    * optionally performed.  Only square systems are handled.
575    */
576   namespace detail {
577 
578     template <typename MatrA, typename VecX, typename VecB>
579     inline
solve(tag::compressed_sparse,int sys,MatrA const & A,VecX & X,VecB const & B,void * Numeric,double const * Control=0,double * Info=0)580     int solve (tag::compressed_sparse, int sys,
581                MatrA const& A, VecX& X, VecB const& B,
582                void *Numeric, double const* Control = 0, double* Info = 0)
583     {
584       return detail::solve (sys, bindings::size_row (A),
585                             bindings::begin_compressed_index_major (A),
586                             bindings::begin_index_minor (A),
587                             bindings::begin_value (A),
588                             bindings::begin_value (X),
589                             bindings::begin_value (B),
590                             Numeric, Control, Info);
591     }
592 
593     template <typename MatrA, typename VecX, typename VecB>
594     inline
solve(tag::coordinate_sparse,int sys,MatrA const & A,VecX & X,VecB const & B,void * Numeric,double const * Control=0,double * Info=0)595     int solve (tag::coordinate_sparse, int sys,
596                MatrA const& A, VecX& X, VecB const& B,
597                void *Numeric, double const* Control = 0, double* Info = 0)
598     {
599 
600       int n = bindings::size_row (A);
601       int nnz = bindings::end_value (A) - bindings::begin_value (A);
602 
603       typedef typename bindings::value_type<MatrA>::type val_t;
604 
605       int const* Ti = bindings::begin_index_minor (A);
606       int const* Tj = bindings::begin_index_major (A);
607       bindings::detail::array<int> Ap (n+1);
608       if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory;
609       bindings::detail::array<int> Ai (nnz);
610       if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory;
611 
612       int status = detail::triplet_to_col (n, n, nnz,
613                                            Ti, Tj, static_cast<val_t*> (0),
614                                            Ap.storage(), Ai.storage(),
615                                            static_cast<val_t*> (0), 0);
616       if (status != UMFPACK_OK) return status;
617 
618       return detail::solve (sys, n, Ap.storage(), Ai.storage(),
619                             bindings::begin_value (A),
620                             bindings::begin_value (X),
621                             bindings::begin_value (B),
622                             Numeric, Control, Info);
623     }
624 
625   } // detail
626 
627   template <typename MatrA, typename VecX, typename VecB>
628   inline
solve(int sys,MatrA const & A,VecX & X,VecB const & B,numeric_type<typename bindings::value_type<MatrA>::type> const & Numeric,double const * Control=0,double * Info=0)629   int solve (int sys, MatrA const& A, VecX& X, VecB const& B,
630              numeric_type<
631                typename bindings::value_type<MatrA>::type
632              > const& Numeric,
633              double const* Control = 0, double* Info = 0)
634   {
635 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
636     check_umfpack_structure<MatrA>();
637 #endif
638     typedef typename bindings::detail::property_at<
639       MatrA, tag::data_structure >::type storage_f;
640 
641     assert (bindings::size_row (A) == bindings::size_row (A));
642     assert (bindings::size_column (A) == bindings::size (X));
643     assert (bindings::size_column (A) == bindings::size (B));
644 
645     return detail::solve (storage_f(), sys, A, X, B,
646                           Numeric.ptr, Control, Info);
647   }
648 
649   template <typename MatrA, typename VecX, typename VecB>
650   inline
solve(int sys,MatrA const & A,VecX & X,VecB const & B,numeric_type<typename bindings::value_type<MatrA>::type> const & Numeric,control_type<typename bindings::value_type<MatrA>::type> const & Control,info_type<typename bindings::value_type<MatrA>::type> & Info)651   int solve (int sys, MatrA const& A, VecX& X, VecB const& B,
652              numeric_type<
653                typename bindings::value_type<MatrA>::type
654              > const& Numeric,
655              control_type<
656                typename bindings::value_type<MatrA>::type
657              > const& Control,
658              info_type<
659                typename bindings::value_type<MatrA>::type
660              >& Info)
661   {
662     return solve (sys, A, X, B, Numeric, Control.ptr, Info.ptr);
663   }
664 
665   template <typename MatrA, typename VecX, typename VecB>
666   inline
solve(int sys,MatrA const & A,VecX & X,VecB const & B,numeric_type<typename bindings::value_type<MatrA>::type> const & Numeric,control_type<typename bindings::value_type<MatrA>::type> const & Control)667   int solve (int sys, MatrA const& A, VecX& X, VecB const& B,
668              numeric_type<
669                typename bindings::value_type<MatrA>::type
670              > const& Numeric,
671              control_type<
672                typename bindings::value_type<MatrA>::type
673              > const& Control)
674   {
675     return solve (sys, A, X, B, Numeric, Control.ptr);
676   }
677 
678   template <typename MatrA, typename VecX, typename VecB>
679   inline
solve(MatrA const & A,VecX & X,VecB const & B,numeric_type<typename bindings::value_type<MatrA>::type> const & Numeric,double const * Control=0,double * Info=0)680   int solve (MatrA const& A, VecX& X, VecB const& B,
681              numeric_type<
682                typename bindings::value_type<MatrA>::type
683              > const& Numeric,
684              double const* Control = 0, double* Info = 0)
685   {
686     return solve (UMFPACK_A, A, X, B, Numeric, Control, Info);
687   }
688 
689   template <typename MatrA, typename VecX, typename VecB>
690   inline
solve(MatrA const & A,VecX & X,VecB const & B,numeric_type<typename bindings::value_type<MatrA>::type> const & Numeric,control_type<typename bindings::value_type<MatrA>::type> const & Control,info_type<typename bindings::value_type<MatrA>::type> & Info)691   int solve (MatrA const& A, VecX& X, VecB const& B,
692              numeric_type<
693                typename bindings::value_type<MatrA>::type
694              > const& Numeric,
695              control_type<
696                typename bindings::value_type<MatrA>::type
697              > const& Control,
698              info_type<
699                typename bindings::value_type<MatrA>::type
700              >& Info)
701   {
702     return solve (UMFPACK_A, A, X, B, Numeric,
703                   Control.ptr, Info.ptr);
704   }
705 
706   template <typename MatrA, typename VecX, typename VecB>
707   inline
solve(MatrA const & A,VecX & X,VecB const & B,numeric_type<typename bindings::value_type<MatrA>::type> const & Numeric,control_type<typename bindings::value_type<MatrA>::type> const & Control)708   int solve (MatrA const& A, VecX& X, VecB const& B,
709              numeric_type<
710                typename bindings::value_type<MatrA>::type
711              > const& Numeric,
712              control_type<
713                typename bindings::value_type<MatrA>::type
714              > const& Control)
715   {
716     return solve (UMFPACK_A, A, X, B, Numeric, Control.ptr);
717   }
718 
719 
720   // umf_solve
721   /*
722    * symbolic, numeric and solve
723    */
724   namespace detail {
725 
726     template <typename MatrA, typename VecX, typename VecB>
727     inline
umf_solve(tag::compressed_sparse,MatrA const & A,VecX & X,VecB const & B,double const * Control=0,double * Info=0)728     int umf_solve (tag::compressed_sparse,
729                    MatrA const& A, VecX& X, VecB const& B,
730                    double const* Control = 0, double* Info = 0)
731     {
732 #ifdef CHECK_TEST_COVERAGE
733       typedef typename MatrA::not_yet_tested i_m_still_here;
734 #endif
735       symbolic_type<typename bindings::value_type<MatrA>::type>
736         Symbolic;
737       numeric_type<typename bindings::value_type<MatrA>::type>
738         Numeric;
739 
740       int status;
741       status = detail::symbolic (bindings::size_row (A),
742                                  bindings::size_column (A),
743                                  bindings::begin_compressed_index_major (A),
744                                  bindings::begin_index_minor (A),
745                                  bindings::begin_value (A),
746                                  &Symbolic.ptr, Control, Info);
747       if (status != UMFPACK_OK) return status;
748 
749       status = detail::numeric (bindings::size_row (A),
750                                 bindings::size_column (A),
751                                 bindings::begin_compressed_index_major (A),
752                                 bindings::begin_index_minor (A),
753                                 bindings::begin_value (A),
754                                 Symbolic.ptr, &Numeric.ptr, Control, Info);
755       if (status != UMFPACK_OK) return status;
756 
757       return detail::solve (UMFPACK_A, bindings::size_row (A),
758                             bindings::begin_compressed_index_major (A),
759                             bindings::begin_index_minor (A),
760                             bindings::begin_value (A),
761                             bindings::begin_value (X),
762                             bindings::begin_value (B),
763                             Numeric.ptr, Control, Info);
764     }
765 
766     template <typename MatrA, typename VecX, typename VecB>
767     inline
umf_solve(tag::coordinate_sparse,MatrA const & A,VecX & X,VecB const & B,double const * Control=0,double * Info=0)768     int umf_solve (tag::coordinate_sparse,
769                    MatrA const& A, VecX& X, VecB const& B,
770                    double const* Control = 0, double* Info = 0)
771     {
772 #ifdef CHECK_TEST_COVERAGE
773       typedef typename MatrAA::not_yet_tested i_m_still_here;
774 #endif
775       int n_row = bindings::size_row (A);
776       int n_col = bindings::size_column (A);
777       int nnz = bindings::end_value (A) - bindings::begin_value (A);
778 
779       typedef typename bindings::value_type<MatrA>::type val_t;
780 
781       int const* Ti = bindings::begin_index_minor (A);
782       int const* Tj = bindings::begin_index_major (A);
783       bindings::detail::array<int> Ap (n_col+1);
784       if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory;
785       bindings::detail::array<int> Ai (nnz);
786       if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory;
787 
788       int status = detail::triplet_to_col (n_row, n_col, nnz,
789                                            Ti, Tj, static_cast<val_t*> (0),
790                                            Ap.storage(), Ai.storage(),
791                                            static_cast<val_t*> (0), 0);
792       if (status != UMFPACK_OK) return status;
793 
794       symbolic_type<typename bindings::value_type<MatrA>::type>
795         Symbolic;
796       numeric_type<typename bindings::value_type<MatrA>::type>
797         Numeric;
798 
799       status = detail::symbolic (n_row, n_col,
800                                  Ap.storage(), Ai.storage(),
801                                  bindings::begin_value (A),
802                                  &Symbolic.ptr, Control, Info);
803       if (status != UMFPACK_OK) return status;
804 
805       status = detail::numeric (n_row, n_col,
806                                 Ap.storage(), Ai.storage(),
807                                 bindings::begin_value (A),
808                                 Symbolic.ptr, &Numeric.ptr, Control, Info);
809       if (status != UMFPACK_OK) return status;
810 
811       return detail::solve (UMFPACK_A, n_row, Ap.storage(), Ai.storage(),
812                             bindings::begin_value (A),
813                             bindings::begin_value (X),
814                             bindings::begin_value (B),
815                             Numeric.ptr, Control, Info);
816     }
817 
818   } // detail
819 
820   template <typename MatrA, typename VecX, typename VecB>
821   inline
umf_solve(MatrA const & A,VecX & X,VecB const & B,double const * Control=0,double * Info=0)822   int umf_solve (MatrA const& A, VecX& X, VecB const& B,
823                  double const* Control = 0, double* Info = 0)
824   {
825 #ifdef CHECK_TEST_COVERAGE
826       typedef typename MatrA::not_yet_tested i_m_still_here;
827 #endif
828 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
829     check_umfpack_structure<MatrA>();
830 #endif
831     typedef typename bindings::detail::property_at<
832       MatrA, tag::data_structure >::type storage_f;
833 
834     assert (bindings::size_row (A) == bindings::size_row (A));
835     assert (bindings::size_column (A) == bindings::size (X));
836     assert (bindings::size_column (A) == bindings::size (B));
837 
838     return detail::umf_solve (storage_f(), A, X, B, Control, Info);
839   }
840 
841   template <typename MatrA, typename VecX, typename VecB>
842   inline
umf_solve(MatrA const & A,VecX & X,VecB const & B,control_type<typename bindings::value_type<MatrA>::type> const & Control,info_type<typename bindings::value_type<MatrA>::type> & Info)843   int umf_solve (MatrA const& A, VecX& X, VecB const& B,
844                  control_type<
845                    typename bindings::value_type<MatrA>::type
846                  > const& Control,
847                  info_type<
848                    typename bindings::value_type<MatrA>::type
849                  >& Info)
850   {
851     return umf_solve (A, X, B, Control.ptr, Info.ptr);
852   }
853 
854   template <typename MatrA, typename VecX, typename VecB>
855   inline
umf_solve(MatrA const & A,VecX & X,VecB const & B,control_type<typename bindings::value_type<MatrA>::type> const & Control)856   int umf_solve (MatrA const& A, VecX& X, VecB const& B,
857                  control_type<
858                    typename bindings::value_type<MatrA>::type
859                  > const& Control)
860   {
861     return umf_solve (A, X, B, Control.ptr);
862   }
863 
864 
865   ///////////////////////
866   // matrix manipulations
867   ///////////////////////
868 
869 
870   // scale
871 
872   template <typename VecX, typename VecB>
873   inline
scale(VecX & X,VecB const & B,numeric_type<typename bindings::value_type<VecB>::type> const & Numeric)874   int scale (VecX& X, VecB const& B,
875              numeric_type<
876                typename bindings::value_type<VecB>::type
877              > const& Numeric)
878   {
879     return detail::scale (bindings::size (B),
880                           bindings::begin_value (X),
881                           bindings::begin_value (B),
882                           Numeric.ptr);
883   }
884 
885 
886   ////////////
887   // reporting
888   ////////////
889 
890 
891   // report status
892 
893   template <typename T>
894   inline
report_status(control_type<T> const & Control,int status)895   void report_status (control_type<T> const& Control, int status) {
896     detail::report_status (T(), 0, Control.ptr, status);
897   }
898 
899 #if 0
900   template <typename T>
901   inline
902   void report_status (int printing_level, int status) {
903     control_type<T> Control;
904     Control[UMFPACK_PRL] = printing_level;
905     detail::report_status (T(), 0, Control.ptr, status);
906   }
907   template <typename T>
908   inline
909   void report_status (int status) {
910     control_type<T> Control;
911     detail::report_status (T(), 0, Control.ptr, status);
912   }
913 #endif
914 
915 
916   // report control
917 
918   template <typename T>
919   inline
report_control(control_type<T> const & Control)920   void report_control (control_type<T> const& Control) {
921     detail::report_control (T(), 0, Control.ptr);
922   }
923 
924 
925   // report info
926 
927   template <typename T>
928   inline
report_info(control_type<T> const & Control,info_type<T> const & Info)929   void report_info (control_type<T> const& Control, info_type<T> const& Info) {
930     detail::report_info (T(), 0, Control.ptr, Info.ptr);
931   }
932 
933 #if 0
934   template <typename T>
935   inline
936   void report_info (int printing_level, info_type<T> const& Info) {
937     control_type<T> Control;
938     Control[UMFPACK_PRL] = printing_level;
939     detail::report_info (T(), 0, Control.ptr, Info.ptr);
940   }
941   template <typename T>
942   inline
943   void report_info (info_type<T> const& Info) {
944     control_type<T> Control;
945     detail::report_info (T(), 0, Control.ptr, Info.ptr);
946   }
947 #endif
948 
949 
950   // report matrix (compressed column and coordinate)
951 
952   namespace detail {
953 
954     template <typename MatrA>
955     inline
report_matrix(tag::compressed_sparse,MatrA const & A,double const * Control)956     int report_matrix (tag::compressed_sparse, MatrA const& A,
957                        double const* Control)
958     {
959       return detail::report_matrix (bindings::size_row (A),
960                                     bindings::size_column (A),
961                                     bindings::begin_compressed_index_major (A),
962                                     bindings::begin_index_minor (A),
963                                     bindings::begin_value (A),
964                                     1, Control);
965     }
966 
967     template <typename MatrA>
968     inline
report_matrix(tag::coordinate_sparse,MatrA const & A,double const * Control)969     int report_matrix (tag::coordinate_sparse, MatrA const& A,
970                        double const* Control)
971     {
972       return detail::report_triplet (bindings::size_row (A),
973                                      bindings::size_column (A),
974                                      bindings::end_value (A) - bindings::begin_value (A),
975                                      bindings::begin_index_major (A),
976                                      bindings::begin_index_minor (A),
977                                      bindings::begin_value (A),
978                                      Control);
979     }
980 
981   } // detail
982 
983   template <typename MatrA>
984   inline
report_matrix(MatrA const & A,control_type<typename bindings::value_type<MatrA>::type> const & Control)985   int report_matrix (MatrA const& A,
986                      control_type<
987                        typename bindings::value_type<MatrA>::type
988                      > const& Control)
989   {
990 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
991     check_umfpack_structure<MatrA>();
992 #endif
993     typedef typename bindings::detail::property_at<
994       MatrA, tag::data_structure >::type storage_f;
995 
996     return detail::report_matrix (storage_f(), A, Control.ptr);
997   }
998 
999 
1000   // report vector
1001 
1002   template <typename VecX>
1003   inline
report_vector(VecX const & X,control_type<typename bindings::value_type<VecX>::type> const & Control)1004   int report_vector (VecX const& X,
1005                      control_type<
1006                        typename bindings::value_type<VecX>::type
1007                      > const& Control)
1008   {
1009     return detail::report_vector (bindings::size (X),
1010                                   bindings::begin_value (X),
1011                                   Control.ptr);
1012   }
1013 
1014 
1015   // report numeric
1016 
1017   template <typename T>
1018   inline
report_numeric(numeric_type<T> const & Numeric,control_type<T> const & Control)1019   int report_numeric (numeric_type<T> const& Numeric,
1020                       control_type<T> const& Control)
1021   {
1022     return detail::report_numeric (T(), 0, Numeric.ptr, Control.ptr);
1023   }
1024 
1025 
1026   // report symbolic
1027 
1028   template <typename T>
1029   inline
report_symbolic(symbolic_type<T> const & Symbolic,control_type<T> const & Control)1030   int report_symbolic (symbolic_type<T> const& Symbolic,
1031                        control_type<T> const& Control)
1032   {
1033     return detail::report_symbolic (T(), 0, Symbolic.ptr, Control.ptr);
1034   }
1035 
1036 
1037   // report permutation vector
1038 
1039   template <typename VecP, typename T>
1040   inline
report_permutation(VecP const & Perm,control_type<T> const & Control)1041   int report_permutation (VecP const& Perm, control_type<T> const& Control) {
1042 #ifdef CHECK_TEST_COVERAGE
1043       typedef typename T::not_yet_tested i_m_still_here;
1044 #endif
1045     return detail::report_perm (T(), 0,
1046                                 bindings::begin_value (Perm),
1047                                 Control.ptr);
1048   }
1049 
1050 
1051 }}}}
1052 
1053 #endif // BOOST_NUMERIC_BINDINGS_UMFPACK_HPP
1054