diff --git a/.gitignore b/.gitignore index 3c6d13b..723e9f7 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,5 @@ *~ *tar.gz + +TEST/* \ No newline at end of file diff --git a/Makefile b/Makefile index a36357b..b634f10 100644 --- a/Makefile +++ b/Makefile @@ -22,9 +22,17 @@ SRCDIR = $(BASEDIR)src/ EXECSDIR = $(BASEDIR)src/EXECS/ BINDIR = $(BASEDIR)bin/ +# Set the compiler choice +#CXX = g++ +CXX = /usr/local/bin/g++ +MPICXX = mpicxx + + #COMPILE = g++ -Wall -I$(BASEDIR)include/ -L$(LIBDIR) -O3 -stdlib=libstdc++ #COMPILE = g++ -Wall -I$(BASEDIR)include/ -L$(LIBDIR) -O3 -COMPILE = g++ -I$(BASEDIR)include/ -L$(LIBDIR) -O3 -w -fopenmp +#COMPILE = g++ -I$(BASEDIR)include/ -L$(LIBDIR) -O3 -w -fopenmp +COMPILE = $(CXX) -I$(BASEDIR)include/ -L$(LIBDIR) -O3 -w -fopenmp +#COMPILE = clang++ -std=c++17 -Weverything -I$(BASEDIR)include/ -L$(LIBDIR) -O3 -w -fopenmp #COMPILE_MPI = mpicxx -I$(BASEDIR)include/ -L$(LIBDIR) -O3 -stdlib=libstdc++ COMPILE_MPI = mpicxx -I$(BASEDIR)include/ -L$(LIBDIR) -O3 -w -fopenmp #COMPILE_OMP = g++ -I$(BASEDIR)include/ -L$(LIBDIR) -O3 -w -fopenmp @@ -120,8 +128,8 @@ lib$(VERSION).a : $(Objects_ALL) $(COMPILE) $(EXECSDIR)Heis_DSF_GeneralState.cc -o $(BINDIR)Heis_DSF_GeneralState -l$(VERSION) $(COMPILE) $(EXECSDIR)Smoothen_Heis_DSF.cc -o $(BINDIR)Smoothen_Heis_DSF -l$(VERSION) $(COMPILE) $(EXECSDIR)XXZ_gpd_StagSz_h0.cc -o $(BINDIR)XXZ_gpd_StagSz_h0 -l$(VERSION) - $(COMPILE) $(EXECSDIR)ODSLF_DSF.cc -o $(BINDIR)ODSLF_DSF -l$(VERSION) - $(COMPILE) $(EXECSDIR)Smoothen_ODSLF_DSF.cc -o $(BINDIR)Smoothen_ODSLF_DSF -l$(VERSION) +# $(COMPILE) $(EXECSDIR)ODSLF_DSF.cc -o $(BINDIR)ODSLF_DSF -l$(VERSION) +# $(COMPILE) $(EXECSDIR)Smoothen_ODSLF_DSF.cc -o $(BINDIR)Smoothen_ODSLF_DSF -l$(VERSION) $(COMPILE) $(EXECSDIR)Check_RAW_File.cc -o $(BINDIR)Check_RAW_File -l$(VERSION) $(COMPILE) $(EXECSDIR)Analyze_RAW_File.cc -o $(BINDIR)Analyze_RAW_File -l$(VERSION) $(COMPILE) $(EXECSDIR)RAW_File_Stats.cc -o $(BINDIR)RAW_File_Stats -l$(VERSION) @@ -382,17 +390,9 @@ $(OBJDIR)General_Scan_Parallel.o : General_Scan_Parallel.cc $(Headers_all) $(OBJDIR)Particle_Hole_Excitation_Cost.o : Particle_Hole_Excitation_Cost.cc $(Headers_all) $(COMPILE) -c $< -o $@ -#$(OBJDIR)Offsets.o : Offsets.cc $(Headers_all) -# $(COMPILE) -c $< -o $@ - $(OBJDIR)Scan_Info.o : Scan_Info.cc $(Headers_all) $(COMPILE) -c $< -o $@ -#$(OBJDIR)Scan_State_List.o : Scan_State_List.cc $(Headers_all) -# $(COMPILE) -c $< -o $@ - -#$(OBJDIR)Scan_Thread_List.o : Scan_Thread_List.cc $(Headers_all) -#$(OBJDIR)Scan_Thread_Set.o : Scan_Thread_Set.cc $(Headers_all) $(OBJDIR)Scan_Thread_Data.o : Scan_Thread_Data.cc $(Headers_all) $(COMPILE) -c $< -o $@ @@ -451,14 +451,6 @@ $(OBJDIR)Young_Tableau.o : Young_Tableau.cc $(Headers_all) $(COMPILE) -c $< -o $@ -############################################ - -#executables : -# $(COMPILE) $(EXECSDIR)LiebLin_DSF.cc -o $(BINDIR)LiebLin_DSF -l$(VERSION) -# $(COMPILE) $(EXECSDIR)Smoothen_LiebLin_DSF.cc -o $(BINDIR)Smoothen_LiebLin_DSF -l$(VERSION) -# $(COMPILE) $(EXECSDIR)Heis_DSF.cc -o $(BINDIR)Heis_DSF -l$(VERSION) -# $(COMPILE) $(EXECSDIR)Smoothen_Heis_DSF.cc -o $(BINDIR)Smoothen_Heis_DSF -l$(VERSION) - ########################################### # Executables (parallel versions) diff --git a/README.md b/README.md index 04a0740..5e6cac8 100644 --- a/README.md +++ b/README.md @@ -1,16 +1,41 @@ # ABACUS -Copyright (c) J.-S. Caux +Copyright (c) J.-S. Caux. +ABACUS is a general set of algorithms for dealing with Bethe Ansatz-solvable systems including: + +* Lieb-Liniger (delta-)interacting bosons +* Heisenberg spin chains + * Isotropic antiferromagnet (XXX) + * Anisotropic gapless antiferromagnet (XXZ, 0 < Delta < 1) + * Anisotropic gapped antiferromagnet (XXZ_gpd, Delta > 1) +* Interacting fermions (Jordan-Wigner'd XXZ) + +The principal purposes of ABACUS are: +* to construct Bethe eigenstates +* to compute (dynamical) correlation functions of basic observables + + +## Installation +From the base directory, simply run +``` +$ make +``` +This will produce all executables, together with a library `ABACUS_[vn]` where vn is of the form [digit][character]. + + --- Previous versions ================= +Development of ABACUS started around 2004 and has known an embarrassingly large number of versions. ------------------------------------------------------------------- -## ABACUS++G_8 [started 2015 02 20] [published 2015 04 04] +Here under are some notes on major versions and overhauls of the codebase. + +#### ABACUS++G_8 +_[started 2015 02 20] [published 2015 04 04]_ Building up on ++G_7. @@ -19,54 +44,57 @@ Objectives: - facilitate use of OpenMP, by removing the recursive aspect of descending. DONE. - improve memory use. DONE. - Replace in-memory Scan_Thread_Set by on-disk Scan_Thread_Data. - Uses Scan_Thread structure. + Replace in-memory `Scan_Thread_Set` by on-disk `Scan_Thread_Data`. + Uses `Scan_Thread` structure. Features: - OpenMP is implemented. - Scanning on spin chains is now implemented. -- Scan_Info's Nfull now contains the sub-Hilbert space dimension (for spin chains; type == double) -- Added the .stat file in General_Scan, containing info about expected data value and actually obtained one -- Timing is more strictly enforced, using OpenMP wtime(). +- `Scan_Info`'s `Nfull` now contains the sub-Hilbert space dimension (for spin chains; `type == double`) +- Added the .stat file in `General_Scan`, containing info about expected data value and actually obtained one +- Timing is more strictly enforced, using OpenMP `wtime()`. Notes: - a first implementation using OpenMP is archived as ABACUS++G_8_v1.tar.gz package. ------------------------------------------------------------------- -ABACUS++G_7 [published 2015 01 17] +#### ABACUS++G_7 +_[published 2015 01 17]_ Building up on ++G_6. - Now using 15 types of descendents, enabling fixed iK scanning with increasing ph nrs. -- Scan_Thread_List replaced by Scan_Thread_Set (to optimize performance and threads storage) +- `Scan_Thread_List` replaced by `Scan_Thread_Set` (to optimize performance and threads storage) Works well for LiebLin, including at finite T. TODO: -- implement scanning for spin chains (scanning over bases still missing for generic AveragingState) +- implement scanning for spin chains (scanning over bases still missing for generic `AveragingState`) ------------------------------------------------------------------- -ABACUS++G_6 [started 2015 01] -Simple scanning: iK stepped up, iK stepped down (from ++G_5). +#### ABACUS++G_6 +_[started 2015 01]_ + +Simple scanning: `iK` stepped up, `iK` stepped down (from ++G_5). Fixed momentum scanning is thus implemented. -Version 6.0: [published 2015 01 16] +###### Version 6.0: _[published 2015 01 16]_ - uses 8 types of descendents - Works well for LiebLin, including finite T. -Version 6.1: proto version of ++G_7 +###### Version 6.1: proto version of ++G_7 - uses greater nr of descendents, 15, enforcing strict ph nr increase ------------------------------------------------------------------- -ABACUS++G_5 [started 2014 12 11] [abandoned 2015 01, except for version 5.0 (great for GS of LiebLin)] + +#### ABACUS++G_5 +_[started 2014 12 11] [abandoned 2015 01, except for version 5.0 (great for GS of LiebLin)]_ Fundamental rewriting of scanning protocol. Version 5.0 works very well for ground state correlations of Lieb-Liniger. Heis: still bugged. ------------------------------------------------------------------- -ABACUS++G_4 [started 2014 12 08] [abandoned 2014 12 11] + +#### ABACUS++G_4 +_[started 2014 12 08] [abandoned 2014 12 11]_ Nontrivial bug: descending with inner and outer skeletons does not always preserve number of p-h excitations. Fatal for large c in LiebLin. @@ -79,33 +107,32 @@ TODO: - For XXX: inclusion of infinite rapidities as genuine particle type ------------------------------------------------------------------- -ABACUS++G_3 [started 2014 11 10] [closed 2014 12 08] +#### ABACUS++G_3 +_[started 2014 11 10] [closed 2014 12 08]_ - Change of naming convention: LiebLin instead of Bose or 1DBG. -- For LiebLin: changed rapidity naming convention: now lambdaoc, to make rescaling \lambda/c explicit +- For LiebLin: changed rapidity naming convention: now `lambdaoc`, to make rescaling lambda/c explicit - Improved small c algorithm for LiebLin - Fixed momentum scanning implemented using inner and outer skeleton logic ------------------------------------------------------------------- -ABACUS++G_2 [closed 2014 11 10] -[2013 09 20] +#### ABACUS++G_2 +_[closed 2014 11 10] [2013 09 20]_ - Scanning is now automatically over remaining pair of excitations; facilitates fixed iK scans - Threads are now over states with NScan-2 particles, intermediate states then include all Nscan states fulfilling the momentum constraints - Scanning for spin chains implemented (general states; also parallel). Needs further testing. Seems to work for ground state. ------------------------------------------------------------------- -ABACUS++G_1.1 -[2013 09 06] +#### ABACUS++G_1.1 +_[2013 09 06]_ + Changed parallel version for Bose: - supercycle time is now an argument to functions -- split the 'Prepare', 'Run' and 'Wrapup' into separate executables -- added the 'filenamesuffix' argument to parallel functions +- split the `Prepare`, `Run` and `Wrapup` into separate executables +- added the `filenamesuffix` argument to parallel functions ------------------------------------------------------------------- -ABACUS++G_1 + +#### ABACUS++G_1 New labelling and descending implemented for arbitrary states with strings, e.g. spin chains. @@ -119,53 +146,82 @@ Include string deviation value in Bethe_State objects and in RAW files. Form fac Changed sum rules for LiebLin density-density: now using f-sumrule from iKmin to iKmax. -Earlier notes for ABACUS++T: -############################################################## +#### Earlier notes for ABACUS++T Notable changes from ABACUS++: -- The state labelling is now done with a single `string' label, detailing the base, nr of particle-hole excitations (as compared to a reference scan state) and Ix2 changes. The number of particle-hole excitations is therefore not limited anymore, and (for LiebLin) there is no upper momentum limit. See the src/State_Label.cc file for the implementation details. +- The state labelling is now done with a single `string` label, detailing the base, nr of particle-hole excitations (as compared to a reference scan state) and `Ix2` changes. The number of particle-hole excitations is therefore not limited anymore, and (for LiebLin) there is no upper momentum limit. See the `src/State_Label.cc` file for the implementation details. - The state scanning procedure (i.e. the descendents logic) is now implemented solely at the level of the quantum numbers, following a recursive logic in which excitations are created at the right, then left Fermi boundary and allowed to propagate deeper outside/inside as the tree is climbed. - Since the labelling and descending are now completely general, dynamical correlations not only on ground states but also on arbitrary excited states can be computed efficiently. This thus allows to deal with finite temperature correlations. -############################################################## - -IMPROVEMENTS YET TO BE IMPLEMENTED: - -- Some form of binning to make very large runs possible. -Possibility: post-processing of the .raw file, joining together contributions at given iK and (about) same energy into effective ones and saving those into a pre-defined matrix format. -Difficulties: the size of the .raw file is of the order of the .thr one, which must be kept. - -############################################################## - -VERSIONS: NOTES - Important conventions: - Versions are numbered with two integers, [VERSION].[SUBVERSION] -- Changes in subversion number indicate `internal' revisions which do not change any of the external conventions, data file structures, etc. Data produced with version V.S can be `polished' with version V.S' with S' != S. +- Changes in subversion number indicate internal revisions which do not change any of the external conventions, data file structures, etc. Data produced with version V.S can be polished with version V.S' with S' != S. - Changes in version number indicate a larger scale revision affecting the conventions. Data produced with earlier versions are then deprecated. ---------------------------------------------------------------- - -ABACUS++T_1.0 [published 23 June 2011] - -First version of new setup and logic for ABACUS. Implementation of all basic new ideas for state labelling and descending. ------------------------------------------------------------------- +#### ABACUS++T_9 +_[in development][abandoned]_ -ABACUS++T_2.0 [in progress][*ABANDONED*] - -Attempt at implementing fixed-momentum-based scanning. +Scanning over ensembles. ------------------------------------------------------------------- +#### ABACUS++T_8.0 -ABACUS++T_3.0 [published 1 Nov 2011] +Changed labels slightly, to avoid identifying empty string with integer 0. +All labels with at least one excitation are now strictly nonempty. + + +#### ABACUS++T_7.0 +_[published 1 Dec 2011]_ + +Optimization of runtime memory use and of output files sizes. + +Labels: introduction of compressed labels (used in raw and thr files). +Removed conv boolean from raw file (all states in this file have converged). +Threads: removed omega and iK (changed `General_Scan` accordingly), use compressed labels. + +Scanning function: it's now possible to give a default file name when invoking +`General_Scan` (and thus `Scan_Bose`, etc). This is done to avoid stupidly long file +names when calculating correlators over `AveragingState`s which are far from +the ground state. + + +#### ABACUS++T_6.0 +_[published 21 Nov 2011]_ + +The scanning is now also done over the hole positions. +Upon the creation of a particle-hole, only hole positions at the edges of +blocks of contiguous `Ix2` in `OriginState` are used. +The holes are then scanned moving towards the center of the blocks. + +There are now thus 3 types of scanning: +0: over holes +1: over particles +2: adding a p-h pair + + +#### ABACUS++T_5.0 +_[published 4 Nov 2011]_ + +Introduced two types of threads, separating scanning over fixed particle-hole numbers +and adding a p-h. + + +#### ABACUS++T_4.0 [published 1 Nov 2011] + +Modification to the labeling logic: +states are now always labeled using the AveragingState's quantum numbers, +even if their bases are different. The form of the `State_Label` is thus changed. + + +#### ABACUS++T_3.0 +_[published 1 Nov 2011]_ New, simpler logic for descendents: the hole positions are scanned immediately upon creation of a new particle-hole pair; the hole positions are then kept fixed @@ -176,60 +232,17 @@ Problems: the extra particle in the AveragingState is not scanned to. Version 4 will address this problem. ------------------------------------------------------------------- -ABACUS++T_4.0 [published 1 Nov 2011] +#### ABACUS++T_2.0 +_[in progress][abandoned]_ -Modification to the labeling logic: -states are now always labeled using the AveragingState's quantum numbers, -even if their bases are different. The form of the State_Label is thus changed. +Attempt at implementing fixed-momentum-based scanning. ------------------------------------------------------------------- +#### ABACUS++T_1.0 +_[published 23 June 2011]_ -ABACUS++T_5.0 [published 4 Nov 2011] - -Introduced two types of threads, separating scanning over fixed particle-hole numbers -and adding a p-h. - ------------------------------------------------------------------- - -ABACUS++T_6.0 [published 21 Nov 2011] - -The scanning is now also done over the hole positions. -Upon the creation of a particle-hole, only hole positions at the edges of -blocks of contiguous Ix2 in OriginState are used. -The holes are then scanned moving towards the center of the blocks. - -There are now thus 3 types of scanning: -0: over holes -1: over particles -2: adding a p-h pair - ------------------------------------------------------------------- - -ABACUS++T_7.0 [published 1 Dec 2011]] - -Optimization of runtime memory use and of output files sizes. - -Labels: introduction of compressed labels (used in raw and thr files). -Removed conv boolean from raw file (all states in this file have converged). -Threads: removed omega and iK (changed General_Scan accordingly), use compressed labels. - -Scanning function: it's now possible to give a default file name when invoking -General_Scan (and thus Scan_Bose, etc). This is done to avoid stupidly long file -names when calculating correlators over AveragingStates which are far from -the ground state. - ------------------------------------------------------------------- - -ABACUS++T_8.0 - -Changed labels slightly, to avoid identifying empty string with integer 0. -All labels with at least one excitation are now strictly nonempty. +First version of new setup and logic for ABACUS. Implementation of all basic new ideas for state labelling and descending. ------------------------------------------------------------------- - -ABACUS++T_9 [in development][abandoned] - -Scanning over ensembles. +### Earlier versions +Earlier version of ABACUS are not documented here. The whole history is available in J.-S. Caux's private backfiles. diff --git a/src/HEIS/Heis_Matrix_Element_Contrib.cc b/src/HEIS/Heis_Matrix_Element_Contrib.cc index bf7821e..0e2799d 100644 --- a/src/HEIS/Heis_Matrix_Element_Contrib.cc +++ b/src/HEIS/Heis_Matrix_Element_Contrib.cc @@ -20,11 +20,11 @@ using namespace JSC; namespace JSC { - //DP Compute_Matrix_Element_Contrib (char whichDSF, bool fixed_iK, XXZ_Bethe_State& LeftState, + //DP Compute_Matrix_Element_Contrib (char whichDSF, bool fixed_iK, XXZ_Bethe_State& LeftState, //XXZ_Bethe_State& RightState, DP Chem_Pot, fstream& DAT_outfile) - //DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, XXZ_Bethe_State& LeftState, + //DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, XXZ_Bethe_State& LeftState, // XXZ_Bethe_State& RightState, DP Chem_Pot, fstream& DAT_outfile) - DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, XXZ_Bethe_State& LeftState, + DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, XXZ_Bethe_State& LeftState, XXZ_Bethe_State& RightState, DP Chem_Pot, stringstream& DAT_outfile) { // This function prints the matrix information to the fstream, @@ -61,20 +61,21 @@ namespace JSC { while(iKout < 0) iKout += RightState.chain.Nsites; while(iKout >= RightState.chain.Nsites) iKout -= RightState.chain.Nsites; + DAT_outfile << setprecision(16); // Print information to fstream: if (iKout >= iKmin && iKout <= iKmax) { if (whichDSF == 'Z') { - DAT_outfile << endl << setprecision(16) << LeftState.E - RightState.E - (LeftState.base.Mdown - RightState.base.Mdown) * Chem_Pot << "\t" + DAT_outfile << endl << setprecision(16) << LeftState.E - RightState.E - (LeftState.base.Mdown - RightState.base.Mdown) * Chem_Pot << "\t" << iKout << "\t" - //<< LeftState.conv << "\t" + //<< LeftState.conv << "\t" << setprecision(3) << LeftState.dev << "\t" << LeftState.label; } else { - DAT_outfile << endl << setprecision(16) << LeftState.E - RightState.E - (LeftState.base.Mdown - RightState.base.Mdown) * Chem_Pot << "\t" + DAT_outfile << endl << setprecision(16) << LeftState.E - RightState.E - (LeftState.base.Mdown - RightState.base.Mdown) * Chem_Pot << "\t" << iKout << "\t" << ME << "\t" - //<< LeftState.conv << "\t" + //<< LeftState.conv << "\t" << setprecision(3) << LeftState.dev << "\t" << LeftState.label; } @@ -92,11 +93,11 @@ namespace JSC { return(data_value); } - //DP Compute_Matrix_Element_Contrib (char whichDSF, bool fixed_iK, XXX_Bethe_State& LeftState, + //DP Compute_Matrix_Element_Contrib (char whichDSF, bool fixed_iK, XXX_Bethe_State& LeftState, //XXX_Bethe_State& RightState, DP Chem_Pot, fstream& DAT_outfile) - //DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, XXX_Bethe_State& LeftState, + //DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, XXX_Bethe_State& LeftState, //XXX_Bethe_State& RightState, DP Chem_Pot, fstream& DAT_outfile) - DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, XXX_Bethe_State& LeftState, + DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, XXX_Bethe_State& LeftState, XXX_Bethe_State& RightState, DP Chem_Pot, stringstream& DAT_outfile) { // This function prints the matrix element information to the fstream, @@ -137,7 +138,7 @@ namespace JSC { if (LeftState.base.Mdown == RightState.base.Mdown - 1) { // two infinite rapidities, use rescaled S^- matrix element instead of S^+ nrinfrap = 2; // Correction factor for MEsq: Smffsq to Spffsq = 2/((N - 2M + 2) (N - 2M + 1)) - ME = sqrt(2.0/((RightState.chain.Nsites - 2*RightState.base.Mdown + 2.0) * (RightState.chain.Nsites - 2*RightState.base.Mdown + 1.0))) + ME = sqrt(2.0/((RightState.chain.Nsites - 2*RightState.base.Mdown + 2.0) * (RightState.chain.Nsites - 2*RightState.base.Mdown + 1.0))) * exp(real(ln_Smin_ME (RightState, LeftState))); } else if (LeftState.base.Mdown == RightState.base.Mdown) { // one infinite rapidity, use rescaled S^z matrix element instead of S^+ @@ -170,29 +171,30 @@ namespace JSC { while(iKout < 0) iKout += RightState.chain.Nsites; while(iKout >= RightState.chain.Nsites) iKout -= RightState.chain.Nsites; + DAT_outfile << setprecision(16); // Print information to fstream: if (iKout >= iKmin && iKout <= iKmax) { if (whichDSF == 'Z') { - DAT_outfile << endl << setprecision(16) << LeftState.E - RightState.E - (LeftState.base.Mdown + nrinfrap - RightState.base.Mdown) * Chem_Pot << "\t" + DAT_outfile << endl << setprecision(16) << LeftState.E - RightState.E - (LeftState.base.Mdown + nrinfrap - RightState.base.Mdown) * Chem_Pot << "\t" << iKout << "\t" - //<< LeftState.conv << "\t" + //<< LeftState.conv << "\t" << setprecision(3) << LeftState.dev << "\t" << LeftState.label; } else if (whichDSF == 'q') { - DAT_outfile << endl << setprecision(16) << LeftState.E - RightState.E - (LeftState.base.Mdown + nrinfrap - RightState.base.Mdown) * Chem_Pot << "\t" + DAT_outfile << endl << setprecision(16) << LeftState.E - RightState.E - (LeftState.base.Mdown + nrinfrap - RightState.base.Mdown) * Chem_Pot << "\t" << iKout << "\t" << real(ME_CX) << "\t" << imag(ME_CX) - twoPI * int(imag(ME_CX)/twoPI + 1.0e-10) << "\t" - //<< LeftState.conv << "\t" + //<< LeftState.conv << "\t" << setprecision(3) << LeftState.dev << "\t" << LeftState.label; } else { - DAT_outfile << endl << setprecision(16) << LeftState.E - RightState.E - (LeftState.base.Mdown + nrinfrap - RightState.base.Mdown) * Chem_Pot << "\t" + DAT_outfile << endl << setprecision(16) << LeftState.E - RightState.E - (LeftState.base.Mdown + nrinfrap - RightState.base.Mdown) * Chem_Pot << "\t" << iKout << "\t" << ME << "\t" - //<< LeftState.conv << "\t" + //<< LeftState.conv << "\t" << setprecision(3) << LeftState.dev << "\t" << LeftState.label; } @@ -203,7 +205,7 @@ namespace JSC { //DP data_value = (iKout == 0 ? 1.0 : 2.0) * MEsq; if (whichDSF == 'Z') // use 1/(1 + omega) data_value = 1.0/(1.0 + LeftState.E - RightState.E - (LeftState.base.Mdown - RightState.base.Mdown) * Chem_Pot); - else if (whichDSF == 'q') + else if (whichDSF == 'q') data_value = exp(2.0 * real(ME_CX)); else if (fixed_iK) // data value is MEsq * omega: data_value = ME * ME * (LeftState.E - RightState.E - (LeftState.base.Mdown - RightState.base.Mdown) * Chem_Pot); @@ -211,11 +213,11 @@ namespace JSC { return(data_value); } - //DP Compute_Matrix_Element_Contrib (char whichDSF, bool fixed_iK, XXZ_gpd_Bethe_State& LeftState, + //DP Compute_Matrix_Element_Contrib (char whichDSF, bool fixed_iK, XXZ_gpd_Bethe_State& LeftState, //XXZ_gpd_Bethe_State& RightState, DP Chem_Pot, fstream& DAT_outfile) - //DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, XXZ_gpd_Bethe_State& LeftState, + //DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, XXZ_gpd_Bethe_State& LeftState, // XXZ_gpd_Bethe_State& RightState, DP Chem_Pot, fstream& DAT_outfile) - DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, XXZ_gpd_Bethe_State& LeftState, + DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, XXZ_gpd_Bethe_State& LeftState, XXZ_gpd_Bethe_State& RightState, DP Chem_Pot, stringstream& DAT_outfile) { // This function prints the matrix element information to the fstream, @@ -233,7 +235,7 @@ namespace JSC { */ bool fixed_iK = (iKmin == iKmax); - + // If any of the rapidities outside of fundamental interval -pi/2 < lambda <= pi/2, set matrix element to zero. bool rap_in_fundamental = true; @@ -254,9 +256,9 @@ namespace JSC { for (int k = 0; k < LeftState.chain.Nstrings; ++k) sum1 += LeftState.base.Nrap[k] * (2 * JSC::min(LeftState.chain.Str_L[j], LeftState.chain.Str_L[k]) - ((j == k) ? 1 : 0)); // This almost does it: only missing are the states with one on -PI/2 and one on PI/2 - if (LeftState.base.Nrap[j] >= 1 - && (LeftState.Ix2[j][0] <= -(LeftState.chain.Nsites - sum1) - || (LeftState.Ix2[j][LeftState.base.Nrap[j] - 1] - LeftState.Ix2[j][0]) + if (LeftState.base.Nrap[j] >= 1 + && (LeftState.Ix2[j][0] <= -(LeftState.chain.Nsites - sum1) + || (LeftState.Ix2[j][LeftState.base.Nrap[j] - 1] - LeftState.Ix2[j][0]) > 2*(LeftState.chain.Nsites - sum1))) rap_in_fundamental = false; */ @@ -267,7 +269,7 @@ namespace JSC { // rap_in_fundamental = false; //} /* - if (LeftState.base.Nrap[j] > 0 && + if (LeftState.base.Nrap[j] > 0 && ((LeftState.lambda[j][LeftState.base.Nrap[j] - 1] - LeftState.lambda[j][0] >= PI) || LeftState.lambda[j][0] < -0.5*PI + 1.0e-10 || LeftState.lambda[j][LeftState.base.Nrap[j] - 1] > 0.5*PI @@ -279,7 +281,7 @@ namespace JSC { rap_in_fundamental = false; */ /* - if (LeftState.base.Nrap[j] > 0 && + if (LeftState.base.Nrap[j] > 0 && ((LeftState.lambda[j][LeftState.base.Nrap[j] - 1] - LeftState.lambda[j][0] >= PI) //|| (LeftState.base.Nrap[j] == 1 && fabs(LeftState.lambda[j][0]) > 0.5*PI) )) @@ -287,10 +289,10 @@ namespace JSC { */ // Logic: allow rapidities -PI/2 <= lambda <= PI/2 (including boundaries) - if (LeftState.base.Nrap[j] > 0 && + if (LeftState.base.Nrap[j] > 0 && (LeftState.lambda[j][0] < -PI/2 || LeftState.lambda[j][LeftState.base.Nrap[j] - 1] > PI/2)) rap_in_fundamental = false; - if (RightState.base.Nrap[j] > 0 && + if (RightState.base.Nrap[j] > 0 && (RightState.lambda[j][0] < -PI/2 || RightState.lambda[j][RightState.base.Nrap[j] - 1] > PI/2)) rap_in_fundamental = false; @@ -334,20 +336,21 @@ namespace JSC { while(iKout < 0) iKout += RightState.chain.Nsites; while(iKout >= RightState.chain.Nsites) iKout -= RightState.chain.Nsites; + DAT_outfile << setprecision(16); // Print information to fstream: if (iKout >= iKmin && iKout <= iKmax) { if (whichDSF == 'Z') { - DAT_outfile << endl << setprecision(16) << LeftState.E - RightState.E - (LeftState.base.Mdown - RightState.base.Mdown) * Chem_Pot << "\t" + DAT_outfile << endl << setprecision(16) << LeftState.E - RightState.E - (LeftState.base.Mdown - RightState.base.Mdown) * Chem_Pot << "\t" << iKout << "\t" - //<< LeftState.conv << "\t" + //<< LeftState.conv << "\t" << setprecision(3) << LeftState.dev << "\t" << LeftState.label; } else { - DAT_outfile << endl << setprecision(16) << LeftState.E - RightState.E - (LeftState.base.Mdown - RightState.base.Mdown) * Chem_Pot << "\t" + DAT_outfile << endl << setprecision(16) << LeftState.E - RightState.E - (LeftState.base.Mdown - RightState.base.Mdown) * Chem_Pot << "\t" << iKout << "\t" << ME << "\t" - //<< LeftState.conv << "\t" + //<< LeftState.conv << "\t" << setprecision(3) << LeftState.dev << "\t" << LeftState.label; diff --git a/src/LIEBLIN/LiebLin_Matrix_Element_Contrib.cc b/src/LIEBLIN/LiebLin_Matrix_Element_Contrib.cc index 7f7ee52..92278ec 100644 --- a/src/LIEBLIN/LiebLin_Matrix_Element_Contrib.cc +++ b/src/LIEBLIN/LiebLin_Matrix_Element_Contrib.cc @@ -20,11 +20,11 @@ using namespace JSC; namespace JSC { - //DP Compute_Matrix_Element_Contrib (char whichDSF, bool fixed_iK, LiebLin_Bethe_State& LeftState, + //DP Compute_Matrix_Element_Contrib (char whichDSF, bool fixed_iK, LiebLin_Bethe_State& LeftState, // LiebLin_Bethe_State& RefState, DP Chem_Pot, fstream& DAT_outfile) - //DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, LiebLin_Bethe_State& LeftState, + //DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, LiebLin_Bethe_State& LeftState, // LiebLin_Bethe_State& RefState, DP Chem_Pot, fstream& DAT_outfile) - DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, LiebLin_Bethe_State& LeftState, + DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, LiebLin_Bethe_State& LeftState, LiebLin_Bethe_State& RefState, DP Chem_Pot, stringstream& DAT_outfile) { // This function prints the matrix element information to the fstream, @@ -45,7 +45,7 @@ namespace JSC { ME = real(exp(ln_Psi_ME (LeftState, RefState))); else if (whichDSF == 'g') ME = real(exp(ln_Psi_ME (RefState, LeftState))); - else if (whichDSF == 'q') // geometrical quench + else if (whichDSF == 'q') // geometrical quench //ME_CX = (LiebLin_Twisted_ln_Overlap(1.0, RefState.lambdaoc, RefState.lnnorm, LeftState)); //ME_CX = (LiebLin_ln_Overlap(RefState.lambdaoc, RefState.lnnorm, LeftState)); ME_CX = LiebLin_ln_Overlap(LeftState.lambdaoc, LeftState.lnnorm, RefState); @@ -56,66 +56,67 @@ namespace JSC { // The product is ME_CX * ME = e^{-\delta S_e} \langle \rho_{sp} | g2 (x=0) | \rho_{sp} + e \rangle // and is the first half of the t-dep expectation value formula coming from the Quench Action formalism } - else if (whichDSF == 'C') { // BEC to finite c quench: overlap + else if (whichDSF == 'C') { // BEC to finite c quench: overlap ME_CX = exp(2.0* ln_Overlap_with_BEC (LeftState)); // overlap sq part ME = 0.0; } else JSCerror("Wrong whichDSF in Compute_Matrix_Element_Contrib."); - if (is_nan(ME)) ME = (whichDSF == 'Z') ? 1.0e+200 : 0.0; + if (is_nan(ME)) ME = (whichDSF == 'Z') ? 1.0e+200 : 0.0; if (is_nan(norm(ME_CX))) ME_CX = -100.0; - + + DAT_outfile << setprecision(16); // Print information to fstream: if (LeftState.iK - RefState.iK >= iKmin && LeftState.iK - RefState.iK <= iKmax) { if (whichDSF == 'Z') { - DAT_outfile << endl << LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot << "\t" - << LeftState.iK - RefState.iK - //<< "\t" << LeftState.conv + DAT_outfile << endl << LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot << "\t" + << LeftState.iK - RefState.iK + //<< "\t" << LeftState.conv << 0 << "\t" // This is the deviation, here always 0 << "\t" << LeftState.label; } else if (whichDSF == 'q') { - DAT_outfile << endl << LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot << "\t" + DAT_outfile << endl << LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot << "\t" << LeftState.iK - RefState.iK << "\t" << real(ME_CX) << "\t" << imag(ME_CX) - twoPI * int(imag(ME_CX)/twoPI + 1.0e-10) << "\t" - //<< LeftState.conv << "\t" + //<< LeftState.conv << "\t" << 0 << "\t" // This is the deviation, here always 0 << LeftState.label; } else if (whichDSF == '1') { - DAT_outfile << endl << LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot << "\t" + DAT_outfile << endl << LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot << "\t" << LeftState.iK - RefState.iK << "\t" << ME << "\t" - //<< LeftState.conv << "\t" + //<< LeftState.conv << "\t" << 0 << "\t" // This is the deviation, here always 0 << LeftState.label; DAT_outfile << "\t" << Momentum_Right_Excitations(LeftState) << "\t" << Momentum_Left_Excitations(LeftState); //cout << Momentum_Right_Excitations(LeftState) << "\t" << Momentum_Left_Excitations(LeftState); } else if (whichDSF == 'B') { // BEC to finite c > 0 quench; g2 (x=0) - if (fabs(real(ME_CX) * ME) > 1.0e-100) + if (fabs(real(ME_CX) * ME) > 1.0e-100) DAT_outfile << endl << LeftState.E - RefState.E << "\t" << LeftState.iK - RefState.iK << "\t" << real(ME_CX) << "\t" // the overlap is always real - << ME << "\t" + << ME << "\t" //<< 0 << "\t" // This is the deviation, here always 0 // omit this here << LeftState.label; } else if (whichDSF == 'C') { // BEC to finite c, overlap sq - if (fabs(real(ME_CX)) > 1.0e-100) + if (fabs(real(ME_CX)) > 1.0e-100) DAT_outfile << endl << LeftState.E - RefState.E << "\t" << LeftState.iK - RefState.iK << "\t" << real(ME_CX) << "\t" // the overlap is always real - << ME << "\t" + << ME << "\t" //<< 0 << "\t" // This is the deviation, here always 0 // omit this here << LeftState.label; } else { - DAT_outfile << endl << LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot << "\t" + DAT_outfile << endl << LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot << "\t" << LeftState.iK - RefState.iK << "\t" << ME << "\t" - //<< LeftState.conv << "\t" + //<< LeftState.conv << "\t" << 0 << "\t" // This is the deviation, here always 0 << LeftState.label; } @@ -126,36 +127,36 @@ namespace JSC { if (whichDSF == 'Z') // use 1/(1 + omega) data_value = 1.0/(1.0 + LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot); else if (whichDSF == 'd' || whichDSF == '1') { - if (fixed_iK) + if (fixed_iK) /* - // use omega * MEsq/iK^2 - //data_value = (LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot) + // use omega * MEsq/iK^2 + //data_value = (LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot) // MEsq/JSC::max(1, (LeftState.iK - RefState.iK) * (LeftState.iK - RefState.iK)) - //: 0.0; + //: 0.0; */ // Careful: use fabs(E) since this must also work with Tgt0 or arbitrary RefState DEPRECATED ++G_1, USE abs_data_value data_value = (LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot) * ME * ME; else if (!fixed_iK) // use ME if momentum in fundamental window -iK_UL <= iK <= iK_UL - //data_value = abs(LeftState.iK - RefState.iK) <= RefState.Tableau[0].Ncols ? // this last nr is iK_UL + //data_value = abs(LeftState.iK - RefState.iK) <= RefState.Tableau[0].Ncols ? // this last nr is iK_UL //MEsq : 0.0; //data_value = (LeftState.iK - RefState.iK == 0 ? 1.0 : 2.0) * MEsq; //data_value = ME * ME; - // use omega * MEsq/iK^2 + // use omega * MEsq/iK^2 // Careful: use fabs(E) since this must also work with Tgt0 or arbitrary RefState DEPRECATED ++G_1, USE abs_data_value data_value = (LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot) * ME * ME/JSC::max(1, (LeftState.iK - RefState.iK) * (LeftState.iK - RefState.iK)); } else if (whichDSF == 'g' || whichDSF == 'o') { - if (fixed_iK) + if (fixed_iK) /* // use omega * MEsq/((k^2 - mu + 4 c N/L)/L) - data_value = (LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot) + data_value = (LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot) * MEsq/((((LeftState.K - RefState.K) * (LeftState.K - RefState.K) - Chem_Pot) + 4.0 * RefState.c_int * RefState.N/RefState.L)/RefState.L) : 0.0; */ // Careful: use fabs(E) since this must also work with Tgt0 or arbitrary RefState data_value = (LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot) * ME * ME; else if (!fixed_iK) { // simply use MEsq if momentum in fundamental window -iK_UL <= iK <= iK_UL - //data_value = abs(LeftState.iK - RefState.iK) <= RefState.Tableau[0].Ncols ? // this last nr is iK_UL + //data_value = abs(LeftState.iK - RefState.iK) <= RefState.Tableau[0].Ncols ? // this last nr is iK_UL //MEsq : 0.0; //data_value = (LeftState.iK - RefState.iK == 0 ? 1.0 : 2.0) * MEsq; data_value = ME * ME; @@ -165,9 +166,9 @@ namespace JSC { } } //else if (whichDSF == 'o') - //data_value = abs(LeftState.iK - RefState.iK) <= RefState.Tableau[0].Ncols ? // this last nr is iK_UL + //data_value = abs(LeftState.iK - RefState.iK) <= RefState.Tableau[0].Ncols ? // this last nr is iK_UL //MEsq : 0.0; - else if (whichDSF == 'q') + else if (whichDSF == 'q') data_value = exp(2.0 * real(ME_CX)); else if (whichDSF == 'B') data_value = abs(ME_CX * ME)/(1.0 + sqrt(fabs(LeftState.E - RefState.E))); diff --git a/src/ODSLF/ODSLF_Matrix_Element_Contrib.cc b/src/ODSLF/ODSLF_Matrix_Element_Contrib.cc index 29afa3c..1f18729 100644 --- a/src/ODSLF/ODSLF_Matrix_Element_Contrib.cc +++ b/src/ODSLF/ODSLF_Matrix_Element_Contrib.cc @@ -20,9 +20,9 @@ using namespace JSC; namespace JSC { - //DP Compute_Matrix_Element_Contrib (char whichDSF, bool fixed_iK, ODSLF_XXZ_Bethe_State& LeftState, + //DP Compute_Matrix_Element_Contrib (char whichDSF, bool fixed_iK, ODSLF_XXZ_Bethe_State& LeftState, //ODSLF_XXZ_Bethe_State& RefState, DP Chem_Pot, fstream& DAT_outfile) - DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, ODSLF_XXZ_Bethe_State& LeftState, + DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, ODSLF_XXZ_Bethe_State& LeftState, ODSLF_XXZ_Bethe_State& RefState, DP Chem_Pot, fstream& DAT_outfile) { // This function prints the matrix element information to the fstream, @@ -55,19 +55,20 @@ namespace JSC { while(iKout < 0) iKout += RefState.chain.Nsites; while(iKout >= RefState.chain.Nsites) iKout -= RefState.chain.Nsites; + DAT_outfile << setprecision(16); // Print information to fstream: if (iKout >= iKmin && iKout <= iKmax) { if (whichDSF == 'Z') { - DAT_outfile << endl << LeftState.E - RefState.E - (LeftState.base.Mdown - RefState.base.Mdown) * Chem_Pot << "\t" + DAT_outfile << endl << LeftState.E - RefState.E - (LeftState.base.Mdown - RefState.base.Mdown) * Chem_Pot << "\t" << iKout << "\t" - //<< LeftState.conv << "\t" + //<< LeftState.conv << "\t" << LeftState.base_id << "\t" << LeftState.type_id << "\t" << LeftState.id; } else { - DAT_outfile << endl << LeftState.E - RefState.E - (LeftState.base.Mdown - RefState.base.Mdown) * Chem_Pot << "\t" + DAT_outfile << endl << LeftState.E - RefState.E - (LeftState.base.Mdown - RefState.base.Mdown) * Chem_Pot << "\t" << iKout << "\t" << ME << "\t" - //<< LeftState.conv << "\t" + //<< LeftState.conv << "\t" << LeftState.base_id << "\t" << LeftState.type_id << "\t" << LeftState.id; } } // if iKmin <= iKout <= iKmax diff --git a/src/SCAN/General_Scan.cc b/src/SCAN/General_Scan.cc index f272805..822e740 100644 --- a/src/SCAN/General_Scan.cc +++ b/src/SCAN/General_Scan.cc @@ -8,10 +8,10 @@ Copyright (c). File: src/SCAN/General_Scan.cc -Purpose: universal implementation of state scanning: +Purpose: universal implementation of state scanning: functions to descend down hierarchy of intermediate states. -NOTE: since templated functions have to be in the same file, +NOTE: since templated functions have to be in the same file, we put all scanning functions here. The externally-accessible functions are defined at the end of this file. @@ -61,7 +61,7 @@ namespace JSC { // - there exists an OriginIx2 between the active Ix2 and the next Ix2 (to right or left depending on type of descendent) Vect > ScanIx2 = Return_Ix2_from_Label (ScanIx2_label, OriginIx2); - + // Determine the level and index of the bottom-most left-most right-moving quantum number sits: int exclevel = -1; int excindex = 0; @@ -69,9 +69,9 @@ namespace JSC { do { exclevel++; if (exclevel == ScanIx2.size()) { // there isn't a single right-moving quantum number in ScanIx2 - break; + break; } - for (int alpha = 0; alpha < ScanIx2[exclevel].size(); ++alpha) + for (int alpha = 0; alpha < ScanIx2[exclevel].size(); ++alpha) if (ScanIx2[exclevel][alpha] > BaseScanIx2[exclevel][alpha]) { excindex = alpha; excfound = true; @@ -81,10 +81,10 @@ namespace JSC { // If we haven't found an excitation, then exclevel == ScanIx2.size() and excindex = 0; if (excfound && !BaseScanIx2[exclevel].includes(ScanIx2[exclevel][excindex])) { // there exists an already dispersing excitation which isn't in Origin - // Is there a possible recombination? + // Is there a possible recombination? if (excindex < ScanIx2[exclevel].size() - 1) { // a particle to the right of excitation has already move right, so there is a hole // check that there exists an occupied Ix2 in Origin sitting between the excitation and the next Ix2 to its right in ScanIx2 - for (int alpha = BaseScanIx2[exclevel].size() - 1; alpha >= 0; --alpha) + for (int alpha = BaseScanIx2[exclevel].size() - 1; alpha >= 0; --alpha) if (BaseScanIx2[exclevel][alpha] > ScanIx2[exclevel][excindex] && BaseScanIx2[exclevel][alpha] < ScanIx2[exclevel][excindex + 1]) { return(true); } @@ -125,14 +125,14 @@ namespace JSC { int exclevel = -1; int excindex = 0; bool excfound = false; - + //cout << "Looking for exclevel and excindex for " << endl << "\tBaseIx2 = " << BaseScanIx2 << endl << "\tScanIx2 = " << ScanIx2 << endl; do { exclevel++; if (exclevel == ScanIx2.size()) { // there isn't a single left-moving quantum number in ScanIx2 - break; + break; } - for (int alpha = ScanIx2[exclevel].size() - 1; alpha >= 0; --alpha) { + for (int alpha = ScanIx2[exclevel].size() - 1; alpha >= 0; --alpha) { //cout << exclevel << "\t" << alpha << "\t" << ScanIx2[exclevel][alpha] << "\t" << BaseScanIx2[exclevel][alpha] << "\t" << (ScanIx2[exclevel][alpha] < BaseScanIx2[exclevel][alpha]) << endl; if (ScanIx2[exclevel][alpha] < BaseScanIx2[exclevel][alpha]) { excindex = alpha; @@ -142,13 +142,13 @@ namespace JSC { } } while (!excfound); // If we haven't found an excitation, then exclevel == ScanIx2.size() and excindex = 0; - if (!excfound) excindex = ScanIx2[exclevel].size() - 1; - + if (!excfound) excindex = ScanIx2[exclevel].size() - 1; + if (excfound && !BaseScanIx2[exclevel].includes(ScanIx2[exclevel][excindex])) { // there exists an already dispersing excitation which isn't in Origin - // Is there a possible recombination? + // Is there a possible recombination? if (excindex > 0) { // a particle to the left of excitation has already moved left, so there is a hole // check that there exists an occupied Ix2 in Origin sitting between the excitation and the next Ix2 to its left in ScanIx2 - for (int alpha = 0; alpha < BaseScanIx2[exclevel].size(); ++alpha) + for (int alpha = 0; alpha < BaseScanIx2[exclevel].size(); ++alpha) if (BaseScanIx2[exclevel][alpha] > ScanIx2[exclevel][excindex - 1] && BaseScanIx2[exclevel][alpha] < ScanIx2[exclevel][excindex]) { return(true); } @@ -180,12 +180,12 @@ namespace JSC { template void Descend_and_Compute_for_Fixed_Base (char whichDSF, Tstate& AveragingState, Tstate& BaseScanState, Tstate& ScanState, int type_required, int iKmin, int iKmax, int iKmod, - //Scan_Thread_List& paused_thread_list, - //Scan_Thread_Set& paused_thread_set, - Scan_Thread_Data& paused_thread_data, - //thresholdremoved DP& running_scan_threshold, //DP ref_abs_data_value, - DP& ph_cost, int Max_Secs, DP sumrule_factor, DP Chem_Pot, Scan_Info& scan_info, - fstream& RAW_outfile, fstream& INADM_outfile, int& ninadm, + //Scan_Thread_List& paused_thread_list, + //Scan_Thread_Set& paused_thread_set, + Scan_Thread_Data& paused_thread_data, + //thresholdremoved DP& running_scan_threshold, //DP ref_abs_data_value, + DP& ph_cost, int Max_Secs, DP sumrule_factor, DP Chem_Pot, Scan_Info& scan_info, + fstream& RAW_outfile, fstream& INADM_outfile, int& ninadm, fstream& CONV0_outfile, int& nconv0, fstream& STAT_outfile) { @@ -207,11 +207,11 @@ namespace JSC { if (whichDSF == 'B') { // symmetric state scanning if (type_required >= 9 && type_required <= 11) desc_label = Descendent_States_with_iK_Stepped_Down_rightIx2only (ScanState.label, BaseScanState, disperse_only_current_exc_down, preserve_nexc_down); - else if (type_required >= 12 && type_required <= 14) + else if (type_required >= 12 && type_required <= 14) desc_label = Descendent_States_with_iK_Stepped_Up_rightIx2only (ScanState.label, BaseScanState, disperse_only_current_exc_up, preserve_nexc_up); } else { - if (type_required >= 0 && type_required <= 8) { + if (type_required >= 0 && type_required <= 8) { desc_label = Descendent_States_with_iK_Preserved(ScanState.label, BaseScanState, disperse_only_current_exc_up, preserve_nexc_up, disperse_only_current_exc_down, preserve_nexc_down); } else if (type_required >= 9 && type_required <= 11) @@ -226,16 +226,16 @@ namespace JSC { //cout << "OK for descend on " << ScanState.label << " with type_required = " << type_required << endl; //cout << desc_label << endl; - + //Vect desc_label = desc_data.label; //Vect desc_type = desc_data.type; - + //bool disp_OK = false; - + //if (ScanState.label == "7|2:1_0|1_|0@2") { //if (ScanState.label == "64_1_63@319") { //if (ScanState.label == "32_1_-29@-33") { - + if (ScanState.label == LABEL_TO_CHECK) { cout << "Called Descend on state " << ScanState << endl; cout << "For type_required == " << type_required << ", found " << desc_label.size() << " descendents, "; @@ -252,19 +252,19 @@ namespace JSC { cin >> LABEL_TO_CHECK; } } - - + + string label_here = ScanState.label; //int ScanState_iK = ScanState.iK; - + for (int idesc = 0; idesc < desc_label.size(); ++idesc) { - + //cout << "\tDealing with descendent " << idesc << " out of " << desc_label.size() << " with label " << desc_label[idesc] << endl; //cout << "\tfrom state with label " << label_here << " and of type_required " << type_required << endl; clock_t start_time_here = clock(); - //if (desc_label[idesc] == "64_2_0yvv7") { + //if (desc_label[idesc] == "64_2_0yvv7") { if (false) { cout << "Found " << desc_label[idesc] << " as descendent of type " << type_required << " of " << label_here << endl; ScanState.Set_to_Label (label_here, BaseScanState.Ix2); @@ -284,22 +284,22 @@ namespace JSC { ScanState.Set_to_Label (desc_label[idesc], BaseScanState.Ix2); bool admissible = ScanState.Check_Admissibility(whichDSF); - + DP data_value = 0.0; - + scan_info.Ndata++; - + ScanState.conv = false; ScanState.Compute_Momentum(); // since momentum is used as forced descent criterion - + if (admissible) { - + ScanState.Compute_All (idesc == 0); //ScanState.Compute_All (true); //scan_info.Ndata++; - + if (ScanState.conv) { scan_info.Ndata_conv++; @@ -307,7 +307,7 @@ namespace JSC { int iKexc = ScanState.iK - AveragingState.iK; while (iKexc > iKmax && iKexc - iKmod >= iKmin) iKexc -= iKmod; while (iKexc < iKmin && iKexc + iKmod <= iKmax) iKexc += iKmod; - + data_value = Compute_Matrix_Element_Contrib (whichDSF, iKmin, iKmax, ScanState, AveragingState, Chem_Pot, RAW_outfile); if (iKexc >= iKmin && iKexc <= iKmax) scan_info.sumrule_obtained += data_value*sumrule_factor; //cout << "data_value found = " << data_value * sumrule_factor << endl; @@ -317,14 +317,14 @@ namespace JSC { } // if (ScanState.conv) else { - if (nconv0++ < 1000) - CONV0_outfile << setw(25) << ScanState.label << setw(25) << ScanState.diffsq << setw(5) << ScanState.Check_Rapidities() + if (nconv0++ < 1000) + CONV0_outfile << setw(25) << ScanState.label << setw(25) << ScanState.diffsq << setw(5) << ScanState.Check_Rapidities() << setw(25) << ScanState.String_delta() << endl; scan_info.Ndata_conv0++; //cout << "State did not converge." << endl; } } // if (admissible) - + else { if (ninadm++ < 1000000) INADM_outfile << ScanState.label << endl; scan_info.Ninadm++; @@ -334,7 +334,7 @@ namespace JSC { } clock_t stop_time_here = clock(); - + scan_info.CPU_ticks += stop_time_here - start_time_here; Tstate state_to_descend; state_to_descend = ScanState; // for checking @@ -353,7 +353,7 @@ namespace JSC { if (whichDSF == 'B') { // We scan over symmetric states. Only types 14 down to 9 are allowed. if (type_required >= 9 && type_required <= 11) { // iK stepped down on rightIx2; step further up or down - allowed[9] = true; allowed[10] = true; allowed[11] = true; + allowed[9] = true; allowed[10] = true; allowed[11] = true; allowed[12] = true; allowed[13] = true; allowed[14] = true; } else if (type_required >= 12 && type_required <= 14) { // iK stepped up on rightIx2; only step further up @@ -377,7 +377,7 @@ namespace JSC { allowed[12] = (iKexc < iKmax); } // The others are just copies of the ones above: - allowed[1] = allowed[0]; allowed[2] = allowed[0]; allowed[3] = allowed[0]; allowed[4] = allowed[0]; allowed[5] = allowed[0]; allowed[6] = allowed[0]; allowed[7] = allowed[0]; allowed[8] = allowed[0]; + allowed[1] = allowed[0]; allowed[2] = allowed[0]; allowed[3] = allowed[0]; allowed[4] = allowed[0]; allowed[5] = allowed[0]; allowed[6] = allowed[0]; allowed[7] = allowed[0]; allowed[8] = allowed[0]; allowed[10] = allowed[9]; allowed[11] = allowed[9]; allowed[13] = allowed[12]; allowed[14] = allowed[12]; } @@ -400,7 +400,7 @@ namespace JSC { //if (iKmin != iKmax && (ScanState.iK - AveragingState.iK > iKmax || ScanState.iK - AveragingState.iK < iKmin)) data_value = 0.0; //if (abs(data_value) * (type_required_here != 2 ? 1.0 : ph_cost) > running_scan_threshold - //if ((abs(data_value) > running_scan_threshold + //if ((abs(data_value) > running_scan_threshold //|| Nr_ph_Recombinations_Possible (ScanState.label, BaseScanState, type_required_here) > 0) //DP expected_abs_data_value = abs(data_value)/pow(ph_cost, DP(Nr_ph_Recombinations_Possible (ScanState.label, BaseScanState, type_required_here))); @@ -414,14 +414,14 @@ namespace JSC { if ((type_required_here == 11 || type_required_here == 8 || type_required_here == 5 || type_required_here == 2) && Expect_ph_Recombination_iK_Down (ScanState.label, BaseScanState)) expected_abs_data_value /= ph_cost; if (type_required_here == 9 || type_required_here == 6 || type_required_here == 3 || type_required_here == 0) - expected_abs_data_value *= ph_cost; + expected_abs_data_value *= ph_cost; //cout << "\tIncluding thread " << expected_abs_data_value << "\t" << ScanState.label << "\t" << type_required_here << endl; //paused_thread_set.Include_Thread (expected_abs_data_value, ScanState.label, type_required_here); paused_thread_data.Include_Thread (expected_abs_data_value, ScanState.label, type_required_here); //cout << "\tDone including thread." << endl; } // for type_required_here - + //cout << "\tFinished with descendent " << idesc << " out of " << desc_label.size() << " with label " << desc_label[idesc] << endl; //cout << "\tfrom state with label " << label_here << endl; } // for idesc @@ -433,7 +433,7 @@ namespace JSC { */ template - Scan_Info General_Scan (char whichDSF, int iKmin, int iKmax, int iKmod, DP kBT, Tstate& AveragingState, Tstate& SeedScanState, + Scan_Info General_Scan (char whichDSF, int iKmin, int iKmax, int iKmod, DP kBT, Tstate& AveragingState, Tstate& SeedScanState, string defaultScanStatename, int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect rank, Vect nr_processors) { // Performs the scan over excited states, writing data to file. @@ -472,14 +472,14 @@ namespace JSC { //clock_t start_time_local = clock(); stringstream filenameprefix; - Data_File_Name (filenameprefix, whichDSF, iKmin, iKmax, kBT, AveragingState, SeedScanState, defaultScanStatename); + Data_File_Name (filenameprefix, whichDSF, iKmin, iKmax, kBT, AveragingState, SeedScanState, defaultScanStatename); if (in_parallel) for (int r = 0; r < paralevel; ++r) filenameprefix << "_" << rank[r] << "_" << nr_processors[r]; - string prefix = filenameprefix.str(); + string prefix = filenameprefix.str(); stringstream filenameprefix_prevparalevel; // without the rank and nr_processors of the highest paralevel - Data_File_Name (filenameprefix_prevparalevel, whichDSF, iKmin, iKmax, kBT, AveragingState, SeedScanState, defaultScanStatename); + Data_File_Name (filenameprefix_prevparalevel, whichDSF, iKmin, iKmax, kBT, AveragingState, SeedScanState, defaultScanStatename); if (in_parallel) for (int r = 0; r < paralevel - 1; ++r) filenameprefix << "_" << rank[r] << "_" << nr_processors[r]; string prefix_prevparalevel = filenameprefix_prevparalevel.str(); @@ -526,48 +526,48 @@ namespace JSC { else RAW_outfile.open(RAW_Cstr, fstream::out | fstream::app); if (RAW_outfile.fail()) { cout << RAW_Cstr << endl; - JSCerror("Could not open RAW_outfile... "); + JSCerror("Could not open RAW_outfile... "); } RAW_outfile.precision(16); fstream INADM_outfile; if (!refine || in_parallel) INADM_outfile.open(INADM_Cstr, fstream::out | fstream::trunc); else INADM_outfile.open(INADM_Cstr, fstream::out | fstream::app); - if (INADM_outfile.fail()) JSCerror("Could not open INADM_outfile... "); + if (INADM_outfile.fail()) JSCerror("Could not open INADM_outfile... "); INADM_outfile.precision(16); fstream CONV0_outfile; if (!refine || in_parallel) CONV0_outfile.open(CONV0_Cstr, fstream::out | fstream::trunc); else CONV0_outfile.open(CONV0_Cstr, fstream::out | fstream::app); - if (CONV0_outfile.fail()) JSCerror("Could not open CONV0_outfile... "); + if (CONV0_outfile.fail()) JSCerror("Could not open CONV0_outfile... "); CONV0_outfile.precision(16); fstream STAT_outfile; if (!refine || in_parallel) STAT_outfile.open(STAT_Cstr, fstream::out | fstream::trunc); else STAT_outfile.open(STAT_Cstr, fstream::out | fstream::app); - if (STAT_outfile.fail()) JSCerror("Could not open STAT_outfile... "); + if (STAT_outfile.fail()) JSCerror("Could not open STAT_outfile... "); STAT_outfile.precision(8); ofstream LOG_outfile; if (!in_parallel) { if (!refine) LOG_outfile.open(LOG_Cstr, fstream::out | fstream::trunc); else LOG_outfile.open(LOG_Cstr, fstream::out | fstream::app); - if (LOG_outfile.fail()) JSCerror("Could not open LOG_outfile... "); + if (LOG_outfile.fail()) JSCerror("Could not open LOG_outfile... "); LOG_outfile.precision(16); } else { // in_parallel LOG_outfile.open(LOG_Cstr, fstream::out | fstream::trunc); - if (LOG_outfile.fail()) JSCerror("Could not open LOG_outfile... "); + if (LOG_outfile.fail()) JSCerror("Could not open LOG_outfile... "); LOG_outfile.precision(16); //LOG_outfile << endl; } Scan_Info scan_info; - + //Scan_Thread_Set paused_thread_set; //Scan_Thread_Set paused_thread_set_this_run; if (!refine) mkdir(THRDIR_string.c_str(), S_IRWXU | S_IRWXG | S_IRWXO); - Scan_Thread_Data paused_thread_data (THRDIR_string, refine); + Scan_Thread_Data paused_thread_data (THRDIR_string, refine); if (refine) { //paused_thread_set.Load(THR_Cstr); @@ -580,7 +580,7 @@ namespace JSC { Scan_Info scan_info_obtained_in_descent; Scan_State_List ScanStateList (whichDSF, SeedScanState); - ScanStateList.Populate_List(whichDSF, SeedScanState); + ScanStateList.Populate_List(whichDSF, SeedScanState); if (refine && !in_parallel) ScanStateList.Load_Info (SUM_Cstr); else if (in_parallel && rank.sum() == 0) {}; // do nothing, keep the info in the higher .sum file! @@ -608,7 +608,7 @@ namespace JSC { int omp_thread_nr = omp_get_thread_num(); - if ((paused_thread_data.lowest_il_with_nthreads_neq_0 == paused_thread_data.nlists - 1) + if ((paused_thread_data.lowest_il_with_nthreads_neq_0 == paused_thread_data.nlists - 1) && omp_thread_nr > 0) { double start_time_wait = omp_get_wtime(); //cout << "omp_thread " << omp_thread_nr << " sleeping for 5 seconds... " << endl; @@ -640,7 +640,7 @@ namespace JSC { ScanStateList.Raise_Scanning_Flags (exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0)); //cout << "flags: " << endl << ScanStateList.flag_for_scan << endl; - + // Get these base/type started: for (int i = 0; i < ScanStateList.ndef; ++i) { @@ -652,7 +652,7 @@ namespace JSC { ScanStateList.scan_attempted[i] = true; - Tstate ScanState; + Tstate ScanState; //scan_info_before_descent = scan_info; @@ -665,18 +665,18 @@ namespace JSC { bool admissible = ScanState.Check_Admissibility(whichDSF); if (admissible) { - + ScanState.Compute_All(true); - + if (ScanState.conv) { // Put momentum in fundamental window, if possible: int iKexc = ScanState.iK - AveragingState.iK; while (iKexc > iKmax && iKexc - iKmod >= iKmin) iKexc -= iKmod; while (iKexc < iKmin && iKexc + iKmod <= iKmax) iKexc += iKmod; - + //if (iKexc >= iKmin && iKexc <= iKmax) RAW_outfile << endl; - + //data_value = Compute_Matrix_Element_Contrib (whichDSF, iKmin, iKmax, ScanState, AveragingState, Chem_Pot, RAW_outfile); stringstream rawfile_entry; data_value = Compute_Matrix_Element_Contrib (whichDSF, iKmin, iKmax, ScanState, AveragingState, Chem_Pot, rawfile_entry); @@ -694,9 +694,9 @@ namespace JSC { } } //cout << data_value * sumrule_factor << endl; - + // If we force descent: modify data_value by hand so that descent is forced on next scanning pass - //if (Force_Descent(whichDSF, ScanState, AveragingState, iKmod, Chem_Pot) && ScanState.iK - AveragingState.iK < iKmax && Sca nState.iK - AveragingState.iK > iKmin) + //if (Force_Descent(whichDSF, ScanState, AveragingState, iKmod, Chem_Pot) && ScanState.iK - AveragingState.iK < iKmax && Sca nState.iK - AveragingState.iK > iKmin) //if (Force_Descent(whichDSF, ScanState, AveragingState, iKmod, Chem_Pot)) for (int itype = 0; itype < 15; ++itype) { DP data_value_used = 0.1* exp(-paused_thread_data.logscale * JSC::min(0, paused_thread_data.lowest_il_with_nthreads_neq_0)); @@ -707,12 +707,12 @@ namespace JSC { // ++G_7 logic Vect allowed(false, 15); if (whichDSF == 'B') { // symmetric state scanning - allowed[9] = true; allowed[10] = true; allowed[11] = true; + allowed[9] = true; allowed[10] = true; allowed[11] = true; allowed[12] = true; allowed[13] = true; allowed[14] = true; } else { allowed[0] = (iKexc >= iKmin && iKexc <= iKmax); - allowed[1] = allowed[0]; allowed[2] = allowed[0]; allowed[3] = allowed[0]; allowed[4] = allowed[0]; allowed[5] = allowed[0]; allowed[6] = allowed[0]; allowed[7] = allowed[0]; allowed[8] = allowed[0]; + allowed[1] = allowed[0]; allowed[2] = allowed[0]; allowed[3] = allowed[0]; allowed[4] = allowed[0]; allowed[5] = allowed[0]; allowed[6] = allowed[0]; allowed[7] = allowed[0]; allowed[8] = allowed[0]; //allowed[9] = (iKexc <= 0 && iKexc > iKmin); allowed[9] = (iKexc > iKmin); allowed[10] = allowed[9]; allowed[11] = allowed[9]; @@ -730,18 +730,18 @@ namespace JSC { } } } // if (ScanState.conv) - + else { - if (nconv0++ < 1000) - CONV0_outfile << setw(25) << ScanState.label << setw(25) << ScanState.diffsq << setw(5) << ScanState.Check_Rapidities() + if (nconv0++ < 1000) + CONV0_outfile << setw(25) << ScanState.label << setw(25) << ScanState.diffsq << setw(5) << ScanState.Check_Rapidities() << setw(25) << ScanState.String_delta() << endl; scan_info_flags.Ndata++; scan_info_flags.Ndata_conv0++; } } // if admissible - + else { // if inadmissible, modify data_value by hand so that descent is forced on next scanning pass - + if (ninadm++ < 10000000) INADM_outfile << ScanState.label << endl; scan_info_flags.Ndata++; scan_info_flags.Ninadm++; @@ -763,12 +763,12 @@ namespace JSC { Vect allowed(false, 15); if (whichDSF == 'B') { // We scan over symmetric states. Only types 14 down to 9 are allowed. - allowed[9] = true; allowed[10] = true; allowed[11] = true; + allowed[9] = true; allowed[10] = true; allowed[11] = true; allowed[12] = true; allowed[13] = true; allowed[14] = true; } else { allowed[0] = (iKexc >= iKmin && iKexc <= iKmax); - allowed[1] = allowed[0]; allowed[2] = allowed[0]; allowed[3] = allowed[0]; allowed[4] = allowed[0]; allowed[5] = allowed[0]; allowed[6] = allowed[0]; allowed[7] = allowed[0]; allowed[8] = allowed[0]; + allowed[1] = allowed[0]; allowed[2] = allowed[0]; allowed[3] = allowed[0]; allowed[4] = allowed[0]; allowed[5] = allowed[0]; allowed[6] = allowed[0]; allowed[7] = allowed[0]; allowed[8] = allowed[0]; //allowed[9] = (iKexc <= 0 && iKexc > iKmin); allowed[9] = (iKexc > iKmin); allowed[10] = allowed[9]; allowed[11] = allowed[9]; @@ -785,13 +785,13 @@ namespace JSC { } } } // inadmissible - + //scan_info_obtained_in_descent = scan_info; //scan_info_obtained_in_descent -= scan_info_before_descent; scan_info_flags.TT += omp_get_wtime() - start_time_flags; - // Put this info into the appropriate ScanStateList.info + // Put this info into the appropriate ScanStateList.info { #pragma omp critical //ScanStateList.Include_Info(scan_info_obtained_in_descent, ScanStateList.base_label[i]); @@ -801,12 +801,12 @@ namespace JSC { } //cout << "Done with state " << ScanState.label << endl; - + } // if flag_for_scan } // for i - + //clock_t stop_time_flags = clock(); - + //scan_info.CPU_ticks += stop_time_flags - start_time_flags; //scan_info.TT += (stop_time_flags - start_time_flags)/CLOCKS_PER_SEC; //scan_info.TT += omp_get_wtime() - start_time_flags; @@ -837,9 +837,9 @@ namespace JSC { //omp1#pragma omp parallel { - //omp1#pragma omp for + //omp1#pragma omp for for (ithread = 0; ithread < threads_to_do.size(); ++ithread) { - + //cout << "\tithread = " << ithread << endl; //scan_info_before_descent = scan_info; @@ -849,7 +849,7 @@ namespace JSC { //cout << "Thread " << tid << " handling ithread " << ithread << " out of " << threads_to_do.size() << "\t" << threads_to_do[ithread].label << "\t" << threads_to_do[ithread].type << endl; //} - Scan_Info scan_info_this_ithread; + Scan_Info scan_info_this_ithread; double start_time_this_ithread = omp_get_wtime(); // If we don't have time anymore, resave the threads instead of computing them: @@ -860,38 +860,38 @@ namespace JSC { } break; // jump out of ithread loop } - - Tstate ScanState; + + Tstate ScanState; { -#pragma omp critical +#pragma omp critical ScanState = ScanStateList.Return_State(Extract_Base_Label(threads_to_do[ithread].label)); } Tstate BaseScanState; BaseScanState = ScanState; - + //cout << "Setting to label = " << threads_to_do[ithread].label << ", descending type " << threads_to_do[ithread].type << endl; ScanState.Set_to_Label(threads_to_do[ithread].label, BaseScanState.Ix2); //cout << "ScanState after setting label: " << threads_to_do[ithread].label << endl << ScanState << endl; - - + + //cout << "Calling Descend_and_Compute with type " << paused_thread_list.type[ithread] << " on state" << endl << ScanState << endl; /* Descend_and_Compute_for_Fixed_Base (whichDSF, AveragingState, BaseScanState, ScanState, //paused_thread_set.type[ilist][ithread], iKmin, iKmax, iKmod, threads_to_do[ithread].type, iKmin, iKmax, iKmod, - //paused_thread_set_this_run, running_scan_threshold, - paused_thread_data, //thresholdremoved running_scan_threshold, - //paused_thread_set[ilist].abs_data_value[ithread], + //paused_thread_set_this_run, running_scan_threshold, + paused_thread_data, //thresholdremoved running_scan_threshold, + //paused_thread_set[ilist].abs_data_value[ithread], ph_cost, - Max_Secs_used, sumrule_factor, Chem_Pot, scan_info, RAW_outfile, + Max_Secs_used, sumrule_factor, Chem_Pot, scan_info, RAW_outfile, INADM_outfile, ninadm, CONV0_outfile, nconv0, STAT_outfile); - */ - + */ + // STARTING Descend_and_Compute block: int type_required = threads_to_do[ithread].type; - + ScanState.Compute_Momentum(); Vect desc_label; - + // ++G_7 logic bool disperse_only_current_exc_up = false; if (type_required == 14 || type_required == 8 || type_required == 7 || type_required == 6) disperse_only_current_exc_up = true; @@ -901,15 +901,15 @@ namespace JSC { if (type_required == 11 || type_required == 8 || type_required == 5 || type_required == 2) disperse_only_current_exc_down = true; bool preserve_nexc_down = false; if (type_required == 10 || type_required == 7 || type_required == 4 || type_required == 1) preserve_nexc_down = true; - + if (whichDSF == 'B') { // symmetric state scanning if (type_required >= 9 && type_required <= 11) desc_label = Descendent_States_with_iK_Stepped_Down_rightIx2only (ScanState.label, BaseScanState, disperse_only_current_exc_down, preserve_nexc_down); - else if (type_required >= 12 && type_required <= 14) + else if (type_required >= 12 && type_required <= 14) desc_label = Descendent_States_with_iK_Stepped_Up_rightIx2only (ScanState.label, BaseScanState, disperse_only_current_exc_up, preserve_nexc_up); } else { - if (type_required >= 0 && type_required <= 8) { + if (type_required >= 0 && type_required <= 8) { desc_label = Descendent_States_with_iK_Preserved(ScanState.label, BaseScanState, disperse_only_current_exc_up, preserve_nexc_up, disperse_only_current_exc_down, preserve_nexc_down); } else if (type_required >= 9 && type_required <= 11) @@ -917,46 +917,46 @@ namespace JSC { else if (type_required >= 12 && type_required <= 14) desc_label = Descendent_States_with_iK_Stepped_Up (ScanState.label, BaseScanState, disperse_only_current_exc_up, preserve_nexc_up); } - + string label_here = ScanState.label; //int ScanState_iK = ScanState.iK; - + for (int idesc = 0; idesc < desc_label.size(); ++idesc) { - + //clock_t start_time_here = clock(); - + ScanState.Set_to_Label (desc_label[idesc], BaseScanState.Ix2); - + bool admissible = ScanState.Check_Admissibility(whichDSF); - + DP data_value = 0.0; - + //scan_info.Ndata++; //scan_info_this_ithread.Ndata++; - + ScanState.conv = false; ScanState.Compute_Momentum(); // since momentum is used as forced descent criterion - - + + if (admissible) { - + ScanState.Compute_All (idesc == 0); //ScanState.Compute_All (true); - + //scan_info.Ndata++; - + if (ScanState.conv) { //scan_info_this_ithread.Ndata_conv++; - + // Put momentum in fundamental window, if possible: int iKexc = ScanState.iK - AveragingState.iK; while (iKexc > iKmax && iKexc - iKmod >= iKmin) iKexc -= iKmod; while (iKexc < iKmin && iKexc + iKmod <= iKmax) iKexc += iKmod; - + //data_value = Compute_Matrix_Element_Contrib (whichDSF, iKmin, iKmax, ScanState, AveragingState, Chem_Pot, RAW_outfile); stringstream rawfile_entry; data_value = Compute_Matrix_Element_Contrib (whichDSF, iKmin, iKmax, ScanState, AveragingState, Chem_Pot, rawfile_entry); - + { #pragma omp critical RAW_outfile << rawfile_entry.str(); @@ -966,17 +966,17 @@ namespace JSC { scan_info_this_ithread.sumrule_obtained += data_value*sumrule_factor; } } - + //if (iKexc >= iKmin && iKexc <= iKmax) scan_info_this_ithread.sumrule_obtained += data_value*sumrule_factor; //cout << "data_value found = " << data_value * sumrule_factor << endl; - + // Uncomment line below if .stat file is desired: //STAT_outfile << setw(20) << label_here << "\t" << setw(5) << type_required << "\t" << setw(16) << std::scientific << running_scan_threshold << "\t" << setw(20) << ScanState.label << "\t" << setw(16) << data_value << "\t" << setw(16) << std::fixed << setprecision(8) << data_value/running_scan_threshold << endl; - + } // if (ScanState.conv) else { - if (nconv0++ < 1000) - CONV0_outfile << setw(25) << ScanState.label << setw(25) << ScanState.diffsq << setw(5) << ScanState.Check_Rapidities() + if (nconv0++ < 1000) + CONV0_outfile << setw(25) << ScanState.label << setw(25) << ScanState.diffsq << setw(5) << ScanState.Check_Rapidities() << setw(25) << ScanState.String_delta() << endl; //scan_info.Ndata_conv0++; scan_info_this_ithread.Ndata++; @@ -984,7 +984,7 @@ namespace JSC { //cout << "State did not converge." << endl; } } // if (admissible) - + else { if (ninadm++ < 1000000) INADM_outfile << ScanState.label << endl; //scan_info.Ninadm++; @@ -994,22 +994,22 @@ namespace JSC { // Set data_value to enable continued scanning later on: //thresholdremoved data_value = 0.1* running_scan_threshold; } - + //clock_t stop_time_here = clock(); - + //scan_info.CPU_ticks += stop_time_here - start_time_here; //scan_info_this_ithread.CPU_ticks += stop_time_here - start_time_here; //scan_info_this_ithread.TT += (stop_time_here - start_time_here)/CLOCKS_PER_SEC; - + Tstate state_to_descend; state_to_descend = ScanState; // for checking - + ScanState.Compute_Momentum(); // Put momentum in fundamental window, if possible: int iKexc = ScanState.iK - AveragingState.iK; while (iKexc > iKmax && iKexc - iKmod >= iKmin) iKexc -= iKmod; while (iKexc < iKmin && iKexc + iKmod <= iKmax) iKexc += iKmod; - + // ++G_7 logic // Momentum-preserving are only descended to momentum-preserving. // Momentum-increasing are only descended to momentum-preserving and momentum-increasing. @@ -1018,7 +1018,7 @@ namespace JSC { if (whichDSF == 'B') { // We scan over symmetric states. Only types 14 down to 9 are allowed. if (type_required >= 9 && type_required <= 11) { // iK stepped down on rightIx2; step further up or down - allowed[9] = true; allowed[10] = true; allowed[11] = true; + allowed[9] = true; allowed[10] = true; allowed[11] = true; allowed[12] = true; allowed[13] = true; allowed[14] = true; } else if (type_required >= 12 && type_required <= 14) { // iK stepped up on rightIx2; only step further up @@ -1042,19 +1042,19 @@ namespace JSC { allowed[12] = (iKexc < iKmax); } // The others are just copies of the ones above: - allowed[1] = allowed[0]; allowed[2] = allowed[0]; allowed[3] = allowed[0]; allowed[4] = allowed[0]; allowed[5] = allowed[0]; allowed[6] = allowed[0]; allowed[7] = allowed[0]; allowed[8] = allowed[0]; + allowed[1] = allowed[0]; allowed[2] = allowed[0]; allowed[3] = allowed[0]; allowed[4] = allowed[0]; allowed[5] = allowed[0]; allowed[6] = allowed[0]; allowed[7] = allowed[0]; allowed[8] = allowed[0]; allowed[10] = allowed[9]; allowed[11] = allowed[9]; allowed[13] = allowed[12]; allowed[14] = allowed[12]; } - - + + for (int type_required_here = 0; type_required_here < 15; ++type_required_here) { - + if (!allowed[type_required_here]) continue; - + // Reset ScanState to what it was, if change on first pass if (type_required_here > 0) ScanState = state_to_descend; - + // We determine if we carry on scanning based on the data_value obtained, or forcing conditions: // Forcing conditions: //if (!admissible || Force_Descent(whichDSF, ScanState, AveragingState, type_required_here, iKmod, Chem_Pot)) @@ -1063,14 +1063,14 @@ namespace JSC { ////data_value = 1.0; // force for all types of desc // If we're sitting out of the iKmin & iKmax window, stop: //if (iKmin != iKmax && (ScanState.iK - AveragingState.iK > iKmax || ScanState.iK - AveragingState.iK < iKmin)) data_value = 0.0; - + //if (abs(data_value) * (type_required_here != 2 ? 1.0 : ph_cost) > running_scan_threshold - //if ((abs(data_value) > running_scan_threshold + //if ((abs(data_value) > running_scan_threshold //|| Nr_ph_Recombinations_Possible (ScanState.label, BaseScanState, type_required_here) > 0) - + //DP expected_abs_data_value = abs(data_value)/pow(ph_cost, DP(Nr_ph_Recombinations_Possible (ScanState.label, BaseScanState, type_required_here))); DP expected_abs_data_value = abs(data_value); - + //++G_7 logic if ((type_required_here == 14 || type_required_here == 8 || type_required_here == 7 || type_required_here == 6) && Expect_ph_Recombination_iK_Up (ScanState.label, BaseScanState)) expected_abs_data_value /= ph_cost; @@ -1079,8 +1079,8 @@ namespace JSC { if ((type_required_here == 11 || type_required_here == 8 || type_required_here == 5 || type_required_here == 2) && Expect_ph_Recombination_iK_Down (ScanState.label, BaseScanState)) expected_abs_data_value /= ph_cost; if (type_required_here == 9 || type_required_here == 6 || type_required_here == 3 || type_required_here == 0) - expected_abs_data_value *= ph_cost; - + expected_abs_data_value *= ph_cost; + //paused_thread_set.Include_Thread (expected_abs_data_value, ScanState.label, type_required_here); { #pragma omp critical @@ -1089,23 +1089,23 @@ namespace JSC { } //cout << "\tDone including thread." << endl; } // for type_required_here - + //cout << "\tFinished with descendent " << idesc << " out of " << desc_label.size() << " with label " << desc_label[idesc] << endl; //cout << "\tfrom state with label " << label_here << endl; } // for idesc - + //cout << "Finished Descend on state " << endl << ScanState.label << endl; - + // FINISHED Descend_and_Compute block - - - + + + //cout << "Finished descending." << endl; //cout << "\tFinished descending ithread = " << ithread << endl; - + //scan_info_obtained_in_descent = scan_info; //scan_info_obtained_in_descent -= scan_info_before_descent; - //// Put this info into the appropriate ScanStateList.info + //// Put this info into the appropriate ScanStateList.info //ScanStateList.Include_Info(scan_info_obtained_in_descent, Extract_Base_Label(threads_to_do[ithread].label)); scan_info_this_ithread.TT += omp_get_wtime() - start_time_this_ithread; @@ -1115,13 +1115,13 @@ namespace JSC { scan_info += scan_info_this_ithread; //cout << "Including info_ihtread: " << scan_info_this_ithread << endl; ScanStateList.Include_Info(scan_info_this_ithread, Extract_Base_Label(threads_to_do[ithread].label)); - } + } } // for ithread } // omp parallel region // Resynchronize all compute threads: - //omp1#pragma omp barrier + //omp1#pragma omp barrier //cout << "Done with threads_to_do." << endl; @@ -1130,10 +1130,10 @@ namespace JSC { //start_time_local = clock(); /* - if (!in_parallel) - LOG_outfile << "Ndata handled up to now: " << scan_info.Ndata_conv + if (!in_parallel) + LOG_outfile << "Ndata handled up to now: " << scan_info.Ndata_conv << ". Threshold level " << paused_thread_data.lowest_il_with_nthreads_neq_0 << " " << setprecision(6) << exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0) - << ". " << scan_info.Ndata - Ndata_previous_cycle << " new data points. Number of threads: " + << ". " << scan_info.Ndata - Ndata_previous_cycle << " new data points. Number of threads: " << paused_thread_data.nthreads_total.sum() << ". Saturation: " << setprecision(12) << scan_info.sumrule_obtained << endl; */ @@ -1141,10 +1141,10 @@ namespace JSC { //int tid = omp_get_thread_num(); #pragma omp master { - if (!in_parallel) - LOG_outfile << "Master cycling. Ndata_conv " << scan_info.Ndata_conv + if (!in_parallel) + LOG_outfile << "Master cycling. Ndata_conv " << scan_info.Ndata_conv << ". Threshold " << paused_thread_data.lowest_il_with_nthreads_neq_0 << " " << setw(9) << setprecision(3) << exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0) - << ". " << setw(12) << scan_info.Ndata - Ndata_previous_cycle << " new data. Nr of threads: " + << ". " << setw(12) << scan_info.Ndata - Ndata_previous_cycle << " new data. Nr of threads: " << setw(14) << paused_thread_data.nthreads_total.sum() << ". Saturation: " << setprecision(12) << scan_info.sumrule_obtained << endl; @@ -1171,9 +1171,9 @@ namespace JSC { //} while (scan_info.CPU_ticks < ((long long int) Max_Secs_used) * ((long long int) CLOCKS_PER_SEC) //} while (current_time - start_time < ((long long int) Max_Secs_used) * ((long long int) CLOCKS_PER_SEC) - } while (current_time_omp - start_time_omp < Max_Secs_used + } while (current_time_omp - start_time_omp < Max_Secs_used && scan_info.sumrule_obtained < target_sumrule - //&& paused_thread_list.Highest_abs_data_value(0.0, 1.0e+10) > 1.0e-30 + //&& paused_thread_list.Highest_abs_data_value(0.0, 1.0e+10) > 1.0e-30 //&& !(all_threads_zero_previous_cycle && all_threads_zero_current_cycle && !at_least_one_new_flag_raised) //thresholdremoved && running_scan_threshold > 1.0e-10*MACHINE_EPS ); @@ -1193,11 +1193,11 @@ namespace JSC { if (!in_parallel) { - if (scan_info.sumrule_obtained >= target_sumrule) - LOG_outfile << endl << "Achieved sumrule saturation of " << scan_info.sumrule_obtained + if (scan_info.sumrule_obtained >= target_sumrule) + LOG_outfile << endl << "Achieved sumrule saturation of " << scan_info.sumrule_obtained << "\t(target was " << target_sumrule << ")." << endl << endl; - //thresholdremoved + //thresholdremoved //if (running_scan_threshold < MACHINE_EPS) //LOG_outfile << endl << "Stopping because threshold lower than machine precision. " << endl << endl; @@ -1218,12 +1218,12 @@ namespace JSC { else { // in_parallel - //thresholdremoved + //thresholdremoved //if (running_scan_threshold < MACHINE_EPS) //LOG_outfile << "rank " << rank << " out of " << nr_processors << " processors: " // << "Stopping because threshold lower than machine precision. " << endl << endl; - LOG_outfile << "rank " << rank << " out of " << nr_processors << " processors: " + LOG_outfile << "rank " << rank << " out of " << nr_processors << " processors: " //thresholdremoved << "run info: " << scan_info << endl << "Latest running_scan_threshold = " << running_scan_threshold << endl; << "run info: " << scan_info << endl << "Latest threshold = " << exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0) << endl; } @@ -1277,13 +1277,13 @@ namespace JSC { // whichDSF == 'o': one-body function < \Psi^{\dagger} \Psi > // Delta is the number of sites involved in the smoothing of the entropy - //int Delta = int(sqrt(N))/2;//6;//N/20; + //int Delta = int(sqrt(N))/2;//6;//N/20; //DP epsilon = log(L)/L; // using Gaussian for density in entropy. //DP epsilon = 1.0/L; // using Lorentzian for density in entropy. // Construct the finite-size saddle-point state: // if we refine, read the quantum numbers of the saddle point state (and seed sps) from the sps file: - + stringstream SPS_stringstream; string SPS_string; //SPS_stringstream << "Tgt0_"; //Data_File_Name (SPS_stringstream, whichDSF, iKmin, iKmax, kBT, spstate, SeedScanState, ""); @@ -1292,7 +1292,7 @@ namespace JSC { SPS_string = SPS_stringstream.str(); const char* SPS_Cstr = SPS_string.c_str(); fstream spsfile; - if (refine) spsfile.open(SPS_Cstr, fstream::in); + if (refine) spsfile.open(SPS_Cstr, fstream::in); else spsfile.open(SPS_Cstr, fstream::out | fstream::trunc); if (spsfile.fail()) { cout << SPS_Cstr << endl; JSCerror("Could not open spsfile."); @@ -1306,7 +1306,7 @@ namespace JSC { else { // read it from the sps file // Check that the sps has the right number of Ix2: int Nspsread; - spsfile >> Nspsread; + spsfile >> Nspsread; if (Nspsread != N) { cout << Nspsread << "\t" << N << endl; JSCerror("Wrong number of Ix2 in saddle-point state."); @@ -1338,7 +1338,7 @@ namespace JSC { else { // read it from the sps file // Check that the sps has the right number of Ix2: int Nsspsread; - spsfile >> Nsspsread; + spsfile >> Nsspsread; if (Nsspsread != Nscan) { cout << Nsspsread << "\t" << Nscan << endl; JSCerror("Wrong number of Ix2 in scan saddle-point state."); @@ -1363,11 +1363,11 @@ namespace JSC { spsfile << endl << spstate << endl << endl; for (int i = 1; i < spstate.N - 2; ++i) - spsfile << 0.5 * (spstate.lambdaoc[i] + spstate.lambdaoc[i+1]) + spsfile << 0.5 * (spstate.lambdaoc[i] + spstate.lambdaoc[i+1]) //<< "\t" << twoPI/(spstate.L * (spstate.lambdaoc[i+1] - spstate.lambdaoc[i])) << endl; - << "\t" << 1.0/spstate.L * (0.25/(spstate.lambdaoc[i] - spstate.lambdaoc[i-1]) - + 0.5/(spstate.lambdaoc[i+1] - spstate.lambdaoc[i]) - + 0.25/(spstate.lambdaoc[i+2] - spstate.lambdaoc[i+1])) + << "\t" << 1.0/spstate.L * (0.25/(spstate.lambdaoc[i] - spstate.lambdaoc[i-1]) + + 0.5/(spstate.lambdaoc[i+1] - spstate.lambdaoc[i]) + + 0.25/(spstate.lambdaoc[i+2] - spstate.lambdaoc[i+1])) << "\t" << rho_of_lambdaoc_1 (spstate, 0.5 * (spstate.lambdaoc[i] + spstate.lambdaoc[i+1]), delta) << "\t" << rho_of_lambdaoc_2 (spstate, 0.5 * (spstate.lambdaoc[i] + spstate.lambdaoc[i+1]), delta) << endl; @@ -1394,7 +1394,7 @@ namespace JSC { // Scanning on an excited state defined by a set of Ix2: - void Scan_LiebLin (char whichDSF, LiebLin_Bethe_State AveragingState, string defaultScanStatename, int iKmin, int iKmax, + void Scan_LiebLin (char whichDSF, LiebLin_Bethe_State AveragingState, string defaultScanStatename, int iKmin, int iKmax, int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect rank, Vect nr_processors) { // This function is as Scan_LiebLin for generic T defined above, except that the @@ -1427,14 +1427,14 @@ namespace JSC { //if (whichDSF == 'g') { //for (int i = 0; i < N; ++i) SeedScanState.Ix2[i] = AveragingState.Ix2[i] - 1; //SeedScanState.Ix2[N] = SeedScanState.Ix2[N-1] + 2; - //} + //} // If 'o', remove midmost and shift quantum numbers by half-integer towards removed one: if (whichDSF == 'o') { for (int i = 0; i < N-1; ++i) SeedScanState.Ix2[i] = AveragingState.Ix2[i + (i >= N/2)] + 1 - 2*(i >= N/2); } - // If 'g', add a quantum number in middle (explicitly: to right of index N/2) + // If 'g', add a quantum number in middle (explicitly: to right of index N/2) // and shift quantum numbers by half-integer away from added one: if (whichDSF == 'g') { SeedScanState.Ix2[N/2] = AveragingState.Ix2[N/2] - 1; @@ -1459,7 +1459,7 @@ namespace JSC { } // Simplified function call of the above: - void Scan_LiebLin (char whichDSF, LiebLin_Bethe_State AveragingState, string defaultScanStatename, int iKmin, int iKmax, + void Scan_LiebLin (char whichDSF, LiebLin_Bethe_State AveragingState, string defaultScanStatename, int iKmin, int iKmax, int Max_Secs, DP target_sumrule, bool refine) { int paralevel = 0; @@ -1472,7 +1472,7 @@ namespace JSC { } // Scanning on a previously-defined AveragingState - void Scan_Heis (char whichDSF, XXZ_Bethe_State& AveragingState, string defaultScanStatename, int iKmin, int iKmax, + void Scan_Heis (char whichDSF, XXZ_Bethe_State& AveragingState, string defaultScanStatename, int iKmin, int iKmax, int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect rank, Vect nr_processors) { // General state scanning for Heisenberg chains @@ -1487,7 +1487,7 @@ namespace JSC { if (whichDSF == 'Z' || whichDSF == 'z') SeedScanState = AveragingState; else if (whichDSF == 'm') SeedScanState = Remove_Particle_at_Center (AveragingState); else if (whichDSF == 'p') SeedScanState = Add_Particle_at_Center (AveragingState); - else JSCerror("Unknown whichDSF in Scan_Heis."); + else JSCerror("Unknown whichDSF in Scan_Heis."); //cout << "In General_Scan: SeedScanState = " << SeedScanState << endl; @@ -1497,7 +1497,7 @@ namespace JSC { } // Scanning on a previously-defined AveragingState - void Scan_Heis (char whichDSF, XXX_Bethe_State& AveragingState, string defaultScanStatename, int iKmin, int iKmax, + void Scan_Heis (char whichDSF, XXX_Bethe_State& AveragingState, string defaultScanStatename, int iKmin, int iKmax, int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect rank, Vect nr_processors) { // General state scanning for Heisenberg chains @@ -1512,7 +1512,7 @@ namespace JSC { if (whichDSF == 'Z' || whichDSF == 'z') SeedScanState = AveragingState; else if (whichDSF == 'm') SeedScanState = Remove_Particle_at_Center (AveragingState); else if (whichDSF == 'p') SeedScanState = Add_Particle_at_Center (AveragingState); - else JSCerror("Unknown whichDSF in Scan_Heis."); + else JSCerror("Unknown whichDSF in Scan_Heis."); // Now the scan itself General_Scan (whichDSF, iKmin, iKmax, AveragingState.chain.Nsites, 0.0, AveragingState, SeedScanState, defaultScanStatename, Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors); @@ -1520,7 +1520,7 @@ namespace JSC { } // Scanning on a previously-defined AveragingState - void Scan_Heis (char whichDSF, XXZ_gpd_Bethe_State& AveragingState, string defaultScanStatename, int iKmin, int iKmax, + void Scan_Heis (char whichDSF, XXZ_gpd_Bethe_State& AveragingState, string defaultScanStatename, int iKmin, int iKmax, int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect rank, Vect nr_processors) { // General state scanning for Heisenberg chains @@ -1535,7 +1535,7 @@ namespace JSC { if (whichDSF == 'Z' || whichDSF == 'z') SeedScanState = AveragingState; else if (whichDSF == 'm') SeedScanState = Remove_Particle_at_Center (AveragingState); else if (whichDSF == 'p') SeedScanState = Add_Particle_at_Center (AveragingState); - else JSCerror("Unknown whichDSF in Scan_Heis."); + else JSCerror("Unknown whichDSF in Scan_Heis."); // Now the scan itself General_Scan (whichDSF, iKmin, iKmax, AveragingState.chain.Nsites, 0.0, AveragingState, SeedScanState, defaultScanStatename, Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors); @@ -1543,7 +1543,7 @@ namespace JSC { } - //void Scan_Heis (char whichDSF, DP Delta, DP N, int M, bool fixed_iK, int iKneeded, + //void Scan_Heis (char whichDSF, DP Delta, DP N, int M, bool fixed_iK, int iKneeded, void Scan_Heis (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax, int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect rank, Vect nr_processors) { @@ -1554,48 +1554,48 @@ namespace JSC { // whichDSF == 'm': S^{-+} // whichDSF == 'z': S^{zz} // whichDSF == 'p': S^{+-} - // whichDSF == 'a': < S^z_j S^z_{j+1} S^z_l S^z_{l+1} > for RIXS - // whichDSF == 'b': < S^z_j S^-_{j+1} S^-_l S^z_{l+1} > + (m <-> z) for RIXS - // whichDSF == 'c': < S^-_j S^-_{j+1} S^-_l S^-_{l+1} > for RIXS + // whichDSF == 'a': < S^z_j S^z_{j+1} S^z_l S^z_{l+1} > for RIXS + // whichDSF == 'b': < S^z_j S^-_{j+1} S^-_l S^z_{l+1} > + (m <-> z) for RIXS + // whichDSF == 'c': < S^-_j S^-_{j+1} S^-_l S^-_{l+1} > for RIXS Heis_Chain BD1(1.0, Delta, 0.0, N); Vect_INT Nrapidities_groundstate(0, BD1.Nstrings); - + Nrapidities_groundstate[0] = M; - + Heis_Base baseconfig_groundstate(BD1, Nrapidities_groundstate); - + if ((Delta > 0.0) && (Delta < 1.0)) { - + XXZ_Bethe_State GroundState(BD1, baseconfig_groundstate); GroundState.Compute_All(true); - // The ground state is now fully defined. + // The ground state is now fully defined. XXZ_Bethe_State SeedScanState; if (whichDSF == 'Z' || whichDSF == 'z') SeedScanState = GroundState; else if (whichDSF == 'm') SeedScanState = XXZ_Bethe_State(GroundState.chain, M - 1); else if (whichDSF == 'p') SeedScanState = XXZ_Bethe_State(GroundState.chain, M + 1); - else JSCerror("Unknown whichDSF in Scan_Heis."); + else JSCerror("Unknown whichDSF in Scan_Heis."); // Now the scan itself General_Scan (whichDSF, iKmin, iKmax, N, 0.0, GroundState, SeedScanState, "", Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors); } - + else if (Delta == 1.0) { XXX_Bethe_State GroundState(BD1, baseconfig_groundstate); GroundState.Compute_All(true); - // The ground state is now fully defined. + // The ground state is now fully defined. XXX_Bethe_State SeedScanState; if (whichDSF == 'Z' || whichDSF == 'z' || whichDSF == 'a' || whichDSF == 'q') SeedScanState = GroundState; else if (whichDSF == 'm') SeedScanState = XXX_Bethe_State(GroundState.chain, M - 1); else if (whichDSF == 'p') SeedScanState = XXX_Bethe_State(GroundState.chain, M + 1); else if (whichDSF == 'c') SeedScanState = XXX_Bethe_State(GroundState.chain, M - 2); - else JSCerror("Unknown whichDSF in Scan_Heis."); + else JSCerror("Unknown whichDSF in Scan_Heis."); // Now the scan itself General_Scan (whichDSF, iKmin, iKmax, N, 0.0, GroundState, SeedScanState, "", Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors); @@ -1605,13 +1605,13 @@ namespace JSC { XXZ_gpd_Bethe_State GroundState(BD1, baseconfig_groundstate); GroundState.Compute_All(true); - // The ground state is now fully defined. + // The ground state is now fully defined. XXZ_gpd_Bethe_State SeedScanState; if (whichDSF == 'Z' || whichDSF == 'z') SeedScanState = GroundState; else if (whichDSF == 'm') SeedScanState = XXZ_gpd_Bethe_State(GroundState.chain, M - 1); else if (whichDSF == 'p') SeedScanState = XXZ_gpd_Bethe_State(GroundState.chain, M + 1); - else JSCerror("Unknown whichDSF in Scan_Heis."); + else JSCerror("Unknown whichDSF in Scan_Heis."); // Now the scan itself General_Scan (whichDSF, iKmin, iKmax, N, 0.0, GroundState, SeedScanState, "", Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors); @@ -1666,20 +1666,20 @@ namespace JSC { // whichDSF == 'm': S^{-+} // whichDSF == 'z': S^{zz} // whichDSF == 'p': S^{+-} - // whichDSF == 'a': < S^z_j S^_{j+1} S^z_l S^z_{l+1} > for RIXS + // whichDSF == 'a': < S^z_j S^_{j+1} S^z_l S^z_{l+1} > for RIXS Heis_Chain BD1(1.0, Delta, 0.0, N); Vect_INT Nrapidities_groundstate(0, BD1.Nstrings); - + Nrapidities_groundstate[0] = M; - + ODSLF_Base baseconfig_groundstate(BD1, Nrapidities_groundstate); - + ODSLF_Ix2_Offsets baseoffsets(baseconfig_groundstate, 0ULL); if ((Delta > 0.0) && (Delta < 1.0)) { - + ODSLF_XXZ_Bethe_State GroundState(BD1, baseconfig_groundstate); GroundState.Compute_All(true); @@ -1689,7 +1689,7 @@ namespace JSC { } */ - /* + /* else if (Delta == 1.0) { XXX_Bethe_State GroundState(BD1, baseconfig_groundstate); @@ -1738,9 +1738,9 @@ namespace JSC { - // Geometric quenches + // Geometric quenches - void Scan_LiebLin_Geometric_Quench (DP c_int, DP L_1, int type_id_1, long long int id_1, DP L_2, int N, int iK_UL, + void Scan_LiebLin_Geometric_Quench (DP c_int, DP L_1, int type_id_1, long long int id_1, DP L_2, int N, int iK_UL, int Max_Secs, DP target_sumrule, bool refine) { // We decompose the wavefunction of state 1 (living on length L_1) into @@ -1768,7 +1768,7 @@ namespace JSC { } - void Scan_Heis_Geometric_Quench (DP Delta, int N_1, int M, long long int base_id_1, long long int type_id_1, long long int id_1, + void Scan_Heis_Geometric_Quench (DP Delta, int N_1, int M, long long int base_id_1, long long int type_id_1, long long int id_1, int N_2, int iKmin, int iKmax, int Max_Secs, DP target_sumrule, bool refine) { // We decompose the wavefunction of state 1 (living on length L_1) into @@ -1781,7 +1781,7 @@ namespace JSC { if ((Delta > 0.0) && (Delta < 1.0)) { JSCerror("Geometric quench not yet implemented for XXZ."); } - + else if (Delta == 1.0) { XXX_Bethe_State BasicState_1(BD_1, base_id_1, type_id_1); diff --git a/src/SCAN/Scan_Thread_List.cc b/src/SCAN/Scan_Thread_List.cc deleted file mode 100644 index e7b97f8..0000000 --- a/src/SCAN/Scan_Thread_List.cc +++ /dev/null @@ -1,466 +0,0 @@ -/********************************************************** - -This software is part of J.-S. Caux's ABACUS library. - -Copyright (c). - ------------------------------------------------------------ - -File: Scan_Thread_List.cc - -Purpose: defines all functions for Scan_Thread_List class. - - -***********************************************************/ - -#include "JSC.h" - -using namespace std; -using namespace JSC; - -namespace JSC { - - Scan_Thread_List::Scan_Thread_List() - { - dim = 2; - nthreads_tot = 0; - nthreads_done = 0; - //omega = Vect(1.0e+6, dim); - //iK = Vect_INT(0, dim); - //abs_data_value = Vect_DP (0.0, Nentries); - abs_data_value = Vect (0.0, dim); - label = Vect (dim); - type = Vect (dim); - isdone = Vect (false, dim); - } - - Scan_Thread_List::Scan_Thread_List (int Nentries) - { - dim = Nentries; - nthreads_tot = 0; - nthreads_done = 0; - ////omega = Vect_DP(1.0e+6, Nentries); - //omega = Vect(1.0e+6, Nentries); - //iK = Vect_INT(0, Nentries); - //abs_data_value = Vect_DP (0.0, Nentries); - abs_data_value = Vect (0.0, Nentries); - label = Vect (Nentries); - type = Vect (Nentries); - isdone = Vect (false, Nentries); - } - - bool Scan_Thread_List::Increase_Size (int nr_to_add) // resizes the vectors to accommodate up to nr_to_add additional entries - { - dim += nr_to_add; - // Nentries stays the same - - //cout << "Called Increase_Size with nr_to_add = " << nr_to_add << endl; - try { - //omega.Increase_Size (nr_to_add); - //iK.Increase_Size (nr_to_add); - abs_data_value.Increase_Size (nr_to_add); - //cout << "Increasing size of label" << endl; - label.Increase_Size (nr_to_add, ""); - //cout << "Done increasing size of label" << endl; - type.Increase_Size (nr_to_add); - isdone.Increase_Size (nr_to_add, false); - } - - catch (bad_alloc) { - cout << "dim " << dim << "\tnr_to_add " << nr_to_add << endl; - JSCerror("Memory allocation failed in Scan_Thread_List::Increase_Size."); - } - - return(true); - } - - //void Scan_Thread_List::Include_Thread (DP omega_ref, int iK_ref, DP abs_data_value_ref, string label_ref, int type_ref) - void Scan_Thread_List::Include_Thread (DP data_value_ref, string label_ref, int type_ref) - { - if (nthreads_tot > dim - 10) { // Resize the Scan_Thread_List, by doubling its size - //cout << "Resizing threads list" << endl; - (*this).Increase_Size (dim); - //cout << "Done resizing" << endl; - } - - //cout << "Including thread for label " << label_ref << " and type " << type_ref << endl; - - //omega[nthreads] = omega_ref; - //iK[nthreads] = iK_ref; - abs_data_value[nthreads_tot] = abs(data_value_ref); - label[nthreads_tot] = label_ref; - type[nthreads_tot] = type_ref; - isdone[nthreads_tot] = false; - - nthreads_tot++; - - if (nthreads_tot >= dim) { - (*this).Save("Threads_stopped.thr"); - JSCerror("nthreads_tot >= dim in Scan_Thread_List."); - } - - } - - void Scan_Thread_List::Order_in_abs_data_value () - { - // Reorders all entries in decreasing order in abs_data_value - - if (nthreads_tot > 1) { - - Vect_INT index(dim); - for (int i = 0; i < nthreads_tot; ++i) index[i] = i; - - abs_data_value.QuickSort(index, 0, nthreads_tot - 1); - - ////Vect_DP omega_ordered(nthreads); - //Vect omega_ordered(nthreads); - //Vect_INT iK_ordered(nthreads); - //Vect_DP abs_data_value_ordered(nthreads); - Vect abs_data_value_ordered(nthreads_tot); - Vect label_ordered(nthreads_tot); - Vect type_ordered(nthreads_tot); - Vect isdone_ordered(nthreads_tot); - - // Put data in proper order - for (int i = 0; i < nthreads_tot; ++i) { - //omega_ordered[i] = omega[index[nthreads - 1 - i] ]; - //iK_ordered[i] = iK[index[nthreads - 1 - i] ]; - abs_data_value_ordered[i] = abs_data_value[nthreads_tot - 1 - i]; - label_ordered[i] = label[index[nthreads_tot - 1 - i] ]; - type_ordered[i] = type[index[nthreads_tot - 1 - i] ]; - //sector_lowest_raisable_ordered[i] = sector_lowest_raisable[index[nthreads - 1 - i] ]; - isdone_ordered[i] = isdone[index[nthreads_tot - 1 - i] ]; - } - - // Put back in *this object: - for (int i = 0; i < nthreads_tot; ++i) { - //omega[i] = omega_ordered[i]; - //iK[i] = iK_ordered[i]; - abs_data_value[i] = abs_data_value_ordered[i]; - label[i] = label_ordered[i]; - type[i] = type_ordered[i]; - isdone[i] = isdone_ordered[i]; - } // The rest are all simply 0. - - // Done ! - } - } - /* - void Scan_Thread_List::Order_in_omega () - { - // Reorders all entries in increasing order in omega - - if (nthreads > 1) { - - Vect_INT index(0, dim); - for (int i = 0; i < nthreads; ++i) index[i] = i; - - omega.QuickSort(index, 0, nthreads - 1); - - //Vect_DP omega_ordered(nthreads); - Vect omega_ordered(nthreads); - Vect_INT iK_ordered(nthreads); - //Vect_DP abs_data_value_ordered(nthreads); - Vect abs_data_value_ordered(nthreads); - Vect label_ordered(nthreads); - Vect type_ordered(nthreads); - - // Put data in proper order - for (int i = 0; i < nthreads; ++i) { - omega_ordered[i] = omega[i]; - iK_ordered[i] = iK[index[i] ]; - abs_data_value_ordered[i] = abs_data_value[index[i] ]; - label_ordered[i] = label[index[i] ]; - type_ordered[i] = type[index[i] ]; - } - - // Put back in *this object: - for (int i = 0; i < nthreads; ++i) { - omega[i] = omega_ordered[i]; - iK[i] = iK_ordered[i]; - abs_data_value[i] = abs_data_value_ordered[i]; - label[i] = label_ordered[i]; - type[i] = type_ordered[i]; - } // The rest are all simply 0. - - // Done ! - } - } - */ - /* - DP Scan_Thread_List::Highest_abs_data_value_ordered (DP omega_MIN, DP omega_MAX) // returns highest abs_data_value in window omega_MIN to omega_MAX - { - // Assume that the threads have been Order_in_abs_data_value - for (int i = 0; i < nthreads; ++i) - if (omega[i] > omega_MIN && omega[i] < omega_MAX) return(abs_data_value[i]); - - return(0.0); - } - */ - //DP Scan_Thread_List::Highest_abs_data_value (DP omega_MIN, DP omega_MAX) // returns highest abs_data_value in window omega_MIN to omega_MAX - DP Scan_Thread_List::Highest_abs_data_value () // returns highest abs_data_value - { - // Works even if list is not ordered - DP maxdatavalue = 0.0; - for (int i = 0; i < nthreads_tot; ++i) - //if (omega[i] > omega_MIN && omega[i] < omega_MAX && abs_data_value[i] > maxdatavalue) - if (!isdone[i] && (abs_data_value[i] > maxdatavalue)) - maxdatavalue = abs_data_value[i]; - - return(maxdatavalue); - } - - DP Scan_Thread_List::kth_highest_abs_data_value (int k) // returns the k-th highest abs_data_value in window omega_MIN to omega_MAX - { - // Works even if list is not ordered - if (k < 1) JSCerror("Give k > 1 in Scan_Thread_List::kth_Highest_abs_data_value."); - else if (k == 1 || k >= nthreads_tot) // Threads list not long enough, return the top one - //return ((*this).Highest_abs_data_value (-1.0e+10, 1.0e+10)); - return ((*this).Highest_abs_data_value ()); - - Vect topk (0.0, k); - - for (int i = 0; i < nthreads_tot; ++i) - if (!isdone[i] && (abs_data_value[i] > topk[0])) { - topk[0] = abs_data_value[i]; - topk.QuickSort(0, k-1); - } - - return(topk[0]); - } - - bool Scan_Thread_List::Exists_data_value_greater_than (DP value) - { - for (int i = 0; i < nthreads_tot; ++i) - if (abs_data_value[i] > value) return(true); - return(false); - } - - /* - DP Scan_Thread_List::Lowest_omega () // returns lowest omega - { - // Assume that the threads have been Order_in_omega - return(omega[0]); - } - */ - void Scan_Thread_List::Merge (const Scan_Thread_List& reflist) - { - if (nthreads_tot + reflist.nthreads_tot >= dim) // JSCerror("Scan_Threads_List: too big to merge."); - //(*this).Increase_Size (reflist.nthreads); - (*this).Increase_Size (nthreads_tot/10 + reflist.nthreads_tot); - - for (int i = 0; i < reflist.nthreads_tot; ++i) { - //omega[nthreads] = reflist.omega[i]; - //iK[nthreads] = reflist.iK[i]; - abs_data_value[nthreads_tot] = reflist.abs_data_value[i]; - label[nthreads_tot] = reflist.label[i]; - type[nthreads_tot] = reflist.type[i]; - isdone[nthreads_tot] = reflist.isdone[i]; - nthreads_tot++; - } - } - - /* // DEACTIVATED IN ABACUS++G_2 - void Scan_Thread_List::Remove_Threads (int ithread_down, int ithread_up) - { - if (ithread_down < 0 || ithread_up > nthreads) { - JSCerror("Trying to remove inexistent entries in Remove_Threads."); - } - - for (int i = 0; i < nthreads - ithread_up - 1; ++i) { - //omega[ithread_down + i] = omega[ithread_up + i + 1]; - //iK[ithread_down + i] = iK[ithread_up + 1 + i]; - abs_data_value[ithread_down + i] = abs_data_value[ithread_up + i + 1]; - label[ithread_down + i] = label[ithread_up + i + 1]; - type[ithread_down + i] = type[ithread_up + i + 1]; - } - // Zero the other entries: - for (int i = nthreads - ithread_up - 1; i < nthreads; ++i) { - //omega[i] = 1.0e+6; - //iK[i] = 0; - abs_data_value[i] = 0.0;//-1.0e-30; // give a small negative value, so these are put at the very bottom of the list - label[i] = ""; - type[i] = 0; - } - - nthreads -= ithread_up - ithread_down + 1; // removed that many entries. - } - */ - /* // DEACTIVATED IN ABACUS++G_2 - int Scan_Thread_List::Remove_Threads (const Vect& threads_done) - { - // If threads_done[ithread] == true, remove from list - - if (threads_done.size() != nthreads) JSCerror("Wrong size boolean vector in Scan_Thread_List::Remove_Threads."); - - int nr_removed = 0; - - for (int i = 0; i < nthreads; ++i) { - if (!threads_done[i]) { - //omega[i - nr_removed] = omega[i]; - //iK[i - nr_removed] = iK[i]; - abs_data_value[i - nr_removed] = abs_data_value[i]; - label[i - nr_removed] = label[i]; - type[i - nr_removed] = type[i]; - - } - else nr_removed++; - } - // Zero other entries: - for (int i = nthreads - nr_removed; i < nthreads; ++i) { - //omega[i] = 1.0e+6; - //iK[i] = 0; - abs_data_value[i] = 0.0; - label[i] = ""; - type[i] = 0; - } - - nthreads -= nr_removed; - - return(nr_removed); - } - */ - - void Scan_Thread_List::Remove_Done_Threads () - { - // If isdone[ithread] == true, remove from list - - int nr_removed = 0; - - for (int i = 0; i < nthreads_tot; ++i) { - if (!isdone[i]) { - //omega[i - nr_removed] = omega[i]; - //iK[i - nr_removed] = iK[i]; - abs_data_value[i - nr_removed] = abs_data_value[i]; - label[i - nr_removed] = label[i]; - type[i - nr_removed] = type[i]; - isdone[i - nr_removed] = isdone[i]; - } - else nr_removed++; - } - // Zero other entries: - for (int i = nthreads_tot - nr_removed; i < nthreads_tot; ++i) { - //omega[i] = 1.0e+6; - //iK[i] = 0; - abs_data_value[i] = 0.0; - label[i] = ""; - type[i] = 0; - isdone[i] = false; - } - - nthreads_tot -= nr_removed; - if (nthreads_done != nr_removed) { - cout << nthreads_done << "\t" << nr_removed << endl; - JSCerror("Miscount in removing threads during Remove_Done_Threads."); - } - nthreads_done = 0; - - return; - } - - void Scan_Thread_List::Clear () - { - // dim is preserved - - nthreads_tot = 0; - nthreads_done = 0; - ////omega = Vect_DP(0.0, dim); - //omega = Vect (1.0e+6, dim); - //iK = Vect_INT(0, dim); - //abs_data_value = Vect_DP (0.0, dim); - abs_data_value = Vect (0.0, dim); - label = Vect (dim); - type = Vect (dim); - isdone = Vect (false, dim); - } - - void Scan_Thread_List::Save(const char* outfile_Cstr) - { - // We save only the undone threads, so after cleanup. - (*this).Remove_Done_Threads(); - - ofstream outfile; - - outfile.open(outfile_Cstr); - if (outfile.fail()) JSCerror("Could not open outfile... "); - outfile.precision(3); - - if (nthreads_tot > 0) { - //outfile << omega[0] << "\t" << iK[0] << "\t" << abs_data_value[0] << "\t" << label[0] << "\t" << type[0]; - outfile << abs_data_value[0] << "\t" << label[0] << "\t" << type[0]; - for (int i = 1; i < nthreads_tot; ++i) { - //outfile << endl << omega[i] << "\t" << iK[i] << "\t" << abs_data_value[i] << "\t" << label[i] << "\t" << type[i]; - outfile << endl << abs_data_value[i] << "\t" << label[i] << "\t" << type[i]; - } - } - - outfile.close(); - - return; - } - - void Scan_Thread_List::Load (const char* thrfile_Cstr) - { - /* - // Check that the data fits within limits: - struct stat statbuf; - - stat (thrfile_Cstr, &statbuf); - int filesize = statbuf.st_size; - - // Determine the number of entries approximately - // Don't overestimate entry_size: pretend label takes no space... - int entry_size = 2* sizeof(float) + sizeof(int); - - int estimate_nr_entries = filesize/entry_size; - */ - - // Get number of lines in the file: - int nrlines = 0; - ifstream infile1(thrfile_Cstr); - if(infile1.fail()) { - cout << "Could not open file " << thrfile_Cstr << " in Scan_Thread_List::Load." << endl; - JSCerror("Terminating."); - } - string dummy; - while (getline(infile1, dummy)) ++nrlines; - infile1.close(); - - int estimate_nr_entries = JSC::max(10, (11*nrlines)/10); - - if (estimate_nr_entries > dim) { - (*this).Increase_Size (estimate_nr_entries - dim/2); // give dim/2 safety - } - - ifstream infile; - infile.open(thrfile_Cstr); - if(infile.fail()) { - cout << "Could not open file " << thrfile_Cstr << " in Scan_Thread_List::Load." << endl; - JSCerror("Terminating."); - } - - (*this).Clear(); - - DP abs_data_value_DP; - while (infile.peek() != EOF) { - //infile >> omega[nthreads]; - //infile >> iK[nthreads]; - //infile >> abs_data_value[nthreads]; // nasty bug: reading into float bugs when value less or more than float limits - infile >> abs_data_value_DP; // cure: first read into double, then cast to float: - abs_data_value[nthreads_tot] = float(abs_data_value_DP); - infile >> label[nthreads_tot]; - infile >> type[nthreads_tot]; - isdone[nthreads_tot] = false; - nthreads_tot++; - if (nthreads_tot >= dim) { - cout << "file " << thrfile_Cstr << endl << "estimate_nr_entries = " << estimate_nr_entries << "\tnrlines = " << nrlines << endl << "nthreads_tot = " << nthreads_tot << "\tdim = " << dim << endl; - JSCerror("Too many threads in input file."); - } - } - - infile.close(); - - return; - } - -} // namespace JSC diff --git a/src/SCAN/Scan_Thread_Set.cc b/src/SCAN/Scan_Thread_Set.cc deleted file mode 100644 index 0c3663e..0000000 --- a/src/SCAN/Scan_Thread_Set.cc +++ /dev/null @@ -1,254 +0,0 @@ -/********************************************************** - -This software is part of J.-S. Caux's ABACUS library. - -Copyright (c). - ------------------------------------------------------------ - -File: Scan_Thread_Set.cc - -Purpose: defines all functions for Scan_Thread_Set class. - - -***********************************************************/ - -#include "JSC.h" - -using namespace std; -using namespace JSC; - -namespace JSC { - - Scan_Thread_Set::Scan_Thread_Set() - { - dim = Vect (nlists); - nthreads_tot = Vect (nlists); - nthreads_done = Vect (nlists); - - label = Vect > (nlists); - type = Vect > (nlists); - isdone = Vect > (nlists); - - - // Give starting values to all: - for (int il = 0; il < nlists; ++il) { - dim[il] = 100; - nthreads_tot[il] = 0; - nthreads_done[il] = 0; - label[il] = Vect (dim[il]); - type[il] = Vect (dim[il]); - isdone[il] = Vect (false, dim[il]); - } - } - - - bool Scan_Thread_Set::Increase_Size (int il, int nr_to_add) // resizes the vectors to accommodate up to nr_to_add additional entries - { - if (il < 0 || il > nlists) JSCerror("ilist out of bounds in Scan_Thread_Set::Increase_Size"); - - dim[il] += nr_to_add; - - try { - label[il].Increase_Size (nr_to_add, ""); - type[il].Increase_Size (nr_to_add); - isdone[il].Increase_Size (nr_to_add, false); - } - - catch (bad_alloc) { - cout << "dim[il] " << dim[il] << "\tnr_to_add " << nr_to_add << endl; - JSCerror("Memory allocation failed in Scan_Thread_Set::Increase_Size."); - } - - return(true); - } - - void Scan_Thread_Set::Include_Thread (DP abs_data_value_ref, string label_ref, int type_ref) - { - // Determine which ilist index is to be used: - int il = int(-log(abs_data_value_ref)/logscale); - if (il < 0) il = 0; - if (il >= nlists) il = nlists - 1; - - if (nthreads_tot[il] > dim[il] - 10) { // Resize the Scan_Thread_Set list, by doubling its size - //cout << "Resizing threads list" << endl; - (*this).Increase_Size (il, dim[il]); - //cout << "Done resizing" << endl; - } - - //cout << "Including thread for label " << label_ref << " and type " << type_ref << endl; - - label[il][nthreads_tot[il] ] = label_ref; - type[il][nthreads_tot[il] ] = type_ref; - isdone[il][nthreads_tot[il] ] = false; - - nthreads_tot[il]++; - - } - - void Scan_Thread_Set::Merge (const Scan_Thread_Set& refset) - { - if (nlists != refset.nlists) JSCerror("nlists don't match in Scan_Thread_Set::Merge"); - - for (int il = 0; il < nlists; ++il) { - - if (nthreads_tot[il] + refset.nthreads_tot[il] >= dim[il]) - (*this).Increase_Size (il, nthreads_tot[il]/10 + refset.nthreads_tot[il]); - - for (int i = 0; i < refset.nthreads_tot[il]; ++i) { - label[il][nthreads_tot[il] ] = refset.label[il][i]; - type[il][nthreads_tot[il] ] = refset.type[il][i]; - isdone[il][nthreads_tot[il] ] = refset.isdone[il][i]; - nthreads_tot[il]++; - } - } - } - - void Scan_Thread_Set::Remove_Done_Threads (int il) - { - // If isdone[ithread] == true, remove from list - - int nr_removed = 0; - - for (int i = 0; i < nthreads_tot[il]; ++i) { - if (!isdone[il][i]) { - label[il][i - nr_removed] = label[il][i]; - type[il][i - nr_removed] = type[il][i]; - isdone[il][i - nr_removed] = isdone[il][i]; - } - else nr_removed++; - } - // Zero other entries: - for (int i = nthreads_tot[il] - nr_removed; i < nthreads_tot[il]; ++i) { - label[il][i] = ""; - type[il][i] = 0; - isdone[il][i] = false; - } - - nthreads_tot[il] -= nr_removed; - if (nthreads_done[il] != nr_removed) { - cout << nthreads_done[il] << "\t" << nr_removed << endl; - JSCerror("Miscount in removing threads during Scan_Thread_Set::Remove_Done_Threads."); - } - nthreads_done[il] = 0; - - return; - } - - void Scan_Thread_Set::Remove_Done_Threads () - { - // If isdone[ithread] == true, remove from list - - for (int il = 0; il < nlists; ++il) - (*this).Remove_Done_Threads(il); - - return; - } - - void Scan_Thread_Set::Clear () - { - for (int il = 0; il < nlists; ++il) { - - nthreads_tot[il] = 0; - nthreads_done[il] = 0; - label[il] = Vect (dim[il]); - type[il] = Vect (dim[il]); - isdone[il] = Vect (false, dim[il]); - } - } - - void Scan_Thread_Set::Save(const char* outfile_Cstr) - { - // We save only the undone threads, so after cleanup. - (*this).Remove_Done_Threads(); - - ofstream outfile; - - outfile.open(outfile_Cstr); - if (outfile.fail()) JSCerror("Could not open outfile... "); - outfile.precision(3); - - //cout << "Saving threads: nthreads_tot vector is" << endl; - //for (int il = 0; il < nlists; ++il) - //if (nthreads_tot[il] > 0) cout << il << "\t" << nthreads_tot[il] << "\t"; - //cout << endl; - - bool started = false; - for (int il = 0; il < nlists; ++il) { - - if (nthreads_tot[il] > 0) { - if (started) outfile << endl; - outfile << il << "\t" << label[il][0] << "\t" << type[il][0]; - started = true; - for (int i = 1; i < nthreads_tot[il]; ++i) { - outfile << endl << il << "\t" << label[il][i] << "\t" << type[il][i]; - } - } - } - - outfile.close(); - - return; - } - - void Scan_Thread_Set::Load (const char* thrfile_Cstr) - { - ifstream infile; - infile.open(thrfile_Cstr); - if(infile.fail()) { - cout << "Could not open file " << thrfile_Cstr << " in Scan_Thread_Set::Load." << endl; - JSCerror("Terminating."); - } - - (*this).Clear(); - - // First count the number of elements in each list: - Vect nthreads_read(0, nlists); - int il_read; - string label_read; - int type_read; - - while (infile.peek() != EOF) { - infile >> il_read; - infile >> label_read; - infile >> type_read; - nthreads_read[il_read]++; - } - - infile.close(); - - //cout << "nthreads_read vector: " << endl; - //for (int il = 0; il < nlists; ++il) - //if (nthreads_read[il] > 0) cout << il << "\t" << nthreads_read[il] << "\t"; - //cout << endl; - - // Now allocate the proper sizes: - for (int il = 0; il < nlists; ++il) { - (*this).Increase_Size (il, JSC::max(0, nthreads_read[il] - dim[il] + 10)); - nthreads_read[il] = 0; - } - - // Read threads data in: - ifstream infile2; - infile2.open(thrfile_Cstr); - infile.seekg(0); - while (infile2.peek() != EOF) { - infile2 >> il_read; - infile2 >> label[il_read][nthreads_read[il_read] ]; - infile2 >> type[il_read][nthreads_read[il_read] ]; - //isdone[il_read][nthreads_read[il_read] ] = false; // no need - nthreads_tot[il_read]++; - nthreads_read[il_read]++; - } - - infile2.close(); - - //cout << "Loading threads: nthreads_tot vector is" << endl; - //for (int il = 0; il < nlists; ++il) - //if (nthreads_tot[il] > 0) cout << il << "\t" << nthreads_tot[il] << "\t"; - //cout << endl; - - return; - } - -} // namespace JSC