TBTK
Need a break? Support the development by playing Polarity Puzzles
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/BlockStructureDescriptor.h"
27 #include "TBTK/CArray.h"
28 #include "TBTK/Communicator.h"
29 #include "TBTK/Model.h"
30 #include "TBTK/Solver/Solver.h"
31 #include "TBTK/Timer.h"
32 
33 #include <complex>
34 
35 namespace TBTK{
36 namespace Solver{
37 
68 class BlockDiagonalizer : public Solver, public Communicator{
69 public:
72  public:
77  virtual bool selfConsistencyCallback(
78  BlockDiagonalizer &blockDiagonalizer
79  ) = 0;
80  };
81 
84 
96  );
97 
106  void setMaxIterations(int maxIterations);
107 
111  void run();
112 
122  const double getEigenValue(int state);
123 
132  const double getEigenValue(const Index &blockIndex, int state);
133 
141  const std::complex<double> getAmplitude(int state, const Index &index);
142 
151  const std::complex<double> getAmplitude(
152  const Index &blockIndex,
153  int state,
154  const Index &intraBlockIndex
155  );
156 
163  unsigned int getFirstStateInBlock(const Index &index) const;
164 
171  unsigned int getLastStateInBlock(const Index &index) const;
172 
177  void setParallelExecution(bool parallelExecution);
178 private:
180  CArray<std::complex<double>> hamiltonian;
181 
183  CArray<double> eigenValues;
184 
186  CArray<std::complex<double>> eigenVectors;
187 
189  BlockStructureDescriptor blockStructureDescriptor;
190 
192  std::vector<unsigned int> blockSizes;
193 
195  std::vector<unsigned int> blockOffsets;
196 
198  std::vector<unsigned int> eigenVectorSizes;
199 
201  std::vector<unsigned int> eigenVectorOffsets;
202 
204  int maxIterations;
205 
207  bool parallelExecution;
208 
212 
214  void init();
215 
217  void update();
218 
220  void solve();
221 };
222 
225 ){
226  this->selfConsistencyCallback = &selfConsistencyCallback;
227 }
228 
229 inline void BlockDiagonalizer::setMaxIterations(int maxIterations){
230  this->maxIterations = maxIterations;
231 }
232 
233 inline const std::complex<double> BlockDiagonalizer::getAmplitude(
234  int state,
235  const Index &index
236 ){
237  const Model &model = getModel();
238  unsigned int block = blockStructureDescriptor.getBlockIndex(state);
239  unsigned int offset = eigenVectorOffsets.at(block);
240  unsigned int linearIndex = model.getBasisIndex(index);
241  unsigned int firstStateInBlock
242  = blockStructureDescriptor.getFirstStateInBlock(block);
243  unsigned int lastStateInBlock = firstStateInBlock
244  + blockStructureDescriptor.getNumStatesInBlock(block)-1;
245  offset += (
246  state - firstStateInBlock
247  )*blockStructureDescriptor.getNumStatesInBlock(block);
248  if(
249  linearIndex >= firstStateInBlock
250  && linearIndex <= lastStateInBlock
251  ){
252  return eigenVectors[
253  offset + (linearIndex - firstStateInBlock)
254  ];
255  }
256  else{
257  return 0;
258  }
259 }
260 
261 inline const std::complex<double> BlockDiagonalizer::getAmplitude(
262  const Index &blockIndex,
263  int state,
264  const Index &intraBlockIndex
265 ){
266  int firstStateInBlock = getModel().getHoppingAmplitudeSet(
267  ).getFirstIndexInBlock(blockIndex);
268  unsigned int block = blockStructureDescriptor.getBlockIndex(
269  firstStateInBlock
270  );
271  TBTKAssert(
272  state >= 0
273  && state < (int)blockStructureDescriptor.getNumStatesInBlock(
274  block
275  ),
276  "BlockDiagonalizer::getAmplitude()",
277  "Out of bound error. The block with block Index "
278  << blockIndex.toString() << " has "
279  << blockStructureDescriptor.getNumStatesInBlock(block)
280  << " states, but state " << state << " was requested.",
281  ""
282  );
283  unsigned int offset = eigenVectorOffsets.at(block)
284  + state*blockStructureDescriptor.getNumStatesInBlock(block);
285  unsigned int linearIndex = getModel().getBasisIndex(
286  Index(blockIndex, intraBlockIndex)
287  );
288 
289  return eigenVectors[offset + (linearIndex - firstStateInBlock)];
290 }
291 
292 inline const double BlockDiagonalizer::getEigenValue(int state){
293  return eigenValues[state];
294 }
295 
297  const Index &blockIndex,
298  int state
299 ){
301  blockIndex
302  );
303 
304  return eigenValues[offset + state];
305 }
306 
308  const Index &index
309 ) const{
310  unsigned int linearIndex = getModel().getBasisIndex(index);
311  unsigned int block
312  = blockStructureDescriptor.getBlockIndex(linearIndex);
313 
314  return blockStructureDescriptor.getFirstStateInBlock(block);
315 }
316 
318  const Index &index
319 ) const{
320  unsigned int linearIndex = getModel().getBasisIndex(index);
321  unsigned int block
322  = blockStructureDescriptor.getBlockIndex(linearIndex);
323 
324  return getFirstStateInBlock(index)
325  + blockStructureDescriptor.getNumStatesInBlock(block) - 1;
326 }
327 
329  bool parallelExecution
330 ){
331  this->parallelExecution = parallelExecution;
332 }
333 
334 }; //End of namespace Solver
335 }; //End of namespace TBTK
336 
337 #endif
Base class for Solvers.
Definition: Solver.h:42
Model & getModel()
Definition: Solver.h:73
void setParallelExecution(bool parallelExecution)
Definition: BlockDiagonalizer.h:328
const std::complex< double > getAmplitude(int state, const Index &index)
Definition: BlockDiagonalizer.h:233
int getBasisIndex(const Index &index) const
Definition: Model.h:331
const double getEigenValue(int state)
Definition: BlockDiagonalizer.h:292
Container of Model related information.
int getFirstIndexInBlock(const Index &blockIndex) const
Definition: HoppingAmplitudeSet.h:386
std::string toString() const
Definition: Index.h:349
Container for a C style array.
virtual bool selfConsistencyCallback(BlockDiagonalizer &blockDiagonalizer)=0
A Timer for measuring execution time.
void setMaxIterations(int maxIterations)
Definition: BlockDiagonalizer.h:229
void setSelfConsistencyCallback(SelfConsistencyCallback &selfConsistencyCallback)
Definition: BlockDiagonalizer.h:223
Base class for Solvers.
unsigned int getFirstStateInBlock(const Index &index) const
Definition: BlockDiagonalizer.h:307
Physical index.
Definition: Index.h:44
Container for a C style array.
Definition: CArray.h:44
Definition: Boolean.h:32
Solves a block diagonal Model using diagonalization.
Definition: BlockDiagonalizer.h:68
Base class for classes that can communicate their status during execution.
unsigned int getLastStateInBlock(const Index &index) const
Definition: BlockDiagonalizer.h:317
Base class for classes that can communicate their status during execution.
Definition: Communicator.h:56
Container of Model related information.
Definition: Model.h:57
const HoppingAmplitudeSet & getHoppingAmplitudeSet() const
Definition: Model.h:367