Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
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 DGEQPF_F77 F77_BLAS_MANGLE(dgeqpf,DGEQPF)
53 #define DGEQP3_F77 F77_BLAS_MANGLE(dgeqp3,DGEQP3)
54 extern "C" {
55 void DGEQPF_F77(const int*, const int*, double*, const int*, int*, double*,
56  double*, int*);
57 void DGEQP3_F77(const int*, const int*, double*, const int*, int*,
58  double*, double*, const int*, int*);
59 }
60 
61 #include "Stokhos_ConfigDefs.h"
62 
63 #ifdef HAVE_STOKHOS_MATLABLIB
64 extern "C" {
65 #include "mat.h"
66 #include "matrix.h"
67 }
68 #endif
69 
70 namespace Stokhos {
71 
72  namespace detail {
73 
75  template <typename ordinal_type, typename scalar_type>
77  const Teuchos::SerialDenseVector<ordinal_type,scalar_type>& x,
78  const Teuchos::SerialDenseVector<ordinal_type,scalar_type>& y,
79  const Teuchos::Array<scalar_type>& w)
80  {
81  ordinal_type n = x.length();
82  scalar_type t = 0;
83  for (ordinal_type i=0; i<n; i++)
84  t += x[i]*w[i]*y[i];
85  return t;
86  }
87 
89  template <typename ordinal_type, typename scalar_type>
90  void saxpy(
91  const scalar_type& alpha,
92  Teuchos::SerialDenseVector<ordinal_type, scalar_type>& x,
93  const scalar_type& beta,
94  const Teuchos::SerialDenseVector<ordinal_type, scalar_type>& y)
95  {
96  ordinal_type n = x.length();
97  for (ordinal_type i=0; i<n; i++)
98  x[i] = alpha*x[i] + beta*y[i];
99  }
100 
101  }
102 
103  // Print a matrix in matlab-esque format
104  template <typename ordinal_type, typename scalar_type>
105  void
106  print_matlab(std::ostream& os,
107  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A)
108  {
109  os << "[ ";
110  for (ordinal_type i=0; i<A.numRows(); i++) {
111  for (ordinal_type j=0; j<A.numCols(); j++)
112  os << A(i,j) << " ";
113  if (i < A.numRows()-1)
114  os << ";" << std::endl << " ";
115  }
116  os << "];" << std::endl;
117  }
118 
119 #ifdef HAVE_STOKHOS_MATLABLIB
120  // Save a matrix to matlab MAT file
121  template <typename ordinal_type>
122  void
123  save_matlab(const std::string& file_name, const std::string& matrix_name,
124  const Teuchos::SerialDenseMatrix<ordinal_type,double>& A,
125  bool overwrite = true)
126  {
127  // Open matfile
128  const char *mode;
129  if (overwrite)
130  mode = "w";
131  else
132  mode = "u";
133  MATFile *file = matOpen(file_name.c_str(), mode);
134  TEUCHOS_ASSERT(file != NULL);
135 
136  // Convert SDM to Matlab matrix
137  mxArray *mat = mxCreateDoubleMatrix(0, 0, mxREAL);
138  TEUCHOS_ASSERT(mat != NULL);
139  mxSetPr(mat, A.values());
140  mxSetM(mat, A.numRows());
141  mxSetN(mat, A.numCols());
142 
143  // Save matrix to file
144  ordinal_type ret = matPutVariable(file, matrix_name.c_str(), mat);
145  TEUCHOS_ASSERT(ret == 0);
146 
147  // Close file
148  ret = matClose(file);
149  TEUCHOS_ASSERT(ret == 0);
150 
151  // Free matlab matrix
152  mxSetPr(mat, NULL);
153  mxSetM(mat, 0);
154  mxSetN(mat, 0);
155  mxDestroyArray(mat);
156  }
157 #endif
158 
160  template <typename ordinal_type, typename scalar_type>
162  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A) {
163  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> vec_A(
164  Teuchos::View, A.values(), 1, A.numRows()*A.numCols(), 1);
165  return vec_A.normInf();
166  }
167 
169 
173  template <typename ordinal_type, typename scalar_type>
174  void QR_CGS(
175  ordinal_type k,
176  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
177  const Teuchos::Array<scalar_type>& w,
178  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
179  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R)
180  {
181  using Teuchos::getCol;
182  typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type> SDV;
183  typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> SDM;
184 
185  // getCol() requires non-const SDM
186  SDM& Anc = const_cast<SDM&>(A);
187 
188  // Make sure Q is the right size
189  ordinal_type m = A.numRows();
190  if (Q.numRows() != m || Q.numCols() != k)
191  Q.shape(m,k);
192  if (R.numRows() != k || R.numCols() != k)
193  R.shape(k,k);
194 
195  for (ordinal_type j=0; j<k; j++) {
196  SDV Aj = getCol(Teuchos::View, Anc, j);
197  SDV Qj = getCol(Teuchos::View, Q, j);
198  Qj.assign(Aj);
199  for (ordinal_type i=0; i<j; i++) {
200  SDV Qi = getCol(Teuchos::View, Q, i);
201  R(i,j) = detail::weighted_inner_product(Qi, Aj, w);
202  detail::saxpy(1.0, Qj, -R(i,j), Qi); // Q(:,j) = 1.0*Q(:,j) - R(i,j)*Q(:,i)
203  }
204  R(j,j) = std::sqrt(detail::weighted_inner_product(Qj, Qj, w));
205  Qj.scale(1.0/R(j,j));
206  }
207 
208  }
209 
211 
215  template <typename ordinal_type, typename scalar_type>
216  void QR_MGS(
217  ordinal_type k,
218  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
219  const Teuchos::Array<scalar_type>& w,
220  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
221  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R)
222  {
223  using Teuchos::getCol;
224  typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type> SDV;
225  typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> SDM;
226 
227  // getCol() requires non-const SDM
228  SDM& Anc = const_cast<SDM&>(A);
229 
230  // Make sure Q is the right size
231  ordinal_type m = A.numRows();
232  if (Q.numRows() != m || Q.numCols() != k)
233  Q.shape(m,k);
234  if (R.numRows() != k || R.numCols() != k)
235  R.shape(k,k);
236 
237  for (ordinal_type i=0; i<k; i++) {
238  SDV Ai = getCol(Teuchos::View, Anc, i);
239  SDV Qi = getCol(Teuchos::View, Q, i);
240  Qi.assign(Ai);
241  }
242  for (ordinal_type i=0; i<k; i++) {
243  SDV Qi = getCol(Teuchos::View, Q, i);
244  for (ordinal_type j=0; j<i; j++) {
245  SDV Qj = getCol(Teuchos::View, Q, j);
247  detail::saxpy(1.0, Qi, -v, Qj); // Q(:,i) = 1.0*Q(:,i) - v*Q(:,j)
248  R(j,i) += v;
249  }
250  R(i,i) = std::sqrt(detail::weighted_inner_product(Qi, Qi, w));
251  Qi.scale(1.0/R(i,i));
252  }
253  }
254 
256 
260  template <typename ordinal_type, typename scalar_type>
261  void QR_MGS2(
262  ordinal_type k,
263  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
264  const Teuchos::Array<scalar_type>& w,
265  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
266  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R)
267  {
268  using Teuchos::getCol;
269  typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type> SDV;
270  typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> SDM;
271 
272  // getCol() requires non-const SDM
273  SDM& Anc = const_cast<SDM&>(A);
274 
275  // Make sure Q is the right size
276  ordinal_type m = A.numRows();
277  if (Q.numRows() != m || Q.numCols() != k)
278  Q.shape(m,k);
279  if (R.numRows() != k || R.numCols() != k)
280  R.shape(k,k);
281 
282  for (ordinal_type i=0; i<k; i++) {
283  SDV Ai = getCol(Teuchos::View, Anc, i);
284  SDV Qi = getCol(Teuchos::View, Q, i);
285  Qi.assign(Ai);
286  }
287  for (ordinal_type i=0; i<k; i++) {
288  SDV Qi = getCol(Teuchos::View, Q, i);
289  for (ordinal_type j=0; j<i; j++) {
290  SDV Qj = getCol(Teuchos::View, Q, j);
292  detail::saxpy(1.0, Qi, -v, Qj); // Q(:,i) = 1.0*Q(:,i) - v*Q(:,j)
293  R(j,i) += v;
294  }
295  for (ordinal_type j=i-1; j>=0; j--) {
296  SDV Qj = getCol(Teuchos::View, Q, j);
298  detail::saxpy(1.0, Qi, -v, Qj); // Q(:,i) = 1.0*Q(:,i) - v*Q(:,j)
299  R(j,i) += v;
300  }
301  R(i,i) = std::sqrt(detail::weighted_inner_product(Qi, Qi, w));
302  Qi.scale(1.0/R(i,i));
303  }
304  }
305 
307 
313  template <typename ordinal_type, typename scalar_type>
314  void
316  ordinal_type k,
317  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
318  const Teuchos::Array<scalar_type>& w,
319  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
320  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R)
321  {
323  ordinal_type m = A.numRows();
324  ordinal_type n = A.numCols();
325  ordinal_type kk = std::min(m,n);
326  if (k > kk)
327  k = kk;
328 
329  // Check that each component of w is 1
330  for (ordinal_type i=0; i<w.size(); i++)
331  TEUCHOS_TEST_FOR_EXCEPTION(
332  w[i] != 1.0, std::logic_error,
333  "CPQR_Householder_threshold() requires unit weight vector!");
334 
335  // Lapack routine overwrites A, so copy into temporary matrix
336  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
337  Teuchos::Copy, A, m, n);
338 
339  // QR
340  ordinal_type lda = AA.stride();
341  Teuchos::Array<scalar_type> tau(k);
342  Teuchos::Array<scalar_type> work(1);
343  ordinal_type info;
344  ordinal_type lwork = -1;
345  lapack.GEQRF(m, n, AA.values(), lda, &tau[0], &work[0], lwork, &info);
346  TEUCHOS_TEST_FOR_EXCEPTION(
347  info < 0, std::logic_error, "geqrf returned info = " << info);
348  lwork = work[0];
349  work.resize(lwork);
350  lapack.GEQRF(m, n, AA.values(), lda, &tau[0], &work[0], lwork, &info);
351  TEUCHOS_TEST_FOR_EXCEPTION(
352  info < 0, std::logic_error, "geqrf returned info = " << info);
353 
354  // Extract R
355  if (R.numRows() != k || R.numCols() != n)
356  R.shape(k,n);
357  R.putScalar(0.0);
358  for (ordinal_type i=0; i<k; i++)
359  for (ordinal_type j=i; j<n; j++)
360  R(i,j) = AA(i,j);
361 
362  // Extract Q
363  if (Q.numRows() != m || Q.numCols() != k)
364  Q.shape(m,k);
365  lwork = -1;
366  lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
367  TEUCHOS_TEST_FOR_EXCEPTION(
368  info < 0, std::logic_error, "orgqr returned info = " << info);
369  lwork = work[0];
370  work.resize(lwork);
371  lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
372  TEUCHOS_TEST_FOR_EXCEPTION(
373  info < 0, std::logic_error, "orgqr returned info = " << info);
374  if (Q.numRows() != m || Q.numCols() != k)
375  Q.shape(m, k);
376  for (ordinal_type i=0; i<m; i++)
377  for (ordinal_type j=0; j<k; j++)
378  Q(i,j) = AA(i,j);
379  }
380 
382 
393  template <typename ordinal_type, typename scalar_type>
394  void
396  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
397  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
398  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
399  Teuchos::Array<ordinal_type>& piv)
400  {
402  ordinal_type m = A.numRows();
403  ordinal_type n = A.numCols();
404  ordinal_type k = std::min(m,n);
405 
406  // Lapack routine overwrites A, so copy into temporary matrix
407  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
408  Teuchos::Copy, A, m, n);
409  if (Q.numRows() != m || Q.numCols() != k)
410  Q.shape(m,k);
411 
412  // Teuchos LAPACK interface doesn't have dgeqpf, so call it directly
413 
414  // Column pivoted QR
415  ordinal_type lda = AA.stride();
416  piv.resize(n);
417  Teuchos::Array<scalar_type> tau(k);
418  Teuchos::Array<scalar_type> work(3*n);
419  ordinal_type info;
420  DGEQPF_F77(&m, &n, AA.values(), &lda, &piv[0], &tau[0], &work[0], &info);
421  TEUCHOS_TEST_FOR_EXCEPTION(
422  info < 0, std::logic_error, "dgeqp3 returned info = " << info);
423 
424  // Extract R
425  if (R.numRows() != k || R.numCols() != n)
426  R.shape(k,n);
427  R.putScalar(0.0);
428  for (ordinal_type i=0; i<k; i++)
429  for (ordinal_type j=i; j<n; j++)
430  R(i,j) = AA(i,j);
431 
432  // Extract Q
433  ordinal_type lwork = -1;
434  lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
435  TEUCHOS_TEST_FOR_EXCEPTION(
436  info < 0, std::logic_error, "orgqr returned info = " << info);
437  lwork = work[0];
438  work.resize(lwork);
439  lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
440  TEUCHOS_TEST_FOR_EXCEPTION(
441  info < 0, std::logic_error, "orgqr returned info = " << info);
442  if (Q.numRows() != m || Q.numCols() != k)
443  Q.shape(m, k);
444  for (ordinal_type i=0; i<m; i++)
445  for (ordinal_type j=0; j<k; j++)
446  Q(i,j) = AA(i,j);
447 
448  // Transform piv to zero-based indexing
449  for (ordinal_type i=0; i<n; i++)
450  piv[i] -= 1;
451  }
452 
454 
466  template <typename ordinal_type, typename scalar_type>
467  void
469  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
470  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
471  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
472  Teuchos::Array<ordinal_type>& piv)
473  {
475  ordinal_type m = A.numRows();
476  ordinal_type n = A.numCols();
477  ordinal_type k = std::min(m,n);
478 
479  // Lapack routine overwrites A, so copy into temporary matrix
480  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
481  Teuchos::Copy, A, m, n);
482  if (Q.numRows() != m || Q.numCols() != k)
483  Q.shape(m,k);
484 
485  // Teuchos LAPACK interface doesn't have dgeqp3, so call it directly
486 
487  // Workspace query
488  ordinal_type lda = AA.stride();
489  piv.resize(n);
490  Teuchos::Array<scalar_type> tau(k);
491  Teuchos::Array<scalar_type> work(1);
492  ordinal_type lwork = -1;
493  ordinal_type info;
494  DGEQP3_F77(&m, &n, AA.values(), &lda, &piv[0], &tau[0], &work[0], &lwork,
495  &info);
496  TEUCHOS_TEST_FOR_EXCEPTION(
497  info < 0, std::logic_error, "dgeqp3 returned info = " << info);
498 
499  // Column pivoted QR
500  lwork = work[0];
501  work.resize(lwork);
502  DGEQP3_F77(&m, &n, AA.values(), &lda, &piv[0], &tau[0], &work[0], &lwork,
503  &info);
504  TEUCHOS_TEST_FOR_EXCEPTION(
505  info < 0, std::logic_error, "dgeqp3 returned info = " << info);
506 
507  // Extract R
508  if (R.numRows() != k || R.numCols() != n)
509  R.shape(k,n);
510  R.putScalar(0.0);
511  for (ordinal_type i=0; i<k; i++)
512  for (ordinal_type j=i; j<n; j++)
513  R(i,j) = AA(i,j);
514 
515  // Extract Q
516  lwork = -1;
517  lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
518  TEUCHOS_TEST_FOR_EXCEPTION(
519  info < 0, std::logic_error, "orgqr returned info = " << info);
520  lwork = work[0];
521  work.resize(lwork);
522  lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
523  TEUCHOS_TEST_FOR_EXCEPTION(
524  info < 0, std::logic_error, "orgqr returned info = " << info);
525  if (Q.numRows() != m || Q.numCols() != k)
526  Q.shape(m, k);
527  for (ordinal_type i=0; i<m; i++)
528  for (ordinal_type j=0; j<k; j++)
529  Q(i,j) = AA(i,j);
530 
531  // Transform piv to zero-based indexing
532  for (ordinal_type i=0; i<n; i++)
533  piv[i] -= 1;
534  }
535 
537 
555  template <typename ordinal_type, typename scalar_type>
558  const scalar_type& rank_threshold,
559  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
560  const Teuchos::Array<scalar_type>& w,
561  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
562  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
563  Teuchos::Array<ordinal_type>& piv)
564  {
565  // Check that each component of w is 1
566  for (ordinal_type i=0; i<w.size(); i++)
567  TEUCHOS_TEST_FOR_EXCEPTION(
568  w[i] != 1.0, std::logic_error,
569  "CPQR_Householder_threshold() requires unit weight vector!");
570 
571  // Compute full QR
572  CPQR_Householder(A, Q, R, piv);
573 
574  // Find leading block of R such that cond(R) <= 1/tau
575  ordinal_type rank = 0;
576  ordinal_type m = R.numRows();
577  scalar_type r_max = std::abs(R(rank,rank));
578  scalar_type r_min = std::abs(R(rank,rank));
579  for (rank=1; rank<m; rank++) {
580  if (std::abs(R(rank,rank)) > r_max)
581  r_max = std::abs(R(rank,rank));
582  if (std::abs(R(rank,rank)) < r_min)
583  r_min = std::abs(R(rank,rank));
584  if (r_min / r_max < rank_threshold)
585  break;
586  }
587 
588  // Extract blocks from Q and R
589  R.reshape(rank,rank);
590  Q.reshape(Q.numRows(), rank);
591 
592  return rank;
593  }
594 
596 
607  template <typename ordinal_type, typename scalar_type>
610  const scalar_type& rank_threshold,
611  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
612  const Teuchos::Array<scalar_type>& w,
613  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
614  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
615  Teuchos::Array<ordinal_type>& piv)
616  {
617  using Teuchos::getCol;
618  typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type> SDV;
619  ordinal_type m = A.numRows();
620  ordinal_type n = A.numCols();
621  ordinal_type k = std::min(m,n);
622 
623  // Make sure Q is the right size
624  if (Q.numRows() != m || Q.numCols() != n)
625  Q.shape(m,n);
626  if (R.numRows() != m || R.numCols() != n)
627  R.shape(m,n);
628  if (piv.size() != n)
629  piv.resize(n);
630  Q.assign(A);
631 
632  // Compute column norms
633  SDV nrm(n);
634  for (ordinal_type j=0; j<n; j++) {
635  SDV Qj = getCol(Teuchos::View, Q, j);
636  nrm(j) = detail::weighted_inner_product(Qj, Qj, w);
637  }
638  SDV Qtmp(m), Rtmp(m);
639 
640  Teuchos::Array<ordinal_type> piv_orig(piv);
641  for (ordinal_type i=0; i<n; i++)
642  piv[i] = i;
643 
644  // Prepivot any columns requested by setting piv[i] != 0
645  ordinal_type nfixed = 0;
646  for (ordinal_type j=0; j<n; j++) {
647  if (piv_orig[j] != 0) {
648  if (j != nfixed) {
649  scalar_type tmp = nrm(j);
650  nrm(j) = nrm(nfixed);
651  nrm(nfixed) = tmp;
652 
653  SDV Qpvt = getCol(Teuchos::View, Q, j);
654  SDV Qj = getCol(Teuchos::View, Q, nfixed);
655  Qtmp.assign(Qpvt);
656  Qpvt.assign(Qj);
657  Qj.assign(Qtmp);
658 
659  ordinal_type t = piv[j];
660  piv[j] = piv[nfixed];
661  piv[nfixed] = t;
662  }
663  ++nfixed;
664  }
665  }
666 
667  scalar_type sigma = 1.0 + rank_threshold;
668  scalar_type r_max = 0.0;
669  ordinal_type j=0;
670  R.putScalar(0.0);
671  while (j < k && sigma >= rank_threshold) {
672 
673  SDV Qj = getCol(Teuchos::View, Q, j);
674 
675  // Find pivot column if we are passed the pre-pivot stage
676  if (j >= nfixed) {
677  ordinal_type pvt = j;
678  for (ordinal_type l=j+1; l<n; l++)
679  if (nrm(l) > nrm(pvt))
680  pvt = l;
681 
682  // Interchange column j and pvt
683  if (pvt != j) {
684  SDV Qpvt = getCol(Teuchos::View, Q, pvt);
685  Qtmp.assign(Qpvt);
686  Qpvt.assign(Qj);
687  Qj.assign(Qtmp);
688 
689  SDV Rpvt = getCol(Teuchos::View, R, pvt);
690  SDV Rj = getCol(Teuchos::View, R, j);
691  Rtmp.assign(Rpvt);
692  Rpvt.assign(Rj);
693  Rj.assign(Rtmp);
694 
695  ordinal_type t = piv[pvt];
696  piv[pvt] = piv[j];
697  piv[j] = t;
698  }
699  }
700 
701  R(j,j) = std::sqrt(detail::weighted_inner_product(Qj, Qj, w));
702  Qj.scale(1.0/R(j,j));
703  for (ordinal_type l=j+1; l<n; l++) {
704  SDV Ql = getCol(Teuchos::View, Q, l);
706  R(j,l) = t;
707  detail::saxpy(1.0, Ql, -t, Qj);
708 
709  // Update norms
710  nrm(l) = detail::weighted_inner_product(Ql, Ql, w);
711  }
712 
713  if (std::abs(R(j,j)) > r_max)
714  r_max = R(j,j);
715  sigma = std::abs(R(j,j)/r_max);
716  j++;
717  }
718 
719  ordinal_type rank = j;
720  if (sigma < rank_threshold)
721  rank--;
722  Q.reshape(m, rank);
723  R.reshape(rank, rank);
724 
725  return rank;
726  }
727 
729 
740  template <typename ordinal_type, typename scalar_type>
743  const scalar_type& rank_threshold,
744  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
745  const Teuchos::Array<scalar_type>& w,
746  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
747  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
748  Teuchos::Array<ordinal_type>& piv)
749  {
750  // First do standard column-pivoted QR
751  ordinal_type rank = CPQR_MGS_threshold(rank_threshold, A, w, Q, R, piv);
752 
753  // Now apply standard MGS to Q
754  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> A2(Q), R2;
755  QR_MGS(rank, A2, w, Q, R2);
756 
757  return rank;
758  }
759 
761  template <typename ordinal_type, typename scalar_type>
763  cond_R(const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R)
764  {
765  ordinal_type k = R.numRows();
766  if (R.numCols() < k)
767  k = R.numCols();
768  scalar_type r_max = std::abs(R(0,0));
769  scalar_type r_min = std::abs(R(0,0));
770  for (ordinal_type i=1; i<k; i++) {
771  if (std::abs(R(i,i)) > r_max)
772  r_max = R(i,i);
773  if (std::abs(R(i,i)) < r_min)
774  r_min = R(i,i);
775  }
776  scalar_type cond_r = r_max / r_min;
777  return cond_r;
778  }
779 
781 
785  template <typename ordinal_type, typename scalar_type>
788  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
789  const Teuchos::Array<scalar_type>& w)
790  {
791  ordinal_type m = Q.numRows();
792  ordinal_type n = Q.numCols();
793  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> Qt(m,n);
794  for (ordinal_type i=0; i<m; i++)
795  for (ordinal_type j=0; j<n; j++)
796  Qt(i,j) = w[i]*Q(i,j);
797  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(n,n);
798  err.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q, Qt, 0.0);
799  for (ordinal_type i=0; i<n; i++)
800  err(i,i) -= 1.0;
801  return err.normInf();
802  }
803 
805 
808  template <typename ordinal_type, typename scalar_type>
811  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q)
812  {
813  ordinal_type n = Q.numCols();
814  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(n,n);
815  err.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q, Q, 0.0);
816  for (ordinal_type i=0; i<n; i++)
817  err(i,i) -= 1.0;
818  return err.normInf();
819  }
820 
822 
827  template <typename ordinal_type, typename scalar_type>
830  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
831  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
832  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R)
833  {
834  ordinal_type k = Q.numCols();
835  ordinal_type m = Q.numRows();
836  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> AA(
837  Teuchos::View, A, m, k, 0, 0);
838  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(m,k);
839  ordinal_type ret =
840  err.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
841  TEUCHOS_ASSERT(ret == 0);
842  err -= AA;
843  return err.normInf();
844  }
845 
847 
853  template <typename ordinal_type, typename scalar_type>
856  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
857  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
858  const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R,
859  const Teuchos::Array<ordinal_type>& piv)
860  {
861  ordinal_type k = Q.numCols();
862  ordinal_type m = Q.numRows();
863  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> AP(m, k);
864  for (ordinal_type j=0; j<k; j++)
865  for (ordinal_type i=0; i<m; i++)
866  AP(i,j) = A(i,piv[j]);
867  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(m,k);
868  ordinal_type ret =
869  err.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
870  TEUCHOS_ASSERT(ret == 0);
871  err -= AP;
872 
873  return err.normInf();
874  }
875 
877 
880  template <typename ordinal_type, typename scalar_type>
881  void
882  svd(const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
883  Teuchos::Array<scalar_type>& s,
884  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& U,
885  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Vt)
886  {
888  ordinal_type m = A.numRows();
889  ordinal_type n = A.numCols();
890  ordinal_type k = std::min(m,n);
891  ordinal_type lda = A.stride();
892 
893  // Copy A since GESVD overwrites it (always)
894  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
895  Teuchos::Copy, A, m, n);
896 
897  // Resize appropriately
898  if (U.numRows() != m || U.numCols() != m)
899  U.shape(m,m);
900  if (Vt.numRows() != n || Vt.numCols() != n)
901  Vt.shape(n,n);
902  if (s.size() != k)
903  s.resize(k);
904  ordinal_type ldu = U.stride();
905  ordinal_type ldv = Vt.stride();
906 
907  // Workspace query
908  Teuchos::Array<scalar_type> work(1);
909  ordinal_type lwork = -1;
910  ordinal_type info;
911  lapack.GESVD('A', 'A', m, n, AA.values(), lda, &s[0], U.values(), ldu,
912  Vt.values(), ldv, &work[0], lwork, NULL, &info);
913  TEUCHOS_TEST_FOR_EXCEPTION(
914  info < 0, std::logic_error, "dgesvd returned info = " << info);
915 
916  // Do SVD
917  lwork = work[0];
918  work.resize(lwork);
919  lapack.GESVD('A', 'A', m, n, AA.values(), lda, &s[0], U.values(), ldu,
920  Vt.values(), ldv, &work[0], lwork, NULL, &info);
921  TEUCHOS_TEST_FOR_EXCEPTION(
922  info < 0, std::logic_error, "dgesvd returned info = " << info);
923 
924  }
925 
926  template <typename ordinal_type, typename scalar_type>
929  const scalar_type& rank_threshold,
930  const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
931  Teuchos::Array<scalar_type>& s,
932  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& U,
933  Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Vt)
934  {
935  // Compute full SVD
936  svd(A, s, U, Vt);
937 
938  // Find leading block of S such that cond(S) <= 1/tau
939  ordinal_type rank = 0;
940  ordinal_type m = s.size();
941  while (rank < m && s[rank]/s[0] > rank_threshold) rank++;
942 
943  // Extract blocks from s, U and Vt
944  s.resize(rank);
945  U.reshape(U.numRows(),rank);
946  Vt.reshape(rank, Vt.numCols());
947 
948  return rank;
949  }
950 
951 }
952 
953 #endif
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.
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
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.
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.
void print_matlab(std::ostream &os, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A)
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.
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.
scalar_type cond_R(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R)
Compute condition number of upper-triangular R.
#define DGEQP3_F77
#define DGEQPF_F77
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.
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.
Teuchos::SerialDenseMatrix< int, pce_type > SDM
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)
KOKKOS_INLINE_FUNCTION PCE< Storage > min(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
Teuchos::SerialDenseVector< int, pce_type > SDV
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.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:260
scalar_type weightedQROrthogonalizationError(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, const Teuchos::Array< scalar_type > &w)
Compute weighted QR orthogonalization error.
Top-level namespace for Stokhos classes and functions.
Specialization for Sacado::UQ::PCE< Storage<...> >
scalar_type QROrthogonalizationError(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q)
Compute QR orthogonalization error.
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
expr expr expr expr j
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.
scalar_type vec_norm_inf(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A)
Vector-infinity norm of a matrix.
const IndexType const IndexType const IndexType const IndexType * Aj
Definition: csr_vector.h:260
void CPQR_Householder(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.
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.
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.
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.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
Definition: csr_vector.h:267
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.