42#include "Teuchos_TestForException.hpp"
43#include "Teuchos_SerialDenseVector.hpp"
49#ifdef HAVE_STOKHOS_GLPK
55#ifdef HAVE_STOKHOS_CLP
56#include "coin/ClpSimplex.hpp"
57#include "coin/CoinBuild.hpp"
58#include "coin/ClpInterior.hpp"
61#ifdef HAVE_STOKHOS_QPOASES
65#ifdef HAVE_STOKHOS_DAKOTA
66#include "LinearSolver.hpp"
69template <
typename ordinal_type,
typename value_type>
72 const Teuchos::ParameterList& 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))
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
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;
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);
102 reducedQuadrature_Q_Squared(Q, F, weights,
103 red_weights, red_points, red_values);
105 else if (reduction_method ==
"Q Squared2") {
106 reducedQuadrature_Q_Squared_CPQR2(Q, F, weights,
107 red_weights, red_points, red_values);
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);
114 reducedQuadrature_Q2(Q,
Q2, F, weights,
115 red_weights, red_points, red_values);
117 else if (reduction_method ==
"None") {
118 ordinal_type sz = Q.numCols();
119 ordinal_type nqp = Q.numRows();
120 ordinal_type d = F.numCols();
122 Teuchos::rcp(
new Teuchos::Array<value_type>(weights));
124 Teuchos::rcp(
new Teuchos::Array< Teuchos::Array<value_type> >(nqp));
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);
137 TEUCHOS_TEST_FOR_EXCEPTION(
138 true, std::logic_error,
139 "Invalid dimension reduction method " << reduction_method);
142 Teuchos::RCP< const Teuchos::Array<value_type> > cred_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;
147 Teuchos::RCP<const Stokhos::UserDefinedQuadrature<ordinal_type, value_type> >
149 Teuchos::rcp(
new UserDefinedQuadrature<ordinal_type,value_type>(
157template <
typename ordinal_type,
typename value_type>
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
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();
175 Teuchos::SerialDenseMatrix<ordinal_type, value_type>
Q2(nqp, sz2);
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);
184 TEUCHOS_ASSERT(jdx == sz2);
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);
191 e1.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0,
Q2, w, 0.0);
192 TEUCHOS_ASSERT(ret == 0);
195 Teuchos::SerialDenseVector<ordinal_type,value_type> u(nqp);
196 underdetermined_solver(
Q2, e1, u, Teuchos::TRANS, Teuchos::UNDEF_TRI);
199 std::cout <<
"sz = " << sz << std::endl;
200 std::cout <<
"nqp = " << nqp << std::endl;
201 std::cout <<
"sz2 = " << sz2 << std::endl;
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;
212 Teuchos::SerialDenseMatrix<ordinal_type,value_type> err2(sz, sz);
214 for (ordinal_type i=0; i<sz; i++)
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)
227 ordinal_type rank = 0;
228 for (ordinal_type i=0; i<nqp; i++)
229 if (std::abs(u[i]) > reduction_tol) ++rank;
233 Teuchos::rcp(
new Teuchos::Array<value_type>(rank));
235 Teuchos::rcp(
new Teuchos::Array< Teuchos::Array<value_type> >(rank));
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);
254 TEUCHOS_ASSERT(idx <= rank);
257 reduced_weights->resize(rank);
258 reduced_points->resize(rank);
259 reduced_values->resize(rank);
263template <
typename ordinal_type,
typename value_type>
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
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();
281 Teuchos::SerialDenseMatrix<ordinal_type, value_type>
Q2(nqp, sz2);
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);
290 TEUCHOS_ASSERT(jdx == sz2);
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);
297 e1.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0,
Q2, w, 0.0);
298 TEUCHOS_ASSERT(ret == 0);
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;
305 ordinal_type r = computeRank(R, reduction_tol);
306 R.reshape(r, R.numCols());
307 Z.reshape(Z.numRows(), r);
310 std::cout <<
"Q2 rank = " << r << std::endl;
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);
319 Teuchos::SerialDenseVector<ordinal_type,value_type> ut(nqp);
320 underdetermined_solver(R, e1t, ut, Teuchos::NO_TRANS, Teuchos::UPPER_TRI);
323 Teuchos::SerialDenseVector<ordinal_type,value_type> u(nqp);
324 for (ordinal_type i=0; i<nqp; i++)
328 std::cout <<
"sz = " << sz << std::endl;
329 std::cout <<
"nqp = " << nqp << std::endl;
330 std::cout <<
"sz2 = " << sz2 << std::endl;
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;
341 Teuchos::SerialDenseMatrix<ordinal_type,value_type> err2(sz, sz);
343 for (ordinal_type i=0; i<sz; i++)
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)
356 ordinal_type rank = 0;
357 for (ordinal_type i=0; i<nqp; i++)
358 if (std::abs(u[i]) > reduction_tol) ++rank;
362 Teuchos::rcp(
new Teuchos::Array<value_type>(rank));
364 Teuchos::rcp(
new Teuchos::Array< Teuchos::Array<value_type> >(rank));
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);
383 TEUCHOS_ASSERT(idx <= rank);
386 reduced_weights->resize(rank);
387 reduced_points->resize(rank);
388 reduced_values->resize(rank);
392template <
typename ordinal_type,
typename value_type>
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
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();
410 Teuchos::SerialDenseMatrix<ordinal_type, value_type>
Q2(nqp, sz2);
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);
419 TEUCHOS_ASSERT(jdx == sz2);
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,
433 bool restrict_r = params.get(
"Restrict Rank",
false);
435 ordinal_type d = F.numCols();
436 ordinal_type p = params.get<ordinal_type>(
"Order Restriction");
440 std::cout <<
"Restricting rank from " << r <<
" to " << n << std::endl;
446 bool use_Z = params.get(
"Use Q in LP",
true);
447 Teuchos::SerialDenseMatrix<ordinal_type, value_type> QQ2(nqp, r);
449 for (ordinal_type
j=0;
j<r;
j++)
450 for (ordinal_type i=0; i<nqp; i++)
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]);
460 std::cout <<
"Q2 rank = " << r << std::endl;
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);
468 ee1.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, QQ2, w, 0.0);
469 TEUCHOS_ASSERT(ret == 0);
475 Teuchos::SerialDenseVector<ordinal_type,value_type> u(nqp);
476 underdetermined_solver(QQ2, ee1, u, Teuchos::TRANS, Teuchos::UNDEF_TRI);
479 std::cout <<
"sz = " << sz << std::endl;
480 std::cout <<
"nqp = " << nqp << std::endl;
481 std::cout <<
"sz2 = " << sz2 << std::endl;
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;
494 Teuchos::SerialDenseMatrix<ordinal_type,value_type> err2(sz, sz);
496 for (ordinal_type i=0; i<sz; i++)
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)
509 ordinal_type rank = 0;
510 for (ordinal_type i=0; i<nqp; i++)
511 if (std::abs(u[i]) > reduction_tol) ++rank;
515 Teuchos::rcp(
new Teuchos::Array<value_type>(rank));
517 Teuchos::rcp(
new Teuchos::Array< Teuchos::Array<value_type> >(rank));
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);
536 TEUCHOS_ASSERT(idx <= rank);
539 reduced_weights->resize(rank);
540 reduced_points->resize(rank);
541 reduced_values->resize(rank);
545template <
typename ordinal_type,
typename value_type>
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
559 ordinal_type sz = Q.numCols();
560 ordinal_type nqp = Q.numRows();
561 ordinal_type sz2 =
Q2.numCols();
562 ordinal_type d = F.numCols();
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);
569 e1.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0,
Q2, w, 0.0);
570 TEUCHOS_ASSERT(ret == 0);
573 Teuchos::SerialDenseVector<ordinal_type,value_type> u(nqp);
574 underdetermined_solver(
Q2, e1, u, Teuchos::TRANS, Teuchos::UNDEF_TRI);
577 std::cout <<
"sz = " << sz << std::endl;
578 std::cout <<
"nqp = " << nqp << std::endl;
579 std::cout <<
"sz2 = " << sz2 << std::endl;
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;
590 Teuchos::SerialDenseMatrix<ordinal_type,value_type> err2(sz, sz);
592 for (ordinal_type i=0; i<sz; i++)
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)
604 ordinal_type rank = 0;
605 for (ordinal_type i=0; i<nqp; i++)
606 if (std::abs(u[i]) > reduction_tol) ++rank;
610 Teuchos::rcp(
new Teuchos::Array<value_type>(rank));
612 Teuchos::rcp(
new Teuchos::Array< Teuchos::Array<value_type> >(rank));
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);
631 TEUCHOS_ASSERT(idx <= rank);
634 reduced_weights->resize(rank);
635 reduced_points->resize(rank);
636 reduced_values->resize(rank);
640template <
typename ordinal_type,
typename value_type>
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
654 ordinal_type sz = Q.numCols();
655 ordinal_type nqp = Q.numRows();
656 ordinal_type sz2 =
Q2.numCols();
657 ordinal_type d = F.numCols();
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);
664 e1.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0,
Q2, w, 0.0);
665 TEUCHOS_ASSERT(ret == 0);
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;
672 ordinal_type r = computeRank(R, reduction_tol);
673 R.reshape(r, R.numCols());
674 Z.reshape(Z.numRows(), r);
677 std::cout <<
"Q2 rank = " << r << std::endl;
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);
686 Teuchos::SerialDenseVector<ordinal_type,value_type> ut(nqp);
687 underdetermined_solver(R, e1t, ut, Teuchos::NO_TRANS, Teuchos::UPPER_TRI);
690 Teuchos::SerialDenseVector<ordinal_type,value_type> u(nqp);
691 for (ordinal_type i=0; i<nqp; i++)
695 std::cout <<
"sz = " << sz << std::endl;
696 std::cout <<
"nqp = " << nqp << std::endl;
697 std::cout <<
"sz2 = " << sz2 << std::endl;
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;
708 Teuchos::SerialDenseMatrix<ordinal_type,value_type> err2(sz, sz);
710 for (ordinal_type i=0; i<sz; i++)
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)
723 ordinal_type rank = 0;
724 for (ordinal_type i=0; i<nqp; i++)
725 if (std::abs(u[i]) > reduction_tol) ++rank;
729 Teuchos::rcp(
new Teuchos::Array<value_type>(rank));
731 Teuchos::rcp(
new Teuchos::Array< Teuchos::Array<value_type> >(rank));
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);
750 TEUCHOS_ASSERT(idx <= rank);
753 reduced_weights->resize(rank);
754 reduced_points->resize(rank);
755 reduced_values->resize(rank);
759template <
typename ordinal_type,
typename value_type>
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
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);
781 TEUCHOS_TEST_FOR_EXCEPTION(
782 true, std::logic_error,
783 "Invalid solver method " << solver_method);
786template <
typename ordinal_type,
typename value_type>
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
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();
801 if (transa == Teuchos::NO_TRANS) {
809 if (x.length() < n_cols)
812 blas.COPY(n_rows, b.values(), 1, x.values(), 1);
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());
820template <
typename ordinal_type,
typename value_type>
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
829#ifdef HAVE_STOKHOS_GLPK
830 ordinal_type m = A.numRows();
831 ordinal_type n = A.numCols();
834 if (transa == Teuchos::NO_TRANS) {
842 if (x.length() < n_cols)
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);
852 for (ordinal_type i=0; i<n_rows; i++)
853 lpx_set_row_bnds(lp, i+1, LPX_FX, b[i], b[i]);
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);
865 Teuchos::SerialDenseMatrix<ordinal_type, value_type> AA(n_rows+1,n_cols+1);
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);
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;
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]);
891 else if (AA_uplo == Teuchos::UPPER_TRI) {
894 for (ordinal_type
j=0;
j<n_cols;
j++) {
895 ordinal_type nr = n_rows;
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]);
904 else if (AA_uplo == Teuchos::LOWER_TRI) {
907 for (ordinal_type
j=0;
j<n_cols;
j++) {
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);
919 if (params.get(
"Write MPS File",
false) ==
true)
920 lpx_write_mps(lp, params.get(
"MPS Filename",
"lp.mps").c_str());
924 int status = lpx_get_status(lp);
926 std::cout <<
"glpk status = " << status << std::endl;
927 double z = lpx_get_obj_val(lp);
928 std::cout <<
"glpk objective = " << z << std::endl;
932 for (ordinal_type i=0; i<n_cols; i++)
933 x[i] = lpx_get_col_prim(lp, i+1);
935 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
936 "GLPK solver called but not enabled!");
940template <
typename ordinal_type,
typename value_type>
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
949#ifdef HAVE_STOKHOS_CLP
950 ordinal_type m = A.numRows();
951 ordinal_type n = A.numCols();
954 if (transa == Teuchos::NO_TRANS) {
962 if (x.length() < n_cols)
967 model.resize(n_rows, 0);
970 for (ordinal_type i=0; i<n_rows; i++) {
971 model.setRowLower(i, b[i]);
972 model.setRowUpper(i, b[i]);
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);
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++)
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;
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);
1003 else if (AA_uplo == Teuchos::UPPER_TRI) {
1006 for (ordinal_type
j=0;
j<n_cols;
j++) {
1007 ordinal_type nr = n_rows;
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);
1016 else if (AA_uplo == Teuchos::LOWER_TRI) {
1019 for (ordinal_type
j=0;
j<n_cols;
j++) {
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);
1029 buildObject.addColumn(0, NULL, NULL, 0.0, DBL_MAX, objective_value);
1032 model.addColumns(buildObject);
1038 double *solution = model.primalColumnSolution();
1039 for (ordinal_type i=0; i<n_cols; i++)
1042 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
1043 "CLP solver called but not enabled!");
1047template <
typename ordinal_type,
typename value_type>
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
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) {
1069 if (x.length() < n_cols)
1074 model.resize(n_rows, 0);
1077 for (ordinal_type i=0; i<n_rows; i++) {
1078 model.setRowLower(i, b[i]);
1079 model.setRowUpper(i, b[i]);
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);
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++)
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;
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);
1110 else if (AA_uplo == Teuchos::UPPER_TRI) {
1113 for (ordinal_type
j=0;
j<n_cols;
j++) {
1114 ordinal_type nr = n_rows;
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);
1123 else if (AA_uplo == Teuchos::LOWER_TRI) {
1126 for (ordinal_type
j=0;
j<n_cols;
j++) {
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);
1136 buildObject.addColumn(0, NULL, NULL, 0.0, DBL_MAX, objective_value);
1139 model.addColumns(buildObject);
1145 double *solution = model.primalColumnSolution();
1146 for (ordinal_type i=0; i<n_cols; i++)
1149 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
1150 "CLP solver called but not enabled!");
1154template <
typename ordinal_type,
typename value_type>
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
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) {
1176 if (x.length() < n_cols)
1180 Teuchos::SerialDenseVector<ordinal_type,value_type> c(n_cols);
1181 c.putScalar(objective_value);
1184 Teuchos::SerialDenseVector<ordinal_type,value_type> lb(n_cols);
1188 qpOASES::QProblem lp(n_cols, n_rows, qpOASES::HST_ZERO);
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);
1221 lp.getPrimalSolution(x.values());
1223 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
1224 "qpOASES solver called but not enabled!");
1228template <
typename ordinal_type,
typename value_type>
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
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) {
1247 if (x.length() < n_cols)
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);
1262 TEUCHOS_TEST_FOR_EXCEPTION(
1263 true, std::logic_error,
1264 "CompressedSensing solver called but not enabled!");
1268template <
typename ordinal_type,
typename value_type>
1272 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& R,
1273 const value_type tol)
const
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)
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;
1302template <
typename ordinal_type,
typename value_type>
1305n_choose_k(
const ordinal_type& n,
const ordinal_type& k)
const {
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++) {
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 ¶ms)
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)! )