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