Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_ReducedQuadratureFactoryImp.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) 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 Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42#include "Teuchos_TestForException.hpp"
43#include "Teuchos_SerialDenseVector.hpp"
44#include "Stokhos_SDMUtils.hpp"
47
48#include "Stokhos_ConfigDefs.h"
49#ifdef HAVE_STOKHOS_GLPK
50extern "C" {
51#include "glpk.h"
52}
53#endif
54
55#ifdef HAVE_STOKHOS_CLP
56#include "coin/ClpSimplex.hpp"
57#include "coin/CoinBuild.hpp"
58#include "coin/ClpInterior.hpp"
59#endif
60
61#ifdef HAVE_STOKHOS_QPOASES
62#include "qpOASES.hpp"
63#endif
64
65#ifdef HAVE_STOKHOS_DAKOTA
66#include "LinearSolver.hpp"
67#endif
68
69template <typename ordinal_type, typename value_type>
72 const Teuchos::ParameterList& params_) :
73 params(params_),
74 reduction_method(params.get("Reduced Quadrature Method", "Q Squared")),
75 solver_method(params.get("Solver Method", "TRSM")),
76 eliminate_dependent_rows(params.get("Eliminate Dependent Rows", true)),
77 verbose(params.get("Verbose", false)),
78 reduction_tol(params.get("Reduction Tolerance", 1.0e-12)),
79 objective_value(params.get("Objective Value", 0.0))
80{
81}
82
83template <typename ordinal_type, typename value_type>
84Teuchos::RCP<const Stokhos::UserDefinedQuadrature<ordinal_type, value_type> >
87 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Q,
88 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Q2,
89 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& F,
90 const Teuchos::Array<value_type>& weights) const
91{
92 Teuchos::RCP< Teuchos::Array<value_type> > red_weights;
93 Teuchos::RCP< Teuchos::Array< Teuchos::Array<value_type> > > red_points;
94 Teuchos::RCP< Teuchos::Array< Teuchos::Array<value_type> > > red_values;
95
96 // Compute reduced quadrature rule
97 if (reduction_method == "Q Squared") {
98 if (eliminate_dependent_rows)
99 reducedQuadrature_Q_Squared_CPQR(Q, F, weights,
100 red_weights, red_points, red_values);
101 else
102 reducedQuadrature_Q_Squared(Q, F, weights,
103 red_weights, red_points, red_values);
104 }
105 else if (reduction_method == "Q Squared2") {
106 reducedQuadrature_Q_Squared_CPQR2(Q, F, weights,
107 red_weights, red_points, red_values);
108 }
109 else if (reduction_method == "Q2") {
110 if (eliminate_dependent_rows)
111 reducedQuadrature_Q2_CPQR(Q, Q2, F, weights,
112 red_weights, red_points, red_values);
113 else
114 reducedQuadrature_Q2(Q, Q2, F, weights,
115 red_weights, red_points, red_values);
116 }
117 else if (reduction_method == "None") {
118 ordinal_type sz = Q.numCols();
119 ordinal_type nqp = Q.numRows();
120 ordinal_type d = F.numCols();
121 red_weights =
122 Teuchos::rcp(new Teuchos::Array<value_type>(weights));
123 red_points =
124 Teuchos::rcp(new Teuchos::Array< Teuchos::Array<value_type> >(nqp));
125 red_values =
126 Teuchos::rcp(new Teuchos::Array< Teuchos::Array<value_type> >(nqp));
127 for (ordinal_type i=0; i<nqp; i++) {
128 (*red_points)[i].resize(d);
129 for (ordinal_type j=0; j<d; j++)
130 (*red_points)[i][j] = F(i,j);
131 (*red_values)[i].resize(sz);
132 for (ordinal_type j=0; j<sz; j++)
133 (*red_values)[i][j] = Q(i,j);
134 }
135 }
136 else
137 TEUCHOS_TEST_FOR_EXCEPTION(
138 true, std::logic_error,
139 "Invalid dimension reduction method " << reduction_method);
140
141 // Build reduced quadrature object
142 Teuchos::RCP< const Teuchos::Array<value_type> > cred_weights =
143 red_weights;
144 Teuchos::RCP< const Teuchos::Array< Teuchos::Array<value_type> > > cred_points = red_points;
145 Teuchos::RCP< const Teuchos::Array< Teuchos::Array<value_type> > > cred_values = red_values;
146
147 Teuchos::RCP<const Stokhos::UserDefinedQuadrature<ordinal_type, value_type> >
148 reduced_quad =
149 Teuchos::rcp(new UserDefinedQuadrature<ordinal_type,value_type>(
150 cred_points,
151 cred_weights,
152 cred_values));
153
154 return reduced_quad;
155}
156
157template <typename ordinal_type, typename value_type>
158void
161 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Q,
162 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& F,
163 const Teuchos::Array<value_type>& weights,
164 Teuchos::RCP< Teuchos::Array<value_type> >& reduced_weights,
165 Teuchos::RCP< Teuchos::Array< Teuchos::Array<value_type> > >& reduced_points,
166 Teuchos::RCP< Teuchos::Array< Teuchos::Array<value_type> > >& reduced_values
167 ) const
168{
169 ordinal_type sz = Q.numCols();
170 ordinal_type sz2 = sz*(sz+1)/2;
171 ordinal_type nqp = Q.numRows();
172 ordinal_type d = F.numCols();
173
174 // Compute Q-squared matrix with all possible products
175 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Q2(nqp, sz2);
176 ordinal_type jdx=0;
177 for (ordinal_type j=0; j<sz; j++) {
178 for (ordinal_type k=j; k<sz; k++) {
179 for (ordinal_type i=0; i<nqp; i++)
180 Q2(i,jdx) = Q(i,j)*Q(i,k);
181 jdx++;
182 }
183 }
184 TEUCHOS_ASSERT(jdx == sz2);
185
186 // Compute e_1 = Q2^T*w which is our constraint
187 Teuchos::SerialDenseVector<ordinal_type,value_type> e1(sz2);
188 Teuchos::SerialDenseVector<ordinal_type,value_type> w(
189 Teuchos::View, const_cast<value_type*>(&weights[0]), nqp);
190 ordinal_type ret =
191 e1.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q2, w, 0.0);
192 TEUCHOS_ASSERT(ret == 0);
193
194 // Solve problem
195 Teuchos::SerialDenseVector<ordinal_type,value_type> u(nqp);
196 underdetermined_solver(Q2, e1, u, Teuchos::TRANS, Teuchos::UNDEF_TRI);
197
198 if (verbose) {
199 std::cout << "sz = " << sz << std::endl;
200 std::cout << "nqp = " << nqp << std::endl;
201 std::cout << "sz2 = " << sz2 << std::endl;
202
203 //std::cout << "reduced weights = " << u << std::endl;
204
205 // Check residual error ||e1-Q2^T*u||
206 Teuchos::SerialDenseVector<ordinal_type,value_type> err1(e1);
207 ret = err1.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, -1.0, Q2, u, 1.0);
208 TEUCHOS_ASSERT(ret == 0);
209 std::cout << "||e1-Q2^T*u||_infty = " << err1.normInf() << std::endl;
210
211 // Check discrete orthogonality error ||I - Q^T*diag(u)*Q||
212 Teuchos::SerialDenseMatrix<ordinal_type,value_type> err2(sz, sz);
213 err2.putScalar(0.0);
214 for (ordinal_type i=0; i<sz; i++)
215 err2(i,i) = 1.0;
216 Teuchos::SerialDenseMatrix<ordinal_type,value_type> WQ(nqp, sz);
217 for (ordinal_type i=0; i<nqp; i++)
218 for (ordinal_type j=0; j<sz; j++)
219 WQ(i,j) = u[i]*Q(i,j);
220 ret = err2.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, -1.0, Q, WQ, 1.0);
221 TEUCHOS_ASSERT(ret == 0);
222 std::cout << "||vec(I-Q^T*diag(u)*Q)||_infty = " << vec_norm_inf(err2)
223 << std::endl;
224 //print_matlab(std::cout, err2);
225 }
226
227 ordinal_type rank = 0;
228 for (ordinal_type i=0; i<nqp; i++)
229 if (std::abs(u[i]) > reduction_tol) ++rank;
230
231 // Get reduced weights, points and values
232 reduced_weights =
233 Teuchos::rcp(new Teuchos::Array<value_type>(rank));
234 reduced_points =
235 Teuchos::rcp(new Teuchos::Array< Teuchos::Array<value_type> >(rank));
236 reduced_values =
237 Teuchos::rcp(new Teuchos::Array< Teuchos::Array<value_type> >(rank));
238 ordinal_type idx = 0;
239 for (ordinal_type i=0; i<nqp; i++) {
240 if (std::abs(u[i]) > reduction_tol) {
241 (*reduced_weights)[idx] = u[i];
242 (*reduced_points)[idx].resize(d);
243 for (ordinal_type j=0; j<d; j++)
244 (*reduced_points)[idx][j] = F(i,j);
245 (*reduced_values)[idx].resize(sz);
246 for (ordinal_type j=0; j<sz; j++)
247 (*reduced_values)[idx][j] = Q(i,j);
248 idx++;
249 }
250 }
251
252 // idx may be less than rank if we obtained zero weights in solving
253 // the least-squares problem
254 TEUCHOS_ASSERT(idx <= rank);
255 if (idx < rank) {
256 rank = idx;
257 reduced_weights->resize(rank);
258 reduced_points->resize(rank);
259 reduced_values->resize(rank);
260 }
261}
262
263template <typename ordinal_type, typename value_type>
264void
267 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Q,
268 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& F,
269 const Teuchos::Array<value_type>& weights,
270 Teuchos::RCP< Teuchos::Array<value_type> >& reduced_weights,
271 Teuchos::RCP< Teuchos::Array< Teuchos::Array<value_type> > >& reduced_points,
272 Teuchos::RCP< Teuchos::Array< Teuchos::Array<value_type> > >& reduced_values
273 ) const
274{
275 ordinal_type sz = Q.numCols();
276 ordinal_type sz2 = sz*(sz+1)/2;
277 ordinal_type nqp = Q.numRows();
278 ordinal_type d = F.numCols();
279
280 // Compute Q-squared matrix with all possible products
281 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Q2(nqp, sz2);
282 ordinal_type jdx=0;
283 for (ordinal_type j=0; j<sz; j++) {
284 for (ordinal_type k=j; k<sz; k++) {
285 for (ordinal_type i=0; i<nqp; i++)
286 Q2(i,jdx) = Q(i,j)*Q(i,k);
287 jdx++;
288 }
289 }
290 TEUCHOS_ASSERT(jdx == sz2);
291
292 // Compute e_1 = Q2^T*w which is our constraint
293 Teuchos::SerialDenseVector<ordinal_type,value_type> e1(sz2);
294 Teuchos::SerialDenseVector<ordinal_type,value_type> w(
295 Teuchos::View, const_cast<value_type*>(&weights[0]), nqp);
296 ordinal_type ret =
297 e1.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q2, w, 0.0);
298 TEUCHOS_ASSERT(ret == 0);
299
300 // Compute QR decomposition of Q2^T: Q2^T*P = Z*R
301 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Q2t(Q2, Teuchos::TRANS);
302 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Z, R;
303 Teuchos::Array<ordinal_type> piv;
304 CPQR_Householder3(Q2t, Z, R, piv);
305 ordinal_type r = computeRank(R, reduction_tol);
306 R.reshape(r, R.numCols());
307 Z.reshape(Z.numRows(), r);
308
309 if (verbose) {
310 std::cout << "Q2 rank = " << r << std::endl;
311 }
312
313 // Compute new constraint vector
314 Teuchos::SerialDenseVector<ordinal_type,value_type> e1t(r);
315 ret = e1t.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Z, e1, 0.0);
316 TEUCHOS_ASSERT(ret == 0);
317
318 // Solve problem
319 Teuchos::SerialDenseVector<ordinal_type,value_type> ut(nqp);
320 underdetermined_solver(R, e1t, ut, Teuchos::NO_TRANS, Teuchos::UPPER_TRI);
321
322 // Transform solution back to original
323 Teuchos::SerialDenseVector<ordinal_type,value_type> u(nqp);
324 for (ordinal_type i=0; i<nqp; i++)
325 u[piv[i]] = ut[i];
326
327 if (verbose) {
328 std::cout << "sz = " << sz << std::endl;
329 std::cout << "nqp = " << nqp << std::endl;
330 std::cout << "sz2 = " << sz2 << std::endl;
331
332 //std::cout << "reduced weights = " << u << std::endl;
333
334 // Check residual error ||e1-Q2^T*u||
335 Teuchos::SerialDenseVector<ordinal_type,value_type> err1(e1);
336 ret = err1.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, -1.0, Q2, u, 1.0);
337 TEUCHOS_ASSERT(ret == 0);
338 std::cout << "||e1-Q2^T*u||_infty = " << err1.normInf() << std::endl;
339
340 // Check discrete orthogonality error ||I - Q^T*diag(u)*Q||
341 Teuchos::SerialDenseMatrix<ordinal_type,value_type> err2(sz, sz);
342 err2.putScalar(0.0);
343 for (ordinal_type i=0; i<sz; i++)
344 err2(i,i) = 1.0;
345 Teuchos::SerialDenseMatrix<ordinal_type,value_type> WQ(nqp, sz);
346 for (ordinal_type i=0; i<nqp; i++)
347 for (ordinal_type j=0; j<sz; j++)
348 WQ(i,j) = u[i]*Q(i,j);
349 ret = err2.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, -1.0, Q, WQ, 1.0);
350 TEUCHOS_ASSERT(ret == 0);
351 std::cout << "||vec(I-Q^T*diag(u)*Q)||_infty = " << vec_norm_inf(err2)
352 << std::endl;
353 //print_matlab(std::cout, err2);
354 }
355
356 ordinal_type rank = 0;
357 for (ordinal_type i=0; i<nqp; i++)
358 if (std::abs(u[i]) > reduction_tol) ++rank;
359
360 // Get reduced weights, points and values
361 reduced_weights =
362 Teuchos::rcp(new Teuchos::Array<value_type>(rank));
363 reduced_points =
364 Teuchos::rcp(new Teuchos::Array< Teuchos::Array<value_type> >(rank));
365 reduced_values =
366 Teuchos::rcp(new Teuchos::Array< Teuchos::Array<value_type> >(rank));
367 ordinal_type idx = 0;
368 for (ordinal_type i=0; i<nqp; i++) {
369 if (std::abs(u[i]) > reduction_tol) {
370 (*reduced_weights)[idx] = u[i];
371 (*reduced_points)[idx].resize(d);
372 for (ordinal_type j=0; j<d; j++)
373 (*reduced_points)[idx][j] = F(i,j);
374 (*reduced_values)[idx].resize(sz);
375 for (ordinal_type j=0; j<sz; j++)
376 (*reduced_values)[idx][j] = Q(i,j);
377 idx++;
378 }
379 }
380
381 // idx may be less than rank if we obtained zero weights in solving
382 // the least-squares problem
383 TEUCHOS_ASSERT(idx <= rank);
384 if (idx < rank) {
385 rank = idx;
386 reduced_weights->resize(rank);
387 reduced_points->resize(rank);
388 reduced_values->resize(rank);
389 }
390}
391
392template <typename ordinal_type, typename value_type>
393void
396 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Q,
397 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& F,
398 const Teuchos::Array<value_type>& weights,
399 Teuchos::RCP< Teuchos::Array<value_type> >& reduced_weights,
400 Teuchos::RCP< Teuchos::Array< Teuchos::Array<value_type> > >& reduced_points,
401 Teuchos::RCP< Teuchos::Array< Teuchos::Array<value_type> > >& reduced_values
402 ) const
403{
404 ordinal_type sz = Q.numCols();
405 ordinal_type sz2 = sz*(sz+1)/2;
406 ordinal_type nqp = Q.numRows();
407 ordinal_type d = F.numCols();
408
409 // Compute Q-squared matrix with all possible products
410 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Q2(nqp, sz2);
411 ordinal_type jdx=0;
412 for (ordinal_type j=0; j<sz; j++) {
413 for (ordinal_type k=j; k<sz; k++) {
414 for (ordinal_type i=0; i<nqp; i++)
415 Q2(i,jdx) = Q(i,j)*Q(i,k);
416 jdx++;
417 }
418 }
419 TEUCHOS_ASSERT(jdx == sz2);
420
421 // Compute QR decomposition of Q2: Q2*P = Z*R
422 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Z, R;
423 Teuchos::Array<ordinal_type> piv;
424 std::string orthogonalization_method =
425 params.get("Orthogonalization Method", "Householder");
426 Teuchos::Array<value_type> ww(weights);
427 if (orthogonalization_method == "Householder")
428 for (ordinal_type i=0; i<nqp; ++i) ww[i] = 1.0;
430 ordinal_type r = SOF::createOrthogonalBasis(
431 orthogonalization_method, reduction_tol, verbose, Q2, ww,
432 Z, R, piv);
433 bool restrict_r = params.get("Restrict Rank", false);
434 if (restrict_r) {
435 ordinal_type d = F.numCols();
436 ordinal_type p = params.get<ordinal_type>("Order Restriction");
437 ordinal_type n = n_choose_k(p+d,d);
438 if (r > n) {
439 if (verbose)
440 std::cout << "Restricting rank from " << r << " to " << n << std::endl;
441 r = n;
442 }
443 }
444
445 // Get the first r columns of Q2*P or Z
446 bool use_Z = params.get("Use Q in LP", true);
447 Teuchos::SerialDenseMatrix<ordinal_type, value_type> QQ2(nqp, r);
448 if (use_Z) {
449 for (ordinal_type j=0; j<r; j++)
450 for (ordinal_type i=0; i<nqp; i++)
451 QQ2(i,j) = Z(i,j);
452 }
453 else {
454 for (ordinal_type j=0; j<r; j++)
455 for (ordinal_type i=0; i<nqp; i++)
456 QQ2(i,j) = Q2(i,piv[j]);
457 }
458
459 if (verbose) {
460 std::cout << "Q2 rank = " << r << std::endl;
461 }
462
463 // Compute e_1 = QQ2^T*w which is our constraint
464 Teuchos::SerialDenseVector<ordinal_type,value_type> ee1(r);
465 Teuchos::SerialDenseVector<ordinal_type,value_type> w(
466 Teuchos::View, const_cast<value_type*>(&weights[0]), nqp);
467 ordinal_type ret =
468 ee1.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, QQ2, w, 0.0);
469 TEUCHOS_ASSERT(ret == 0);
470
471 // save_matlab("lp.mat", "A", QQ2, true);
472 // save_matlab("lp.mat", "b", ee1, false);
473
474 // Solve problem
475 Teuchos::SerialDenseVector<ordinal_type,value_type> u(nqp);
476 underdetermined_solver(QQ2, ee1, u, Teuchos::TRANS, Teuchos::UNDEF_TRI);
477
478 if (verbose) {
479 std::cout << "sz = " << sz << std::endl;
480 std::cout << "nqp = " << nqp << std::endl;
481 std::cout << "sz2 = " << sz2 << std::endl;
482
483 //std::cout << "reduced weights = " << u << std::endl;
484
485 // Check residual error ||e1-Q2^T*u||
486 Teuchos::SerialDenseVector<ordinal_type,value_type> err1(sz2);
487 ret = err1.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q2, w, 0.0);
488 TEUCHOS_ASSERT(ret == 0);
489 ret = err1.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, -1.0, Q2, u, 1.0);
490 TEUCHOS_ASSERT(ret == 0);
491 std::cout << "||e1-Q2^T*u||_infty = " << err1.normInf() << std::endl;
492
493 // Check discrete orthogonality error ||I - Q^T*diag(u)*Q||
494 Teuchos::SerialDenseMatrix<ordinal_type,value_type> err2(sz, sz);
495 err2.putScalar(0.0);
496 for (ordinal_type i=0; i<sz; i++)
497 err2(i,i) = 1.0;
498 Teuchos::SerialDenseMatrix<ordinal_type,value_type> WQ(nqp, sz);
499 for (ordinal_type i=0; i<nqp; i++)
500 for (ordinal_type j=0; j<sz; j++)
501 WQ(i,j) = u[i]*Q(i,j);
502 ret = err2.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, -1.0, Q, WQ, 1.0);
503 TEUCHOS_ASSERT(ret == 0);
504 std::cout << "||vec(I-Q^T*diag(u)*Q)||_infty = " << vec_norm_inf(err2)
505 << std::endl;
506 //print_matlab(std::cout, err2);
507 }
508
509 ordinal_type rank = 0;
510 for (ordinal_type i=0; i<nqp; i++)
511 if (std::abs(u[i]) > reduction_tol) ++rank;
512
513 // Get reduced weights, points and values
514 reduced_weights =
515 Teuchos::rcp(new Teuchos::Array<value_type>(rank));
516 reduced_points =
517 Teuchos::rcp(new Teuchos::Array< Teuchos::Array<value_type> >(rank));
518 reduced_values =
519 Teuchos::rcp(new Teuchos::Array< Teuchos::Array<value_type> >(rank));
520 ordinal_type idx = 0;
521 for (ordinal_type i=0; i<nqp; i++) {
522 if (std::abs(u[i]) > reduction_tol) {
523 (*reduced_weights)[idx] = u[i];
524 (*reduced_points)[idx].resize(d);
525 for (ordinal_type j=0; j<d; j++)
526 (*reduced_points)[idx][j] = F(i,j);
527 (*reduced_values)[idx].resize(sz);
528 for (ordinal_type j=0; j<sz; j++)
529 (*reduced_values)[idx][j] = Q(i,j);
530 idx++;
531 }
532 }
533
534 // idx may be less than rank if we obtained zero weights in solving
535 // the least-squares problem
536 TEUCHOS_ASSERT(idx <= rank);
537 if (idx < rank) {
538 rank = idx;
539 reduced_weights->resize(rank);
540 reduced_points->resize(rank);
541 reduced_values->resize(rank);
542 }
543}
544
545template <typename ordinal_type, typename value_type>
546void
549 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Q,
550 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Q2,
551 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& F,
552 const Teuchos::Array<value_type>& weights,
553 Teuchos::RCP< Teuchos::Array<value_type> >& reduced_weights,
554 Teuchos::RCP< Teuchos::Array< Teuchos::Array<value_type> > >& reduced_points,
555 Teuchos::RCP< Teuchos::Array< Teuchos::Array<value_type> > >& reduced_values
556 ) const
557{
558
559 ordinal_type sz = Q.numCols();
560 ordinal_type nqp = Q.numRows();
561 ordinal_type sz2 = Q2.numCols();
562 ordinal_type d = F.numCols();
563
564 // Compute e_1 = Q2^T*w which is our constraint
565 Teuchos::SerialDenseVector<ordinal_type,value_type> w(
566 Teuchos::View, const_cast<value_type*>(&weights[0]), nqp);
567 Teuchos::SerialDenseVector<ordinal_type,value_type> e1(sz2);
568 ordinal_type ret =
569 e1.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q2, w, 0.0);
570 TEUCHOS_ASSERT(ret == 0);
571
572 // Solve problem
573 Teuchos::SerialDenseVector<ordinal_type,value_type> u(nqp);
574 underdetermined_solver(Q2, e1, u, Teuchos::TRANS, Teuchos::UNDEF_TRI);
575
576 if (verbose) {
577 std::cout << "sz = " << sz << std::endl;
578 std::cout << "nqp = " << nqp << std::endl;
579 std::cout << "sz2 = " << sz2 << std::endl;
580
581 //std::cout << "reduced weights = " << u << std::endl;
582
583 // Check residual error ||e1-B2^T*u||
584 Teuchos::SerialDenseVector<ordinal_type,value_type> err1(e1);
585 ret = err1.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, -1.0, Q2, u, 1.0);
586 TEUCHOS_ASSERT(ret == 0);
587 std::cout << "||e1-Q2^T*u||_infty = " << err1.normInf() << std::endl;
588
589 // Check discrete orthogonality error ||I - Q^T*diag(u)*Q||
590 Teuchos::SerialDenseMatrix<ordinal_type,value_type> err2(sz, sz);
591 err2.putScalar(0.0);
592 for (ordinal_type i=0; i<sz; i++)
593 err2(i,i) = 1.0;
594 Teuchos::SerialDenseMatrix<ordinal_type,value_type> WQ(nqp, sz);
595 for (ordinal_type i=0; i<nqp; i++)
596 for (ordinal_type j=0; j<sz; j++)
597 WQ(i,j) = u[i]*Q(i,j);
598 ret = err2.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, -1.0, Q, WQ, 1.0);
599 TEUCHOS_ASSERT(ret == 0);
600 std::cout << "||vec(I-Q^T*diag(u)*Q)||_infty = " << vec_norm_inf(err2)
601 << std::endl;
602 }
603
604 ordinal_type rank = 0;
605 for (ordinal_type i=0; i<nqp; i++)
606 if (std::abs(u[i]) > reduction_tol) ++rank;
607
608 // Get reduced weights, points and values
609 reduced_weights =
610 Teuchos::rcp(new Teuchos::Array<value_type>(rank));
611 reduced_points =
612 Teuchos::rcp(new Teuchos::Array< Teuchos::Array<value_type> >(rank));
613 reduced_values =
614 Teuchos::rcp(new Teuchos::Array< Teuchos::Array<value_type> >(rank));
615 ordinal_type idx = 0;
616 for (ordinal_type i=0; i<nqp; i++) {
617 if (std::abs(u[i]) > reduction_tol) {
618 (*reduced_weights)[idx] = u[i];
619 (*reduced_points)[idx].resize(d);
620 for (ordinal_type j=0; j<d; j++)
621 (*reduced_points)[idx][j] = F(i,j);
622 (*reduced_values)[idx].resize(sz);
623 for (ordinal_type j=0; j<sz; j++)
624 (*reduced_values)[idx][j] = Q(i,j);
625 idx++;
626 }
627 }
628
629 // idx may be less than rank if we obtained zero weights in solving
630 // the least-squares problem
631 TEUCHOS_ASSERT(idx <= rank);
632 if (idx < rank) {
633 rank = idx;
634 reduced_weights->resize(rank);
635 reduced_points->resize(rank);
636 reduced_values->resize(rank);
637 }
638}
639
640template <typename ordinal_type, typename value_type>
641void
644 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Q,
645 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Q2,
646 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& F,
647 const Teuchos::Array<value_type>& weights,
648 Teuchos::RCP< Teuchos::Array<value_type> >& reduced_weights,
649 Teuchos::RCP< Teuchos::Array< Teuchos::Array<value_type> > >& reduced_points,
650 Teuchos::RCP< Teuchos::Array< Teuchos::Array<value_type> > >& reduced_values
651 ) const
652{
653
654 ordinal_type sz = Q.numCols();
655 ordinal_type nqp = Q.numRows();
656 ordinal_type sz2 = Q2.numCols();
657 ordinal_type d = F.numCols();
658
659 // Compute e_1 = Q2^T*w which is our constraint
660 Teuchos::SerialDenseVector<ordinal_type,value_type> w(
661 Teuchos::View, const_cast<value_type*>(&weights[0]), nqp);
662 Teuchos::SerialDenseVector<ordinal_type,value_type> e1(sz2);
663 ordinal_type ret =
664 e1.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q2, w, 0.0);
665 TEUCHOS_ASSERT(ret == 0);
666
667 // Compute QR decomposition of Q2^T: Q2^T*P = Z*R
668 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Q2t(Q2, Teuchos::TRANS);
669 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Z, R;
670 Teuchos::Array<ordinal_type> piv;
671 CPQR_Householder3(Q2t, Z, R, piv);
672 ordinal_type r = computeRank(R, reduction_tol);
673 R.reshape(r, R.numCols());
674 Z.reshape(Z.numRows(), r);
675
676 if (verbose) {
677 std::cout << "Q2 rank = " << r << std::endl;
678 }
679
680 // Compute new constraint vector
681 Teuchos::SerialDenseVector<ordinal_type,value_type> e1t(r);
682 ret = e1t.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Z, e1, 0.0);
683 TEUCHOS_ASSERT(ret == 0);
684
685 // Solve problem
686 Teuchos::SerialDenseVector<ordinal_type,value_type> ut(nqp);
687 underdetermined_solver(R, e1t, ut, Teuchos::NO_TRANS, Teuchos::UPPER_TRI);
688
689 // Transform solution back to original
690 Teuchos::SerialDenseVector<ordinal_type,value_type> u(nqp);
691 for (ordinal_type i=0; i<nqp; i++)
692 u[piv[i]] = ut[i];
693
694 if (verbose) {
695 std::cout << "sz = " << sz << std::endl;
696 std::cout << "nqp = " << nqp << std::endl;
697 std::cout << "sz2 = " << sz2 << std::endl;
698
699 //std::cout << "reduced weights = " << u << std::endl;
700
701 // Check residual error ||e1-B2^T*u||
702 Teuchos::SerialDenseVector<ordinal_type,value_type> err1(e1);
703 ret = err1.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, -1.0, Q2, u, 1.0);
704 TEUCHOS_ASSERT(ret == 0);
705 std::cout << "||e1-Q2^T*u||_infty = " << err1.normInf() << std::endl;
706
707 // Check discrete orthogonality error ||I - Q^T*diag(u)*Q||
708 Teuchos::SerialDenseMatrix<ordinal_type,value_type> err2(sz, sz);
709 err2.putScalar(0.0);
710 for (ordinal_type i=0; i<sz; i++)
711 err2(i,i) = 1.0;
712 Teuchos::SerialDenseMatrix<ordinal_type,value_type> WQ(nqp, sz);
713 for (ordinal_type i=0; i<nqp; i++)
714 for (ordinal_type j=0; j<sz; j++)
715 WQ(i,j) = u[i]*Q(i,j);
716 ret = err2.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, -1.0, Q, WQ, 1.0);
717 TEUCHOS_ASSERT(ret == 0);
718 std::cout << "||vec(I-Q^T*diag(u)*Q)||_infty = " << vec_norm_inf(err2)
719 << std::endl;
720 //print_matlab(std::cout, err2);
721 }
722
723 ordinal_type rank = 0;
724 for (ordinal_type i=0; i<nqp; i++)
725 if (std::abs(u[i]) > reduction_tol) ++rank;
726
727 // Get reduced weights, points and values
728 reduced_weights =
729 Teuchos::rcp(new Teuchos::Array<value_type>(rank));
730 reduced_points =
731 Teuchos::rcp(new Teuchos::Array< Teuchos::Array<value_type> >(rank));
732 reduced_values =
733 Teuchos::rcp(new Teuchos::Array< Teuchos::Array<value_type> >(rank));
734 ordinal_type idx = 0;
735 for (ordinal_type i=0; i<nqp; i++) {
736 if (std::abs(u[i]) > reduction_tol) {
737 (*reduced_weights)[idx] = u[i];
738 (*reduced_points)[idx].resize(d);
739 for (ordinal_type j=0; j<d; j++)
740 (*reduced_points)[idx][j] = F(i,j);
741 (*reduced_values)[idx].resize(sz);
742 for (ordinal_type j=0; j<sz; j++)
743 (*reduced_values)[idx][j] = Q(i,j);
744 idx++;
745 }
746 }
747
748 // idx may be less than rank if we obtained zero weights in solving
749 // the least-squares problem
750 TEUCHOS_ASSERT(idx <= rank);
751 if (idx < rank) {
752 rank = idx;
753 reduced_weights->resize(rank);
754 reduced_points->resize(rank);
755 reduced_values->resize(rank);
756 }
757}
758
759template <typename ordinal_type, typename value_type>
760void
763 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& A,
764 const Teuchos::SerialDenseVector<ordinal_type, value_type>& b,
765 Teuchos::SerialDenseVector<ordinal_type, value_type>& x,
766 Teuchos::ETransp transa, Teuchos::EUplo uplo) const
767{
768 if (solver_method == "TRSM")
769 solver_TRSM(A, b, x, transa, uplo);
770 else if (solver_method == "Clp")
771 solver_CLP(A, b, x, transa, uplo);
772 else if (solver_method == "Clp-IP")
773 solver_CLP_IP(A, b, x, transa, uplo);
774 else if (solver_method == "GLPK")
775 solver_GLPK(A, b, x, transa, uplo);
776 else if (solver_method == "qpOASES")
777 solver_qpOASES(A, b, x, transa, uplo);
778 else if (solver_method == "Compressed Sensing")
779 solver_CompressedSensing(A, b, x, transa, uplo);
780 else
781 TEUCHOS_TEST_FOR_EXCEPTION(
782 true, std::logic_error,
783 "Invalid solver method " << solver_method);
784}
785
786template <typename ordinal_type, typename value_type>
787void
790 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& A,
791 const Teuchos::SerialDenseVector<ordinal_type, value_type>& b,
792 Teuchos::SerialDenseVector<ordinal_type, value_type>& x,
793 Teuchos::ETransp transa, Teuchos::EUplo uplo) const
794{
795 TEUCHOS_TEST_FOR_EXCEPTION(uplo == Teuchos::UNDEF_TRI, std::logic_error,
796 "TRSM solver requires triangular matrix!");
797 ordinal_type m = A.numRows();
798 ordinal_type n = A.numCols();
799 ordinal_type n_rows;
800 ordinal_type n_cols;
801 if (transa == Teuchos::NO_TRANS) {
802 n_rows = m;
803 n_cols = n;
804 }
805 else {
806 n_rows = n;
807 n_cols = m;
808 }
809 if (x.length() < n_cols)
810 x.shape(n_cols, 1);
811 x.putScalar(0.0);
812 blas.COPY(n_rows, b.values(), 1, x.values(), 1);
813
814 // Here we assume A is upper triangular
815 blas.TRSM(Teuchos::LEFT_SIDE, uplo, transa,
816 Teuchos::NON_UNIT_DIAG, n_rows, 1, 1.0, A.values(), A.stride(),
817 x.values(), x.stride());
818}
819
820template <typename ordinal_type, typename value_type>
821void
824 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& A,
825 const Teuchos::SerialDenseVector<ordinal_type, value_type>& b,
826 Teuchos::SerialDenseVector<ordinal_type, value_type>& x,
827 Teuchos::ETransp transa, Teuchos::EUplo uplo) const
828{
829#ifdef HAVE_STOKHOS_GLPK
830 ordinal_type m = A.numRows();
831 ordinal_type n = A.numCols();
832 ordinal_type n_rows;
833 ordinal_type n_cols;
834 if (transa == Teuchos::NO_TRANS) {
835 n_rows = m;
836 n_cols = n;
837 }
838 else {
839 n_rows = n;
840 n_cols = m;
841 }
842 if (x.length() < n_cols)
843 x.shape(n_cols, 1);
844
845 LPX *lp = lpx_create_prob();
846 lpx_set_prob_name(lp, "Reduced Quadrature");
847 lpx_set_obj_dir(lp, LPX_MIN);
848 lpx_add_rows(lp, n_rows);
849 lpx_add_cols(lp, n_cols);
850
851 // Set row bounds
852 for (ordinal_type i=0; i<n_rows; i++)
853 lpx_set_row_bnds(lp, i+1, LPX_FX, b[i], b[i]);
854
855 // Set columns bounds and object coefficients
856 for (ordinal_type j=0; j<n_cols; j++) {
857 lpx_set_col_bnds(lp, j+1, LPX_LO, 0.0, 0.0);
858 lpx_set_obj_coef(lp, j+1, objective_value);
859 }
860
861 // Set constraint matrix
862 // GLPK uses this silly 1-based indexing scheme which essentially
863 // requires us to copy the matrix. While copying we will also transpose
864 // to make loading in GLPK easier
865 Teuchos::SerialDenseMatrix<ordinal_type, value_type> AA(n_rows+1,n_cols+1);
866 AA.putScalar(0.0);
867 Teuchos::EUplo AA_uplo = uplo;
868 if (transa == Teuchos::NO_TRANS) {
869 Teuchos::SerialDenseMatrix<ordinal_type, value_type> AAA(
870 Teuchos::View, AA, n_rows, n_cols, 1, 1);
871 AAA.assign(A);
872 }
873 else {
874 for (ordinal_type i=0; i<n_rows; i++)
875 for (ordinal_type j=0; j<n_cols; j++)
876 AA(i+1,j+1) = A(j,i);
877 if (uplo == Teuchos::UPPER_TRI)
878 AA_uplo = Teuchos::LOWER_TRI;
879 else if (uplo == Teuchos::LOWER_TRI)
880 AA_uplo = Teuchos::UPPER_TRI;
881 }
882 Teuchos::Array< Teuchos::Array<int> > row_indices(n_cols);
883 if (AA_uplo == Teuchos::UNDEF_TRI) {
884 for (ordinal_type j=0; j<n_cols; j++) {
885 row_indices[j].resize(n_rows+1);
886 for (ordinal_type i=0; i<n_rows; i++)
887 row_indices[j][i+1] = i+1;
888 lpx_set_mat_col(lp, j+1, n_rows, &row_indices[j][0], AA[j+1]);
889 }
890 }
891 else if (AA_uplo == Teuchos::UPPER_TRI) {
892 // For AA upper-triangular, only need to include leading
893 // min(n_rows,n_cols) x n_cols block (remaining rows must be zero)
894 for (ordinal_type j=0; j<n_cols; j++) {
895 ordinal_type nr = n_rows;
896 if (j < n_rows)
897 nr = j+1;
898 row_indices[j].resize(nr+1);
899 for (ordinal_type i=0; i<nr; i++)
900 row_indices[j][i+1] = i+1;
901 lpx_set_mat_col(lp, j+1, nr, &row_indices[j][0], AA[j+1]);
902 }
903 }
904 else if (AA_uplo == Teuchos::LOWER_TRI) {
905 // For AA lower-triangular, only need to include leading
906 // n_rows x min(n_rows,n_cols) block (remaining columns must be zero)
907 for (ordinal_type j=0; j<n_cols; j++) {
908 if (j < n_rows) {
909 ordinal_type nr = n_rows-j;
910 row_indices[j].resize(nr+1);
911 for (ordinal_type i=0; i<nr; i++)
912 row_indices[j][i+1] = j+i+1;
913 lpx_set_mat_col(lp, j+1, nr, &row_indices[j][0], AA[j+1]+j);
914 }
915 }
916 }
917
918 // Write MPS-file if requested
919 if (params.get("Write MPS File", false) == true)
920 lpx_write_mps(lp, params.get("MPS Filename", "lp.mps").c_str());
921
922 // Solve linear program
923 lpx_simplex(lp);
924 int status = lpx_get_status(lp);
925 if (verbose) {
926 std::cout << "glpk status = " << status << std::endl;
927 double z = lpx_get_obj_val(lp);
928 std::cout << "glpk objective = " << z << std::endl;
929 }
930
931 // Get solution
932 for (ordinal_type i=0; i<n_cols; i++)
933 x[i] = lpx_get_col_prim(lp, i+1);
934#else
935 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
936 "GLPK solver called but not enabled!");
937#endif
938}
939
940template <typename ordinal_type, typename value_type>
941void
944 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& A,
945 const Teuchos::SerialDenseVector<ordinal_type, value_type>& b,
946 Teuchos::SerialDenseVector<ordinal_type, value_type>& x,
947 Teuchos::ETransp transa, Teuchos::EUplo uplo) const
948{
949#ifdef HAVE_STOKHOS_CLP
950 ordinal_type m = A.numRows();
951 ordinal_type n = A.numCols();
952 ordinal_type n_rows;
953 ordinal_type n_cols;
954 if (transa == Teuchos::NO_TRANS) {
955 n_rows = m;
956 n_cols = n;
957 }
958 else {
959 n_rows = n;
960 n_cols = m;
961 }
962 if (x.length() < n_cols)
963 x.shape(n_cols, 1);
964
965 // Setup linear program
966 ClpSimplex model;
967 model.resize(n_rows, 0);
968
969 // Set row bounds
970 for (ordinal_type i=0; i<n_rows; i++) {
971 model.setRowLower(i, b[i]);
972 model.setRowUpper(i, b[i]);
973 }
974
975 // Set constraint matrix, columns bounds and objective coefficients
976 Teuchos::SerialDenseMatrix<ordinal_type, value_type> AA;
977 Teuchos::EUplo AA_uplo = uplo;
978 if (transa == Teuchos::NO_TRANS) {
979 AA = Teuchos::SerialDenseMatrix<ordinal_type, value_type>(
980 Teuchos::View, A, n_rows, n_cols);
981 }
982 else {
983 AA.reshape(n_rows, n_cols);
984 for (ordinal_type i=0; i<n_rows; i++)
985 for (ordinal_type j=0; j<n_cols; j++)
986 AA(i,j) = A(j,i);
987 if (uplo == Teuchos::UPPER_TRI)
988 AA_uplo = Teuchos::LOWER_TRI;
989 else if (uplo == Teuchos::LOWER_TRI)
990 AA_uplo = Teuchos::UPPER_TRI;
991 }
992 Teuchos::Array< Teuchos::Array<int> > row_indices(n_cols);
993 CoinBuild buildObject;
994 if (AA_uplo == Teuchos::UNDEF_TRI) {
995 for (ordinal_type j=0; j<n_cols; j++) {
996 row_indices[j].resize(n_rows);
997 for (ordinal_type i=0; i<n_rows; i++)
998 row_indices[j][i] = i;
999 buildObject.addColumn(
1000 n_rows, &row_indices[j][0], AA[j], 0.0, DBL_MAX, objective_value);
1001 }
1002 }
1003 else if (AA_uplo == Teuchos::UPPER_TRI) {
1004 // For AA upper-triangular, only need to include leading
1005 // min(n_rows,n_cols) x n_cols block (remaining rows must be zero)
1006 for (ordinal_type j=0; j<n_cols; j++) {
1007 ordinal_type nr = n_rows;
1008 if (j < n_rows)
1009 nr = j+1;
1010 row_indices[j].resize(nr);
1011 for (ordinal_type i=0; i<nr; i++)
1012 row_indices[j][i] = i;
1013 buildObject.addColumn(nr, &row_indices[j][0], AA[j], 0.0, DBL_MAX, objective_value);
1014 }
1015 }
1016 else if (AA_uplo == Teuchos::LOWER_TRI) {
1017 // For AA lower-triangular, only need to include leading
1018 // n_rows x min(n_rows,n_cols) block (remaining columns must be zero)
1019 for (ordinal_type j=0; j<n_cols; j++) {
1020 if (j < n_rows) {
1021 ordinal_type nr = n_rows-j;
1022 row_indices[j].resize(nr);
1023 for (ordinal_type i=0; i<nr; i++)
1024 row_indices[j][i] = j+i;
1025 buildObject.addColumn(
1026 nr, &row_indices[j][0], AA[j]+j, 0.0, DBL_MAX, objective_value);
1027 }
1028 else
1029 buildObject.addColumn(0, NULL, NULL, 0.0, DBL_MAX, objective_value);
1030 }
1031 }
1032 model.addColumns(buildObject);
1033
1034 // Solve linear program
1035 model.primal();
1036
1037 // Get solution
1038 double *solution = model.primalColumnSolution();
1039 for (ordinal_type i=0; i<n_cols; i++)
1040 x[i] = solution[i];
1041#else
1042 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
1043 "CLP solver called but not enabled!");
1044#endif
1045}
1046
1047template <typename ordinal_type, typename value_type>
1048void
1051 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& A,
1052 const Teuchos::SerialDenseVector<ordinal_type, value_type>& b,
1053 Teuchos::SerialDenseVector<ordinal_type, value_type>& x,
1054 Teuchos::ETransp transa, Teuchos::EUplo uplo) const
1055{
1056#ifdef HAVE_STOKHOS_CLP
1057 ordinal_type m = A.numRows();
1058 ordinal_type n = A.numCols();
1059 ordinal_type n_rows;
1060 ordinal_type n_cols;
1061 if (transa == Teuchos::NO_TRANS) {
1062 n_rows = m;
1063 n_cols = n;
1064 }
1065 else {
1066 n_rows = n;
1067 n_cols = m;
1068 }
1069 if (x.length() < n_cols)
1070 x.shape(n_cols, 1);
1071
1072 // Setup linear program
1073 ClpInterior model;
1074 model.resize(n_rows, 0);
1075
1076 // Set row bounds
1077 for (ordinal_type i=0; i<n_rows; i++) {
1078 model.setRowLower(i, b[i]);
1079 model.setRowUpper(i, b[i]);
1080 }
1081
1082 // Set constraint matrix, columns bounds and objective coefficients
1083 Teuchos::SerialDenseMatrix<ordinal_type, value_type> AA;
1084 Teuchos::EUplo AA_uplo = uplo;
1085 if (transa == Teuchos::NO_TRANS) {
1086 AA = Teuchos::SerialDenseMatrix<ordinal_type, value_type>(
1087 Teuchos::View, A, n_rows, n_cols);
1088 }
1089 else {
1090 AA.reshape(n_rows, n_cols);
1091 for (ordinal_type i=0; i<n_rows; i++)
1092 for (ordinal_type j=0; j<n_cols; j++)
1093 AA(i,j) = A(j,i);
1094 if (uplo == Teuchos::UPPER_TRI)
1095 AA_uplo = Teuchos::LOWER_TRI;
1096 else if (uplo == Teuchos::LOWER_TRI)
1097 AA_uplo = Teuchos::UPPER_TRI;
1098 }
1099 Teuchos::Array< Teuchos::Array<int> > row_indices(n_cols);
1100 CoinBuild buildObject;
1101 if (AA_uplo == Teuchos::UNDEF_TRI) {
1102 for (ordinal_type j=0; j<n_cols; j++) {
1103 row_indices[j].resize(n_rows);
1104 for (ordinal_type i=0; i<n_rows; i++)
1105 row_indices[j][i] = i;
1106 buildObject.addColumn(
1107 n_rows, &row_indices[j][0], AA[j], 0.0, DBL_MAX, objective_value);
1108 }
1109 }
1110 else if (AA_uplo == Teuchos::UPPER_TRI) {
1111 // For AA upper-triangular, only need to include leading
1112 // min(n_rows,n_cols) x n_cols block (remaining rows must be zero)
1113 for (ordinal_type j=0; j<n_cols; j++) {
1114 ordinal_type nr = n_rows;
1115 if (j < n_rows)
1116 nr = j+1;
1117 row_indices[j].resize(nr);
1118 for (ordinal_type i=0; i<nr; i++)
1119 row_indices[j][i] = i;
1120 buildObject.addColumn(nr, &row_indices[j][0], AA[j], 0.0, DBL_MAX, objective_value);
1121 }
1122 }
1123 else if (AA_uplo == Teuchos::LOWER_TRI) {
1124 // For AA lower-triangular, only need to include leading
1125 // n_rows x min(n_rows,n_cols) block (remaining columns must be zero)
1126 for (ordinal_type j=0; j<n_cols; j++) {
1127 if (j < n_rows) {
1128 ordinal_type nr = n_rows-j;
1129 row_indices[j].resize(nr);
1130 for (ordinal_type i=0; i<nr; i++)
1131 row_indices[j][i] = j+i;
1132 buildObject.addColumn(
1133 nr, &row_indices[j][0], AA[j]+j, 0.0, DBL_MAX, objective_value);
1134 }
1135 else
1136 buildObject.addColumn(0, NULL, NULL, 0.0, DBL_MAX, objective_value);
1137 }
1138 }
1139 model.addColumns(buildObject);
1140
1141 // Solve linear program
1142 model.primalDual();
1143
1144 // Get solution
1145 double *solution = model.primalColumnSolution();
1146 for (ordinal_type i=0; i<n_cols; i++)
1147 x[i] = solution[i];
1148#else
1149 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
1150 "CLP solver called but not enabled!");
1151#endif
1152}
1153
1154template <typename ordinal_type, typename value_type>
1155void
1158 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& A,
1159 const Teuchos::SerialDenseVector<ordinal_type, value_type>& b,
1160 Teuchos::SerialDenseVector<ordinal_type, value_type>& x,
1161 Teuchos::ETransp transa, Teuchos::EUplo uplo) const
1162{
1163#ifdef HAVE_STOKHOS_QPOASES
1164 ordinal_type m = A.numRows();
1165 ordinal_type n = A.numCols();
1166 ordinal_type n_rows;
1167 ordinal_type n_cols;
1168 if (transa == Teuchos::NO_TRANS) {
1169 n_rows = m;
1170 n_cols = n;
1171 }
1172 else {
1173 n_rows = n;
1174 n_cols = m;
1175 }
1176 if (x.length() < n_cols)
1177 x.shape(n_cols, 1);
1178
1179 // Compute objective vector
1180 Teuchos::SerialDenseVector<ordinal_type,value_type> c(n_cols);
1181 c.putScalar(objective_value);
1182
1183 // Compute lower bounds
1184 Teuchos::SerialDenseVector<ordinal_type,value_type> lb(n_cols);
1185 lb.putScalar(0.0);
1186
1187 // Setup linear program
1188 qpOASES::QProblem lp(n_cols, n_rows, qpOASES::HST_ZERO);
1189
1190 // Solve linear program -- qpOASES expects row-wise storage
1191 // so transpose the matrix if necessary. Since qpOASES uses dense
1192 // linear algebra, can't take advantage of upper/lower triangular info
1193 int nWSR = params.get("Max Working Set Recalculations", 10000);
1194 if (transa == Teuchos::NO_TRANS) {
1195 Teuchos::SerialDenseMatrix<ordinal_type, value_type> AA(A, Teuchos::TRANS);
1196 lp.init(NULL, // zero Hessian
1197 c.values(), // objective
1198 AA.values(), // constrait matrix
1199 lb.values(), // variable lower bounds
1200 NULL, // no upper bounds
1201 b.values(), // constraint lower bounds
1202 b.values(), // constraint upper bounds
1203 nWSR, // maximum number of working set recalculations
1204 NULL // maximum CPU time
1205 );
1206 }
1207 else {
1208 lp.init(NULL, // zero Hessian
1209 c.values(), // objective
1210 A.values(), // constrait matrix
1211 lb.values(), // variable lower bounds
1212 NULL, // no upper bounds
1213 b.values(), // constraint lower bounds
1214 b.values(), // constraint upper bounds
1215 nWSR, // maximum number of working set recalculations
1216 NULL // maximum CPU time
1217 );
1218 }
1219
1220 // Get solution
1221 lp.getPrimalSolution(x.values());
1222#else
1223 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
1224 "qpOASES solver called but not enabled!");
1225#endif
1226}
1227
1228template <typename ordinal_type, typename value_type>
1229void
1232 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& A,
1233 const Teuchos::SerialDenseVector<ordinal_type, value_type>& b,
1234 Teuchos::SerialDenseVector<ordinal_type, value_type>& x,
1235 Teuchos::ETransp transa, Teuchos::EUplo uplo) const
1236{
1237#ifdef HAVE_STOKHOS_DAKOTA
1238 ordinal_type m = A.numRows();
1239 ordinal_type n = A.numCols();
1240 ordinal_type n_cols;
1241 if (transa == Teuchos::NO_TRANS) {
1242 n_cols = n;
1243 }
1244 else {
1245 n_cols = m;
1246 }
1247 if (x.length() < n_cols)
1248 x.shape(n_cols, 1);
1249
1250 // Setup L1 minimization problem
1251 // Todo: parse options from parameter list
1252 Pecos::CompressedSensingTool CS;
1253 Teuchos::SerialDenseMatrix<ordinal_type, value_type> AA(A, transa);
1254 Teuchos::SerialDenseVector<ordinal_type, value_type> bb(b);
1255 Pecos::RealMatrixArray xx;
1256 Pecos::CompressedSensingOptions opts;
1257 Pecos::CompressedSensingOptionsList opts_list;
1258 CS.solve(AA, bb, xx, opts, opts_list);
1259 x.assign(xx[0]);
1260
1261#else
1262 TEUCHOS_TEST_FOR_EXCEPTION(
1263 true, std::logic_error,
1264 "CompressedSensing solver called but not enabled!");
1265#endif
1266}
1267
1268template <typename ordinal_type, typename value_type>
1269ordinal_type
1272 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& R,
1273 const value_type tol) const
1274{
1275 // Compute rank -- since we constrain the first d+1 columns of Bp to lie
1276 // in the basis, the diagonal elements of Rp may not be monotonically
1277 // decreasing until we get past d+1
1278 ordinal_type rank = 0;
1279 ordinal_type m = R.numRows();
1280 value_type r_max = std::abs(R(rank,rank));
1281 value_type r_min = std::abs(R(rank,rank));
1282 for (rank=1; rank<m; rank++) {
1283 if (std::abs(R(rank,rank)) > r_max)
1284 r_max = std::abs(R(rank,rank));
1285 if (std::abs(R(rank,rank)) < r_min)
1286 r_min = std::abs(R(rank,rank));
1287 if (r_min / r_max < tol)
1288 break;
1289 }
1290
1291 // Print condition number of R
1292 if (verbose) {
1293 value_type cond_r = r_max / r_min;
1294 std::cout << "r_max = " << r_max << std::endl;
1295 std::cout << "r_min = " << r_min << std::endl;
1296 std::cout << "Condition number of R = " << cond_r << std::endl;
1297 }
1298
1299 return rank;
1300}
1301
1302template <typename ordinal_type, typename value_type>
1303ordinal_type
1305n_choose_k(const ordinal_type& n, const ordinal_type& k) const {
1306 // Use formula
1307 // n!/(k!(n-k)!) = n(n-1)...(k+1) / ( (n-k)(n-k-1)...1 ) ( n-k terms )
1308 // = n(n-1)...(n-k+1) / ( k(k-1)...1 ) ( k terms )
1309 // which ever has fewer terms
1310 if (k > n)
1311 return 0;
1312 ordinal_type num = 1;
1313 ordinal_type den = 1;
1314 ordinal_type l = std::min(n-k,k);
1315 for (ordinal_type i=0; i<l; i++) {
1316 num *= n-i;
1317 den *= i+1;
1318 }
1319 return num / den;
1320}
Encapsulate various orthogonalization (ie QR) methods.
void solver_CLP_IP(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &A, const Teuchos::SerialDenseVector< ordinal_type, value_type > &b, Teuchos::SerialDenseVector< ordinal_type, value_type > &x, Teuchos::ETransp transa, Teuchos::EUplo uplo) const
void reducedQuadrature_Q2_CPQR(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q2, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &F, const Teuchos::Array< value_type > &weights, Teuchos::RCP< Teuchos::Array< value_type > > &red_weights, Teuchos::RCP< Teuchos::Array< Teuchos::Array< value_type > > > &red_points, Teuchos::RCP< Teuchos::Array< Teuchos::Array< value_type > > > &red_values) const
void reducedQuadrature_Q_Squared_CPQR2(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &F, const Teuchos::Array< value_type > &weights, Teuchos::RCP< Teuchos::Array< value_type > > &red_weights, Teuchos::RCP< Teuchos::Array< Teuchos::Array< value_type > > > &red_points, Teuchos::RCP< Teuchos::Array< Teuchos::Array< value_type > > > &red_values) const
ordinal_type n_choose_k(const ordinal_type &n, const ordinal_type &k) const
Compute bionomial coefficient (n ; k) = n!/( k! (n-k)! )
void reducedQuadrature_Q_Squared(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &F, const Teuchos::Array< value_type > &weights, Teuchos::RCP< Teuchos::Array< value_type > > &red_weights, Teuchos::RCP< Teuchos::Array< Teuchos::Array< value_type > > > &red_points, Teuchos::RCP< Teuchos::Array< Teuchos::Array< value_type > > > &red_values) const
void reducedQuadrature_Q_Squared_CPQR(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &F, const Teuchos::Array< value_type > &weights, Teuchos::RCP< Teuchos::Array< value_type > > &red_weights, Teuchos::RCP< Teuchos::Array< Teuchos::Array< value_type > > > &red_points, Teuchos::RCP< Teuchos::Array< Teuchos::Array< value_type > > > &red_values) const
void underdetermined_solver(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &A, const Teuchos::SerialDenseVector< ordinal_type, value_type > &b, Teuchos::SerialDenseVector< ordinal_type, value_type > &x, Teuchos::ETransp transa, Teuchos::EUplo uplo) const
ordinal_type computeRank(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &R, const value_type tol) const
void solver_TRSM(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &A, const Teuchos::SerialDenseVector< ordinal_type, value_type > &b, Teuchos::SerialDenseVector< ordinal_type, value_type > &x, Teuchos::ETransp transa, Teuchos::EUplo uplo) const
virtual Teuchos::RCP< const Stokhos::UserDefinedQuadrature< ordinal_type, value_type > > createReducedQuadrature(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q2, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &F, const Teuchos::Array< value_type > &weights) const
Get reduced quadrature object.
void reducedQuadrature_Q2(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q2, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &F, const Teuchos::Array< value_type > &weights, Teuchos::RCP< Teuchos::Array< value_type > > &red_weights, Teuchos::RCP< Teuchos::Array< Teuchos::Array< value_type > > > &red_points, Teuchos::RCP< Teuchos::Array< Teuchos::Array< value_type > > > &red_values) const
void solver_CompressedSensing(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &A, const Teuchos::SerialDenseVector< ordinal_type, value_type > &b, Teuchos::SerialDenseVector< ordinal_type, value_type > &x, Teuchos::ETransp transa, Teuchos::EUplo uplo) const
ReducedQuadratureFactory(const Teuchos::ParameterList &params)
Constructor.
void solver_GLPK(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &A, const Teuchos::SerialDenseVector< ordinal_type, value_type > &b, Teuchos::SerialDenseVector< ordinal_type, value_type > &x, Teuchos::ETransp transa, Teuchos::EUplo uplo) const
void solver_CLP(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &A, const Teuchos::SerialDenseVector< ordinal_type, value_type > &b, Teuchos::SerialDenseVector< ordinal_type, value_type > &x, Teuchos::ETransp transa, Teuchos::EUplo uplo) const
void solver_qpOASES(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &A, const Teuchos::SerialDenseVector< ordinal_type, value_type > &b, Teuchos::SerialDenseVector< ordinal_type, value_type > &x, Teuchos::ETransp transa, Teuchos::EUplo uplo) const
scalar_type vec_norm_inf(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A)
Vector-infinity norm of a matrix.
void CPQR_Householder3(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R, Teuchos::Array< ordinal_type > &piv)
Compute column-pivoted QR using Householder reflections.
ordinal_type n_choose_k(const ordinal_type &n, const ordinal_type &k)
Compute bionomial coefficient (n ; k) = n!/( k! (n-k)! )