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 62 #include "qpOASES.hpp" 65 #ifdef HAVE_STOKHOS_DAKOTA 66 #include "LinearSolver.hpp" 69 template <
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))
83 template <
typename ordinal_type,
typename value_type>
84 Teuchos::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") {
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));
128 (*red_points)[i].resize(d);
130 (*red_points)[i][
j] = F(i,
j);
131 (*red_values)[i].resize(sz);
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> >
157 template <
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
175 Teuchos::SerialDenseMatrix<ordinal_type, value_type>
Q2(nqp, sz2);
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);
216 Teuchos::SerialDenseMatrix<ordinal_type,value_type> WQ(nqp, sz);
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)
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));
240 if (
std::abs(u[i]) > reduction_tol) {
241 (*reduced_weights)[idx] = u[i];
242 (*reduced_points)[idx].resize(d);
244 (*reduced_points)[idx][
j] = F(i,
j);
245 (*reduced_values)[idx].resize(sz);
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);
263 template <
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
281 Teuchos::SerialDenseMatrix<ordinal_type, value_type>
Q2(nqp, sz2);
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;
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);
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);
345 Teuchos::SerialDenseMatrix<ordinal_type,value_type> WQ(nqp, sz);
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)
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));
369 if (
std::abs(u[i]) > reduction_tol) {
370 (*reduced_weights)[idx] = u[i];
371 (*reduced_points)[idx].resize(d);
373 (*reduced_points)[idx][
j] = F(i,
j);
374 (*reduced_values)[idx].resize(sz);
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);
392 template <
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
410 Teuchos::SerialDenseMatrix<ordinal_type, value_type>
Q2(nqp, sz2);
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")
431 orthogonalization_method, reduction_tol, verbose,
Q2, ww,
433 bool restrict_r = params.get(
"Restrict Rank",
false);
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);
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);
498 Teuchos::SerialDenseMatrix<ordinal_type,value_type> WQ(nqp, sz);
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)
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));
522 if (
std::abs(u[i]) > reduction_tol) {
523 (*reduced_weights)[idx] = u[i];
524 (*reduced_points)[idx].resize(d);
526 (*reduced_points)[idx][
j] = F(i,
j);
527 (*reduced_values)[idx].resize(sz);
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);
545 template <
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
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);
594 Teuchos::SerialDenseMatrix<ordinal_type,value_type> WQ(nqp, sz);
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)
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));
617 if (
std::abs(u[i]) > reduction_tol) {
618 (*reduced_weights)[idx] = u[i];
619 (*reduced_points)[idx].resize(d);
621 (*reduced_points)[idx][
j] = F(i,
j);
622 (*reduced_values)[idx].resize(sz);
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);
640 template <
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
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;
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);
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);
712 Teuchos::SerialDenseMatrix<ordinal_type,value_type> WQ(nqp, sz);
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)
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));
736 if (
std::abs(u[i]) > reduction_tol) {
737 (*reduced_weights)[idx] = u[i];
738 (*reduced_points)[idx].resize(d);
740 (*reduced_points)[idx][
j] = F(i,
j);
741 (*reduced_values)[idx].resize(sz);
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);
759 template <
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);
786 template <
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!");
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());
820 template <
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 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);
853 lpx_set_row_bnds(lp, i+1, LPX_FX, b[i], b[i]);
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);
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) {
885 row_indices[
j].resize(n_rows+1);
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) {
898 row_indices[
j].resize(nr+1);
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) {
910 row_indices[
j].resize(nr+1);
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;
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!");
940 template <
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 954 if (transa == Teuchos::NO_TRANS) {
962 if (
x.length() < n_cols)
967 model.resize(n_rows, 0);
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);
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) {
996 row_indices[
j].resize(n_rows);
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) {
1010 row_indices[
j].resize(nr);
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) {
1022 row_indices[
j].resize(nr);
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();
1042 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
1043 "CLP solver called but not enabled!");
1047 template <
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 1061 if (transa == Teuchos::NO_TRANS) {
1069 if (
x.length() < n_cols)
1074 model.resize(n_rows, 0);
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);
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) {
1103 row_indices[
j].resize(n_rows);
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) {
1117 row_indices[
j].resize(nr);
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) {
1129 row_indices[
j].resize(nr);
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();
1149 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
1150 "CLP solver called but not enabled!");
1154 template <
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 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!");
1228 template <
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 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!");
1268 template <
typename ordinal_type,
typename value_type>
1272 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& R,
1282 for (rank=1; rank<m; rank++) {
1283 if (
std::abs(R(rank,rank)) > r_max)
1285 if (
std::abs(R(rank,rank)) < r_min)
1287 if (r_min / r_max < tol)
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;
1302 template <
typename ordinal_type,
typename value_type>
ordinal_type n_choose_k(const ordinal_type &n, const ordinal_type &k)
Compute bionomial coefficient (n ; k) = n!/( k! (n-k)! )
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.
KOKKOS_INLINE_FUNCTION PCE< Storage > min(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
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
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
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
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_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
ReducedQuadratureFactory(const Teuchos::ParameterList ¶ms)
Constructor.
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
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
scalar_type vec_norm_inf(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A)
Vector-infinity norm of a matrix.
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
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_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
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
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
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 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
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_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
Encapsulate various orthogonalization (ie QR) methods.