umfpack.hpp

Go to the documentation of this file.
00001 /*
00002  * 
00003  * Copyright (c) Kresimir Fresl 2003 
00004  *
00005  * Distributed under the Boost Software License, Version 1.0.
00006  * (See accompanying file LICENSE_1_0.txt or copy at
00007  * http://www.boost.org/LICENSE_1_0.txt)
00008  *
00009  * Author acknowledges the support of the Faculty of Civil Engineering, 
00010  * University of Zagreb, Croatia.
00011  *
00012  */
00013 
00014 /* for UMFPACK Copyright, License and Availability see umfpack_inc.hpp */ 
00015 
00016 
00017 #ifndef BOOST_NUMERIC_BINDINGS_UMFPACK_HPP
00018 #define BOOST_NUMERIC_BINDINGS_UMFPACK_HPP
00019 
00020 
00021 #include <boost/noncopyable.hpp> 
00022 #include <boost/numeric/bindings/traits/traits.hpp>
00023 #include <boost/numeric/bindings/traits/sparse_traits.hpp>
00024 #include <boost/numeric/bindings/umfpack/umfpack_overloads.hpp>
00025 
00026 
00027 namespace boost { namespace numeric { namespace bindings {  namespace umfpack {
00028 
00029 
00030   template <typename T = double>
00031   struct symbolic_type : private noncopyable { 
00032     void *ptr; 
00033     ~symbolic_type() { 
00034       if (ptr)
00035         detail::free_symbolic (T(), 0, &ptr); 
00036     }
00037     void free() {
00038       if (ptr)
00039         detail::free_symbolic (T(), 0, &ptr); 
00040       ptr = 0; 
00041     }
00042   }; 
00043 
00044   template <typename T>
00045   void free (symbolic_type<T>& s) { s.free(); }
00046 
00047   template <typename T = double>
00048   struct numeric_type : private noncopyable { 
00049     void *ptr; 
00050     ~numeric_type() { 
00051       if (ptr)
00052         detail::free_numeric (T(), 0, &ptr); 
00053     }
00054     void free() { 
00055       if (ptr)
00056         detail::free_numeric (T(), 0, &ptr); 
00057       ptr = 0; 
00058     }
00059   }; 
00060 
00061   template <typename T>
00062   void free (numeric_type<T>& n) { n.free(); }
00063 
00064 
00065   template <typename T = double>
00066   struct control_type : private noncopyable {
00067     double ptr[UMFPACK_CONTROL]; 
00068     control_type() { detail::defaults (T(), 0, ptr); }
00069     double operator[] (int i) const { return ptr[i]; }
00070     double& operator[] (int i) { return ptr[i]; }
00071     void defaults() { detail::defaults (T(), 0, ptr); }
00072   }; 
00073 
00074   template <typename T>
00075   void defaults (control_type<T>& c) { c.defaults(); } 
00076 
00077   template <typename T = double>
00078   struct info_type : private noncopyable {
00079     double ptr[UMFPACK_INFO]; 
00080     double operator[] (int i) const { return ptr[i]; }
00081     double& operator[] (int i) { return ptr[i]; }
00082   }; 
00083 
00084 
00086   // solving system of linear equations
00088 
00089 
00090   // symbolic 
00091   /* 
00092    * Given nonzero pattern of a sparse matrix A in column-oriented form,
00093    * umfpack_*_symbolic performs a column pre-ordering to reduce fill-in
00094    * (using COLAMD or AMD) and a symbolic factorisation.  This is required
00095    * before the matrix can be numerically factorised with umfpack_*_numeric.
00096    */
00097   namespace detail {
00098 
00099     template <typename MatrA>
00100     inline
00101     int symbolic (traits::compressed_t, 
00102                   MatrA const& A, void **Symbolic, 
00103                   double const* Control = 0, double* Info = 0) 
00104     {
00105       return detail::symbolic (traits::spmatrix_size1 (A),
00106                                traits::spmatrix_size2 (A),
00107                                traits::spmatrix_index1_storage (A),
00108                                traits::spmatrix_index2_storage (A),
00109                                traits::spmatrix_value_storage (A),
00110                                Symbolic, Control, Info); 
00111     }
00112 
00113     template <typename MatrA, typename QVec>
00114     inline
00115     int symbolic (traits::compressed_t, 
00116                   MatrA const& A, QVec const& Qinit, void **Symbolic, 
00117                   double const* Control = 0, double* Info = 0) 
00118     {
00119       return detail::qsymbolic (traits::spmatrix_size1 (A),
00120                                 traits::spmatrix_size2 (A),
00121                                 traits::spmatrix_index1_storage (A),
00122                                 traits::spmatrix_index2_storage (A),
00123                                 traits::spmatrix_value_storage (A),
00124                                 traits::vector_storage (Qinit), 
00125                                 Symbolic, Control, Info); 
00126     }
00127 
00128     template <typename MatrA>
00129     inline
00130     int symbolic (traits::coordinate_t, 
00131                   MatrA const& A, void **Symbolic, 
00132                   double const* Control = 0, double* Info = 0) 
00133     {
00134       int n_row = traits::spmatrix_size1 (A); 
00135       int n_col = traits::spmatrix_size2 (A); 
00136       int nnz = traits::spmatrix_num_nonzeros (A); 
00137 
00138       typedef typename traits::sparse_matrix_traits<MatrA>::value_type val_t; 
00139 
00140       int const* Ti = traits::spmatrix_index2_storage (A);
00141       int const* Tj = traits::spmatrix_index1_storage (A); 
00142       traits::detail::array<int> Ap (n_col+1); 
00143       if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory;
00144       traits::detail::array<int> Ai (nnz); 
00145       if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory;
00146 
00147       int status = detail::triplet_to_col (n_row, n_col, nnz, 
00148                                            Ti, Tj, static_cast<val_t*> (0),
00149                                            Ap.storage(), Ai.storage(), 
00150                                            static_cast<val_t*> (0), 0); 
00151       if (status != UMFPACK_OK) return status; 
00152 
00153       return detail::symbolic (n_row, n_col, 
00154                                Ap.storage(), Ai.storage(),
00155                                traits::spmatrix_value_storage (A),
00156                                Symbolic, Control, Info); 
00157     }
00158 
00159     template <typename MatrA, typename QVec>
00160     inline
00161     int symbolic (traits::coordinate_t, 
00162                   MatrA const& A, QVec const& Qinit, void **Symbolic, 
00163                   double const* Control = 0, double* Info = 0) 
00164     {
00165       int n_row = traits::spmatrix_size1 (A); 
00166       int n_col = traits::spmatrix_size2 (A); 
00167       int nnz = traits::spmatrix_num_nonzeros (A); 
00168 
00169       typedef typename traits::sparse_matrix_traits<MatrA>::value_type val_t; 
00170 
00171       int const* Ti = traits::spmatrix_index2_storage (A);
00172       int const* Tj = traits::spmatrix_index1_storage (A); 
00173       traits::detail::array<int> Ap (n_col+1); 
00174       if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory;
00175       traits::detail::array<int> Ai (nnz); 
00176       if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory;
00177 
00178       int status = detail::triplet_to_col (n_row, n_col, nnz, 
00179                                            Ti, Tj, static_cast<val_t*> (0),
00180                                            Ap.storage(), Ai.storage(), 
00181                                            static_cast<val_t*> (0), 0); 
00182       if (status != UMFPACK_OK) return status; 
00183 
00184       return detail::qsymbolic (n_row, n_col, 
00185                                 Ap.storage(), Ai.storage(),
00186                                 traits::spmatrix_value_storage (A),
00187                                 traits::vector_storage (Qinit), 
00188                                 Symbolic, Control, Info); 
00189     }
00190 
00191   } // detail 
00192  
00193   template <typename MatrA>
00194   inline
00195   int symbolic (MatrA const& A, 
00196                 symbolic_type<
00197                   typename traits::sparse_matrix_traits<MatrA>::value_type
00198                 >& Symbolic, 
00199                 double const* Control = 0, double* Info = 0) 
00200   {
00201 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00202     BOOST_STATIC_ASSERT((boost::is_same<
00203       typename traits::sparse_matrix_traits<MatrA>::matrix_structure, 
00204       traits::general_t
00205     >::value)); 
00206     BOOST_STATIC_ASSERT((boost::is_same<
00207       typename traits::sparse_matrix_traits<MatrA>::ordering_type,
00208       traits::column_major_t
00209     >::value)); 
00210     BOOST_STATIC_ASSERT(traits::sparse_matrix_traits<MatrA>::index_base == 0);
00211 #endif 
00212 
00213     typedef 
00214       typename traits::sparse_matrix_traits<MatrA>::storage_format storage_f; 
00215 
00216 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00217     BOOST_STATIC_ASSERT(
00218       (boost::is_same<storage_f, traits::compressed_t>::value
00219        || 
00220        boost::is_same<storage_f, traits::coordinate_t>::value
00221        )); 
00222 #endif 
00223 
00224     return detail::symbolic (storage_f(), A, &Symbolic.ptr, Control, Info); 
00225   }
00226 
00227   template <typename MatrA>
00228   inline
00229   int symbolic (MatrA const& A, 
00230                 symbolic_type<
00231                   typename traits::sparse_matrix_traits<MatrA>::value_type
00232                 >& Symbolic, 
00233                 control_type<
00234                   typename traits::sparse_matrix_traits<MatrA>::value_type
00235                 > const& Control, 
00236                 info_type<
00237                   typename traits::sparse_matrix_traits<MatrA>::value_type
00238                 >& Info) 
00239   {
00240     return symbolic (A, Symbolic, Control.ptr, Info.ptr); 
00241   }
00242 
00243   template <typename MatrA>
00244   inline
00245   int symbolic (MatrA const& A, 
00246                 symbolic_type<
00247                   typename traits::sparse_matrix_traits<MatrA>::value_type
00248                 >& Symbolic, 
00249                 control_type<
00250                   typename traits::sparse_matrix_traits<MatrA>::value_type
00251                 > const& Control)
00252   {
00253     return symbolic (A, Symbolic, Control.ptr); 
00254   }
00255 
00256   template <typename MatrA, typename QVec>
00257   inline
00258   int symbolic (MatrA const& A, QVec const& Qinit, 
00259                 symbolic_type<
00260                   typename traits::sparse_matrix_traits<MatrA>::value_type
00261                 >& Symbolic, 
00262                 double const* Control = 0, double* Info = 0) 
00263   {
00264 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00265     BOOST_STATIC_ASSERT((boost::is_same<
00266       typename traits::sparse_matrix_traits<MatrA>::matrix_structure, 
00267       traits::general_t
00268     >::value)); 
00269     BOOST_STATIC_ASSERT((boost::is_same<
00270       typename traits::sparse_matrix_traits<MatrA>::ordering_type,
00271       traits::column_major_t
00272     >::value)); 
00273     BOOST_STATIC_ASSERT(traits::sparse_matrix_traits<MatrA>::index_base == 0);
00274 #endif 
00275 
00276     typedef 
00277       typename traits::sparse_matrix_traits<MatrA>::storage_format storage_f; 
00278 
00279 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00280     BOOST_STATIC_ASSERT(
00281       (boost::is_same<storage_f, traits::compressed_t>::value
00282        || 
00283        boost::is_same<storage_f, traits::coordinate_t>::value
00284        )); 
00285 #endif 
00286 
00287     assert (traits::spmatrix_size2 (A) == traits::vector_size (Qinit)); 
00288 
00289     return detail::symbolic (storage_f(), A, Qinit, 
00290                              &Symbolic.ptr, Control, Info); 
00291   }
00292 
00293   template <typename MatrA, typename QVec>
00294   inline
00295   int symbolic (MatrA const& A, QVec const& Qinit, 
00296                 symbolic_type<
00297                   typename traits::sparse_matrix_traits<MatrA>::value_type
00298                 >& Symbolic, 
00299                 control_type<
00300                   typename traits::sparse_matrix_traits<MatrA>::value_type
00301                 > const& Control, 
00302                 info_type<
00303                   typename traits::sparse_matrix_traits<MatrA>::value_type
00304                 >& Info) 
00305   {
00306     return symbolic (A, Qinit, Symbolic, Control.ptr, Info.ptr); 
00307   }
00308 
00309   template <typename MatrA, typename QVec>
00310   inline
00311   int symbolic (MatrA const& A, QVec const& Qinit,  
00312                 symbolic_type<
00313                   typename traits::sparse_matrix_traits<MatrA>::value_type
00314                 >& Symbolic, 
00315                 control_type<
00316                   typename traits::sparse_matrix_traits<MatrA>::value_type
00317                 > const& Control)
00318   {
00319     return symbolic (A, Qinit, Symbolic, Control.ptr); 
00320   }
00321 
00322 
00323   // numeric 
00324   /*
00325    * Given a sparse matrix A in column-oriented form, and a symbolic analysis
00326    * computed by umfpack_*_*symbolic, the umfpack_*_numeric routine performs 
00327    * the numerical factorisation, PAQ=LU, PRAQ=LU, or P(R\A)Q=LU, where P 
00328    * and Q are permutation matrices (represented as permutation vectors), 
00329    * R is the row scaling, L is unit-lower triangular, and U is upper 
00330    * triangular.  This is required before the system Ax=b (or other related 
00331    * linear systems) can be solved.  
00332    */
00333   namespace detail {
00334 
00335     template <typename MatrA>
00336     inline
00337     int numeric (traits::compressed_t, MatrA const& A, 
00338                  void *Symbolic, void** Numeric, 
00339                  double const* Control = 0, double* Info = 0) 
00340     {
00341       return detail::numeric (traits::spmatrix_size1 (A),
00342                               traits::spmatrix_size2 (A),
00343                               traits::spmatrix_index1_storage (A),
00344                               traits::spmatrix_index2_storage (A),
00345                               traits::spmatrix_value_storage (A),
00346                               Symbolic, Numeric, Control, Info); 
00347     }
00348 
00349     template <typename MatrA>
00350     inline
00351     int numeric (traits::coordinate_t, MatrA const& A, 
00352                  void *Symbolic, void** Numeric, 
00353                  double const* Control = 0, double* Info = 0) 
00354     {
00355       int n_row = traits::spmatrix_size1 (A); 
00356       int n_col = traits::spmatrix_size2 (A); 
00357       int nnz = traits::spmatrix_num_nonzeros (A); 
00358 
00359       typedef typename traits::sparse_matrix_traits<MatrA>::value_type val_t; 
00360 
00361       int const* Ti = traits::spmatrix_index2_storage (A);
00362       int const* Tj = traits::spmatrix_index1_storage (A); 
00363       traits::detail::array<int> Ap (n_col+1); 
00364       if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory;
00365       traits::detail::array<int> Ai (nnz); 
00366       if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory;
00367 
00368       int status = detail::triplet_to_col (n_row, n_col, nnz, 
00369                                            Ti, Tj, static_cast<val_t*> (0),
00370                                            Ap.storage(), Ai.storage(), 
00371                                            static_cast<val_t*> (0), 0); 
00372       if (status != UMFPACK_OK) return status; 
00373 
00374       return detail::numeric (n_row, n_col, 
00375                               Ap.storage(), Ai.storage(),
00376                               traits::spmatrix_value_storage (A),
00377                               Symbolic, Numeric, Control, Info); 
00378     }
00379 
00380   } // detail 
00381 
00382   template <typename MatrA>
00383   inline
00384   int numeric (MatrA const& A, 
00385                symbolic_type<
00386                  typename traits::sparse_matrix_traits<MatrA>::value_type
00387                > const& Symbolic, 
00388                numeric_type<
00389                  typename traits::sparse_matrix_traits<MatrA>::value_type
00390                >& Numeric, 
00391                double const* Control = 0, double* Info = 0) 
00392   {
00393 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00394     BOOST_STATIC_ASSERT((boost::is_same<
00395       typename traits::sparse_matrix_traits<MatrA>::matrix_structure, 
00396       traits::general_t
00397     >::value)); 
00398     BOOST_STATIC_ASSERT((boost::is_same<
00399       typename traits::sparse_matrix_traits<MatrA>::ordering_type,
00400       traits::column_major_t
00401     >::value)); 
00402     BOOST_STATIC_ASSERT(traits::sparse_matrix_traits<MatrA>::index_base == 0);
00403 #endif 
00404 
00405     typedef 
00406       typename traits::sparse_matrix_traits<MatrA>::storage_format storage_f; 
00407 
00408 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00409     BOOST_STATIC_ASSERT(
00410       (boost::is_same<storage_f, traits::compressed_t>::value
00411        || 
00412        boost::is_same<storage_f, traits::coordinate_t>::value
00413        )); 
00414 #endif 
00415 
00416     return detail::numeric (storage_f(), A, 
00417                             Symbolic.ptr, &Numeric.ptr, Control, Info); 
00418   }
00419 
00420   template <typename MatrA>
00421   inline
00422   int numeric (MatrA const& A, 
00423                symbolic_type<
00424                  typename traits::sparse_matrix_traits<MatrA>::value_type
00425                > const& Symbolic, 
00426                numeric_type<
00427                  typename traits::sparse_matrix_traits<MatrA>::value_type
00428                >& Numeric, 
00429                control_type<
00430                  typename traits::sparse_matrix_traits<MatrA>::value_type
00431                > const& Control, 
00432                info_type<
00433                  typename traits::sparse_matrix_traits<MatrA>::value_type
00434                >& Info) 
00435 
00436   {
00437     // g++ (3.2) is unable to distinguish 
00438     //           function numeric() and namespace boost::numeric ;o) 
00439     return umfpack::numeric (A, Symbolic, Numeric, Control.ptr, Info.ptr); 
00440   }
00441     
00442   template <typename MatrA>
00443   inline
00444   int numeric (MatrA const& A, 
00445                symbolic_type<
00446                  typename traits::sparse_matrix_traits<MatrA>::value_type
00447                > const& Symbolic, 
00448                numeric_type<
00449                  typename traits::sparse_matrix_traits<MatrA>::value_type
00450                >& Numeric, 
00451                control_type<
00452                  typename traits::sparse_matrix_traits<MatrA>::value_type
00453                > const& Control)
00454   {
00455     return umfpack::numeric (A, Symbolic, Numeric, Control.ptr); 
00456   }
00457     
00458 
00459   // factor 
00460   /* 
00461    * symbolic and numeric
00462    */
00463   namespace detail {
00464 
00465     template <typename MatrA>
00466     inline
00467     int factor (traits::compressed_t, MatrA const& A, 
00468                 void** Numeric, double const* Control = 0, double* Info = 0) 
00469     {
00470       symbolic_type<typename traits::sparse_matrix_traits<MatrA>::value_type>
00471         Symbolic; 
00472 
00473       int status;
00474       status = detail::symbolic (traits::spmatrix_size1 (A),
00475                                  traits::spmatrix_size2 (A),
00476                                  traits::spmatrix_index1_storage (A),
00477                                  traits::spmatrix_index2_storage (A),
00478                                  traits::spmatrix_value_storage (A),
00479                                  &Symbolic.ptr, Control, Info); 
00480       if (status != UMFPACK_OK) return status; 
00481 
00482       return detail::numeric (traits::spmatrix_size1 (A),
00483                               traits::spmatrix_size2 (A),
00484                               traits::spmatrix_index1_storage (A),
00485                               traits::spmatrix_index2_storage (A),
00486                               traits::spmatrix_value_storage (A),
00487                               Symbolic.ptr, Numeric, Control, Info); 
00488     }
00489 
00490     template <typename MatrA>
00491     inline
00492     int factor (traits::coordinate_t, MatrA const& A, 
00493                 void** Numeric, double const* Control = 0, double* Info = 0) 
00494     {
00495       int n_row = traits::spmatrix_size1 (A); 
00496       int n_col = traits::spmatrix_size2 (A); 
00497       int nnz = traits::spmatrix_num_nonzeros (A); 
00498 
00499       typedef typename traits::sparse_matrix_traits<MatrA>::value_type val_t; 
00500 
00501       int const* Ti = traits::spmatrix_index2_storage (A);
00502       int const* Tj = traits::spmatrix_index1_storage (A); 
00503       traits::detail::array<int> Ap (n_col+1); 
00504       if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory;
00505       traits::detail::array<int> Ai (nnz); 
00506       if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory;
00507 
00508       int status = detail::triplet_to_col (n_row, n_col, nnz, 
00509                                            Ti, Tj, static_cast<val_t*> (0),
00510                                            Ap.storage(), Ai.storage(), 
00511                                            static_cast<val_t*> (0), 0); 
00512       if (status != UMFPACK_OK) return status; 
00513 
00514       symbolic_type<typename traits::sparse_matrix_traits<MatrA>::value_type>
00515         Symbolic; 
00516 
00517       status = detail::symbolic (n_row, n_col, 
00518                                  Ap.storage(), Ai.storage(),
00519                                  traits::spmatrix_value_storage (A),
00520                                  &Symbolic.ptr, Control, Info); 
00521       if (status != UMFPACK_OK) return status; 
00522 
00523       return detail::numeric (n_row, n_col, 
00524                               Ap.storage(), Ai.storage(),
00525                               traits::spmatrix_value_storage (A),
00526                               Symbolic.ptr, Numeric, Control, Info); 
00527     }
00528 
00529   } // detail 
00530 
00531   template <typename MatrA>
00532   inline
00533   int factor (MatrA const& A, 
00534               numeric_type<
00535                 typename traits::sparse_matrix_traits<MatrA>::value_type
00536               >& Numeric, 
00537               double const* Control = 0, double* Info = 0) 
00538   {
00539 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00540     BOOST_STATIC_ASSERT((boost::is_same<
00541       typename traits::sparse_matrix_traits<MatrA>::matrix_structure, 
00542       traits::general_t
00543     >::value)); 
00544     BOOST_STATIC_ASSERT((boost::is_same<
00545       typename traits::sparse_matrix_traits<MatrA>::ordering_type,
00546       traits::column_major_t
00547     >::value)); 
00548     BOOST_STATIC_ASSERT(traits::sparse_matrix_traits<MatrA>::index_base == 0);
00549 #endif 
00550 
00551     typedef 
00552       typename traits::sparse_matrix_traits<MatrA>::storage_format storage_f; 
00553 
00554 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00555     BOOST_STATIC_ASSERT(
00556       (boost::is_same<storage_f, traits::compressed_t>::value
00557        || 
00558        boost::is_same<storage_f, traits::coordinate_t>::value
00559        )); 
00560 #endif 
00561 
00562     return detail::factor (storage_f(), A, &Numeric.ptr, Control, Info); 
00563   }
00564 
00565   template <typename MatrA>
00566   inline
00567   int factor (MatrA const& A, 
00568               numeric_type<
00569                 typename traits::sparse_matrix_traits<MatrA>::value_type
00570               >& Numeric, 
00571               control_type<
00572                 typename traits::sparse_matrix_traits<MatrA>::value_type
00573               > const& Control, 
00574               info_type<
00575                 typename traits::sparse_matrix_traits<MatrA>::value_type
00576               >& Info) 
00577   {
00578     return factor (A, Numeric, Control.ptr, Info.ptr); 
00579   }
00580     
00581   template <typename MatrA>
00582   inline
00583   int factor (MatrA const& A, 
00584               numeric_type<
00585                 typename traits::sparse_matrix_traits<MatrA>::value_type
00586               >& Numeric, 
00587               control_type<
00588                 typename traits::sparse_matrix_traits<MatrA>::value_type
00589               > const& Control)
00590   {
00591     return factor (A, Numeric, Control.ptr); 
00592   }
00593     
00594   
00595   // solve
00596   /*
00597    * Given LU factors computed by umfpack_*_numeric and the right-hand-side, 
00598    * B, solve a linear system for the solution X.  Iterative refinement is 
00599    * optionally performed.  Only square systems are handled. 
00600    */
00601   namespace detail {
00602 
00603     template <typename MatrA, typename VecX, typename VecB> 
00604     inline 
00605     int solve (traits::compressed_t, int sys, 
00606                MatrA const& A, VecX& X, VecB const& B, 
00607                void *Numeric, double const* Control = 0, double* Info = 0) 
00608     {
00609       return detail::solve (sys, traits::spmatrix_size1 (A),
00610                             traits::spmatrix_index1_storage (A),
00611                             traits::spmatrix_index2_storage (A),
00612                             traits::spmatrix_value_storage (A),
00613                             traits::vector_storage (X),
00614                             traits::vector_storage (B),
00615                             Numeric, Control, Info); 
00616     }
00617 
00618     template <typename MatrA, typename VecX, typename VecB> 
00619     inline 
00620     int solve (traits::coordinate_t, int sys, 
00621                MatrA const& A, VecX& X, VecB const& B, 
00622                void *Numeric, double const* Control = 0, double* Info = 0) 
00623     {
00624 
00625       int n = traits::spmatrix_size1 (A); 
00626       int nnz = traits::spmatrix_num_nonzeros (A); 
00627 
00628       typedef typename traits::sparse_matrix_traits<MatrA>::value_type val_t; 
00629 
00630       int const* Ti = traits::spmatrix_index2_storage (A);
00631       int const* Tj = traits::spmatrix_index1_storage (A); 
00632       traits::detail::array<int> Ap (n+1); 
00633       if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory;
00634       traits::detail::array<int> Ai (nnz); 
00635       if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory;
00636 
00637       int status = detail::triplet_to_col (n, n, nnz, 
00638                                            Ti, Tj, static_cast<val_t*> (0),
00639                                            Ap.storage(), Ai.storage(), 
00640                                            static_cast<val_t*> (0), 0); 
00641       if (status != UMFPACK_OK) return status; 
00642  
00643       return detail::solve (sys, n, Ap.storage(), Ai.storage(),
00644                             traits::spmatrix_value_storage (A),
00645                             traits::vector_storage (X),
00646                             traits::vector_storage (B),
00647                             Numeric, Control, Info); 
00648     }
00649 
00650   } // detail 
00651 
00652   template <typename MatrA, typename VecX, typename VecB> 
00653   inline 
00654   int solve (int sys, MatrA const& A, VecX& X, VecB const& B, 
00655              numeric_type<
00656                typename traits::sparse_matrix_traits<MatrA>::value_type
00657              > const& Numeric, 
00658              double const* Control = 0, double* Info = 0) 
00659   {
00660 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00661     BOOST_STATIC_ASSERT((boost::is_same<
00662       typename traits::sparse_matrix_traits<MatrA>::matrix_structure, 
00663       traits::general_t
00664     >::value)); 
00665     BOOST_STATIC_ASSERT((boost::is_same<
00666       typename traits::sparse_matrix_traits<MatrA>::ordering_type,
00667       traits::column_major_t
00668     >::value)); 
00669     BOOST_STATIC_ASSERT(traits::sparse_matrix_traits<MatrA>::index_base == 0);
00670 #endif 
00671 
00672     typedef 
00673       typename traits::sparse_matrix_traits<MatrA>::storage_format storage_f; 
00674 
00675 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00676     BOOST_STATIC_ASSERT(
00677       (boost::is_same<storage_f, traits::compressed_t>::value
00678        || 
00679        boost::is_same<storage_f, traits::coordinate_t>::value
00680        )); 
00681 #endif 
00682 
00683     assert (traits::spmatrix_size1 (A) == traits::spmatrix_size1 (A)); 
00684     assert (traits::spmatrix_size2 (A) == traits::vector_size (X)); 
00685     assert (traits::spmatrix_size2 (A) == traits::vector_size (B)); 
00686 
00687     return detail::solve (storage_f(), sys, A, X, B, 
00688                           Numeric.ptr, Control, Info); 
00689   }
00690 
00691   template <typename MatrA, typename VecX, typename VecB> 
00692   inline 
00693   int solve (int sys, MatrA const& A, VecX& X, VecB const& B, 
00694              numeric_type<
00695                typename traits::sparse_matrix_traits<MatrA>::value_type
00696              > const& Numeric, 
00697              control_type<
00698                typename traits::sparse_matrix_traits<MatrA>::value_type
00699              > const& Control, 
00700              info_type<
00701                typename traits::sparse_matrix_traits<MatrA>::value_type
00702              >& Info) 
00703   {
00704     return solve (sys, A, X, B, Numeric, Control.ptr, Info.ptr); 
00705   }
00706 
00707   template <typename MatrA, typename VecX, typename VecB> 
00708   inline 
00709   int solve (int sys, MatrA const& A, VecX& X, VecB const& B, 
00710              numeric_type<
00711                typename traits::sparse_matrix_traits<MatrA>::value_type
00712              > const& Numeric, 
00713              control_type<
00714                typename traits::sparse_matrix_traits<MatrA>::value_type
00715              > const& Control)
00716   {
00717     return solve (sys, A, X, B, Numeric, Control.ptr); 
00718   }
00719 
00720   template <typename MatrA, typename VecX, typename VecB> 
00721   inline 
00722   int solve (MatrA const& A, VecX& X, VecB const& B, 
00723              numeric_type<
00724                typename traits::sparse_matrix_traits<MatrA>::value_type
00725              > const& Numeric, 
00726              double const* Control = 0, double* Info = 0) 
00727   {
00728     return solve (UMFPACK_A, A, X, B, Numeric, Control, Info); 
00729   }
00730 
00731   template <typename MatrA, typename VecX, typename VecB> 
00732   inline 
00733   int solve (MatrA const& A, VecX& X, VecB const& B, 
00734              numeric_type<
00735                typename traits::sparse_matrix_traits<MatrA>::value_type
00736              > const& Numeric, 
00737              control_type<
00738                typename traits::sparse_matrix_traits<MatrA>::value_type
00739              > const& Control, 
00740              info_type<
00741                typename traits::sparse_matrix_traits<MatrA>::value_type
00742              >& Info) 
00743   {
00744     return solve (UMFPACK_A, A, X, B, Numeric, 
00745                   Control.ptr, Info.ptr); 
00746   }
00747 
00748   template <typename MatrA, typename VecX, typename VecB> 
00749   inline 
00750   int solve (MatrA const& A, VecX& X, VecB const& B, 
00751              numeric_type<
00752                typename traits::sparse_matrix_traits<MatrA>::value_type
00753              > const& Numeric, 
00754              control_type<
00755                typename traits::sparse_matrix_traits<MatrA>::value_type
00756              > const& Control)
00757   {
00758     return solve (UMFPACK_A, A, X, B, Numeric, Control.ptr); 
00759   }
00760 
00761 
00762   // umf_solve 
00763   /* 
00764    * symbolic, numeric and solve 
00765    */
00766   namespace detail {
00767 
00768     template <typename MatrA, typename VecX, typename VecB>
00769     inline
00770     int umf_solve (traits::compressed_t, 
00771                    MatrA const& A, VecX& X, VecB const& B, 
00772                    double const* Control = 0, double* Info = 0) 
00773     {
00774       symbolic_type<typename traits::sparse_matrix_traits<MatrA>::value_type>
00775         Symbolic; 
00776       numeric_type<typename traits::sparse_matrix_traits<MatrA>::value_type>
00777         Numeric; 
00778 
00779       int status;
00780       status = detail::symbolic (traits::spmatrix_size1 (A),
00781                                  traits::spmatrix_size2 (A),
00782                                  traits::spmatrix_index1_storage (A),
00783                                  traits::spmatrix_index2_storage (A),
00784                                  traits::spmatrix_value_storage (A),
00785                                  &Symbolic.ptr, Control, Info); 
00786       if (status != UMFPACK_OK) return status; 
00787 
00788       status = detail::numeric (traits::spmatrix_size1 (A),
00789                                 traits::spmatrix_size2 (A),
00790                                 traits::spmatrix_index1_storage (A),
00791                                 traits::spmatrix_index2_storage (A),
00792                                 traits::spmatrix_value_storage (A),
00793                                 Symbolic.ptr, &Numeric.ptr, Control, Info); 
00794       if (status != UMFPACK_OK) return status; 
00795 
00796       return detail::solve (UMFPACK_A, traits::spmatrix_size1 (A),
00797                             traits::spmatrix_index1_storage (A),
00798                             traits::spmatrix_index2_storage (A),
00799                             traits::spmatrix_value_storage (A),
00800                             traits::vector_storage (X),
00801                             traits::vector_storage (B),
00802                             Numeric.ptr, Control, Info); 
00803     }
00804 
00805     template <typename MatrA, typename VecX, typename VecB>
00806     inline
00807     int umf_solve (traits::coordinate_t, 
00808                    MatrA const& A, VecX& X, VecB const& B, 
00809                    double const* Control = 0, double* Info = 0) 
00810     {
00811       int n_row = traits::spmatrix_size1 (A); 
00812       int n_col = traits::spmatrix_size2 (A); 
00813       int nnz = traits::spmatrix_num_nonzeros (A); 
00814 
00815       typedef typename traits::sparse_matrix_traits<MatrA>::value_type val_t; 
00816 
00817       int const* Ti = traits::spmatrix_index2_storage (A);
00818       int const* Tj = traits::spmatrix_index1_storage (A); 
00819       traits::detail::array<int> Ap (n_col+1); 
00820       if (!Ap.valid()) return UMFPACK_ERROR_out_of_memory;
00821       traits::detail::array<int> Ai (nnz); 
00822       if (!Ai.valid()) return UMFPACK_ERROR_out_of_memory;
00823 
00824       int status = detail::triplet_to_col (n_row, n_col, nnz, 
00825                                            Ti, Tj, static_cast<val_t*> (0),
00826                                            Ap.storage(), Ai.storage(), 
00827                                            static_cast<val_t*> (0), 0); 
00828       if (status != UMFPACK_OK) return status; 
00829 
00830       symbolic_type<typename traits::sparse_matrix_traits<MatrA>::value_type>
00831         Symbolic; 
00832       numeric_type<typename traits::sparse_matrix_traits<MatrA>::value_type>
00833         Numeric; 
00834 
00835       status = detail::symbolic (n_row, n_col, 
00836                                  Ap.storage(), Ai.storage(),
00837                                  traits::spmatrix_value_storage (A),
00838                                  &Symbolic.ptr, Control, Info); 
00839       if (status != UMFPACK_OK) return status; 
00840 
00841       status = detail::numeric (n_row, n_col, 
00842                                 Ap.storage(), Ai.storage(),
00843                                 traits::spmatrix_value_storage (A),
00844                                 Symbolic.ptr, &Numeric.ptr, Control, Info); 
00845       if (status != UMFPACK_OK) return status; 
00846 
00847       return detail::solve (UMFPACK_A, n_row, Ap.storage(), Ai.storage(),
00848                             traits::spmatrix_value_storage (A),
00849                             traits::vector_storage (X),
00850                             traits::vector_storage (B),
00851                             Numeric.ptr, Control, Info); 
00852     }
00853 
00854   } // detail 
00855 
00856   template <typename MatrA, typename VecX, typename VecB>
00857   inline
00858   int umf_solve (MatrA const& A, VecX& X, VecB const& B, 
00859                  double const* Control = 0, double* Info = 0) 
00860   {
00861 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00862     BOOST_STATIC_ASSERT((boost::is_same<
00863       typename traits::sparse_matrix_traits<MatrA>::matrix_structure, 
00864       traits::general_t
00865     >::value)); 
00866     BOOST_STATIC_ASSERT((boost::is_same<
00867       typename traits::sparse_matrix_traits<MatrA>::ordering_type,
00868       traits::column_major_t
00869     >::value)); 
00870     BOOST_STATIC_ASSERT(traits::sparse_matrix_traits<MatrA>::index_base == 0);
00871 #endif 
00872 
00873     typedef 
00874       typename traits::sparse_matrix_traits<MatrA>::storage_format storage_f; 
00875 
00876 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00877     BOOST_STATIC_ASSERT(
00878       (boost::is_same<storage_f, traits::compressed_t>::value
00879        || 
00880        boost::is_same<storage_f, traits::coordinate_t>::value
00881        )); 
00882 #endif 
00883 
00884     assert (traits::spmatrix_size1 (A) == traits::spmatrix_size1 (A)); 
00885     assert (traits::spmatrix_size2 (A) == traits::vector_size (X)); 
00886     assert (traits::spmatrix_size2 (A) == traits::vector_size (B)); 
00887 
00888     return detail::umf_solve (storage_f(), A, X, B, Control, Info); 
00889   }
00890 
00891   template <typename MatrA, typename VecX, typename VecB>
00892   inline
00893   int umf_solve (MatrA const& A, VecX& X, VecB const& B, 
00894                  control_type<
00895                    typename traits::sparse_matrix_traits<MatrA>::value_type
00896                  > const& Control, 
00897                  info_type<
00898                    typename traits::sparse_matrix_traits<MatrA>::value_type
00899                  >& Info) 
00900   {
00901     return umf_solve (A, X, B, Control.ptr, Info.ptr); 
00902   }    
00903 
00904   template <typename MatrA, typename VecX, typename VecB>
00905   inline
00906   int umf_solve (MatrA const& A, VecX& X, VecB const& B, 
00907                  control_type<
00908                    typename traits::sparse_matrix_traits<MatrA>::value_type
00909                  > const& Control)
00910   {
00911     return umf_solve (A, X, B, Control.ptr); 
00912   }    
00913 
00914 
00916   // matrix manipulations
00918 
00919 
00920   // scale 
00921   
00922   template <typename VecX, typename VecB> 
00923   inline 
00924   int scale (VecX& X, VecB const& B, 
00925              numeric_type<
00926                typename traits::vector_traits<VecB>::value_type
00927              > const& Numeric) 
00928   {
00929     return detail::scale (traits::vector_size (B),
00930                           traits::vector_storage (X),
00931                           traits::vector_storage (B),
00932                           Numeric.ptr);
00933   }
00934 
00935 
00937   // reporting
00939 
00940 
00941   // report status
00942 
00943   template <typename T>
00944   inline
00945   void report_status (control_type<T> const& Control, int status) {
00946     detail::report_status (T(), 0, Control.ptr, status); 
00947   }
00948 
00949 #if 0
00950   template <typename T>
00951   inline
00952   void report_status (int printing_level, int status) {
00953     control_type<T> Control; 
00954     Control[UMFPACK_PRL] = printing_level; 
00955     detail::report_status (T(), 0, Control.ptr, status); 
00956   }
00957   template <typename T>
00958   inline
00959   void report_status (int status) {
00960     control_type<T> Control; 
00961     detail::report_status (T(), 0, Control.ptr, status); 
00962   }
00963 #endif 
00964   
00965 
00966   // report control
00967 
00968   template <typename T>
00969   inline
00970   void report_control (control_type<T> const& Control) {
00971     detail::report_control (T(), 0, Control.ptr); 
00972   }
00973   
00974 
00975   // report info 
00976 
00977   template <typename T>
00978   inline
00979   void report_info (control_type<T> const& Control, info_type<T> const& Info) {
00980     detail::report_info (T(), 0, Control.ptr, Info.ptr); 
00981   }
00982 
00983 #if 0
00984   template <typename T>
00985   inline
00986   void report_info (int printing_level, info_type<T> const& Info) {
00987     control_type<T> Control; 
00988     Control[UMFPACK_PRL] = printing_level; 
00989     detail::report_info (T(), 0, Control.ptr, Info.ptr); 
00990   }
00991   template <typename T>
00992   inline
00993   void report_info (info_type<T> const& Info) {
00994     control_type<T> Control; 
00995     detail::report_info (T(), 0, Control.ptr, Info.ptr); 
00996   }
00997 #endif 
00998 
00999 
01000   // report matrix (compressed column and coordinate) 
01001 
01002   namespace detail {
01003 
01004     template <typename MatrA>
01005     inline
01006     int report_matrix (traits::compressed_t, MatrA const& A, 
01007                        double const* Control)
01008     {
01009       return detail::report_matrix (traits::spmatrix_size1 (A),
01010                                     traits::spmatrix_size2 (A),
01011                                     traits::spmatrix_index1_storage (A),
01012                                     traits::spmatrix_index2_storage (A),
01013                                     traits::spmatrix_value_storage (A),
01014                                     1, Control); 
01015     }
01016     
01017     template <typename MatrA>
01018     inline
01019     int report_matrix (traits::coordinate_t, MatrA const& A, 
01020                        double const* Control)
01021     {
01022       return detail::report_triplet (traits::spmatrix_size1 (A),
01023                                      traits::spmatrix_size2 (A),
01024                                      traits::spmatrix_num_nonzeros (A), 
01025                                      traits::spmatrix_index1_storage (A),
01026                                      traits::spmatrix_index2_storage (A),
01027                                      traits::spmatrix_value_storage (A),
01028                                      Control); 
01029     }
01030     
01031   } // detail 
01032 
01033   template <typename MatrA>
01034   inline
01035   int report_matrix (MatrA const& A, 
01036                      control_type<
01037                        typename traits::sparse_matrix_traits<MatrA>::value_type
01038                      > const& Control) 
01039   {
01040 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
01041     BOOST_STATIC_ASSERT((boost::is_same<
01042       typename traits::sparse_matrix_traits<MatrA>::matrix_structure, 
01043       traits::general_t
01044     >::value)); 
01045     BOOST_STATIC_ASSERT((boost::is_same<
01046       typename traits::sparse_matrix_traits<MatrA>::ordering_type,
01047       traits::column_major_t
01048     >::value)); 
01049     BOOST_STATIC_ASSERT(traits::sparse_matrix_traits<MatrA>::index_base == 0);
01050 #endif 
01051 
01052     typedef 
01053       typename traits::sparse_matrix_traits<MatrA>::storage_format storage_f; 
01054 
01055 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
01056     BOOST_STATIC_ASSERT(
01057       (boost::is_same<storage_f, traits::compressed_t>::value
01058        || 
01059        boost::is_same<storage_f, traits::coordinate_t>::value
01060        )); 
01061 #endif 
01062 
01063     return detail::report_matrix (storage_f(), A, Control.ptr); 
01064   }
01065 
01066 
01067   // report vector 
01068 
01069   template <typename VecX>
01070   inline
01071   int report_vector (VecX const& X, 
01072                      control_type<
01073                        typename traits::vector_traits<VecX>::value_type
01074                      > const& Control) 
01075   {
01076     return detail::report_vector (traits::vector_size (X), 
01077                                   traits::vector_storage (X), 
01078                                   Control.ptr);
01079   }
01080 
01081 
01082   // report numeric 
01083 
01084   template <typename T> 
01085   inline
01086   int report_numeric (numeric_type<T> const& Numeric, 
01087                       control_type<T> const& Control)
01088   {
01089     return detail::report_numeric (T(), 0, Numeric.ptr, Control.ptr); 
01090   }
01091 
01092 
01093   // report symbolic 
01094 
01095   template <typename T> 
01096   inline
01097   int report_symbolic (symbolic_type<T> const& Symbolic, 
01098                        control_type<T> const& Control)
01099   {
01100     return detail::report_symbolic (T(), 0, Symbolic.ptr, Control.ptr); 
01101   }
01102 
01103 
01104   // report permutation vector 
01105 
01106   template <typename VecP, typename T> 
01107   inline
01108   int report_permutation (VecP const& Perm, control_type<T> const& Control) {
01109     return detail::report_perm (T(), 0, 
01110                                 traits::vector_storage (Perm),
01111                                 Control.ptr); 
01112   }
01113 
01114 
01115 }}}} 
01116 
01117 #endif // BOOST_NUMERIC_BINDINGS_UMFPACK_HPP