@@ -612,7 +612,7 @@ namespace BasisClasses
612612 Eigen::VectorXd workE;
613613
614614 // M loop
615- for (int m=0 , moffset=0 , moffE= 0 ; m<=l; m++, L++) {
615+ for (int m=0 , moffset=0 ; m<=l; m++, L++) {
616616
617617 fac = factorial (l, m) * legs[tid](l, m);
618618
@@ -4055,7 +4055,7 @@ namespace BasisClasses
40554055 expcoef (ix, iy, iz) += - mass * curr (0 )*curr (1 )*curr (2 ) * norm;
40564056
40574057 if (pcavar)
4058- g[index1D (ix, iy, iz)] = - mass * curr (0 )*curr (1 )*curr (2 ) * norm;
4058+ g[index1D (ix, iy, iz)] = - curr (0 )*curr (1 )*curr (2 ) * norm;
40594059 }
40604060 }
40614061 }
@@ -4068,7 +4068,7 @@ namespace BasisClasses
40684068 sampleMasses (T) += mass;
40694069
40704070 meanV[T].noalias () += g * mass;
4071- covrV[T].noalias () += g * g.adjoint () * mass;
4071+ if (covar) covrV[T].noalias () += g * g.adjoint () * mass;
40724072 }
40734073 }
40744074
@@ -4099,8 +4099,9 @@ namespace BasisClasses
40994099 MPI_Allreduce (MPI_IN_PLACE , meanV[T].data (), meanV[T].size (),
41004100 MPI_DOUBLE_COMPLEX , MPI_SUM , MPI_COMM_WORLD );
41014101
4102- MPI_Allreduce (MPI_IN_PLACE , covrV[T].data (), covrV[T].size (),
4103- MPI_DOUBLE_COMPLEX , MPI_SUM , MPI_COMM_WORLD );
4102+ if (covar)
4103+ MPI_Allreduce (MPI_IN_PLACE , covrV[T].data (), covrV[T].size (),
4104+ MPI_DOUBLE_COMPLEX , MPI_SUM , MPI_COMM_WORLD );
41044105 }
41054106 // END: sample loop
41064107 }
@@ -4120,7 +4121,7 @@ namespace BasisClasses
41204121 for (int T=0 ; T<sampT; T++) {
41214122 ret[T].resize (1 );
41224123 std::get<0 >(ret[T][0 ]) = meanV[T];
4123- std::get<1 >(ret[T][0 ]) = covrV[T];
4124+ if (covar) std::get<1 >(ret[T][0 ]) = covrV[T];
41244125 }
41254126 }
41264127
@@ -4243,10 +4244,12 @@ namespace BasisClasses
42434244 v.resize (Itot);
42444245 }
42454246
4246- covrV.resize (sampT);
4247- for (auto & v : covrV) {
4248- v.resize (Itot, Itot);
4249- }
4247+ if (covar) {
4248+ covrV.resize (sampT);
4249+ for (auto & v : covrV) {
4250+ v.resize (Itot, Itot);
4251+ }
4252+ } else covrV.clear ();
42504253
42514254 sampleCounts.resize (sampT);
42524255 sampleMasses.resize (sampT);
@@ -4260,7 +4263,7 @@ namespace BasisClasses
42604263 {
42614264 for (int T=0 ; T<sampT; T++) {
42624265 meanV[T].setZero ();
4263- covrV[T].setZero ();
4266+ if (covar) covrV[T].setZero ();
42644267 }
42654268
42664269 sampleCounts.setZero ();
@@ -5109,6 +5112,11 @@ namespace BasisClasses
51095112 }
51105113 }
51115114
5115+ // Check for existence of a covariance matrix (only Cube can
5116+ // toggle this so far)
5117+ //
5118+ int varsz = std::get<1 >(covar[0 ][0 ]).size ();
5119+
51125120 // Pack the coefficient data
51135121 //
51145122 if (floatType) {
@@ -5131,23 +5139,26 @@ namespace BasisClasses
51315139
51325140 // Pack the covariance data in an upper triangular format
51335141 //
5134- real_part.resize (ltot*diagonalSize*sampleSize);
5135- imag_part.resize (ltot*diagonalSize*sampleSize);
5142+ if (varsz) {
51365143
5137- for (size_t T=0 , c=0 ; T<sampleCounts.size (); T++) {
5138- for (size_t l=0 ; l<ltot; l++) {
5139- for (size_t n1=0 ; n1<nmax; n1++) {
5140- for (size_t n2=n1; n2<nmax; n2++, c++) {
5141- real_part (c) = std::real (std::get<1 >(covar[T][l])(n1, n2));
5142- imag_part (c) = std::imag (std::get<1 >(covar[T][l])(n1, n2));
5144+ real_part.resize (ltot*diagonalSize*sampleSize);
5145+ imag_part.resize (ltot*diagonalSize*sampleSize);
5146+
5147+ for (size_t T=0 , c=0 ; T<sampleCounts.size (); T++) {
5148+ for (size_t l=0 ; l<ltot; l++) {
5149+ for (size_t n1=0 ; n1<nmax; n1++) {
5150+ for (size_t n2=n1; n2<nmax; n2++, c++) {
5151+ real_part (c) = std::real (std::get<1 >(covar[T][l])(n1, n2));
5152+ imag_part (c) = std::imag (std::get<1 >(covar[T][l])(n1, n2));
5153+ }
51435154 }
51445155 }
51455156 }
5146- }
51475157
5148- // Create two separate, compressed datasets
5149- stanza.createDataSet (" covariance_real" , real_part, dcpl3);
5150- stanza.createDataSet (" covariance_imag" , imag_part, dcpl3);
5158+ // Create two separate, compressed datasets
5159+ stanza.createDataSet (" covariance_real" , real_part, dcpl3);
5160+ stanza.createDataSet (" covariance_imag" , imag_part, dcpl3);
5161+ }
51515162
51525163 } else {
51535164 Eigen::VectorXd real_part (ltot*nmax*sampleSize);
@@ -5169,24 +5180,26 @@ namespace BasisClasses
51695180
51705181 // Pack the covariance data in an upper triangular format
51715182 //
5172- real_part.resize (ltot*diagonalSize*sampleSize);
5173- imag_part.resize (ltot*diagonalSize*sampleSize);
5174-
5175- for (size_t T=0 , c=0 ; T<sampleCounts.size (); T++) {
5176- for (size_t l=0 ; l<ltot; l++) {
5177- for (size_t n1=0 ; n1<nmax; n1++) {
5178- for (size_t n2=n1; n2<nmax; n2++, c++) {
5179- real_part (c) = std::real (std::get<1 >(covar[T][l])(n1, n2));
5180- imag_part (c) = std::imag (std::get<1 >(covar[T][l])(n1, n2));
5183+ if (varsz) {
5184+ real_part.resize (ltot*diagonalSize*sampleSize);
5185+ imag_part.resize (ltot*diagonalSize*sampleSize);
5186+
5187+ for (size_t T=0 , c=0 ; T<sampleCounts.size (); T++) {
5188+ for (size_t l=0 ; l<ltot; l++) {
5189+ for (size_t n1=0 ; n1<nmax; n1++) {
5190+ for (size_t n2=n1; n2<nmax; n2++, c++) {
5191+ real_part (c) = std::real (std::get<1 >(covar[T][l])(n1, n2));
5192+ imag_part (c) = std::imag (std::get<1 >(covar[T][l])(n1, n2));
5193+ }
51815194 }
51825195 }
51835196 }
5184- }
51855197
5186- // Create two separate, compressed datasets
5187- //
5188- stanza.createDataSet (" covariance_real" , real_part, dcpl3);
5189- stanza.createDataSet (" covariance_imag" , imag_part, dcpl3);
5198+ // Create two separate, compressed datasets
5199+ //
5200+ stanza.createDataSet (" covariance_real" , real_part, dcpl3);
5201+ stanza.createDataSet (" covariance_imag" , imag_part, dcpl3);
5202+ }
51905203 }
51915204 // END: sample loop
51925205
@@ -5436,16 +5449,21 @@ namespace BasisClasses
54365449 data0.real () = data_real.cast <double >();
54375450 data0.imag () = data_imag.cast <double >();
54385451
5439- data_real =
5440- stanza.getDataSet (" covariance_real" ).read <Eigen::VectorXf>();
5452+ // Check for existence of covariance
5453+ //
5454+ if (stanza.exist (" covariance_real" )) {
5455+
5456+ data_real =
5457+ stanza.getDataSet (" covariance_real" ).read <Eigen::VectorXf>();
54415458
5442- data_imag =
5443- stanza.getDataSet (" covariance_imag" ).read <Eigen::VectorXf>();
5459+ data_imag =
5460+ stanza.getDataSet (" covariance_imag" ).read <Eigen::VectorXf>();
54445461
5445- // Resize the complex array and assign
5446- data1.resize (data_real.size ());
5447- data1.real () = data_real.cast <double >();
5448- data1.imag () = data_imag.cast <double >();
5462+ // Resize the complex array and assign
5463+ data1.resize (data_real.size ());
5464+ data1.real () = data_real.cast <double >();
5465+ data1.imag () = data_imag.cast <double >();
5466+ }
54495467 } else {
54505468 // Get the real and imaginary parts
54515469 Eigen::VectorXd data_real =
@@ -5459,17 +5477,22 @@ namespace BasisClasses
54595477 data0.real () = data_real;
54605478 data0.imag () = data_imag;
54615479
5462- // Get the real and imaginary parts
5463- data_real =
5464- stanza.getDataSet (" covariance_real" ). read <Eigen::VectorXd>();
5480+ // Check for existence of covariance
5481+ //
5482+ if ( stanza.exist (" covariance_real" )) {
54655483
5466- data_imag =
5467- stanza.getDataSet (" covariance_imag" ).read <Eigen::VectorXd>();
5468-
5469- // Resize the complex array and assign
5470- data1.resize (data_real.size ());
5471- data1.real () = data_real;
5472- data1.imag () = data_imag;
5484+ // Get the real and imaginary parts
5485+ data_real =
5486+ stanza.getDataSet (" covariance_real" ).read <Eigen::VectorXd>();
5487+
5488+ data_imag =
5489+ stanza.getDataSet (" covariance_imag" ).read <Eigen::VectorXd>();
5490+
5491+ // Resize the complex array and assign
5492+ data1.resize (data_real.size ());
5493+ data1.real () = data_real;
5494+ data1.imag () = data_imag;
5495+ }
54735496 }
54745497
54755498 // Positions in data stanzas
@@ -5484,7 +5507,7 @@ namespace BasisClasses
54845507 // Coefficients
54855508 std::get<0 >(e).resize (rank);
54865509 // Covariance matrix
5487- std::get<1 >(e).resize (rank, rank);
5510+ if (data1. size ()) std::get<1 >(e).resize (rank, rank);
54885511 }
54895512
54905513 // Pack the coefficient data
@@ -5499,11 +5522,14 @@ namespace BasisClasses
54995522 // Pack the covariance data
55005523 c = 0 ;
55015524 for (size_t l=0 ; l<lSize; l++) {
5502- for (size_t n1=0 ; n1<rank; n1++) {
5503- for (size_t n2=n1; n2<rank; n2++) {
5504- std::get<1 >(elem[l])(n1, n2) = data1 (sCov + c++);
5505- if (n1 != n2)
5506- std::get<1 >(elem[l])(n2, n1) = std::get<1 >(elem[l])(n1, n2);
5525+
5526+ if (data1.size ()) {
5527+ for (size_t n1=0 ; n1<rank; n1++) {
5528+ for (size_t n2=n1; n2<rank; n2++) {
5529+ std::get<1 >(elem[l])(n1, n2) = data1 (sCov + c++);
5530+ if (n1 != n2)
5531+ std::get<1 >(elem[l])(n2, n1) = std::get<1 >(elem[l])(n1, n2);
5532+ }
55075533 }
55085534 }
55095535 }
0 commit comments