SVBackSubstitution.hpp

Go to the documentation of this file.

/* 
 * SVBackSubstitution.hpp 
 *
 * This file is part of the Chemical Data Processing Toolkit
 *
 * Copyright (C) 2003 Thomas Seidel <thomas.seidel@univie.ac.at>
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this library; see the file COPYING. If not, write to
 * the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */
 
#ifndef CDPL_MATH_SVBACKSUBSTITUTION_HPP
#define CDPL_MATH_SVBACKSUBSTITUTION_HPP
 
#include "CDPL/Math/Vector.hpp"
#include "CDPL/Math/Matrix.hpp"
 
 
namespace CDPL
{
 
    namespace Math
    {
 
        template <typename M1, typename V1, typename M2, typename V2, typename V3>
        void svBackSubstitution(const M1& u, const V1& w, const M2& v, const V2& b, V3& x);
    } // namespace Math
} // namespace CDPL
 
 
// Implementation
 
template <typename M1, typename V1, typename M2, typename V2, typename V3>
void CDPL::Math::svBackSubstitution(const M1& u, const V1& w, const M2& v, const V2& b, V3& x)
{
    typedef typename V3::value_type T;
 
    typename M1::size_type n = u.size2();
 
    Vector<T> tmp(n);
 
    for (typename M1::size_type j = 0; j < n; j++) { // Calculate trans(U) * B
        T s = 0;
 
        if (w(j) != 0) // Nonzero result only if w(j) is nonzero
            s = inner_prod(column(u, j), b) / w(j);
 
        tmp(j) = s;
    }
 
    x = prod(v, tmp);
}
 
#endif // CDPL_MATH_SVBACKSUBSTITUTION_HPP