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) 54 void DGEQP3_F77(
const int*,
const int*,
double*,
const int*,
int*,
55 double*,
double*,
const int*,
int*);
60 #ifdef HAVE_STOKHOS_MATLABLIB 72 template <
typename ordinal_type,
typename scalar_type>
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)
86 template <
typename ordinal_type,
typename scalar_type>
89 Teuchos::SerialDenseVector<ordinal_type, scalar_type>&
x,
91 const Teuchos::SerialDenseVector<ordinal_type, scalar_type>&
y)
95 x[i] = alpha*
x[i] + beta*
y[i];
101 template <
typename ordinal_type,
typename scalar_type>
104 const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A)
110 if (i < A.numRows()-1)
111 os <<
";" << std::endl <<
" ";
113 os <<
"];" << std::endl;
116 #ifdef HAVE_STOKHOS_MATLABLIB 118 template <
typename ordinal_type>
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)
130 MATFile *file = matOpen(file_name.c_str(), mode);
131 TEUCHOS_ASSERT(file != NULL);
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());
141 ordinal_type ret = matPutVariable(file, matrix_name.c_str(), mat);
142 TEUCHOS_ASSERT(ret == 0);
145 ret = matClose(file);
146 TEUCHOS_ASSERT(ret == 0);
157 template <
typename ordinal_type,
typename scalar_type>
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();
170 template <
typename ordinal_type,
typename scalar_type>
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)
178 using Teuchos::getCol;
179 typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type>
SDV;
180 typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>
SDM;
183 SDM& Anc =
const_cast<SDM&
>(A);
187 if (Q.numRows() != m || Q.numCols() != k)
189 if (R.numRows() != k || R.numCols() != k)
193 SDV Aj = getCol(Teuchos::View, Anc,
j);
194 SDV Qj = getCol(Teuchos::View, Q,
j);
197 SDV Qi = getCol(Teuchos::View, Q, i);
202 Qj.scale(1.0/R(
j,
j));
212 template <
typename ordinal_type,
typename scalar_type>
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)
220 using Teuchos::getCol;
221 typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type>
SDV;
222 typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>
SDM;
225 SDM& Anc =
const_cast<SDM&
>(A);
229 if (Q.numRows() != m || Q.numCols() != k)
231 if (R.numRows() != k || R.numCols() != k)
235 SDV Ai = getCol(Teuchos::View, Anc, i);
236 SDV Qi = getCol(Teuchos::View, Q, i);
240 SDV Qi = getCol(Teuchos::View, Q, i);
242 SDV Qj = getCol(Teuchos::View, Q,
j);
248 Qi.scale(1.0/R(i,i));
257 template <
typename ordinal_type,
typename scalar_type>
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)
265 using Teuchos::getCol;
266 typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type>
SDV;
267 typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>
SDM;
270 SDM& Anc =
const_cast<SDM&
>(A);
274 if (Q.numRows() != m || Q.numCols() != k)
276 if (R.numRows() != k || R.numCols() != k)
280 SDV Ai = getCol(Teuchos::View, Anc, i);
281 SDV Qi = getCol(Teuchos::View, Q, i);
285 SDV Qi = getCol(Teuchos::View, Q, i);
287 SDV Qj = getCol(Teuchos::View, Q,
j);
293 SDV Qj = getCol(Teuchos::View, Q,
j);
299 Qi.scale(1.0/R(i,i));
310 template <
typename ordinal_type,
typename scalar_type>
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)
328 TEUCHOS_TEST_FOR_EXCEPTION(
329 w[i] != 1.0, std::logic_error,
330 "CPQR_Householder_threshold() requires unit weight vector!");
333 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
334 Teuchos::Copy, A, m, n);
338 Teuchos::Array<scalar_type> tau(k);
339 Teuchos::Array<scalar_type> work(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);
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);
352 if (R.numRows() != k || R.numCols() != n)
360 if (Q.numRows() != m || Q.numCols() != k)
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);
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)
391 template <
typename ordinal_type,
typename scalar_type>
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)
405 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
406 Teuchos::Copy, A, m, n);
407 if (Q.numRows() != m || Q.numCols() != k)
415 Teuchos::Array<scalar_type> tau(k);
416 Teuchos::Array<scalar_type> work(1);
419 DGEQP3_F77(&m, &n, AA.values(), &lda, &piv[0], &tau[0], &work[0], &lwork,
421 TEUCHOS_TEST_FOR_EXCEPTION(
422 info < 0, std::logic_error,
"dgeqp3 returned info = " << info);
427 DGEQP3_F77(&m, &n, AA.values(), &lda, &piv[0], &tau[0], &work[0], &lwork,
429 TEUCHOS_TEST_FOR_EXCEPTION(
430 info < 0, std::logic_error,
"dgeqp3 returned info = " << info);
433 if (R.numRows() != k || R.numCols() != n)
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);
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)
480 template <
typename ordinal_type,
typename scalar_type>
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)
492 TEUCHOS_TEST_FOR_EXCEPTION(
493 w[i] != 1.0, std::logic_error,
494 "CPQR_Householder_threshold() requires unit weight vector!");
504 for (rank=1; rank<m; rank++) {
509 if (r_min / r_max < rank_threshold)
514 R.reshape(rank,rank);
515 Q.reshape(Q.numRows(), rank);
532 template <
typename ordinal_type,
typename scalar_type>
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)
542 using Teuchos::getCol;
543 typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type>
SDV;
549 if (Q.numRows() != m || Q.numCols() != n)
551 if (R.numRows() != m || R.numCols() != n)
560 SDV Qj = getCol(Teuchos::View, Q,
j);
563 SDV Qtmp(m), Rtmp(m);
565 Teuchos::Array<ordinal_type> piv_orig(piv);
572 if (piv_orig[
j] != 0) {
575 nrm(
j) = nrm(nfixed);
578 SDV Qpvt = getCol(Teuchos::View, Q,
j);
579 SDV Qj = getCol(Teuchos::View, Q, nfixed);
585 piv[
j] = piv[nfixed];
596 while (j < k && sigma >= rank_threshold) {
598 SDV Qj = getCol(Teuchos::View, Q,
j);
604 if (nrm(l) > nrm(pvt))
609 SDV Qpvt = getCol(Teuchos::View, Q, pvt);
614 SDV Rpvt = getCol(Teuchos::View, R, pvt);
615 SDV Rj = getCol(Teuchos::View, R,
j);
627 Qj.scale(1.0/R(
j,
j));
629 SDV Ql = getCol(Teuchos::View, Q, l);
645 if (sigma < rank_threshold)
648 R.reshape(rank, rank);
665 template <
typename ordinal_type,
typename scalar_type>
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)
679 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> A2(Q), R2;
680 QR_MGS(rank, A2, w, Q, R2);
686 template <
typename ordinal_type,
typename scalar_type>
688 cond_R(
const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R)
710 template <
typename ordinal_type,
typename scalar_type>
713 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
714 const Teuchos::Array<scalar_type>& w)
718 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> Qt(m,n);
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);
726 return err.normInf();
733 template <
typename ordinal_type,
typename scalar_type>
736 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q)
739 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(n,n);
740 err.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q, Q, 0.0);
743 return err.normInf();
752 template <
typename ordinal_type,
typename scalar_type>
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)
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);
765 err.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
766 TEUCHOS_ASSERT(ret == 0);
768 return err.normInf();
778 template <
typename ordinal_type,
typename scalar_type>
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)
788 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> AP(m, k);
791 AP(i,
j) = A(i,piv[
j]);
792 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> err(m,k);
794 err.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
795 TEUCHOS_ASSERT(ret == 0);
798 return err.normInf();
805 template <
typename ordinal_type,
typename scalar_type>
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)
819 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
820 Teuchos::Copy, A, m, n);
823 if (U.numRows() != m || U.numCols() != m)
825 if (Vt.numRows() != n || Vt.numCols() != n)
833 Teuchos::Array<scalar_type> work(1);
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);
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);
851 template <
typename ordinal_type,
typename scalar_type>
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)
866 while (rank < m && s[rank]/s[0] > rank_threshold) rank++;
870 U.reshape(U.numRows(),rank);
871 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 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.