sysv.hpp

Go to the documentation of this file.
00001 /*
00002  * 
00003  * Copyright (c) Toon Knapen & 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  * KF acknowledges the support of the Faculty of Civil Engineering, 
00010  * University of Zagreb, Croatia.
00011  *
00012  */
00013 
00014 #ifndef BOOST_NUMERIC_BINDINGS_LAPACK_SYSV_HPP
00015 #define BOOST_NUMERIC_BINDINGS_LAPACK_SYSV_HPP
00016 
00017 #include <boost/numeric/bindings/traits/type_traits.hpp>
00018 #include <boost/numeric/bindings/traits/traits.hpp>
00019 #include <boost/numeric/bindings/lapack/lapack.h>
00020 #include <boost/numeric/bindings/lapack/ilaenv.hpp>
00021 #include <boost/numeric/bindings/traits/detail/array.hpp>
00022 
00023 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00024 #  include <boost/static_assert.hpp>
00025 #  include <boost/type_traits/is_same.hpp>
00026 #endif 
00027 
00028 #include <cassert>
00029 
00030 namespace boost { namespace numeric { namespace bindings { 
00031 
00032   namespace lapack {
00033 
00035     //
00036     // system of linear equations A * X = B with A symmetric matrix
00037     //
00039 
00040     namespace detail {
00041 
00042       inline 
00043       int sytrf_block (float, int const ispec, char const ul, int const n) {
00044         char ul2[2] = "x"; ul2[0] = ul; 
00045         return ilaenv (ispec, "SSYTRF", ul2, n); 
00046       }
00047       inline 
00048       int sytrf_block (double, int const ispec, char const ul, int const n) {
00049         char ul2[2] = "x"; ul2[0] = ul; 
00050         return ilaenv (ispec, "DSYTRF", ul2, n); 
00051       }
00052       inline 
00053       int sytrf_block (traits::complex_f, 
00054                        int const ispec, char const ul, int const n) 
00055       {
00056         char ul2[2] = "x"; ul2[0] = ul; 
00057         return ilaenv (ispec, "CSYTRF", ul2, n); 
00058       }
00059       inline 
00060       int sytrf_block (traits::complex_d, 
00061                        int const ispec, char const ul, int const n) 
00062       {
00063         char ul2[2] = "x"; ul2[0] = ul; 
00064         return ilaenv (ispec, "ZSYTRF", ul2, n); 
00065       }
00066 
00067     }
00068 
00069 
00070     template <typename SymmA>
00071     inline
00072     int sytrf_block (char const q, char const ul, SymmA const& a) {
00073 
00074 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
00075       BOOST_STATIC_ASSERT((boost::is_same<
00076         typename traits::matrix_traits<SymmA>::matrix_structure, 
00077         traits::general_t
00078       >::value));
00079 #endif
00080       assert (q == 'O' || q == 'M'); 
00081       assert (ul == 'U' || ul == 'L'); 
00082 
00083       int n = traits::matrix_size1 (a); 
00084       assert (n == traits::matrix_size2 (a)); 
00085 
00086 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00087       typedef typename traits::matrix_traits<SymmA>::value_type val_t; 
00088 #else 
00089       typedef typename SymmA::value_type val_t; 
00090 #endif 
00091       int ispec = (q == 'O' ? 1 : 2); 
00092       return detail::sytrf_block (val_t(), ispec, ul, n); 
00093     }
00094 
00095     template <typename SymmA>
00096     inline
00097     int sytrf_block (char const q, SymmA const& a) {
00098 
00099 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
00100       BOOST_STATIC_ASSERT((boost::is_same<
00101         typename traits::matrix_traits<SymmA>::matrix_structure, 
00102         traits::symmetric_t
00103       >::value));
00104 #endif
00105       assert (q == 'O' || q == 'M'); 
00106 
00107       char ul = traits::matrix_uplo_tag (a);
00108       int n = traits::matrix_size1 (a); 
00109       assert (n == traits::matrix_size2 (a)); 
00110 
00111 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00112       typedef typename traits::matrix_traits<SymmA>::value_type val_t; 
00113 #else 
00114       typedef typename SymmA::value_type val_t; 
00115 #endif 
00116       int ispec = (q == 'O' ? 1 : 2); 
00117       return detail::sytrf_block (val_t(), ispec, ul, n); 
00118     }
00119 
00120     template <typename SymmA>
00121     inline
00122     int sytrf_work (char const q, char const ul, SymmA const& a) {
00123 
00124 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
00125       BOOST_STATIC_ASSERT((boost::is_same<
00126         typename traits::matrix_traits<SymmA>::matrix_structure, 
00127         traits::general_t
00128       >::value));
00129 #endif
00130       assert (q == 'O' || q == 'M'); 
00131       assert (ul == 'U' || ul == 'L'); 
00132 
00133       int n = traits::matrix_size1 (a); 
00134       assert (n == traits::matrix_size2 (a)); 
00135 
00136 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00137       typedef typename traits::matrix_traits<SymmA>::value_type val_t; 
00138 #else 
00139       typedef typename SymmA::value_type val_t; 
00140 #endif 
00141       int lw = -13; 
00142       if (q == 'M') 
00143         lw = 1;
00144       if (q == 'O') 
00145         lw = n * detail::sytrf_block (val_t(), 1, ul, n); 
00146       return lw; 
00147     }
00148 
00149     template <typename SymmA>
00150     inline
00151     int sytrf_work (char const q, SymmA const& a) {
00152 
00153 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
00154       BOOST_STATIC_ASSERT((boost::is_same<
00155         typename traits::matrix_traits<SymmA>::matrix_structure, 
00156         traits::symmetric_t
00157       >::value));
00158 #endif
00159       assert (q == 'O' || q == 'M'); 
00160 
00161       char ul = traits::matrix_uplo_tag (a);
00162       int n = traits::matrix_size1 (a); 
00163       assert (n == traits::matrix_size2 (a)); 
00164 
00165 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00166       typedef typename traits::matrix_traits<SymmA>::value_type val_t; 
00167 #else 
00168       typedef typename SymmA::value_type val_t; 
00169 #endif 
00170       int lw = -13; 
00171       if (q == 'M') 
00172         lw = 1;
00173       if (q == 'O') 
00174         lw = n * detail::sytrf_block (val_t(), 1, ul, n); 
00175       return lw; 
00176     }
00177 
00178 
00179     template <typename SymmA>
00180     inline
00181     int sysv_work (char const q, char const ul, SymmA const& a) {
00182       return sytrf_work (q, ul, a); 
00183     }
00184 
00185     template <typename SymmA>
00186     inline
00187     int sysv_work (char const q, SymmA const& a) { return sytrf_work (q, a); }
00188 
00189 
00190     /*
00191      * sysv() computes the solution to a system of linear equations 
00192      * A * X = B, where A is an N-by-N symmetric matrix and X and B 
00193      * are N-by-NRHS matrices.
00194      *
00195      * The diagonal pivoting method is used to factor A as
00196      *   A = U * D * U^T,  if UPLO = 'U', 
00197      *   A = L * D * L^T,  if UPLO = 'L',
00198      * where  U (or L) is a product of permutation and unit upper 
00199      * (lower) triangular matrices, and D is symmetric and block 
00200      * diagonal with 1-by-1 and 2-by-2 diagonal blocks.  
00201      * The factored form of A is then used to solve the system 
00202      * of equations A * X = B.
00203      */
00204 
00205     namespace detail {
00206 
00207       inline 
00208       void sysv (char const uplo, int const n, int const nrhs,
00209                  float* a, int const lda, int* ipiv, 
00210                  float* b, int const ldb, 
00211                  float* w, int const lw, int* info) 
00212       {
00213         LAPACK_SSYSV (&uplo, &n, &nrhs, a, &lda, ipiv, b, &ldb, w, &lw, info);
00214       }
00215 
00216       inline 
00217       void sysv (char const uplo, int const n, int const nrhs,
00218                  double* a, int const lda, int* ipiv, 
00219                  double* b, int const ldb, 
00220                  double* w, int const lw, int* info) 
00221       {
00222         LAPACK_DSYSV (&uplo, &n, &nrhs, a, &lda, ipiv, b, &ldb, w, &lw, info);
00223       }
00224 
00225       inline 
00226       void sysv (char const uplo, int const n, int const nrhs,
00227                  traits::complex_f* a, int const lda, int* ipiv,  
00228                  traits::complex_f* b, int const ldb, 
00229                  traits::complex_f* w, int const lw, int* info) 
00230       {
00231         LAPACK_CSYSV (&uplo, &n, &nrhs, 
00232                       traits::complex_ptr (a), &lda, ipiv, 
00233                       traits::complex_ptr (b), &ldb, 
00234                       traits::complex_ptr (w), &lw, info);
00235       }
00236 
00237       inline 
00238       void sysv (char const uplo, int const n, int const nrhs,
00239                  traits::complex_d* a, int const lda, int* ipiv, 
00240                  traits::complex_d* b, int const ldb, 
00241                  traits::complex_d* w, int const lw, int* info) 
00242       {
00243         LAPACK_ZSYSV (&uplo, &n, &nrhs, 
00244                       traits::complex_ptr (a), &lda, ipiv, 
00245                       traits::complex_ptr (b), &ldb, 
00246                       traits::complex_ptr (w), &lw, info);
00247       }
00248 
00249       template <typename SymmA, typename MatrB, typename IVec, typename Work>
00250       inline
00251       int sysv (char const ul, SymmA& a, IVec& i, MatrB& b, Work& w) {
00252 
00253         int const n = traits::matrix_size1 (a);
00254         assert (n == traits::matrix_size2 (a)); 
00255         assert (n == traits::matrix_size1 (b)); 
00256         assert (n == traits::vector_size (i)); 
00257 
00258         int info; 
00259         sysv (ul, n, traits::matrix_size2 (b), 
00260               traits::matrix_storage (a), 
00261               traits::leading_dimension (a),
00262               traits::vector_storage (i),  
00263               traits::matrix_storage (b),
00264               traits::leading_dimension (b),
00265               traits::vector_storage (w), 
00266               traits::vector_size (w), 
00267               &info);
00268         return info; 
00269       }
00270 
00271     }
00272 
00273     template <typename SymmA, typename MatrB, typename IVec, typename Work>
00274     inline
00275     int sysv (char const ul, SymmA& a, IVec& i, MatrB& b, Work& w) {
00276 
00277       assert (ul == 'U' || ul == 'L'); 
00278 
00279 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
00280       BOOST_STATIC_ASSERT((boost::is_same<
00281         typename traits::matrix_traits<SymmA>::matrix_structure, 
00282         traits::general_t
00283       >::value));
00284       BOOST_STATIC_ASSERT((boost::is_same<
00285         typename traits::matrix_traits<MatrB>::matrix_structure, 
00286         traits::general_t
00287       >::value));
00288 #endif
00289 
00290       assert (traits::vector_size (w) >= 1); 
00291       return detail::sysv (ul, a, i, b, w); 
00292     }
00293 
00294     template <typename SymmA, typename MatrB, typename IVec, typename Work>
00295     inline
00296     int sysv (SymmA& a, IVec& i, MatrB& b, Work& w) {
00297 
00298 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
00299       BOOST_STATIC_ASSERT((boost::is_same<
00300         typename traits::matrix_traits<SymmA>::matrix_structure, 
00301         traits::symmetric_t
00302       >::value));
00303       BOOST_STATIC_ASSERT((boost::is_same<
00304         typename traits::matrix_traits<MatrB>::matrix_structure, 
00305         traits::general_t
00306       >::value));
00307 #endif
00308 
00309       assert (traits::vector_size (w) >= 1); 
00310       char uplo = traits::matrix_uplo_tag (a);
00311       return detail::sysv (uplo, a, i, b, w); 
00312     }
00313 
00314     template <typename SymmA, typename MatrB>
00315     inline
00316     int sysv (char const ul, SymmA& a, MatrB& b) {
00317       // with 'internal' pivot and work vectors 
00318 
00319       assert (ul == 'U' || ul == 'L'); 
00320 
00321 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
00322       BOOST_STATIC_ASSERT((boost::is_same<
00323         typename traits::matrix_traits<SymmA>::matrix_structure, 
00324         traits::general_t
00325       >::value));
00326       BOOST_STATIC_ASSERT((boost::is_same<
00327         typename traits::matrix_traits<MatrB>::matrix_structure, 
00328         traits::general_t
00329       >::value));
00330 #endif
00331 
00332       int const n = traits::matrix_size1 (a); 
00333       int info = -101; 
00334       traits::detail::array<int> i (n); 
00335 
00336       if (i.valid()) {
00337         info = -102; 
00338         int lw = sytrf_work ('O', ul, a); 
00339         assert (lw >= 1); // paranoia ? 
00340 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00341         typedef typename traits::matrix_traits<SymmA>::value_type val_t; 
00342 #else 
00343         typedef typename SymmA::value_type val_t; 
00344 #endif 
00345         traits::detail::array<val_t> w (lw); 
00346         if (w.valid()) 
00347           info =  detail::sysv (ul, a, i, b, w); 
00348       }
00349       return info; 
00350     }
00351 
00352     template <typename SymmA, typename MatrB>
00353     inline
00354     int sysv (SymmA& a, MatrB& b) {
00355       // with 'internal' pivot and work vectors 
00356 
00357 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
00358       BOOST_STATIC_ASSERT((boost::is_same<
00359         typename traits::matrix_traits<SymmA>::matrix_structure, 
00360         traits::symmetric_t
00361       >::value));
00362       BOOST_STATIC_ASSERT((boost::is_same<
00363         typename traits::matrix_traits<MatrB>::matrix_structure, 
00364         traits::general_t
00365       >::value));
00366 #endif
00367 
00368       int const n = traits::matrix_size1 (a); 
00369       char uplo = traits::matrix_uplo_tag (a);
00370       int info = -101; 
00371       traits::detail::array<int> i (n); 
00372 
00373       if (i.valid()) {
00374         info = -102; 
00375         int lw = sytrf_work ('O', a); 
00376         assert (lw >= 1); // paranoia ? 
00377 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00378         typedef typename traits::matrix_traits<SymmA>::value_type val_t; 
00379 #else 
00380         typedef typename SymmA::value_type val_t; 
00381 #endif 
00382         traits::detail::array<val_t> w (lw); 
00383         if (w.valid()) 
00384           info =  detail::sysv (uplo, a, i, b, w); 
00385       }
00386       return info; 
00387     }
00388 
00389 
00390     /*
00391      * sytrf() computes the factorization of a symmetric matrix A using
00392      * the  Bunch-Kaufman diagonal pivoting method. The form of the 
00393      * factorization is
00394      *    A = U * D * U^T  or  A = L * D * L^T
00395      * where U (or L) is a product of permutation and unit upper (lower)  
00396      * triangular matrices, and D is symmetric and block diagonal with 
00397      * 1-by-1 and 2-by-2 diagonal blocks.
00398      */
00399 
00400     namespace detail {
00401 
00402       inline 
00403       void sytrf (char const uplo, int const n, 
00404                   float* a, int const lda, int* ipiv, 
00405                   float* w, int const lw, int* info) 
00406       {
00407         LAPACK_SSYTRF (&uplo, &n, a, &lda, ipiv, w, &lw, info);
00408       }
00409 
00410       inline 
00411       void sytrf (char const uplo, int const n, 
00412                   double* a, int const lda, int* ipiv, 
00413                   double* w, int const lw, int* info) 
00414       {
00415         LAPACK_DSYTRF (&uplo, &n, a, &lda, ipiv, w, &lw, info);
00416       }
00417 
00418       inline 
00419       void sytrf (char const uplo, int const n, 
00420                   traits::complex_f* a, int const lda, int* ipiv,  
00421                   traits::complex_f* w, int const lw, int* info) 
00422       {
00423         LAPACK_CSYTRF (&uplo, &n, 
00424                        traits::complex_ptr (a), &lda, ipiv, 
00425                        traits::complex_ptr (w), &lw, info);
00426       }
00427 
00428       inline 
00429       void sytrf (char const uplo, int const n, 
00430                   traits::complex_d* a, int const lda, int* ipiv, 
00431                   traits::complex_d* w, int const lw, int* info) 
00432       {
00433         LAPACK_ZSYTRF (&uplo, &n, 
00434                        traits::complex_ptr (a), &lda, ipiv, 
00435                        traits::complex_ptr (w), &lw, info);
00436       }
00437 
00438       template <typename SymmA, typename IVec, typename Work>
00439       inline
00440       int sytrf (char const ul, SymmA& a, IVec& i, Work& w) {
00441 
00442         int const n = traits::matrix_size1 (a);
00443         assert (n == traits::matrix_size2 (a)); 
00444         assert (n == traits::vector_size (i)); 
00445 
00446         int info; 
00447         sytrf (ul, n, traits::matrix_storage (a), 
00448                traits::leading_dimension (a),
00449                traits::vector_storage (i),  
00450                traits::vector_storage (w), 
00451                traits::vector_size (w), 
00452                &info);
00453         return info; 
00454       }
00455 
00456     }
00457 
00458     template <typename SymmA, typename IVec, typename Work>
00459     inline
00460     int sytrf (char const ul, SymmA& a, IVec& i, Work& w) {
00461 
00462       assert (ul == 'U' || ul == 'L'); 
00463 
00464 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
00465       BOOST_STATIC_ASSERT((boost::is_same<
00466         typename traits::matrix_traits<SymmA>::matrix_structure, 
00467         traits::general_t
00468       >::value));
00469 #endif
00470 
00471       assert (traits::vector_size (w) >= 1); 
00472       return detail::sytrf (ul, a, i, w); 
00473     }
00474 
00475     template <typename SymmA, typename IVec, typename Work>
00476     inline
00477     int sytrf (SymmA& a, IVec& i, Work& w) {
00478 
00479 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
00480       BOOST_STATIC_ASSERT((boost::is_same<
00481         typename traits::matrix_traits<SymmA>::matrix_structure, 
00482         traits::symmetric_t
00483       >::value));
00484 #endif
00485 
00486       assert (traits::vector_size (w) >= 1); 
00487       char uplo = traits::matrix_uplo_tag (a);
00488       return detail::sytrf (uplo, a, i, w); 
00489     }
00490 
00491     template <typename SymmA, typename Ivec>
00492     inline
00493     int sytrf (char const ul, SymmA& a, Ivec& i) {
00494       // with 'internal' work vector
00495 
00496       assert (ul == 'U' || ul == 'L'); 
00497 
00498 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
00499       BOOST_STATIC_ASSERT((boost::is_same<
00500         typename traits::matrix_traits<SymmA>::matrix_structure, 
00501         traits::general_t
00502       >::value));
00503 #endif
00504 
00505       int info = -101; 
00506       int lw = sytrf_work ('O', ul, a); 
00507       assert (lw >= 1); // paranoia ? 
00508 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00509       typedef typename traits::matrix_traits<SymmA>::value_type val_t; 
00510 #else 
00511       typedef typename SymmA::value_type val_t; 
00512 #endif 
00513       traits::detail::array<val_t> w (lw); 
00514       if (w.valid()) 
00515         info =  detail::sytrf (ul, a, i, w); 
00516       return info; 
00517     }
00518 
00519     template <typename SymmA, typename Ivec>
00520     inline
00521     int sytrf (SymmA& a, Ivec& i) {
00522       // with 'internal' work vector 
00523 
00524 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
00525       BOOST_STATIC_ASSERT((boost::is_same<
00526         typename traits::matrix_traits<SymmA>::matrix_structure, 
00527         traits::symmetric_t
00528       >::value));
00529 #endif
00530 
00531       char uplo = traits::matrix_uplo_tag (a);
00532       int info = -101; 
00533       int lw = sytrf_work ('O', a); 
00534       assert (lw >= 1); // paranoia ? 
00535 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00536       typedef typename traits::matrix_traits<SymmA>::value_type val_t; 
00537 #else 
00538       typedef typename SymmA::value_type val_t; 
00539 #endif 
00540       traits::detail::array<val_t> w (lw); 
00541       if (w.valid()) 
00542         info =  detail::sytrf (uplo, a, i, w); 
00543       return info; 
00544     }
00545 
00546 
00547     /*
00548      * sytrs() solves a system of linear equations A*X = B with 
00549      * a symmetric matrix A using the factorization 
00550      *    A = U * D * U^T   or  A = L * D * L^T
00551      * computed by sytrf().
00552      */
00553 
00554     namespace detail {
00555 
00556       inline 
00557       void sytrs (char const uplo, int const n, int const nrhs,
00558                   float const* a, int const lda, int const* ipiv, 
00559                   float* b, int const ldb, int* info) 
00560       {
00561         LAPACK_SSYTRS (&uplo, &n, &nrhs, a, &lda, ipiv, b, &ldb, info);
00562       }
00563 
00564       inline 
00565       void sytrs (char const uplo, int const n, int const nrhs,
00566                   double const* a, int const lda, int const* ipiv, 
00567                   double* b, int const ldb, int* info) 
00568       {
00569         LAPACK_DSYTRS (&uplo, &n, &nrhs, a, &lda, ipiv, b, &ldb, info);
00570       }
00571 
00572       inline 
00573       void sytrs (char const uplo, int const n, int const nrhs,
00574                   traits::complex_f const* a, int const lda, 
00575                   int const* ipiv,  
00576                   traits::complex_f* b, int const ldb, int* info) 
00577       {
00578         LAPACK_CSYTRS (&uplo, &n, &nrhs, 
00579                       traits::complex_ptr (a), &lda, ipiv, 
00580                       traits::complex_ptr (b), &ldb, info);
00581       }
00582 
00583       inline 
00584       void sytrs (char const uplo, int const n, int const nrhs,
00585                   traits::complex_d const* a, int const lda, 
00586                   int const* ipiv, 
00587                   traits::complex_d* b, int const ldb, int* info) 
00588       {
00589         LAPACK_ZSYTRS (&uplo, &n, &nrhs, 
00590                        traits::complex_ptr (a), &lda, ipiv, 
00591                        traits::complex_ptr (b), &ldb, info);
00592       }
00593 
00594       template <typename SymmA, typename MatrB, typename IVec>
00595       inline
00596       int sytrs (char const ul, SymmA const& a, IVec const& i, MatrB& b) {
00597 
00598         int const n = traits::matrix_size1 (a);
00599         assert (n == traits::matrix_size2 (a)); 
00600         assert (n == traits::matrix_size1 (b)); 
00601         assert (n == traits::vector_size (i)); 
00602 
00603         int info; 
00604         sytrs (ul, n, traits::matrix_size2 (b), 
00605 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00606                traits::matrix_storage (a), 
00607 #else
00608                traits::matrix_storage_const (a), 
00609 #endif 
00610                traits::leading_dimension (a),
00611 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00612                traits::vector_storage (i),  
00613 #else
00614                traits::vector_storage_const (i),  
00615 #endif
00616                traits::matrix_storage (b),
00617                traits::leading_dimension (b), &info);
00618         return info; 
00619       }
00620 
00621     }
00622 
00623     template <typename SymmA, typename MatrB, typename IVec>
00624     inline
00625     int sytrs (char const ul, SymmA const& a, IVec const& i, MatrB& b) {
00626 
00627       assert (ul == 'U' || ul == 'L'); 
00628 
00629 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
00630       BOOST_STATIC_ASSERT((boost::is_same<
00631         typename traits::matrix_traits<SymmA>::matrix_structure, 
00632         traits::general_t
00633       >::value));
00634       BOOST_STATIC_ASSERT((boost::is_same<
00635         typename traits::matrix_traits<MatrB>::matrix_structure, 
00636         traits::general_t
00637       >::value));
00638 #endif
00639 
00640       return detail::sytrs (ul, a, i, b); 
00641     }
00642 
00643     template <typename SymmA, typename MatrB, typename IVec>
00644     inline
00645     int sytrs (SymmA const& a, IVec const& i, MatrB& b) {
00646 
00647 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
00648       BOOST_STATIC_ASSERT((boost::is_same<
00649         typename traits::matrix_traits<SymmA>::matrix_structure, 
00650         traits::symmetric_t
00651       >::value));
00652       BOOST_STATIC_ASSERT((boost::is_same<
00653         typename traits::matrix_traits<MatrB>::matrix_structure, 
00654         traits::general_t
00655       >::value));
00656 #endif
00657 
00658       char uplo = traits::matrix_uplo_tag (a);
00659       return detail::sytrs (uplo, a, i, b); 
00660     }
00661 
00662 
00663     // TO DO: sytri 
00664 
00665   }
00666 
00667 }}}
00668 
00669 #endif