TBTK
BlockDiagonalizer.h
1 /* Copyright 2017 Kristofer Bj√∂rnson
2  *
3  * Licensed under the Apache License, Version 2.0 (the "License");
4  * you may not use this file except in compliance with the License.
5  * You may obtain a copy of the License at
6  *
7  * http://www.apache.org/licenses/LICENSE-2.0
8  *
9  * Unless required by applicable law or agreed to in writing, software
10  * distributed under the License is distributed on an "AS IS" BASIS,
11  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12  * See the License for the specific language governing permissions and
13  * limitations under the License.
14  */
15 
23 #ifndef COM_DAFER45_TBTK_SOLVER_BLOCK_DIAGONALIZER
24 #define COM_DAFER45_TBTK_SOLVER_BLOCK_DIAGONALIZER
25 
26 #include "TBTK/Communicator.h"
27 #include "TBTK/Model.h"
28 #include "TBTK/Solver/Solver.h"
29 #include "TBTK/Timer.h"
30 
31 #include <complex>
32 
33 namespace TBTK{
34 namespace Solver{
35 
43 class BlockDiagonalizer : public Solver, public Communicator{
44 public:
47 
49  virtual ~BlockDiagonalizer();
50 
54  bool (*selfConsistencyCallback)(
55  BlockDiagonalizer &blockDiagonalizer
56  )
57  );
58 
60  void setMaxIterations(int maxIterations);
61 
65  void run();
66 
68  const double getEigenValue(int state);
69 
73  const double getEigenValue(const Index &blockIndex, int state);
74 
80  const std::complex<double> getAmplitude(int state, const Index &index);
81 
87  const std::complex<double> getAmplitude(
88  const Index &blockIndex,
89  int state,
90  const Index &intraBlockIndex
91  );
92 
94  unsigned int getFirstStateInBlock(const Index &index) const;
95 
97  unsigned int getLastStateInBlock(const Index &index) const;
98 
100  void setParallelExecution(bool parallelExecution);
101 private:
103  std::complex<double> *hamiltonian;
104 
106  double *eigenValues;
107 
109  std::complex<double> *eigenVectors;
110 
112  std::vector<unsigned int> numStatesPerBlock;
113 
115  std::vector<unsigned int> stateToBlockMap;
116 
118  std::vector<unsigned int> blockToStateMap;
119 
121  std::vector<unsigned int> blockSizes;
122 
124  std::vector<unsigned int> blockOffsets;
125 
127  std::vector<unsigned int> eigenVectorSizes;
128 
130  std::vector<unsigned int> eigenVectorOffsets;
131 
133  int numBlocks;
134 
136  int maxIterations;
137 
139  bool parallelExecution;
140 
143  bool (*selfConsistencyCallback)(BlockDiagonalizer &blockDiagonalizer);
144 
146  void init();
147 
149  void update();
150 
152  void solve();
153 };
154 
156  bool (*selfConsistencyCallback)(
157  BlockDiagonalizer &blockDiagonalizer
158  )
159 ){
160  this->selfConsistencyCallback = selfConsistencyCallback;
161 }
162 
163 inline void BlockDiagonalizer::setMaxIterations(int maxIterations){
164  this->maxIterations = maxIterations;
165 }
166 
167 inline const std::complex<double> BlockDiagonalizer::getAmplitude(
168  int state,
169  const Index &index
170 ){
171  const Model &model = getModel();
172  unsigned int block = stateToBlockMap.at(state);
173  unsigned int offset = eigenVectorOffsets.at(block);
174  unsigned int linearIndex = model.getBasisIndex(index);
175  unsigned int firstStateInBlock = blockToStateMap.at(block);
176  unsigned int lastStateInBlock = firstStateInBlock + numStatesPerBlock.at(block)-1;
177  offset += (state - firstStateInBlock)*numStatesPerBlock.at(block);
178  if(linearIndex >= firstStateInBlock && linearIndex <= lastStateInBlock)
179  return eigenVectors[offset + (linearIndex - firstStateInBlock)];
180  else
181  return 0;
182 }
183 
184 inline const std::complex<double> BlockDiagonalizer::getAmplitude(
185  const Index &blockIndex,
186  int state,
187  const Index &intraBlockIndex
188 ){
189  int firstStateInBlock = getModel().getHoppingAmplitudeSet()->getFirstIndexInBlock(
190  blockIndex
191  );
192  unsigned int block = stateToBlockMap.at(firstStateInBlock);
193  TBTKAssert(
194  state >= 0 && state < (int)numStatesPerBlock.at(block),
195  "BlockDiagonalizer::getAmplitude()",
196  "Out of bound error. The block with block Index "
197  << blockIndex.toString() << " has "
198  << numStatesPerBlock.at(block) << " states, but state "
199  << state << " was requested.",
200  ""
201  );
202  unsigned int offset = eigenVectorOffsets.at(block) + state*numStatesPerBlock.at(block);
203  unsigned int linearIndex = getModel().getBasisIndex(
204  Index(blockIndex, intraBlockIndex)
205  );
206 // Streams::out << linearIndex << "\t" << Index(blockIndex, intraBlockIndex).toString() << "\n";
207 
208  return eigenVectors[offset + (linearIndex - firstStateInBlock)];
209 }
210 
211 inline const double BlockDiagonalizer::getEigenValue(int state){
212  return eigenValues[state];
213 }
214 
216  const Index &blockIndex,
217  int state
218 ){
220  blockIndex
221  );
222 
223  return eigenValues[offset + state];
224 }
225 
227  const Index &index
228 ) const{
229  unsigned int linearIndex = getModel().getBasisIndex(index);
230  unsigned int block = stateToBlockMap.at(linearIndex);
231 
232  return blockToStateMap.at(block);
233 }
234 
236  const Index &index
237 ) const{
238  unsigned int linearIndex = getModel().getBasisIndex(index);
239  unsigned int block = stateToBlockMap.at(linearIndex);
240 
241  return getFirstStateInBlock(index) + numStatesPerBlock.at(block)-1;
242 }
243 
245  bool parallelExecution
246 ){
247  this->parallelExecution = parallelExecution;
248 }
249 
250 }; //End of namespace Solver
251 }; //End of namespace TBTK
252 
253 #endif
Base class for Solvers.
Definition: Solver.h:32
Model & getModel()
Definition: Solver.h:57
void setParallelExecution(bool parallelExecution)
Definition: BlockDiagonalizer.h:244
const std::complex< double > getAmplitude(int state, const Index &index)
Definition: BlockDiagonalizer.h:167
const double getEigenValue(int state)
Definition: BlockDiagonalizer.h:211
Container of Model related information.
Timer class providing stacked tick and tock functions.
void setMaxIterations(int maxIterations)
Definition: BlockDiagonalizer.h:163
void setSelfConsistencyCallback(bool(*selfConsistencyCallback)(BlockDiagonalizer &blockDiagonalizer))
Definition: BlockDiagonalizer.h:155
int getBasisIndex(const Index &index) const
Definition: Model.h:214
unsigned int getLastStateInBlock(const Index &index) const
Definition: BlockDiagonalizer.h:235
Base class for Solvers.
Flexible physical index.
Definition: Index.h:70
Definition: ModelFactory.h:35
Solves a block diagonal Model using diagonalization.
Definition: BlockDiagonalizer.h:43
Base class that communicate their current status during execution.
unsigned int getFirstStateInBlock(const Index &index) const
Definition: BlockDiagonalizer.h:226
const HoppingAmplitudeSet * getHoppingAmplitudeSet() const
Definition: Model.h:262
Definition: Communicator.h:28
Container of Model related information.
Definition: Model.h:52
std::string toString() const
Definition: Index.h:265
int getFirstIndexInBlock(const Index &blockIndex) const
Definition: HoppingAmplitudeSet.h:344