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