queso-0.56.0
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-2015 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::leftDiagScaling
GslMatrix leftDiagScaling(const GslVector &vec, const GslMatrix &mat)
Definition: GslMatrix.C:2052
work
the intent is to exercise the right to control the distribution of derivative or collective works based on the Library In mere aggregation of another work not based on the Library with the you must alter all the notices that refer to this so that they refer to the ordinary GNU General Public instead of to this it is irreversible for that so the ordinary GNU General Public License applies to all subsequent copies and derivative works made from that copy This option is useful when you wish to copy part of the code of the Library into a program that is not a library You may copy and distribute the which must be distributed under the terms of Sections and above on a medium customarily used for software interchange If distribution of object code is made by offering access to copy from a designated then offering equivalent access to copy the source code from the same place satisfies the requirement to distribute the source even though third parties are not compelled to copy the source along with the object code A program that contains no derivative of any portion of the but is designed to work with the Library by being compiled or linked with is called a work that uses the Library Such a work
Definition: License.txt:237
QUESO::UQ_MATRIX_IS_NOT_POS_DEFINITE_RC
const int UQ_MATRIX_IS_NOT_POS_DEFINITE_RC
Definition: Defines.h:97
RawValue_MPI_SUM
#define RawValue_MPI_SUM
Definition: MpiComm.h:71
queso_require_greater_equal_msg
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:78
QUESO::operator<<
std::ostream & operator<<(std::ostream &os, const BaseEnvironment &obj)
QUESO::operator+
GslMatrix operator+(const GslMatrix &m1, const GslMatrix &m2)
Definition: GslMatrix.C:2022
copy
and that you are informed that you can do these things To protect your we need to make restrictions that forbid distributors to deny you these rights or to ask you to surrender these rights These restrictions translate to certain responsibilities for you if you distribute copies of the library or if you modify it For if you distribute copies of the whether gratis or for a you must give the recipients all the rights that we gave you You must make sure that receive or can get the source code If you link other code with the you must provide complete object files to the so that they can relink them with the library after making changes to the library and recompiling it And you must show them these terms so they know their rights We protect your rights with a two step which gives you legal permission to copy
Definition: License.txt:60
dim
int dim
Definition: ann2fig.cpp:81
UQ_RC_MACRO
#define UQ_RC_MACRO(macroIRc, givenRank, where, what, retValue)
Definition: Defines.h:137
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:1155
queso_require_equal_to_msg
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:73
queso_require_less_msg
#define queso_require_less_msg(expr1, expr2, msg)
Definition: asserts.h:75
k
int k
Definition: ann_sample.cpp:53
QUESO::GslMatrix::numRowsLocal
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
queso_require_greater_msg
#define queso_require_greater_msg(expr1, expr2, msg)
Definition: asserts.h:76
QUESO::operator-
GslMatrix operator-(const GslMatrix &m1, const GslMatrix &m2)
Definition: GslMatrix.C:2029
queso_require_not_equal_to_msg
#define queso_require_not_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:74
queso_require_less_equal_msg
#define queso_require_less_equal_msg(expr1, expr2, msg)
Definition: asserts.h:77
queso_require_msg
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
QUESO::rightDiagScaling
GslMatrix rightDiagScaling(const GslMatrix &mat, const GslVector &vec)
Definition: GslMatrix.C:2073
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::operator*
GslMatrix operator*(double a, const GslMatrix &mat)
Definition: GslMatrix.C:1982
QUESO::matrixProduct
GslMatrix matrixProduct(const GslVector &v1, const GslVector &v2)
Definition: GslMatrix.C:2036
Generated on Tue Nov 29 2016 10:53:10 for queso-0.56.0 by  doxygen 1.8.5