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 //
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 #ifndef STOKHOS_SDM_UTILS_HPP
43 #define STOKHOS_SDM_UTILS_HPP
44 
45 #include "Teuchos_SerialDenseMatrix.hpp"
46 #include "Teuchos_SerialDenseVector.hpp"
47 #include "Teuchos_SerialDenseHelpers.hpp"
48 #include "Teuchos_Array.hpp"
49 #include "Teuchos_LAPACK.hpp"
50 #include <ostream>
51 
52 #define DGEQP3_F77 F77_BLAS_MANGLE(dgeqp3,DGEQP3)
53 extern "C" {
54 void DGEQP3_F77(const int*, const int*, double*, const int*, int*,
55  double*, double*, const int*, int*);
56 }
57 
58 #include "Stokhos_ConfigDefs.h"
59 
60 #ifdef HAVE_STOKHOS_MATLABLIB
61 extern "C" {
62 #include "mat.h"
63 #include "matrix.h"
64 }
65 #endif
66 
67 namespace Stokhos {
68 
69  namespace detail {
70 
72  template <typename ordinal_type, typename scalar_type>
73  scalar_type weighted_inner_product(
74  const Teuchos::SerialDenseVector<ordinal_type,scalar_type>& x,
75  const Teuchos::SerialDenseVector<ordinal_type,scalar_type>& y,
76  const Teuchos::Array<scalar_type>& w)
77  {
78  ordinal_type n = x.length();
79  scalar_type t = 0;
80  for (ordinal_type i=0; i<n; i++)
81  t += x[i]*w[i]*y[i];
82  return t;
83  }
84 
86  template <typename ordinal_type, typename scalar_type>
87  void saxpy(
88  const scalar_type& alpha,
89  Teuchos::SerialDenseVector<ordinal_type, scalar_type>& x,
90  const scalar_type& beta,
91  const Teuchos::SerialDenseVector<ordinal_type, scalar_type>& y)
92  {
93  ordinal_type n = x.length();
94  for (ordinal_type i=0; i<n; i++)
95  x[i] = alpha*x[i] + beta*y[i];
96  }
97 
98  }
99 
100  // Print a matrix in matlab-esque format
101  template <typename ordinal_type, typename scalar_type>
102  void
103  print_matlab(std::ostream& os,
104  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A)
105  {
106  os << "[ ";
107  for (ordinal_type i=0; i<A.numRows(); i++) {
108  for (ordinal_type j=0; j<A.numCols(); j++)
109  os << A(i,j) << " ";
110  if (i < A.numRows()-1)
111  os << ";" << std::endl << " ";
112  }
113  os << "];" << std::endl;
114  }
115 
116 #ifdef HAVE_STOKHOS_MATLABLIB
117  // Save a matrix to matlab MAT file
118  template <typename ordinal_type>
119  void
120  save_matlab(const std::string& file_name, const std::string& matrix_name,
121  const Teuchos::SerialDenseMatrix<ordinal_type,double>& A,
122  bool overwrite = true)
123  {
124  // Open matfile
125  const char *mode;
126  if (overwrite)
127  mode = "w";
128  else
129  mode = "u";
130  MATFile *file = matOpen(file_name.c_str(), mode);
131  TEUCHOS_ASSERT(file != NULL);
132 
133  // Convert SDM to Matlab matrix
134  mxArray *mat = mxCreateDoubleMatrix(0, 0, mxREAL);
135  TEUCHOS_ASSERT(mat != NULL);
136  mxSetPr(mat, A.values());
137  mxSetM(mat, A.numRows());
138  mxSetN(mat, A.numCols());
139 
140  // Save matrix to file
141  ordinal_type ret = matPutVariable(file, matrix_name.c_str(), mat);
142  TEUCHOS_ASSERT(ret == 0);
143 
144  // Close file
145  ret = matClose(file);
146  TEUCHOS_ASSERT(ret == 0);
147 
148  // Free matlab matrix
149  mxSetPr(mat, NULL);
150  mxSetM(mat, 0);
151  mxSetN(mat, 0);
152  mxDestroyArray(mat);
153  }
154 #endif
155 
157  template <typename ordinal_type, typename scalar_type>
158  scalar_type vec_norm_inf(
159  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A) {
160  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> vec_A(
161  Teuchos::View, A.values(), 1, A.numRows()*A.numCols(), 1);
162  return vec_A.normInf();
163  }
164 
166 
170  template <typename ordinal_type, typename scalar_type>
171  void QR_CGS(
172  ordinal_type k,
173  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
174  const Teuchos::Array<scalar_type>& w,
175  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
176  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R)
177  {
178  using Teuchos::getCol;
179  typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type> SDV;
180  typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> SDM;
181 
182  // getCol() requires non-const SDM
183  SDM& Anc = const_cast<SDM&>(A);
184 
185  // Make sure Q is the right size
186  ordinal_type m = A.numRows();
187  if (Q.numRows() != m || Q.numCols() != k)
188  Q.shape(m,k);
189  if (R.numRows() != k || R.numCols() != k)
190  R.shape(k,k);
191 
192  for (ordinal_type j=0; j<k; j++) {
193  SDV Aj = getCol(Teuchos::View, Anc, j);
194  SDV Qj = getCol(Teuchos::View, Q, j);
195  Qj.assign(Aj);
196  for (ordinal_type i=0; i<j; i++) {
197  SDV Qi = getCol(Teuchos::View, Q, i);
198  R(i,j) = detail::weighted_inner_product(Qi, Aj, w);
199  detail::saxpy(1.0, Qj, -R(i,j), Qi); // Q(:,j) = 1.0*Q(:,j) - R(i,j)*Q(:,i)
200  }
201  R(j,j) = std::sqrt(detail::weighted_inner_product(Qj, Qj, w));
202  Qj.scale(1.0/R(j,j));
203  }
204 
205  }
206 
208 
212  template <typename ordinal_type, typename scalar_type>
213  void QR_MGS(
214  ordinal_type k,
215  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
216  const Teuchos::Array<scalar_type>& w,
217  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
218  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R)
219  {
220  using Teuchos::getCol;
221  typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type> SDV;
222  typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> SDM;
223 
224  // getCol() requires non-const SDM
225  SDM& Anc = const_cast<SDM&>(A);
226 
227  // Make sure Q is the right size
228  ordinal_type m = A.numRows();
229  if (Q.numRows() != m || Q.numCols() != k)
230  Q.shape(m,k);
231  if (R.numRows() != k || R.numCols() != k)
232  R.shape(k,k);
233 
234  for (ordinal_type i=0; i<k; i++) {
235  SDV Ai = getCol(Teuchos::View, Anc, i);
236  SDV Qi = getCol(Teuchos::View, Q, i);
237  Qi.assign(Ai);
238  }
239  for (ordinal_type i=0; i<k; i++) {
240  SDV Qi = getCol(Teuchos::View, Q, i);
241  for (ordinal_type j=0; j<i; j++) {
242  SDV Qj = getCol(Teuchos::View, Q, j);
243  scalar_type v = detail::weighted_inner_product(Qi, Qj, w);
244  detail::saxpy(1.0, Qi, -v, Qj); // Q(:,i) = 1.0*Q(:,i) - v*Q(:,j)
245  R(j,i) += v;
246  }
247  R(i,i) = std::sqrt(detail::weighted_inner_product(Qi, Qi, w));
248  Qi.scale(1.0/R(i,i));
249  }
250  }
251 
253 
257  template <typename ordinal_type, typename scalar_type>
258  void QR_MGS2(
259  ordinal_type k,
260  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
261  const Teuchos::Array<scalar_type>& w,
262  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
263  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R)
264  {
265  using Teuchos::getCol;
266  typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type> SDV;
267  typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> SDM;
268 
269  // getCol() requires non-const SDM
270  SDM& Anc = const_cast<SDM&>(A);
271 
272  // Make sure Q is the right size
273  ordinal_type m = A.numRows();
274  if (Q.numRows() != m || Q.numCols() != k)
275  Q.shape(m,k);
276  if (R.numRows() != k || R.numCols() != k)
277  R.shape(k,k);
278 
279  for (ordinal_type i=0; i<k; i++) {
280  SDV Ai = getCol(Teuchos::View, Anc, i);
281  SDV Qi = getCol(Teuchos::View, Q, i);
282  Qi.assign(Ai);
283  }
284  for (ordinal_type i=0; i<k; i++) {
285  SDV Qi = getCol(Teuchos::View, Q, i);
286  for (ordinal_type j=0; j<i; j++) {
287  SDV Qj = getCol(Teuchos::View, Q, j);
288  scalar_type v = detail::weighted_inner_product(Qi, Qj, w);
289  detail::saxpy(1.0, Qi, -v, Qj); // Q(:,i) = 1.0*Q(:,i) - v*Q(:,j)
290  R(j,i) += v;
291  }
292  for (ordinal_type j=i-1; j>=0; j--) {
293  SDV Qj = getCol(Teuchos::View, Q, j);
294  scalar_type v = detail::weighted_inner_product(Qi, Qj, w);
295  detail::saxpy(1.0, Qi, -v, Qj); // Q(:,i) = 1.0*Q(:,i) - v*Q(:,j)
296  R(j,i) += v;
297  }
298  R(i,i) = std::sqrt(detail::weighted_inner_product(Qi, Qi, w));
299  Qi.scale(1.0/R(i,i));
300  }
301  }
302 
304 
310  template <typename ordinal_type, typename scalar_type>
311  void
312  QR_Householder(
313  ordinal_type k,
314  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
315  const Teuchos::Array<scalar_type>& w,
316  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
317  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R)
318  {
319  Teuchos::LAPACK<ordinal_type,scalar_type> lapack;
320  ordinal_type m = A.numRows();
321  ordinal_type n = A.numCols();
322  ordinal_type kk = std::min(m,n);
323  if (k > kk)
324  k = kk;
325 
326  // Check that each component of w is 1
327  for (ordinal_type i=0; i<w.size(); i++)
328  TEUCHOS_TEST_FOR_EXCEPTION(
329  w[i] != 1.0, std::logic_error,
330  "CPQR_Householder_threshold() requires unit weight vector!");
331 
332  // Lapack routine overwrites A, so copy into temporary matrix
333  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
334  Teuchos::Copy, A, m, n);
335 
336  // QR
337  ordinal_type lda = AA.stride();
338  Teuchos::Array<scalar_type> tau(k);
339  Teuchos::Array<scalar_type> work(1);
340  ordinal_type info;
341  ordinal_type lwork = -1;
342  lapack.GEQRF(m, n, AA.values(), lda, &tau[0], &work[0], lwork, &info);
343  TEUCHOS_TEST_FOR_EXCEPTION(
344  info < 0, std::logic_error, "geqrf returned info = " << info);
345  lwork = work[0];
346  work.resize(lwork);
347  lapack.GEQRF(m, n, AA.values(), lda, &tau[0], &work[0], lwork, &info);
348  TEUCHOS_TEST_FOR_EXCEPTION(
349  info < 0, std::logic_error, "geqrf returned info = " << info);
350 
351  // Extract R
352  if (R.numRows() != k || R.numCols() != n)
353  R.shape(k,n);
354  R.putScalar(0.0);
355  for (ordinal_type i=0; i<k; i++)
356  for (ordinal_type j=i; j<n; j++)
357  R(i,j) = AA(i,j);
358 
359  // Extract Q
360  if (Q.numRows() != m || Q.numCols() != k)
361  Q.shape(m,k);
362  lwork = -1;
363  lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
364  TEUCHOS_TEST_FOR_EXCEPTION(
365  info < 0, std::logic_error, "orgqr returned info = " << info);
366  lwork = work[0];
367  work.resize(lwork);
368  lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
369  TEUCHOS_TEST_FOR_EXCEPTION(
370  info < 0, std::logic_error, "orgqr returned info = " << info);
371  if (Q.numRows() != m || Q.numCols() != k)
372  Q.shape(m, k);
373  for (ordinal_type i=0; i<m; i++)
374  for (ordinal_type j=0; j<k; j++)
375  Q(i,j) = AA(i,j);
376  }
377 
379 
391  template <typename ordinal_type, typename scalar_type>
392  void
393  CPQR_Householder3(
394  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
395  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
396  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
397  Teuchos::Array<ordinal_type>& piv)
398  {
399  Teuchos::LAPACK<ordinal_type,scalar_type> lapack;
400  ordinal_type m = A.numRows();
401  ordinal_type n = A.numCols();
402  ordinal_type k = std::min(m,n);
403 
404  // Lapack routine overwrites A, so copy into temporary matrix
405  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
406  Teuchos::Copy, A, m, n);
407  if (Q.numRows() != m || Q.numCols() != k)
408  Q.shape(m,k);
409 
410  // Teuchos LAPACK interface doesn't have dgeqp3, so call it directly
411 
412  // Workspace query
413  ordinal_type lda = AA.stride();
414  piv.resize(n);
415  Teuchos::Array<scalar_type> tau(k);
416  Teuchos::Array<scalar_type> work(1);
417  ordinal_type lwork = -1;
418  ordinal_type info;
419  DGEQP3_F77(&m, &n, AA.values(), &lda, &piv[0], &tau[0], &work[0], &lwork,
420  &info);
421  TEUCHOS_TEST_FOR_EXCEPTION(
422  info < 0, std::logic_error, "dgeqp3 returned info = " << info);
423 
424  // Column pivoted QR
425  lwork = work[0];
426  work.resize(lwork);
427  DGEQP3_F77(&m, &n, AA.values(), &lda, &piv[0], &tau[0], &work[0], &lwork,
428  &info);
429  TEUCHOS_TEST_FOR_EXCEPTION(
430  info < 0, std::logic_error, "dgeqp3 returned info = " << info);
431 
432  // Extract R
433  if (R.numRows() != k || R.numCols() != n)
434  R.shape(k,n);
435  R.putScalar(0.0);
436  for (ordinal_type i=0; i<k; i++)
437  for (ordinal_type j=i; j<n; j++)
438  R(i,j) = AA(i,j);
439 
440  // Extract Q
441  lwork = -1;
442  lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
443  TEUCHOS_TEST_FOR_EXCEPTION(
444  info < 0, std::logic_error, "orgqr returned info = " << info);
445  lwork = work[0];
446  work.resize(lwork);
447  lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
448  TEUCHOS_TEST_FOR_EXCEPTION(
449  info < 0, std::logic_error, "orgqr returned info = " << info);
450  if (Q.numRows() != m || Q.numCols() != k)
451  Q.shape(m, k);
452  for (ordinal_type i=0; i<m; i++)
453  for (ordinal_type j=0; j<k; j++)
454  Q(i,j) = AA(i,j);
455 
456  // Transform piv to zero-based indexing
457  for (ordinal_type i=0; i<n; i++)
458  piv[i] -= 1;
459  }
460 
462 
480  template <typename ordinal_type, typename scalar_type>
481  ordinal_type
482  CPQR_Householder_threshold(
483  const scalar_type& rank_threshold,
484  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
485  const Teuchos::Array<scalar_type>& w,
486  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
487  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
488  Teuchos::Array<ordinal_type>& piv)
489  {
490  // Check that each component of w is 1
491  for (ordinal_type i=0; i<w.size(); i++)
492  TEUCHOS_TEST_FOR_EXCEPTION(
493  w[i] != 1.0, std::logic_error,
494  "CPQR_Householder_threshold() requires unit weight vector!");
495 
496  // Compute full QR
497  CPQR_Householder3(A, Q, R, piv);
498 
499  // Find leading block of R such that cond(R) <= 1/tau
500  ordinal_type rank = 0;
501  ordinal_type m = R.numRows();
502  scalar_type r_max = std::abs(R(rank,rank));
503  scalar_type r_min = std::abs(R(rank,rank));
504  for (rank=1; rank<m; rank++) {
505  if (std::abs(R(rank,rank)) > r_max)
506  r_max = std::abs(R(rank,rank));
507  if (std::abs(R(rank,rank)) < r_min)
508  r_min = std::abs(R(rank,rank));
509  if (r_min / r_max < rank_threshold)
510  break;
511  }
512 
513  // Extract blocks from Q and R
514  R.reshape(rank,rank);
515  Q.reshape(Q.numRows(), rank);
516 
517  return rank;
518  }
519 
521 
532  template <typename ordinal_type, typename scalar_type>
533  ordinal_type
534  CPQR_MGS_threshold(
535  const scalar_type& rank_threshold,
536  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
537  const Teuchos::Array<scalar_type>& w,
538  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
539  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
540  Teuchos::Array<ordinal_type>& piv)
541  {
542  using Teuchos::getCol;
543  typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type> SDV;
544  ordinal_type m = A.numRows();
545  ordinal_type n = A.numCols();
546  ordinal_type k = std::min(m,n);
547 
548  // Make sure Q is the right size
549  if (Q.numRows() != m || Q.numCols() != n)
550  Q.shape(m,n);
551  if (R.numRows() != m || R.numCols() != n)
552  R.shape(m,n);
553  if (piv.size() != n)
554  piv.resize(n);
555  Q.assign(A);
556 
557  // Compute column norms
558  SDV nrm(n);
559  for (ordinal_type j=0; j<n; j++) {
560  SDV Qj = getCol(Teuchos::View, Q, j);
561  nrm(j) = detail::weighted_inner_product(Qj, Qj, w);
562  }
563  SDV Qtmp(m), Rtmp(m);
564 
565  Teuchos::Array<ordinal_type> piv_orig(piv);
566  for (ordinal_type i=0; i<n; i++)
567  piv[i] = i;
568 
569  // Prepivot any columns requested by setting piv[i] != 0
570  ordinal_type nfixed = 0;
571  for (ordinal_type j=0; j<n; j++) {
572  if (piv_orig[j] != 0) {
573  if (j != nfixed) {
574  scalar_type tmp = nrm(j);
575  nrm(j) = nrm(nfixed);
576  nrm(nfixed) = tmp;
577 
578  SDV Qpvt = getCol(Teuchos::View, Q, j);
579  SDV Qj = getCol(Teuchos::View, Q, nfixed);
580  Qtmp.assign(Qpvt);
581  Qpvt.assign(Qj);
582  Qj.assign(Qtmp);
583 
584  ordinal_type t = piv[j];
585  piv[j] = piv[nfixed];
586  piv[nfixed] = t;
587  }
588  ++nfixed;
589  }
590  }
591 
592  scalar_type sigma = 1.0 + rank_threshold;
593  scalar_type r_max = 0.0;
594  ordinal_type j=0;
595  R.putScalar(0.0);
596  while (j < k && sigma >= rank_threshold) {
597 
598  SDV Qj = getCol(Teuchos::View, Q, j);
599 
600  // Find pivot column if we are passed the pre-pivot stage
601  if (j >= nfixed) {
602  ordinal_type pvt = j;
603  for (ordinal_type l=j+1; l<n; l++)
604  if (nrm(l) > nrm(pvt))
605  pvt = l;
606 
607  // Interchange column j and pvt
608  if (pvt != j) {
609  SDV Qpvt = getCol(Teuchos::View, Q, pvt);
610  Qtmp.assign(Qpvt);
611  Qpvt.assign(Qj);
612  Qj.assign(Qtmp);
613 
614  SDV Rpvt = getCol(Teuchos::View, R, pvt);
615  SDV Rj = getCol(Teuchos::View, R, j);
616  Rtmp.assign(Rpvt);
617  Rpvt.assign(Rj);
618  Rj.assign(Rtmp);
619 
620  ordinal_type t = piv[pvt];
621  piv[pvt] = piv[j];
622  piv[j] = t;
623  }
624  }
625 
626  R(j,j) = std::sqrt(detail::weighted_inner_product(Qj, Qj, w));
627  Qj.scale(1.0/R(j,j));
628  for (ordinal_type l=j+1; l<n; l++) {
629  SDV Ql = getCol(Teuchos::View, Q, l);
630  scalar_type t = detail::weighted_inner_product(Qj, Ql, w);
631  R(j,l) = t;
632  detail::saxpy(1.0, Ql, -t, Qj);
633 
634  // Update norms
635  nrm(l) = detail::weighted_inner_product(Ql, Ql, w);
636  }
637 
638  if (std::abs(R(j,j)) > r_max)
639  r_max = R(j,j);
640  sigma = std::abs(R(j,j)/r_max);
641  j++;
642  }
643 
644  ordinal_type rank = j;
645  if (sigma < rank_threshold)
646  rank--;
647  Q.reshape(m, rank);
648  R.reshape(rank, rank);
649 
650  return rank;
651  }
652 
654 
665  template <typename ordinal_type, typename scalar_type>
666  ordinal_type
667  CPQR_MGS_reorthog_threshold(
668  const scalar_type& rank_threshold,
669  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
670  const Teuchos::Array<scalar_type>& w,
671  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
672  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
673  Teuchos::Array<ordinal_type>& piv)
674  {
675  // First do standard column-pivoted QR
676  ordinal_type rank = CPQR_MGS_threshold(rank_threshold, A, w, Q, R, piv);
677 
678  // Now apply standard MGS to Q
679  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> A2(Q), R2;
680  QR_MGS(rank, A2, w, Q, R2);
681 
682  return rank;
683  }
684 
686  template <typename ordinal_type, typename scalar_type>
687  scalar_type
688  cond_R(const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R)
689  {
690  ordinal_type k = R.numRows();
691  if (R.numCols() < k)
692  k = R.numCols();
693  scalar_type r_max = std::abs(R(0,0));
694  scalar_type r_min = std::abs(R(0,0));
695  for (ordinal_type i=1; i<k; i++) {
696  if (std::abs(R(i,i)) > r_max)
697  r_max = R(i,i);
698  if (std::abs(R(i,i)) < r_min)
699  r_min = R(i,i);
700  }
701  scalar_type cond_r = r_max / r_min;
702  return cond_r;
703  }
704 
706 
710  template <typename ordinal_type, typename scalar_type>
711  scalar_type
712  weightedQROrthogonalizationError(
713  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
714  const Teuchos::Array<scalar_type>& w)
715  {
716  ordinal_type m = Q.numRows();
717  ordinal_type n = Q.numCols();
718  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> Qt(m,n);
719  for (ordinal_type i=0; i<m; i++)
720  for (ordinal_type j=0; j<n; j++)
721  Qt(i,j) = w[i]*Q(i,j);
722  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(n,n);
723  err.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q, Qt, 0.0);
724  for (ordinal_type i=0; i<n; i++)
725  err(i,i) -= 1.0;
726  return err.normInf();
727  }
728 
730 
733  template <typename ordinal_type, typename scalar_type>
734  scalar_type
735  QROrthogonalizationError(
736  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q)
737  {
738  ordinal_type n = Q.numCols();
739  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(n,n);
740  err.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q, Q, 0.0);
741  for (ordinal_type i=0; i<n; i++)
742  err(i,i) -= 1.0;
743  return err.normInf();
744  }
745 
747 
752  template <typename ordinal_type, typename scalar_type>
753  scalar_type
754  residualQRError(
755  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
756  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
757  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R)
758  {
759  ordinal_type k = Q.numCols();
760  ordinal_type m = Q.numRows();
761  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> AA(
762  Teuchos::View, A, m, k, 0, 0);
763  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(m,k);
764  ordinal_type ret =
765  err.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
766  TEUCHOS_ASSERT(ret == 0);
767  err -= AA;
768  return err.normInf();
769  }
770 
772 
778  template <typename ordinal_type, typename scalar_type>
779  scalar_type
780  residualCPQRError(
781  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
782  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
783  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R,
784  const Teuchos::Array<ordinal_type>& piv)
785  {
786  ordinal_type k = Q.numCols();
787  ordinal_type m = Q.numRows();
788  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> AP(m, k);
789  for (ordinal_type j=0; j<k; j++)
790  for (ordinal_type i=0; i<m; i++)
791  AP(i,j) = A(i,piv[j]);
792  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(m,k);
793  ordinal_type ret =
794  err.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
795  TEUCHOS_ASSERT(ret == 0);
796  err -= AP;
797 
798  return err.normInf();
799  }
800 
802 
805  template <typename ordinal_type, typename scalar_type>
806  void
807  svd(const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
808  Teuchos::Array<scalar_type>& s,
809  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& U,
810  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Vt)
811  {
812  Teuchos::LAPACK<ordinal_type,scalar_type> lapack;
813  ordinal_type m = A.numRows();
814  ordinal_type n = A.numCols();
815  ordinal_type k = std::min(m,n);
816  ordinal_type lda = A.stride();
817 
818  // Copy A since GESVD overwrites it (always)
819  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
820  Teuchos::Copy, A, m, n);
821 
822  // Resize appropriately
823  if (U.numRows() != m || U.numCols() != m)
824  U.shape(m,m);
825  if (Vt.numRows() != n || Vt.numCols() != n)
826  Vt.shape(n,n);
827  if (s.size() != k)
828  s.resize(k);
829  ordinal_type ldu = U.stride();
830  ordinal_type ldv = Vt.stride();
831 
832  // Workspace query
833  Teuchos::Array<scalar_type> work(1);
834  ordinal_type lwork = -1;
835  ordinal_type info;
836  lapack.GESVD('A', 'A', m, n, AA.values(), lda, &s[0], U.values(), ldu,
837  Vt.values(), ldv, &work[0], lwork, NULL, &info);
838  TEUCHOS_TEST_FOR_EXCEPTION(
839  info < 0, std::logic_error, "dgesvd returned info = " << info);
840 
841  // Do SVD
842  lwork = work[0];
843  work.resize(lwork);
844  lapack.GESVD('A', 'A', m, n, AA.values(), lda, &s[0], U.values(), ldu,
845  Vt.values(), ldv, &work[0], lwork, NULL, &info);
846  TEUCHOS_TEST_FOR_EXCEPTION(
847  info < 0, std::logic_error, "dgesvd returned info = " << info);
848 
849  }
850 
851  template <typename ordinal_type, typename scalar_type>
852  ordinal_type
853  svd_threshold(
854  const scalar_type& rank_threshold,
855  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
856  Teuchos::Array<scalar_type>& s,
857  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& U,
858  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Vt)
859  {
860  // Compute full SVD
861  svd(A, s, U, Vt);
862 
863  // Find leading block of S such that cond(S) <= 1/tau
864  ordinal_type rank = 0;
865  ordinal_type m = s.size();
866  while (rank < m && s[rank]/s[0] > rank_threshold) rank++;
867 
868  // Extract blocks from s, U and Vt
869  s.resize(rank);
870  U.reshape(U.numRows(),rank);
871  Vt.reshape(rank, Vt.numCols());
872 
873  return rank;
874  }
875 
876 }
877 
878 #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:534
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:878
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:754
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:482
Stokhos::print_matlab
void print_matlab(std::ostream &os, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A)
Definition: Stokhos_SDMUtils.hpp:103
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:171
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:213
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:688
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:52
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:393
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:87
SDM
Teuchos::SerialDenseMatrix< int, pce_type > SDM
Definition: stieltjes_example5.cpp:80
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:853
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
Teuchos_LAPACK.hpp
SDV
Teuchos::SerialDenseVector< int, pce_type > SDV
Definition: stieltjes_example5.cpp:81
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:780
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:712
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:735
Sacado::UQ::abs
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1181
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:807
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
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:258
Teuchos::SerialDenseMatrix::reshape
int reshape(OrdinalType numRows, OrdinalType numCols)
Teuchos::Array::size
size_type size() const
Teuchos_SerialDenseVector.hpp
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:73
Teuchos::SerialDenseMatrix::shape
int shape(OrdinalType numRows, OrdinalType numCols)
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
Stokhos_ConfigDefs.h
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:312
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:667
Teuchos::SerialDenseMatrix::numRows
OrdinalType numRows() const
Generated by   doxygen 1.8.5