include/ewalena/lac/solver_jacobi.h

Go to the documentation of this file.

// -------------------------------------------------------------------
// @author Toby D. Young
// @version $Id: solver_jacobi.h 988 2012-11-26 19:09:39Z oneliefleft $
//
// Copyright 2012 namespace ewalena authors. All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions
// are met:
//
// 1. Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above
//    copyright notice, this list of conditions and the following
//    disclaimer in the documentation and/or other materials provided
//    with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE NAMEPSACE EWALENA AUTHORS ``AS
// IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
// NAMESPACE EWALENA AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
// OF THE POSSIBILITY OF SUCH DAMAGE.
//
// The views and conclusions contained in the software and
// documentation are those of the authors and should not be
// interpreted as representing official policies, either expressed or
// implied, of the namespace ewalena authors.
// -------------------------------------------------------------------

#ifndef __ewalena_solver_jacobi_h
#define __ewalena_solver_jacobi_h

#ifdef EWALENA_HAVE_CONFIG_H
#include <ewalena/base/config.h>
#endif

#include <ewalena/lac/solver_base.h>

namespace ewalena 
{

  template <typename ValueType = double>
    class SolverJacobi 
    : 
    public SolverBase<ValueType>
    {
      
    public:
    
    SolverJacobi (SolverControl<ValueType> &solver_control);
    
    virtual void solve (const Matrix<ValueType> &A,
                        Vector<ValueType>       &lambda,
                        VectorBasis<ValueType>  &X);
    
    protected:
    
    void do_jacobi_rotation (ValueType       *index_i, 
                             ValueType       *index_j, 
                             const ValueType &root);
    
    void sort (Vector<ValueType> &lambda);
    
    void sort (Vector<ValueType>      &lambda,
               VectorBasis<ValueType> &X);

    }; /* SolverJacobi */

  /*-------------- Inline and Other Functions -----------------------*/

  template <typename ValueType>
    inline
    void
    SolverJacobi<ValueType>::do_jacobi_rotation (ValueType       *index_i, 
                                                 ValueType       *index_j, 
                                                 const ValueType &root)
    {
      const ValueType c     = ValueType(1) / sqrt(root*root+ValueType(1));
      const ValueType s     = root*c;
      const ValueType tau   = s / (ValueType(1)+c);
      const ValueType tmp_1 = *index_i;
      const ValueType tmp_2 = *index_j;
      
      *index_i -= s*(tmp_2 + tmp_1*tau);
      *index_j += s*(tmp_1 - tmp_2*tau);
    }

  template <typename ValueType>
    inline
    void
    SolverJacobi<ValueType>::sort (Vector<ValueType>      &lambda)
    {

      /* Track if swapping is taking place. */
      bool swapped = false;

      /* Loop over all eigenvalues. */
      for (unsigned int i=0; i<lambda.size (); ++i)
        {

          /* We haven't swapped anything yet. */
          swapped = false;

          for (unsigned int j=i+1; j<lambda.size (); ++j)
            {
              if (lambda(i)>lambda(j))
                {
                  /* Swap the eigenvalues. */
                  const ValueType tmp_value = lambda(i);
                  lambda(i) = lambda(j);
                  lambda(j) = tmp_value;

                  swapped = true;
                }
            }

          /* If no swapping is taking place then the algorithm is
             done. */
          if (!swapped)
            return;

        }

      /* If the eigenvalues are not in order by now, something went
         deeply wrong... */
      if (swapped)
        assert (false);
    }

  template <typename ValueType>
    inline
    void
    SolverJacobi<ValueType>::sort (Vector<ValueType>      &lambda,
                                   VectorBasis<ValueType> &X)
    {

      /* Track if swapping is taking place. */
      bool swapped = false;

      /* Loop over all eigenvalues. */
      for (unsigned int i=0; i<lambda.size (); ++i)
        {

          /* We haven't swapped anything yet. */
          swapped = false;

          for (unsigned int j=i+1; j<lambda.size (); ++j)
            {
              if (lambda(i)>lambda(j))
                {

                  /* Swap the eigenvalues. */
                  const ValueType tmp_value = lambda(i);
                  lambda(i) = lambda(j);
                  lambda(j) = tmp_value;

                  /* Swap the eigenvectors by column. */
                  Vector<ValueType> tmp_vector (X.n_rows ());

                  for (unsigned int k=0; k<X.n_rows (); ++k)
                    {
                      tmp_vector(k) = X(k,i);
                      X(k,i) = X(k,j);
                      X(k,j) = tmp_vector(k);
                    }

                  swapped = true;
                }
            }

          /* If no swapping is taking place then the algorithm is
             done. */
          if (!swapped)
            return;

        }

      /* If the eigenvalues are not in order by now, something went
         deeply wrong... */
      if (swapped)
        assert (false);
    }

} /* namespace ewalena */

#endif /* __ewalena_solver_jacobi_h */