Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_SDMUtils.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef STOKHOS_SDM_UTILS_HPP
43#define STOKHOS_SDM_UTILS_HPP
44
45#include "Teuchos_SerialDenseMatrix.hpp"
46#include "Teuchos_SerialDenseVector.hpp"
47#include "Teuchos_SerialDenseHelpers.hpp"
48#include "Teuchos_Array.hpp"
49#include "Teuchos_LAPACK.hpp"
50#include <ostream>
51
52#define DGEQP3_F77 F77_BLAS_MANGLE(dgeqp3,DGEQP3)
53
54extern "C" {
55#if defined(HAVE_STOKHOS_MKL)
56 #include "mkl_version.h"
57 #if __INTEL_MKL__ >= 2021
58 #define MKL_NO_EXCEPT noexcept
59 #else
60 #define MKL_NO_EXCEPT
61 #endif
62#else
63 #define MKL_NO_EXCEPT
64#endif
65
66void DGEQP3_F77(const int*, const int*, double*, const int*, int*,
67 double*, double*, const int*, int*) MKL_NO_EXCEPT;
68}
69
70#include "Stokhos_ConfigDefs.h"
71
72#ifdef HAVE_STOKHOS_MATLABLIB
73extern "C" {
74#include "mat.h"
75#include "matrix.h"
76}
77#endif
78
79namespace Stokhos {
80
81 namespace detail {
82
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)
89 {
90 ordinal_type n = x.length();
91 scalar_type t = 0;
92 for (ordinal_type i=0; i<n; i++)
93 t += x[i]*w[i]*y[i];
94 return t;
95 }
96
98 template <typename ordinal_type, typename scalar_type>
99 void saxpy(
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)
104 {
105 ordinal_type n = x.length();
106 for (ordinal_type i=0; i<n; i++)
107 x[i] = alpha*x[i] + beta*y[i];
108 }
109
110 }
111
112 // Print a matrix in matlab-esque format
113 template <typename ordinal_type, typename scalar_type>
114 void
115 print_matlab(std::ostream& os,
116 const Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>& A)
117 {
118 os << "[ ";
119 for (ordinal_type i=0; i<A.numRows(); i++) {
120 for (ordinal_type j=0; j<A.numCols(); j++)
121 os << A(i,j) << " ";
122 if (i < A.numRows()-1)
123 os << ";" << std::endl << " ";
124 }
125 os << "];" << std::endl;
126 }
127
128#ifdef HAVE_STOKHOS_MATLABLIB
129 // Save a matrix to matlab MAT file
130 template <typename ordinal_type>
131 void
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)
135 {
136 // Open matfile
137 const char *mode;
138 if (overwrite)
139 mode = "w";
140 else
141 mode = "u";
142 MATFile *file = matOpen(file_name.c_str(), mode);
143 TEUCHOS_ASSERT(file != NULL);
144
145 // Convert SDM to Matlab matrix
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());
151
152 // Save matrix to file
153 ordinal_type ret = matPutVariable(file, matrix_name.c_str(), mat);
154 TEUCHOS_ASSERT(ret == 0);
155
156 // Close file
157 ret = matClose(file);
158 TEUCHOS_ASSERT(ret == 0);
159
160 // Free matlab matrix
161 mxSetPr(mat, NULL);
162 mxSetM(mat, 0);
163 mxSetN(mat, 0);
164 mxDestroyArray(mat);
165 }
166#endif
167
169 template <typename ordinal_type, typename scalar_type>
170 scalar_type vec_norm_inf(
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();
175 }
176
178
182 template <typename ordinal_type, typename scalar_type>
183 void QR_CGS(
184 ordinal_type k,
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)
189 {
190 using Teuchos::getCol;
191 typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type> SDV;
192 typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> SDM;
193
194 // getCol() requires non-const SDM
195 SDM& Anc = const_cast<SDM&>(A);
196
197 // Make sure Q is the right size
198 ordinal_type m = A.numRows();
199 if (Q.numRows() != m || Q.numCols() != k)
200 Q.shape(m,k);
201 if (R.numRows() != k || R.numCols() != k)
202 R.shape(k,k);
203
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);
207 Qj.assign(Aj);
208 for (ordinal_type i=0; i<j; i++) {
209 SDV Qi = getCol(Teuchos::View, Q, i);
210 R(i,j) = detail::weighted_inner_product(Qi, Aj, w);
211 detail::saxpy(1.0, Qj, -R(i,j), Qi); // Q(:,j) = 1.0*Q(:,j) - R(i,j)*Q(:,i)
212 }
213 R(j,j) = std::sqrt(detail::weighted_inner_product(Qj, Qj, w));
214 Qj.scale(1.0/R(j,j));
215 }
216
217 }
218
220
224 template <typename ordinal_type, typename scalar_type>
225 void QR_MGS(
226 ordinal_type k,
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)
231 {
232 using Teuchos::getCol;
233 typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type> SDV;
234 typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> SDM;
235
236 // getCol() requires non-const SDM
237 SDM& Anc = const_cast<SDM&>(A);
238
239 // Make sure Q is the right size
240 ordinal_type m = A.numRows();
241 if (Q.numRows() != m || Q.numCols() != k)
242 Q.shape(m,k);
243 if (R.numRows() != k || R.numCols() != k)
244 R.shape(k,k);
245
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);
249 Qi.assign(Ai);
250 }
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);
255 scalar_type v = detail::weighted_inner_product(Qi, Qj, w);
256 detail::saxpy(1.0, Qi, -v, Qj); // Q(:,i) = 1.0*Q(:,i) - v*Q(:,j)
257 R(j,i) += v;
258 }
259 R(i,i) = std::sqrt(detail::weighted_inner_product(Qi, Qi, w));
260 Qi.scale(1.0/R(i,i));
261 }
262 }
263
265
269 template <typename ordinal_type, typename scalar_type>
271 ordinal_type k,
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)
276 {
277 using Teuchos::getCol;
278 typedef Teuchos::SerialDenseVector<ordinal_type,scalar_type> SDV;
279 typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> SDM;
280
281 // getCol() requires non-const SDM
282 SDM& Anc = const_cast<SDM&>(A);
283
284 // Make sure Q is the right size
285 ordinal_type m = A.numRows();
286 if (Q.numRows() != m || Q.numCols() != k)
287 Q.shape(m,k);
288 if (R.numRows() != k || R.numCols() != k)
289 R.shape(k,k);
290
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);
294 Qi.assign(Ai);
295 }
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);
300 scalar_type v = detail::weighted_inner_product(Qi, Qj, w);
301 detail::saxpy(1.0, Qi, -v, Qj); // Q(:,i) = 1.0*Q(:,i) - v*Q(:,j)
302 R(j,i) += v;
303 }
304 for (ordinal_type j=i-1; j>=0; j--) {
305 SDV Qj = getCol(Teuchos::View, Q, j);
306 scalar_type v = detail::weighted_inner_product(Qi, Qj, w);
307 detail::saxpy(1.0, Qi, -v, Qj); // Q(:,i) = 1.0*Q(:,i) - v*Q(:,j)
308 R(j,i) += v;
309 }
310 R(i,i) = std::sqrt(detail::weighted_inner_product(Qi, Qi, w));
311 Qi.scale(1.0/R(i,i));
312 }
313 }
314
316
322 template <typename ordinal_type, typename scalar_type>
323 void
325 ordinal_type k,
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)
330 {
332 ordinal_type m = A.numRows();
333 ordinal_type n = A.numCols();
334 ordinal_type kk = std::min(m,n);
335 if (k > kk)
336 k = kk;
337
338 // Check that each component of w is 1
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!");
343
344 // Lapack routine overwrites A, so copy into temporary matrix
345 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
346 Teuchos::Copy, A, m, n);
347
348 // QR
349 ordinal_type lda = AA.stride();
350 Teuchos::Array<scalar_type> tau(k);
351 Teuchos::Array<scalar_type> work(1);
352 ordinal_type info;
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);
357 lwork = work[0];
358 work.resize(lwork);
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);
362
363 // Extract R
364 if (R.numRows() != k || R.numCols() != n)
365 R.shape(k,n);
366 R.putScalar(0.0);
367 for (ordinal_type i=0; i<k; i++)
368 for (ordinal_type j=i; j<n; j++)
369 R(i,j) = AA(i,j);
370
371 // Extract Q
372 if (Q.numRows() != m || Q.numCols() != k)
373 Q.shape(m,k);
374 lwork = -1;
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);
378 lwork = work[0];
379 work.resize(lwork);
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)
384 Q.shape(m, k);
385 for (ordinal_type i=0; i<m; i++)
386 for (ordinal_type j=0; j<k; j++)
387 Q(i,j) = AA(i,j);
388 }
389
391
403 template <typename ordinal_type, typename scalar_type>
404 void
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)
410 {
412 ordinal_type m = A.numRows();
413 ordinal_type n = A.numCols();
414 ordinal_type k = std::min(m,n);
415
416 // Lapack routine overwrites A, so copy into temporary matrix
417 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
418 Teuchos::Copy, A, m, n);
419 if (Q.numRows() != m || Q.numCols() != k)
420 Q.shape(m,k);
421
422 // Teuchos LAPACK interface doesn't have dgeqp3, so call it directly
423
424 // Workspace query
425 ordinal_type lda = AA.stride();
426 piv.resize(n);
427 Teuchos::Array<scalar_type> tau(k);
428 Teuchos::Array<scalar_type> work(1);
429 ordinal_type lwork = -1;
430 ordinal_type info;
431 DGEQP3_F77(&m, &n, AA.values(), &lda, &piv[0], &tau[0], &work[0], &lwork,
432 &info);
433 TEUCHOS_TEST_FOR_EXCEPTION(
434 info < 0, std::logic_error, "dgeqp3 returned info = " << info);
435
436 // Column pivoted QR
437 lwork = work[0];
438 work.resize(lwork);
439 DGEQP3_F77(&m, &n, AA.values(), &lda, &piv[0], &tau[0], &work[0], &lwork,
440 &info);
441 TEUCHOS_TEST_FOR_EXCEPTION(
442 info < 0, std::logic_error, "dgeqp3 returned info = " << info);
443
444 // Extract R
445 if (R.numRows() != k || R.numCols() != n)
446 R.shape(k,n);
447 R.putScalar(0.0);
448 for (ordinal_type i=0; i<k; i++)
449 for (ordinal_type j=i; j<n; j++)
450 R(i,j) = AA(i,j);
451
452 // Extract Q
453 lwork = -1;
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);
457 lwork = work[0];
458 work.resize(lwork);
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)
463 Q.shape(m, k);
464 for (ordinal_type i=0; i<m; i++)
465 for (ordinal_type j=0; j<k; j++)
466 Q(i,j) = AA(i,j);
467
468 // Transform piv to zero-based indexing
469 for (ordinal_type i=0; i<n; i++)
470 piv[i] -= 1;
471 }
472
474
492 template <typename ordinal_type, typename scalar_type>
493 ordinal_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)
501 {
502 // Check that each component of w is 1
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!");
507
508 // Compute full QR
509 CPQR_Householder3(A, Q, R, piv);
510
511 // Find leading block of R such that cond(R) <= 1/tau
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)
522 break;
523 }
524
525 // Extract blocks from Q and R
526 R.reshape(rank,rank);
527 Q.reshape(Q.numRows(), rank);
528
529 return rank;
530 }
531
533
544 template <typename ordinal_type, typename scalar_type>
545 ordinal_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)
553 {
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);
559
560 // Make sure Q is the right size
561 if (Q.numRows() != m || Q.numCols() != n)
562 Q.shape(m,n);
563 if (R.numRows() != m || R.numCols() != n)
564 R.shape(m,n);
565 if (piv.size() != n)
566 piv.resize(n);
567 Q.assign(A);
568
569 // Compute column norms
570 SDV nrm(n);
571 for (ordinal_type j=0; j<n; j++) {
572 SDV Qj = getCol(Teuchos::View, Q, j);
573 nrm(j) = detail::weighted_inner_product(Qj, Qj, w);
574 }
575 SDV Qtmp(m), Rtmp(m);
576
577 Teuchos::Array<ordinal_type> piv_orig(piv);
578 for (ordinal_type i=0; i<n; i++)
579 piv[i] = i;
580
581 // Prepivot any columns requested by setting piv[i] != 0
582 ordinal_type nfixed = 0;
583 for (ordinal_type j=0; j<n; j++) {
584 if (piv_orig[j] != 0) {
585 if (j != nfixed) {
586 scalar_type tmp = nrm(j);
587 nrm(j) = nrm(nfixed);
588 nrm(nfixed) = tmp;
589
590 SDV Qpvt = getCol(Teuchos::View, Q, j);
591 SDV Qj = getCol(Teuchos::View, Q, nfixed);
592 Qtmp.assign(Qpvt);
593 Qpvt.assign(Qj);
594 Qj.assign(Qtmp);
595
596 ordinal_type t = piv[j];
597 piv[j] = piv[nfixed];
598 piv[nfixed] = t;
599 }
600 ++nfixed;
601 }
602 }
603
604 scalar_type sigma = 1.0 + rank_threshold;
605 scalar_type r_max = 0.0;
606 ordinal_type j=0;
607 R.putScalar(0.0);
608 while (j < k && sigma >= rank_threshold) {
609
610 SDV Qj = getCol(Teuchos::View, Q, j);
611
612 // Find pivot column if we are passed the pre-pivot stage
613 if (j >= nfixed) {
614 ordinal_type pvt = j;
615 for (ordinal_type l=j+1; l<n; l++)
616 if (nrm(l) > nrm(pvt))
617 pvt = l;
618
619 // Interchange column j and pvt
620 if (pvt != j) {
621 SDV Qpvt = getCol(Teuchos::View, Q, pvt);
622 Qtmp.assign(Qpvt);
623 Qpvt.assign(Qj);
624 Qj.assign(Qtmp);
625
626 SDV Rpvt = getCol(Teuchos::View, R, pvt);
627 SDV Rj = getCol(Teuchos::View, R, j);
628 Rtmp.assign(Rpvt);
629 Rpvt.assign(Rj);
630 Rj.assign(Rtmp);
631
632 ordinal_type t = piv[pvt];
633 piv[pvt] = piv[j];
634 piv[j] = t;
635 }
636 }
637
638 R(j,j) = std::sqrt(detail::weighted_inner_product(Qj, Qj, w));
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);
642 scalar_type t = detail::weighted_inner_product(Qj, Ql, w);
643 R(j,l) = t;
644 detail::saxpy(1.0, Ql, -t, Qj);
645
646 // Update norms
647 nrm(l) = detail::weighted_inner_product(Ql, Ql, w);
648 }
649
650 if (std::abs(R(j,j)) > r_max)
651 r_max = R(j,j);
652 sigma = std::abs(R(j,j)/r_max);
653 j++;
654 }
655
656 ordinal_type rank = j;
657 if (sigma < rank_threshold)
658 rank--;
659 Q.reshape(m, rank);
660 R.reshape(rank, rank);
661
662 return rank;
663 }
664
666
677 template <typename ordinal_type, typename scalar_type>
678 ordinal_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)
686 {
687 // First do standard column-pivoted QR
688 ordinal_type rank = CPQR_MGS_threshold(rank_threshold, A, w, Q, R, piv);
689
690 // Now apply standard MGS to Q
691 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> A2(Q), R2;
692 QR_MGS(rank, A2, w, Q, R2);
693
694 return rank;
695 }
696
698 template <typename ordinal_type, typename scalar_type>
699 scalar_type
700 cond_R(const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& R)
701 {
702 ordinal_type k = R.numRows();
703 if (R.numCols() < k)
704 k = R.numCols();
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)
709 r_max = R(i,i);
710 if (std::abs(R(i,i)) < r_min)
711 r_min = R(i,i);
712 }
713 scalar_type cond_r = r_max / r_min;
714 return cond_r;
715 }
716
718
722 template <typename ordinal_type, typename scalar_type>
723 scalar_type
725 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q,
726 const Teuchos::Array<scalar_type>& w)
727 {
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++)
737 err(i,i) -= 1.0;
738 return err.normInf();
739 }
740
742
745 template <typename ordinal_type, typename scalar_type>
746 scalar_type
748 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& Q)
749 {
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++)
754 err(i,i) -= 1.0;
755 return err.normInf();
756 }
757
759
764 template <typename ordinal_type, typename scalar_type>
765 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)
770 {
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);
776 ordinal_type ret =
777 err.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
778 TEUCHOS_ASSERT(ret == 0);
779 err -= AA;
780 return err.normInf();
781 }
782
784
790 template <typename ordinal_type, typename scalar_type>
791 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)
797 {
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);
805 ordinal_type ret =
806 err.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
807 TEUCHOS_ASSERT(ret == 0);
808 err -= AP;
809
810 return err.normInf();
811 }
812
814
817 template <typename ordinal_type, typename scalar_type>
818 void
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)
823 {
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();
829
830 // Copy A since GESVD overwrites it (always)
831 Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> AA(
832 Teuchos::Copy, A, m, n);
833
834 // Resize appropriately
835 if (U.numRows() != m || U.numCols() != m)
836 U.shape(m,m);
837 if (Vt.numRows() != n || Vt.numCols() != n)
838 Vt.shape(n,n);
839 if (s.size() != k)
840 s.resize(k);
841 ordinal_type ldu = U.stride();
842 ordinal_type ldv = Vt.stride();
843
844 // Workspace query
845 Teuchos::Array<scalar_type> work(1);
846 ordinal_type lwork = -1;
847 ordinal_type info;
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);
852
853 // Do SVD
854 lwork = work[0];
855 work.resize(lwork);
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);
860
861 }
862
863 template <typename ordinal_type, typename scalar_type>
864 ordinal_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)
871 {
872 // Compute full SVD
873 svd(A, s, U, Vt);
874
875 // Find leading block of S such that cond(S) <= 1/tau
876 ordinal_type rank = 0;
877 ordinal_type m = s.size();
878 while (rank < m && s[rank]/s[0] > rank_threshold) rank++;
879
880 // Extract blocks from s, U and Vt
881 s.resize(rank);
882 U.reshape(U.numRows(),rank);
883 Vt.reshape(rank, Vt.numCols());
884
885 return rank;
886 }
887
888}
889
890#endif
#define DGEQP3_F77
#define MKL_NO_EXCEPT
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