AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects Version of the Day
AbstractLinAlgPack_rank_2_chol_update.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 // 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 "AbstractLinAlgPack_rank_2_chol_update.hpp"
00043 #include "DenseLinAlgPack_DMatrixAsTriSym.hpp"
00044 #include "DenseLinAlgPack_DVectorClass.hpp"
00045 #include "DenseLinAlgPack_DVectorOp.hpp"
00046 #include "DenseLinAlgPack_AssertOp.hpp"
00047 #include "DenseLinAlgPack_BLAS_Cpp.hpp"
00048 
00049 void AbstractLinAlgPack::rank_2_chol_update(
00050   const value_type     a
00051   ,DVectorSlice         *u
00052   ,const DVectorSlice   &v
00053   ,DVectorSlice         *w
00054   ,DMatrixSliceTriEle         *R
00055   ,BLAS_Cpp::Transp    R_trans
00056   )
00057 {
00058   using BLAS_Cpp::no_trans;
00059   using BLAS_Cpp::trans;
00060   using BLAS_Cpp::upper;
00061   using BLAS_Cpp::lower;
00062   using BLAS_Cpp::rotg;
00063   using DenseLinAlgPack::row;
00064   using DenseLinAlgPack::col;
00065   using DenseLinAlgPack::rot;
00066   using DenseLinAlgPack::Vp_StV;
00067   //
00068   // The idea for this routine is to perform a set of givens rotations to return
00069   // op(R) +a *u*v' back to triangular form.  The first set of (n-1) givens
00070   // rotations Q1 eleminates the elements in u to form ||u||2 * e(last_i) and transform:
00071   //
00072   //     Q1 * (op(R) + a*u*v') -> op(R1) + a*||u||2 * e(last_i) * v'
00073   //
00074   // where op(R1) is a upper or lower Hessenberg matrix.  The ordering of the rotations
00075   // and and wether last_i = 1 or n depends on whether op(R) is logically upper or lower
00076   // triangular.
00077   //
00078   // If op(R) logically lower triangular then we have:
00079   //
00080   // [ x          ]       [ x ]
00081   // [ x  x       ]       [ x ]
00082   // [ x  x  x    ]       [ x ]
00083   // [ x  x  x  x ] + a * [ x ] * [ x x x x ]
00084   // 
00085   //       op(R)    + a *   u   *  v'
00086   // 
00087   //  Rotations are applied to zero out u(i), for i = 1..n-1 to form
00088   //  (first_i = 1, di = +1, last_i = n):
00089   // 
00090   //  [ x  x       ]       [   ]
00091   //  [ x  x  x    ]       [   ]
00092   //  [ x  x  x  x ]       [   ]
00093   //  [ x  x  x  x ] + a * [ x ] * [ x x x x ]
00094   // 
00095   //      op(R1)     + a * ||u||2*e(n) * v'
00096   // 
00097   //  If op(R) logically upper triangular then we have:
00098   // 
00099   //  [ x  x  x  x ] + a * [ x ] * [ x x x x ]
00100   //  [    x  x  x ]       [ x ]
00101   //  [       x  x ]       [ x ]
00102   //  [          x ]       [ x ]
00103   // 
00104   //       op(R)     + a *   u   *  v'
00105   // 
00106   //  Rotations are applied to zero out u(i), for i = n..2 to form
00107   //  (first_i = n, di = -1, last_i = 1):
00108   // 
00109   //  [ x  x  x  x ] + a * [ x ] * [ x x x x ]
00110   //  [ x  x  x  x ]       [   ]
00111   //  [    x  x  x ]       [   ]
00112   //  [       x  x ]       [   ]
00113   // 
00114   //      op(R1)     + a * ||u||2*e(1) * v'
00115   // 
00116   //  The off diagonal elements in R created by this set of rotations is not stored
00117   //  in R(*,*) but instead in the workspace vector w(*).  More precisely,
00118   //
00119   //       w(i+wo) = R(i,i+di), for i = first_i...last_i-di
00120   //
00121   //       where wo = 0 if lower, wo = -1 if upper
00122   //
00123   //  This gives the client more flexibility in how workspace is handeled.
00124   //  Other implementations overwrite the off diagonal of R (lazyness I guess).
00125   // 
00126   //  The update a*||u||2*e(last_i)*v' is then added to the row op(R1)(last_i,:) to
00127   //  give op(R2) which is another upper or lower Hessenberg matrix:
00128   //  
00129   //       op(R1) + a*||u||2 * e(last_i) * v' -> op(R2)
00130   // 
00131   //  Then, we apply (n-1) more givens rotations to return op(R2) back to triangular
00132   //  form and we are finished:
00133   // 
00134   //       Q2 * op(R2) -> op(R_new)
00135   // 
00136   const size_type n = R->rows();
00137   DenseLinAlgPack::Mp_M_assert_sizes( n, n, no_trans, u->dim(), v.dim(), no_trans );
00138   // Is op(R) logically upper or lower triangular
00139   const BLAS_Cpp::Uplo
00140     opR_uplo = ( (R->uplo()==upper && R_trans==no_trans) || (R->uplo()==lower && R_trans==trans)
00141            ? upper : lower );
00142   // Get frame of reference
00143   size_type first_i, last_i;
00144   int di, wo;
00145   if( opR_uplo == lower ) {
00146     first_i =  1;
00147     di      = +1;
00148     last_i  =  n;
00149     wo      =  0;
00150   }
00151   else {
00152     first_i =  n;
00153     di      = -1;
00154     last_i  =  1;
00155     wo      = -1;
00156   }
00157   // Zero out off diagonal workspace
00158   if(n > 1)
00159     *w = 0.0;
00160   // Find the first u(k) != 0 in the direction of the forward sweeps
00161   size_type k = first_i;
00162   for( ; k != last_i+di; k+=di )
00163     if( (*u)(k) != 0.0 ) break;
00164   if( k == last_i+di ) return; // All of u is zero, no update to perform!
00165   //
00166   // Transform Q(.)*...*Q(.) * u -> ||u||2 * e(last_i) while applying
00167   // Q1*R forming R1 (upper or lower Hessenberg)
00168   //
00169   {for( size_type i = k; i != last_i; i+=di ) {
00170     // [  c  s ] * [ u(i+di) ] -> [ u(i+id) ]
00171     // [ -s  c ]   [ u(i)    ]    [ 0       ]
00172     value_type c, s;  // Here c and s are computed to zero out u(i)
00173     rotg( &(*u)(i+di), &(*u)(i), &c, &s );
00174     // [  c  s ] * [ op(R)(i+di,:) ] -> [ op(R)(i+di,:) ]
00175     // [ -s  c ]   [ op(R)(i   ,:) ]    [ op(R)(i   ,:) ]
00176     DVectorSlice
00177       opR_row_idi = row(R->gms(),R_trans,i+di),
00178       opR_row_i   = row(R->gms(),R_trans,i);
00179     if( opR_uplo == lower ) {
00180       // op(R)(i   ,:) = [ x  x  x  w(i+wo)     ] i
00181       // op(R)(i+di,:) = [ x  x  x   x          ]
00182       //                         i
00183       rot( c, s, &opR_row_idi(1  ,i  ), &opR_row_i(1,i)  ); // First nonzero columns
00184       rot( c, s, &opR_row_idi(i+1,i+1), &(*w)(i+wo,i+wo) ); // Last nonzero column
00185     }
00186     else {
00187       // op(R)(i+di,:) = [         x  x  x  x ]
00188       // op(R)(i   ,:) = [    w(i+wo) x  x  x ] i
00189       //                              i
00190       rot( c, s, &opR_row_idi(i-1,i-1), &(*w)(i+wo,i+wo) ); // First nonzero column
00191       rot( c, s, &opR_row_idi(i  ,n  ), &opR_row_i(i,n)  ); // Last nonzero columns
00192     }
00193   }}
00194   //
00195   // op(R1) + a * ||u||2 * e(last_i) * v' -> op(R2)
00196   //
00197   Vp_StV( &row(R->gms(),R_trans,last_i), a * (*u)(last_i), v );
00198   //
00199   // Apply (k-1) more givens rotations to return op(R2) back to triangular form:
00200   //
00201   // Q2 * R2 -> R_new
00202   //
00203   {for( size_type i = last_i; i != k; i-=di ) {
00204     DVectorSlice
00205       opR_row_i   = row(R->gms(),R_trans,i),
00206       opR_row_idi = row(R->gms(),R_trans,i-di);
00207     value_type
00208       &w_idiwo = (*w)(i-di+wo);
00209     // [  c  s ]   [ op(R)(i,i)    ]    [ op(R)(i,i) ]
00210     // [ -s  c ] * [ op(R)(i-id,i) ] -> [ 0          ]  w(i-di+wo)
00211     value_type &R_i_i = opR_row_i(i);
00212     value_type c, s;  // Here c and s are computed to zero out w(i-di+wo)
00213     rotg( &R_i_i, &w_idiwo, &c, &s );
00214     // Make sure the diagonal is positive
00215     value_type scale = +1.0;
00216     if( R_i_i < 0.0 ) {
00217       R_i_i    = -R_i_i;
00218       scale    = -1.0;
00219       w_idiwo  = -w_idiwo; // Must also change this stored value
00220     }
00221     // [  c  s ] * [ op(R)(i,:)   ] -> [ op(R)(i,:)    ]
00222     // [ -s  c ]   [ op(R)(i-id,:)]    [ op(R)(i-id,:) ]
00223     if( opR_uplo == lower ) {
00224       // op(R)(i-di,:) = [ x  x  x   0     ]      
00225       // op(R)(i,:)    = [ x  x  x   x     ] i
00226       //                             i
00227       rot( scale*c, scale*s, &opR_row_i(1,i-1), &opR_row_idi(1,i-1) );
00228     }
00229     else {
00230       // op(R)(i,:)     = [      x  x  x  x ] i
00231       // op(R)(i-di,:)  = [      0  x  x  x ]
00232       //                         i
00233       rot( scale*c, scale*s, &opR_row_i(i+1,n), &opR_row_idi(i+1,n) );
00234     }
00235   }}
00236   // When we get here note that u contains the first set of rotations Q1 and
00237   // w contains the second set of rotations Q2.
00238 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends