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