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