umfpack_overloads.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_OVERLOADS_HPP
00018 #define BOOST_NUMERIC_BINDINGS_UMFPACK_OVERLOADS_HPP
00019 
00020 #include <boost/numeric/bindings/umfpack/umfpack_inc.hpp>
00021 #include <boost/numeric/bindings/traits/type.hpp>
00022 #include <boost/numeric/bindings/traits/detail/array.hpp>
00023 #include <boost/numeric/bindings/traits/detail/utils.hpp>
00024 
00025 namespace boost { namespace numeric { namespace bindings { 
00026 
00027   namespace umfpack { namespace detail {
00028 
00030     // UMFPACK primary routines:
00032 
00033     // symbolic 
00034 
00035     inline
00036     int symbolic (int n_row, int n_col, 
00037                   int const* Ap, int const* Ai, double const* Ax, 
00038                   void **Symbolic, double const* Control, double* Info) 
00039     {
00040       return umfpack_di_symbolic (n_row, n_col, Ap, Ai, Ax, 
00041                                   Symbolic, Control, Info); 
00042     }
00043 
00044     inline
00045     int symbolic (int n_row, int n_col,
00046                   int const* Ap, int const* Ai, traits::complex_d const* Ax, 
00047                   void **Symbolic, double const* Control, double* Info)
00048     {
00049       int nnz = Ap[n_col];
00050       traits::detail::array<double> Axr (nnz); 
00051       if (!Axr.valid()) return UMFPACK_ERROR_out_of_memory;
00052       traits::detail::array<double> Axi (nnz); 
00053       if (!Axi.valid()) return UMFPACK_ERROR_out_of_memory;
00054       traits::detail::disentangle (Ax, Ax+nnz, 
00055                                    Axr.storage(), Axi.storage()); 
00056       return umfpack_zi_symbolic (n_row, n_col, Ap, Ai, 
00057                                   Axr.storage(), Axi.storage(), 
00058                                   Symbolic, Control, Info); 
00059     }
00060 
00061 
00062     // numeric 
00063 
00064     inline
00065     int numeric (int, int, 
00066                  int const* Ap, int const* Ai, double const* Ax, 
00067                  void *Symbolic, void **Numeric, 
00068                  double const* Control, double* Info) 
00069     {
00070       return umfpack_di_numeric (Ap, Ai, Ax, Symbolic, Numeric, Control, Info);
00071     }
00072 
00073     inline
00074     int numeric (int, int n_col, 
00075                  int const* Ap, int const* Ai, traits::complex_d const* Ax, 
00076                  void *Symbolic, void **Numeric,
00077                  double const* Control, double* Info)
00078     {
00079       int nnz = Ap[n_col];
00080       traits::detail::array<double> Axr (nnz); 
00081       if (!Axr.valid()) return UMFPACK_ERROR_out_of_memory;
00082       traits::detail::array<double> Axi (nnz); 
00083       if (!Axi.valid()) return UMFPACK_ERROR_out_of_memory;
00084       traits::detail::disentangle (Ax, Ax+nnz, 
00085                                    Axr.storage(), Axi.storage()); 
00086       return umfpack_zi_numeric (Ap, Ai, Axr.storage(), Axi.storage(), 
00087                                  Symbolic, Numeric, Control, Info); 
00088     }
00089 
00090 
00091     // solve 
00092 
00093     inline
00094     int solve (int sys, int, 
00095                int const* Ap, int const* Ai, double const* Ax,
00096                double* X, double const* B, void *Numeric,
00097                double const* Control, double* Info)
00098     {
00099       return umfpack_di_solve (sys, Ap, Ai, Ax, X, B, Numeric, Control, Info);
00100     }
00101 
00102     inline
00103     int solve (int sys, int n, 
00104                int const* Ap, int const* Ai, traits::complex_d const* Ax, 
00105                traits::complex_d* X, traits::complex_d const* B, 
00106                void *Numeric, double const* Control, double* Info)
00107     {
00108       int nnz = Ap[n];
00109       traits::detail::array<double> Axr (nnz); 
00110       if (!Axr.valid()) return UMFPACK_ERROR_out_of_memory;
00111       traits::detail::array<double> Axi (nnz); 
00112       if (!Axi.valid()) return UMFPACK_ERROR_out_of_memory;
00113       traits::detail::disentangle (Ax, Ax+nnz, 
00114                                    Axr.storage(), Axi.storage()); 
00115       traits::detail::array<double> Br (n); 
00116       if (!Br.valid()) return UMFPACK_ERROR_out_of_memory;
00117       traits::detail::array<double> Bi (n); 
00118       if (!Bi.valid()) return UMFPACK_ERROR_out_of_memory;
00119       traits::detail::disentangle (B, B+n, 
00120                                    Br.storage(), Bi.storage()); 
00121       traits::detail::array<double> Xr (n); 
00122       if (!Xr.valid()) return UMFPACK_ERROR_out_of_memory;
00123       traits::detail::array<double> Xi (n); 
00124       if (!Xi.valid()) return UMFPACK_ERROR_out_of_memory;
00125 
00126       int status = umfpack_zi_solve (sys, Ap, Ai, 
00127                                      Axr.storage(), Axi.storage(), 
00128                                      Xr.storage(), Xi.storage(), 
00129                                      Br.storage(), Bi.storage(), 
00130                                      Numeric, Control, Info);
00131       if (status != UMFPACK_OK) return status; 
00132       traits::detail::interlace (Xr.storage(), Xr.storage() + n,
00133                                  Xi.storage(), X); 
00134       return status; 
00135     }
00136 
00137 
00138     // free_symbolic 
00139 
00140     inline
00141     void free_symbolic (double, int, void **Symbolic) {
00142       umfpack_di_free_symbolic (Symbolic); 
00143     }
00144     inline
00145     void free_symbolic (traits::complex_d const&, int, void **Symbolic) {
00146       umfpack_zi_free_symbolic (Symbolic); 
00147     }
00148 
00149 
00150     // free_numeric
00151 
00152     inline
00153     void free_numeric (double, int, void **Numeric) {
00154       umfpack_di_free_numeric (Numeric); 
00155     }
00156     inline
00157     void free_numeric (traits::complex_d const&, int, void **Numeric) {
00158       umfpack_zi_free_numeric (Numeric); 
00159     }
00160 
00161 
00163     // UMFPACK alternative routines:
00165 
00166     // default control 
00167 
00168     inline
00169     void defaults (double, int, double* Control) {
00170       umfpack_di_defaults (Control); 
00171     }
00172     inline
00173     void defaults (traits::complex_d const&, int, double* Control) {
00174       umfpack_zi_defaults (Control); 
00175     }
00176 
00177 
00178     // symbolic with specified column preordering
00179 
00180     inline
00181     int qsymbolic (int n_row, int n_col, 
00182                    int const* Ap, int const* Ai, double const* Ax, 
00183                    int const* Qinit, 
00184                    void **Symbolic, double const* Control, double* Info) 
00185     {
00186       return umfpack_di_qsymbolic (n_row, n_col, Ap, Ai, Ax, Qinit, 
00187                                   Symbolic, Control, Info); 
00188     }
00189 
00190     inline
00191     int qsymbolic (int n_row, int n_col,
00192                    int const* Ap, int const* Ai, traits::complex_d const* Ax, 
00193                    int const* Qinit, 
00194                    void **Symbolic, double const* Control, double* Info)
00195     {
00196       int nnz = Ap[n_col];
00197       traits::detail::array<double> Axr (nnz); 
00198       if (!Axr.valid()) return UMFPACK_ERROR_out_of_memory;
00199       traits::detail::array<double> Axi (nnz); 
00200       if (!Axi.valid()) return UMFPACK_ERROR_out_of_memory;
00201       traits::detail::disentangle (Ax, Ax+nnz, 
00202                                    Axr.storage(), Axi.storage()); 
00203       return umfpack_zi_qsymbolic (n_row, n_col, Ap, Ai, 
00204                                    Axr.storage(), Axi.storage(), 
00205                                    Qinit, Symbolic, Control, Info); 
00206     }
00207 
00208 
00210     // UMFPACK matrix manipulation routines
00212 
00213     // triplet (coordinate) to compressed column 
00214 
00215     inline
00216     int triplet_to_col (int n_row, int n_col, int nz,
00217                         int const* Ti, int const* Tj, double const* Tx, 
00218                         int* Ap, int* Ai, double* Ax, int* Map) 
00219     {
00220       return umfpack_di_triplet_to_col (n_row, n_col, nz, 
00221                                         Ti, Tj, Tx,
00222                                         Ap, Ai, Ax, Map);
00223     }
00224 
00225     inline 
00226     int triplet_to_col (int n_row, int n_col, int nz,
00227                         int const* Ti, int const* Tj, 
00228                         traits::complex_d const* Tx, 
00229                         int* Ap, int* Ai, traits::complex_d* Ax, 
00230                         int* Map)
00231     {
00232       assert (Tx == 0 && Ax == 0 || Tx != 0 && Ax != 0); 
00233       double *Txr = 0, *Txi = 0;
00234       if (Tx != 0) {
00235         traits::detail::array<double> ATxr (nz); 
00236         if (!ATxr.valid()) return UMFPACK_ERROR_out_of_memory;
00237         traits::detail::array<double> ATxi (nz); 
00238         if (!ATxi.valid()) return UMFPACK_ERROR_out_of_memory;
00239         Txr = ATxr.storage();
00240         Txi = ATxi.storage(); 
00241         traits::detail::disentangle (Tx, Tx+nz, Txr, Txi); 
00242       }
00243       double *Axr = 0, *Axi = 0;
00244       if (Ax != 0) {
00245         traits::detail::array<double> AAxr (nz); 
00246         if (!AAxr.valid()) return UMFPACK_ERROR_out_of_memory;
00247         traits::detail::array<double> AAxi (nz); 
00248         if (!AAxi.valid()) return UMFPACK_ERROR_out_of_memory;
00249         Axr = AAxr.storage();
00250         Axi = AAxi.storage(); 
00251       } 
00252       int status; 
00253       status = umfpack_zi_triplet_to_col (n_row, n_col, nz, 
00254                                           Ti, Tj, Txr, Txi,
00255                                           Ap, Ai, Axr, Axi, Map);
00256       if (Ax != 0) {
00257         if (status != UMFPACK_OK) return status; 
00258         traits::detail::interlace (Axr, Axr + nz, Axi, Ax); 
00259       }
00260       return status; 
00261     }
00262 
00263 
00264     // scale
00265 
00266     inline
00267     int scale (int, double* X, double const* B, void* Numeric) {
00268       return umfpack_di_scale (X, B, Numeric); 
00269     } 
00270 
00271     inline
00272     int scale (int n, traits::complex_d* X, 
00273                traits::complex_d const* B, void* Numeric) 
00274     {
00275       traits::detail::array<double> Br (n); 
00276       if (!Br.valid()) return UMFPACK_ERROR_out_of_memory;
00277       traits::detail::array<double> Bi (n); 
00278       if (!Bi.valid()) return UMFPACK_ERROR_out_of_memory;
00279       traits::detail::disentangle (B, B+n, 
00280                                    Br.storage(), Bi.storage()); 
00281       traits::detail::array<double> Xr (n); 
00282       if (!Xr.valid()) return UMFPACK_ERROR_out_of_memory;
00283       traits::detail::array<double> Xi (n); 
00284       if (!Xi.valid()) return UMFPACK_ERROR_out_of_memory;
00285 
00286       int status = umfpack_zi_scale (Xr.storage(), Xi.storage(), 
00287                                      Br.storage(), Bi.storage(), 
00288                                      Numeric); 
00289       if (status != UMFPACK_OK) return status; 
00290       traits::detail::interlace (Xr.storage(), Xr.storage() + n,
00291                                  Xi.storage(), X); 
00292       return status; 
00293     } 
00294 
00295 
00297     // UMFPACK reporting routines:
00299 
00300     // report status
00301 
00302     inline
00303     void report_status (double, int, double const* Control, int status) {
00304       umfpack_di_report_status (Control, status); 
00305     }
00306     inline
00307     void report_status (traits::complex_d const&, int, 
00308                         double const* Control, int status) 
00309     {
00310       umfpack_zi_report_status (Control, status); 
00311     }
00312 
00313 
00314     // report control
00315 
00316     inline
00317     void report_control (double, int, double const* Control) {
00318       umfpack_di_report_control (Control); 
00319     }
00320     inline
00321     void 
00322     report_control (traits::complex_d const&, int, double const* Control) {
00323       umfpack_zi_report_control (Control); 
00324     }
00325 
00326 
00327     // report info
00328 
00329     inline
00330     void report_info (double, int, double const* Control, double const* Info) {
00331       umfpack_di_report_info (Control, Info); 
00332     }
00333     inline
00334     void report_info (traits::complex_d const&, int, 
00335                       double const* Control, double const* Info) 
00336     {
00337       umfpack_zi_report_info (Control, Info); 
00338     }
00339 
00340 
00341     // report matrix 
00342 
00343     inline 
00344     int report_matrix (int n_row, int n_col, 
00345                        int const* Ap, int const* Ai, double const* Ax,
00346                        int col_form, double const* Control) 
00347     {
00348       return umfpack_di_report_matrix (n_row, n_col, Ap, Ai, Ax, 
00349                                        col_form, Control); 
00350     }
00351 
00352     inline 
00353     int report_matrix (int n_row, int n_col, 
00354                        int const* Ap, int const* Ai, 
00355                        traits::complex_d const* Ax, 
00356                        int col_form, double const* Control) 
00357     {
00358       int nnz = Ap[n_col];
00359       traits::detail::array<double> Axr (nnz); 
00360       if (!Axr.valid()) return UMFPACK_ERROR_out_of_memory;
00361       traits::detail::array<double> Axi (nnz); 
00362       if (!Axi.valid()) return UMFPACK_ERROR_out_of_memory;
00363       traits::detail::disentangle (Ax, Ax+nnz, 
00364                                    Axr.storage(), Axi.storage()); 
00365       return umfpack_zi_report_matrix (n_row, n_col, Ap, Ai, 
00366                                        Axr.storage(), Axi.storage(), 
00367                                        col_form, Control); 
00368     }
00369 
00370 
00371     // report triplet (coordinate)
00372 
00373     int report_triplet (int n_row, int n_col, int nz,
00374                         int const* Ti, int const* Tj, double const* Tx,
00375                         double const* Control) 
00376     {
00377       return umfpack_di_report_triplet (n_row, n_col, nz, Ti, Tj, Tx, Control);
00378     }
00379 
00380     int report_triplet (int n_row, int n_col, int nz,
00381                         int const* Ti, int const* Tj, 
00382                         traits::complex_d const* Tx, 
00383                         double const* Control) 
00384     {
00385       traits::detail::array<double> Txr (nz); 
00386       if (!Txr.valid()) return UMFPACK_ERROR_out_of_memory;
00387       traits::detail::array<double> Txi (nz); 
00388       if (!Txi.valid()) return UMFPACK_ERROR_out_of_memory;
00389       traits::detail::disentangle (Tx, Tx+nz, 
00390                                    Txr.storage(), Txi.storage()); 
00391       return umfpack_zi_report_triplet (n_row, n_col, nz, Ti, Tj, 
00392                                         Txr.storage(), Txi.storage(), 
00393                                         Control);
00394     }
00395 
00396 
00397     // report vector 
00398 
00399     int report_vector (int n, double const* X, double const* Control) {
00400       return umfpack_di_report_vector (n, X, Control);
00401     }
00402 
00403     int report_vector (int n, traits::complex_d const* X, 
00404                        double const* Control) 
00405     {
00406 #if 0
00407       // see UMFPACK v 4.1 User Guide
00408       traits::detail::array<double> Xr (n); 
00409       if (!Xr.valid()) return UMFPACK_ERROR_out_of_memory;
00410       traits::detail::array<double> Xi (n); 
00411       if (!Xi.valid()) return UMFPACK_ERROR_out_of_memory;
00412       traits::detail::disentangle (X, X+n, 
00413                                    Xr.storage(), Xi.storage()); 
00414       return umfpack_zi_report_vector (n, Xr.storage(), Xi.storage(), Control);
00415 #endif 
00416       return umfpack_zi_report_vector (n, 
00417                                        reinterpret_cast<double const*> (X), 
00418                                        reinterpret_cast<double const*> (0), 
00419                                        Control);
00420     }
00421 
00422 
00423     // report numeric 
00424 
00425     inline
00426     int report_numeric (double, int, void* Numeric, double const* Control) {
00427       return umfpack_di_report_numeric (Numeric, Control); 
00428     }
00429     inline
00430     int report_numeric (traits::complex_d const&, int, 
00431                         void* Numeric, double const* Control) 
00432     {
00433       return umfpack_zi_report_numeric (Numeric, Control); 
00434     }
00435 
00436 
00437     // report symbolic 
00438 
00439     inline
00440     int report_symbolic (double, int, void* Symbolic, double const* Control) {
00441       return umfpack_di_report_symbolic (Symbolic, Control); 
00442     }
00443     inline
00444     int report_symbolic (traits::complex_d const&, int, 
00445                          void* Symbolic, double const* Control) 
00446     {
00447       return umfpack_zi_report_symbolic (Symbolic, Control); 
00448     }
00449 
00450 
00451     // report permutation vector
00452 
00453     inline
00454     int report_perm (double, int, int np, 
00455                      int const* Perm, double const* Control) {
00456       return umfpack_di_report_perm (np, Perm, Control); 
00457     }
00458     inline
00459     int report_perm (traits::complex_d const&, int, int np, 
00460                      int const* Perm, double const* Control) {
00461       return umfpack_zi_report_perm (np, Perm, Control); 
00462     }
00463 
00464   }}
00465 
00466 }}}
00467 
00468 #endif // BOOST_NUMERIC_BINDINGS_UMFPACK_OVERLOADS_HPP