DenseLinAlgPack_InvCholUpdate.cpp

// @HEADER
// ***********************************************************************
// 
// Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
//                  Copyright (2003) Sandia Corporation
// 
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
// 
// 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.1 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; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
// 
// ***********************************************************************
// @HEADER

#include <math.h>

#include "DenseLinAlgPack_InvCholUpdate.hpp"
#include "DenseLinAlgPack_DMatrixClass.hpp"
#include "DenseLinAlgPack_DVectorOp.hpp"

namespace {
// sign function
inline int sign(double v) {
  if(v > 0.0)
    return 1;
  else if(v == 0.0)
    return 0;
  else
    return -1;
}

}

void DenseLinAlgPack::update_chol_factor(DMatrixSlice* pM, DVectorSlice* pu
  , const DVectorSlice& v)
{
  DMatrixSlice  &M = *pM;
  DVectorSlice    &u = *pu;

  assert_gms_square(M);
  if(M.rows() != u.dim() || u.dim() != v.dim())
    throw std::length_error("update_chol_factor(): Sizes of matrix and vectors do not match");

  // Algorithm A3.4.1a in Dennis and Schabel
  
  // 0. Set lower diagonal to zero
  M.diag(-1) = 0.0;

  // 1,2 Find the largest k such that u(k) != 0.0
  size_type n = M.rows(), k = n;
  {
    DVectorSlice::reverse_iterator r_itr_u = u.rbegin();
    while(*r_itr_u == 0.0 && k > 1) --k;
  }
  
  // 3.
  {
    DVectorSlice::reverse_iterator
      r_itr_u_i   = u(1,k-1).rbegin(),  // iterator for u(i), i = k-1,...,1
      r_itr_u_ip1   = u(2,k).rbegin();    // iterator for u(i), i = k,...,2
    for(size_type i = k-1; i > 0 ; ++r_itr_u_i, ++r_itr_u_ip1, --i) {
      value_type u_i = *r_itr_u_i, u_ip1 = *r_itr_u_ip1;
      jacobi_rotate(&M, i, u_i, - u_ip1); // 3.1
      *r_itr_u_i = (u_i == 0.0) ? ::fabs(u_ip1) : ::sqrt(u_i * u_i + u_ip1 * u_ip1); // 3.2
    }
  }

  // 4. M.row(1) += u(1) * v 
  DenseLinAlgPack::Vp_StV(&M.row(1), u(1), v);

  // 5.
  {
    DVectorSlice::const_iterator
      itr_M_i_i = M.diag().begin(),   // iterator for M(i,i), for i = 1,...,k-1
      itr_M_ip1_i = M.diag(-1).begin(); // iterator for M(i+1,i), for i = 1,...,k-1
    for(size_type i = 1; i < k; ++i)
      jacobi_rotate(&M, i, *itr_M_i_i++, - *itr_M_ip1_i++);
  }

/*  This code does not work

  size_type n = M.rows();

  // 0.
  {
    for(size_type i = 2; i <= n; ++i)
      M(i,i-1) = 0;
  }

  // 1.
  size_type k = n;
  
  // 2.
  while(u(k) == 0.0 && k > 1) --k;

  // 3.
  {
    for(size_type i = k-1; i >= 1; --i) {
      jacobi_rotate(M, i, u(i), - u(i+1));
      if(u(i) == 0)
        u(i) = ::fabs(u(i+1));
      else
        u(i) = ::sqrt(u(i) * u(i) + u(i+1) * u(i+1));
    }
  }

  // 4.
  {
    for(size_type j = 1; j <= n; ++j)
      M(1,j) = M(1,j) + dot(u(1),v(j));
  }

  // 5.
  {
    for(size_type i = 1; i <= k-1; ++i)
      jacobi_rotate(M, i, M(i,i), - M(i+1,i));
  }



*/


}

void DenseLinAlgPack::jacobi_rotate(DMatrixSlice* pM, size_type row_i, value_type alpha
  , value_type beta)
{
  DMatrixSlice  &M = *pM;

  assert_gms_square(M);
    
  // Algorithm A3.4.1a in Dennis and Schabel
  
  // 1.
  value_type c, s;
  if(alpha == 0.0) {
    c = 0.0;
    s = sign(beta);
  }
  else {
    value_type den = ::sqrt(alpha * alpha + beta * beta);
    c = alpha / den;
    s = beta / den;
  }

  // 2.
  size_type i = row_i, n = M.rows();

  // Use iterators instead of element access
  DVectorSlice::iterator
    itr_M_i   = M.row(i)(i,n).begin(),  // iterator for M(i,j), for j = i,...,n
    itr_M_i_end = M.row(i)(i,n).end(),
    itr_M_ip1 = M.row(i+1)(i,n).begin();  // iterator for M(i+1,j), for j = i,...,n

  while(itr_M_i != itr_M_i_end) {
    value_type y = *itr_M_i, w = *itr_M_ip1;
    *itr_M_i++    = c * y - s * w;
    *itr_M_ip1++  = s * y + c * w;
  }

/*  This code does not work

  size_type n = M.rows(), i = row_i;

  // 1.
  value_type c, s;
  if(alpha == 0) {
    c = 0.0;
    s = ::fabs(beta);
  }
  else {
    value_type den = ::sqrt(alpha*alpha + beta*beta);
    c = alpha / den;
    s = beta / den;
  }

  // 2.
  {
    for(size_type j = i; j <= n; ++j) {
      size_type y = M(i,j), w = M(i+1,j);
      M(i,j) = c*y - s*w;
      M(i+1,j) = s*y + c*w;
    }
  }



*/

}