42#ifndef STOKHOS_SDM_UTILS_HPP
43#define STOKHOS_SDM_UTILS_HPP
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"
52#define DGEQP3_F77 F77_BLAS_MANGLE(dgeqp3,DGEQP3)
55#if defined(HAVE_STOKHOS_MKL)
56 #include "mkl_version.h"
57 #if __INTEL_MKL__ >= 2021
58 #define MKL_NO_EXCEPT noexcept
66void DGEQP3_F77(
const int*,
const int*,
double*,
const int*,
int*,
72#ifdef HAVE_STOKHOS_MATLABLIB
84 template <
typename ordinal_type,
typename scalar_type>
86 const Teuchos::SerialDenseVector<ordinal_type,scalar_type>& x,
87 const Teuchos::SerialDenseVector<ordinal_type,scalar_type>& y,
88 const Teuchos::Array<scalar_type>& w)
90 ordinal_type n = x.length();
92 for (ordinal_type i=0; i<n; i++)
98 template <
typename ordinal_type,
typename scalar_type>
100 const scalar_type& alpha,
101 Teuchos::SerialDenseVector<ordinal_type, scalar_type>& x,
102 const scalar_type& beta,
103 const Teuchos::SerialDenseVector<ordinal_type, scalar_type>& y)
105 ordinal_type n = x.length();
106 for (ordinal_type i=0; i<n; i++)
107 x[i] = alpha*x[i] + beta*y[i];
113 template <
typename ordinal_type,
typename scalar_type>
116 const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A)
119 for (ordinal_type i=0; i<A.numRows(); i++) {
120 for (ordinal_type
j=0;
j<A.numCols();
j++)
122 if (i < A.numRows()-1)
123 os <<
";" << std::endl <<
" ";
125 os <<
"];" << std::endl;
128#ifdef HAVE_STOKHOS_MATLABLIB
130 template <
typename ordinal_type>
132 save_matlab(
const std::string& file_name,
const std::string& matrix_name,
133 const Teuchos::SerialDenseMatrix<ordinal_type,double>& A,
134 bool overwrite =
true)
142 MATFile *file = matOpen(file_name.c_str(), mode);
143 TEUCHOS_ASSERT(file != NULL);
146 mxArray *mat = mxCreateDoubleMatrix(0, 0, mxREAL);
147 TEUCHOS_ASSERT(mat != NULL);
148 mxSetPr(mat, A.values());
149 mxSetM(mat, A.numRows());
150 mxSetN(mat, A.numCols());
153 ordinal_type ret = matPutVariable(file, matrix_name.c_str(), mat);
154 TEUCHOS_ASSERT(ret == 0);
157 ret = matClose(file);
158 TEUCHOS_ASSERT(ret == 0);
169 template <
typename ordinal_type,
typename scalar_type>
171 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A) {
172 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> vec_A(
173 Teuchos::View, A.values(), 1, A.numRows()*A.numCols(), 1);
174 return vec_A.normInf();
182 template <
typename ordinal_type,
typename scalar_type>
185 const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
186 const Teuchos::Array<scalar_type>& w,
187 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
188 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R)
190 using Teuchos::getCol;
191 typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type>
SDV;
192 typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>
SDM;
195 SDM& Anc =
const_cast<SDM&
>(A);
198 ordinal_type m = A.numRows();
199 if (Q.numRows() != m || Q.numCols() != k)
201 if (R.numRows() != k || R.numCols() != k)
204 for (ordinal_type
j=0;
j<k;
j++) {
205 SDV Aj = getCol(Teuchos::View, Anc,
j);
206 SDV Qj = getCol(Teuchos::View, Q,
j);
208 for (ordinal_type i=0; i<
j; i++) {
209 SDV Qi = getCol(Teuchos::View, Q, i);
214 Qj.scale(1.0/R(
j,
j));
224 template <
typename ordinal_type,
typename scalar_type>
227 const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
228 const Teuchos::Array<scalar_type>& w,
229 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
230 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R)
232 using Teuchos::getCol;
233 typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type>
SDV;
234 typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>
SDM;
237 SDM& Anc =
const_cast<SDM&
>(A);
240 ordinal_type m = A.numRows();
241 if (Q.numRows() != m || Q.numCols() != k)
243 if (R.numRows() != k || R.numCols() != k)
246 for (ordinal_type i=0; i<k; i++) {
247 SDV Ai = getCol(Teuchos::View, Anc, i);
248 SDV Qi = getCol(Teuchos::View, Q, i);
251 for (ordinal_type i=0; i<k; i++) {
252 SDV Qi = getCol(Teuchos::View, Q, i);
253 for (ordinal_type
j=0;
j<i;
j++) {
254 SDV Qj = getCol(Teuchos::View, Q,
j);
260 Qi.scale(1.0/R(i,i));
269 template <
typename ordinal_type,
typename scalar_type>
272 const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
273 const Teuchos::Array<scalar_type>& w,
274 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
275 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R)
277 using Teuchos::getCol;
278 typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type>
SDV;
279 typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>
SDM;
282 SDM& Anc =
const_cast<SDM&
>(A);
285 ordinal_type m = A.numRows();
286 if (Q.numRows() != m || Q.numCols() != k)
288 if (R.numRows() != k || R.numCols() != k)
291 for (ordinal_type i=0; i<k; i++) {
292 SDV Ai = getCol(Teuchos::View, Anc, i);
293 SDV Qi = getCol(Teuchos::View, Q, i);
296 for (ordinal_type i=0; i<k; i++) {
297 SDV Qi = getCol(Teuchos::View, Q, i);
298 for (ordinal_type
j=0;
j<i;
j++) {
299 SDV Qj = getCol(Teuchos::View, Q,
j);
304 for (ordinal_type
j=i-1;
j>=0;
j--) {
305 SDV Qj = getCol(Teuchos::View, Q,
j);
311 Qi.scale(1.0/R(i,i));
322 template <
typename ordinal_type,
typename scalar_type>
326 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
327 const Teuchos::Array<scalar_type>& w,
328 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
329 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R)
332 ordinal_type m = A.numRows();
333 ordinal_type n = A.numCols();
334 ordinal_type kk = std::min(m,n);
339 for (ordinal_type i=0; i<w.size(); i++)
340 TEUCHOS_TEST_FOR_EXCEPTION(
341 w[i] != 1.0, std::logic_error,
342 "CPQR_Householder_threshold() requires unit weight vector!");
345 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
346 Teuchos::Copy, A, m, n);
349 ordinal_type lda = AA.stride();
350 Teuchos::Array<scalar_type> tau(k);
351 Teuchos::Array<scalar_type> work(1);
353 ordinal_type lwork = -1;
354 lapack.GEQRF(m, n, AA.values(), lda, &tau[0], &work[0], lwork, &info);
355 TEUCHOS_TEST_FOR_EXCEPTION(
356 info < 0, std::logic_error,
"geqrf returned info = " << info);
359 lapack.GEQRF(m, n, AA.values(), lda, &tau[0], &work[0], lwork, &info);
360 TEUCHOS_TEST_FOR_EXCEPTION(
361 info < 0, std::logic_error,
"geqrf returned info = " << info);
364 if (R.numRows() != k || R.numCols() != n)
367 for (ordinal_type i=0; i<k; i++)
368 for (ordinal_type
j=i;
j<n;
j++)
372 if (Q.numRows() != m || Q.numCols() != k)
375 lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
376 TEUCHOS_TEST_FOR_EXCEPTION(
377 info < 0, std::logic_error,
"orgqr returned info = " << info);
380 lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
381 TEUCHOS_TEST_FOR_EXCEPTION(
382 info < 0, std::logic_error,
"orgqr returned info = " << info);
383 if (Q.numRows() != m || Q.numCols() != k)
385 for (ordinal_type i=0; i<m; i++)
386 for (ordinal_type
j=0;
j<k;
j++)
403 template <
typename ordinal_type,
typename scalar_type>
406 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
407 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
408 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
409 Teuchos::Array<ordinal_type>& piv)
412 ordinal_type m = A.numRows();
413 ordinal_type n = A.numCols();
414 ordinal_type k = std::min(m,n);
417 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
418 Teuchos::Copy, A, m, n);
419 if (Q.numRows() != m || Q.numCols() != k)
425 ordinal_type lda = AA.stride();
427 Teuchos::Array<scalar_type> tau(k);
428 Teuchos::Array<scalar_type> work(1);
429 ordinal_type lwork = -1;
431 DGEQP3_F77(&m, &n, AA.values(), &lda, &piv[0], &tau[0], &work[0], &lwork,
433 TEUCHOS_TEST_FOR_EXCEPTION(
434 info < 0, std::logic_error,
"dgeqp3 returned info = " << info);
439 DGEQP3_F77(&m, &n, AA.values(), &lda, &piv[0], &tau[0], &work[0], &lwork,
441 TEUCHOS_TEST_FOR_EXCEPTION(
442 info < 0, std::logic_error,
"dgeqp3 returned info = " << info);
445 if (R.numRows() != k || R.numCols() != n)
448 for (ordinal_type i=0; i<k; i++)
449 for (ordinal_type
j=i;
j<n;
j++)
454 lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
455 TEUCHOS_TEST_FOR_EXCEPTION(
456 info < 0, std::logic_error,
"orgqr returned info = " << info);
459 lapack.ORGQR(m, k, k, AA.values(), lda, &tau[0], &work[0], lwork, &info);
460 TEUCHOS_TEST_FOR_EXCEPTION(
461 info < 0, std::logic_error,
"orgqr returned info = " << info);
462 if (Q.numRows() != m || Q.numCols() != k)
464 for (ordinal_type i=0; i<m; i++)
465 for (ordinal_type
j=0;
j<k;
j++)
469 for (ordinal_type i=0; i<n; i++)
492 template <
typename ordinal_type,
typename scalar_type>
495 const scalar_type& rank_threshold,
496 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
497 const Teuchos::Array<scalar_type>& w,
498 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
499 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
500 Teuchos::Array<ordinal_type>& piv)
503 for (ordinal_type i=0; i<w.size(); i++)
504 TEUCHOS_TEST_FOR_EXCEPTION(
505 w[i] != 1.0, std::logic_error,
506 "CPQR_Householder_threshold() requires unit weight vector!");
512 ordinal_type rank = 0;
513 ordinal_type m = R.numRows();
514 scalar_type r_max = std::abs(R(rank,rank));
515 scalar_type r_min = std::abs(R(rank,rank));
516 for (rank=1; rank<m; rank++) {
517 if (std::abs(R(rank,rank)) > r_max)
518 r_max = std::abs(R(rank,rank));
519 if (std::abs(R(rank,rank)) < r_min)
520 r_min = std::abs(R(rank,rank));
521 if (r_min / r_max < rank_threshold)
526 R.reshape(rank,rank);
527 Q.reshape(Q.numRows(), rank);
544 template <
typename ordinal_type,
typename scalar_type>
547 const scalar_type& rank_threshold,
548 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
549 const Teuchos::Array<scalar_type>& w,
550 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
551 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
552 Teuchos::Array<ordinal_type>& piv)
554 using Teuchos::getCol;
555 typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type>
SDV;
556 ordinal_type m = A.numRows();
557 ordinal_type n = A.numCols();
558 ordinal_type k = std::min(m,n);
561 if (Q.numRows() != m || Q.numCols() != n)
563 if (R.numRows() != m || R.numCols() != n)
571 for (ordinal_type
j=0;
j<n;
j++) {
572 SDV Qj = getCol(Teuchos::View, Q,
j);
575 SDV Qtmp(m), Rtmp(m);
577 Teuchos::Array<ordinal_type> piv_orig(piv);
578 for (ordinal_type i=0; i<n; i++)
582 ordinal_type nfixed = 0;
583 for (ordinal_type
j=0;
j<n;
j++) {
584 if (piv_orig[
j] != 0) {
586 scalar_type tmp = nrm(
j);
587 nrm(
j) = nrm(nfixed);
590 SDV Qpvt = getCol(Teuchos::View, Q,
j);
591 SDV Qj = getCol(Teuchos::View, Q, nfixed);
596 ordinal_type t = piv[
j];
597 piv[
j] = piv[nfixed];
604 scalar_type sigma = 1.0 + rank_threshold;
605 scalar_type r_max = 0.0;
608 while (j < k && sigma >= rank_threshold) {
610 SDV Qj = getCol(Teuchos::View, Q,
j);
614 ordinal_type pvt =
j;
615 for (ordinal_type l=
j+1; l<n; l++)
616 if (nrm(l) > nrm(pvt))
621 SDV Qpvt = getCol(Teuchos::View, Q, pvt);
626 SDV Rpvt = getCol(Teuchos::View, R, pvt);
627 SDV Rj = getCol(Teuchos::View, R,
j);
632 ordinal_type t = piv[pvt];
639 Qj.scale(1.0/R(
j,
j));
640 for (ordinal_type l=
j+1; l<n; l++) {
641 SDV Ql = getCol(Teuchos::View, Q, l);
650 if (std::abs(R(
j,
j)) > r_max)
652 sigma = std::abs(R(
j,
j)/r_max);
656 ordinal_type rank =
j;
657 if (sigma < rank_threshold)
660 R.reshape(rank, rank);
677 template <
typename ordinal_type,
typename scalar_type>
680 const scalar_type& rank_threshold,
681 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
682 const Teuchos::Array<scalar_type>& w,
683 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
684 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R,
685 Teuchos::Array<ordinal_type>& piv)
691 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> A2(Q), R2;
692 QR_MGS(rank, A2, w, Q, R2);
698 template <
typename ordinal_type,
typename scalar_type>
700 cond_R(
const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R)
702 ordinal_type k = R.numRows();
705 scalar_type r_max = std::abs(R(0,0));
706 scalar_type r_min = std::abs(R(0,0));
707 for (ordinal_type i=1; i<k; i++) {
708 if (std::abs(R(i,i)) > r_max)
710 if (std::abs(R(i,i)) < r_min)
713 scalar_type cond_r = r_max / r_min;
722 template <
typename ordinal_type,
typename scalar_type>
725 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
726 const Teuchos::Array<scalar_type>& w)
728 ordinal_type m = Q.numRows();
729 ordinal_type n = Q.numCols();
730 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> Qt(m,n);
731 for (ordinal_type i=0; i<m; i++)
732 for (ordinal_type
j=0;
j<n;
j++)
733 Qt(i,
j) = w[i]*Q(i,
j);
734 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(n,n);
735 err.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q, Qt, 0.0);
736 for (ordinal_type i=0; i<n; i++)
738 return err.normInf();
745 template <
typename ordinal_type,
typename scalar_type>
748 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q)
750 ordinal_type n = Q.numCols();
751 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(n,n);
752 err.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q, Q, 0.0);
753 for (ordinal_type i=0; i<n; i++)
755 return err.normInf();
764 template <
typename ordinal_type,
typename scalar_type>
767 const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
768 const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
769 const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R)
771 ordinal_type k = Q.numCols();
772 ordinal_type m = Q.numRows();
773 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> AA(
774 Teuchos::View, A, m, k, 0, 0);
775 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(m,k);
777 err.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
778 TEUCHOS_ASSERT(ret == 0);
780 return err.normInf();
790 template <
typename ordinal_type,
typename scalar_type>
793 const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A,
794 const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& Q,
795 const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& R,
796 const Teuchos::Array<ordinal_type>& piv)
798 ordinal_type k = Q.numCols();
799 ordinal_type m = Q.numRows();
800 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> AP(m, k);
801 for (ordinal_type
j=0;
j<k;
j++)
802 for (ordinal_type i=0; i<m; i++)
803 AP(i,
j) = A(i,piv[
j]);
804 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(m,k);
806 err.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
807 TEUCHOS_ASSERT(ret == 0);
810 return err.normInf();
817 template <
typename ordinal_type,
typename scalar_type>
819 svd(
const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
820 Teuchos::Array<scalar_type>& s,
821 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& U,
822 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Vt)
825 ordinal_type m = A.numRows();
826 ordinal_type n = A.numCols();
827 ordinal_type k = std::min(m,n);
828 ordinal_type lda = A.stride();
831 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
832 Teuchos::Copy, A, m, n);
835 if (U.numRows() != m || U.numCols() != m)
837 if (Vt.numRows() != n || Vt.numCols() != n)
841 ordinal_type ldu = U.stride();
842 ordinal_type ldv = Vt.stride();
845 Teuchos::Array<scalar_type> work(1);
846 ordinal_type lwork = -1;
848 lapack.GESVD(
'A',
'A', m, n, AA.values(), lda, &s[0], U.values(), ldu,
849 Vt.values(), ldv, &work[0], lwork, NULL, &info);
850 TEUCHOS_TEST_FOR_EXCEPTION(
851 info < 0, std::logic_error,
"dgesvd returned info = " << info);
856 lapack.GESVD(
'A',
'A', m, n, AA.values(), lda, &s[0], U.values(), ldu,
857 Vt.values(), ldv, &work[0], lwork, NULL, &info);
858 TEUCHOS_TEST_FOR_EXCEPTION(
859 info < 0, std::logic_error,
"dgesvd returned info = " << info);
863 template <
typename ordinal_type,
typename scalar_type>
866 const scalar_type& rank_threshold,
867 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
868 Teuchos::Array<scalar_type>& s,
869 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& U,
870 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Vt)
876 ordinal_type rank = 0;
877 ordinal_type m = s.size();
878 while (rank < m && s[rank]/s[0] > rank_threshold) rank++;
882 U.reshape(U.numRows(),rank);
883 Vt.reshape(rank, Vt.numCols());
Specialization for Sacado::UQ::PCE< Storage<...> >
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 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.
Top-level namespace for Stokhos classes and functions.
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.
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.
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.
scalar_type vec_norm_inf(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A)
Vector-infinity norm of a matrix.
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)
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.
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.
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.
void print_matlab(std::ostream &os, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A)
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.
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 cond_R(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R)
Compute condition number of upper-triangular R.
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.
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_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.
scalar_type QROrthogonalizationError(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q)
Compute QR orthogonalization error.
scalar_type weightedQROrthogonalizationError(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, const Teuchos::Array< scalar_type > &w)
Compute weighted QR orthogonalization error.
Teuchos::SerialDenseMatrix< int, pce_type > SDM
Teuchos::SerialDenseVector< int, pce_type > SDV