From cd1565c1447f5fe3c11325a8b131f50b44d39fff Mon Sep 17 00:00:00 2001 From: Wiarda <wiardada@ornl.gov> Date: Thu, 19 Aug 2021 17:27:26 -0400 Subject: [PATCH] Add some changes to GridList, DerivativeList to allow accelerated access Add a reserve function to the GridData class so that we can resize the vector more easily --- sammy/src/avg/mavg0.f | 5 +- sammy/src/clq/mclq0.f | 4 +- sammy/src/cro/mcro0.f | 2 +- sammy/src/dat/mdat0.f90 | 13 ++ sammy/src/dat/mdat3.f90 | 1 + sammy/src/dop/mdop0.f90 | 3 +- sammy/src/dop/mdop1.f90 | 2 + sammy/src/dop/mdop2.f90 | 3 + sammy/src/int/mint1.f | 6 + sammy/src/ipq/FillFitData_m.f90 | 3 +- sammy/src/mlb/mmlb0.f | 4 +- sammy/src/mxw/mmxw5.f | 1 + sammy/src/ntg/mntg4.f | 2 + sammy/src/ntg/mntg5.f | 10 +- sammy/src/ntg/mntg6.f | 18 +- sammy/src/rec/mrec5.f | 2 + sammy/src/rec/mrec8.f90 | 1 + sammy/src/rpt/mrpt.f | 1 + sammy/src/salmon/DerivativeList.cpp | 179 +++++++++++++++--- sammy/src/salmon/DerivativeList.h | 90 ++++++++- sammy/src/salmon/GridData.cpp | 61 +++++- sammy/src/salmon/GridData.h | 62 +++++- .../cix/DerivativeListHolder.cpp2f.xml | 20 +- .../salmon/interface/cix/GridData.cpp2f.xml | 10 +- .../cpp/DerivativeListHolderInterface.cpp | 22 ++- .../cpp/DerivativeListHolderInterface.h | 6 +- .../interface/cpp/GridDataInterface.cpp | 12 +- .../salmon/interface/cpp/GridDataInterface.h | 4 +- .../fortran/DerivativeListHolder_I.f90 | 30 ++- .../fortran/DerivativeListHolder_M.f90 | 54 ++++++ .../salmon/interface/fortran/GridData_I.f90 | 14 +- .../salmon/interface/fortran/GridData_M.f90 | 16 +- sammy/src/salmon/tests/DerivativeListTest.cpp | 154 ++++++++++++++- sammy/src/salmon/tests/GridDataTest.cpp | 82 ++++++++ sammy/src/smc/msmc0.f | 6 +- sammy/src/ssm/mssm18.f90 | 2 + sammy/src/the/DerivativeHandler.cpp | 97 ++++++++-- sammy/src/the/DerivativeHandler.h | 15 +- .../interface/cix/DerivativeHandler.cpp2f.xml | 7 + .../cpp/DerivativeHandlerInterface.cpp | 19 +- .../cpp/DerivativeHandlerInterface.h | 9 +- .../interface/fortran/DerivativeHandler_I.f90 | 14 +- .../interface/fortran/DerivativeHandler_M.f90 | 37 +++- sammy/src/xct/mxct02.f90 | 4 +- 44 files changed, 990 insertions(+), 117 deletions(-) diff --git a/sammy/src/avg/mavg0.f b/sammy/src/avg/mavg0.f index 69a346acc..e5ee6449e 100755 --- a/sammy/src/avg/mavg0.f +++ b/sammy/src/avg/mavg0.f @@ -14,7 +14,7 @@ C use oopsch_common_m, only : Nowwww, Segmen use AllocateFunctions_m, only : allocate_real_data, * allocate_integer_data - use EndfData_common_m, only : expData, resParData + use EndfData_common_m, only : expData, resParData, radFitFlags use DerivativeHandler_M use GridData_M use sammy_ipq_common_m, only : resultData, derivStart, nimplict @@ -96,7 +96,7 @@ C oldCalc%instance_ptr = calcData%instance_ptr call ourCalc%initialize() - call ourCalc%setUpList(resParData, 1) + call ourCalc%setUpList(resParData, radFitFlags, 1) call ourCalc%setNnsig(1) calcData%instance_ptr = ourCalc%instance_ptr C @@ -222,6 +222,7 @@ C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -< C C Copy theoretical data to auxillary grid C + call grid%reserve(resultData%getLength(), ipos) do i = 1, resultData%getLength() dd = resultData%getData(i, derivStart) call grid%addData(i, ipos, dd) diff --git a/sammy/src/clq/mclq0.f b/sammy/src/clq/mclq0.f index f09d19b81..2bc390003 100644 --- a/sammy/src/clq/mclq0.f +++ b/sammy/src/clq/mclq0.f @@ -17,7 +17,7 @@ C use templc_common_m use cbro_common_m use lbro_common_m - use EndfData_common_m, only : resParData + use EndfData_common_m, only : resParData,radFitFlags use rsl7_m use AuxGridHelper_M, only : setAuxGridOffset, setAuxGridRowMax IMPLICIT DOUBLE PRECISION (a-h,o-z) @@ -46,7 +46,7 @@ C *** GUESSTIMATE SIZE OF ARRAY NEEDED FOR SAMMY-CLQ C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - < CALL Set_Kws_Xct - call calcData%setUpList(resParData, Iq_Iso) + call calcData%setUpList(resParData, radFitFlags, Iq_Iso) CALL Work_Clq (A_Isigxx ) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > C diff --git a/sammy/src/cro/mcro0.f b/sammy/src/cro/mcro0.f index 18632f938..c3430c57d 100644 --- a/sammy/src/cro/mcro0.f +++ b/sammy/src/cro/mcro0.f @@ -72,7 +72,7 @@ C *** GUESSTIMATE SIZE OF ARRAY NEEDED FOR SAMMY-Cross section C C *** one *** CALL Set_Kws_Xct - call calcData%setUpList(resParData, Iq_Iso) + call calcData%setUpList(resParData, radFitFlags, Iq_Iso) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - < Krext = Nrext diff --git a/sammy/src/dat/mdat0.f90 b/sammy/src/dat/mdat0.f90 index cbf97f77e..bc26a117c 100644 --- a/sammy/src/dat/mdat0.f90 +++ b/sammy/src/dat/mdat0.f90 @@ -47,6 +47,7 @@ module Samdat_0_M type(EndfData)::reader type(GridData)::grid, auxGrid logical::hasAuxGrid + integer::maxRow, maxCol ! Segmen(1) = 'D' Segmen(2) = 'A' @@ -141,6 +142,12 @@ module Samdat_0_M if(Ktzero.NE.0) minDataIndex = minDataIndex + 1 if(Numcro.gt.1) minDataIndex = minDataIndex + 1 call grid%setDataIndex(minDataIndex) + maxRow = Ndat * Numcro + maxCol = 2 + if (Numcro.gt.1) maxCol = 3 + if (Ktzero.NE.0) maxCol = maxCol + 1 + call grid%reserve(maxRow, maxCol) + do i = 1, Ndat ee = A_Ienerg(i) etzero = ee @@ -333,6 +340,9 @@ module Samdat_0_M CALL Artificial_Energy (A_Ienerg, Ndat, Keveng) call expData%getGrid(grid, 1) ii = 1 + maxRow = Ndat * Numcro + maxCol = 1 + call grid%reserve(maxRow, maxCol) do i = 1, Ndat do j = 1, Numcro call grid%addData(ii, 1, A_Ienerg(i)) @@ -349,6 +359,9 @@ module Samdat_0_M if (hasAuxGrid) then call auxGrid%initialize() + maxRow = Ndatb * Numcro + maxCol = 1 + call auxGrid%reserve(maxRow, maxCol) ipos = 1 do i = 1, Ndatb do ii = 1, numcro diff --git a/sammy/src/dat/mdat3.f90 b/sammy/src/dat/mdat3.f90 index 08a4ddbf4..cf8475618 100755 --- a/sammy/src/dat/mdat3.f90 +++ b/sammy/src/dat/mdat3.f90 @@ -323,6 +323,7 @@ module mdat3_M end do CALL Convert_To_Caps (keyword, 10, Kpound) IF (keyword.EQ.'PARTIAL DE') THEN + call derivs%reserve(ener%getLength(), Numidc) DO K=1,ener%getLength() READ (30,10200) (X(1, I),I=1,Numidc) Do I = 1, Numidc diff --git a/sammy/src/dop/mdop0.f90 b/sammy/src/dop/mdop0.f90 index 7fcf7f8d7..740c0a487 100644 --- a/sammy/src/dop/mdop0.f90 +++ b/sammy/src/dop/mdop0.f90 @@ -58,7 +58,7 @@ module dop_m Left = Msize - Icoef - 2*I2pls1 Ldatb = Ndatb - ! ensure we have an auxillary energy grid + ! ensure we have an auxillary energy grid if ( expData%getLength().lt.2) then call expData%getGrid(grid, 1) ipos = 1 @@ -70,6 +70,7 @@ module dop_m end if end if call auxGrid%initialize() + call auxGrid%reserve(grid%getLength(),1) do i = 1, grid%getLength() call auxGrid%addData(i, 1, grid%getData(i, ipos)) end do diff --git a/sammy/src/dop/mdop1.f90 b/sammy/src/dop/mdop1.f90 index c2b58cada..411b65d71 100644 --- a/sammy/src/dop/mdop1.f90 +++ b/sammy/src/dop/mdop1.f90 @@ -263,12 +263,14 @@ module dop1_m ! iposJJ = (Jj - 1) * numcro + 1 ! all energies for each numcro are the same ee = grid%getData(iposJJ, 1) + call grid%reserve(iposJJ + Numcro,1) do Jjnum = iposJJ, iposJJ + Numcro -1 call grid%addData( Jjnum, 1, ee) end do J = J - 1 END DO iposJJ = (Ngtvv - 1) * numcro + 1 ! all energies for each numcro are the same + call grid%reserve(iposJJ + Numcro,1) do Jjnum = iposJJ, iposJJ + Numcro -1 call grid%addData( Jjnum, 1, 0.0d0) end do diff --git a/sammy/src/dop/mdop2.f90 b/sammy/src/dop/mdop2.f90 index 5e3b6bdfd..290e00f2d 100755 --- a/sammy/src/dop/mdop2.f90 +++ b/sammy/src/dop/mdop2.f90 @@ -26,6 +26,7 @@ module dop2_m Mstopx = 0 Np1 = N + 1 Itimes = 0 + call Coef%reserve(Np1, 1) DO I=1,Np1 IF((I/1000)*1000.EQ.I) WRITE (6,66666) I, Coef%getData(I-1,1) 66666 FORMAT (' I, Coef(I-1)=', i10, 1PG20.10) @@ -102,6 +103,7 @@ module dop2_m ! *** rearrange storage of Coef ... 1 to m+1 now goes to m+1 to 2m+1 C = Coef%getData(Np1, 1) Ii = Np1 + call Coef%reserve(Ii+Mmmmmm,1) DO I=1,Mmmmmm tmp = Coef%getData(I,1) call Coef%addData(Ii,1,tmp) @@ -111,6 +113,7 @@ module dop2_m ! ! *** now put M+2 to 2M+1 into [1 to M] in reverse order Ii = M2pls1 + call Coef%reserve(Ii+Mmmmmm,1) DO I=1,Mmmmmm tmp = Coef%getData(Ii,1) call Coef%addData(I, 1, tmp) diff --git a/sammy/src/int/mint1.f b/sammy/src/int/mint1.f index 30461a1dc..559f1ca44 100644 --- a/sammy/src/int/mint1.f +++ b/sammy/src/int/mint1.f @@ -162,6 +162,12 @@ C C DO Iso=1,calcData%getNumberIsotopes() ! update cross section (ipar=0) as well as derivatives + DO N=1,calcData%getNnnsig() + call calcData%reserveColumnsNs(Idat,N, Ndasig+Ndbsig) + end do + if (Ksindi.ne.0) then + call calcDataSelf%reserveColumns(Idat,Ndasig+Ndbsig) + end if DO Ipar=0,Ndasig+Ndbsig DO N=1,calcData%getNnnsig() val1 = calcData%getDataNs(Iw1, N, Ipar, Iso) diff --git a/sammy/src/ipq/FillFitData_m.f90 b/sammy/src/ipq/FillFitData_m.f90 index 427c8076f..bd2082ae3 100644 --- a/sammy/src/ipq/FillFitData_m.f90 +++ b/sammy/src/ipq/FillFitData_m.f90 @@ -117,8 +117,9 @@ contains ! copy the user supplied deriviatives as needed call grid%getImplicitDerivs(implGrid) if (C_ASSOCIATED(implGrid%instance_ptr)) then + call resultData%reserve(implGrid%getLength(), nimplgiven+npars+derivStart) do i = 1, implGrid%getLength() - do j = 1, nimplgiven + do j = 1, nimplgiven call resultData%addData(i, j + npars + derivStart, implGrid%getData(i,j)) end do end do diff --git a/sammy/src/mlb/mmlb0.f b/sammy/src/mlb/mmlb0.f index 16a102b3f..a719821fe 100644 --- a/sammy/src/mlb/mmlb0.f +++ b/sammy/src/mlb/mmlb0.f @@ -22,7 +22,7 @@ C use rsl7_m use xct_m use mxct27_m - use EndfData_common_m, only : covData, resParData + use EndfData_common_m, only : covData, resParData, radFitFlags use AuxGridHelper_M, only : setAuxGridOffset, setAuxGridRowMax IMPLICIT DOUBLE PRECISION (a-h,o-z) real(kind=8),allocatable,dimension(:)::A_Idum @@ -74,7 +74,7 @@ C C ### one ### Ks_Res = Ksolve CALL Set_Kws_Xct - call calcData%setUpList(resParData, Iq_Iso) + call calcData%setUpList(resParData, radFitFlags, Iq_Iso) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - < C C ### three ### diff --git a/sammy/src/mxw/mmxw5.f b/sammy/src/mxw/mmxw5.f index ef60f8946..dc73df8ef 100644 --- a/sammy/src/mxw/mmxw5.f +++ b/sammy/src/mxw/mmxw5.f @@ -170,6 +170,7 @@ C Gj1 = Gj1 + Aa*val2 val = calcData%getDataNs(J, 1, 0, 1) val = val + Aa*val2 + call calcData%reserveColumnsNs(J, 1, Ndasig+Ndbsig) call calcData%addDataNs(J, 1, 0, 1, val) IF (Ksolve.NE.2) then DO Ipar=1,Ndasig + Ndbsig diff --git a/sammy/src/ntg/mntg4.f b/sammy/src/ntg/mntg4.f index edd0039a9..7125fa1b6 100644 --- a/sammy/src/ntg/mntg4.f +++ b/sammy/src/ntg/mntg4.f @@ -462,6 +462,8 @@ C Aaold = Bb val1 = val1 + Alfa*Aa val2 = val2 + Etaa*Aa*Xnu + call calcData%reserveColumns(Kountr, Ndasig+Ndbsig) + call calcData%reserveColumns(Kountr+1, Ndasig+Ndbsig) call calcData%addData(Kountr, 0, 1, val1) call calcData%addData(Kountr+1, 0, 1, val2) DO Ipar=1,Ndasig + Ndbsig diff --git a/sammy/src/ntg/mntg5.f b/sammy/src/ntg/mntg5.f index 75bb258a8..68d682741 100644 --- a/sammy/src/ntg/mntg5.f +++ b/sammy/src/ntg/mntg5.f @@ -110,6 +110,7 @@ C W12 = calcData%getData(Kountr-12, 0, 1) W13 = calcData%getData(Kountr-13, 0, 1) val = (Xnu*W12-W13)*Tosp + call calcData%reserveColumns(Kountr, Ndasig+Ndbsig) call calcData%addData(Kountr, 0, 1, val) IF (Ndasig.NE.0) THEN DO Ipar=1,Ndasig @@ -137,9 +138,7 @@ C Kountr = Kountr + 1 if (calcData%getLength().lt.Kountr) then ! update cross section (ipar=0) as well as derivatives - do Ipar = 0, Ndasig + Ndbsig - call calcData%setSharedDataVal(Kountr, Ipar, 0.0d0) - end do + call calcData%reserveColumns(Kountr, Ndasig+Ndbsig) end if RETURN END @@ -164,6 +163,7 @@ C IF (W16.EQ.0.0d0) return W15 = calcData%getData(Kountr-15,0,1) WK = W15 / W16 + call calcData%reserveColumns(Kountr, Ndasig+Ndbsig) call calcData%addData(Kountr,0,1,WK) IF (Ndasig.GE.0) THEN DO Ipar=1,Ndasig @@ -178,9 +178,7 @@ C Kountr = Kountr + 1 if (calcData%getLength().lt.Kountr) then ! update cross section (ipar=0) as well as derivatives - do Ipar = 0, Ndasig + Ndbsig - call calcData%setSharedDataVal(Kountr, Ipar, 0.0d0) - end do + call calcData%reserveColumns(Kountr, Ndasig+Ndbsig) end if RETURN END diff --git a/sammy/src/ntg/mntg6.f b/sammy/src/ntg/mntg6.f index 34f104271..c23b02618 100644 --- a/sammy/src/ntg/mntg6.f +++ b/sammy/src/ntg/mntg6.f @@ -152,6 +152,7 @@ C C IF (Nnnsig.EQ.1) THEN val = calcData%getData(Kountr, 0, 1) + Const(2,Kountr) + call calcData%reserveColumns(Kountr, Ndasig+Ndbsig) call calcData%addData(Kountr, 0, 1, val) WRITE (21,10400) val Kountr = Kountr + 1 @@ -165,16 +166,13 @@ C ELSE eIp0 = grid%getEnergy(JEquit , expData) val = calcData%getData(Kountr, 0, 1) + Const(2,Kountr) + call calcData%reserveColumns(Kountr, Ndasig+Ndbsig) call calcData%addData(Kountr, 0, 1, val) WRITE (21,10200) val, eIp0, Equit 10200 FORMAT (' RIABS = Resonance integrl Absrptn =', 1P4G14.6) Kountr = Kountr + 1 - if (calcData%getLength().lt.Kountr) then - ! update cross section (ipar=0) as well as derivatives - do Ipar = 0, Ndasig + Ndbsig - call calcData%setSharedDataVal(Kountr, Ipar, 0.0d0) - end do - end if + ! update cross section (ipar=0) as well as derivatives + call calcData%reserveColumns(Kountr, Ndasig+Ndbsig) val = calcData%getData(Kountr, 0, 1) + Const(2,Kountr) call calcData%addData(Kountr, 0, 1, val) WRITE (21,10300) val, eIp0, Equit @@ -191,12 +189,8 @@ C WRITE (21,10400) val, eIp0, Equit 10400 FORMAT (' RICAP = Resonance integrl Capture =', 1P4G14.6) Kountr = Kountr + 1 - if (calcData%getLength().lt.Kountr) then - ! update cross section (ipar=0) as well as derivatives - do Ipar = 0, Ndasig + Ndbsig - call calcData%setSharedDataVal(Kountr, Ipar, 0.0d0) - end do - end if + ! update cross section (ipar=0) as well as derivatives + call calcData%reserveColumns(Kountr, Ndasig+Ndbsig) END IF call grid%destroy() RETURN diff --git a/sammy/src/rec/mrec5.f b/sammy/src/rec/mrec5.f index b56729918..1b7d400dd 100644 --- a/sammy/src/rec/mrec5.f +++ b/sammy/src/rec/mrec5.f @@ -101,6 +101,7 @@ C *** Sig2(2,I) is capture cross section (absorption - fission) call expData%initialize() call grid%initialize() index = 1 + call grid%reserve(Ndatb*Numcro, 1) do i = 1, Ndatb do j = 1, Numcro call grid%addData(index, 1, Etab2(i)) @@ -115,6 +116,7 @@ C *** Sig2(2,I) is capture cross section (absorption - fission) end if if (calcData%getNumberIsotopes().eq.0) then call calcData%setUpList(resParData, + * radFitFlags, * resParData%getNumIso()) end if Ksolve = 2 diff --git a/sammy/src/rec/mrec8.f90 b/sammy/src/rec/mrec8.f90 index e738ff0fc..449b8673b 100755 --- a/sammy/src/rec/mrec8.f90 +++ b/sammy/src/rec/mrec8.f90 @@ -87,6 +87,7 @@ module mrec8_m 101 CONTINUE Npurr = L - 1 call reallocate_real_data(X, Npurr) + call Energb%reserve(Nemax+Npurr, 1) DO N=1,Npurr call Energb%addData(Nemax+N, 1, X(N)) END DO diff --git a/sammy/src/rpt/mrpt.f b/sammy/src/rpt/mrpt.f index 11b4de815..26d927caf 100644 --- a/sammy/src/rpt/mrpt.f +++ b/sammy/src/rpt/mrpt.f @@ -294,6 +294,7 @@ C T = Tup Dt = (Tup-Tlow)/(dFLOAT(N-1)) Ii = 0 + call grid%reserve(N, 1) DO I=1,N Ii = Ii + 1 Tt(Ii) = T - Dt*dFLOAT(I-1) diff --git a/sammy/src/salmon/DerivativeList.cpp b/sammy/src/salmon/DerivativeList.cpp index e0e04167a..1caac597f 100644 --- a/sammy/src/salmon/DerivativeList.cpp +++ b/sammy/src/salmon/DerivativeList.cpp @@ -9,7 +9,7 @@ namespace sammy { void DerivativeList::addGrid(){ std::unique_ptr<GridData> data = nemesis::make_unique<GridData>(); - gridData.addGrid(data); + gridData.addGrid(data); // gridData automatically sets notSetReturnsZero } int DerivativeList::getLength() const{ @@ -24,30 +24,35 @@ namespace sammy { return gridData.getGrid(iso)->getLength(); } - void DerivativeList::updateRealCol(int col){ - if (col < (int)realColIndices.size()) return; // already done + void DerivativeList::updateRealCol(int col, bool force){ + if (col < (int)realColIndices.size() && !force) return; // already done - size_t start = realColIndices.size(); - for (int c = start; c <= col; c++){ - if (isSharedColumn(c)){ - for (size_t i = 0; i < sharedIndices.size(); i++){ - if ( c == sharedIndices[i].first){ - realColIndices.push_back(-1*(i+1)); - break; - - } - } - } - else{ - int pos = -1; - for ( int c = 0; c <= col; c++){ - if ( isSharedColumn(c)) continue; - pos++; - } - realColIndices.push_back(pos+1); + size_t start = realColIndices.size(); + if (force) { + start = 0; + realColIndices.clear(); } - } - } + for (int c = start; c <= col; c++){ + int index = -1; + if (isSharedColumn(c)){ + for (size_t i = 0; i < sharedIndices.size(); i++){ + if ( c == sharedIndices[i].first){ + index = -1*(i+1); + + } + } + } + else{ + int pos = -1; + for ( int i = 0; i <= c; i++){ + if ( isSharedColumn(i)) continue; + pos++; + } + index = pos+1; + } + realColIndices.push_back(index); + } + } void DerivativeList::addData(int row, int col, int iso, double val){ @@ -56,6 +61,11 @@ namespace sammy { throw std::runtime_error(std::to_string(iso) + " Isotope out of range"); } + if (sharedIndices.size() == 0){ + gridData.getGrid(iso)->addData(row,col, val); + return; + } + updateRealCol(col); int ourCol = realColIndices[col]; @@ -75,20 +85,89 @@ namespace sammy { gridData.getGrid(iso)->addData(row,ourCol, val); } + namespace detail{ + void getMaxColumnData(int & maxShared, int & maxUn, const std::vector<int> & cols){ + maxShared = maxUn = 0; + auto min = std::min_element(cols.begin(), cols.end()); + int minVal = (*min); + if (minVal < 0){ + maxShared = -1 * minVal; + } + auto max = std::max_element(cols.begin(), cols.end()); + maxUn = (*max); + } + } + + void DerivativeList::reserveColumns(int row, int maxCol){ + if(maxCol == 0) return; + + if (sharedIndices.size() == 0){ + for (int iso = 0; iso < gridData.getLength(); iso++){ + gridData.getGrid(iso)->reserveColumn(row, maxCol); + } + return; + } + + updateRealCol(maxCol-1); + + int maxShared, maxUn; + detail::getMaxColumnData(maxShared, maxUn, realColIndices); + if (maxShared > 0){ + sharedData.reserveColumn(row, maxShared); + } + if (maxUn > 0){ + for (int iso = 0; iso < gridData.getLength(); iso++){ + gridData.getGrid(iso)->reserveColumn(row, maxUn); + } + } + } + + + void DerivativeList::reserve(int maxRow, int maxCol){ + if(maxCol == 0 || maxRow == 0) return; + + if (sharedIndices.size() == 0){ + for (int iso = 0; iso < gridData.getLength(); iso++){ + gridData.getGrid(iso)->reserve(maxRow, maxCol); + } + return; + } + + + updateRealCol(maxCol-1); + + int maxShared, maxUn; + detail::getMaxColumnData(maxShared, maxUn, realColIndices); + if (maxShared > 0){ + sharedData.reserve(maxRow, maxShared); + } + if (maxUn > 0){ + for (int iso = 0; iso < gridData.getLength(); iso++){ + gridData.getGrid(iso)->reserve(maxRow, maxUn); + } + } + } + double DerivativeList::getData(int row, int col, int iso) const{ // check whether it's a valid isotope if (iso < 0 || iso >= gridData.getLength()){ throw std::runtime_error(std::to_string(iso) + " Isotope out of range"); } + + if (sharedIndices.size() == 0){ + return gridData.getGrid(iso)->getData(row,col); + } + if (col >= (int)realColIndices.size()){ + if( notSetReturnsZero) return 0.0; throw std::runtime_error(std::to_string(col) + " column out of range"); } // if it is a shared column, just get it - int ourCol = realColIndices[col]; + int ourCol = realColIndices[col]; if (ourCol < 0){ ourCol = -1*ourCol - 1; - if (sharedIndices[ourCol].second == iso){ + if (sharedIndices[ourCol].second == iso){ return sharedData.getData(row, ourCol); } return 0.0; @@ -131,7 +210,7 @@ namespace sammy { // Therefore the C++ setup get's a chance to reset the // isotope consistent with the isotope ID in the // resonance data - auto it = std::find_if( sharedIndices.begin(), sharedIndices.end(), detail::comp(col)); + auto it = std::find_if( sharedIndices.begin(), sharedIndices.end(), detail::comp(col)); it->second = iso; return; } @@ -142,7 +221,9 @@ namespace sammy { throw std::runtime_error("Can't add shared column if we already have data"); } sharedIndices.push_back(std::make_pair(col, iso)); - updateRealCol(col); + bool force = false; + if (col < (int)realColIndices.size() ) force = true; + updateRealCol(col, force); } void DerivativeList::clear(){ @@ -155,10 +236,10 @@ namespace sammy { } void DerivativeList::nullify(){ - sharedData.clearGrid(); + sharedData.nullify(); for ( int i = 0; i < gridData.getLength(); i++){ - gridData.getGrid(i)->clearGrid(); + gridData.getGrid(i)->nullify(); } } @@ -185,8 +266,32 @@ namespace sammy { return -1; } + void DerivativeList::setNotSetReturnsZero(bool empty){ + notSetReturnsZero = empty; + gridData.setNotSetReturnsZero(empty); + sharedData.setNotSetReturnsZero(empty); + } + + void DerivativeList::setAccumulate(bool add){ + accumulate = add; + gridData.setAccumulate(add); + sharedData.setAccumulate(add); + } + + void DerivativeListHolder::setNotSetReturnsZero(bool empty){ + notSetReturnsZero = empty; + list[0].setNotSetReturnsZero(empty); + list[1].setNotSetReturnsZero(empty); + } + + void DerivativeListHolder::setAccumulate(bool add){ + accumulate = add; + list[0].setAccumulate(add); + list[1].setAccumulate(add); + } + void DerivativeListHolder::addGrid(){ - list[currentList].addGrid(); + list[currentList].addGrid(); // list updates notSetReturnsZero internally } int DerivativeListHolder::getLength(bool current) const{ @@ -202,6 +307,20 @@ namespace sammy { return list[c].getLength(iso); } + void DerivativeListHolder::reserveColumns(int row, int maxCol, bool current){ + int c = currentList; + if( !current) c = parentList; + + list[c].reserveColumns(row, maxCol); + } + + void DerivativeListHolder::reserve(int maxRow, int maxCol, bool current){ + int c = currentList; + if( !current) c = parentList; + + list[c].reserve(maxRow, maxCol); + } + void DerivativeListHolder::addData(int row, int col, int iso, double val, bool current){ int c = currentList; if( !current) c = parentList; diff --git a/sammy/src/salmon/DerivativeList.h b/sammy/src/salmon/DerivativeList.h index fbafd20c7..e1725cd67 100644 --- a/sammy/src/salmon/DerivativeList.h +++ b/sammy/src/salmon/DerivativeList.h @@ -47,7 +47,7 @@ namespace sammy{ */ class DerivativeList{ public: - DerivativeList(){}; + DerivativeList():notSetReturnsZero(false),accumulate(false){}; DerivativeList(const DerivativeList & orig) = delete; virtual ~DerivativeList(){} @@ -148,10 +148,47 @@ namespace sammy{ * a shared column */ int getIsotopeForShared(int col) const; + + /** + * Indicate whether getData should return 0 instead of throw if the + * corresponding addData was not called. + * + * This is sometime convenient if calculating derivatives, where + * derivatives are only set for paramters where the matter, + * but the array needs to be over all flagged parameters. + * + * @param empty if true, getData returns 0 if addData was not called + */ + void setNotSetReturnsZero(bool empty); + + /** + * Should addData accumulate instead of set + * @param add + */ + void setAccumulate(bool add); + + /** + * Make sure the indicated row as at least maxCol data. + * If not set, additional columns are set to zero + * + * @param row the row for which to set the data + * @param maxCol the maximum number of columns + */ + void reserveColumns(int row, int maxCol); + + /** + * Reserve maxRow and maxCol. + * If not set, additional rows columns are set to zero + * + * @param maxRow the maximum number of rows + * @param maxCol the maximum number of columns + * @return + */ + void reserve(int maxRow, int maxCol); private: /* For look-up acceleration set up the linking between columns and * isotope to either shared or "normal colum */ - void updateRealCol(int col); + void updateRealCol(int col, bool force = false); /** the list of "normal" columns */ GridDataList gridData; @@ -162,6 +199,10 @@ namespace sammy{ /** The isotope to link to for shared column */ std::vector<std::pair<int,int>> sharedIndices; + bool notSetReturnsZero; + + bool accumulate; + /* For look-up acceleration set up the linking between columns and * isotope to either shared or "normal colum */ std::vector<int> realColIndices; @@ -184,7 +225,7 @@ namespace sammy{ */ class DerivativeListHolder{ public: - DerivativeListHolder():currentList(0),parentList(0){} + DerivativeListHolder():currentList(0),parentList(0), notSetReturnsZero(false), accumulate(false){} DerivativeListHolder(const DerivativeListHolder & orig) = delete; virtual ~DerivativeListHolder(){} @@ -256,7 +297,7 @@ namespace sammy{ /** * Set the length of the desired DerivativeList objects to 0 - * @param current true if using the current grif, false for the parent grid + * @param current true if using the current grid, false for the parent grid */ void nullify(bool current = true); @@ -282,6 +323,43 @@ namespace sammy{ * @return isotope index for the shared column */ int getIsotopeForShared(int col) const; + + /** + * Indicate whether getData should return 0 instead of throw if the + * corresponding addData was not called. + * + * This is sometime convenient if calculating derivatives, where + * derivatives are only set for paramters where the matter, + * but the array needs to be over all flagged parameters. + * + * @param empty if true, getData returns 0 if addData was not called + */ + void setNotSetReturnsZero(bool empty); + + /** + * Should addData accumulate instead of just set + * @param add + */ + void setAccumulate(bool add); + + /** + * Make sure the indicated row as at least maxCol data. + * If not set, additional columns are set to zero + * + * @param row the row for which to set the data + * @param maxCol the maximum number of columns + */ + void reserveColumns(int row, int maxCol, bool current = true); + + /** + * Reserve maxRow and maxCol. + * If not set, additional rows columns are set to zero + * + * @param maxRow the maximum number of rows + * @param maxCol the maximum number of columns + * @return + */ + void reserve(int maxRow, int maxCol, bool current = true); private: /** Index of the current index in list used to accumulate the broadened data */ int currentList; @@ -289,6 +367,10 @@ namespace sammy{ /** Index of the index to the list containing the unbroadened data */ int parentList; + bool notSetReturnsZero; + + bool accumulate; + /** The two grids to ping-pong between if chaining broadening data */ DerivativeList list[2]; diff --git a/sammy/src/salmon/GridData.cpp b/sammy/src/salmon/GridData.cpp index 1f7380cdd..940f6591e 100644 --- a/sammy/src/salmon/GridData.cpp +++ b/sammy/src/salmon/GridData.cpp @@ -1,11 +1,14 @@ #include "GridData.h" #include "Nemesis/utils/Future.hh" +#include <algorithm> namespace sammy{ GridData::GridData(const GridData & orig):data(orig.data), dataIndex(orig.dataIndex), rowOffset(orig.rowOffset), - rowMax(orig.rowMax){ + rowMax(orig.rowMax), + notSetReturnsZero(orig.notSetReturnsZero), + accumulate(orig.accumulate){ if (orig.implicitParCov != nullptr){ implicitParCov = nemesis::make_unique<endf::ResonanceCovariance>(*(orig.implicitParCov)); @@ -23,6 +26,9 @@ namespace sammy{ if (row < 0) throw std::runtime_error("row index " + std::to_string(row) + " out of range in addData"); if (col < 0) throw std::runtime_error("col index " + std::to_string(col) + " out of range in addData"); + if (accumulate && val == 0.0) return; + + if ( row >= getLength() ) { data.resize(row+1); @@ -31,14 +37,40 @@ namespace sammy{ data[row].resize(col+1); } - data[row][col] = val; + + if (accumulate) data[row][col] += val; + else data[row][col] = val; + } + + + void GridData::reserve(int maxRow, int maxCol){ + if (maxRow > (int)data.size()){ + data.resize(maxRow); + } + + for( auto & d : data){ + if (maxCol > (int)d.size()){ + d.resize(maxCol); + } + } + } + + void GridData::reserveColumn(int row, int maxCol){ + if (row >= (int)data.size()){ + data.resize(row+1); + } + if( maxCol > (int)data[row].size()){ + data[row].resize(maxCol); + } } double GridData::getData(int row, int col) const{ if ( row < 0 || row >= (int)data.size()) { + if (notSetReturnsZero) return 0.0; throw std::runtime_error("Row index out of range in GridData::getData: " + std::to_string(row)); } if (col < 0 || col >= (int)data[row].size()){ + if (notSetReturnsZero) return 0.0; throw std::runtime_error("Column index out of range in GridData::getData: " + std::to_string(col)); } @@ -73,10 +105,31 @@ namespace sammy{ data.clear(); } - GridDataList::GridDataList():auxGridIndex(1), expGridIndex(0){ // normally auxillary grid index is 1, and experimental grid 0 + void GridData::nullify(){ + for (auto & d: data){ + std::fill(d.begin(), d.end(), 0.0); + } + } + + GridDataList::GridDataList():auxGridIndex(1), expGridIndex(0), notSetReturnsZero(false),accumulate(false){ // normally auxillary grid index is 1, and experimental grid 0 expCov = std::unique_ptr<endf::ResonanceCovariance>( new endf::ResonanceCovariance()); } + void GridDataList::setNotSetReturnsZero(bool empty){ + notSetReturnsZero = empty; + for (auto & g : grids){ + g->setNotSetReturnsZero(empty); + } + } + + void GridDataList::setAccumulate(bool add){ + accumulate = add; + + for (auto & g : grids){ + g->setAccumulate(add); + } + } + int GridDataList::getLength() const{ return (int)grids.size(); @@ -100,6 +153,8 @@ namespace sammy{ void GridDataList::addGrid(std::unique_ptr<GridData> & grid){ grids.push_back( std::unique_ptr<GridData>()); grids[grids.size() -1 ] = std::move(grid); + grids[grids.size() -1 ]->setNotSetReturnsZero(notSetReturnsZero); + grids[grids.size() -1 ]->setAccumulate(accumulate); } void GridDataList::addExperimentalCov(int row, int col, double cov){ diff --git a/sammy/src/salmon/GridData.h b/sammy/src/salmon/GridData.h index c4f39383b..cc9ed0524 100644 --- a/sammy/src/salmon/GridData.h +++ b/sammy/src/salmon/GridData.h @@ -31,7 +31,7 @@ namespace sammy{ */ class GridData { public: - GridData():dataIndex(-1), rowOffset(0), rowMax(0){} + GridData():dataIndex(-1), rowOffset(0), rowMax(0), notSetReturnsZero(false),accumulate(false){} GridData(const GridData & orig); virtual ~GridData(){} @@ -64,6 +64,30 @@ namespace sammy{ */ double getData(int row, int col) const; + /** + * Indicate whether getData should return 0 instead of throw if the + * corresponding addData was not called. + * + * This is sometime convenient if calculating derivatives, where + * derivatives are only set for paramters where the matter, + * but the array needs to be over all flagged parameters. + * + * @param empty if true, getData returns 0 if addData was not called + */ + void setNotSetReturnsZero(bool empty){ + notSetReturnsZero = empty; + } + + + /** + * If true, addData accumulates data instead of set data + * + * @param add set to true to accumulate data with addData + */ + void setAccumulate(bool add){ + accumulate = add; + } + /** * @brief getDataColumn Get the index at which the data start * @return index at which the data start @@ -101,6 +125,7 @@ namespace sammy{ int getRowMax() const; void clearGrid(); + void nullify(); /** * Add a covariance for implicit parameters for @@ -136,6 +161,15 @@ namespace sammy{ std::unique_ptr<GridData> & getImplicitDerivs(){ return implicitDerivs; } + + /** + * Reserve (and fill with zero), maxRow and maxCol data + * + * @param maxRow the number of rows desired + * @param maxCol the number of columns desired + */ + void reserve(int maxRow, int maxCol); + void reserveColumn(int row, int maxCol); protected: std::vector< std::vector<double> > data; @@ -145,6 +179,9 @@ namespace sammy{ int rowOffset; int rowMax; + bool notSetReturnsZero; + bool accumulate; + // arrays for the implicit data covariance // covariance for the implicit parameters @@ -187,6 +224,26 @@ namespace sammy{ */ void addGrid(std::unique_ptr<GridData> & grid); + /** + * Indicate whether getData should return 0 instead of throw if the + * corresponding addData was not called. + * + * This is sometime convenient if calculating derivatives, where + * derivatives are only set for paramters where the matter, + * but the array needs to be over all flagged parameters. + * + * @param empty if true, getData returns 0 if addData was not called + */ + void setNotSetReturnsZero(bool empty); + + /** + * If true, addData for all grids accumulates data, instead of + * just updating them + * + * @param add set to true to accumulate in addData + */ + void setAccumulate(bool add); + /** * @brief addExperimentalCov Add or set covariance data * Currently I am thinking covariance data are for all grids, and index is starting at grid index 0, gridsize of 1. grid + [0, 1, ..] for @@ -267,6 +324,9 @@ namespace sammy{ int auxGridIndex; int expGridIndex; + + bool notSetReturnsZero; + bool accumulate; }; } diff --git a/sammy/src/salmon/interface/cix/DerivativeListHolder.cpp2f.xml b/sammy/src/salmon/interface/cix/DerivativeListHolder.cpp2f.xml index 74c3fe0e7..4318755e6 100644 --- a/sammy/src/salmon/interface/cix/DerivativeListHolder.cpp2f.xml +++ b/sammy/src/salmon/interface/cix/DerivativeListHolder.cpp2f.xml @@ -19,6 +19,16 @@ <param name="iso" type="int" offset="-1"/> <param name="current" type="bool"/> </method> + <method name="reserveColumns"> + <param name="row" type="int" offset="-1"/> + <param name="maxCol" type="int" offset="+1"/> + <param name="current" type="bool"/> + </method> + <method name="reserve"> + <param name="maxRow" type="int"/> + <param name="maxCol" type="int" offset="+1"/> + <param name="current" type="bool"/> + </method> <method name="clear"/> <method name="nullify"> <param name="current" type="bool"/> @@ -28,8 +38,14 @@ <param name="col" type="int"/> </method> <method name="addSharedColumn"> - <param name="col" type="int"/> - <param name="iso" type="int" offset="-1"/> + <param name="col" type="int"/> + <param name="iso" type="int" offset="-1"/> </method> + <method name="setNotSetReturnsZero"> + <param name="empty" type="bool"/> + </method> + <method name="setAccumulate"> + <param name="add" type="bool"/> + </method> </class> </generate> diff --git a/sammy/src/salmon/interface/cix/GridData.cpp2f.xml b/sammy/src/salmon/interface/cix/GridData.cpp2f.xml index 4c659b249..f47920359 100644 --- a/sammy/src/salmon/interface/cix/GridData.cpp2f.xml +++ b/sammy/src/salmon/interface/cix/GridData.cpp2f.xml @@ -26,6 +26,12 @@ <param name="offset" type="int"/> </method> <method name="clearGrid"/> + <method name="nullify"/> + + <method name="reserve"> + <param name="row" type="int"/> + <param name="col" type="int"/> + </method> <method name="addImplicitDerivs"> <param name="grid" type="GridData*"/> @@ -33,7 +39,7 @@ <method name="getImplicitDerivs" return_type="GridData*"/> <method name="addImplicitParCov"> <param name="cov" type="ResonanceCovariance*"/> - </method> + </method> <method name="getImplicitParCov" return_type="ResonanceCovariance*"/> </class> @@ -53,7 +59,7 @@ <method name="getExperimentalCov" return_type="double"> <param name="row" type="int" offset="-1"/> <param name="col" type="int" offset="-1"/> - </method> + </method> <method name="getAuxGridIndex" return_type="int"/> <method name="getExpGridIndex" return_type="int"/> <method name="setAuxGridIndex"> diff --git a/sammy/src/salmon/interface/cpp/DerivativeListHolderInterface.cpp b/sammy/src/salmon/interface/cpp/DerivativeListHolderInterface.cpp index e4edba3b7..b1b749ac5 100644 --- a/sammy/src/salmon/interface/cpp/DerivativeListHolderInterface.cpp +++ b/sammy/src/salmon/interface/cpp/DerivativeListHolderInterface.cpp @@ -2,7 +2,7 @@ * This file has been dynamically generated by Class Interface Xml (CIX) * DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION * If changes need to occur, modify the appropriate CIX xml file -* Date Generated: Tue Apr 06 18:09:02 EDT 2021 +* Date Generated: Fri Aug 20 14:10:56 EDT 2021 * If any issues are experiences with this generated file that cannot be fixed * with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov */ @@ -24,6 +24,16 @@ double DerivativeListHolder_getData(void * DerivativeListHolder_ptr,int * row,in return ((DerivativeListHolder*)DerivativeListHolder_ptr)->getData(*row,*col,*iso,*current); } +void DerivativeListHolder_reserveColumns(void * DerivativeListHolder_ptr,int * row,int * maxCol,bool * current) +{ + ((DerivativeListHolder*)DerivativeListHolder_ptr)->reserveColumns(*row,*maxCol,*current); +} + +void DerivativeListHolder_reserve(void * DerivativeListHolder_ptr,int * maxRow,int * maxCol,bool * current) +{ + ((DerivativeListHolder*)DerivativeListHolder_ptr)->reserve(*maxRow,*maxCol,*current); +} + void DerivativeListHolder_clear(void * DerivativeListHolder_ptr) { ((DerivativeListHolder*)DerivativeListHolder_ptr)->clear(); @@ -49,6 +59,16 @@ void DerivativeListHolder_addSharedColumn(void * DerivativeListHolder_ptr,int * ((DerivativeListHolder*)DerivativeListHolder_ptr)->addSharedColumn(*col,*iso); } +void DerivativeListHolder_setNotSetReturnsZero(void * DerivativeListHolder_ptr,bool * empty) +{ + ((DerivativeListHolder*)DerivativeListHolder_ptr)->setNotSetReturnsZero(*empty); +} + +void DerivativeListHolder_setAccumulate(void * DerivativeListHolder_ptr,bool * add) +{ + ((DerivativeListHolder*)DerivativeListHolder_ptr)->setAccumulate(*add); +} + void* DerivativeListHolder_initialize() { return new DerivativeListHolder(); diff --git a/sammy/src/salmon/interface/cpp/DerivativeListHolderInterface.h b/sammy/src/salmon/interface/cpp/DerivativeListHolderInterface.h index 4a7cdeaf8..e3e70dc98 100644 --- a/sammy/src/salmon/interface/cpp/DerivativeListHolderInterface.h +++ b/sammy/src/salmon/interface/cpp/DerivativeListHolderInterface.h @@ -2,7 +2,7 @@ * This file has been dynamically generated by Class Interface Xml (CIX) * DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION * If changes need to occur, modify the appropriate CIX xml file -* Date Generated: Tue Apr 06 18:09:02 EDT 2021 +* Date Generated: Fri Aug 20 14:10:56 EDT 2021 * If any issues are experiences with this generated file that cannot be fixed * with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov */ @@ -16,11 +16,15 @@ extern "C" { int DerivativeListHolder_getLength(void * DerivativeListHolder_ptr,bool * current); void DerivativeListHolder_addData(void * DerivativeListHolder_ptr,int * row,int * col,int * iso,double * val,bool * current); double DerivativeListHolder_getData(void * DerivativeListHolder_ptr,int * row,int * col,int * iso,bool * current); +void DerivativeListHolder_reserveColumns(void * DerivativeListHolder_ptr,int * row,int * maxCol,bool * current); +void DerivativeListHolder_reserve(void * DerivativeListHolder_ptr,int * maxRow,int * maxCol,bool * current); void DerivativeListHolder_clear(void * DerivativeListHolder_ptr); void DerivativeListHolder_nullify(void * DerivativeListHolder_ptr,bool * current); void DerivativeListHolder_switchGrid(void * DerivativeListHolder_ptr); int DerivativeListHolder_getIsotopeForShared(void * DerivativeListHolder_ptr,int * col); void DerivativeListHolder_addSharedColumn(void * DerivativeListHolder_ptr,int * col,int * iso); +void DerivativeListHolder_setNotSetReturnsZero(void * DerivativeListHolder_ptr,bool * empty); +void DerivativeListHolder_setAccumulate(void * DerivativeListHolder_ptr,bool * add); void* DerivativeListHolder_initialize(); void DerivativeListHolder_destroy(void * DerivativeListHolder_ptr); #ifdef __cplusplus diff --git a/sammy/src/salmon/interface/cpp/GridDataInterface.cpp b/sammy/src/salmon/interface/cpp/GridDataInterface.cpp index 155b538e0..2c18d8328 100644 --- a/sammy/src/salmon/interface/cpp/GridDataInterface.cpp +++ b/sammy/src/salmon/interface/cpp/GridDataInterface.cpp @@ -2,7 +2,7 @@ * This file has been dynamically generated by Class Interface Xml (CIX) * DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION * If changes need to occur, modify the appropriate CIX xml file -* Date Generated: Thu Jan 07 15:23:53 EST 2021 +* Date Generated: Fri Aug 20 09:26:07 EDT 2021 * If any issues are experiences with this generated file that cannot be fixed * with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov */ @@ -59,6 +59,16 @@ void GridData_clearGrid(void * GridData_ptr) ((GridData*)GridData_ptr)->clearGrid(); } +void GridData_nullify(void * GridData_ptr) +{ + ((GridData*)GridData_ptr)->nullify(); +} + +void GridData_reserve(void * GridData_ptr,int * row,int * col) +{ + ((GridData*)GridData_ptr)->reserve(*row,*col); +} + void GridData_addImplicitDerivs(void * GridData_ptr,GridData* grid) { std::unique_ptr<GridData> gridPtr(grid); diff --git a/sammy/src/salmon/interface/cpp/GridDataInterface.h b/sammy/src/salmon/interface/cpp/GridDataInterface.h index 00f0bac8e..a85b0ebea 100644 --- a/sammy/src/salmon/interface/cpp/GridDataInterface.h +++ b/sammy/src/salmon/interface/cpp/GridDataInterface.h @@ -2,7 +2,7 @@ * This file has been dynamically generated by Class Interface Xml (CIX) * DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION * If changes need to occur, modify the appropriate CIX xml file -* Date Generated: Thu Jan 07 15:23:53 EST 2021 +* Date Generated: Fri Aug 20 09:26:07 EDT 2021 * If any issues are experiences with this generated file that cannot be fixed * with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov */ @@ -24,6 +24,8 @@ void GridData_setRowOffset(void * GridData_ptr,int * offset); int GridData_getRowMax(void * GridData_ptr); void GridData_setRowMax(void * GridData_ptr,int * offset); void GridData_clearGrid(void * GridData_ptr); +void GridData_nullify(void * GridData_ptr); +void GridData_reserve(void * GridData_ptr,int * row,int * col); void GridData_addImplicitDerivs(void * GridData_ptr,GridData* grid); void* GridData_getImplicitDerivs(void * GridData_ptr); void GridData_addImplicitParCov(void * GridData_ptr,endf::ResonanceCovariance* cov); diff --git a/sammy/src/salmon/interface/fortran/DerivativeListHolder_I.f90 b/sammy/src/salmon/interface/fortran/DerivativeListHolder_I.f90 index 78cd7d2a7..550254f21 100644 --- a/sammy/src/salmon/interface/fortran/DerivativeListHolder_I.f90 +++ b/sammy/src/salmon/interface/fortran/DerivativeListHolder_I.f90 @@ -2,7 +2,7 @@ !! This file has been dynamically generated by Class Interface Xml (CIX) !! DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION !! If changes need to occur, modify the appropriate CIX xml file -!! Date Generated: Tue Apr 06 18:09:02 EDT 2021 +!! Date Generated: Fri Aug 20 14:10:56 EDT 2021 !! If any issues are experiences with this generated file that cannot be fixed !! with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov !!/ @@ -34,6 +34,22 @@ real(C_DOUBLE) function f_DerivativeListHolder_getData(DerivativeListHolder_ptr, integer(C_INT) :: iso; logical(C_BOOL) :: current; end function +subroutine f_DerivativeListHolder_reserveColumns(DerivativeListHolder_ptr, row,maxCol,current ) BIND(C,name="DerivativeListHolder_reserveColumns") + use,intrinsic :: ISO_C_BINDING + implicit none + type(C_PTR), value :: DerivativeListHolder_ptr; + integer(C_INT) :: row; + integer(C_INT) :: maxCol; + logical(C_BOOL) :: current; +end subroutine +subroutine f_DerivativeListHolder_reserve(DerivativeListHolder_ptr, maxRow,maxCol,current ) BIND(C,name="DerivativeListHolder_reserve") + use,intrinsic :: ISO_C_BINDING + implicit none + type(C_PTR), value :: DerivativeListHolder_ptr; + integer(C_INT) :: maxRow; + integer(C_INT) :: maxCol; + logical(C_BOOL) :: current; +end subroutine subroutine f_DerivativeListHolder_clear(DerivativeListHolder_ptr ) BIND(C,name="DerivativeListHolder_clear") use,intrinsic :: ISO_C_BINDING implicit none @@ -63,6 +79,18 @@ subroutine f_DerivativeListHolder_addSharedColumn(DerivativeListHolder_ptr, col, integer(C_INT) :: col; integer(C_INT) :: iso; end subroutine +subroutine f_DerivativeListHolder_setNotSetReturnsZero(DerivativeListHolder_ptr, empty ) BIND(C,name="DerivativeListHolder_setNotSetReturnsZero") + use,intrinsic :: ISO_C_BINDING + implicit none + type(C_PTR), value :: DerivativeListHolder_ptr; + logical(C_BOOL) :: empty; +end subroutine +subroutine f_DerivativeListHolder_setAccumulate(DerivativeListHolder_ptr, add ) BIND(C,name="DerivativeListHolder_setAccumulate") + use,intrinsic :: ISO_C_BINDING + implicit none + type(C_PTR), value :: DerivativeListHolder_ptr; + logical(C_BOOL) :: add; +end subroutine type(C_PTR) function f_DerivativeListHolder_initialize( )BIND(C,name="DerivativeListHolder_initialize") use,intrinsic :: ISO_C_BINDING implicit none diff --git a/sammy/src/salmon/interface/fortran/DerivativeListHolder_M.f90 b/sammy/src/salmon/interface/fortran/DerivativeListHolder_M.f90 index 25d69068f..6cb8811c8 100644 --- a/sammy/src/salmon/interface/fortran/DerivativeListHolder_M.f90 +++ b/sammy/src/salmon/interface/fortran/DerivativeListHolder_M.f90 @@ -18,12 +18,18 @@ type DerivativeListHolder procedure, pass(this) :: addDataOld => DerivativeListHolder_addDataOld procedure, pass(this) :: getData => DerivativeListHolder_getData procedure, pass(this) :: getDataOld => DerivativeListHolder_getDataOld + procedure, pass(this) :: reserveColumns => DerivativeListHolder_reserveColumns + procedure, pass(this) :: reserveColumnsOld => DerivativeListHolder_reserveColumnsOld + procedure, pass(this) :: reserve => DerivativeListHolder_reserve + procedure, pass(this) :: reserveOld => DerivativeListHolder_reserveOld procedure, pass(this) :: clear => DerivativeListHolder_clear procedure, pass(this) :: nullify => DerivativeListHolder_nullify procedure, pass(this) :: nullifyOld => DerivativeListHolder_nullifyOld procedure, pass(this) :: switchGrid => DerivativeListHolder_switchGrid procedure, pass(this) :: getIsotopeForShared => DerivativeListHolder_getIsotopeForShared procedure, pass(this) :: addSharedColumn => DerivativeListHolder_addSharedColumn + procedure, pass(this) :: setNotSetReturnsZero => DerivativeListHolder_setNotSetReturnsZero + procedure, pass(this) :: setAccumulate => DerivativeListHolder_setAccumulate procedure, pass(this) :: initialize => DerivativeListHolder_initialize procedure, pass(this) :: destroy => DerivativeListHolder_destroy end type DerivativeListHolder @@ -88,6 +94,42 @@ function DerivativeListHolder_getDataOld(this, row, col, iso) result(result2Retu current = .false. result2Return=f_DerivativeListHolder_getData(this%instance_ptr, row-1,col,iso-1,current) end function +subroutine DerivativeListHolder_reserveColumns(this, row, maxCol) + implicit none + class(DerivativeListHolder)::this + integer(C_INT)::row + integer(C_INT)::maxCol + logical(C_BOOL)::current + current=.true. + call f_DerivativeListHolder_reserveColumns(this%instance_ptr, row-1,maxCol+1,current) +end subroutine +subroutine DerivativeListHolder_reserveColumnsOld(this, row, maxCol) + implicit none + class(DerivativeListHolder)::this + integer(C_INT)::row + integer(C_INT)::maxCol + logical(C_BOOL)::current + current=.false. + call f_DerivativeListHolder_reserveColumns(this%instance_ptr, row-1,maxCol+1,current) +end subroutine +subroutine DerivativeListHolder_reserve(this, maxRow, maxCol) + implicit none + class(DerivativeListHolder)::this + integer(C_INT)::maxRow + integer(C_INT)::maxCol + logical(C_BOOL)::current + current=.true. + call f_DerivativeListHolder_reserve(this%instance_ptr, maxRow,maxCol+1,current) +end subroutine +subroutine DerivativeListHolder_reserveOld(this, maxRow, maxCol) + implicit none + class(DerivativeListHolder)::this + integer(C_INT)::maxRow + integer(C_INT)::maxCol + logical(C_BOOL)::current + current=.false. + call f_DerivativeListHolder_reserve(this%instance_ptr, maxRow,maxCol+1,current) +end subroutine subroutine DerivativeListHolder_clear(this) implicit none class(DerivativeListHolder)::this @@ -126,6 +168,18 @@ subroutine DerivativeListHolder_addSharedColumn(this, col, iso) integer(C_INT)::iso call f_DerivativeListHolder_addSharedColumn(this%instance_ptr, col,iso-1) end subroutine +subroutine DerivativeListHolder_setNotSetReturnsZero(this, empty) + implicit none + class(DerivativeListHolder)::this + logical(C_BOOL)::empty + call f_DerivativeListHolder_setNotSetReturnsZero(this%instance_ptr, empty) +end subroutine +subroutine DerivativeListHolder_setAccumulate(this, add) + implicit none + class(DerivativeListHolder)::this + logical(C_BOOL)::add + call f_DerivativeListHolder_setAccumulate(this%instance_ptr, add) +end subroutine subroutine DerivativeListHolder_initialize(this) implicit none class(DerivativeListHolder) :: this diff --git a/sammy/src/salmon/interface/fortran/GridData_I.f90 b/sammy/src/salmon/interface/fortran/GridData_I.f90 index 3b6109be7..96ed8dc0c 100644 --- a/sammy/src/salmon/interface/fortran/GridData_I.f90 +++ b/sammy/src/salmon/interface/fortran/GridData_I.f90 @@ -2,7 +2,7 @@ !! This file has been dynamically generated by Class Interface Xml (CIX) !! DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION !! If changes need to occur, modify the appropriate CIX xml file -!! Date Generated: Thu Jan 07 15:23:53 EST 2021 +!! Date Generated: Fri Aug 20 09:26:07 EDT 2021 !! If any issues are experiences with this generated file that cannot be fixed !! with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov !!/ @@ -67,6 +67,18 @@ subroutine f_GridData_clearGrid(GridData_ptr ) BIND(C,name="GridData_clearGrid") implicit none type(C_PTR), value :: GridData_ptr; end subroutine +subroutine f_GridData_nullify(GridData_ptr ) BIND(C,name="GridData_nullify") + use,intrinsic :: ISO_C_BINDING + implicit none + type(C_PTR), value :: GridData_ptr; +end subroutine +subroutine f_GridData_reserve(GridData_ptr, row,col ) BIND(C,name="GridData_reserve") + use,intrinsic :: ISO_C_BINDING + implicit none + type(C_PTR), value :: GridData_ptr; + integer(C_INT) :: row; + integer(C_INT) :: col; +end subroutine subroutine f_GridData_addImplicitDerivs(GridData_ptr, grid ) BIND(C,name="GridData_addImplicitDerivs") use,intrinsic :: ISO_C_BINDING implicit none diff --git a/sammy/src/salmon/interface/fortran/GridData_M.f90 b/sammy/src/salmon/interface/fortran/GridData_M.f90 index 040ca2b04..5e0d99d5c 100644 --- a/sammy/src/salmon/interface/fortran/GridData_M.f90 +++ b/sammy/src/salmon/interface/fortran/GridData_M.f90 @@ -2,7 +2,7 @@ !! This file has been dynamically generated by Class Interface Xml (CIX) !! DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION !! If changes need to occur, modify the appropriate CIX xml file -!! Date Generated: Thu Jan 07 15:23:53 EST 2021 +!! Date Generated: Fri Aug 20 09:26:07 EDT 2021 !! If any issues are experiences with this generated file that cannot be fixed !! with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov !!/ @@ -23,6 +23,8 @@ type GridData procedure, pass(this) :: getRowMax => GridData_getRowMax procedure, pass(this) :: setRowMax => GridData_setRowMax procedure, pass(this) :: clearGrid => GridData_clearGrid + procedure, pass(this) :: nullify => GridData_nullify + procedure, pass(this) :: reserve => GridData_reserve procedure, pass(this) :: addImplicitDerivs => GridData_addImplicitDerivs procedure, pass(this) :: getImplicitDerivs => GridData_getImplicitDerivs procedure, pass(this) :: addImplicitParCov => GridData_addImplicitParCov @@ -109,6 +111,18 @@ subroutine GridData_clearGrid(this) class(GridData)::this call f_GridData_clearGrid(this%instance_ptr) end subroutine +subroutine GridData_nullify(this) + implicit none + class(GridData)::this + call f_GridData_nullify(this%instance_ptr) +end subroutine +subroutine GridData_reserve(this, row, col) + implicit none + class(GridData)::this + integer(C_INT)::row + integer(C_INT)::col + call f_GridData_reserve(this%instance_ptr, row,col) +end subroutine subroutine GridData_addImplicitDerivs(this, grid) implicit none class(GridData)::this diff --git a/sammy/src/salmon/tests/DerivativeListTest.cpp b/sammy/src/salmon/tests/DerivativeListTest.cpp index 5759268e9..a32f41e7d 100644 --- a/sammy/src/salmon/tests/DerivativeListTest.cpp +++ b/sammy/src/salmon/tests/DerivativeListTest.cpp @@ -87,7 +87,6 @@ TEST(LinkedGridDataList, basic_data){ ASSERT_TRUE(list.isSharedColumn(i)); ASSERT_TRUE(!list.isSharedColumn(i+1)); } - ASSERT_EQ(0, list.getLength()); // clear clears everything list.clear(); @@ -186,6 +185,82 @@ TEST(LinkedGridDataList, copySharedInfo){ ASSERT_NO_THROW(listCopy.addData(2, 5, 2, 2)); } + +TEST(LinkedGridDataList, setToZeroOnEmpty){ + sammy::DerivativeList list; + + // shared columns are added and repored correctly + int j = 0; + for (int i = 0; i < 10; i+= 2){ + list.addSharedColumn(i, j); + j++; + if (j > 1) j = 0; + } + + // add three grids + list.addGrid(); + list.addGrid(); + + j = 0; + for (int i = 8; i >= 0; i-= 2){ + for( int r = 5; r < 10; r++){ + list.addData(r, i, j, i+r+20); + list.addData(r, i+1, 0, i+r+30); + list.addData(r, i+1, 1, i+r+40); + } + j++; + if (j > 1) j = 0; + } + list.addData(5, 40, 0, 40); + list.addData(5, 35, 0, 35); + ASSERT_NEAR(list.getData(5, 40,0), 40, 1e-5); + ASSERT_NEAR(list.getData(5, 35,0), 35, 1e-5); + + + + j = 0; + for (int i = 0; i < 10; i+= 2){ + for( int r = 5; r < 10; r++){ + ASSERT_NEAR(list.getData(r, i, j), i+r+20, 1e-5); + ASSERT_NEAR(list.getData(r, i+1, 0), i+r+30, 1e-5); + ASSERT_NEAR(list.getData(r, i+1, 1), i+r+40, 1e-5); + } + for( int r = 0; r < 5; r++){ + ASSERT_THROW(list.getData(r, i, j), std::runtime_error); + ASSERT_THROW(list.getData(r,i+1,0), std::runtime_error); + ASSERT_THROW(list.getData(r,i+1,1), std::runtime_error); + } + j++; + if (j > 1) j = 0; + } + + // now allow zero + list.setNotSetReturnsZero(true); + j = 0; + for (int i = 0; i < 10; i+= 2){ + for( int r = 0; r < 5; r++){ + ASSERT_EQ(list.getData(r, i, j), 0.0); + ASSERT_EQ(list.getData(r,i+1,0), 0.0); + ASSERT_EQ(list.getData(r,i+1,1), 0.0); + } + j++; + if (j > 1) j = 0; + } + + // nullify and test + list.nullify(); + j = 0; + for (int i = 0; i < 10; i+= 2){ + for( int r = 5; r < 10; r++){ + ASSERT_EQ(list.getData(r, i, j), 0.0); + ASSERT_EQ(list.getData(r,i+1,0), 0.0); + ASSERT_EQ(list.getData(r,i+1,1), 0.0); + } + j++; + if (j > 1) j = 0; + } +} + TEST(LinkedGridDataListHolder, updateShared){ sammy::DerivativeListHolder list; @@ -394,7 +469,7 @@ TEST(LinkedGridDataList, DerivativeListHolder){ ASSERT_TRUE(list.isSharedColumn(i)); ASSERT_TRUE(!list.isSharedColumn(i+1)); } - ASSERT_EQ(0, list.getLength()); + ASSERT_EQ(5, list.getLength()); ASSERT_EQ(5, list.getLength(false)); list.nullify(false); @@ -402,8 +477,8 @@ TEST(LinkedGridDataList, DerivativeListHolder){ ASSERT_TRUE(list.isSharedColumn(i)); ASSERT_TRUE(!list.isSharedColumn(i+1)); } - ASSERT_EQ(0, list.getLength()); - ASSERT_EQ(0, list.getLength(false)); + ASSERT_EQ(5, list.getLength()); + ASSERT_EQ(5, list.getLength(false)); // clear list.clear(); @@ -415,5 +490,76 @@ TEST(LinkedGridDataList, DerivativeListHolder){ ASSERT_EQ(0, list.getLength(false)); } +TEST(LinkedGridDataListHolder, setToZeroOnEmpty){ + sammy::DerivativeListHolder list; + + // shared columns are added and repored correctly + int j = 0; + for (int i = 0; i < 10; i+= 2){ + list.addSharedColumn(i, j); + j++; + if (j > 1) j = 0; + } + + // add three grids + list.addGrid(); + list.addGrid(); + + for (int ig = 0; ig < 1; ig++){ + j = 0; + for (int i = 0; i < 10; i+= 2){ + for( int r = 5; r < 10; r++){ + list.addData(r, i, j, i+r+20); + list.addData(r, i+1, 0, i+r+30); + list.addData(r, i+1, 1, i+r+40); + } + j++; + if (j > 1) j = 0; + } + + if (ig == 0){ + j = 0; + for (int i = 0; i < 10; i+= 2){ + for( int r = 0; r < 5; r++){ + ASSERT_THROW(list.getData(r, i, j), std::runtime_error); + ASSERT_THROW(list.getData(r,i+1,0), std::runtime_error); + ASSERT_THROW(list.getData(r,i+1,1), std::runtime_error); + } + j++; + if (j > 1) j = 0; + } + // now allow zero + list.setNotSetReturnsZero(true); + } + + j = 0; + for (int i = 0; i < 10; i+= 2){ + for( int r = 0; r < 5; r++){ + ASSERT_EQ(list.getData(r, i, j), 0.0); + ASSERT_EQ(list.getData(r,i+1,0), 0.0); + ASSERT_EQ(list.getData(r,i+1,1), 0.0); + } + j++; + if (j > 1) j = 0; + } + + // nullify and test + list.nullify(); + j = 0; + for (int i = 0; i < 10; i+= 2){ + for( int r = 5; r < 10; r++){ + ASSERT_EQ(list.getData(r, i, j), 0.0); + ASSERT_EQ(list.getData(r,i+1,0), 0.0); + ASSERT_EQ(list.getData(r,i+1,1), 0.0); + } + j++; + if (j > 1) j = 0; + } + + list.switchGrid(); + } +} + + diff --git a/sammy/src/salmon/tests/GridDataTest.cpp b/sammy/src/salmon/tests/GridDataTest.cpp index f063f7d41..1fdcbfd33 100644 --- a/sammy/src/salmon/tests/GridDataTest.cpp +++ b/sammy/src/salmon/tests/GridDataTest.cpp @@ -46,6 +46,14 @@ TEST(GridDataSetter,gridData){ } } + // now nullify + grid.nullify(); + for ( int i = 10; i < 20; i++ ){ + for ( int j = 5; j < 20; j++){ + ASSERT_NEAR(0.0, grid.getData(i,j), 1e-3); + } + } + // clear all data grid.clearGrid(); for ( int i = 10; i < 20; i++ ){ @@ -55,6 +63,27 @@ TEST(GridDataSetter,gridData){ } + // and a reserve + sammy::GridData gridReserve; + gridReserve.reserve(20, 2); + for ( int i = 0; i < 20; i++ ){ + for ( int j = 0; j < 2; j++){ + ASSERT_NEAR(0.0, gridReserve.getData(i,j), 1e-3); + gridReserve.addData(i, j, 1.0); + } + } + gridReserve.reserve(25, 3); + for ( int i = 0; i < 20; i++ ){ + for ( int j = 0; j < 2; j++){ + ASSERT_NEAR(1.0, gridReserve.getData(i,j), 1e-3); + } + ASSERT_NEAR(0.0, gridReserve.getData(i,2), 1e-3); + } + for ( int i = 20; i < 25; i++ ){ + for ( int j = 0; j < 3; j++){ + ASSERT_NEAR(0.0, gridReserve.getData(i,j), 1e-3); + } + } } @@ -98,6 +127,32 @@ TEST(DataInfo,gridData){ ASSERT_THROW(grid.addData(1, -5, 4.4), std::runtime_error); } +TEST(notSetReturnsZero,gridData){ + sammy::GridData grid; + + ASSERT_THROW(grid.getData(5,7), std::runtime_error); + grid.setNotSetReturnsZero(true); + ASSERT_EQ(0.0, grid.getData(5,7)); + + grid.addData(5,7, 6.0); + ASSERT_NEAR(6.0, grid.getData(5,7), 1e-5); + grid.setNotSetReturnsZero(false); + ASSERT_NEAR(6.0, grid.getData(5,7), 1e-5); + + ASSERT_THROW(grid.getData(5,8), std::runtime_error); + + // and a copy + grid.setNotSetReturnsZero(true); + sammy::GridData gridCopy(grid); + ASSERT_EQ(0.0, grid.getData(5,8)); + + grid.setNotSetReturnsZero(false); + sammy::GridData gridCopy2(grid); + ASSERT_THROW(grid.getData(5,8), std::runtime_error); +} + + + TEST(GridList, gridDataList){ sammy::GridDataList list; @@ -196,4 +251,31 @@ TEST(GridData, getImplicitParCov){ checkImplictParCov(gridDataCopy.getImplicitParCov()); } +TEST(notSetReturnsZero,gridDataList){ + sammy::GridDataList gridList; + + std::unique_ptr<sammy::GridData> grid = nemesis::make_unique<sammy::GridData>(); + grid->setNotSetReturnsZero(true); + gridList.addGrid(grid); + + ASSERT_THROW(gridList.getGrid(0)->getData(5,7), std::runtime_error); + + gridList.setNotSetReturnsZero(true); + ASSERT_EQ(0.0, gridList.getGrid(0)->getData(5,7)); + + grid = nemesis::make_unique<sammy::GridData>(); + grid->setNotSetReturnsZero(false); + gridList.addGrid(grid); + ASSERT_EQ(0.0, gridList.getGrid(1)->getData(5,7)); + + gridList.setNotSetReturnsZero(false); + ASSERT_THROW(gridList.getGrid(0)->getData(5,7), std::runtime_error); + ASSERT_THROW(gridList.getGrid(1)->getData(5,7), std::runtime_error); + + gridList.setNotSetReturnsZero(true); + ASSERT_EQ(0.0, gridList.getGrid(0)->getData(5,7)); + ASSERT_EQ(0.0, gridList.getGrid(1)->getData(5,7)); +} + + diff --git a/sammy/src/smc/msmc0.f b/sammy/src/smc/msmc0.f index c40bb99ea..8bb96e973 100644 --- a/sammy/src/smc/msmc0.f +++ b/sammy/src/smc/msmc0.f @@ -51,7 +51,7 @@ C *** Construct the energy scale I2 = 1 + (Eh-El)/De if (expData%getLength().gt.0) then call expData%getGrid(grid, 1) - call grid%clearGrid() + call grid%nullify() else call grid%initialize() call expData%addGrid(grid) @@ -72,7 +72,7 @@ C *** Read the energy scale from a file 11100 FORMAT (I6) if (expData%getLength().gt.0) then call expData%getGrid(grid, 1) - call grid%clearGrid() + call grid%nullify() else call grid%initialize() call expData%addGrid(grid) @@ -248,7 +248,7 @@ C *** set energy scale if needed IF (Kkkuse.EQ.1) THEN if (expData%getLength().gt.0) then call expData%getGrid(grid, 1) - call grid%clearGrid() + call grid%nullify() else call grid%initialize() call expData%addGrid(grid) diff --git a/sammy/src/ssm/mssm18.f90 b/sammy/src/ssm/mssm18.f90 index 191454a15..dd3ebdb21 100644 --- a/sammy/src/ssm/mssm18.f90 +++ b/sammy/src/ssm/mssm18.f90 @@ -332,6 +332,7 @@ module ssm_18_m Kk = Kk + 1 IF (dABS(auxGrid%getEnergy(K, expData)-Em).LT.De) THEN val = derivs%getData(Kk, 0, 1) + call derivs%reserveColumns(Jj, Ndasig+Ndbsig) call derivs%addData(Jj, 0, 1, val) Theory(Jj) = val IF (Ndasig.GT.0) THEN @@ -360,6 +361,7 @@ module ssm_18_m A = A/D B = B/D val = derivs%getData(Kk, 0, 1) * B + derivs%getData(Kk-1, 0, 1) * A + call derivs%reserveColumns(Jj, Ndasig+Ndbsig) call derivs%addData(Jj, 0, 1, val) Theory(Jj) = val IF (Ndasig.GT.0) THEN diff --git a/sammy/src/the/DerivativeHandler.cpp b/sammy/src/the/DerivativeHandler.cpp index 93b07b700..cdc52c745 100644 --- a/sammy/src/the/DerivativeHandler.cpp +++ b/sammy/src/the/DerivativeHandler.cpp @@ -1,8 +1,9 @@ #include "DerivativeHandler.h" namespace sammy{ - void DerivativeHandler::setUpList(SammyRMatrixParameters & params, int num){ - int numIso = params.getNumIso(); + + void DerivativeHandler::setUpList(SammyRMatrixParameters & params, AdjustedRadiusData & radii, int num){ + int numIso = params.getNumIso(); if (numIso < num) numIso = num; if (num == 1) numIso = 1; // numIso is sometimes redefined - why oh why if (numIso <= 1) { @@ -12,29 +13,98 @@ namespace sammy{ // add indices for resonance parameters // they are shared across isotopes - for ( int ires = 0; ires < params.getNumResonances(); ires++){ + for ( int ires = 0; ires < params.getNumResonances(); ires++){ sammy::SammyResonanceInfo * resInfo = params.getResonanceInfo(ires); - if (!resInfo->getIncludeInCalc()) continue; - if (!resInfo->hasAnyVariedParams()) continue; + bool includeInCalc = true; + + if (!resInfo->getIncludeInCalc()) includeInCalc = false; + if (!resInfo->hasAnyVariedParams()) includeInCalc = false; + + int igr = resInfo->getSpinGroupIndex(); + sammy::SammySpinGroupInfo *spinInfo = params.getSpinGroupInfo(igr); + if (!spinInfo->getIncludeInCalc()) includeInCalc = false; - int igr = resInfo->getSpinGroupIndex(); - sammy::SammySpinGroupInfo *spinInfo = params.getSpinGroupInfo(igr); - if (!spinInfo->getIncludeInCalc()) continue; + int iso = spinInfo->getIsotopeIndex(); + + // this allows legacy SAMMY to work. This shouldn't have been flagged + if (!includeInCalc) iso = 0; + + int iflGam = spinInfo->getGammaWidthIndex(); + if (iflGam > 0){ + addSharedColumn(iflGam, iso); + } int ifl = resInfo->getEnergyFitOption(); if (ifl > 0){ - addSharedColumn(ifl, spinInfo->getIsotopeIndex()); + addSharedColumn(ifl, iso); } int ichan = spinInfo->getAllChannels(); for (int ic = 0; ic < ichan; ic++){ ifl = resInfo->getChannelFitOption(ic); + if (ifl == iflGam) continue; + if (ifl > 0){ + addSharedColumn(ifl, iso); + } + } + } + + // r-external + for (int iext = 0; iext < params.getNumRext(); iext++){ + sammy::SammyRExternalInfo * resInfo = params.getRextInfo(iext); + int igr = resInfo->getGrp(); + sammy::SammySpinGroupInfo * spinInfo = params.getSpinGroupInfo(igr); + int iso = spinInfo->getIsotopeIndex(); + + for ( int i = 1; i <= 7; i++){ + int ifl = resInfo->getIflSammyIndex(i); if (ifl > 0){ - addSharedColumn(ifl, spinInfo->getIsotopeIndex()); + addSharedColumn(ifl, iso); } } } + // radius info + int nn = radii.getNumRadInfo(); + for (int ir = 0; ir < nn; ir++){ + int iflTrue, iflEff; + iflTrue = radii.getTrueFitFlagByIndex(ir); + iflEff = radii.getEffFitFlagByIndex(ir); + if (radii.trueEqualToEff(ir)) iflTrue = 0; + if (radii.hasMatchingRadius()) iflTrue = 0; + + if (iflEff <= 0 && iflTrue <= 0) continue; + + + if (radii.getNumGroups(ir) == 0) continue; + int igr = radii.getGroupIndex(ir, 0); // take istope from first group + sammy::SammySpinGroupInfo * spinInfo = params.getSpinGroupInfo(igr); + int iso = spinInfo->getIsotopeIndex(); + + if (iflTrue > 0){ + addSharedColumn(iflTrue, iso); + } + if (iflEff > 0){ + addSharedColumn(iflEff, iso); + } + } + + // matching radius, add to first isotope + if (nn == 0 || radii.hasMatchingRadius()){ + int ifl = radii.matchFitFlag(); + if (ifl > 0){ + addSharedColumn(ifl, 0); + } + } + + // abundances + for (int iso = 0; iso < params.getNumIso(); iso++){ + int ifl = params.getIsoInfo(iso)->getFitOption(); + if (ifl > 0){ + addSharedColumn(ifl, iso); + } + } + while ( getGridNumber() < numIso) addGrid(); } @@ -54,6 +124,11 @@ namespace sammy{ } } + void DerivativeHandler::reserveColumnsNs(int row, int ns, int maxCol, bool current){ + int ipos = row * nnsig + ns; + reserveColumns(ipos, maxCol, current); + } + double DerivativeHandler::getSharedDataVal(int row, int col, bool current) const{ int posOur = getIsotopeForShared(col); if (posOur < 0) posOur = 0; @@ -120,7 +195,7 @@ namespace sammy{ } double DerivativeHandler::getDataNs(int row, int ns, int col, int iso, bool current) const{ - int ipos = row * nnsig + ns; + int ipos = row * nnsig + ns; return getData(ipos, col, iso, current); } } diff --git a/sammy/src/the/DerivativeHandler.h b/sammy/src/the/DerivativeHandler.h index b1d0e711f..7ee9d4c18 100644 --- a/sammy/src/the/DerivativeHandler.h +++ b/sammy/src/the/DerivativeHandler.h @@ -3,6 +3,7 @@ #include "salmon/DerivativeList.h" #include "endf/SammyRMatrixParameters.h" +#include "endf/AdjustedRadiusData.h" namespace sammy{ @@ -57,8 +58,8 @@ namespace sammy{ * * @param params the resonance parameter information * @param num the number of isotopes - */ - void setUpList(SammyRMatrixParameters & params, int num=0); + */ + void setUpList(SammyRMatrixParameters & params, AdjustedRadiusData & radii, int num=0); /** * Set the data for the indicated row to zero @@ -132,6 +133,16 @@ namespace sammy{ */ void setNnsig(int nnsig); + + /** + * Make sure the indicated row as at least maxCol data. + * If not set, additional columns are set to zero + * + * @param row the row for which to set the data + * @param maxCol the maximum number of columns + */ + void reserveColumnsNs(int row, int ns, int maxCol, bool current = true); + /** * Same as getSharedDataVal, except return the value for row=((row-1)*getNnnsig() +ns) * diff --git a/sammy/src/the/interface/cix/DerivativeHandler.cpp2f.xml b/sammy/src/the/interface/cix/DerivativeHandler.cpp2f.xml index 2119ed6fb..600de7fb9 100644 --- a/sammy/src/the/interface/cix/DerivativeHandler.cpp2f.xml +++ b/sammy/src/the/interface/cix/DerivativeHandler.cpp2f.xml @@ -5,6 +5,7 @@ <class name="DerivativeHandler" parent="DerivativeListHolder"> <method name="setUpList"> <param name="params" type="SammyRMatrixParameters"/> + <param name="radii" type="AdjustedRadiusData"/> <param name="num" type="int"/> </method> <method name="setToZero"> @@ -52,6 +53,12 @@ <param name="value" type="double"/> <param name="current" type="bool"/> </method> + <method name="reserveColumnsNs"> + <param name="row" type="int" offset="-1"/> + <param name="ns" type="int" offset="-1"/> + <param name="maxCol" type="int" offset="+1"/> + <param name="current" type="bool"/> + </method> <method name="addDataNs"> <param name="row" type="int" offset="-1"/> <param name="ns" type="int" offset="-1"/> diff --git a/sammy/src/the/interface/cpp/DerivativeHandlerInterface.cpp b/sammy/src/the/interface/cpp/DerivativeHandlerInterface.cpp index 054ace377..bf20c2439 100644 --- a/sammy/src/the/interface/cpp/DerivativeHandlerInterface.cpp +++ b/sammy/src/the/interface/cpp/DerivativeHandlerInterface.cpp @@ -2,16 +2,16 @@ * This file has been dynamically generated by Class Interface Xml (CIX) * DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION * If changes need to occur, modify the appropriate CIX xml file -* Date Generated: Mon Apr 26 18:06:30 EDT 2021 +* Date Generated: Fri Aug 20 14:16:52 EDT 2021 * If any issues are experiences with this generated file that cannot be fixed * with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov */ #include <string.h> #include "DerivativeHandlerInterface.h" using namespace sammy; -void DerivativeHandler_setUpList(void * DerivativeHandler_ptr,SammyRMatrixParameters * params,int * num) +void DerivativeHandler_setUpList(void * DerivativeHandler_ptr,SammyRMatrixParameters * params,AdjustedRadiusData * radii,int * num) { - ((DerivativeHandler*)DerivativeHandler_ptr)->setUpList(*params,*num); + ((DerivativeHandler*)DerivativeHandler_ptr)->setUpList(*params,*radii,*num); } void DerivativeHandler_setToZero(void * DerivativeHandler_ptr,int * row,int * maxIfl,bool * current) @@ -64,14 +64,19 @@ void DerivativeHandler_setSharedValNs(void * DerivativeHandler_ptr,int * row,int ((DerivativeHandler*)DerivativeHandler_ptr)->setSharedValNs(*row,*ns,*col,*value,*current); } -void DerivativeHandler_addDataNs(void * DerivativeHandler_ptr,int * row,int * ns,int * col,int * index,double * value,bool * current) +void DerivativeHandler_reserveColumnsNs(void * DerivativeHandler_ptr,int * row,int * ns,int * maxCol,bool * current) { - ((DerivativeHandler*)DerivativeHandler_ptr)->addDataNs(*row,*ns,*col,*index,*value,*current); + ((DerivativeHandler*)DerivativeHandler_ptr)->reserveColumnsNs(*row,*ns,*maxCol,*current); } -double DerivativeHandler_getDataNs(void * DerivativeHandler_ptr,int * row,int * ns,int * col,int * index,bool * current) +void DerivativeHandler_addDataNs(void * DerivativeHandler_ptr,int * row,int * ns,int * col,int * iso,double * value,bool * current) { - return ((DerivativeHandler*)DerivativeHandler_ptr)->getDataNs(*row,*ns,*col,*index,*current); + ((DerivativeHandler*)DerivativeHandler_ptr)->addDataNs(*row,*ns,*col,*iso,*value,*current); +} + +double DerivativeHandler_getDataNs(void * DerivativeHandler_ptr,int * row,int * ns,int * col,int * iso,bool * current) +{ + return ((DerivativeHandler*)DerivativeHandler_ptr)->getDataNs(*row,*ns,*col,*iso,*current); } void* DerivativeHandler_initialize() diff --git a/sammy/src/the/interface/cpp/DerivativeHandlerInterface.h b/sammy/src/the/interface/cpp/DerivativeHandlerInterface.h index 243e1e5f7..808fde4fc 100644 --- a/sammy/src/the/interface/cpp/DerivativeHandlerInterface.h +++ b/sammy/src/the/interface/cpp/DerivativeHandlerInterface.h @@ -2,7 +2,7 @@ * This file has been dynamically generated by Class Interface Xml (CIX) * DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION * If changes need to occur, modify the appropriate CIX xml file -* Date Generated: Mon Apr 26 18:06:30 EDT 2021 +* Date Generated: Fri Aug 20 14:16:51 EDT 2021 * If any issues are experiences with this generated file that cannot be fixed * with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov */ @@ -13,7 +13,7 @@ using namespace sammy; #ifdef __cplusplus extern "C" { #endif -void DerivativeHandler_setUpList(void * DerivativeHandler_ptr,SammyRMatrixParameters * params,int * num); +void DerivativeHandler_setUpList(void * DerivativeHandler_ptr,SammyRMatrixParameters * params,AdjustedRadiusData * radii,int * num); void DerivativeHandler_setToZero(void * DerivativeHandler_ptr,int * row,int * maxIfl,bool * current); double DerivativeHandler_getSharedDataVal(void * DerivativeHandler_ptr,int * row,int * col,bool * current); void DerivativeHandler_setSharedDataVal(void * DerivativeHandler_ptr,int * row,int * col,double * value,bool * current); @@ -24,8 +24,9 @@ int DerivativeHandler_getNnnsig(void * DerivativeHandler_ptr); void DerivativeHandler_setNnsig(void * DerivativeHandler_ptr,int * nnsig); double DerivativeHandler_getSharedValNs(void * DerivativeHandler_ptr,int * row,int * ns,int * col,bool * current); void DerivativeHandler_setSharedValNs(void * DerivativeHandler_ptr,int * row,int * ns,int * col,double * value,bool * current); -void DerivativeHandler_addDataNs(void * DerivativeHandler_ptr,int * row,int * ns,int * col,int * index,double * value,bool * current); -double DerivativeHandler_getDataNs(void * DerivativeHandler_ptr,int * row,int * ns,int * col,int * index,bool * current); +void DerivativeHandler_reserveColumnsNs(void * DerivativeHandler_ptr,int * row,int * ns,int * maxCol,bool * current); +void DerivativeHandler_addDataNs(void * DerivativeHandler_ptr,int * row,int * ns,int * col,int * iso,double * value,bool * current); +double DerivativeHandler_getDataNs(void * DerivativeHandler_ptr,int * row,int * ns,int * col,int * iso,bool * current); void* DerivativeHandler_initialize(); void DerivativeHandler_destroy(void * DerivativeHandler_ptr); #ifdef __cplusplus diff --git a/sammy/src/the/interface/fortran/DerivativeHandler_I.f90 b/sammy/src/the/interface/fortran/DerivativeHandler_I.f90 index 88e3a4467..00e30e462 100644 --- a/sammy/src/the/interface/fortran/DerivativeHandler_I.f90 +++ b/sammy/src/the/interface/fortran/DerivativeHandler_I.f90 @@ -2,18 +2,19 @@ !! This file has been dynamically generated by Class Interface Xml (CIX) !! DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION !! If changes need to occur, modify the appropriate CIX xml file -!! Date Generated: Mon Apr 26 18:06:30 EDT 2021 +!! Date Generated: Fri Aug 20 14:16:52 EDT 2021 !! If any issues are experiences with this generated file that cannot be fixed !! with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov !!/ module DerivativeHandler_I use, intrinsic :: ISO_C_BINDING interface -subroutine f_DerivativeHandler_setUpList(DerivativeHandler_ptr, params,num ) BIND(C,name="DerivativeHandler_setUpList") +subroutine f_DerivativeHandler_setUpList(DerivativeHandler_ptr, params,radii,num ) BIND(C,name="DerivativeHandler_setUpList") use,intrinsic :: ISO_C_BINDING implicit none type(C_PTR), value :: DerivativeHandler_ptr; type(C_PTR), value :: params; + type(C_PTR), value :: radii; integer(C_INT) :: num; end subroutine subroutine f_DerivativeHandler_setToZero(DerivativeHandler_ptr, row,maxIfl,current ) BIND(C,name="DerivativeHandler_setToZero") @@ -92,6 +93,15 @@ subroutine f_DerivativeHandler_setSharedValNs(DerivativeHandler_ptr, row,ns,col, real(C_DOUBLE) :: value; logical(C_BOOL) :: current; end subroutine +subroutine f_DerivativeHandler_reserveColumnsNs(DerivativeHandler_ptr, row,ns,maxCol,current ) BIND(C,name="DerivativeHandler_reserveColumnsNs") + use,intrinsic :: ISO_C_BINDING + implicit none + type(C_PTR), value :: DerivativeHandler_ptr; + integer(C_INT) :: row; + integer(C_INT) :: ns; + integer(C_INT) :: maxCol; + logical(C_BOOL) :: current; +end subroutine subroutine f_DerivativeHandler_addDataNs(DerivativeHandler_ptr, row,ns,col,iso,value,current ) BIND(C,name="DerivativeHandler_addDataNs") use,intrinsic :: ISO_C_BINDING implicit none diff --git a/sammy/src/the/interface/fortran/DerivativeHandler_M.f90 b/sammy/src/the/interface/fortran/DerivativeHandler_M.f90 index c220f1d8f..37c2a80b9 100644 --- a/sammy/src/the/interface/fortran/DerivativeHandler_M.f90 +++ b/sammy/src/the/interface/fortran/DerivativeHandler_M.f90 @@ -11,6 +11,7 @@ use, intrinsic :: ISO_C_BINDING use DerivativeHandler_I use DerivativeListHolder_M use SammyRMatrixParameters_M +use AdjustedRadiusData_M type , extends(DerivativeListHolder) :: DerivativeHandler contains procedure, pass(this) :: setUpList => DerivativeHandler_setUpList @@ -29,6 +30,8 @@ type , extends(DerivativeListHolder) :: DerivativeHandler procedure, pass(this) :: getSharedValNsOld => DerivativeHandler_getSharedValNsOld procedure, pass(this) :: setSharedValNs => DerivativeHandler_setSharedValNs procedure, pass(this) :: setSharedValNsOld => DerivativeHandler_setSharedValNsOld + procedure, pass(this) :: reserveColumnsNs => DerivativeHandler_reserveColumnsNs + procedure, pass(this) :: reserveColumnsNsOld => DerivativeHandler_reserveColumnsNsOld procedure, pass(this) :: addDataNs => DerivativeHandler_addDataNs procedure, pass(this) :: addDataNsOld => DerivativeHandler_addDataNsOld procedure, pass(this) :: getDataNs => DerivativeHandler_getDataNs @@ -40,12 +43,13 @@ type , extends(DerivativeListHolder) :: DerivativeHandler procedure, pass(this) :: destroy => DerivativeHandler_destroy end type DerivativeHandler contains -subroutine DerivativeHandler_setUpList(this, params, num) +subroutine DerivativeHandler_setUpList(this, params, radii, num) implicit none class(DerivativeHandler)::this class(SammyRMatrixParameters)::params + class(AdjustedRadiusData)::radii integer(C_INT)::num - call f_DerivativeHandler_setUpList(this%instance_ptr, params%instance_ptr,num) + call f_DerivativeHandler_setUpList(this%instance_ptr, params%instance_ptr, radii%instance_ptr, num) end subroutine subroutine DerivativeHandler_setToZero(this, row, maxIfl) implicit none @@ -187,6 +191,26 @@ subroutine DerivativeHandler_setSharedValNsOld(this, row, ns, col, value) current=.false. call f_DerivativeHandler_setSharedValNs(this%instance_ptr, row-1,ns-1,col,value,current) end subroutine +subroutine DerivativeHandler_reserveColumnsNs(this, row, ns, maxCol) + implicit none + class(DerivativeHandler)::this + integer(C_INT)::row + integer(C_INT)::ns + integer(C_INT)::maxCol + logical(C_BOOL)::current + current = .true. + call f_DerivativeHandler_reserveColumnsNs(this%instance_ptr, row-1,ns-1,maxCol+1,current) +end subroutine +subroutine DerivativeHandler_reserveColumnsNsOld(this, row, ns, maxCol) + implicit none + class(DerivativeHandler)::this + integer(C_INT)::row + integer(C_INT)::ns + integer(C_INT)::maxCol + logical(C_BOOL)::current + current = .false. + call f_DerivativeHandler_reserveColumnsNs(this%instance_ptr, row-1,ns-1,maxCol+1,current) +end subroutine subroutine DerivativeHandler_addDataNs(this, row, ns, col, iso, value) implicit none class(DerivativeHandler)::this @@ -264,18 +288,19 @@ subroutine DerivativeHandler_addCalculatedData(this, row, nnsig, ndasig, ndbsig, ourIso = iso if (ourIso.lt.0) ourIso = 1 do Jsig = 1,nnsig + call this%reserveColumnsNs(row, Jsig, Ndasig+Ndbsig) call this%addDataNs(row, Jsig, 0, ourIso, Sigx(Jsig)) DO Iipar=1,Ndasig isop = -1 ! indicate the shared value is to be set if (iso.ge.0) then isop = this%getIsotopeForShared(Iipar) if (isop.le.0) isop = 1 ! for one-isotope problems we don't use shared column - end if + end if if (iso.eq.isop) then call this%setSharedValNs(row, Jsig, Iipar, Dasigx(jsig,Iipar)) end if end do - DO Iipar=1,Ndbsig + DO Iipar=1,Ndbsig call this%addDataNs(row, Jsig, Iipar + Ndasig, ourIso, Dbsigx(Jsig,Iipar)) end do end do @@ -297,7 +322,7 @@ subroutine DerivativeHandler_getCalculatedData(this, row, nnsig, ndasig, ndbsig, DO Iipar=0,Ndasig + Ndbsig sharedIso = this%getIsotopeForShared(Iipar) if (sharedIso.le.0) sharedIso = 1 - if (iso.gt.0) sharedIso = iso + if (iso.gt.0) sharedIso = iso val = this%getDataNs(row, Jsig, Iipar, sharedIso) if (Iipar.eq.0) then Sigx(Jsig) = val @@ -325,7 +350,7 @@ subroutine DerivativeHandler_getCalculatedDataOld(this, row, nnsig, ndasig, ndbs DO Iipar=0,Ndasig + Ndbsig sharedIso = this%getIsotopeForShared(Iipar) if (sharedIso.le.0) sharedIso = 1 - if (iso.gt.0) sharedIso = iso + if (iso.gt.0) sharedIso = iso val = this%getDataNsOld(row, Jsig, Iipar, sharedIso) if (Iipar.eq.0) then Sigx(Jsig) = val diff --git a/sammy/src/xct/mxct02.f90 b/sammy/src/xct/mxct02.f90 index 22aa1f17a..96969d90f 100644 --- a/sammy/src/xct/mxct02.f90 +++ b/sammy/src/xct/mxct02.f90 @@ -99,7 +99,7 @@ module xct2_m call derivs%addSharedColumn(Iipar, 1) end do end if - call derivs%setUpList(resParData, Iq_Iso) + call derivs%setUpList(resParData, radFitFlags, Iq_Iso) IF (Iw.EQ.1.or.Ksitmp.gt.0) THEN call derivsSelf%nullify() @@ -109,7 +109,7 @@ module xct2_m call derivsSelf%addSharedColumn(Iipar, 1) end do end if - call derivsSelf%setUpList(resParData, Iq_Iso) + call derivsSelf%setUpList(resParData, radFitFlags, Iq_Iso) end if CALL Zero_Integer (Isopar, Ndasig) -- GitLab