queso-0.55.0
GcmJointInfo.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2015 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <queso/GcmJointInfo.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
28 
29 namespace QUESO {
30 
31 template <class S_V,class S_M,class D_V,class D_M,class P_V,class P_M,class Q_V,class Q_M>
32 GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::GcmJointInfo(
33  const GpmsaComputerModelOptions& gcmOptionsObj,
34  bool allOutputsAreScalar,
35  const GcmSimulationInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>& s,
36  const GcmExperimentInfo<S_V,S_M,D_V,D_M,P_V,P_M>& e)
37  :
38  m_env (s.m_env),
39  m_unique_u_space (m_env, "unique_u_", s.m_paper_p_eta, NULL),
40  m_Smat_u_asterisk_u_asterisk (m_unique_u_space.zeroVector()),
41  m_u_size (e.m_paper_n * s.m_paper_p_eta),
42  m_u_space (m_env, "u_", m_u_size, NULL),
43  m_Rmat_u_is (s.m_paper_p_eta, (D_M*) NULL), // to be deleted on destructor
44  m_Smat_u_is (s.m_paper_p_eta, (D_M*) NULL), // to be deleted on destructor
45  m_Rmat_uw_is (s.m_paper_p_eta, (D_M*) NULL), // to be deleted on destructor
46  m_Smat_uw_is (s.m_paper_p_eta, (D_M*) NULL), // to be deleted on destructor
47  m_Smat_uw (m_env, m_u_space.map(),s.m_w_size),
48  m_Smat_uw_t (m_env,s.m_w_space.map(), m_u_size),
49  m_Rmat_u_hat_u_asterisk_is (s.m_paper_p_eta, (D_M*) NULL), // to be deleted on destructor
50  m_Smat_u_hat_u_asterisk_is (s.m_paper_p_eta, (D_M*) NULL), // to be deleted on destructor
51  m_Smat_u_hat_u_asterisk (m_env, m_u_space.map(), s.m_paper_p_eta),
52  m_Smat_u_hat_u_asterisk_t (m_env, m_unique_u_space.map(), m_u_size),
53  m_Rmat_w_hat_u_asterisk_is (s.m_paper_p_eta, (D_M*) NULL), // to be deleted on destructor
54  m_Smat_w_hat_u_asterisk_is (s.m_paper_p_eta, (D_M*) NULL), // to be deleted on destructor
55  m_Smat_w_hat_u_asterisk (m_env, s.m_w_space.map(), s.m_paper_p_eta),
56  m_Smat_w_hat_u_asterisk_t (m_env, m_unique_u_space.map(), s.m_w_size),
57  m_vu_size (e.m_v_size + m_u_size),
58  m_vu_space (m_env, "vu_", m_vu_size, NULL),
59  m_unique_vu_space (m_env, "unique_vu_", e.m_paper_p_delta + s.m_paper_p_eta, NULL),
60  m_predVU_counter (MiscUintDebugMessage(0,NULL)),
61  m_predVU_summingRVs_unique_vu_meanVec (m_unique_vu_space.zeroVector()),
62  m_predVU_summingRVs_mean_of_unique_vu_covMatrices(m_unique_vu_space.zeroVector()),
63  m_predVU_summingRVs_covMatrix_of_unique_vu_means (m_unique_vu_space.zeroVector()),
64  m_predVU_summingRVs_corrMatrix_of_unique_vu_means(m_unique_vu_space.zeroVector()),
65  m_omega_size (e.m_paper_n_y + s.m_eta_size),
66  m_omega_space (m_env, "omega_", m_omega_size, NULL),
67  m_Zvec_hat_vu (m_vu_space.zeroVector()),
68  m_Smat_u (m_u_space.zeroVector()),
69  m_Bmat_with_permut (NULL), // to be deleted on destructor
70  m_Bmat_without_permut (NULL), // to be deleted on destructor
71  m_Bmat_rank (0),
72  m_Bwp_t__Wy__Bwp (NULL), // to be deleted on destructor
73  m_Bop_t__Wy__Bop (NULL), // to be deleted on destructor
74  m_Bwp_t__Wy__Bwp__inv (NULL), // to be deleted on destructor
75  m_Bop_t__Wy__Bop__inv (NULL), // to be deleted on destructor
76  m_a_y_modifier (0.),
77  m_b_y_modifier (0.)
78 {
79  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
80  *m_env.subDisplayFile() << "Entering GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
81  << ": allOutputsAreScalar = " << allOutputsAreScalar
82  << ": key-debug"
83  << ", some entities just created (not yet populated)"
84  << ", m_Smat_uw.numRowsLocal() = " << m_Smat_uw.numRowsLocal()
85  << ", m_Smat_uw.numCols() = " << m_Smat_uw.numCols()
86  << ", m_Smat_uw_t.numRowsLocal() = " << m_Smat_uw_t.numRowsLocal()
87  << ", m_Smat_uw_t.numCols() = " << m_Smat_uw_t.numCols()
88  << ", m_Zvec_hat_vu.sizeLocal() = " << m_Zvec_hat_vu.sizeLocal()
89  << ", m_Smat_u.numRowsLocal() = " << m_Smat_u.numRowsLocal()
90  << ", m_Smat_u.numCols() = " << m_Smat_u.numCols()
91  << std::endl;
92  }
93 
94  std::set<unsigned int> tmpSet;
95  tmpSet.insert(m_env.subId());
96 
97  if (allOutputsAreScalar) {
98  // Do nothing
99  // ppp: set m_Zvec_hat_vu?
100  }
101  else {
102  //********************************************************************************
103  // Form 'P_K' matrix
104  //********************************************************************************
105  D_M PK(m_u_space.zeroVector());
106  for (unsigned int i = 0; i < s.m_paper_p_eta; ++i) {
107  for (unsigned int j = 0; j < e.m_paper_n; ++j) {
108  unsigned int row = j + (e.m_paper_n*i);
109  unsigned int col = (j*s.m_paper_p_eta)+i;
110  PK(row,col) = 1.;
111  }
112  }
113 
114  if (m_env.checkingLevel() >= 1) {
115  // Check transpose operation
116  D_M PKt(PK.transpose());
117  if (m_env.subDisplayFile()) {
118  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
119  << ", tests on PK"
120  << ": PK.numRowsLocal() = " << PK.numRowsLocal()
121  << ", PK.numCols() = " << PK.numCols()
122  << ": PKt.numRowsLocal() = " << PKt.numRowsLocal()
123  << ", PKt.numCols() = " << PKt.numCols()
124  << std::endl;
125  }
126 
127  D_M matShouldBeI1( PK * PKt );
128  D_M matI1 (m_u_space.zeroVector());
129  for (unsigned int i = 0; i < matI1.numRowsLocal(); ++i) {
130  matI1(i,i) = 1.;
131  }
132  matShouldBeI1 -= matI1;
133  double auxNorm1 = matShouldBeI1.normFrob();
134 
135  D_M matShouldBeI2( PKt * PK );
136  D_M matI2 (m_u_space.zeroVector());
137  for (unsigned int i = 0; i < matI2.numRowsLocal(); ++i) {
138  matI2(i,i) = 1.;
139  }
140  matShouldBeI2 -= matI2;
141  double auxNorm2 = matShouldBeI2.normFrob();
142 
143  if (m_env.subDisplayFile()) {
144  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
145  << ", tests on PK"
146  << ": matShouldBeI1.numRowsLocal() = " << matShouldBeI1.numRowsLocal()
147  << ", ||matI1||_2^2 = " << matI1.normFrob() * matI1.normFrob()
148  << ", ||matShouldBeI1 - matI1||_2^2 = " << auxNorm1 * auxNorm1
149  << "; matShouldBeI2.numRowsLocal() = " << matShouldBeI2.numRowsLocal()
150  << ", ||matI2||_2^2 = " << matI2.normFrob() * matI2.normFrob()
151  << ", ||matShouldBeI2 - matI2||_2^2 = " << auxNorm2 * auxNorm2
152  << std::endl;
153  }
154  }
155 
156  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
157  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
158  << ": finished forming 'P_K'"
159  << std::endl;
160  }
161 
162  //********************************************************************************
163  // Form 'Kmat_interp_BlockDiag' matrix
164  //********************************************************************************
165  D_M Kmat_interp_BlockDiag (m_env,e.m_y_space.map(),m_u_size); // Formed with 'experimentModel.Kmats_interp()' // rr0: use D_M& and copy from 'experimentModel' object
166  D_M Kmat_interp_BlockDiag_permut(m_env,e.m_y_space.map(),m_u_size);
167 
168  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
169  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
170  << ": Kmat_interp_BlockDiag.numRowsLocal() = " << Kmat_interp_BlockDiag.numRowsLocal()
171  << ", Kmat_interp_BlockDiag.numCols() = " << Kmat_interp_BlockDiag.numCols()
172  << std::endl;
173  }
174 
175  Kmat_interp_BlockDiag.fillWithBlocksDiagonally(0,0,e.m_experimentModel.Kmats_interp(),true,true);
176 
177  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
178  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
179  << ": finished forming 'Kmat_interp_BlockDiag'"
180  << std::endl;
181  }
182 
183  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
184  Kmat_interp_BlockDiag.subWriteContents("Kmat_interp_BlockDiag",
185  "mat_Kmat_interp_BlockDiag",
186  "m",
187  tmpSet);
188  }
189 
190  //********************************************************************************
191  // Compute 'Kmat_interp_BlockDiag_permut' matrix
192  //********************************************************************************
193  Kmat_interp_BlockDiag_permut = Kmat_interp_BlockDiag * (PK.transpose());
194 
195  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
196  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
197  << ": finished computing 'Kmat_interp_BlockDiag_permut'"
198  << std::endl;
199  }
200 
201  //********************************************************************************
202  // Form 'B' matrix
203  //********************************************************************************
204  m_Bmat_with_permut = new D_M(m_env,e.m_y_space.map(),m_vu_size); // to be deleted on destructor
205  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
206  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
207  << ": key-debug"
208  << ", m_Bmat_with_permut just created (not yet populated)"
209  << ", numRowsLocal() = " << m_Bmat_with_permut->numRowsLocal()
210  << ", numCols() = " << m_Bmat_with_permut->numCols()
211  << std::endl;
212  }
213  m_Bmat_without_permut = new D_M(m_env,e.m_y_space.map(),m_vu_size); // to be deleted on destructor
214  m_Bmat_rank = std::min(m_Bmat_with_permut->numRowsGlobal(),m_Bmat_with_permut->numCols()); // Might be smaller
215  m_Bwp_t__Wy__Bwp = new D_M(m_vu_space.zeroVector()); // to be deleted on destructor
216  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
217  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
218  << ": key-debug"
219  << ", m_Bwp_t__Wy__Bwp just created (not yet populated)"
220  << ", numRowsLocal() = " << m_Bwp_t__Wy__Bwp->numRowsLocal()
221  << ", numCols() = " << m_Bwp_t__Wy__Bwp->numCols()
222  << std::endl;
223  }
224  m_Bop_t__Wy__Bop = new D_M(m_vu_space.zeroVector()); // to be deleted on destructor
225  m_Bwp_t__Wy__Bwp__inv = new D_M(m_vu_space.zeroVector()); // to be deleted on destructor
226  m_Bop_t__Wy__Bop__inv = new D_M(m_vu_space.zeroVector()); // to be deleted on destructor
227 
228  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
229  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor():"
230  << "\n m_Bmat_with_permut->numRowsLocal() = " << m_Bmat_with_permut->numRowsLocal()
231  << "\n m_Bmat_with_permut->numCols() = " << m_Bmat_with_permut->numCols()
232  << "\n m_Dmat_BlockDiag_permut->numRowsLocal() = " << e.m_Dmat_BlockDiag_permut->numRowsLocal()
233  << "\n m_Dmat_BlockDiag_permut->numCols() = " << e.m_Dmat_BlockDiag_permut->numCols()
234  << "\n m_PD->numRowsLocal() = " << e.m_PD->numRowsLocal()
235  << "\n m_PD->numCols() = " << e.m_PD->numCols()
236  << "\n Kmat_interp_BlockDiag_permut.numRowsLocal() = " << Kmat_interp_BlockDiag_permut.numRowsLocal()
237  << "\n Kmat_interp_BlockDiag_permut.numCols() = " << Kmat_interp_BlockDiag_permut.numCols()
238  << "\n PK.numRowsLocal() = " << PK.numRowsLocal()
239  << "\n PK.numCols() = " << PK.numCols()
240  << "\n gcmOptionsObj.m_ov.m_nuggetValueForBtWyB = " << gcmOptionsObj.m_ov.m_nuggetValueForBtWyB
241  << "\n gcmOptionsObj.m_ov.m_nuggetValueForBtWyBInv = " << gcmOptionsObj.m_ov.m_nuggetValueForBtWyBInv
242  << std::endl;
243  }
244 
245  std::vector<const P_M* > twoMats(2, (P_M*) NULL);
246 
247  twoMats[0] = e.m_Dmat_BlockDiag_permut;
248  twoMats[1] = &Kmat_interp_BlockDiag_permut;
249  m_Bmat_with_permut->fillWithBlocksHorizontally(0,0,twoMats,true,true);
250 
251  twoMats[0] = e.m_Dmat_BlockDiag;
252  twoMats[1] = &Kmat_interp_BlockDiag;
253  m_Bmat_without_permut->fillWithBlocksHorizontally(0,0,twoMats,true,true);
254 
255  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
256  m_Bmat_without_permut->subWriteContents("B_op",
257  "mat_B_op",
258  "m",
259  tmpSet);
260  }
261 
262  //********************************************************************************
263  // Form 'Bwp^T' matrix
264  //********************************************************************************
265  D_M Bwp_t(m_env,m_vu_space.map(),e.m_paper_n_y);
266  Bwp_t.fillWithTranspose(0,0,*m_Bmat_with_permut,true,true);
267 
268  if (m_env.checkingLevel() >= 1) {
269  // Check transpose operation
270  D_M Btt(m_env,e.m_y_space.map(),m_vu_size);
271  Btt.fillWithTranspose(0,0,Bwp_t,true,true);
272  Btt -= *m_Bmat_with_permut;
273  double btDiffNorm = Btt.normFrob();
274  if (m_env.subDisplayFile()) {
275  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
276  << ": ||Btt - B||_2 = " << btDiffNorm
277  << std::endl;
278  }
279  }
280 
281  double bRank14 = 0.;
282  if (m_Bmat_with_permut->numRowsGlobal() >= m_Bmat_with_permut->numCols()) {
283  m_Bmat_rank = m_Bmat_with_permut->rank(0.,1.e-8 ); // todo: should be an option
284  bRank14 = m_Bmat_with_permut->rank(0.,1.e-14);
285  }
286  else {
287  m_Bmat_rank = Bwp_t.rank(0.,1.e-8); // todo: should be an option
288  bRank14 = Bwp_t.rank(0.,1.e-14);
289  }
290 
291  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
292  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
293  << ": finished forming 'm_Bmat'"
294  << ", m_Bmat_with_permut->numRowsLocal() = " << m_Bmat_with_permut->numRowsLocal()
295  << ", m_Bmat_with_permut->numCols() = " << m_Bmat_with_permut->numCols()
296  << ", m_Bmat_with_permut->rank(0.,1.e-8) = " << m_Bmat_rank
297  << ", m_Bmat_with_permut->rank(0.,1.e-14) = " << bRank14
298  << std::endl;
299  }
300  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
301  m_Bmat_with_permut->subWriteContents("B_wp",
302  "mat_B_wp",
303  "m",
304  tmpSet);
305  }
306  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
307  D_M Bmat_filtered(*m_Bmat_with_permut);
308  Bmat_filtered.setPrintHorizontally(false);
309  Bmat_filtered.filterSmallValues(1.e-6);
310  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
311  << ": Bmat_filtered.numRowsLocal() = " << Bmat_filtered.numRowsLocal()
312  << ", Bmat_filtered.numCols() = " << Bmat_filtered.numCols()
313  << ", Bmat_filtered contents =\n" << Bmat_filtered
314  << std::endl;
315  }
316 
317  //********************************************************************************
318  // Compute 'Bwp^T W_y Bwp' matrix
319  //********************************************************************************
320  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
321  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor():"
322  << "\n m_Bwp_t__Wy__Bwp->numRowsLocal() = " << m_Bwp_t__Wy__Bwp->numRowsLocal()
323  << "\n m_Bwp_t__Wy__Bwp->numCols() = " << m_Bwp_t__Wy__Bwp->numCols()
324  << "\n Bwp_t.numRowsLocal() = " << Bwp_t.numRowsLocal()
325  << "\n Bwp_t.numCols() = " << Bwp_t.numCols()
326  << "\n m_Wy->numRowsLocal() = " << e.m_Wy->numRowsLocal()
327  << "\n m_Wy->numCols() = " << e.m_Wy->numCols()
328  << "\n m_Bmat_with_permut->numRowsLocal() = " << m_Bmat_with_permut->numRowsLocal()
329  << "\n m_Bmat_with_permut->numCols() = " << m_Bmat_with_permut->numCols()
330  << std::endl;
331  }
332  *m_Bwp_t__Wy__Bwp = Bwp_t * (*e.m_Wy * *m_Bmat_with_permut);
333 
334  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
335  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
336  << ": finished computing 'm_Bwp_t__Wy__Bwp'"
337  << std::endl;
338  }
339  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
340  m_Bwp_t__Wy__Bwp->subWriteContents("Bwp_t__Wy__Bwp",
341  "mat_Bwp_t__Wy__Bwp",
342  "m",
343  tmpSet);
344  }
345 
346  if (m_env.displayVerbosity() >= 4) {
347  double Bwp_t__Wy__Bwp__LnDeterminant = m_Bwp_t__Wy__Bwp->lnDeterminant();
348  unsigned int Bwp_t__Wy__Bwp__Rank = m_Bwp_t__Wy__Bwp->rank(0.,1.e-8 ); // todo: should be an option
349  unsigned int Bwp_t__Wy__Bwp__Rank14 = m_Bwp_t__Wy__Bwp->rank(0.,1.e-14);
350  if (m_env.subDisplayFile()) {
351  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
352  << ": m_Bwp_t__Wy__Bwp->numRowsLocal() = " << m_Bwp_t__Wy__Bwp->numRowsLocal()
353  << ", m_Bwp_t__Wy__Bwp->numCols() = " << m_Bwp_t__Wy__Bwp->numCols()
354  << ", m_Bwp_t__Wy__Bwp->lnDeterminant() = " << Bwp_t__Wy__Bwp__LnDeterminant
355  << ", m_Bwp_t__Wy__Bwp->rank(0.,1.e-8) = " << Bwp_t__Wy__Bwp__Rank
356  << ", m_Bwp_t__Wy__Bwp->rank(0.,1.e-14) = " << Bwp_t__Wy__Bwp__Rank14
357  << std::endl;
358  }
359  }
360 
361  //********************************************************************************
362  // Compute 'Bwp^T W_y Bwp' inverse
363  //********************************************************************************
364  *m_Bwp_t__Wy__Bwp__inv = m_Bwp_t__Wy__Bwp->inverse(); // inversion savings
365  if (m_env.displayVerbosity() >= 4) {
366  double Bwp_t__Wy__Bwp__inv__LnDeterminant = m_Bwp_t__Wy__Bwp__inv->lnDeterminant();
367  if (m_env.subDisplayFile()) {
368  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
369  << ": finished computing 'm_Bwp_t__Wy__Bwp__inv'"
370  << ", m_Bwp_t__Wy__Bwp__inv->lnDeterminant() = " << Bwp_t__Wy__Bwp__inv__LnDeterminant
371  << std::endl;
372  }
373  }
374  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
375  m_Bwp_t__Wy__Bwp__inv->subWriteContents("Bwp_t__Wy__Bwp__inv",
376  "mat_Bwp_t__Wy__Bwp__inv",
377  "m",
378  tmpSet);
379  }
380 
381  if (m_env.displayVerbosity() >= 4) {
382  double Bwp_t__Wy__Bwp__InvLnDeterminant = m_Bwp_t__Wy__Bwp__inv->lnDeterminant();
383  unsigned int Bwp_t__Wy__Bwp__InvRank = m_Bwp_t__Wy__Bwp__inv->rank(0.,1.e-8 ); // todo: should be an option
384  unsigned int Bwp_t__Wy__Bwp__InvRank14 = m_Bwp_t__Wy__Bwp__inv->rank(0.,1.e-14);
385  if (m_env.subDisplayFile()) {
386  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
387  << ", m_Bwp_t__Wy__Bwp__inv->numRowsLocal() = " << m_Bwp_t__Wy__Bwp__inv->numRowsLocal()
388  << ", m_Bwp_t__Wy__Bwp__inv->numCols() = " << m_Bwp_t__Wy__Bwp__inv->numCols()
389  << ": m_Bwp_t__Wy__Bwp__inv->lnDeterminant() = " << Bwp_t__Wy__Bwp__InvLnDeterminant
390  << ": m_Bwp_t__Wy__Bwp__inv->rank(0.,1.e-8) = " << Bwp_t__Wy__Bwp__InvRank
391  << ": m_Bwp_t__Wy__Bwp__inv->rank(0.,1.e-14) = " << Bwp_t__Wy__Bwp__InvRank14
392  << std::endl;
393  }
394  }
395 
396  if (m_env.checkingLevel() >= 1) {
397  // Check transpose operation
398  D_M matShouldBeI1(*m_Bwp_t__Wy__Bwp * *m_Bwp_t__Wy__Bwp__inv);
399  D_M matShouldBeI2(*m_Bwp_t__Wy__Bwp__inv * *m_Bwp_t__Wy__Bwp );
400  D_M matI (m_vu_space.zeroVector());
401  for (unsigned int i = 0; i < matI.numRowsLocal(); ++i) {
402  matI(i,i) = 1.;
403  }
404  matShouldBeI1 -= matI;
405  double auxNorm1 = matShouldBeI1.normFrob();
406  matShouldBeI2 -= matI;
407  double auxNorm2 = matShouldBeI2.normFrob();
408  if (m_env.subDisplayFile()) {
409  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
410  << ", tests on m_Bt_Wy_B"
411  << ": matShouldBeI1.numRowsLocal() = " << matShouldBeI1.numRowsLocal()
412  //<< ", ||matShouldBeI1||_2^2 = " << matShouldBeI1.normFrob() * matShouldBeI1.normFrob()
413  //<< ", ||matShouldBeI2||_2^2 = " << matShouldBeI2.normFrob() * matShouldBeI2.normFrob()
414  << ", ||matI||_2^2 = " << matI.normFrob() * matI.normFrob()
415  << ", ||matShouldBeI1 - matI||_2^2 = " << auxNorm1 * auxNorm1
416  << ", ||matShouldBeI2 - matI||_2^2 = " << auxNorm2 * auxNorm2
417  << std::endl;
418  }
419  }
420 
421  //********************************************************************************
422  // Form 'Bop^T' matrix
423  //********************************************************************************
424  D_M Bop_t(m_env,m_vu_space.map(),e.m_paper_n_y);
425  Bop_t.fillWithTranspose(0,0,*m_Bmat_without_permut,true,true);
426 
427  //********************************************************************************
428  // Compute 'Bop^T W_y Bop' matrix
429  //********************************************************************************
430  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
431  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor():"
432  << "\n m_Bop_t__Wy__Bop->numRowsLocal() = " << m_Bop_t__Wy__Bop->numRowsLocal()
433  << "\n m_Bop_t__Wy__Bop->numCols() = " << m_Bop_t__Wy__Bop->numCols()
434  << "\n Bop_t.numRowsLocal() = " << Bop_t.numRowsLocal()
435  << "\n Bop_t.numCols() = " << Bop_t.numCols()
436  << "\n m_Wy->numRowsLocal() = " << e.m_Wy->numRowsLocal()
437  << "\n m_Wy->numCols() = " << e.m_Wy->numCols()
438  << "\n m_Bmat_with_permut->numRowsLocal() = " << m_Bmat_with_permut->numRowsLocal()
439  << "\n m_Bmat_with_permut->numCols() = " << m_Bmat_with_permut->numCols()
440  << std::endl;
441  }
442 
443  *m_Bop_t__Wy__Bop = Bop_t * (*e.m_Wy * *m_Bmat_without_permut);
444  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
445  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
446  << ": finished computing 'm_Bop_t__Wy__Bop before nugget'"
447  << std::endl;
448  }
449  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
450  m_Bop_t__Wy__Bop->subWriteContents("Bop_t__Wy__Bop__beforePermut_beforeNugget",
451  "mat_Bop_t__Wy__Bop__beforePermut_beforeNugget",
452  "m",
453  tmpSet);
454  }
455 
456  // Just for checking
457  if (m_env.displayVerbosity() >= 4) {
458  double Bop_t__Wy__Bop__beforeNugget__LnDeterminant = m_Bop_t__Wy__Bop->lnDeterminant();
459  unsigned int Bop_t__Wy__Bop__beforeNugget__Rank = m_Bop_t__Wy__Bop->rank(0.,1.e-8 ); // todo: should be an option
460  unsigned int Bop_t__Wy__Bop__beforeNugget__Rank14 = m_Bop_t__Wy__Bop->rank(0.,1.e-14);
461  if (m_env.subDisplayFile()) {
462  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
463  << ", m_Bop_t__Wy__Bop__beforeNugget.numRowsLocal() = " << m_Bop_t__Wy__Bop->numRowsLocal()
464  << ", m_Bop_t__Wy__Bop__beforeNugget.numCols() = " << m_Bop_t__Wy__Bop->numCols()
465  << ": m_Bop_t__Wy__Bop__beforeNugget.lnDeterminant() = " << Bop_t__Wy__Bop__beforeNugget__LnDeterminant
466  << ", m_Bop_t__Wy__Bop__beforeNugget.rank(0.,1.e-8) = " << Bop_t__Wy__Bop__beforeNugget__Rank
467  << ", m_Bop_t__Wy__Bop__beforeNugget.rank(0.,1.e-14) = " << Bop_t__Wy__Bop__beforeNugget__Rank14
468  << std::endl;
469  }
470  }
471 
472  //********************************************************************************
473  // Apply permutation only now
474  //********************************************************************************
475  twoMats[0] = e.m_PD;
476  twoMats[1] = &PK;
477  D_M matP(m_vu_space.zeroVector());
478  matP.fillWithBlocksDiagonally(0,0,twoMats,true,true);
479  *m_Bop_t__Wy__Bop = matP * (*m_Bop_t__Wy__Bop * (matP.transpose()));
480 
481  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
482  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
483  << ": finished computing 'm_Bop_t__Wy__Bop after permutation'"
484  << std::endl;
485  }
486  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
487  m_Bop_t__Wy__Bop->subWriteContents("Bop_t__Wy__Bop__afterPermut_beforeNugget",
488  "mat_Bop_t__Wy__Bop__afterPermut_beforeNugget",
489  "m",
490  tmpSet);
491  }
492 
493  if (m_env.displayVerbosity() >= 4) {
494  double Bop_t__Wy__Bop__LnDeterminant = m_Bop_t__Wy__Bop->lnDeterminant();
495  unsigned int Bop_t__Wy__Bop__Rank = m_Bop_t__Wy__Bop->rank(0.,1.e-8 ); // todo: should be an option
496  unsigned int Bop_t__Wy__Bop__Rank14 = m_Bop_t__Wy__Bop->rank(0.,1.e-14);
497  if (m_env.subDisplayFile()) {
498  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
499  << ": m_Bop_t__Wy__Bop->numRowsLocal() = " << m_Bop_t__Wy__Bop->numRowsLocal()
500  << ", m_Bop_t__Wy__Bop->numCols() = " << m_Bop_t__Wy__Bop->numCols()
501  << ", m_Bop_t__Wy__Bop->lnDeterminant() = " << Bop_t__Wy__Bop__LnDeterminant
502  << ", m_Bop_t__Wy__Bop->rank(0.,1.e-8) = " << Bop_t__Wy__Bop__Rank
503  << ", m_Bop_t__Wy__Bop->rank(0.,1.e-14) = " << Bop_t__Wy__Bop__Rank14
504  << std::endl;
505  }
506  }
507 
508  //********************************************************************************
509  // Add nugget only now
510  //********************************************************************************
511  if (gcmOptionsObj.m_ov.m_nuggetValueForBtWyB != 0.) {
512  for (unsigned int i = 0; i < m_Bop_t__Wy__Bop->numRowsLocal(); ++i) {
513  (*m_Bop_t__Wy__Bop)(i,i) += gcmOptionsObj.m_ov.m_nuggetValueForBtWyB;
514  }
515  }
516  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
517  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
518  << ": finished computing 'm_Bop_t__Wy__Bop after permutation and after nugget'"
519  << std::endl;
520  }
521  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
522  m_Bop_t__Wy__Bop->subWriteContents("Bop_t__Wy__Bop__afterPermut_afterNugget",
523  "mat_Bop_t__Wy__Bop__afterPermut_afterNugget",
524  "m",
525  tmpSet);
526  }
527 
528  // Just for checking
529  if (m_env.displayVerbosity() >= 4) {
530  double Bop_t__Wy__Bop__afterNugget__LnDeterminant = m_Bop_t__Wy__Bop->lnDeterminant();
531  unsigned int Bop_t__Wy__Bop__afterNugget__Rank = m_Bop_t__Wy__Bop->rank(0.,1.e-8 ); // todo: should be an option
532  unsigned int Bop_t__Wy__Bop__afterNugget__Rank14 = m_Bop_t__Wy__Bop->rank(0.,1.e-14);
533  if (m_env.subDisplayFile()) {
534  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
535  << ", m_Bop_t__Wy__Bop__afterNugget.numRowsLocal() = " << m_Bop_t__Wy__Bop->numRowsLocal()
536  << ", m_Bop_t__Wy__Bop__afterNugget.numCols() = " << m_Bop_t__Wy__Bop->numCols()
537  << ": m_Bop_t__Wy__Bop__afterNugget.lnDeterminant() = " << Bop_t__Wy__Bop__afterNugget__LnDeterminant
538  << ", m_Bop_t__Wy__Bop__afterNugget.rank(0.,1.e-8) = " << Bop_t__Wy__Bop__afterNugget__Rank
539  << ", m_Bop_t__Wy__Bop__afterNugget.rank(0.,1.e-14) = " << Bop_t__Wy__Bop__afterNugget__Rank14
540  << std::endl;
541  }
542  }
543 
544  //********************************************************************************
545  // Compute 'Bop^T W_y Bop' inverse
546  //********************************************************************************
547  *m_Bop_t__Wy__Bop__inv = m_Bop_t__Wy__Bop->inverse(); // inversion savings
548 
549  if (m_env.displayVerbosity() >= 4) {
550  double Bop_t__Wy__Bop__inv__beforeNugget__LnDeterminant = m_Bop_t__Wy__Bop__inv->lnDeterminant();
551  unsigned int Bop_t__Wy__Bop__beforeNugget__InvRank = m_Bop_t__Wy__Bop__inv->rank(0.,1.e-8 ); // todo: should be an option
552  unsigned int Bop_t__Wy__Bop__beforeNugget__InvRank14 = m_Bop_t__Wy__Bop__inv->rank(0.,1.e-14);
553  if (m_env.subDisplayFile()) {
554  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
555  << ", m_Bop_t__Wy__Bop__inv__beforeNugget.numRowsLocal() = " << m_Bop_t__Wy__Bop__inv->numRowsLocal()
556  << ", m_Bop_t__Wy__Bop__inv__beforeNugget.numCols() = " << m_Bop_t__Wy__Bop__inv->numCols()
557  << ": m_Bop_t__Wy__Bop__inv__beforeNugget.lnDeterminant() = " << Bop_t__Wy__Bop__inv__beforeNugget__LnDeterminant
558  << ", m_Bop_t__Wy__Bop__inv__beforeNugget.rank(0.,1.e-8) = " << Bop_t__Wy__Bop__beforeNugget__InvRank
559  << ", m_Bop_t__Wy__Bop__inv__beforeNugget.rank(0.,1.e-14) = " << Bop_t__Wy__Bop__beforeNugget__InvRank14
560  << std::endl;
561  }
562  }
563 
564  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
565  m_Bop_t__Wy__Bop__inv->subWriteContents("Bop_t__Wy__Bop__inv__beforeNuggetForInv",
566  "mat_Bop_t__Wy__Bop__inv__beforeNuggetForInv",
567  "m",
568  tmpSet);
569  }
570 
571  // Add nugget
572  if (gcmOptionsObj.m_ov.m_nuggetValueForBtWyBInv != 0.) {
573  for (unsigned int i = 0; i < m_Bop_t__Wy__Bop__inv->numRowsLocal(); ++i) {
574  (*m_Bop_t__Wy__Bop__inv)(i,i) += gcmOptionsObj.m_ov.m_nuggetValueForBtWyBInv;
575  }
576  }
577 
578  if (m_env.displayVerbosity() >= 4) {
579  double Bop_t__Wy__Bop__inv__afterNugget__LnDeterminant = m_Bop_t__Wy__Bop__inv->lnDeterminant();
580  unsigned int Bop_t__Wy__Bop__afterNugget__InvRank = m_Bop_t__Wy__Bop__inv->rank(0.,1.e-8 ); // todo: should be an option
581  unsigned int Bop_t__Wy__Bop__afterNugget__InvRank14 = m_Bop_t__Wy__Bop__inv->rank(0.,1.e-14);
582  if (m_env.subDisplayFile()) {
583  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
584  << ", m_Bop_t__Wy__Bop__inv__afterNugget.numRowsLocal() = " << m_Bop_t__Wy__Bop__inv->numRowsLocal()
585  << ", m_Bop_t__Wy__Bop__inv__afterNugget.numCols() = " << m_Bop_t__Wy__Bop__inv->numCols()
586  << ": m_Bop_t__Wy__Bop__inv__afterNugget.lnDeterminant() = " << Bop_t__Wy__Bop__inv__afterNugget__LnDeterminant
587  << ", m_Bop_t__Wy__Bop__inv__afterNugget.rank(0.,1.e-8) = " << Bop_t__Wy__Bop__afterNugget__InvRank
588  << ", m_Bop_t__Wy__Bop__inv__afterNugget.rank(0.,1.e-14) = " << Bop_t__Wy__Bop__afterNugget__InvRank14
589  << std::endl;
590  }
591  }
592 
593  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
594  m_Bop_t__Wy__Bop__inv->subWriteContents("Bop_t__Wy__Bop__inv__afterNuggetForInv",
595  "mat_Bop_t__Wy__Bop__inv__afterNuggetForInv",
596  "m",
597  tmpSet);
598  }
599 
600  if (m_env.checkingLevel() >= 1) {
601  // Check transpose operation
602  D_M matShouldBeI1(*m_Bop_t__Wy__Bop * *m_Bop_t__Wy__Bop__inv);
603  D_M matShouldBeI2(*m_Bop_t__Wy__Bop__inv * *m_Bop_t__Wy__Bop );
604  D_M matI (m_vu_space.zeroVector());
605  for (unsigned int i = 0; i < matI.numRowsLocal(); ++i) {
606  matI(i,i) = 1.;
607  }
608  matShouldBeI1 -= matI;
609  double auxNorm1 = matShouldBeI1.normFrob();
610  matShouldBeI2 -= matI;
611  double auxNorm2 = matShouldBeI2.normFrob();
612  if (m_env.subDisplayFile()) {
613  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
614  << ", tests on m_Bt_Wy_B"
615  << ": matShouldBeI1.numRowsLocal() = " << matShouldBeI1.numRowsLocal()
616  //<< ", ||matShouldBeI1||_2^2 = " << matShouldBeI1.normFrob() * matShouldBeI1.normFrob()
617  //<< ", ||matShouldBeI2||_2^2 = " << matShouldBeI2.normFrob() * matShouldBeI2.normFrob()
618  << ", ||matI||_2^2 = " << matI.normFrob() * matI.normFrob()
619  << ", ||matShouldBeI1 - matI||_2^2 = " << auxNorm1 * auxNorm1
620  << ", ||matShouldBeI2 - matI||_2^2 = " << auxNorm2 * auxNorm2
621  << std::endl;
622  }
623  }
624 
625  //********************************************************************************
626  // Compute exponent modifiers
627  //********************************************************************************
628  m_a_y_modifier = ((double) (e.m_paper_n_y - m_Bmat_rank)) / 2.;
629  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
630  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
631  << ": m_a_y_modifier = " << m_a_y_modifier
632  << std::endl;
633  }
634 
635  D_V yVec_transformed(e.m_experimentStorage.yVec_transformed());
636  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
637  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
638  << ": m_Zvec_hat_vu.sizeLocal() = " << m_Zvec_hat_vu.sizeLocal()
639  << ", m_Bop_t__Wy__Bop__inv->numRowsLocal() = " << m_Bop_t__Wy__Bop__inv->numRowsLocal()
640  << ", m_Bop_t__Wy__Bop__inv->numCols() = " << m_Bop_t__Wy__Bop__inv->numCols()
641  << ", Bop_t.numRowsLocal() = " << Bop_t.numRowsLocal()
642  << ", Bop_t.numCols() = " << Bop_t.numCols()
643  << ", m_Wy->numRowsLocal() = " << e.m_Wy->numRowsLocal()
644  << ", m_Wy->numCols() = " << e.m_Wy->numCols()
645  << ", yVec_transformed.sizeLocal() = " << yVec_transformed.sizeLocal()
646  << std::endl;
647  }
648  m_Zvec_hat_vu = *m_Bop_t__Wy__Bop__inv * (Bwp_t * (*e.m_Wy * yVec_transformed)); // aqui
649  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
650  m_Zvec_hat_vu.subWriteContents("Zvec_hat_vu",
651  "vec_Zvec_hat_vu",
652  "m",
653  tmpSet);
654  }
655  D_V tmpVec2(yVec_transformed - (*m_Bmat_with_permut * m_Zvec_hat_vu));
656  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
657  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
658  << ", yVec_transformed = " << yVec_transformed.sizeLocal()
659  << ", m_Zvec_hat_vu = " << m_Zvec_hat_vu
660  << ", tmpVec2 (initial) = " << tmpVec2
661  << std::endl;
662  }
663 
664  tmpVec2 = *e.m_Wy * tmpVec2;
665  m_b_y_modifier = scalarProduct(yVec_transformed,tmpVec2) / 2.;
666  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
667  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
668  << ": tmpVec2 (final) = " << tmpVec2
669  << ", m_b_y_modifier = " << m_b_y_modifier
670  << std::endl;
671  }
672  } // if (allOutputsAreScalar)
673 
674  //********************************************************************************
675  // Instantiate Smat spaces
676  //********************************************************************************
677  unsigned int sumDims = 0;
678  for (unsigned int i = 0; i < m_Smat_u_is.size(); ++i) {
679  m_Rmat_u_is[i] = new D_M(e.m_paper_n_space.zeroVector()); // to be deleted on destructor
680  m_Smat_u_is[i] = new D_M(e.m_paper_n_space.zeroVector()); // to be deleted on destructor
681  sumDims += e.m_paper_n;
682  }
683  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
684  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
685  << ": finished instantiating the m_Smat_u_i spaces"
686  << std::endl;
687  }
688  queso_require_equal_to_msg(sumDims, m_u_size, "'sumDims' and 'm_u_size' should be equal");
689 
690  unsigned int sumNumRows = 0;
691  unsigned int sumNumCols = 0;
692  for (unsigned int i = 0; i < m_Smat_uw_is.size(); ++i) {
693  m_Rmat_uw_is[i] = new D_M(m_env,e.m_paper_n_space.map(),s.m_paper_m); // Yes, 'u' only // to be deleted on destructor
694  m_Smat_uw_is[i] = new D_M(m_env,e.m_paper_n_space.map(),s.m_paper_m); // Yes, 'u' only // to be deleted on destructor
695  sumNumRows += e.m_paper_n_space.dimLocal();
696  sumNumCols += s.m_paper_m;
697  }
698  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
699  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
700  << ": finished instantiating the m_Smat_uw_i matrices"
701  << std::endl;
702  }
703  queso_require_equal_to_msg(sumNumRows, m_u_size, "'sumNumRows' and 'm_u_size' should be equal");
704  queso_require_equal_to_msg(sumNumCols, s.m_w_size, "'sumNumCols' and 'm_w_size' should be equal");
705 
706  //********************************************************************************
707  // Instantiate 'u_hat_u_asterisk' matrices
708  //********************************************************************************
709  sumNumRows = 0;
710  sumNumCols = 0;
711  for (unsigned int i = 0; i < m_Smat_u_hat_u_asterisk_is.size(); ++i) {
712  m_Rmat_u_hat_u_asterisk_is[i] = new D_M(m_env, e.m_paper_n_space.map(), (unsigned int) 1); // to be deleted on destructor; Yes, only 1 column
713  m_Smat_u_hat_u_asterisk_is[i] = new D_M(m_env, e.m_paper_n_space.map(), (unsigned int) 1); // to be deleted on destructor; Yes, only 1 column
714  sumNumRows += e.m_paper_n_space.dimLocal();
715  sumNumCols += 1;
716  }
717  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
718  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
719  << ": finished instantiating the m_Smat_u_hat_u_asterisk_i matrices"
720  << std::endl;
721  }
722  queso_require_equal_to_msg(sumNumRows, m_u_size, "'sumNumRows' and 'm_u_size' should be equal");
723  queso_require_equal_to_msg(sumNumCols, s.m_paper_p_eta, "'sumNumCols (1)' and 'm_paper_p_eta' should be equal");
724 
725  //********************************************************************************
726  // Instantiate 'w_hat_u_asterisk' matrices
727  //********************************************************************************
728  sumNumRows = 0;
729  sumNumCols = 0;
730  for (unsigned int i = 0; i < m_Smat_w_hat_u_asterisk_is.size(); ++i) {
731  m_Rmat_w_hat_u_asterisk_is[i] = new D_M(m_env, s.m_paper_m_space.map(), (unsigned int) 1); // to be deleted on destructor; Yes, only 1 column
732  m_Smat_w_hat_u_asterisk_is[i] = new D_M(m_env, s.m_paper_m_space.map(), (unsigned int) 1); // to be deleted on destructor; Yes, only 1 column
733  sumNumRows += s.m_paper_m_space.dimLocal();
734  sumNumCols += 1;
735  }
736  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
737  *m_env.subDisplayFile() << "In GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
738  << ": finished instantiating the m_Smat_w_hat_u_asterisk_i matrices"
739  << std::endl;
740  }
741  queso_require_equal_to_msg(sumNumRows, s.m_w_size, "'sumNumRows' and 'm_w_size' should be equal");
742  queso_require_equal_to_msg(sumNumCols, s.m_paper_p_eta, "'sumNumCols (2)' and 'm_paper_p_eta' should be equal");
743 
744  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
745  *m_env.subDisplayFile() << "Leaving GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
746  << std::endl;
747  }
748 }
749 
750 template <class S_V,class S_M,class D_V,class D_M,class P_V,class P_M,class Q_V,class Q_M>
751 GcmJointInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::~GcmJointInfo()
752 {
753  for (unsigned int i = 0; i < m_Smat_w_hat_u_asterisk_is.size(); ++i) {
754  delete m_Smat_w_hat_u_asterisk_is[i]; // to be deleted on destructor
755  m_Smat_w_hat_u_asterisk_is[i] = NULL;
756  delete m_Rmat_w_hat_u_asterisk_is[i]; // to be deleted on destructor
757  m_Rmat_w_hat_u_asterisk_is[i] = NULL;
758  }
759 
760  for (unsigned int i = 0; i < m_Smat_u_hat_u_asterisk_is.size(); ++i) {
761  delete m_Smat_u_hat_u_asterisk_is[i]; // to be deleted on destructor
762  m_Smat_u_hat_u_asterisk_is[i] = NULL;
763  delete m_Rmat_u_hat_u_asterisk_is[i]; // to be deleted on destructor
764  m_Rmat_u_hat_u_asterisk_is[i] = NULL;
765  }
766 
767  for (unsigned int i = 0; i < m_Smat_uw_is.size(); ++i) {
768  delete m_Smat_uw_is[i]; // to be deleted on destructor
769  m_Smat_uw_is[i] = NULL;
770  delete m_Rmat_uw_is[i]; // to be deleted on destructor
771  m_Rmat_uw_is[i] = NULL;
772  }
773 
774  for (unsigned int i = 0; i < m_Smat_u_is.size(); ++i) {
775  delete m_Smat_u_is[i]; // to be deleted on destructor
776  m_Smat_u_is[i] = NULL;
777  delete m_Rmat_u_is[i]; // to be deleted on destructor
778  m_Rmat_u_is[i] = NULL;
779  }
780 
781  delete m_Bop_t__Wy__Bop__inv;
782  delete m_Bwp_t__Wy__Bwp__inv;
783  delete m_Bop_t__Wy__Bop;
784  delete m_Bwp_t__Wy__Bwp;
785  delete m_Bmat_without_permut;
786  delete m_Bmat_with_permut;
787 }
788 
789 } // End namespace QUESO
790 
791 template class QUESO::GcmJointInfo<QUESO::GslVector, QUESO::GslMatrix, QUESO::GslVector, QUESO::GslMatrix, QUESO::GslVector, QUESO::GslMatrix, QUESO::GslVector, QUESO::GslMatrix>;
QUESO::BaseEnvironment::displayVerbosity
unsigned int displayVerbosity() const
Definition: Environment.C:400
QUESO::GcmJointInfo::GcmJointInfo
GcmJointInfo(const GpmsaComputerModelOptions &gcmOptionsObj, bool allOutputsAreScalar, const GcmSimulationInfo< S_V, S_M, P_V, P_M, Q_V, Q_M > &s, const GcmExperimentInfo< S_V, S_M, D_V, D_M, P_V, P_M > &e)
Definition: GcmJointInfo.C:32
QUESO::GcmJointInfo::m_vu_space
VectorSpace< D_V, D_M > m_vu_space
Definition: GcmJointInfo.h:72
QUESO::GcmJointInfo::m_Smat_w_hat_u_asterisk_is
std::vector< D_M * > m_Smat_w_hat_u_asterisk_is
Definition: GcmJointInfo.h:67
QUESO::GcmJointInfo::m_vu_size
unsigned int m_vu_size
Definition: GcmJointInfo.h:71
QUESO::GcmJointInfo::m_u_space
VectorSpace< D_V, D_M > m_u_space
Definition: GcmJointInfo.h:53
QUESO::VectorSpace::map
const Map & map() const
Map.
Definition: VectorSpace.C:157
QUESO::GcmJointInfo::m_Bmat_rank
unsigned int m_Bmat_rank
Definition: GcmJointInfo.h:111
QUESO::GcmJointInfo::m_env
const BaseEnvironment & m_env
Definition: GcmJointInfo.h:48
QUESO::GcmJointInfo::m_Smat_uw_is
std::vector< D_M * > m_Smat_uw_is
Definition: GcmJointInfo.h:57
QUESO::GcmJointInfo::m_u_size
unsigned int m_u_size
Definition: GcmJointInfo.h:52
queso_require_equal_to_msg
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:85
QUESO::GcmSimulationInfo::m_paper_m_space
VectorSpace< P_V, P_M > m_paper_m_space
Definition: GcmSimulationInfo.h:61
QUESO::GcmJointInfo::m_Rmat_u_is
std::vector< D_M * > m_Rmat_u_is
Definition: GcmJointInfo.h:54
QUESO::GcmExperimentInfo::m_experimentModel
const ExperimentModel< S_V, S_M, D_V, D_M > & m_experimentModel
Definition: GcmExperimentInfo.h:53
QUESO::GcmExperimentInfo::m_y_space
VectorSpace< D_V, D_M > m_y_space
Definition: GcmExperimentInfo.h:108
QUESO::GcmExperimentInfo::m_paper_n_space
VectorSpace< P_V, P_M > m_paper_n_space
Definition: GcmExperimentInfo.h:63
QUESO::GcmJointInfo::m_Rmat_u_hat_u_asterisk_is
std::vector< D_M * > m_Rmat_u_hat_u_asterisk_is
Definition: GcmJointInfo.h:61
QUESO::BaseEnvironment::subDisplayFile
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:278
QUESO::GcmJointInfo::m_Rmat_uw_is
std::vector< D_M * > m_Rmat_uw_is
Definition: GcmJointInfo.h:56
QUESO::BaseEnvironment::subId
unsigned int subId() const
Access function to the number of each sub-environment Id: m_subId.
Definition: Environment.C:299
QUESO::GcmOptionsValues::m_dataOutputAllowedSet
std::set< unsigned int > m_dataOutputAllowedSet
Definition: GpmsaComputerModelOptions.h:79
QUESO::GcmJointInfo::m_Smat_u_hat_u_asterisk_is
std::vector< D_M * > m_Smat_u_hat_u_asterisk_is
Definition: GcmJointInfo.h:62
QUESO::VectorSpace::zeroVector
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:189
QUESO::GcmExperimentInfo::m_experimentStorage
const ExperimentStorage< S_V, S_M, D_V, D_M > & m_experimentStorage
Definition: GcmExperimentInfo.h:52
QUESO::GcmJointInfo::m_Rmat_w_hat_u_asterisk_is
std::vector< D_M * > m_Rmat_w_hat_u_asterisk_is
Definition: GcmJointInfo.h:66
QUESO::MiscUintDebugMessage
unsigned int MiscUintDebugMessage(unsigned int value, const char *message)
Definition: Miscellaneous.C:314
QUESO::BaseEnvironment::checkingLevel
unsigned int checkingLevel() const
Access function to private attribute m_checkingLevel.
Definition: Environment.C:414
QUESO::VectorSpace::dimLocal
unsigned int dimLocal() const
Definition: VectorSpace.C:170
QUESO::GcmJointInfo::m_Smat_u_is
std::vector< D_M * > m_Smat_u_is
Definition: GcmJointInfo.h:55
QUESO::scalarProduct
double scalarProduct(const GslVector &x, const GslVector &y)
Definition: GslVector.C:1136
Generated on Fri Jun 17 2016 14:17:41 for queso-0.55.0 by  doxygen 1.8.5