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