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 DGEQPF_F77 F77_BLAS_MANGLE(dgeqpf,DGEQPF) 53 #define DGEQP3_F77 F77_BLAS_MANGLE(dgeqp3,DGEQP3) 55 void DGEQPF_F77(
const int*,
const int*,
double*,
const int*,
int*,
double*,
57 void DGEQP3_F77(
const int*,
const int*,
double*,
const int*,
int*,
58 double*,
double*,
const int*,
int*);
63 #ifdef HAVE_STOKHOS_MATLABLIB 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)
89 template <
typename ordinal_type,
typename scalar_type>
92 Teuchos::SerialDenseVector<ordinal_type, scalar_type>&
x,
94 const Teuchos::SerialDenseVector<ordinal_type, scalar_type>&
y)
98 x[i] = alpha*
x[i] + beta*
y[i];
104 template <
typename ordinal_type,
typename scalar_type>
107 const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A)
113 if (i < A.numRows()-1)
114 os <<
";" << std::endl <<
" ";
116 os <<
"];" << std::endl;
119 #ifdef HAVE_STOKHOS_MATLABLIB 121 template <
typename ordinal_type>
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)
133 MATFile *file = matOpen(file_name.c_str(), mode);
134 TEUCHOS_ASSERT(file != NULL);
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());
144 ordinal_type ret = matPutVariable(file, matrix_name.c_str(), mat);
145 TEUCHOS_ASSERT(ret == 0);
148 ret = matClose(file);
149 TEUCHOS_ASSERT(ret == 0);
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();
173 template <
typename ordinal_type,
typename scalar_type>
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)
181 using Teuchos::getCol;
182 typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type>
SDV;
183 typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>
SDM;
186 SDM& Anc =
const_cast<SDM&
>(A);
190 if (Q.numRows() != m || Q.numCols() != k)
192 if (R.numRows() != k || R.numCols() != k)
196 SDV Aj = getCol(Teuchos::View, Anc,
j);
197 SDV Qj = getCol(Teuchos::View, Q,
j);
200 SDV Qi = getCol(Teuchos::View, Q, i);
205 Qj.scale(1.0/R(
j,
j));
215 template <
typename ordinal_type,
typename scalar_type>
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)
223 using Teuchos::getCol;
224 typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type>
SDV;
225 typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>
SDM;
228 SDM& Anc =
const_cast<SDM&
>(A);
232 if (Q.numRows() != m || Q.numCols() != k)
234 if (R.numRows() != k || R.numCols() != k)
238 SDV Ai = getCol(Teuchos::View, Anc, i);
239 SDV Qi = getCol(Teuchos::View, Q, i);
243 SDV Qi = getCol(Teuchos::View, Q, i);
245 SDV Qj = getCol(Teuchos::View, Q,
j);
251 Qi.scale(1.0/R(i,i));
260 template <
typename ordinal_type,
typename scalar_type>
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)
268 using Teuchos::getCol;
269 typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type>
SDV;
270 typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>
SDM;
273 SDM& Anc =
const_cast<SDM&
>(A);
277 if (Q.numRows() != m || Q.numCols() != k)
279 if (R.numRows() != k || R.numCols() != k)
283 SDV Ai = getCol(Teuchos::View, Anc, i);
284 SDV Qi = getCol(Teuchos::View, Q, i);
288 SDV Qi = getCol(Teuchos::View, Q, i);
290 SDV Qj = getCol(Teuchos::View, Q,
j);
296 SDV Qj = getCol(Teuchos::View, Q,
j);
302 Qi.scale(1.0/R(i,i));
313 template <
typename ordinal_type,
typename scalar_type>
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)
331 TEUCHOS_TEST_FOR_EXCEPTION(
332 w[i] != 1.0, std::logic_error,
333 "CPQR_Householder_threshold() requires unit weight vector!");
336 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
337 Teuchos::Copy, A, m, n);
341 Teuchos::Array<scalar_type> tau(k);
342 Teuchos::Array<scalar_type> work(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);
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);
355 if (R.numRows() != k || R.numCols() != n)
363 if (Q.numRows() != m || Q.numCols() != k)
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);
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)
393 template <
typename ordinal_type,
typename scalar_type>
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)
407 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
408 Teuchos::Copy, A, m, n);
409 if (Q.numRows() != m || Q.numCols() != k)
417 Teuchos::Array<scalar_type> tau(k);
418 Teuchos::Array<scalar_type> work(3*n);
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);
425 if (R.numRows() != k || R.numCols() != n)
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);
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)
466 template <
typename ordinal_type,
typename scalar_type>
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)
480 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
481 Teuchos::Copy, A, m, n);
482 if (Q.numRows() != m || Q.numCols() != k)
490 Teuchos::Array<scalar_type> tau(k);
491 Teuchos::Array<scalar_type> work(1);
494 DGEQP3_F77(&m, &n, AA.values(), &lda, &piv[0], &tau[0], &work[0], &lwork,
496 TEUCHOS_TEST_FOR_EXCEPTION(
497 info < 0, std::logic_error,
"dgeqp3 returned info = " << info);
502 DGEQP3_F77(&m, &n, AA.values(), &lda, &piv[0], &tau[0], &work[0], &lwork,
504 TEUCHOS_TEST_FOR_EXCEPTION(
505 info < 0, std::logic_error,
"dgeqp3 returned info = " << info);
508 if (R.numRows() != k || R.numCols() != n)
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);
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)
555 template <
typename ordinal_type,
typename scalar_type>
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)
567 TEUCHOS_TEST_FOR_EXCEPTION(
568 w[i] != 1.0, std::logic_error,
569 "CPQR_Householder_threshold() requires unit weight vector!");
579 for (rank=1; rank<m; rank++) {
584 if (r_min / r_max < rank_threshold)
589 R.reshape(rank,rank);
590 Q.reshape(Q.numRows(), rank);
607 template <
typename ordinal_type,
typename scalar_type>
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)
617 using Teuchos::getCol;
618 typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type>
SDV;
624 if (Q.numRows() != m || Q.numCols() != n)
626 if (R.numRows() != m || R.numCols() != n)
635 SDV Qj = getCol(Teuchos::View, Q,
j);
638 SDV Qtmp(m), Rtmp(m);
640 Teuchos::Array<ordinal_type> piv_orig(piv);
647 if (piv_orig[
j] != 0) {
650 nrm(
j) = nrm(nfixed);
653 SDV Qpvt = getCol(Teuchos::View, Q,
j);
654 SDV Qj = getCol(Teuchos::View, Q, nfixed);
660 piv[
j] = piv[nfixed];
671 while (j < k && sigma >= rank_threshold) {
673 SDV Qj = getCol(Teuchos::View, Q,
j);
679 if (nrm(l) > nrm(pvt))
684 SDV Qpvt = getCol(Teuchos::View, Q, pvt);
689 SDV Rpvt = getCol(Teuchos::View, R, pvt);
690 SDV Rj = getCol(Teuchos::View, R,
j);
702 Qj.scale(1.0/R(
j,
j));
704 SDV Ql = getCol(Teuchos::View, Q, l);
720 if (sigma < rank_threshold)
723 R.reshape(rank, rank);
740 template <
typename ordinal_type,
typename scalar_type>
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)
754 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> A2(Q), R2;
755 QR_MGS(rank, A2, w, Q, R2);
761 template <
typename ordinal_type,
typename scalar_type>
763 cond_R(
const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R)
785 template <
typename ordinal_type,
typename scalar_type>
788 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
789 const Teuchos::Array<scalar_type>& w)
793 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> Qt(m,n);
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);
801 return err.normInf();
808 template <
typename ordinal_type,
typename scalar_type>
811 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q)
814 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(n,n);
815 err.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q, Q, 0.0);
818 return err.normInf();
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)
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);
840 err.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
841 TEUCHOS_ASSERT(ret == 0);
843 return err.normInf();
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)
863 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> AP(m, k);
866 AP(i,
j) = A(i,piv[
j]);
867 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(m,k);
869 err.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
870 TEUCHOS_ASSERT(ret == 0);
873 return err.normInf();
880 template <
typename ordinal_type,
typename scalar_type>
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)
894 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
895 Teuchos::Copy, A, m, n);
898 if (U.numRows() != m || U.numCols() != m)
900 if (Vt.numRows() != n || Vt.numCols() != n)
908 Teuchos::Array<scalar_type> work(1);
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);
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);
926 template <
typename ordinal_type,
typename scalar_type>
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)
941 while (rank < m && s[rank]/s[0] > rank_threshold) rank++;
945 U.reshape(U.numRows(),rank);
946 Vt.reshape(rank, Vt.numCols());
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.
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
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)
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
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
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.