DenseLinAlgPack_InvCholUpdate.cpp

00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <math.h>
00030 
00031 #include "DenseLinAlgPack_InvCholUpdate.hpp"
00032 #include "DenseLinAlgPack_DMatrixClass.hpp"
00033 #include "DenseLinAlgPack_DVectorOp.hpp"
00034 
00035 namespace {
00036 // sign function
00037 inline int sign(double v) {
00038   if(v > 0.0)
00039     return 1;
00040   else if(v == 0.0)
00041     return 0;
00042   else
00043     return -1;
00044 }
00045 
00046 }
00047 
00048 void DenseLinAlgPack::update_chol_factor(DMatrixSlice* pM, DVectorSlice* pu
00049   , const DVectorSlice& v)
00050 {
00051   DMatrixSlice  &M = *pM;
00052   DVectorSlice    &u = *pu;
00053 
00054   assert_gms_square(M);
00055   if(M.rows() != u.dim() || u.dim() != v.dim())
00056     throw std::length_error("update_chol_factor(): Sizes of matrix and vectors do not match");
00057 
00058   // Algorithm A3.4.1a in Dennis and Schabel
00059   
00060   // 0. Set lower diagonal to zero
00061   M.diag(-1) = 0.0;
00062 
00063   // 1,2 Find the largest k such that u(k) != 0.0
00064   size_type n = M.rows(), k = n;
00065   {
00066     DVectorSlice::reverse_iterator r_itr_u = u.rbegin();
00067     while(*r_itr_u == 0.0 && k > 1) --k;
00068   }
00069   
00070   // 3.
00071   {
00072     DVectorSlice::reverse_iterator
00073       r_itr_u_i   = u(1,k-1).rbegin(),  // iterator for u(i), i = k-1,...,1
00074       r_itr_u_ip1   = u(2,k).rbegin();    // iterator for u(i), i = k,...,2
00075     for(size_type i = k-1; i > 0 ; ++r_itr_u_i, ++r_itr_u_ip1, --i) {
00076       value_type u_i = *r_itr_u_i, u_ip1 = *r_itr_u_ip1;
00077       jacobi_rotate(&M, i, u_i, - u_ip1); // 3.1
00078       *r_itr_u_i = (u_i == 0.0) ? ::fabs(u_ip1) : ::sqrt(u_i * u_i + u_ip1 * u_ip1); // 3.2
00079     }
00080   }
00081 
00082   // 4. M.row(1) += u(1) * v 
00083   DenseLinAlgPack::Vp_StV(&M.row(1), u(1), v);
00084 
00085   // 5.
00086   {
00087     DVectorSlice::const_iterator
00088       itr_M_i_i = M.diag().begin(),   // iterator for M(i,i), for i = 1,...,k-1
00089       itr_M_ip1_i = M.diag(-1).begin(); // iterator for M(i+1,i), for i = 1,...,k-1
00090     for(size_type i = 1; i < k; ++i)
00091       jacobi_rotate(&M, i, *itr_M_i_i++, - *itr_M_ip1_i++);
00092   }
00093 
00094 /*  This code does not work
00095 
00096   size_type n = M.rows();
00097 
00098   // 0.
00099   {
00100     for(size_type i = 2; i <= n; ++i)
00101       M(i,i-1) = 0;
00102   }
00103 
00104   // 1.
00105   size_type k = n;
00106   
00107   // 2.
00108   while(u(k) == 0.0 && k > 1) --k;
00109 
00110   // 3.
00111   {
00112     for(size_type i = k-1; i >= 1; --i) {
00113       jacobi_rotate(M, i, u(i), - u(i+1));
00114       if(u(i) == 0)
00115         u(i) = ::fabs(u(i+1));
00116       else
00117         u(i) = ::sqrt(u(i) * u(i) + u(i+1) * u(i+1));
00118     }
00119   }
00120 
00121   // 4.
00122   {
00123     for(size_type j = 1; j <= n; ++j)
00124       M(1,j) = M(1,j) + dot(u(1),v(j));
00125   }
00126 
00127   // 5.
00128   {
00129     for(size_type i = 1; i <= k-1; ++i)
00130       jacobi_rotate(M, i, M(i,i), - M(i+1,i));
00131   }
00132 
00133 
00134 
00135 */
00136 
00137 
00138 }
00139 
00140 void DenseLinAlgPack::jacobi_rotate(DMatrixSlice* pM, size_type row_i, value_type alpha
00141   , value_type beta)
00142 {
00143   DMatrixSlice  &M = *pM;
00144 
00145   assert_gms_square(M);
00146     
00147   // Algorithm A3.4.1a in Dennis and Schabel
00148   
00149   // 1.
00150   value_type c, s;
00151   if(alpha == 0.0) {
00152     c = 0.0;
00153     s = sign(beta);
00154   }
00155   else {
00156     value_type den = ::sqrt(alpha * alpha + beta * beta);
00157     c = alpha / den;
00158     s = beta / den;
00159   }
00160 
00161   // 2.
00162   size_type i = row_i, n = M.rows();
00163 
00164   // Use iterators instead of element access
00165   DVectorSlice::iterator
00166     itr_M_i   = M.row(i)(i,n).begin(),  // iterator for M(i,j), for j = i,...,n
00167     itr_M_i_end = M.row(i)(i,n).end(),
00168     itr_M_ip1 = M.row(i+1)(i,n).begin();  // iterator for M(i+1,j), for j = i,...,n
00169 
00170   while(itr_M_i != itr_M_i_end) {
00171     value_type y = *itr_M_i, w = *itr_M_ip1;
00172     *itr_M_i++    = c * y - s * w;
00173     *itr_M_ip1++  = s * y + c * w;
00174   }
00175 
00176 /*  This code does not work
00177 
00178   size_type n = M.rows(), i = row_i;
00179 
00180   // 1.
00181   value_type c, s;
00182   if(alpha == 0) {
00183     c = 0.0;
00184     s = ::fabs(beta);
00185   }
00186   else {
00187     value_type den = ::sqrt(alpha*alpha + beta*beta);
00188     c = alpha / den;
00189     s = beta / den;
00190   }
00191 
00192   // 2.
00193   {
00194     for(size_type j = i; j <= n; ++j) {
00195       size_type y = M(i,j), w = M(i+1,j);
00196       M(i,j) = c*y - s*w;
00197       M(i+1,j) = s*y + c*w;
00198     }
00199   }
00200 
00201 
00202 
00203 */
00204 
00205 }

Generated on Tue Jul 13 09:28:46 2010 for DenseLinAlgPack: Concreate C++ Classes for Dense Blas-Compatible Linear Algebra by  doxygen 1.4.7