queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
TeuchosMatrix.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2017 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <queso/Defines.h>
26 
27 #ifdef QUESO_HAS_TRILINOS
28 #include <queso/TeuchosMatrix.h>
29 #include <queso/TeuchosVector.h>
30 #endif
31 #include <sys/time.h>
32 #include <cmath>
33 
34 #ifdef QUESO_HAS_TRILINOS
35 
36 namespace QUESO {
37 
38 using std:: cout;
39 using std:: endl;
40 
41 // ---------------------------------------------------
42 // shaped constructor (square/rectangular)------------
43 TeuchosMatrix::TeuchosMatrix( // can be a rectangular matrix
44  const BaseEnvironment& env,
45  const Map& map,
46  unsigned int nCols)
47  :
48  Matrix (env,map),
49  m_inverse (NULL),
50  m_svdColMap (NULL),
51  m_svdUmat (NULL),
52  m_svdSvec (NULL),
53  m_svdVmat (NULL),
54  m_svdVTmat (NULL),
55  m_determinant (-INFINITY),
56  m_lnDeterminant(-INFINITY),
57  v_pivoting (NULL),
58  m_signum (0),
59  m_isSingular (false)
60 {
61  m_mat.shape(map.NumGlobalElements(),nCols);
62  m_LU.shape(0,0);
63 }
64 
65 // ---------------------------------------------------
66 // shaped constructor (square) -----------------------
67 TeuchosMatrix::TeuchosMatrix( // square matrix
68  const BaseEnvironment& env,
69  const Map& map,
70  double diagValue)
71  :
72  Matrix (env,map),
73  m_inverse (NULL),
74  m_svdColMap (NULL),
75  m_svdUmat (NULL),
76  m_svdSvec (NULL),
77  m_svdVmat (NULL),
78  m_svdVTmat (NULL),
79  m_determinant (-INFINITY),
80  m_lnDeterminant(-INFINITY),
81  v_pivoting (NULL),
82  m_signum (0),
83  m_isSingular (false)
84 {
85  m_mat.shape (map.NumGlobalElements(),map.NumGlobalElements());
86  m_LU.shape(0,0);
87 
88  for (unsigned int i = 0; i < (unsigned int) m_mat.numRows(); ++i) {
89  m_mat(i,i) = diagValue;
90  }
91 }
92 // ---------------------------------------------------
93 // shaped constructor (diagonal) ---------------------
94 // Kemelli tested on 12/05/12
95 TeuchosMatrix::TeuchosMatrix( // square matrix
96  const TeuchosVector& v,
97  double diagValue)
98  :
99  Matrix (v.env(),v.map()),
100  m_inverse (NULL),
101  m_svdColMap (NULL),
102  m_svdUmat (NULL),
103  m_svdSvec (NULL),
104  m_svdVmat (NULL),
105  m_svdVTmat (NULL),
106  m_determinant (-INFINITY),
107  m_lnDeterminant(-INFINITY),
108  v_pivoting (NULL),
109  m_signum (0),
110  m_isSingular (false)
111 {
112  m_mat.shape (v.sizeLocal(),v.sizeLocal());
113  m_LU.shape(0,0);
114 
115  for (unsigned int i = 0; i < (unsigned int) m_mat.numRows(); ++i) {
116  m_mat(i,i) = diagValue;
117  }
118 }
119 
120 // ---------------------------------------------------
121 // shaped constructor (diagonal, from vector) --------
122 // Kemelli, 12/05/12 - tested
123 TeuchosMatrix::TeuchosMatrix(const TeuchosVector& v) // square matrix
124  :
125  Matrix (v.env(),v.map()),
126  m_inverse (NULL),
127  m_svdColMap (NULL),
128  m_svdUmat (NULL),
129  m_svdSvec (NULL),
130  m_svdVmat (NULL),
131  m_svdVTmat (NULL),
132  m_determinant (-INFINITY),
133  m_lnDeterminant(-INFINITY),
134  v_pivoting (NULL),
135  m_signum (0),
136  m_isSingular (false)
137 {
138  m_mat.shape (v.sizeLocal(),v.sizeLocal());
139  m_LU.shape(0,0);
140 
141  unsigned int dim = v.sizeLocal();
142 
143  for (unsigned int i = 0; i < dim; ++i) {
144  m_mat(i,i) = v[i];
145  }
146 }
147 // ---------------------------------------------------
148 // copy constructor-----------------------------------
149 // Kemelli 12/05/12 - tested
150 TeuchosMatrix::TeuchosMatrix(const TeuchosMatrix& B) // can be a rectangular matrix
151  :
152  Matrix (B.env(),B.map()),
153  m_inverse (NULL),
154  m_svdColMap (NULL),
155  m_svdUmat (NULL),
156  m_svdSvec (NULL),
157  m_svdVmat (NULL),
158  m_svdVTmat (NULL),
159  m_determinant (-INFINITY),
160  m_lnDeterminant(-INFINITY),
161  v_pivoting (NULL),
162  m_signum (0),
163  m_isSingular (false)
164 {
165  m_mat.shape (B.numRowsLocal(),B.numCols());
166  m_LU.shape(0,0);
167 
168  this->Matrix::base_copy(B);
169  this->copy(B);
170 }
171 
172 // ---------------------------------------------------
173 // destructor ----------------------------------------
174 TeuchosMatrix::~TeuchosMatrix()
175 {
176  this->resetLU();
177 }
178 
179 
180 // ---------------------------------------------------
181 // Set methodos --------------------------------------
182 
183 TeuchosMatrix& TeuchosMatrix::operator=(const TeuchosMatrix& obj)
184 {
185  this->resetLU();
186  this->copy(obj);
187  return *this;
188 }
189 
190 // ---------------------------------------------------
191 // Kemelli 12/05/12 - tested
192 TeuchosMatrix& TeuchosMatrix::operator*=(double a)
193 {
194  this->resetLU();
195  int iRC;
196  iRC = m_mat.scale(a);
197  queso_require_msg(!(iRC), "scaling failed");
198  return *this;
199 }
200 
201 // ---------------------------------------------------
202 // Kemelli 12/05/12 - tested
203 TeuchosMatrix& TeuchosMatrix::operator/=(double a)
204 {
205  this->resetLU();
206  m_mat.scale(1./a);
207  return *this;
208 }
209 
210 // ---------------------------------------------------
211 // Kemelli 12/05/12 - tested
212 TeuchosMatrix& TeuchosMatrix::operator+=(const TeuchosMatrix& rhs)
213 {
214  this->resetLU();
215 
216  unsigned int i,j, nrows=rhs.numRowsLocal(), ncols=rhs.numCols();
217 
218  for(i=0; i< nrows ; i++)
219  for (j = 0; j < ncols; j++)
220  (*this)(i,j) += rhs(i,j);
221 
222  return *this;
223 }
224 // ---------------------------------------------------
225 // Kemelli 12/05/12 - tested
226 TeuchosMatrix&
227 TeuchosMatrix::operator-=(const TeuchosMatrix& rhs)
228 {
229  this->resetLU();
230 
231  unsigned int i,j, nrows=rhs.numRowsLocal(), ncols=rhs.numCols();
232 
233  for(i=0; i< nrows ; i++)
234  for (j = 0; j < ncols; j++)
235  (*this)(i,j) -= rhs(i,j);
236 
237  return *this;
238 }
239 
240 // ---------------------------------------------------
241 // Accessor methods ----------------------------------
242 
243 // Kemelli 12/05/12 - tested
244 double& TeuchosMatrix::operator()(unsigned int i, unsigned int j)
245 {
246  this->resetLU();
247  if ((i >= (unsigned int) m_mat.numRows()) ||
248  (j >= (unsigned int) m_mat.numCols())) {
249  std::cerr << "In TeuchosMatrix::operator(i,j) (non-const)"
250  << ": i = " << i
251  << ", j = " << j
252  << ", m_mat.numRows() = " << (unsigned int) m_mat.numRows()
253  << ", m_mat.numCols() = " << (unsigned int) m_mat.numCols()
254  << std::endl;
255  queso_require_less_msg(i, (unsigned int) m_mat.numRows(), "i is too large");
256  queso_require_less_msg(j, (unsigned int) m_mat.numCols(), "j is too large");
257  }
258  return m_mat(i,j);
259 }
260 
261 // ---------------------------------------------------
262 const double& TeuchosMatrix::operator()(unsigned int i, unsigned int j) const
263 {
264  queso_require_less_msg(i, (unsigned int) m_mat.numRows(), "i is too large");
265  queso_require_less_msg(j, (unsigned int) m_mat.numCols(), "j is too large");
266  return m_mat(i,j);
267 }
268 
269 // ---------------------------------------------------
270 // Attribute methods ---------------------------------
271 
272 // Kemelli 12/05/12 - tested
273 unsigned int
274 TeuchosMatrix::numRowsLocal() const
275 {
276  return (unsigned int) m_mat.numRows();
277 }
278 
279 // ---------------------------------------------------
280 // Kemelli 12/05/12 - tested
281 unsigned int
282 TeuchosMatrix::numRowsGlobal() const
283 {
284  return (unsigned int) m_mat.numRows();
285 }
286 
287 // ---------------------------------------------------
288 // Kemelli 12/05/12 - tested
289 unsigned int
290 TeuchosMatrix::numCols() const
291 {
292  return (unsigned int) m_mat.numCols();
293 };
294 
295 // ---------------------------------------------------
296 // Kemelli 12/05/12 - tested
297 double*
298 TeuchosMatrix::values()
299 {
300  return m_mat.values();
301 };
302 
303 // ---------------------------------------------------
304 // Kemelli 12/05/12 - tested
305 int
306 TeuchosMatrix::stride()
307 {
308  return m_mat.stride();
309 };
310 
311 
312 
313 // ---------------------------------------------------
314 double
315 TeuchosMatrix::max() const
316 {
317  double value = -INFINITY;
318 
319  unsigned int nRows = this->numRowsLocal();
320  unsigned int nCols = this->numCols();
321  double aux = 0.;
322  for (unsigned int i = 0; i < nRows; i++) {
323  for (unsigned int j = 0; j < nCols; j++) {
324  aux = (*this)(i,j);
325  if (aux > value) value = aux;
326  }
327  }
328 
329  return value;
330 }
331 
332 
333 // ---------------------------------------------------
334 unsigned int
335 TeuchosMatrix::rank(double absoluteZeroThreshold, double relativeZeroThreshold) const
336 {
337  int iRC = 0;
338  iRC = internalSvd();
339  if (iRC) {}; // just to remove compiler warning
340 
341  TeuchosVector relativeVec(*m_svdSvec);
342  if (relativeVec[0] > 0.) {
343  relativeVec = (1./relativeVec[0])*relativeVec;
344  }
345 
346  unsigned int rankValue = 0;
347  for (unsigned int i = 0; i < relativeVec.sizeLocal(); ++i) {
348  if (( (*m_svdSvec)[i] >= absoluteZeroThreshold ) &&
349  ( relativeVec [i] >= relativeZeroThreshold )) {
350  rankValue += 1;
351  }
352  }
353 
354  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
355  *m_env.subDisplayFile() << "In TeuchosMatrix::rank()"
356  << ": this->numRowsLocal() = " << this->numRowsLocal()
357  << ", this->numCols() = " << this->numCols()
358  << ", absoluteZeroThreshold = " << absoluteZeroThreshold
359  << ", relativeZeroThreshold = " << relativeZeroThreshold
360  << ", rankValue = " << rankValue
361  << ", m_svdSvec = " << *m_svdSvec
362  << ", relativeVec = " << relativeVec
363  << std::endl;
364  }
365 
366  return rankValue;
367 }
368 
369 // ---------------------------------------------------
370 // checked 1/07/13
371 TeuchosMatrix
372 TeuchosMatrix::transpose() const
373 {
374  unsigned int nRows = this->numRowsLocal();
375  unsigned int nCols = this->numCols();
376 
377  queso_require_equal_to_msg(nRows, nCols, "routine works only for square matrices");
378 
379  TeuchosMatrix mat(m_env,m_map,nCols);
380  for (unsigned int row = 0; row < nRows; ++row) {
381  for (unsigned int col = 0; col < nCols; ++col) {
382  mat(row,col) = (*this)(col,row);
383  }
384  }
385 
386  return mat;
387 }
388 
389 // ---------------------------------------------------
390 // tested 1/31/13
391 TeuchosMatrix
392 TeuchosMatrix::inverse() const
393 {
394  unsigned int nRows = this->numRowsLocal();
395  unsigned int nCols = this->numCols();
396 
397  queso_require_equal_to_msg(nRows, nCols, "matrix is not square");
398 
399  if (m_inverse == NULL) {
400  m_inverse = new TeuchosMatrix(m_env,m_map,nCols);
401  TeuchosVector unitVector(m_env,m_map);
402  unitVector.cwSet(0.);
403  TeuchosVector multVector(m_env,m_map);
404  for (unsigned int j = 0; j < nCols; ++j) {
405  if (j > 0) unitVector[j-1] = 0.;
406  unitVector[j] = 1.;
407  this->invertMultiply(unitVector, multVector);
408  for (unsigned int i = 0; i < nRows; ++i) {
409  (*m_inverse)(i,j) = multVector[i];
410  }
411  }
412  }
413  if (m_env.checkingLevel() >= 1) {
414  *m_env.subDisplayFile() << "CHECKING In TeuchosMatrix::inverse()"
415  << ": M.lnDet = " << this->lnDeterminant()
416  << ", M^{-1}.lnDet = " << m_inverse->lnDeterminant()
417  << std::endl;
418  }
419  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
420  *m_env.subDisplayFile() << "In TeuchosMatrix::inverse():"
421  << "\n M = " << *this
422  << "\n M^{-1} = " << *m_inverse
423  << "\n M*M^{-1} = " << (*this)*(*m_inverse)
424  << "\n M^{-1}*M = " << (*m_inverse)*(*this)
425  << std::endl;
426  }
427 
428 /* cout << "Inverse ---------------------\n";
429  for (int i=0; i<m_inverse->numRowsLocal(); i++) {
430  for (int j=0; j<m_inverse->numCols(); j++)
431  cout<< m_inverse->m_mat(i,j) << " ";
432  cout << endl;
433  }
434  cout << endl;
435 */
436  return *m_inverse;
437 }
438 
439 // ----------------------------------------------
440 //checked 12/10/12
441 double
442 TeuchosMatrix::determinant() const
443 {
444  if (m_determinant == -INFINITY)
445  {
446  if(m_LU.numRows() ==0 && m_LU.numCols() ==0) //dummy
447  {
448  TeuchosVector tmpB(m_env,m_map);
449  TeuchosVector tmpX(m_env,m_map);
450  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
451  *m_env.subDisplayFile() << "In TeuchosMatrix::lnDeterminant()"
452  << ": before 'this->invertMultiply()'"
453  << std::endl;
454  }
455  this->invertMultiply(tmpB,tmpX);
456  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
457  *m_env.subDisplayFile() << "In TeuchosMatrix::lnDeterminant()"
458  << ": after 'this->invertMultiply()'"
459  << std::endl;
460  }
461  }
462  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
463  *m_env.subDisplayFile() << "In TeuchosMatrix::lnDeterminant()"
464  << ": before computing det"
465  << std::endl;
466  }
467 
468  double det = 1.0;
469  double lnDet = 0.0;
470  for (int i=0;i<m_LU.numCols();i++) {
471  det *= m_LU(i,i);
472  lnDet += std::log(m_LU(i,i));
473  }
474 
475  m_determinant = det;
476  m_lnDeterminant = lnDet;
477 
478  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
479  *m_env.subDisplayFile() << "In TeuchosMatrix::lnDeterminant()"
480  << ": after computing det"
481  << std::endl;
482  }
483  }
484 
485  return m_determinant;
486 }
487 
488 // ----------------------------------------------
489 //checked 12/10/12
490 double
491 TeuchosMatrix::lnDeterminant() const
492 {
493  if (m_lnDeterminant == -INFINITY)
494  {
495  if(m_LU.numRows() ==0 && m_LU.numCols() ==0) //dummy
496  {
497  TeuchosVector tmpB(m_env,m_map);
498  TeuchosVector tmpX(m_env,m_map);
499  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
500  *m_env.subDisplayFile() << "In TeuchosMatrix::lnDeterminant()"
501  << ": before 'this->invertMultiply()'"
502  << std::endl;
503  }
504  this->invertMultiply(tmpB,tmpX);
505  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
506  *m_env.subDisplayFile() << "In TeuchosMatrix::lnDeterminant()"
507  << ": after 'this->invertMultiply()'"
508  << std::endl;
509  }
510  }
511  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
512  *m_env.subDisplayFile() << "In TeuchosMatrix::lnDeterminant()"
513  << ": before computing lnDet"
514  << std::endl;
515  }
516 
517  double det = 1.0;
518  double lnDet = 0.0;
519  for (int i=0;i<m_LU.numCols();i++) {
520  det *= m_LU(i,i);
521  lnDet += std::log(m_LU(i,i));
522  }
523 
524  m_determinant = det;
525  m_lnDeterminant = lnDet;
526 
527  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
528  *m_env.subDisplayFile() << "In TeuchosMatrix::lnDeterminant()"
529  << ": after computing lnDet"
530  << std::endl;
531  }
532  }
533 
534  return m_lnDeterminant;
535 }
536 
537 // ---------------------------------------------------
538 // Norm methods --------------------------------------
539 
540 double
541 TeuchosMatrix::normFrob() const
542 {
543  return m_mat.normFrobenius ();
544 }
545 
546 // ---------------------------------------------------
547 double
548 TeuchosMatrix::normMax() const
549 {
550  double value = 0.;
551 
552  unsigned int nRows = this->numRowsLocal();
553  unsigned int nCols = this->numCols();
554  double aux = 0.;
555  for (unsigned int i = 0; i < nRows; i++) {
556  for (unsigned int j = 0; j < nCols; j++) {
557  aux = fabs((*this)(i,j));
558  if (aux > value) value = aux;
559  }
560  }
561  return value;
562 }
563 
564 // Mathematical methods ------------------------------
565 // ---------------------------------------------------
566 
567 // Kemelli: added and tested 12/10/12
568 int TeuchosMatrix::chol()
569 {
570  int return_success =0 ;
571 /* If UPLO = 'L', the leading N-by-N lower triangular part of A contains the lower
572  * triangular part of the matrix A (lets call it L), and the strictly upper
573  * triangular part of A is not referenced. Thus, the upper triangular part
574  * of A +++must be manually+++ overwritten with L^T. */
575 
576  Teuchos::LAPACK<int, double> lapack;
577  int info;//= 0: successful exit
578  //< 0: if INFO = -i, the i-th argument had an illegal value
579  //> 0: if INFO = i, the leading minor of order i is not
580  // positive definite, and the factorization could not be
581  // completed.
582  char UPLO = 'U';
583  //= 'U': Upper triangle of A is stored;
584  //= 'L': Lower triangle of A is stored.
585 
586  lapack.POTRF (UPLO, m_mat.numRows(), m_mat.values(), m_mat.stride(), &info);
587 
588  // Overwriting the upper triangular part of the input matrix A with L^T
589  //(the diagonal terms are identical for both L and L^T)
590 
591  for (int i=0;i<m_mat.numRows();i++){
592  for (int j=i+1;j<m_mat.numCols();j++)
593  m_mat(i,j) = m_mat(j,i) ;
594  }
595 
596  if (info != 0) {
597  std::cerr << "In TeuchosMtrix::chol()"
598  << ": INFO = " << info
599  << ",\n INFO < 0: if INFO = -i, the i-th argument had an illegal value."
600  << "\n INFO > 0: if INFO = i, the leading minor of order i is not "
601  << " positive definite, and the factorization could not be completed."
602  << std::endl;
603  return_success =1 ;
604  }
605 
606  UQ_RC_MACRO(info, // Yes, *not* a fatal check on RC
607  m_env.worldRank(),
608  "TeuchosMatrix::chol()",
609  "matrix is not positive definite",
610  UQ_MATRIX_IS_NOT_POS_DEFINITE_RC);
611 
612  return return_success;
613 };
614 
615 // ---------------------------------------------------
616 int
617 TeuchosMatrix::svd(TeuchosMatrix& matU, TeuchosVector& vecS, TeuchosMatrix& matVt) const
618 {
619  unsigned int nRows = this->numRowsLocal();
620  unsigned int nCols = this->numCols();
621 
622  queso_require_msg(!((matU.numRowsLocal() != nRows) || (matU.numCols() != nCols)), "invalid matU");
623 
624  queso_require_equal_to_msg(vecS.sizeLocal(), nCols, "invalid vecS");
625 
626  queso_require_msg(!((matVt.numRowsLocal() != nCols) || (matVt.numCols() != nCols)), "invalid matVt");
627 
628  int iRC = internalSvd();
629 
630  matU = *m_svdUmat;
631  vecS = *m_svdSvec;
632  matVt = *m_svdVTmat;
633 
634  return iRC;
635 }
636 
637 //---------------------------------------------------------------
638 //tested 2/1/13
639 const TeuchosMatrix& TeuchosMatrix::svdMatU() const
640 {
641  int iRC = 0;
642  iRC = internalSvd();
643  if (iRC) {}; // just to remove compiler warning
644  return *m_svdUmat;
645 }
646 
647 //---------------------------------------------------------------
648 //tested 2/1/13
649 const TeuchosMatrix& TeuchosMatrix::svdMatV() const
650 {
651  int iRC = 0;
652  iRC = internalSvd();
653  if (iRC) {}; // just to remove compiler warning
654  return *m_svdVmat;
655 }
656 
657 //---------------------------------------------------------------
658 // checked 2/27/13
659 /* An orthogonal matrix M has a norm-preserving property, i.e.
660  * for any vector v, \f[||Mv|| = ||v|| \f] (1). Then:
661  * \f[ min(||Ax − b||^2) = min(||Ax − b||) = min(||UDVT x − b||) = (1) min(||DV x − U b||) \f].
662  * Substituting \f[ y = VT x \f] and \f[ b' = UT b \f] gives us \f[ Dy = b' \f]
663  * with D a diagonal matrix. Or, \f[ y = inv(D)*UT*b \f] and we only have to
664  * solve the linear system: \f[ VT x = y \f].
665  */
666 int
667 TeuchosMatrix::svdSolve(const TeuchosVector& rhsVec, TeuchosVector& solVec) const
668 {
669  unsigned int nRows = this->numRowsLocal();
670  unsigned int nCols = this->numCols();
671  unsigned int i;
672 
673  queso_require_equal_to_msg(rhsVec.sizeLocal(), nRows, "invalid rhsVec");
674 
675  queso_require_equal_to_msg(solVec.sizeLocal(), nCols, "invalid solVec");
676 
677  int iRC = internalSvd();
678 
679  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
680  *m_env.subDisplayFile() << "In TeuchosMatrix::svdSolve():"
681  << "\n this->numRowsLocal() = " << this->numRowsLocal()
682  << ", this->numCols() = " << this->numCols()
683  << "\n m_svdUmat->numRowsLocal() = " << m_svdUmat->numRowsLocal()
684  << ", m_svdUmat->numCols() = " << m_svdUmat->numCols()
685  << "\n m_svdVmat->numRowsLocal() = " << m_svdVmat->numRowsLocal()
686  << ", m_svdVmat->numCols() = " << m_svdVmat->numCols()
687  << "\n m_svdSvec->sizeLocal() = " << m_svdSvec->sizeLocal()
688  << "\n rhsVec.sizeLocal() = " << rhsVec.sizeLocal()
689  << "\n solVec.sizeLocal() = " << solVec.sizeLocal()
690  << std::endl;
691  }
692 
693  if (iRC == 0)
694  {
695  TeuchosMatrix invD = TeuchosMatrix(solVec);
696  TeuchosMatrix auxMatrix = TeuchosMatrix(solVec);
697 
698  for (i=0; i<nRows; i++){
699  invD(i,i) = 1./(m_svdSvec->values()[i]);
700  }
701 
702  // GESV: For a system Ax=b, on entry, b contains the right-hand side b;
703  // on exit it contains the solution x.
704  // Thus, intead of doing y = inv(D)*UT*rhsVec = auxMatrix*rhsVec, with
705  // auxMatrix = inv(D)*UT; and then assigning solVec=y (requirement for using GESV)
706  // lets do solVec= auxMatrix*rhsVec, and save one step in the calculation
707 
708  auxMatrix = invD * svdMatU().transpose();
709  solVec = auxMatrix*rhsVec;
710 
711  // solve the linear system VT * solVec = y
712  // GESV changes the values of the matrix, so lets make a copy of it and use it.
713 
714  TeuchosMatrix* aux_m_svdVTmat = new TeuchosMatrix(*m_svdVTmat);
715 
716  int ipiv[0], info;
717  Teuchos::LAPACK<int, double> lapack;
718  lapack.GESV(nCols, 1, aux_m_svdVTmat->values(), aux_m_svdVTmat->stride(), ipiv, &solVec[0], solVec.sizeLocal(), &info );
719 
720  /* GESV output INFO: = 0: successful exit
721  * < 0: if INFO = -i, the i-th argument had an illegal value
722  * > 0: if INFO = i, U(i,i) is exactly zero. The factorization
723  * has been completed, but the factor U is exactly
724  * singular, so the solution could not be computed. */
725 
726  iRC = info;
727  delete aux_m_svdVTmat;
728  }
729  return iRC;
730 }
731 
732 // ---------------------------------------------------
733 //checked 2/27/13
734 int
735 TeuchosMatrix::svdSolve(const TeuchosMatrix& rhsMat, TeuchosMatrix& solMat) const
736 {
737  unsigned int nRows = this->numRowsLocal();
738  unsigned int nCols = this->numCols();
739 
740  queso_require_equal_to_msg(rhsMat.numRowsLocal(), nRows, "invalid rhsMat");
741 
742  queso_require_equal_to_msg(solMat.numRowsLocal(), nCols, "invalid solMat");
743 
744  queso_require_equal_to_msg(rhsMat.numCols(), solMat.numCols(), "rhsMat and solMat are not compatible");
745 
746  TeuchosVector rhsVec(m_env,rhsMat.map());
747  TeuchosVector solVec(m_env,solMat.map());
748  int iRC = 0;
749  for (unsigned int j = 0; j < rhsMat.numCols(); ++j) {
750  rhsVec = rhsMat.getColumn(j);
751  iRC = this->svdSolve(rhsVec, solVec);
752  if (iRC) break;
753  solMat.setColumn(j,solVec);
754  }
755 
756  return iRC;
757 }
758 
759 // ---------------------------------------------------
760 // implemented/checked 1/8/13
761 // multiply this matrix by vector x and store in vector y, which is returned
762 // eg: TeuchosVector v1(paramSpace.zeroVector() );
763 // v1=covMatrix.multiply(meanVector);
764 TeuchosVector
765 TeuchosMatrix::multiply(const TeuchosVector& x) const
766 {
767  queso_require_equal_to_msg(this->numCols(), x.sizeLocal(), "matrix and vector have incompatible sizes");
768 
769  TeuchosVector y(m_env,m_map);
770  this->multiply(x,y);
771 
772  return y;
773 }
774 
775 // ---------------------------------------------------
776 //Kemelli checked 12/06/12
777 TeuchosVector
778 TeuchosMatrix::invertMultiply(const TeuchosVector& b) const
779 {
780  queso_require_equal_to_msg(this->numCols(), b.sizeLocal(), "matrix and rhs have incompatible sizes");
781  TeuchosVector x(m_env,m_map);
782 
783  this->invertMultiply(b,x);
784 
785  return x;
786 }
787 
788 // ---------------------------------------------------
789 //Kemelli checked 12/06/12
790 void
791 TeuchosMatrix::invertMultiply(const TeuchosVector& b, TeuchosVector& x) const
792 {
793  queso_require_equal_to_msg(this->numCols(), b.sizeLocal(), "matrix and rhs have incompatible sizes");
794 
795  queso_require_equal_to_msg(x.sizeLocal(), b.sizeLocal(), "solution and rhs have incompatible sizes");
796 
797  if (m_LU.numCols() == 0 && m_LU.numRows() == 0)
798  {
799  queso_require_msg(!(v_pivoting), "v_pivoting should be NULL");
800 
801  //allocate m_LU and v_pivoting
802  m_LU = m_mat;
803  v_pivoting =(int *) malloc(sizeof(int)*m_LU.numCols() );
804 
805  queso_require_not_equal_to_msg(m_LU.numCols(), 0 && m_LU.numRows() == 0, "malloc() failed");
806 
807  queso_require_msg(v_pivoting, "malloc() failed");
808 
809  if (m_inDebugMode) {
810  std::cout << "In TeuchosMatrix::invertMultiply()"
811  << ": before LU decomposition, m_LU = ";
812  for (int i=0;i<3;i++){
813  for (int j=0;j<3;j++)
814  cout << m_LU(i,j) <<"\t" ;
815  cout << endl;
816  }
817  }
818 
819  // Perform an LU factorization of matrix m_LU. Checked 12/06/12
820  Teuchos::LAPACK<int, double> lapack;
821  int info;
822 
823  lapack.GETRF( m_LU.numRows(), m_LU.numCols(), m_LU.values(), m_LU.stride(), v_pivoting, &info );
824 
825  if (info != 0) {
826  std::cerr << "In TeuchosMatrix::invertMultiply()"
827  << ", after lapack.GETRF"
828  << ": INFO = " << info
829  << ",\nINFO < 0: if INFO = -i, the i-th argument had an illegal value.\n"
830  << "INFO > 0: if INFO = i, U(i,i) is exactly zero. The factorization \n"
831  << "has been completed, but the factor U is exactly singular, and division \n"
832  << "by zero will occur if it is used to solve a system of equations."
833  << std::endl;
834  }
835  queso_require_msg(!(info), "GETRF() failed");
836 
837  if (info > 0)
838  m_isSingular = true;
839 
840  if (m_inDebugMode) {
841  std::cout << "In TeuchosMatrix::invertMultiply()"
842  << ": after LU decomposition, m_LU = ";
843  for (int i=0;i<3;i++){
844  for (int j=0;j<3;j++)
845  cout << m_LU(i,j) <<"\t" ;
846  cout << endl;
847  }
848  std::cout << std::endl;
849  }
850  }
851  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
852  *m_env.subDisplayFile() << "In TeuchosMatrix::invertMultiply()"
853  << ": before 'lapack.GETRS()'"
854  << std::endl;
855  }
856 
857  // Solve the linear system.
858  Teuchos::LAPACK<int, double> lapack;
859  int NRHS = 1; // NRHS: number of right hand sides, i.e., the number
860  // of columns of the matrix B. In this case, vector b.
861  char TRANS = 'N'; // 'N': A * x= B (No transpose). Specifies the
862  // form of the system of equations.
863  int info02;
864 
865  //GETRS expects the matrix to be already factored in LU and uses the
866  //same ipiv vector, which are the pivot indices of the LU factorization
867 
868  x=b;
869  lapack.GETRS(TRANS, m_LU.numRows(), NRHS, m_LU.values(), m_LU.stride(), v_pivoting, &x[0],x.sizeLocal(), &info02 );
870 
871  if (info02 != 0) {
872  std::cerr << "In TeuchosMatrix::invertMultiply()"
873  << ", after lapack.GETRS - solve LU system"
874  << ": INFO = " << info02
875  << ",\nINFO < 0: if INFO = -i, the i-th argument had an illegal value.\n"
876  << std::endl;
877  }
878  queso_require_msg(!(info02), "GETRS() failed");
879  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
880  *m_env.subDisplayFile() << "In TeuchosMatrix::invertMultiply()"
881  << ", after lapack.GETRS() - solve LU system."
882  << std::endl;
883  }
884  if (m_inDebugMode) {
885  TeuchosVector tmpVec(b - (*this)*x);
886  std::cout << "In TeuchosMatrix::invertMultiply()"
887  << ": ||b - Ax||_2 = " << tmpVec.norm2()
888  << ": ||b - Ax||_2/||b||_2 = " << tmpVec.norm2()/b.norm2()
889  << std::endl;
890  }
891  return;
892 }
893 
894 // ----------------------------------------------
895 TeuchosMatrix
896 TeuchosMatrix::invertMultiply(const TeuchosMatrix& B) const
897 {
898  TeuchosMatrix X(m_env,m_map,B.numCols());
899  this->invertMultiply(B,X);
900  return X;
901 }
902 
903 // ----------------------------------------------
904 //checked 12/10/12
905 void
906 TeuchosMatrix::invertMultiply(const TeuchosMatrix& B, TeuchosMatrix& X) const
907 {
908  // Sanity Checks
909  queso_require_msg(!((B.numRowsLocal() != X.numRowsLocal()) || (B.numCols() != X.numCols())), "Matrices B and X are incompatible");
910 
911  queso_require_equal_to_msg(this->numRowsLocal(), X.numRowsLocal(), "This and X matrices are incompatible");
912 
913  // Some local variables used within the loop.
914  TeuchosVector b(m_env, m_map);
915  TeuchosVector x(m_env, m_map);
916 
917  for( unsigned int j = 0; j < B.numCols(); ++j ) {
918  b = B.getColumn(j);
919  //invertMultiply will only do the LU factorization once and store it.
920  x = this->invertMultiply(b);
921  X.setColumn(j,x);
922  }
923  return;
924 }
925 
926 //-----------------------------------------------
927 TeuchosVector
928 TeuchosMatrix::invertMultiplyForceLU(const TeuchosVector& b) const
929 {
930  queso_require_equal_to_msg(this->numCols(), b.sizeLocal(), "matrix and rhs have incompatible sizes");
931 
932  TeuchosVector x(m_env,m_map);
933  this->invertMultiplyForceLU(b,x);
934 
935  return x;
936 }
937 
938 // ---------------------------------------------------
939 // implemented and checked on 1/9/13
940 void
941 TeuchosMatrix::invertMultiplyForceLU(const TeuchosVector& b, TeuchosVector& x) const
942 {
943  queso_require_equal_to_msg(this->numCols(), b.sizeLocal(), "matrix and rhs have incompatible sizes");
944 
945  queso_require_equal_to_msg(x.sizeLocal(), b.sizeLocal(), "solution and rhs have incompatible sizes");
946 
947  if (m_LU.numCols() == 0 && m_LU.numRows() == 0) {
948  queso_require_msg(!(v_pivoting), "v_pivoting should be NULL");
949  }
950 
951  //allocate m_LU , yes outside the if above
952  m_LU = m_mat;
953 
954  queso_require_not_equal_to_msg(m_LU.numCols(), 0 && m_LU.numRows() == 0, "Teuchos atttribuition m_LU = m_mat failed");
955 
956  //allocate v_pivoting
957  if ( v_pivoting == NULL )
958  v_pivoting =(int *) malloc(sizeof(int)*m_LU.numCols() );
959 
960  queso_require_msg(v_pivoting, "malloc() for v_pivoting failed");
961 
962  // Perform an LU factorization of matrix m_LU. Checked 12/06/12
963  Teuchos::LAPACK<int, double> lapack;
964  int info;
965 
966  lapack.GETRF( m_LU.numRows(), m_LU.numCols(), m_LU.values(), m_LU.stride(), v_pivoting, &info );
967 
968  if (info != 0) {
969  std::cerr << "In TeuchosMatrix::invertMultiply()"
970  << ", after lapack.GETRF"
971  << ": INFO = " << info
972  << ",\nINFO < 0: if INFO = -i, the i-th argument had an illegal value.\n"
973  << "INFO > 0: if INFO = i, U(i,i) is exactly zero. The factorization \n"
974  << "has been completed, but the factor U is exactly singular, and division \n"
975  << "by zero will occur if it is used to solve a system of equations."
976  << std::endl;
977  }
978 
979  queso_require_msg(!(info), "GETRF() failed");
980 
981  if (info > 0)
982  m_isSingular = true;
983 
984  // Solve the linear system.
985  int NRHS = 1; // NRHS: number of right hand sides, i.e., the number
986  // of columns of the matrix B. In this case, vector b.
987  char TRANS = 'N'; // 'N': A * x= B (No transpose). Specifies the
988  // form of the system of equations.
989  int info02;
990 
991  //GETRS expects the matrix to be already factored in LU and uses the
992  //same ipiv vector, which are the pivot indices of the LU factorization
993 
994  x=b;
995  lapack.GETRS(TRANS, m_LU.numRows(), NRHS, m_LU.values(), m_LU.stride(), v_pivoting, &x[0],x.sizeLocal(), &info02 );
996 
997  if (info02 != 0) {
998  std::cerr << "In TeuchosMatrix::invertMultiply()"
999  << ", after lapack.GETRS - solve LU system"
1000  << ": INFO = " << info02
1001  << ",\nINFO < 0: if INFO = -i, the i-th argument had an illegal value.\n"
1002  << std::endl;
1003  }
1004 
1005  queso_require_msg(!(info02), "GETRS() failed");
1006 
1007  return;
1008 }
1009 // ---------------------------------------------------
1010 // Implemented and checked on 1/7/2013
1011 
1012 /* DSYEV computes all eigenvalues and, optionally, eigenvectors of a real symmetric matrix A.
1013 *
1014 * Arguments
1015 * JOBZ = 'N': Compute eigenvalues only;
1016 * = 'V': Compute eigenvalues and eigenvectors.
1017 * UPLO = 'U': Upper triangle of A is stored;
1018 * = 'L': Lower triangle of A is stored.
1019 * N The order of the matrix A. N >= 0.
1020 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N).On entry, the symmetric
1021 * matrix A. If UPLO = 'U', the leading N-by-N upper triangular part of A contains the
1022 * upper triangular part of the matrix A. If UPLO = 'L', the leading N-by-N lower triangular
1023 * part of A contains the lower triangular part of the matrix A. On exit, if JOBZ = 'V', then
1024 * if INFO = 0, A contains the orthonormal eigenvectors of the matrix A. If JOBZ = 'N', then
1025 * on exit the lower triangle (if UPLO='L') or the upper triangle (if UPLO='U') of A,
1026 * including the diagonal, is destroyed.
1027 * LDA (input) INTEGER - the leading dimension of the array A. LDA >= max(1,N).
1028 * W (output) DOUBLE PRECISION array, dimension (N)
1029 * If INFO = 0, the eigenvalues in ascending order.
1030 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
1031 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
1032 * LWORK (input) INTEGER - the length of the array WORK. LWORK >= max(1,3*N-1).
1033 * INFO (output) INTEGER
1034 * = 0: successful exit
1035 * < 0: if INFO = -i, the i-th argument had an illegal value
1036 * > 0: if INFO = i, the algorithm failed to converge; i off-diagonal elements of
1037 * an intermediate tridiagonal form did not converge to zero.*/
1038 void
1039 TeuchosMatrix::eigen(TeuchosVector& eigenValues, TeuchosMatrix* eigenVectors) const
1040 {
1041  unsigned int n = eigenValues.sizeLocal();
1042 
1043  queso_require_not_equal_to_msg(n, 0, "invalid input vector size");
1044 
1045  if (eigenVectors) {
1046  queso_require_equal_to_msg(eigenValues.sizeLocal(), eigenVectors->numRowsLocal(), "different input vector sizes");
1047  }
1048 
1049  // At the end of execution, Lapack funcion destroys the input matrix
1050  // So, lets not use m_mat, but instead, a copy of it
1051  Teuchos::SerialDenseMatrix<int,double> copy_m_mat;
1052  copy_m_mat = m_mat;
1053 
1054  Teuchos::LAPACK<int, double> lapack;
1055  int lwork = 3*n -1;
1056  int info;
1057  char UPLO = 'L';
1058  double *W, *WORK;
1059 
1060  W = new double[n];
1061  WORK = new double[lwork];
1062 
1063  if (eigenVectors == NULL) {
1064  char JOBZ = 'N';//eigenvalues only
1065 
1066  lapack.SYEV(JOBZ, UPLO, copy_m_mat.numRows(), copy_m_mat.values(), copy_m_mat.stride(), &W[0], &WORK[0], lwork, &info);
1067 
1068  for (unsigned int i=0; i< n; i++)
1069  eigenValues[i] = W[i];
1070 
1071  queso_require_equal_to_msg(info, 0, "invalid input vector size");
1072  }
1073  else {
1074  char JOBZ = 'V';//eigenvalues and eigenvectors
1075 
1076  lapack.SYEV(JOBZ, UPLO, copy_m_mat.numRows(), copy_m_mat.values(), copy_m_mat.stride(), &W[0], &WORK[0], lwork, &info);
1077 
1078  for (unsigned int i=0; i< n; i++)
1079  eigenValues[i] = W[i];
1080 
1081  eigenVectors->m_mat = copy_m_mat;
1082  }
1083 
1084  if (info != 0) {
1085  std::cerr << "In TeuchosMtrix::eigen()"
1086  << ": INFO = " << info
1087  << ",\n INFO < 0: if INFO = -i, the i-th argument had an illegal value."
1088  << "\n INFO > 0: if INFO = i, the algorithm failed to converge; i off-diagonal "
1089  << " elements of an intermediate tridiagonal form did not converge to zero."
1090  << std::endl;
1091  }
1092 
1093  return;
1094 }
1095 // ---------------------------------------------------
1096 void
1097 TeuchosMatrix::largestEigen(double& eigenValue, TeuchosVector& eigenVector) const
1098 {
1099  unsigned int n = eigenVector.sizeLocal();
1100 
1101  queso_require_not_equal_to_msg(n, 0, "invalid input vector size");
1102  Teuchos::LAPACK<int, double> lapack;
1103 
1104  //SYEV destroys the input matrix, so lets operate in a copy
1105  Teuchos::SerialDenseMatrix<int,double> copy_m_mat;
1106  copy_m_mat = m_mat;
1107 
1108  int lwork = 3*n -1;
1109  int info;
1110  char UPLO = 'L';
1111  char JOBZ = 'V'; //eigenvalues and eigenvectors
1112  double *W, *WORK;
1113 
1114  W = new double[n];
1115  WORK = new double[lwork];
1116 
1117  lapack.SYEV(JOBZ, UPLO, copy_m_mat.numRows(), copy_m_mat.values(), copy_m_mat.stride(), &W[0], &WORK[0], lwork, &info);
1118 
1119  if (info != 0) {
1120  std::cerr << "In TeuchosMtrix::largestEigen()"
1121  << ": INFO = " << info
1122  << ",\n INFO < 0: if INFO = -i, the i-th argument had an illegal value."
1123  << "\n INFO > 0: if INFO = i, the algorithm failed to converge; i off-diagonal "
1124  << " elements of an intermediate tridiagonal form did not converge to zero."
1125  << std::endl;
1126  }
1127  // If INFO = 0, W contains the eigenvalues in ascending order.
1128  // Thus the largest eigenvalue is in W[n-1].
1129  eigenValue = W[n-1];
1130 
1131  // Eigenvector associated to the largest eigenvalue.
1132  // Stored in the n-th column of matrix copy_m_mat.
1133  for (int i=0; i< copy_m_mat.numRows(); i++)
1134  eigenVector[i] = copy_m_mat(i,n-1);
1135 
1136  return;
1137 }
1138 
1139 // ---------------------------------------------------
1140 void
1141 TeuchosMatrix::smallestEigen(double& eigenValue, TeuchosVector& eigenVector) const
1142 {
1143  unsigned int n = eigenVector.sizeLocal();
1144 
1145  queso_require_not_equal_to_msg(n, 0, "invalid input vector size");
1146  Teuchos::LAPACK<int, double> lapack;
1147 
1148  //SYEV destroys the input matrix, so lets operate in a copy
1149  Teuchos::SerialDenseMatrix<int,double> copy_m_mat;
1150  copy_m_mat = m_mat;
1151 
1152  int lwork = 3*n -1;
1153  int info;
1154  char UPLO = 'L';
1155  char JOBZ = 'V';//eigenvalues and eigenvectors
1156  double *W, *WORK;
1157 
1158  W = new double[n];
1159  WORK = new double[lwork];
1160 
1161  lapack.SYEV(JOBZ, UPLO, copy_m_mat.numRows(), copy_m_mat.values(), copy_m_mat.stride(), &W[0], &WORK[0], lwork, &info);
1162 
1163  if (info != 0) {
1164  std::cerr << "In TeuchosMtrix::smallestEigen()"
1165  << ": INFO = " << info
1166  << ",\n INFO < 0: if INFO = -i, the i-th argument had an illegal value."
1167  << "\n INFO > 0: if INFO = i, the algorithm failed to converge; i off-diagonal "
1168  << " elements of an intermediate tridiagonal form did not converge to zero."
1169  << std::endl;
1170  }
1171 
1172  // If INFO = 0, W contains the eigenvalues in ascending order.
1173  // Thus the smallest eigenvalue is in W[0].
1174  eigenValue = W[0];
1175 
1176  // Eigenvector associated to the smallest eigenvalue.
1177  // Stored in the n-th column of matrix copy_m_mat.
1178  for (int i=0; i< copy_m_mat.numRows(); i++)
1179  eigenVector[i] = copy_m_mat(i,0);
1180 
1181  return;
1182 }
1183 
1184 // Get/Set methodos ----------------------------------
1185 // ---------------------------------------------------
1186 
1187 // ---------------------------------------------------
1188 void
1189 TeuchosMatrix::cwSet(double value)
1190 {
1191  m_mat.putScalar(value);
1192  return;
1193 }
1194 
1195 // ---------------------------------------------------
1196 // tested on 1/31/13: copies matrix mat to this matrix, starting at row rowId and column colId
1197 void
1198 TeuchosMatrix::cwSet(unsigned int rowId,unsigned int colId,const TeuchosMatrix& mat)
1199 {
1200  queso_require_less_msg(rowId, this->numRowsLocal(), "invalid rowId");
1201 
1202  queso_require_less_equal_msg((rowId + mat.numRowsLocal()), this->numRowsLocal(), "invalid vec.numRowsLocal()");
1203 
1204  queso_require_less_msg(colId, this->numCols(), "invalid colId");
1205 
1206  queso_require_less_equal_msg((colId + mat.numCols()), this->numCols(), "invalid vec.numCols()");
1207 
1208  for (unsigned int i = 0; i < mat.numRowsLocal(); ++i) {
1209  for (unsigned int j = 0; j < mat.numCols(); ++j) {
1210  (*this)(rowId+i,colId+j) = mat(i,j);
1211  }
1212  }
1213 
1214  return;
1215 }
1216 
1217 // ---------------------------------------------------
1218 //Kemelli tested 07/12/12
1219 void
1220 TeuchosMatrix::getColumn(unsigned int column_num, TeuchosVector& column) const
1221 {
1222  // Sanity checks
1223  queso_require_less_msg(column_num, this->numCols(), "Specified row number not within range");
1224 
1225  queso_require_equal_to_msg(column.sizeLocal(), this->numRowsLocal(), "column vector not same size as this matrix");
1226 
1227  // Temporary working pointer
1228  const double* temp_ptr ;
1229 
1230  // get the column_num- of matrix m_mat
1231  temp_ptr = m_mat[column_num];
1232 
1233  // Copy column from Teuchos matrix into our TeuchosVector object
1234  for (unsigned int i=0; i< column.sizeLocal();i++)
1235  column[i] = temp_ptr[i];
1236 
1237  return;
1238 }
1239 
1240 // ---------------------------------------------------
1241 TeuchosVector
1242 TeuchosMatrix::getColumn(unsigned int column_num ) const
1243 {
1244  // We rely on the void getColumn's to do sanity checks
1245  // So we don't do them here.
1246  TeuchosVector column(m_env, m_map);
1247  this->getColumn( column_num, column );
1248  return column;
1249 }
1250 
1251 
1252 // ---------------------------------------------------
1253 //Kemelli tested 07/12/12
1254 void
1255 TeuchosMatrix::setColumn(unsigned int column_num, const TeuchosVector& column)
1256 {
1257  this->resetLU();
1258  // Sanity checks
1259  queso_require_less_msg(column_num, this->numCols(), "Specified column number not within range");
1260 
1261  queso_require_equal_to_msg(column.sizeLocal(), this->numRowsLocal(), "column vector not same size as this matrix");
1262 
1263  for (unsigned int i =0; i < column.sizeLocal(); i++)
1264  m_mat(i,column_num) = column[i];
1265 
1266  return;
1267 }
1268 
1269 // ---------------------------------------------------
1270 // tested 1/31/13
1271 void
1272 TeuchosMatrix::getRow(unsigned int row_num, TeuchosVector& row) const
1273 {
1274  // Sanity checks
1275  queso_require_less_msg(row_num, this->numRowsLocal(), "Specified row number not within range");
1276 
1277  queso_require_equal_to_msg(row.sizeLocal(), this->numRowsLocal(), "row vector not same size as this matrix");
1278 
1279  // Copy row from Teuchos matrix into our TeuchosVector object
1280  for (unsigned int i=0; i< row.sizeLocal();i++)
1281  row[i] = m_mat(row_num,i);
1282 
1283  return;
1284 }
1285 
1286 // ---------------------------------------------------
1287 // tested 1/31/13
1288 TeuchosVector
1289 TeuchosMatrix::getRow(unsigned int row_num ) const
1290 {
1291  // We rely on the void getRow's to do sanity checks
1292  // So we don't do them here.
1293  TeuchosVector row(m_env, m_map);
1294  this->getRow( row_num, row );
1295  return row;
1296 }
1297 
1298 // ---------------------------------------------------
1299 // tested 1/31/13
1300 void
1301 TeuchosMatrix::setRow (const unsigned int row_num, const TeuchosVector& row)
1302 {
1303  // Sanity checks
1304  queso_require_less_msg(row_num, this->numRowsLocal(), "Specified row number not within range");
1305 
1306  queso_require_equal_to_msg(row.sizeLocal(), this->numRowsLocal(), "row vector not same size as this matrix");
1307 
1308  // Copy our TeuchosVector object to our Teuchos Matrix
1309  for (unsigned int i=0; i< row.sizeLocal();i++)
1310  m_mat(row_num,i) = row[i] ;
1311 
1312  return;
1313 }
1314 
1315 // ---------------------------------------------------
1316 // Kemelli 12/05/12 - tested
1317 void
1318 TeuchosMatrix::zeroLower(bool includeDiagonal)
1319 {
1320  unsigned int nRows = this->numRowsLocal();
1321  unsigned int nCols = this->numCols();
1322 
1323  queso_require_equal_to_msg(nRows, nCols, "routine works only for square matrices");
1324  this->resetLU();
1325 
1326  if (includeDiagonal) {
1327  for (unsigned int i = 0; i < nRows; i++) {
1328  for (unsigned int j = 0; j <= i; j++) {
1329  (*this)(i,j) = 0.;
1330  }
1331  }
1332  }
1333  else {
1334  for (unsigned int i = 0; i < nRows; i++) {
1335  for (unsigned int j = 0; j < i; j++) {
1336  (*this)(i,j) = 0.;
1337  }
1338  }
1339  }
1340 
1341  return;
1342 }
1343 
1344 // ---------------------------------------------------
1345 // Kemelli 12/05/12 - tested
1346 void
1347 TeuchosMatrix::zeroUpper(bool includeDiagonal)
1348 {
1349  unsigned int nRows = this->numRowsLocal();
1350  unsigned int nCols = this->numCols();
1351 
1352  queso_require_equal_to_msg(nRows, nCols, "routine works only for square matrices");
1353  this->resetLU();
1354 
1355  if (includeDiagonal) {
1356  for (unsigned int i = 0; i < nRows; i++) {
1357  for (unsigned int j = i; j < nCols; j++) {
1358  (*this)(i,j) = 0.;
1359  }
1360  }
1361  }
1362  else {
1363  for (unsigned int i = 0; i < nRows; i++) {
1364  for (unsigned int j = (i+1); j < nCols; j++) {
1365  (*this)(i,j) = 0.;
1366  }
1367  }
1368  }
1369 
1370  return;
1371 }
1372 
1373 // ---------------------------------------------------
1374 void
1375 TeuchosMatrix::filterSmallValues(double thresholdValue)
1376 {
1377  unsigned int nRows = this->numRowsLocal();
1378  unsigned int nCols = this->numCols();
1379 
1380  for (unsigned int i = 0; i < nRows; ++i) {
1381  for (unsigned int j = 0; j < nCols; ++j) {
1382  double aux = (*this)(i,j);
1383  // If 'thresholdValue' is negative, no values will be filtered
1384  if ((aux < 0. ) && (-thresholdValue < aux)) {
1385  (*this)(i,j) = 0.;
1386  }
1387  if ((aux > 0. ) && (thresholdValue > aux)) {
1388  (*this)(i,j) = 0.;
1389  }
1390  }
1391  }
1392  return;
1393 }
1394 
1395 // ---------------------------------------------------
1396 void
1397 TeuchosMatrix::filterLargeValues(double thresholdValue)
1398 {
1399  unsigned int nRows = this->numRowsLocal();
1400  unsigned int nCols = this->numCols();
1401 
1402  for (unsigned int i = 0; i < nRows; ++i) {
1403  for (unsigned int j = 0; j < nCols; ++j) {
1404  double aux = (*this)(i,j);
1405  // If 'thresholdValue' is negative, no values will be filtered
1406  if ( (aux < 0. ) && (-thresholdValue > aux))
1407  (*this)(i,j) = 0.;
1408 
1409  if ((aux > 0. ) && (thresholdValue < aux))
1410  (*this)(i,j) = 0.;
1411  }
1412  }
1413  return;
1414 }
1415 
1416 
1417 // ----------------------------------------------
1418 // tested 1/31/13
1419 void
1420 TeuchosMatrix::fillWithTranspose(const TeuchosMatrix& mat)
1421 {
1422  unsigned int nRows = mat.numRowsLocal();
1423  unsigned int nCols = mat.numCols();
1424  queso_require_equal_to_msg(this->numRowsLocal(), nCols, "inconsistent local number of rows");
1425  queso_require_equal_to_msg(this->numCols(), nRows, "inconsistent number of cols");
1426 
1427  for (unsigned int row = 0; row < nRows; ++row) {
1428  for (unsigned int col = 0; col < nCols; ++col) {
1429  (*this)(col,row) = mat(row,col);
1430  }
1431  }
1432 
1433  return;
1434 }
1435 
1436 // ---------------------------------------------------
1437 // added 2/28/13
1438 void
1439 TeuchosMatrix::fillWithBlocksDiagonally(const std::vector<const TeuchosMatrix* >& matrices)
1440 {
1441  unsigned int sumNumRowsLocals = 0;
1442  unsigned int sumNumCols = 0;
1443  for (unsigned int i = 0; i < matrices.size(); ++i) {
1444  sumNumRowsLocals += matrices[i]->numRowsLocal();
1445  sumNumCols += matrices[i]->numCols();
1446  }
1447  queso_require_equal_to_msg(this->numRowsLocal(), sumNumRowsLocals, "inconsistent local number of rows");
1448  queso_require_equal_to_msg(this->numCols(), sumNumCols, "inconsistent number of cols");
1449 
1450  unsigned int cumulativeRowId = 0;
1451  unsigned int cumulativeColId = 0;
1452  for (unsigned int i = 0; i < matrices.size(); ++i) {
1453  unsigned int nRows = matrices[i]->numRowsLocal();
1454  unsigned int nCols = matrices[i]->numCols();
1455  for (unsigned int rowId = 0; rowId < nRows; ++rowId) {
1456  for (unsigned int colId = 0; colId < nCols; ++colId) {
1457  (*this)(cumulativeRowId + rowId, cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
1458  }
1459  }
1460  cumulativeRowId += nRows;
1461  cumulativeColId += nCols;
1462  }
1463 
1464  return;
1465 }
1466 // ---------------------------------------------------
1467 // added 2/28/13
1468 void
1469 TeuchosMatrix::fillWithBlocksDiagonally(const std::vector<TeuchosMatrix* >& matrices)
1470 {
1471  unsigned int sumNumRowsLocals = 0;
1472  unsigned int sumNumCols = 0;
1473  for (unsigned int i = 0; i < matrices.size(); ++i) {
1474  sumNumRowsLocals += matrices[i]->numRowsLocal();
1475  sumNumCols += matrices[i]->numCols();
1476  }
1477  queso_require_equal_to_msg(this->numRowsLocal(), sumNumRowsLocals, "inconsistent local number of rows");
1478  queso_require_equal_to_msg(this->numCols(), sumNumCols, "inconsistent number of cols");
1479 
1480  unsigned int cumulativeRowId = 0;
1481  unsigned int cumulativeColId = 0;
1482  for (unsigned int i = 0; i < matrices.size(); ++i) {
1483  unsigned int nRows = matrices[i]->numRowsLocal();
1484  unsigned int nCols = matrices[i]->numCols();
1485  for (unsigned int rowId = 0; rowId < nRows; ++rowId) {
1486  for (unsigned int colId = 0; colId < nCols; ++colId) {
1487  (*this)(cumulativeRowId + rowId, cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
1488  }
1489  }
1490  cumulativeRowId += nRows;
1491  cumulativeColId += nCols;
1492  }
1493 
1494  return;
1495 }
1496 // ---------------------------------------------------
1497 // added 2/28/13
1498 void
1499 TeuchosMatrix::fillWithBlocksHorizontally(const std::vector<const TeuchosMatrix* >& matrices)
1500 {
1501  unsigned int sumNumCols = 0;
1502  for (unsigned int i = 0; i < matrices.size(); ++i) {
1503  queso_require_equal_to_msg(this->numRowsLocal(), matrices[i]->numRowsLocal(), "inconsistent local number of rows");
1504  sumNumCols += matrices[i]->numCols();
1505  }
1506  queso_require_equal_to_msg(this->numCols(), sumNumCols, "inconsistent number of cols");
1507 
1508  unsigned int cumulativeColId = 0;
1509  for (unsigned int i = 0; i < matrices.size(); ++i) {
1510  unsigned int nRows = matrices[i]->numRowsLocal();
1511  unsigned int nCols = matrices[i]->numCols();
1512  for (unsigned int rowId = 0; rowId < nRows; ++rowId) {
1513  for (unsigned int colId = 0; colId < nCols; ++colId) {
1514  (*this)(rowId, cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
1515  }
1516  }
1517  cumulativeColId += nCols;
1518  }
1519 
1520  return;
1521 }
1522 
1523 // ---------------------------------------------------
1524 // added 2/28/13
1525 void
1526 TeuchosMatrix::fillWithBlocksHorizontally(const std::vector<TeuchosMatrix* >& matrices)
1527 {
1528  unsigned int sumNumCols = 0;
1529  for (unsigned int i = 0; i < matrices.size(); ++i) {
1530  queso_require_equal_to_msg(this->numRowsLocal(), matrices[i]->numRowsLocal(), "inconsistent local number of rows");
1531  sumNumCols += matrices[i]->numCols();
1532  }
1533  queso_require_equal_to_msg(this->numCols(), sumNumCols, "inconsistent number of cols");
1534 
1535  unsigned int cumulativeColId = 0;
1536  for (unsigned int i = 0; i < matrices.size(); ++i) {
1537  unsigned int nRows = matrices[i]->numRowsLocal();
1538  unsigned int nCols = matrices[i]->numCols();
1539  for (unsigned int rowId = 0; rowId < nRows; ++rowId) {
1540  for (unsigned int colId = 0; colId < nCols; ++colId) {
1541  (*this)(rowId, cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
1542  }
1543  }
1544  cumulativeColId += nCols;
1545  }
1546 
1547  return;
1548 }
1549 
1550 // ---------------------------------------------------
1551 // added 2/28/13
1552 void
1553 TeuchosMatrix::fillWithBlocksVertically(const std::vector<const TeuchosMatrix* >& matrices)
1554 {
1555  unsigned int sumNumRows = 0;
1556  for (unsigned int i = 0; i < matrices.size(); ++i) {
1557  queso_require_equal_to_msg(this->numCols(), matrices[i]->numCols(), "inconsistent local number of cols");
1558  sumNumRows += matrices[i]->numRowsLocal();
1559  }
1560  queso_require_equal_to_msg(this->numRowsLocal(), sumNumRows, "inconsistent number of rows");
1561 
1562  unsigned int cumulativeRowId = 0;
1563  for (unsigned int i = 0; i < matrices.size(); ++i) {
1564  unsigned int nRows = matrices[i]->numRowsLocal();
1565  unsigned int nCols = matrices[i]->numCols();
1566  for (unsigned int rowId = 0; rowId < nRows; ++rowId) {
1567  for (unsigned int colId = 0; colId < nCols; ++colId) {
1568  (*this)(cumulativeRowId + rowId, colId) = (*(matrices[i]))(rowId,colId);
1569  }
1570  }
1571  cumulativeRowId += nRows;
1572  }
1573  return;
1574 }
1575 
1576 // ---------------------------------------------------
1577 // added 2/28/13
1578 void
1579 TeuchosMatrix::fillWithBlocksVertically(const std::vector<TeuchosMatrix* >& matrices)
1580 {
1581  unsigned int sumNumRows = 0;
1582  for (unsigned int i = 0; i < matrices.size(); ++i) {
1583  queso_require_equal_to_msg(this->numCols(), matrices[i]->numCols(), "inconsistent local number of cols");
1584  sumNumRows += matrices[i]->numRowsLocal();
1585  }
1586  queso_require_equal_to_msg(this->numRowsLocal(), sumNumRows, "inconsistent number of rows");
1587 
1588  unsigned int cumulativeRowId = 0;
1589  for (unsigned int i = 0; i < matrices.size(); ++i) {
1590  unsigned int nRows = matrices[i]->numRowsLocal();
1591  unsigned int nCols = matrices[i]->numCols();
1592  for (unsigned int rowId = 0; rowId < nRows; ++rowId) {
1593  for (unsigned int colId = 0; colId < nCols; ++colId) {
1594  (*this)(cumulativeRowId + rowId, colId) = (*(matrices[i]))(rowId,colId);
1595  }
1596  }
1597  cumulativeRowId += nRows;
1598  }
1599  return;
1600 }
1601 
1602 // ---------------------------------------------------
1603 // added 2/28/13
1604 void
1605 TeuchosMatrix::fillWithTensorProduct(const TeuchosMatrix& mat1, const TeuchosMatrix& mat2)
1606 {
1607  queso_require_equal_to_msg(this->numRowsLocal(), (mat1.numRowsLocal() * mat2.numRowsLocal()), "inconsistent local number of rows");
1608  queso_require_equal_to_msg(this->numCols(), (mat1.numCols() * mat2.numCols()), "inconsistent number of columns");
1609 
1610  for (unsigned int rowId1 = 0; rowId1 < mat1.numRowsLocal(); ++rowId1) {
1611  for (unsigned int colId1 = 0; colId1 < mat1.numCols(); ++colId1) {
1612  double multiplicativeFactor = mat1(rowId1,colId1);
1613  unsigned int targetRowId = rowId1 * mat2.numRowsLocal();
1614  unsigned int targetColId = colId1 * mat2.numCols();
1615  for (unsigned int rowId2 = 0; rowId2 < mat2.numRowsLocal(); ++rowId2) {
1616  for (unsigned int colId2 = 0; colId2 < mat2.numCols(); ++colId2) {
1617  (*this)(targetRowId + rowId2, targetColId + colId2) = multiplicativeFactor * mat2(rowId2,colId2);
1618  }
1619  }
1620  }
1621  }
1622  return;
1623 }
1624 
1625 // ---------------------------------------------------
1626 // added 2/28/13
1627 void
1628 TeuchosMatrix::fillWithTensorProduct(const TeuchosMatrix& mat1, const TeuchosVector& vec2)
1629 {
1630  queso_require_equal_to_msg(this->numRowsLocal(), (mat1.numRowsLocal() * vec2.sizeLocal()), "inconsistent local number of rows");
1631  queso_require_equal_to_msg(this->numCols(), (mat1.numCols() * 1), "inconsistent number of columns");
1632 
1633  for (unsigned int rowId1 = 0; rowId1 < mat1.numRowsLocal(); ++rowId1) {
1634  for (unsigned int colId1 = 0; colId1 < mat1.numCols(); ++colId1) {
1635  double multiplicativeFactor = mat1(rowId1,colId1);
1636  unsigned int targetRowId = rowId1 * vec2.sizeLocal();
1637  unsigned int targetColId = colId1 * 1;
1638  for (unsigned int rowId2 = 0; rowId2 < vec2.sizeLocal(); ++rowId2) {
1639  for (unsigned int colId2 = 0; colId2 < 1; ++colId2) {
1640  (*this)(targetRowId + rowId2, targetColId + colId2) = multiplicativeFactor * vec2[rowId2];
1641  }
1642  }
1643  }
1644  }
1645  return;
1646 }
1647 
1648 
1649 // Miscellaneous methodos ----------------------------
1650 // ---------------------------------------------------
1651 
1652 void
1653 TeuchosMatrix::mpiSum( const MpiComm& comm, TeuchosMatrix& M_global ) const
1654 {
1655  // Sanity Checks
1656  queso_require_equal_to_msg(this->numRowsLocal(),
1657  M_global.numRowsLocal(),
1658  "local and global matrices incompatible");
1659  queso_require_equal_to_msg(this->numCols(),
1660  M_global.numCols(),
1661  "local and global matrices incompatible");
1662 
1663  /* TODO: Probably a better way to handle this unpacking/packing of data */
1664  int size = M_global.numRowsLocal()*M_global.numCols();
1665  std::vector<double> local( size, 0.0 );
1666  std::vector<double> global( size, 0.0 );
1667 
1668  int k;
1669  for( unsigned int i = 0; i < this->numRowsLocal(); i++ ) {
1670  for( unsigned int j = 0; j < this->numCols(); j++ ) {
1671  k = i + j*M_global.numCols();
1672  local[k] = (*this)(i,j);
1673  }
1674  }
1675 
1676  comm.Allreduce<double>(&local[0], &global[0], size, RawValue_MPI_SUM,
1677  "TeuchosMatrix::mpiSum()",
1678  "failed MPI.Allreduce()");
1679 
1680  for( unsigned int i = 0; i < this->numRowsLocal(); i++ ) {
1681  for( unsigned int j = 0; j < this->numCols(); j++ ) {
1682  k = i + j*M_global.numCols();
1683  M_global(i,j) = global[k];
1684  }
1685  }
1686 
1687  return;
1688 }
1689 
1690 //--------------------------------------------------------
1691 // tested 2/28/13
1692 void
1693 TeuchosMatrix::matlabLinearInterpExtrap(
1694  const TeuchosVector& x1Vec,
1695  const TeuchosMatrix& y1Mat,
1696  const TeuchosVector& x2Vec)
1697 {
1698  queso_require_greater_msg(x1Vec.sizeLocal(), 1, "invalid 'x1' size");
1699 
1700  queso_require_equal_to_msg(x1Vec.sizeLocal(), y1Mat.numRowsLocal(), "invalid 'x1' and 'y1' sizes");
1701 
1702  queso_require_equal_to_msg(x2Vec.sizeLocal(), this->numRowsLocal(), "invalid 'x2' and 'this' sizes");
1703 
1704  queso_require_equal_to_msg(y1Mat.numCols(), this->numCols(), "invalid 'y1' and 'this' sizes");
1705 
1706  TeuchosVector y1Vec(x1Vec);
1707  TeuchosVector y2Vec(x2Vec);
1708  for (unsigned int colId = 0; colId < y1Mat.numCols(); ++colId) {
1709  y1Mat.getColumn(colId,y1Vec);
1710  y2Vec.matlabLinearInterpExtrap(x1Vec,y1Vec,x2Vec);
1711  this->setColumn(colId,y2Vec);
1712  }
1713 
1714  return;
1715 }
1716 
1717 // I/O methodos --------------------------------------
1718 // ---------------------------------------------------
1719 // Kemelli 12/05/12 - tested
1720 void
1721 TeuchosMatrix::print(std::ostream& os) const
1722 {
1723  unsigned int nRows = this->numRowsLocal();
1724  unsigned int nCols = this->numCols();
1725 
1726  if (m_printHorizontally) {
1727  for (unsigned int i = 0; i < nRows; ++i) {
1728  for (unsigned int j = 0; j < nCols; ++j) {
1729  os << (*this)(i,j)
1730  << " ";
1731  }
1732  if (i != (nRows-1)) os << "; ";
1733  }
1734  //os << std::endl;
1735  }
1736  else {
1737  for (unsigned int i = 0; i < nRows; ++i) {
1738  for (unsigned int j = 0; j < nCols; ++j) {
1739  os << (*this)(i,j)
1740  << " ";
1741  }
1742  os << std::endl;
1743  }
1744  }
1745 
1746  return;
1747 }
1748 // ---------------------------------------------------
1749 
1750 void
1751 TeuchosMatrix::subReadContents(
1752  const std::string& fileName,
1753  const std::string& fileType,
1754  const std::set<unsigned int>& allowedSubEnvIds)
1755 {
1756  queso_require_greater_equal_msg(m_env.subRank(), 0, "unexpected subRank");
1757 
1758  queso_require_less_equal_msg(this->numOfProcsForStorage(), 1, "implemented just for sequential vectors for now");
1759 
1760  FilePtrSetStruct filePtrSet;
1761  if (m_env.openInputFile(fileName,
1762  fileType, // "m or hdf"
1763  allowedSubEnvIds,
1764  filePtrSet)) {
1765 
1766  // palms
1767  unsigned int nRowsLocal = this->numRowsLocal();
1768 
1769  // In the logic below, the id of a line' begins with value 0 (zero)
1770  unsigned int idOfMyFirstLine = 1;
1771  unsigned int idOfMyLastLine = nRowsLocal;
1772  unsigned int nCols = this->numCols();
1773 
1774  // Read number of matrix rows in the file by taking care of the first line,
1775  // which resembles something like 'variable_name = zeros(n_rows,n_cols);'
1776  std::string tmpString;
1777 
1778  // Read 'variable name' string
1779  *filePtrSet.ifsVar >> tmpString;
1780 
1781  // Read '=' sign
1782  *filePtrSet.ifsVar >> tmpString;
1783 
1784  queso_require_equal_to_msg(tmpString, "=", "string should be the '=' sign");
1785 
1786  // Read 'zeros(n_rows,n_cols)' string
1787  *filePtrSet.ifsVar >> tmpString;
1788 
1789  unsigned int posInTmpString = 6;
1790 
1791  // Isolate 'n_rows' in a string
1792  char nRowsString[tmpString.size()-posInTmpString+1];
1793  unsigned int posInRowsString = 0;
1794  do {
1795  queso_require_less_msg(posInTmpString, tmpString.size(), "symbol ',' not found in first line of file");
1796  nRowsString[posInRowsString++] = tmpString[posInTmpString++];
1797  } while (tmpString[posInTmpString] != ',');
1798  nRowsString[posInRowsString] = '\0';
1799 
1800  // Isolate 'n_cols' in a string
1801  posInTmpString++; // Avoid reading ',' char
1802  char nColsString[tmpString.size()-posInTmpString+1];
1803  unsigned int posInColsString = 0;
1804  do {
1805  queso_require_less_msg(posInTmpString, tmpString.size(), "symbol ')' not found in first line of file");
1806  nColsString[posInColsString++] = tmpString[posInTmpString++];
1807  } while (tmpString[posInTmpString] != ')');
1808  nColsString[posInColsString] = '\0';
1809 
1810  // Convert 'n_rows' and 'n_cols' strings to numbers
1811  unsigned int numRowsInFile = (unsigned int) strtod(nRowsString,NULL);
1812  unsigned int numColsInFile = (unsigned int) strtod(nColsString,NULL);
1813  if (m_env.subDisplayFile()) {
1814  *m_env.subDisplayFile() << "In TeuchosMatrix::subReadContents()"
1815  << ": fullRank " << m_env.fullRank()
1816  << ", numRowsInFile = " << numRowsInFile
1817  << ", numColsInFile = " << numColsInFile
1818  << ", nRowsLocal = " << nRowsLocal
1819  << ", nCols = " << nCols
1820  << std::endl;
1821  }
1822 
1823  // Check if [num of rows in file] == [requested matrix row size]
1824  queso_require_equal_to_msg(numRowsInFile, nRowsLocal, "size of vec in file is not big enough");
1825 
1826  // Check if [num of cols in file] == [num cols in current matrix]
1827  queso_require_equal_to_msg(numColsInFile, nCols, "number of parameters of vec in file is different than number of parameters in this vec object");
1828 
1829  // Code common to any core in a communicator
1830  unsigned int maxCharsPerLine = 64*nCols; // Up to about 60 characters to represent each parameter value
1831 
1832  unsigned int lineId = 0;
1833  while (lineId < idOfMyFirstLine) {
1834  filePtrSet.ifsVar->ignore(maxCharsPerLine,'\n');
1835  lineId++;
1836  };
1837 
1838  if (m_env.subDisplayFile()) {
1839  *m_env.subDisplayFile() << "In TeuchosMatrix::subReadContents()"
1840  << ": beginning to read input actual data"
1841  << std::endl;
1842  }
1843 
1844  // Take care of initial part of the first data line,
1845  // which resembles something like 'variable_name = [value1 value2 ...'
1846  // Read 'variable name' string
1847  *filePtrSet.ifsVar >> tmpString;
1848  //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
1849 
1850  // Read '=' sign
1851  *filePtrSet.ifsVar >> tmpString;
1852  //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
1853  queso_require_equal_to_msg(tmpString, "=", "in core 0, string should be the '=' sign");
1854 
1855  // Take into account the ' [' portion
1856  std::streampos tmpPos = filePtrSet.ifsVar->tellg();
1857  filePtrSet.ifsVar->seekg(tmpPos+(std::streampos)2);
1858 
1859  if (m_env.subDisplayFile()) {
1860  *m_env.subDisplayFile() << "In TeuchosMatrix::subReadContents()"
1861  << ": beginning to read lines with numbers only"
1862  << ", lineId = " << lineId
1863  << ", idOfMyFirstLine = " << idOfMyFirstLine
1864  << ", idOfMyLastLine = " << idOfMyLastLine
1865  << std::endl;
1866  }
1867 
1868  double tmpRead;
1869  while (lineId <= idOfMyLastLine) {
1870  for (unsigned int j = 0; j < nCols; ++j) {
1871  *filePtrSet.ifsVar >> tmpRead;
1872  (*this)(lineId-idOfMyFirstLine,j) = tmpRead;
1873  }
1874  lineId++;
1875  };
1876 
1877  m_env.closeFile(filePtrSet,fileType);
1878  }
1879 
1880  return;
1881 }
1882 
1883 // ---------------------------------------------------
1884 void
1885 TeuchosMatrix::subWriteContents(
1886  const std::string& varNamePrefix,
1887  const std::string& fileName,
1888  const std::string& fileType,
1889  const std::set<unsigned int>& allowedSubEnvIds) const
1890 {
1891  queso_require_greater_equal_msg(m_env.subRank(), 0, "unexpected subRank");
1892 
1893  queso_require_less_equal_msg(this->numOfProcsForStorage(), 1, "implemented just for sequential vectors for now");
1894 
1895  FilePtrSetStruct filePtrSet;
1896  if (m_env.openOutputFile(fileName,
1897  fileType, // "m or hdf"
1898  allowedSubEnvIds,
1899  false,
1900  filePtrSet)) {
1901  unsigned int nRows = this->numRowsLocal();
1902  unsigned int nCols = this->numCols();
1903  *filePtrSet.ofsVar << varNamePrefix << "_sub" << m_env.subIdString() << " = zeros(" << nRows
1904  << "," << nCols
1905  << ");"
1906  << std::endl;
1907  *filePtrSet.ofsVar << varNamePrefix << "_sub" << m_env.subIdString() << " = [";
1908 
1909  for (unsigned int i = 0; i < nRows; ++i) {
1910  for (unsigned int j = 0; j < nCols; ++j) {
1911  *filePtrSet.ofsVar << (*this)(i,j)
1912  << " ";
1913  }
1914  *filePtrSet.ofsVar << "\n";
1915  }
1916  *filePtrSet.ofsVar << "];\n";
1917 
1918  m_env.closeFile(filePtrSet,fileType);
1919  }
1920 
1921  return;
1922 }
1923 
1924 //====================================================
1925 // Private members
1926 //====================================================
1927 
1928 // ---------------------------------------------------
1929 // Kemelli 12/05/12 - tested
1930 void
1931 TeuchosMatrix::copy(const TeuchosMatrix& src) //dummy
1932 {
1933  // FIXME - should we be calling Matrix::base_copy here? - RHS
1934  this->resetLU();
1935  unsigned int i,j, nrows=src.numRowsLocal(), ncols=src.numCols();
1936 
1937  for(i=0; i< nrows ; i++)
1938  for (j = 0; j < ncols; j++)
1939  m_mat(i,j) = src(i,j);
1940 
1941  return;
1942 }
1943 
1944 // ---------------------------------------------------
1945 void
1946 TeuchosMatrix::resetLU()
1947 {
1948  if (m_LU.numCols() >0 || m_LU.numRows() > 0) {
1949  m_LU.reshape(0,0); //Kemelli, 12/06/12, dummy
1950  }
1951  if (m_inverse) {
1952  delete m_inverse;
1953  m_inverse = NULL;
1954  }
1955  if (m_svdColMap) {
1956  delete m_svdColMap;
1957  m_svdColMap = NULL;
1958  }
1959  if (m_svdUmat) {
1960  delete m_svdUmat;
1961  m_svdUmat = NULL;
1962  }
1963  if (m_svdSvec) {
1964  delete m_svdSvec;
1965  m_svdSvec = NULL;
1966  }
1967  if (m_svdVmat) {
1968  delete m_svdVmat;
1969  m_svdVmat = NULL;
1970  }
1971  if (m_svdVTmat) {
1972  delete m_svdVTmat;
1973  m_svdVTmat = NULL;
1974  }
1975  m_determinant = -INFINITY;
1976  m_lnDeterminant = -INFINITY;
1977 
1978  if (v_pivoting) { //Kemelli added 12/09/12
1979  free(v_pivoting);
1980  v_pivoting = NULL;
1981 
1982  }
1983  m_signum = 0;
1984  m_isSingular = false;
1985 
1986  return;
1987 }
1988 
1989 // ---------------------------------------------------
1990 // multiply this matrix by vector x and store in vector y
1991 // checked 12/10/12
1992 void
1993 TeuchosMatrix::multiply(const TeuchosVector& x, TeuchosVector& y) const
1994 {
1995  queso_require_equal_to_msg(this->numCols(), x.sizeLocal(), "matrix and x have incompatible sizes");
1996 
1997  queso_require_equal_to_msg(this->numRowsLocal(), y.sizeLocal(), "matrix and y have incompatible sizes");
1998 
1999  unsigned int sizeX = this->numCols();
2000  unsigned int sizeY = this->numRowsLocal();
2001  for (unsigned int i = 0; i < sizeY; ++i) {
2002  double value = 0.;
2003  for (unsigned int j = 0; j < sizeX; ++j) {
2004  value += (*this)(i,j)*x[j];
2005  }
2006  y[i] = value;
2007  }
2008 
2009  return;
2010 }
2011 
2012 // ---------------------------------------------------
2013 // Implemented(finally) and checked 1/10/13
2014 int
2015 TeuchosMatrix::internalSvd() const
2016 {
2017  if (m_svdColMap == NULL) {
2018  int nRows = (int) this->numRowsLocal();
2019  int nCols = (int) this->numCols();
2020  queso_require_greater_equal_msg(nRows, nCols, "LAPACK/Teuchos only supports cases where nRows >= nCols");
2021 
2022  m_svdColMap = new Map(this->numCols(),0,this->map().Comm()); // see 'VectorSpace<.,.>::newMap()'
2023  //in src/basic/src/TeuchosVectorSpace.C //old comment already existent in GslMatrix
2024  m_svdUmat = new TeuchosMatrix(*this); // Yes, 'this'
2025  m_svdSvec = new TeuchosVector(m_env,*m_svdColMap);
2026  m_svdVmat = new TeuchosMatrix(*m_svdSvec);
2027  m_svdVTmat = new TeuchosMatrix(*m_svdSvec);
2028 
2029  int minRowsCols, maxRowsCols;
2030 
2031  if (nRows>=nCols) { minRowsCols = nCols; maxRowsCols = nRows; } else { minRowsCols = nRows; maxRowsCols = nCols; }
2032 
2033  char jobu, jobvt;
2034  int lwork, info;
2035  double *work, *rwork;
2036 
2037  jobu = 'S';
2038  jobvt = 'S';
2039 
2040  lwork = 15*maxRowsCols; // Set up the work array, larger than needed.
2041  work = new double[lwork];
2042 
2043  int aux1= 5*minRowsCols+7, aux2= 2*maxRowsCols+2*minRowsCols+1;
2044  int aux_dim;
2045 
2046  if (aux1>=aux2) { aux_dim = minRowsCols*aux1; } else {aux_dim = minRowsCols*aux2; }
2047 
2048  rwork = new double[aux_dim];
2049 
2050  Teuchos::LAPACK<int, double> lapack;
2051 
2052  lapack.GESVD(jobu,jobvt,m_mat.numRows(),m_mat.numCols(),m_mat.values(),m_mat.stride(),
2053  m_svdSvec->values(),m_svdUmat->values(),m_svdUmat->stride(),m_svdVTmat->values(),
2054  m_svdVTmat->stride(),&work[0],lwork,&rwork[0],&info);
2055  }
2056  return 0;
2057 }
2058 
2059 
2060 //++++++++++++++++++++++++++++++++++++++++++++++++++++
2061 // Operators outside class definition
2062 //++++++++++++++++++++++++++++++++++++++++++++++++++++
2063 
2064 TeuchosMatrix operator*(double a, const TeuchosMatrix& mat)
2065 {
2066  TeuchosMatrix answer(mat);
2067  answer *= a;
2068  return answer;
2069 }
2070 
2071 // ---------------------------------------------------
2072 TeuchosVector operator*(const TeuchosMatrix& mat, const TeuchosVector& vec)
2073 {
2074  return mat.multiply(vec);
2075 }
2076 
2077 // ---------------------------------------------------
2078 TeuchosMatrix operator*(const TeuchosMatrix& m1, const TeuchosMatrix& m2)
2079 {
2080  unsigned int m1Rows = m1.numRowsLocal();
2081  unsigned int m1Cols = m1.numCols();
2082  unsigned int m2Rows = m2.numRowsLocal();
2083  unsigned int m2Cols = m2.numCols();
2084 
2085  queso_require_equal_to_msg(m1Cols, m2Rows, "different sizes m1Cols and m2Rows");
2086 
2087  TeuchosMatrix mat(m1.env(),m1.map(),m2Cols);
2088  unsigned int commonSize = m1Cols;
2089 
2090  for (unsigned int row1 = 0; row1 < m1Rows; ++row1) {
2091  for (unsigned int col2 = 0; col2 < m2Cols; ++col2) {
2092  double result = 0.;
2093  for (unsigned int k = 0; k < commonSize; ++k) {
2094  result += m1(row1,k)*m2(k,col2);
2095  }
2096  mat(row1,col2) = result;
2097  }
2098  }
2099 
2100  return mat;
2101 }
2102 
2103 // ---------------------------------------------------
2104 TeuchosMatrix operator+(const TeuchosMatrix& m1, const TeuchosMatrix& m2)
2105 {
2106  TeuchosMatrix answer(m1);
2107  answer += m2;
2108  return answer;
2109 }
2110 
2111 // ---------------------------------------------------
2112 TeuchosMatrix operator-(const TeuchosMatrix& m1, const TeuchosMatrix& m2)
2113 {
2114  TeuchosMatrix answer(m1);
2115  answer -= m2;
2116  return answer;
2117 }
2118 
2119 // ---------------------------------------------------
2120 TeuchosMatrix matrixProduct(const TeuchosVector& v1, const TeuchosVector& v2)
2121 {
2122  unsigned int nRows = v1.sizeLocal();
2123  unsigned int nCols = v2.sizeLocal();
2124  TeuchosMatrix answer(v1.env(),v1.map(),nCols);
2125 
2126  for (unsigned int i = 0; i < nRows; ++i) {
2127  double value1 = v1[i];
2128  for (unsigned int j = 0; j < nCols; ++j) {
2129  answer(i,j) = value1*v2[j];
2130  }
2131  }
2132 
2133  return answer;
2134 }
2135 
2136 // ---------------------------------------------------
2137 TeuchosMatrix leftDiagScaling(const TeuchosVector& vec, const TeuchosMatrix& mat)
2138 {
2139  unsigned int vSize = vec.sizeLocal();
2140  unsigned int mRows = mat.numRowsLocal();
2141  unsigned int mCols = mat.numCols();
2142 
2143  queso_require_equal_to_msg(vSize, mRows, "size of vector is different from the number of rows in matrix");
2144 
2145  queso_require_equal_to_msg(mCols, mRows, "routine currently works for square matrices only");
2146 
2147  TeuchosMatrix answer(mat);
2148  for (unsigned int i = 0; i < mRows; ++i) {
2149  double vecValue = vec[i];
2150  for (unsigned int j = 0; j < mCols; ++j) {
2151  answer(i,j) *= vecValue;
2152  }
2153  }
2154 
2155  return answer;
2156 }
2157 
2158 // ---------------------------------------------------
2159 TeuchosMatrix rightDiagScaling(const TeuchosMatrix& mat, const TeuchosVector& vec)
2160 {
2161  unsigned int vSize = vec.sizeLocal();
2162  unsigned int mRows = mat.numRowsLocal();
2163  unsigned int mCols = mat.numCols();
2164 
2165  queso_require_equal_to_msg(vSize, mCols, "size of vector is different from the number of cols in matrix");
2166 
2167  queso_require_equal_to_msg(mCols, mRows, "routine currently works for square matrices only");
2168 
2169  TeuchosMatrix answer(mat);
2170  for (unsigned int j = 0; j < mCols; ++j) {
2171  double vecValue = vec[j];
2172  for (unsigned int i = 0; i < mRows; ++i) {
2173  answer(i,j) *= vecValue;
2174  }
2175  }
2176 
2177  return answer;
2178 }
2179 
2180 // ---------------------------------------------------
2181 std::ostream&
2182 operator<<(std::ostream& os, const TeuchosMatrix& obj)
2183 {
2184  obj.print(os);
2185 
2186  return os;
2187 }
2188 
2189 } // End namespace QUESO
2190 
2191 #endif // ifdef QUESO_HAS_TRILINOS
QUESO::TeuchosMatrix::cwSet
void cwSet(double value)
Component-wise set all values to this with value.
Definition: TeuchosMatrix.C:1189
QUESO::TeuchosMatrix::operator-=
TeuchosMatrix & operator-=(const TeuchosMatrix &rhs)
Stores in this the coordinate-wise subtraction of this by rhs.
Definition: TeuchosMatrix.C:227
QUESO::TeuchosMatrix::numRowsLocal
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: TeuchosMatrix.C:274
QUESO::TeuchosMatrix::fillWithBlocksHorizontally
void fillWithBlocksHorizontally(const std::vector< const TeuchosMatrix * > &matrices)
This function fills this matrix horizontally with const block matrices.
Definition: TeuchosMatrix.C:1499
QUESO::TeuchosMatrix::TeuchosMatrix
TeuchosMatrix()
Default Constructor.
QUESO::TeuchosMatrix::normFrob
double normFrob() const
Returns the Frobenius norm of this matrix.
Definition: TeuchosMatrix.C:541
QUESO::TeuchosMatrix::svd
int svd(TeuchosMatrix &matU, TeuchosVector &vecS, TeuchosMatrix &matVt) const
Checks for the dimension of this matrix, matU, VecS and matVt, and calls the protected routine intern...
Definition: TeuchosMatrix.C:617
QUESO::TeuchosMatrix::stride
int stride()
Returns the stride between the columns of this matrix in memory.
Definition: TeuchosMatrix.C:306
QUESO::Matrix::m_map
const Map m_map
Mapping variable.
Definition: Matrix.h:121
QUESO::Map
A class for partitioning vectors and matrices.
Definition: Map.h:49
QUESO::Matrix::env
const BaseEnvironment & env() const
Definition: Matrix.C:47
QUESO::TeuchosMatrix::svdMatV
const TeuchosMatrix & svdMatV() const
This function calls private member TeuchosMatrix::internalSvd() to set a N-by-N orthogonal square mat...
Definition: TeuchosMatrix.C:649
QUESO::TeuchosMatrix::m_svdColMap
Map * m_svdColMap
Mapping for matrices involved in the singular value decomposition (svd) routine.
Definition: TeuchosMatrix.h:361
QUESO::TeuchosMatrix::invertMultiplyForceLU
void invertMultiplyForceLU(const TeuchosVector &b, TeuchosVector &x) const
This function calculates the inverse of this matrix, multiplies it with vector b and stores the resul...
Definition: TeuchosMatrix.C:941
QUESO::TeuchosVector::sizeLocal
unsigned int sizeLocal() const
Returns the length of this vector.
Definition: TeuchosVector.C:249
QUESO::TeuchosMatrix::eigen
void eigen(TeuchosVector &eigenValues, TeuchosMatrix *eigenVectors) const
This function computes the eigenvalues of a real symmetric matrix.
Definition: TeuchosMatrix.C:1039
QUESO::BaseEnvironment
This (virtual) class sets up the environment underlying the use of the QUESO library by an executable...
Definition: Environment.h:198
QUESO::TeuchosMatrix::m_LU
Teuchos::SerialDenseMatrix< int, double > m_LU
Teuchos matrix for the LU decomposition of m_mat.
Definition: TeuchosMatrix.h:355
QUESO::TeuchosMatrix::m_svdVmat
TeuchosMatrix * m_svdVmat
m_svdVmat stores the N-by-N orthogonal square matrix V after the singular value decomposition of a ma...
Definition: TeuchosMatrix.h:370
QUESO::TeuchosVector::matlabLinearInterpExtrap
void matlabLinearInterpExtrap(const TeuchosVector &xVec, const TeuchosVector &yVec, const TeuchosVector &xiVec)
Reproduces MATLAB linear inter/extra-polation.
Definition: TeuchosVector.C:834
QUESO::TeuchosMatrix::m_svdUmat
TeuchosMatrix * m_svdUmat
m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix...
Definition: TeuchosMatrix.h:364
QUESO::TeuchosMatrix::m_lnDeterminant
double m_lnDeterminant
The natural logarithm of the determinant of this matrix.
Definition: TeuchosMatrix.h:379
QUESO::TeuchosMatrix::~TeuchosMatrix
~TeuchosMatrix()
Destructor.
Definition: TeuchosMatrix.C:174
QUESO::TeuchosMatrix::svdMatU
const TeuchosMatrix & svdMatU() const
This function calls private member TeuchosMatrix::internalSvd() to set a M-by-N orthogonal matrix U o...
Definition: TeuchosMatrix.C:639
QUESO::TeuchosMatrix::chol
int chol()
Computes Cholesky factorization of this, a real symmetric positive definite matrix.
Definition: TeuchosMatrix.C:568
QUESO::TeuchosMatrix::setColumn
void setColumn(const unsigned int column_num, const TeuchosVector &column)
This function copies vector column into the column_num-th column of this matrix.
Definition: TeuchosMatrix.C:1255
QUESO::TeuchosMatrix::operator*=
TeuchosMatrix & operator*=(double a)
Stores in this the coordinate-wise multiplication of this and a.
Definition: TeuchosMatrix.C:192
QUESO::Matrix::numOfProcsForStorage
unsigned int numOfProcsForStorage() const
Definition: Matrix.C:62
QUESO::operator<<
std::ostream & operator<<(std::ostream &os, const SequenceStatisticalOptions &obj)
Definition: SequenceStatisticalOptions.C:1457
QUESO::BaseEnvironment::closeFile
void closeFile(FilePtrSetStruct &filePtrSet, const std::string &fileType) const
Closes the file.
Definition: Environment.C:1084
QUESO::queso_require_not_equal_to_msg
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix queso_require_not_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the absence of an options input file"))
Definition: MonteCarloSGOptions.C:417
QUESO::TeuchosMatrix::inverse
TeuchosMatrix inverse() const
This function calculates the inverse of this matrix (square).
Definition: TeuchosMatrix.C:392
QUESO::BaseEnvironment::subRank
int subRank() const
Returns the rank of the MPI process in the sub-communicator subComm()
Definition: Environment.C:287
QUESO::TeuchosMatrix::numRowsGlobal
unsigned int numRowsGlobal() const
Returns the global row dimension of this matrix.
Definition: TeuchosMatrix.C:282
QUESO::operator-
GslMatrix operator-(const GslMatrix &m1, const GslMatrix &m2)
Definition: GslMatrix.C:2029
QUESO::Matrix::m_printHorizontally
bool m_printHorizontally
Flag for either or not print this matrix.
Definition: Matrix.h:131
QUESO::TeuchosMatrix::subWriteContents
void subWriteContents(const std::string &varNamePrefix, const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds) const
Write contents of subenvironment in file fileName.
Definition: TeuchosMatrix.C:1885
QUESO::TeuchosMatrix::fillWithTranspose
void fillWithTranspose(const TeuchosMatrix &mat)
This function stores the transpose of this matrix into this matrix.
Definition: TeuchosMatrix.C:1420
QUESO::TeuchosMatrix::m_inverse
TeuchosMatrix * m_inverse
Stores the inverse of this matrix.
Definition: TeuchosMatrix.h:358
QUESO::BaseEnvironment::checkingLevel
unsigned int checkingLevel() const
Access function to private attribute m_checkingLevel.
Definition: Environment.C:464
QUESO::UQ_MATRIX_IS_NOT_POS_DEFINITE_RC
const int UQ_MATRIX_IS_NOT_POS_DEFINITE_RC
Definition: Defines.h:100
QUESO::TeuchosMatrix::setRow
void setRow(const unsigned int row_num, const TeuchosVector &row)
This function copies vector row into the row_num-th column of this matrix.
Definition: TeuchosMatrix.C:1301
QUESO::TeuchosMatrix::zeroUpper
void zeroUpper(bool includeDiagonal=false)
This function sets all the entries bellow the main diagonal (inclusive or not) of this matrix to zero...
Definition: TeuchosMatrix.C:1347
QUESO::TeuchosMatrix::largestEigen
void largestEigen(double &eigenValue, TeuchosVector &eigenVector) const
This function finds largest eigenvalue, namely eigenValue, of this matrix and its corresponding eigen...
Definition: TeuchosMatrix.C:1097
QUESO::TeuchosMatrix::multiply
TeuchosVector multiply(const TeuchosVector &x) const
This function multiplies this matrix by vector x and returns a vector.
Definition: TeuchosMatrix.C:765
QUESO::TeuchosMatrix::getRow
void getRow(const unsigned int row_num, TeuchosVector &row) const
This function gets the row_num-th column of this matrix and stores it into vector row...
Definition: TeuchosMatrix.C:1272
QUESO::BaseEnvironment::worldRank
int worldRank() const
Returns the same thing as fullRank()
Definition: Environment.C:262
QUESO::TeuchosMatrix::determinant
double determinant() const
Calculates the determinant of this matrix.
Definition: TeuchosMatrix.C:442
QUESO::TeuchosMatrix::getColumn
void getColumn(const unsigned int column_num, TeuchosVector &column) const
This function gets the column_num-th column of this matrix and stores it into vector column...
Definition: TeuchosMatrix.C:1220
QUESO::BaseEnvironment::fullRank
int fullRank() const
Returns the rank of the MPI process in QUESO&#39;s full communicator.
Definition: Environment.C:268
QUESO::TeuchosMatrix::normMax
double normMax() const
Returns the Frobenius norm of this matrix.
Definition: TeuchosMatrix.C:548
QUESO::TeuchosMatrix::internalSvd
int internalSvd() const
This function factorizes the M-by-N matrix A into the singular value decomposition A = U S V^T for M ...
Definition: TeuchosMatrix.C:2015
QUESO::BaseEnvironment::openInputFile
bool openInputFile(const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds, FilePtrSetStruct &filePtrSet) const
Opens an input file.
Definition: Environment.C:896
QUESO::TeuchosMatrix::zeroLower
void zeroLower(bool includeDiagonal=false)
This function sets all the entries above the main diagonal (inclusive or not) of this matrix to zero...
Definition: TeuchosMatrix.C:1318
QUESO::TeuchosMatrix::values
double * values()
Returns a pointer to the first element of this matrix.
Definition: TeuchosMatrix.C:298
QUESO::TeuchosMatrix::operator+=
TeuchosMatrix & operator+=(const TeuchosMatrix &rhs)
Stores in this the coordinate-wise addition of this and rhs.
Definition: TeuchosMatrix.C:212
QUESO::TeuchosMatrix::v_pivoting
int * v_pivoting
The pivoting vector of a LU decomposition.
Definition: TeuchosMatrix.h:384
QUESO::Matrix
Class for matrix operations (virtual).
Definition: Matrix.h:46
QUESO::TeuchosMatrix::m_svdSvec
TeuchosVector * m_svdSvec
m_svdSvec stores the diagonal of the N-by-N diagonal matrix of singular values S after the singular v...
Definition: TeuchosMatrix.h:367
QUESO::operator*
GslMatrix operator*(double a, const GslMatrix &mat)
Definition: GslMatrix.C:1982
QUESO::TeuchosMatrix::matlabLinearInterpExtrap
void matlabLinearInterpExtrap(const TeuchosVector &x1Vec, const TeuchosMatrix &y1Mat, const TeuchosVector &x2Vec)
Definition: TeuchosMatrix.C:1693
QUESO::TeuchosMatrix::filterSmallValues
void filterSmallValues(double thresholdValue)
This function sets to zero (filters) all entries of this matrix which are smaller than thresholdValue...
Definition: TeuchosMatrix.C:1375
QUESO::TeuchosMatrix::m_mat
Teuchos::SerialDenseMatrix< int, double > m_mat
Teuchos matrix, also referred to as this matrix.
Definition: TeuchosMatrix.h:352
QUESO::TeuchosMatrix::operator=
TeuchosMatrix & operator=(const TeuchosMatrix &rhs)
Copies values from matrix rhs to this.
Definition: TeuchosMatrix.C:183
QUESO::TeuchosMatrix::transpose
TeuchosMatrix transpose() const
This function calculates the transpose of this matrix (square).
Definition: TeuchosMatrix.C:372
QUESO::TeuchosVector
Class for vector operations using Teuchos (Trilinos).
Definition: TeuchosVector.h:55
QUESO::Vector::env
const BaseEnvironment & env() const
Definition: Vector.C:54
QUESO::TeuchosMatrix::resetLU
void resetLU()
In this function resets the LU decomposition of this matrix, as well as deletes the private member po...
Definition: TeuchosMatrix.C:1946
QUESO::TeuchosMatrix::m_svdVTmat
TeuchosMatrix * m_svdVTmat
m_svdVmatT stores the transpose of N-by-N orthogonal square matrix V, namely V^T, after the singular ...
Definition: TeuchosMatrix.h:373
QUESO::MpiComm::Allreduce
void Allreduce(void *sendbuf, void *recvbuf, int count, RawType_MPI_Datatype datatype, RawType_MPI_Op op, const char *whereMsg, const char *whatMsg) const
Combines values from all processes and distributes the result back to all processes.
Definition: MpiComm.C:143
QUESO::Map::NumGlobalElements
int NumGlobalElements() const
Returns the total number of elements across all processors.
Definition: Map.C:85
QUESO::TeuchosMatrix::m_isSingular
bool m_isSingular
Indicates whether or not this matrix is singular.
Definition: TeuchosMatrix.h:390
QUESO::FilePtrSetStruct::ifsVar
std::ifstream * ifsVar
Provides a stream interface to read data from files.
Definition: Environment.h:89
QUESO::Vector::map
const Map & map() const
Definition: Vector.C:61
QUESO::TeuchosMatrix::fillWithBlocksVertically
void fillWithBlocksVertically(const std::vector< const TeuchosMatrix * > &matrices)
This function fills this matrix vertically with const block matrices.
Definition: TeuchosMatrix.C:1553
QUESO::TeuchosMatrix::m_determinant
double m_determinant
The determinant of this matrix.
Definition: TeuchosMatrix.h:376
QUESO::TeuchosVector::cwSet
void cwSet(double value)
Component-wise sets all values to this with value.
Definition: TeuchosVector.C:472
QUESO::matrixProduct
GslMatrix matrixProduct(const GslVector &v1, const GslVector &v2)
Definition: GslMatrix.C:2036
QUESO::TeuchosMatrix::fillWithTensorProduct
void fillWithTensorProduct(const TeuchosMatrix &mat1, const TeuchosMatrix &mat2)
This function calculates the tensor product of matrices mat1 and mat2 and stores it in this matrix...
Definition: TeuchosMatrix.C:1605
QUESO::TeuchosMatrix::max
double max() const
Returns the maximum element value of the matrix.
Definition: TeuchosMatrix.C:315
QUESO::TeuchosMatrix::mpiSum
void mpiSum(const MpiComm &comm, TeuchosMatrix &M_global) const
Definition: TeuchosMatrix.C:1653
dim
int dim
Definition: ann_test.cpp:472
QUESO::queso_require_equal_to_msg
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
Definition: MonteCarloSGOptions.C:460
QUESO::TeuchosMatrix::lnDeterminant
double lnDeterminant() const
Calculates the ln(determinant) of this matrix.
Definition: TeuchosMatrix.C:491
QUESO::TeuchosVector::norm2
double norm2() const
Returns the 2-norm (Euclidean norm) of the vector.
Definition: TeuchosVector.C:367
QUESO::Matrix::base_copy
virtual void base_copy(const Matrix &src)
Copies base data from matrix src to this matrix.
Definition: Matrix.C:99
QUESO::TeuchosMatrix::numCols
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: TeuchosMatrix.C:290
QUESO::MpiComm
The QUESO MPI Communicator Class.
Definition: MpiComm.h:203
QUESO::TeuchosMatrix::print
void print(std::ostream &os) const
Print method. Defines the behavior of the ostream &lt;&lt; operator inherited from the Object class...
Definition: TeuchosMatrix.C:1721
QUESO::TeuchosMatrix::subReadContents
void subReadContents(const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds)
Read contents of subenvironment from file fileName.
Definition: TeuchosMatrix.C:1751
QUESO::TeuchosMatrix::operator/=
TeuchosMatrix & operator/=(double a)
Stores in this the coordinate-wise division of this by a.
Definition: TeuchosMatrix.C:203
QUESO::TeuchosMatrix::copy
void copy(const TeuchosMatrix &src)
In this function this matrix receives a copy of matrix src.
Definition: TeuchosMatrix.C:1931
QUESO::TeuchosMatrix::operator()
double & operator()(unsigned int i, unsigned int j)
Element access method (non-const).
Definition: TeuchosMatrix.C:244
QUESO::BaseEnvironment::displayVerbosity
unsigned int displayVerbosity() const
Definition: Environment.C:450
QUESO::BaseEnvironment::subIdString
const std::string & subIdString() const
Access to the attribute m_subIdString; which stores the string for the sub-environment, and it will be used, for instance, to create the output files for each sub-environment.
Definition: Environment.C:348
QUESO::FilePtrSetStruct::ofsVar
std::ofstream * ofsVar
Provides a stream interface to write data to files.
Definition: Environment.h:86
QUESO::leftDiagScaling
GslMatrix leftDiagScaling(const GslVector &vec, const GslMatrix &mat)
Definition: GslMatrix.C:2052
QUESO::Matrix::m_env
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
QUESO::TeuchosMatrix::svdSolve
int svdSolve(const TeuchosVector &rhsVec, TeuchosVector &solVec) const
This function solves the system A x = b using the singular value decomposition (U, S, V) of A which must have been computed previously with TeuchosMatrix::svd (x=solVec, b=rhsVec).
Definition: TeuchosMatrix.C:667
QUESO::Matrix::map
const Map & map() const
Definition: Matrix.C:54
QUESO::FilePtrSetStruct
Struct for handling data input and output from files.
Definition: Environment.h:77
QUESO::TeuchosMatrix::rank
unsigned int rank(double absoluteZeroThreshold, double relativeZeroThreshold) const
This function returns the number of singular values of this matrix (rank).
Definition: TeuchosMatrix.C:335
QUESO::rightDiagScaling
GslMatrix rightDiagScaling(const GslMatrix &mat, const GslVector &vec)
Definition: GslMatrix.C:2073
QUESO::TeuchosMatrix::smallestEigen
void smallestEigen(double &eigenValue, TeuchosVector &eigenVector) const
This function finds smallest eigenvalue, namely eigenValue, of this matrix and its corresponding eige...
Definition: TeuchosMatrix.C:1141
QUESO::TeuchosMatrix::invertMultiply
void invertMultiply(const TeuchosVector &b, TeuchosVector &x) const
This function calculates the inverse of this matrix, multiplies it with vector b and stores the resul...
Definition: TeuchosMatrix.C:791
QUESO::Matrix::m_inDebugMode
bool m_inDebugMode
Flag for either or not QUESO is in debug mode.
Definition: Matrix.h:134
QUESO::BaseEnvironment::openOutputFile
bool openOutputFile(const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds, bool writeOver, FilePtrSetStruct &filePtrSet) const
Opens an output file for each sub-environment that was chosen to send data to the file...
Definition: Environment.C:521
QUESO::TeuchosMatrix
Class for matrix operations using Teuchos (Trilinos).
Definition: TeuchosMatrix.h:53
QUESO::operator+
GslMatrix operator+(const GslMatrix &m1, const GslMatrix &m2)
Definition: GslMatrix.C:2022
QUESO::BaseEnvironment::subDisplayFile
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
QUESO::TeuchosMatrix::filterLargeValues
void filterLargeValues(double thresholdValue)
This function sets to zero (filters) all entries of this matrix which are greater than thresholdValue...
Definition: TeuchosMatrix.C:1397
QUESO::TeuchosMatrix::fillWithBlocksDiagonally
void fillWithBlocksDiagonally(const std::vector< const TeuchosMatrix * > &matrices)
This function fills this matrix diagonally with const block matrices.
Definition: TeuchosMatrix.C:1439
Generated on Tue Jun 5 2018 19:48:55 for queso-0.57.1 by  doxygen 1.8.5