AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AbstractLinAlgPack_rank_2_chol_update.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include "AbstractLinAlgPack_rank_2_chol_update.hpp"
43 #include "DenseLinAlgPack_DMatrixAsTriSym.hpp"
44 #include "DenseLinAlgPack_DVectorClass.hpp"
45 #include "DenseLinAlgPack_DVectorOp.hpp"
46 #include "DenseLinAlgPack_AssertOp.hpp"
47 #include "DenseLinAlgPack_BLAS_Cpp.hpp"
48 
49 void AbstractLinAlgPack::rank_2_chol_update(
50  const value_type a
51  ,DVectorSlice *u
52  ,const DVectorSlice &v
53  ,DVectorSlice *w
54  ,DMatrixSliceTriEle *R
55  ,BLAS_Cpp::Transp R_trans
56  )
57 {
58  using BLAS_Cpp::no_trans;
59  using BLAS_Cpp::trans;
60  using BLAS_Cpp::upper;
61  using BLAS_Cpp::lower;
62  using BLAS_Cpp::rotg;
63  using DenseLinAlgPack::row;
64  using DenseLinAlgPack::col;
65  using DenseLinAlgPack::rot;
67  //
68  // The idea for this routine is to perform a set of givens rotations to return
69  // op(R) +a *u*v' back to triangular form. The first set of (n-1) givens
70  // rotations Q1 eleminates the elements in u to form ||u||2 * e(last_i) and transform:
71  //
72  // Q1 * (op(R) + a*u*v') -> op(R1) + a*||u||2 * e(last_i) * v'
73  //
74  // where op(R1) is a upper or lower Hessenberg matrix. The ordering of the rotations
75  // and and wether last_i = 1 or n depends on whether op(R) is logically upper or lower
76  // triangular.
77  //
78  // If op(R) logically lower triangular then we have:
79  //
80  // [ x ] [ x ]
81  // [ x x ] [ x ]
82  // [ x x x ] [ x ]
83  // [ x x x x ] + a * [ x ] * [ x x x x ]
84  //
85  // op(R) + a * u * v'
86  //
87  // Rotations are applied to zero out u(i), for i = 1..n-1 to form
88  // (first_i = 1, di = +1, last_i = n):
89  //
90  // [ x x ] [ ]
91  // [ x x x ] [ ]
92  // [ x x x x ] [ ]
93  // [ x x x x ] + a * [ x ] * [ x x x x ]
94  //
95  // op(R1) + a * ||u||2*e(n) * v'
96  //
97  // If op(R) logically upper triangular then we have:
98  //
99  // [ x x x x ] + a * [ x ] * [ x x x x ]
100  // [ x x x ] [ x ]
101  // [ x x ] [ x ]
102  // [ x ] [ x ]
103  //
104  // op(R) + a * u * v'
105  //
106  // Rotations are applied to zero out u(i), for i = n..2 to form
107  // (first_i = n, di = -1, last_i = 1):
108  //
109  // [ x x x x ] + a * [ x ] * [ x x x x ]
110  // [ x x x x ] [ ]
111  // [ x x x ] [ ]
112  // [ x x ] [ ]
113  //
114  // op(R1) + a * ||u||2*e(1) * v'
115  //
116  // The off diagonal elements in R created by this set of rotations is not stored
117  // in R(*,*) but instead in the workspace vector w(*). More precisely,
118  //
119  // w(i+wo) = R(i,i+di), for i = first_i...last_i-di
120  //
121  // where wo = 0 if lower, wo = -1 if upper
122  //
123  // This gives the client more flexibility in how workspace is handeled.
124  // Other implementations overwrite the off diagonal of R (lazyness I guess).
125  //
126  // The update a*||u||2*e(last_i)*v' is then added to the row op(R1)(last_i,:) to
127  // give op(R2) which is another upper or lower Hessenberg matrix:
128  //
129  // op(R1) + a*||u||2 * e(last_i) * v' -> op(R2)
130  //
131  // Then, we apply (n-1) more givens rotations to return op(R2) back to triangular
132  // form and we are finished:
133  //
134  // Q2 * op(R2) -> op(R_new)
135  //
136  const size_type n = R->rows();
137  DenseLinAlgPack::Mp_M_assert_sizes( n, n, no_trans, u->dim(), v.dim(), no_trans );
138  // Is op(R) logically upper or lower triangular
139  const BLAS_Cpp::Uplo
140  opR_uplo = ( (R->uplo()==upper && R_trans==no_trans) || (R->uplo()==lower && R_trans==trans)
141  ? upper : lower );
142  // Get frame of reference
143  size_type first_i, last_i;
144  int di, wo;
145  if( opR_uplo == lower ) {
146  first_i = 1;
147  di = +1;
148  last_i = n;
149  wo = 0;
150  }
151  else {
152  first_i = n;
153  di = -1;
154  last_i = 1;
155  wo = -1;
156  }
157  // Zero out off diagonal workspace
158  if(n > 1)
159  *w = 0.0;
160  // Find the first u(k) != 0 in the direction of the forward sweeps
161  size_type k = first_i;
162  for( ; k != last_i+di; k+=di )
163  if( (*u)(k) != 0.0 ) break;
164  if( k == last_i+di ) return; // All of u is zero, no update to perform!
165  //
166  // Transform Q(.)*...*Q(.) * u -> ||u||2 * e(last_i) while applying
167  // Q1*R forming R1 (upper or lower Hessenberg)
168  //
169  {for( size_type i = k; i != last_i; i+=di ) {
170  // [ c s ] * [ u(i+di) ] -> [ u(i+id) ]
171  // [ -s c ] [ u(i) ] [ 0 ]
172  value_type c, s; // Here c and s are computed to zero out u(i)
173  rotg( &(*u)(i+di), &(*u)(i), &c, &s );
174  // [ c s ] * [ op(R)(i+di,:) ] -> [ op(R)(i+di,:) ]
175  // [ -s c ] [ op(R)(i ,:) ] [ op(R)(i ,:) ]
176  DVectorSlice
177  opR_row_idi = row(R->gms(),R_trans,i+di),
178  opR_row_i = row(R->gms(),R_trans,i);
179  if( opR_uplo == lower ) {
180  // op(R)(i ,:) = [ x x x w(i+wo) ] i
181  // op(R)(i+di,:) = [ x x x x ]
182  // i
183  rot( c, s, &opR_row_idi(1 ,i ), &opR_row_i(1,i) ); // First nonzero columns
184  rot( c, s, &opR_row_idi(i+1,i+1), &(*w)(i+wo,i+wo) ); // Last nonzero column
185  }
186  else {
187  // op(R)(i+di,:) = [ x x x x ]
188  // op(R)(i ,:) = [ w(i+wo) x x x ] i
189  // i
190  rot( c, s, &opR_row_idi(i-1,i-1), &(*w)(i+wo,i+wo) ); // First nonzero column
191  rot( c, s, &opR_row_idi(i ,n ), &opR_row_i(i,n) ); // Last nonzero columns
192  }
193  }}
194  //
195  // op(R1) + a * ||u||2 * e(last_i) * v' -> op(R2)
196  //
197  Vp_StV( &row(R->gms(),R_trans,last_i), a * (*u)(last_i), v );
198  //
199  // Apply (k-1) more givens rotations to return op(R2) back to triangular form:
200  //
201  // Q2 * R2 -> R_new
202  //
203  {for( size_type i = last_i; i != k; i-=di ) {
204  DVectorSlice
205  opR_row_i = row(R->gms(),R_trans,i),
206  opR_row_idi = row(R->gms(),R_trans,i-di);
207  value_type
208  &w_idiwo = (*w)(i-di+wo);
209  // [ c s ] [ op(R)(i,i) ] [ op(R)(i,i) ]
210  // [ -s c ] * [ op(R)(i-id,i) ] -> [ 0 ] w(i-di+wo)
211  value_type &R_i_i = opR_row_i(i);
212  value_type c, s; // Here c and s are computed to zero out w(i-di+wo)
213  rotg( &R_i_i, &w_idiwo, &c, &s );
214  // Make sure the diagonal is positive
215  value_type scale = +1.0;
216  if( R_i_i < 0.0 ) {
217  R_i_i = -R_i_i;
218  scale = -1.0;
219  w_idiwo = -w_idiwo; // Must also change this stored value
220  }
221  // [ c s ] * [ op(R)(i,:) ] -> [ op(R)(i,:) ]
222  // [ -s c ] [ op(R)(i-id,:)] [ op(R)(i-id,:) ]
223  if( opR_uplo == lower ) {
224  // op(R)(i-di,:) = [ x x x 0 ]
225  // op(R)(i,:) = [ x x x x ] i
226  // i
227  rot( scale*c, scale*s, &opR_row_i(1,i-1), &opR_row_idi(1,i-1) );
228  }
229  else {
230  // op(R)(i,:) = [ x x x x ] i
231  // op(R)(i-di,:) = [ 0 x x x ]
232  // i
233  rot( scale*c, scale*s, &opR_row_i(i+1,n), &opR_row_idi(i+1,n) );
234  }
235  }}
236  // When we get here note that u contains the first set of rotations Q1 and
237  // w contains the second set of rotations Q2.
238 }
Uplo
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
size_t size_type
Transp