Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Stokhos_SDMUtils.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 #ifndef STOKHOS_SDM_UTILS_HPP
11 #define STOKHOS_SDM_UTILS_HPP
12 
13 #include "Teuchos_SerialDenseMatrix.hpp"
14 #include "Teuchos_SerialDenseVector.hpp"
15 #include "Teuchos_SerialDenseHelpers.hpp"
16 #include "Teuchos_Array.hpp"
17 #include "Teuchos_LAPACK.hpp"
18 #include <ostream>
19 
20 #define DGEQP3_F77 F77_BLAS_MANGLE(dgeqp3,DGEQP3)
21 
22 extern "C" {
23 #if defined(HAVE_STOKHOS_MKL)
24  #include "mkl_version.h"
25  #if __INTEL_MKL__ >= 2021
26  #define MKL_NO_EXCEPT noexcept
27  #else
28  #define MKL_NO_EXCEPT
29  #endif
30 #else
31  #define MKL_NO_EXCEPT
32 #endif
33 
34 void DGEQP3_F77(const int*, const int*, double*, const int*, int*,
35  double*, double*, const int*, int*) MKL_NO_EXCEPT;
36 }
37 
38 #include "Stokhos_ConfigDefs.h"
39 
40 #ifdef HAVE_STOKHOS_MATLABLIB
41 extern "C" {
42 #include "mat.h"
43 #include "matrix.h"
44 }
45 #endif
46 
47 namespace Stokhos {
48 
49  namespace detail {
50 
52  template <typename ordinal_type, typename scalar_type>
53  scalar_type weighted_inner_product(
54  const Teuchos::SerialDenseVector<ordinal_type,scalar_type>& x,
55  const Teuchos::SerialDenseVector<ordinal_type,scalar_type>& y,
56  const Teuchos::Array<scalar_type>& w)
57  {
58  ordinal_type n = x.length();
59  scalar_type t = 0;
60  for (ordinal_type i=0; i<n; i++)
61  t += x[i]*w[i]*y[i];
62  return t;
63  }
64 
66  template <typename ordinal_type, typename scalar_type>
67  void saxpy(
68  const scalar_type& alpha,
69  Teuchos::SerialDenseVector<ordinal_type, scalar_type>& x,
70  const scalar_type& beta,
71  const Teuchos::SerialDenseVector<ordinal_type, scalar_type>& y)
72  {
73  ordinal_type n = x.length();
74  for (ordinal_type i=0; i<n; i++)
75  x[i] = alpha*x[i] + beta*y[i];
76  }
77 
78  }
79 
80  // Print a matrix in matlab-esque format
81  template <typename ordinal_type, typename scalar_type>
82  void
83  print_matlab(std::ostream& os,
84  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A)
85  {
86  os << "[ ";
87  for (ordinal_type i=0; i<A.numRows(); i++) {
88  for (ordinal_type j=0; j<A.numCols(); j++)
89  os << A(i,j) << " ";
90  if (i < A.numRows()-1)
91  os << ";" << std::endl << " ";
92  }
93  os << "];" << std::endl;
94  }
95 
96 #ifdef HAVE_STOKHOS_MATLABLIB
97  // Save a matrix to matlab MAT file
98  template <typename ordinal_type>
99  void
100  save_matlab(const std::string& file_name, const std::string& matrix_name,
101  const Teuchos::SerialDenseMatrix<ordinal_type,double>& A,
102  bool overwrite = true)
103  {
104  // Open matfile
105  const char *mode;
106  if (overwrite)
107  mode = "w";
108  else
109  mode = "u";
110  MATFile *file = matOpen(file_name.c_str(), mode);
111  TEUCHOS_ASSERT(file != NULL);
112 
113  // Convert SDM to Matlab matrix
114  mxArray *mat = mxCreateDoubleMatrix(0, 0, mxREAL);
115  TEUCHOS_ASSERT(mat != NULL);
116  mxSetPr(mat, A.values());
117  mxSetM(mat, A.numRows());
118  mxSetN(mat, A.numCols());
119 
120  // Save matrix to file
121  ordinal_type ret = matPutVariable(file, matrix_name.c_str(), mat);
122  TEUCHOS_ASSERT(ret == 0);
123 
124  // Close file
125  ret = matClose(file);
126  TEUCHOS_ASSERT(ret == 0);
127 
128  // Free matlab matrix
129  mxSetPr(mat, NULL);
130  mxSetM(mat, 0);
131  mxSetN(mat, 0);
132  mxDestroyArray(mat);
133  }
134 #endif
135 
137  template <typename ordinal_type, typename scalar_type>
138  scalar_type vec_norm_inf(
139  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A) {
140  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> vec_A(
141  Teuchos::View, A.values(), 1, A.numRows()*A.numCols(), 1);
142  return vec_A.normInf();
143  }
144 
146 
150  template <typename ordinal_type, typename scalar_type>
151  void QR_CGS(
152  ordinal_type k,
153  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
154  const Teuchos::Array<scalar_type>& w,
155  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
156  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R)
157  {
158  using Teuchos::getCol;
159  typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type> SDV;
160  typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> SDM;
161 
162  // getCol() requires non-const SDM
163  SDM& Anc = const_cast<SDM&>(A);
164 
165  // Make sure Q is the right size
166  ordinal_type m = A.numRows();
167  if (Q.numRows() != m || Q.numCols() != k)
168  Q.shape(m,k);
169  if (R.numRows() != k || R.numCols() != k)
170  R.shape(k,k);
171 
172  for (ordinal_type j=0; j<k; j++) {
173  SDV Aj = getCol(Teuchos::View, Anc, j);
174  SDV Qj = getCol(Teuchos::View, Q, j);
175  Qj.assign(Aj);
176  for (ordinal_type i=0; i<j; i++) {
177  SDV Qi = getCol(Teuchos::View, Q, i);
178  R(i,j) = detail::weighted_inner_product(Qi, Aj, w);
179  detail::saxpy(1.0, Qj, -R(i,j), Qi); // Q(:,j) = 1.0*Q(:,j) - R(i,j)*Q(:,i)
180  }
181  R(j,j) = std::sqrt(detail::weighted_inner_product(Qj, Qj, w));
182  Qj.scale(1.0/R(j,j));
183  }
184 
185  }
186 
188 
192  template <typename ordinal_type, typename scalar_type>
193  void QR_MGS(
194  ordinal_type k,
195  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
196  const Teuchos::Array<scalar_type>& w,
197  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
198  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R)
199  {
200  using Teuchos::getCol;
201  typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type> SDV;
202  typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> SDM;
203 
204  // getCol() requires non-const SDM
205  SDM& Anc = const_cast<SDM&>(A);
206 
207  // Make sure Q is the right size
208  ordinal_type m = A.numRows();
209  if (Q.numRows() != m || Q.numCols() != k)
210  Q.shape(m,k);
211  if (R.numRows() != k || R.numCols() != k)
212  R.shape(k,k);
213 
214  for (ordinal_type i=0; i<k; i++) {
215  SDV Ai = getCol(Teuchos::View, Anc, i);
216  SDV Qi = getCol(Teuchos::View, Q, i);
217  Qi.assign(Ai);
218  }
219  for (ordinal_type i=0; i<k; i++) {
220  SDV Qi = getCol(Teuchos::View, Q, i);
221  for (ordinal_type j=0; j<i; j++) {
222  SDV Qj = getCol(Teuchos::View, Q, j);
223  scalar_type v = detail::weighted_inner_product(Qi, Qj, w);
224  detail::saxpy(1.0, Qi, -v, Qj); // Q(:,i) = 1.0*Q(:,i) - v*Q(:,j)
225  R(j,i) += v;
226  }
227  R(i,i) = std::sqrt(detail::weighted_inner_product(Qi, Qi, w));
228  Qi.scale(1.0/R(i,i));
229  }
230  }
231 
233 
237  template <typename ordinal_type, typename scalar_type>
238  void QR_MGS2(
239  ordinal_type k,
240  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
241  const Teuchos::Array<scalar_type>& w,
242  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
243  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R)
244  {
245  using Teuchos::getCol;
246  typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type> SDV;
247  typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> SDM;
248 
249  // getCol() requires non-const SDM
250  SDM& Anc = const_cast<SDM&>(A);
251 
252  // Make sure Q is the right size
253  ordinal_type m = A.numRows();
254  if (Q.numRows() != m || Q.numCols() != k)
255  Q.shape(m,k);
256  if (R.numRows() != k || R.numCols() != k)
257  R.shape(k,k);
258 
259  for (ordinal_type i=0; i<k; i++) {
260  SDV Ai = getCol(Teuchos::View, Anc, i);
261  SDV Qi = getCol(Teuchos::View, Q, i);
262  Qi.assign(Ai);
263  }
264  for (ordinal_type i=0; i<k; i++) {
265  SDV Qi = getCol(Teuchos::View, Q, i);
266  for (ordinal_type j=0; j<i; j++) {
267  SDV Qj = getCol(Teuchos::View, Q, j);
268  scalar_type v = detail::weighted_inner_product(Qi, Qj, w);
269  detail::saxpy(1.0, Qi, -v, Qj); // Q(:,i) = 1.0*Q(:,i) - v*Q(:,j)
270  R(j,i) += v;
271  }
272  for (ordinal_type j=i-1; j>=0; j--) {
273  SDV Qj = getCol(Teuchos::View, Q, j);
274  scalar_type v = detail::weighted_inner_product(Qi, Qj, w);
275  detail::saxpy(1.0, Qi, -v, Qj); // Q(:,i) = 1.0*Q(:,i) - v*Q(:,j)
276  R(j,i) += v;
277  }
278  R(i,i) = std::sqrt(detail::weighted_inner_product(Qi, Qi, w));
279  Qi.scale(1.0/R(i,i));
280  }
281  }
282 
284 
290  template <typename ordinal_type, typename scalar_type>
291  void
292  QR_Householder(
293  ordinal_type k,
294  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
295  const Teuchos::Array<scalar_type>& w,
296  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
297  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R)
298  {
299  Teuchos::LAPACK<ordinal_type,scalar_type> lapack;
300  ordinal_type m = A.numRows();
301  ordinal_type n = A.numCols();
302  ordinal_type kk = std::min(m,n);
303  if (k > kk)
304  k = kk;
305 
306  // Check that each component of w is 1
307  for (ordinal_type i=0; i<w.size(); i++)
308  TEUCHOS_TEST_FOR_EXCEPTION(
309  w[i] != 1.0, std::logic_error,
310  "CPQR_Householder_threshold() requires unit weight vector!");
311 
312  // Lapack routine overwrites A, so copy into temporary matrix
313  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
314  Teuchos::Copy, A, m, n);
315 
316  // QR
317  ordinal_type lda = AA.stride();
318  Teuchos::Array<scalar_type> tau(k);
319  Teuchos::Array<scalar_type> work(1);
320  ordinal_type info;
321  ordinal_type lwork = -1;
322  lapack.GEQRF(m, n, AA.values(), lda, &tau[0], &work[0], lwork, &info);
323  TEUCHOS_TEST_FOR_EXCEPTION(
324  info < 0, std::logic_error, "geqrf returned info = " << info);
325  lwork = work[0];
326  work.resize(lwork);
327  lapack.GEQRF(m, n, AA.values(), lda, &tau[0], &work[0], lwork, &info);
328  TEUCHOS_TEST_FOR_EXCEPTION(
329  info < 0, std::logic_error, "geqrf returned info = " << info);
330 
331  // Extract R
332  if (R.numRows() != k || R.numCols() != n)
333  R.shape(k,n);
334  R.putScalar(0.0);
335  for (ordinal_type i=0; i<k; i++)
336  for (ordinal_type j=i; j<n; j++)
337  R(i,j) = AA(i,j);
338 
339  // Extract Q
340  if (Q.numRows() != m || Q.numCols() != k)
341  Q.shape(m,k);
342  lwork = -1;
343  lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
344  TEUCHOS_TEST_FOR_EXCEPTION(
345  info < 0, std::logic_error, "orgqr returned info = " << info);
346  lwork = work[0];
347  work.resize(lwork);
348  lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
349  TEUCHOS_TEST_FOR_EXCEPTION(
350  info < 0, std::logic_error, "orgqr returned info = " << info);
351  if (Q.numRows() != m || Q.numCols() != k)
352  Q.shape(m, k);
353  for (ordinal_type i=0; i<m; i++)
354  for (ordinal_type j=0; j<k; j++)
355  Q(i,j) = AA(i,j);
356  }
357 
359 
371  template <typename ordinal_type, typename scalar_type>
372  void
373  CPQR_Householder3(
374  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
375  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
376  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
377  Teuchos::Array<ordinal_type>& piv)
378  {
379  Teuchos::LAPACK<ordinal_type,scalar_type> lapack;
380  ordinal_type m = A.numRows();
381  ordinal_type n = A.numCols();
382  ordinal_type k = std::min(m,n);
383 
384  // Lapack routine overwrites A, so copy into temporary matrix
385  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
386  Teuchos::Copy, A, m, n);
387  if (Q.numRows() != m || Q.numCols() != k)
388  Q.shape(m,k);
389 
390  // Teuchos LAPACK interface doesn't have dgeqp3, so call it directly
391 
392  // Workspace query
393  ordinal_type lda = AA.stride();
394  piv.resize(n);
395  Teuchos::Array<scalar_type> tau(k);
396  Teuchos::Array<scalar_type> work(1);
397  ordinal_type lwork = -1;
398  ordinal_type info;
399  DGEQP3_F77(&m, &n, AA.values(), &lda, &piv[0], &tau[0], &work[0], &lwork,
400  &info);
401  TEUCHOS_TEST_FOR_EXCEPTION(
402  info < 0, std::logic_error, "dgeqp3 returned info = " << info);
403 
404  // Column pivoted QR
405  lwork = work[0];
406  work.resize(lwork);
407  DGEQP3_F77(&m, &n, AA.values(), &lda, &piv[0], &tau[0], &work[0], &lwork,
408  &info);
409  TEUCHOS_TEST_FOR_EXCEPTION(
410  info < 0, std::logic_error, "dgeqp3 returned info = " << info);
411 
412  // Extract R
413  if (R.numRows() != k || R.numCols() != n)
414  R.shape(k,n);
415  R.putScalar(0.0);
416  for (ordinal_type i=0; i<k; i++)
417  for (ordinal_type j=i; j<n; j++)
418  R(i,j) = AA(i,j);
419 
420  // Extract Q
421  lwork = -1;
422  lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
423  TEUCHOS_TEST_FOR_EXCEPTION(
424  info < 0, std::logic_error, "orgqr returned info = " << info);
425  lwork = work[0];
426  work.resize(lwork);
427  lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
428  TEUCHOS_TEST_FOR_EXCEPTION(
429  info < 0, std::logic_error, "orgqr returned info = " << info);
430  if (Q.numRows() != m || Q.numCols() != k)
431  Q.shape(m, k);
432  for (ordinal_type i=0; i<m; i++)
433  for (ordinal_type j=0; j<k; j++)
434  Q(i,j) = AA(i,j);
435 
436  // Transform piv to zero-based indexing
437  for (ordinal_type i=0; i<n; i++)
438  piv[i] -= 1;
439  }
440 
442 
460  template <typename ordinal_type, typename scalar_type>
461  ordinal_type
462  CPQR_Householder_threshold(
463  const scalar_type& rank_threshold,
464  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
465  const Teuchos::Array<scalar_type>& w,
466  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
467  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
468  Teuchos::Array<ordinal_type>& piv)
469  {
470  // Check that each component of w is 1
471  for (ordinal_type i=0; i<w.size(); i++)
472  TEUCHOS_TEST_FOR_EXCEPTION(
473  w[i] != 1.0, std::logic_error,
474  "CPQR_Householder_threshold() requires unit weight vector!");
475 
476  // Compute full QR
477  CPQR_Householder3(A, Q, R, piv);
478 
479  // Find leading block of R such that cond(R) <= 1/tau
480  ordinal_type rank = 0;
481  ordinal_type m = R.numRows();
482  scalar_type r_max = std::abs(R(rank,rank));
483  scalar_type r_min = std::abs(R(rank,rank));
484  for (rank=1; rank<m; rank++) {
485  if (std::abs(R(rank,rank)) > r_max)
486  r_max = std::abs(R(rank,rank));
487  if (std::abs(R(rank,rank)) < r_min)
488  r_min = std::abs(R(rank,rank));
489  if (r_min / r_max < rank_threshold)
490  break;
491  }
492 
493  // Extract blocks from Q and R
494  R.reshape(rank,rank);
495  Q.reshape(Q.numRows(), rank);
496 
497  return rank;
498  }
499 
501 
512  template <typename ordinal_type, typename scalar_type>
513  ordinal_type
514  CPQR_MGS_threshold(
515  const scalar_type& rank_threshold,
516  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
517  const Teuchos::Array<scalar_type>& w,
518  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
519  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
520  Teuchos::Array<ordinal_type>& piv)
521  {
522  using Teuchos::getCol;
523  typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type> SDV;
524  ordinal_type m = A.numRows();
525  ordinal_type n = A.numCols();
526  ordinal_type k = std::min(m,n);
527 
528  // Make sure Q is the right size
529  if (Q.numRows() != m || Q.numCols() != n)
530  Q.shape(m,n);
531  if (R.numRows() != m || R.numCols() != n)
532  R.shape(m,n);
533  if (piv.size() != n)
534  piv.resize(n);
535  Q.assign(A);
536 
537  // Compute column norms
538  SDV nrm(n);
539  for (ordinal_type j=0; j<n; j++) {
540  SDV Qj = getCol(Teuchos::View, Q, j);
541  nrm(j) = detail::weighted_inner_product(Qj, Qj, w);
542  }
543  SDV Qtmp(m), Rtmp(m);
544 
545  Teuchos::Array<ordinal_type> piv_orig(piv);
546  for (ordinal_type i=0; i<n; i++)
547  piv[i] = i;
548 
549  // Prepivot any columns requested by setting piv[i] != 0
550  ordinal_type nfixed = 0;
551  for (ordinal_type j=0; j<n; j++) {
552  if (piv_orig[j] != 0) {
553  if (j != nfixed) {
554  scalar_type tmp = nrm(j);
555  nrm(j) = nrm(nfixed);
556  nrm(nfixed) = tmp;
557 
558  SDV Qpvt = getCol(Teuchos::View, Q, j);
559  SDV Qj = getCol(Teuchos::View, Q, nfixed);
560  Qtmp.assign(Qpvt);
561  Qpvt.assign(Qj);
562  Qj.assign(Qtmp);
563 
564  ordinal_type t = piv[j];
565  piv[j] = piv[nfixed];
566  piv[nfixed] = t;
567  }
568  ++nfixed;
569  }
570  }
571 
572  scalar_type sigma = 1.0 + rank_threshold;
573  scalar_type r_max = 0.0;
574  ordinal_type j=0;
575  R.putScalar(0.0);
576  while (j < k && sigma >= rank_threshold) {
577 
578  SDV Qj = getCol(Teuchos::View, Q, j);
579 
580  // Find pivot column if we are passed the pre-pivot stage
581  if (j >= nfixed) {
582  ordinal_type pvt = j;
583  for (ordinal_type l=j+1; l<n; l++)
584  if (nrm(l) > nrm(pvt))
585  pvt = l;
586 
587  // Interchange column j and pvt
588  if (pvt != j) {
589  SDV Qpvt = getCol(Teuchos::View, Q, pvt);
590  Qtmp.assign(Qpvt);
591  Qpvt.assign(Qj);
592  Qj.assign(Qtmp);
593 
594  SDV Rpvt = getCol(Teuchos::View, R, pvt);
595  SDV Rj = getCol(Teuchos::View, R, j);
596  Rtmp.assign(Rpvt);
597  Rpvt.assign(Rj);
598  Rj.assign(Rtmp);
599 
600  ordinal_type t = piv[pvt];
601  piv[pvt] = piv[j];
602  piv[j] = t;
603  }
604  }
605 
606  R(j,j) = std::sqrt(detail::weighted_inner_product(Qj, Qj, w));
607  Qj.scale(1.0/R(j,j));
608  for (ordinal_type l=j+1; l<n; l++) {
609  SDV Ql = getCol(Teuchos::View, Q, l);
610  scalar_type t = detail::weighted_inner_product(Qj, Ql, w);
611  R(j,l) = t;
612  detail::saxpy(1.0, Ql, -t, Qj);
613 
614  // Update norms
615  nrm(l) = detail::weighted_inner_product(Ql, Ql, w);
616  }
617 
618  if (std::abs(R(j,j)) > r_max)
619  r_max = R(j,j);
620  sigma = std::abs(R(j,j)/r_max);
621  j++;
622  }
623 
624  ordinal_type rank = j;
625  if (sigma < rank_threshold)
626  rank--;
627  Q.reshape(m, rank);
628  R.reshape(rank, rank);
629 
630  return rank;
631  }
632 
634 
645  template <typename ordinal_type, typename scalar_type>
646  ordinal_type
647  CPQR_MGS_reorthog_threshold(
648  const scalar_type& rank_threshold,
649  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
650  const Teuchos::Array<scalar_type>& w,
651  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
652  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
653  Teuchos::Array<ordinal_type>& piv)
654  {
655  // First do standard column-pivoted QR
656  ordinal_type rank = CPQR_MGS_threshold(rank_threshold, A, w, Q, R, piv);
657 
658  // Now apply standard MGS to Q
659  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> A2(Q), R2;
660  QR_MGS(rank, A2, w, Q, R2);
661 
662  return rank;
663  }
664 
666  template <typename ordinal_type, typename scalar_type>
667  scalar_type
668  cond_R(const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R)
669  {
670  ordinal_type k = R.numRows();
671  if (R.numCols() < k)
672  k = R.numCols();
673  scalar_type r_max = std::abs(R(0,0));
674  scalar_type r_min = std::abs(R(0,0));
675  for (ordinal_type i=1; i<k; i++) {
676  if (std::abs(R(i,i)) > r_max)
677  r_max = R(i,i);
678  if (std::abs(R(i,i)) < r_min)
679  r_min = R(i,i);
680  }
681  scalar_type cond_r = r_max / r_min;
682  return cond_r;
683  }
684 
686 
690  template <typename ordinal_type, typename scalar_type>
691  scalar_type
692  weightedQROrthogonalizationError(
693  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
694  const Teuchos::Array<scalar_type>& w)
695  {
696  ordinal_type m = Q.numRows();
697  ordinal_type n = Q.numCols();
698  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> Qt(m,n);
699  for (ordinal_type i=0; i<m; i++)
700  for (ordinal_type j=0; j<n; j++)
701  Qt(i,j) = w[i]*Q(i,j);
702  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(n,n);
703  err.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q, Qt, 0.0);
704  for (ordinal_type i=0; i<n; i++)
705  err(i,i) -= 1.0;
706  return err.normInf();
707  }
708 
710 
713  template <typename ordinal_type, typename scalar_type>
714  scalar_type
715  QROrthogonalizationError(
716  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q)
717  {
718  ordinal_type n = Q.numCols();
719  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(n,n);
720  err.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q, Q, 0.0);
721  for (ordinal_type i=0; i<n; i++)
722  err(i,i) -= 1.0;
723  return err.normInf();
724  }
725 
727 
732  template <typename ordinal_type, typename scalar_type>
733  scalar_type
734  residualQRError(
735  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
736  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
737  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R)
738  {
739  ordinal_type k = Q.numCols();
740  ordinal_type m = Q.numRows();
741  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> AA(
742  Teuchos::View, A, m, k, 0, 0);
743  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(m,k);
744  ordinal_type ret =
745  err.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
746  TEUCHOS_ASSERT(ret == 0);
747  err -= AA;
748  return err.normInf();
749  }
750 
752 
758  template <typename ordinal_type, typename scalar_type>
759  scalar_type
760  residualCPQRError(
761  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
762  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
763  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R,
764  const Teuchos::Array<ordinal_type>& piv)
765  {
766  ordinal_type k = Q.numCols();
767  ordinal_type m = Q.numRows();
768  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> AP(m, k);
769  for (ordinal_type j=0; j<k; j++)
770  for (ordinal_type i=0; i<m; i++)
771  AP(i,j) = A(i,piv[j]);
772  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(m,k);
773  ordinal_type ret =
774  err.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
775  TEUCHOS_ASSERT(ret == 0);
776  err -= AP;
777 
778  return err.normInf();
779  }
780 
782 
785  template <typename ordinal_type, typename scalar_type>
786  void
787  svd(const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
788  Teuchos::Array<scalar_type>& s,
789  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& U,
790  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Vt)
791  {
792  Teuchos::LAPACK<ordinal_type,scalar_type> lapack;
793  ordinal_type m = A.numRows();
794  ordinal_type n = A.numCols();
795  ordinal_type k = std::min(m,n);
796  ordinal_type lda = A.stride();
797 
798  // Copy A since GESVD overwrites it (always)
799  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
800  Teuchos::Copy, A, m, n);
801 
802  // Resize appropriately
803  if (U.numRows() != m || U.numCols() != m)
804  U.shape(m,m);
805  if (Vt.numRows() != n || Vt.numCols() != n)
806  Vt.shape(n,n);
807  if (s.size() != k)
808  s.resize(k);
809  ordinal_type ldu = U.stride();
810  ordinal_type ldv = Vt.stride();
811 
812  // Workspace query
813  Teuchos::Array<scalar_type> work(1);
814  ordinal_type lwork = -1;
815  ordinal_type info;
816  lapack.GESVD('A', 'A', m, n, AA.values(), lda, &s[0], U.values(), ldu,
817  Vt.values(), ldv, &work[0], lwork, NULL, &info);
818  TEUCHOS_TEST_FOR_EXCEPTION(
819  info < 0, std::logic_error, "dgesvd returned info = " << info);
820 
821  // Do SVD
822  lwork = work[0];
823  work.resize(lwork);
824  lapack.GESVD('A', 'A', m, n, AA.values(), lda, &s[0], U.values(), ldu,
825  Vt.values(), ldv, &work[0], lwork, NULL, &info);
826  TEUCHOS_TEST_FOR_EXCEPTION(
827  info < 0, std::logic_error, "dgesvd returned info = " << info);
828 
829  }
830 
831  template <typename ordinal_type, typename scalar_type>
832  ordinal_type
833  svd_threshold(
834  const scalar_type& rank_threshold,
835  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
836  Teuchos::Array<scalar_type>& s,
837  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& U,
838  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Vt)
839  {
840  // Compute full SVD
841  svd(A, s, U, Vt);
842 
843  // Find leading block of S such that cond(S) <= 1/tau
844  ordinal_type rank = 0;
845  ordinal_type m = s.size();
846  while (rank < m && s[rank]/s[0] > rank_threshold) rank++;
847 
848  // Extract blocks from s, U and Vt
849  s.resize(rank);
850  U.reshape(U.numRows(),rank);
851  Vt.reshape(rank, Vt.numCols());
852 
853  return rank;
854  }
855 
856 }
857 
858 #endif
Stokhos::CPQR_MGS_threshold
ordinal_type CPQR_MGS_threshold(const scalar_type &rank_threshold, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::Array< scalar_type > &w, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R, Teuchos::Array< ordinal_type > &piv)
Compute column-pivoted QR using modified Gram-Schmidt.
Definition: Stokhos_SDMUtils.hpp:514
Teuchos::SerialDenseMatrix::values
ScalarType * values() const
Sacado::UQ::sqrt
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:847
Teuchos::LAPACK::ORGQR
void ORGQR(const OrdinalType &m, const OrdinalType &n, const OrdinalType &k, ScalarType *A, const OrdinalType &lda, const ScalarType *TAU, ScalarType *WORK, const OrdinalType &lwork, OrdinalType *info) const
Stokhos::residualQRError
scalar_type residualQRError(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R)
Compute QR residual error.
Definition: Stokhos_SDMUtils.hpp:734
Teuchos_SerialDenseHelpers.hpp
Teuchos_SerialDenseMatrix.hpp
Stokhos::CPQR_Householder_threshold
ordinal_type CPQR_Householder_threshold(const scalar_type &rank_threshold, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::Array< scalar_type > &w, 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:462
Stokhos::print_matlab
void print_matlab(std::ostream &os, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A)
Definition: Stokhos_SDMUtils.hpp:83
Stokhos::QR_CGS
void QR_CGS(ordinal_type k, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::Array< scalar_type > &w, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R)
Compute thin QR using classical Gram-Schmidt.
Definition: Stokhos_SDMUtils.hpp:151
Stokhos::QR_MGS
void QR_MGS(ordinal_type k, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::Array< scalar_type > &w, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R)
Compute thin QR using modified Gram-Schmidt.
Definition: Stokhos_SDMUtils.hpp:193
Stokhos::cond_R
scalar_type cond_R(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R)
Compute condition number of upper-triangular R.
Definition: Stokhos_SDMUtils.hpp:668
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)
DGEQP3_F77
#define DGEQP3_F77
Definition: Stokhos_SDMUtils.hpp:20
Teuchos::SerialDenseMatrix::scale
int scale(const ScalarType alpha)
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::detail::saxpy
void saxpy(const scalar_type &alpha, Teuchos::SerialDenseVector< ordinal_type, scalar_type > &x, const scalar_type &beta, const Teuchos::SerialDenseVector< ordinal_type, scalar_type > &y)
Overwrite x with alpha*x + beta*y.
Definition: Stokhos_SDMUtils.hpp:67
SDM
Teuchos::SerialDenseMatrix< int, pce_type > SDM
Definition: stieltjes_example5.cpp:46
Stokhos::svd_threshold
ordinal_type svd_threshold(const scalar_type &rank_threshold, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, Teuchos::Array< scalar_type > &s, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &U, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Vt)
Definition: Stokhos_SDMUtils.hpp:833
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
Teuchos_LAPACK.hpp
SDV
Teuchos::SerialDenseVector< int, pce_type > SDV
Definition: stieltjes_example5.cpp:47
Stokhos::residualCPQRError
scalar_type residualCPQRError(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R, const Teuchos::Array< ordinal_type > &piv)
Compute column-pivoted QR residual error.
Definition: Stokhos_SDMUtils.hpp:760
Teuchos::SerialDenseMatrix::putScalar
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Stokhos::weightedQROrthogonalizationError
scalar_type weightedQROrthogonalizationError(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, const Teuchos::Array< scalar_type > &w)
Compute weighted QR orthogonalization error.
Definition: Stokhos_SDMUtils.hpp:692
Teuchos::LAPACK
Specialization for Sacado::UQ::PCE&lt; Storage&lt;...&gt; &gt;
Teuchos::SerialDenseVector::length
OrdinalType length() const
Stokhos::QROrthogonalizationError
scalar_type QROrthogonalizationError(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q)
Compute QR orthogonalization error.
Definition: Stokhos_SDMUtils.hpp:715
Sacado::UQ::abs
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1150
Teuchos::Array::resize
void resize(size_type new_size, const value_type &x=value_type())
Stokhos::svd
void svd(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, Teuchos::Array< scalar_type > &s, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &U, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Vt)
Compute SVD of matrix.
Definition: Stokhos_SDMUtils.hpp:787
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
Teuchos::SerialDenseMatrix::numCols
OrdinalType numCols() const
Teuchos_Array.hpp
Teuchos::SerialDenseMatrix::normInf
ScalarTraits< ScalarType >::magnitudeType normInf() const
Teuchos::LAPACK::GEQRF
void GEQRF(const OrdinalType &m, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, ScalarType *TAU, ScalarType *WORK, const OrdinalType &lwork, OrdinalType *info) const
Teuchos::LAPACK::GESVD
void GESVD(const char &JOBU, const char &JOBVT, const OrdinalType &m, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, MagnitudeType *S, ScalarType *U, const OrdinalType &ldu, ScalarType *V, const OrdinalType &ldv, ScalarType *WORK, const OrdinalType &lwork, MagnitudeType *RWORK, OrdinalType *info) const
Stokhos::QR_MGS2
void QR_MGS2(ordinal_type k, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::Array< scalar_type > &w, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R)
Compute thin QR using modified Gram-Schmidt with reorthogonalization.
Definition: Stokhos_SDMUtils.hpp:238
Teuchos::SerialDenseMatrix::reshape
int reshape(OrdinalType numRows, OrdinalType numCols)
Teuchos::Array::size
size_type size() const
Teuchos_SerialDenseVector.hpp
MKL_NO_EXCEPT
#define MKL_NO_EXCEPT
Definition: Stokhos_SDMUtils.hpp:31
Stokhos::detail::weighted_inner_product
scalar_type weighted_inner_product(const Teuchos::SerialDenseVector< ordinal_type, scalar_type > &x, const Teuchos::SerialDenseVector< ordinal_type, scalar_type > &y, const Teuchos::Array< scalar_type > &w)
Compute weighted inner product between vectors x and y.
Definition: Stokhos_SDMUtils.hpp:53
Teuchos::SerialDenseMatrix::shape
int shape(OrdinalType numRows, OrdinalType numCols)
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
Stokhos::QR_Householder
void QR_Householder(ordinal_type k, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::Array< scalar_type > &w, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R)
Compute thin QR using Householder reflections.
Definition: Stokhos_SDMUtils.hpp:292
Teuchos::SerialDenseMatrix::assign
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
n
int n
Teuchos::SerialDenseMatrix::stride
OrdinalType stride() const
Stokhos::CPQR_MGS_reorthog_threshold
ordinal_type CPQR_MGS_reorthog_threshold(const scalar_type &rank_threshold, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::Array< scalar_type > &w, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R, Teuchos::Array< ordinal_type > &piv)
Compute column-pivoted QR using modified Gram-Schmidt and reorthogonalization.
Definition: Stokhos_SDMUtils.hpp:647
Teuchos::SerialDenseMatrix::numRows
OrdinalType numRows() const
Generated on Fri Nov 22 2024 09:18:14 for Stokhos Package Browser (Single Doxygen Collection) by   doxygen 1.8.5