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