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

Generated on Tue Jul 13 09:30:50 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7