MOOCHO (Single Doxygen Collection) Version of the Day
DenseLinAlgPack_InvCholUpdate.cpp
Go to the documentation of this file.
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <math.h>
00043 
00044 #include "DenseLinAlgPack_InvCholUpdate.hpp"
00045 #include "DenseLinAlgPack_DMatrixClass.hpp"
00046 #include "DenseLinAlgPack_DVectorOp.hpp"
00047 
00048 namespace {
00049 // sign function
00050 inline int sign(double v) {
00051   if(v > 0.0)
00052     return 1;
00053   else if(v == 0.0)
00054     return 0;
00055   else
00056     return -1;
00057 }
00058 
00059 }
00060 
00061 void DenseLinAlgPack::update_chol_factor(DMatrixSlice* pM, DVectorSlice* pu
00062   , const DVectorSlice& v)
00063 {
00064   DMatrixSlice  &M = *pM;
00065   DVectorSlice    &u = *pu;
00066 
00067   assert_gms_square(M);
00068   if(M.rows() != u.dim() || u.dim() != v.dim())
00069     throw std::length_error("update_chol_factor(): Sizes of matrix and vectors do not match");
00070 
00071   // Algorithm A3.4.1a in Dennis and Schabel
00072   
00073   // 0. Set lower diagonal to zero
00074   M.diag(-1) = 0.0;
00075 
00076   // 1,2 Find the largest k such that u(k) != 0.0
00077   size_type n = M.rows(), k = n;
00078   {
00079     DVectorSlice::reverse_iterator r_itr_u = u.rbegin();
00080     while(*r_itr_u == 0.0 && k > 1) --k;
00081   }
00082   
00083   // 3.
00084   {
00085     DVectorSlice::reverse_iterator
00086       r_itr_u_i   = u(1,k-1).rbegin(),  // iterator for u(i), i = k-1,...,1
00087       r_itr_u_ip1   = u(2,k).rbegin();    // iterator for u(i), i = k,...,2
00088     for(size_type i = k-1; i > 0 ; ++r_itr_u_i, ++r_itr_u_ip1, --i) {
00089       value_type u_i = *r_itr_u_i, u_ip1 = *r_itr_u_ip1;
00090       jacobi_rotate(&M, i, u_i, - u_ip1); // 3.1
00091       *r_itr_u_i = (u_i == 0.0) ? ::fabs(u_ip1) : ::sqrt(u_i * u_i + u_ip1 * u_ip1); // 3.2
00092     }
00093   }
00094 
00095   // 4. M.row(1) += u(1) * v 
00096   DenseLinAlgPack::Vp_StV(&M.row(1), u(1), v);
00097 
00098   // 5.
00099   {
00100     DVectorSlice::const_iterator
00101       itr_M_i_i = M.diag().begin(),   // iterator for M(i,i), for i = 1,...,k-1
00102       itr_M_ip1_i = M.diag(-1).begin(); // iterator for M(i+1,i), for i = 1,...,k-1
00103     for(size_type i = 1; i < k; ++i)
00104       jacobi_rotate(&M, i, *itr_M_i_i++, - *itr_M_ip1_i++);
00105   }
00106 
00107 /*  This code does not work
00108 
00109   size_type n = M.rows();
00110 
00111   // 0.
00112   {
00113     for(size_type i = 2; i <= n; ++i)
00114       M(i,i-1) = 0;
00115   }
00116 
00117   // 1.
00118   size_type k = n;
00119   
00120   // 2.
00121   while(u(k) == 0.0 && k > 1) --k;
00122 
00123   // 3.
00124   {
00125     for(size_type i = k-1; i >= 1; --i) {
00126       jacobi_rotate(M, i, u(i), - u(i+1));
00127       if(u(i) == 0)
00128         u(i) = ::fabs(u(i+1));
00129       else
00130         u(i) = ::sqrt(u(i) * u(i) + u(i+1) * u(i+1));
00131     }
00132   }
00133 
00134   // 4.
00135   {
00136     for(size_type j = 1; j <= n; ++j)
00137       M(1,j) = M(1,j) + dot(u(1),v(j));
00138   }
00139 
00140   // 5.
00141   {
00142     for(size_type i = 1; i <= k-1; ++i)
00143       jacobi_rotate(M, i, M(i,i), - M(i+1,i));
00144   }
00145 
00146 
00147 
00148 */
00149 
00150 
00151 }
00152 
00153 void DenseLinAlgPack::jacobi_rotate(DMatrixSlice* pM, size_type row_i, value_type alpha
00154   , value_type beta)
00155 {
00156   DMatrixSlice  &M = *pM;
00157 
00158   assert_gms_square(M);
00159     
00160   // Algorithm A3.4.1a in Dennis and Schabel
00161   
00162   // 1.
00163   value_type c, s;
00164   if(alpha == 0.0) {
00165     c = 0.0;
00166     s = sign(beta);
00167   }
00168   else {
00169     value_type den = ::sqrt(alpha * alpha + beta * beta);
00170     c = alpha / den;
00171     s = beta / den;
00172   }
00173 
00174   // 2.
00175   size_type i = row_i, n = M.rows();
00176 
00177   // Use iterators instead of element access
00178   DVectorSlice::iterator
00179     itr_M_i   = M.row(i)(i,n).begin(),  // iterator for M(i,j), for j = i,...,n
00180     itr_M_i_end = M.row(i)(i,n).end(),
00181     itr_M_ip1 = M.row(i+1)(i,n).begin();  // iterator for M(i+1,j), for j = i,...,n
00182 
00183   while(itr_M_i != itr_M_i_end) {
00184     value_type y = *itr_M_i, w = *itr_M_ip1;
00185     *itr_M_i++    = c * y - s * w;
00186     *itr_M_ip1++  = s * y + c * w;
00187   }
00188 
00189 /*  This code does not work
00190 
00191   size_type n = M.rows(), i = row_i;
00192 
00193   // 1.
00194   value_type c, s;
00195   if(alpha == 0) {
00196     c = 0.0;
00197     s = ::fabs(beta);
00198   }
00199   else {
00200     value_type den = ::sqrt(alpha*alpha + beta*beta);
00201     c = alpha / den;
00202     s = beta / den;
00203   }
00204 
00205   // 2.
00206   {
00207     for(size_type j = i; j <= n; ++j) {
00208       size_type y = M(i,j), w = M(i+1,j);
00209       M(i,j) = c*y - s*w;
00210       M(i+1,j) = s*y + c*w;
00211     }
00212   }
00213 
00214 
00215 
00216 */
00217 
00218 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines