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 
61  bool (*selfConsistencyCallback)(
62  BlockDiagonalizer &blockDiagonalizer
63  )
64  );
65 
74  void setMaxIterations(int maxIterations);
75 
79  void run();
80 
90  const double getEigenValue(int state);
91 
100  const double getEigenValue(const Index &blockIndex, int state);
101 
109  const std::complex<double> getAmplitude(int state, const Index &index);
110 
119  const std::complex<double> getAmplitude(
120  const Index &blockIndex,
121  int state,
122  const Index &intraBlockIndex
123  );
124 
131  unsigned int getFirstStateInBlock(const Index &index) const;
132 
139  unsigned int getLastStateInBlock(const Index &index) const;
140 
145  void setParallelExecution(bool parallelExecution);
146 private:
148  std::complex<double> *hamiltonian;
149 
151  double *eigenValues;
152 
154  std::complex<double> *eigenVectors;
155 
157  std::vector<unsigned int> numStatesPerBlock;
158 
160  std::vector<unsigned int> stateToBlockMap;
161 
163  std::vector<unsigned int> blockToStateMap;
164 
166  std::vector<unsigned int> blockSizes;
167 
169  std::vector<unsigned int> blockOffsets;
170 
172  std::vector<unsigned int> eigenVectorSizes;
173 
175  std::vector<unsigned int> eigenVectorOffsets;
176 
178  int numBlocks;
179 
181  int maxIterations;
182 
184  bool parallelExecution;
185 
188  bool (*selfConsistencyCallback)(BlockDiagonalizer &blockDiagonalizer);
189 
191  void init();
192 
194  void update();
195 
197  void solve();
198 };
199 
201  bool (*selfConsistencyCallback)(
202  BlockDiagonalizer &blockDiagonalizer
203  )
204 ){
205  this->selfConsistencyCallback = selfConsistencyCallback;
206 }
207 
208 inline void BlockDiagonalizer::setMaxIterations(int maxIterations){
209  this->maxIterations = maxIterations;
210 }
211 
212 inline const std::complex<double> BlockDiagonalizer::getAmplitude(
213  int state,
214  const Index &index
215 ){
216  const Model &model = getModel();
217  unsigned int block = stateToBlockMap.at(state);
218  unsigned int offset = eigenVectorOffsets.at(block);
219  unsigned int linearIndex = model.getBasisIndex(index);
220  unsigned int firstStateInBlock = blockToStateMap.at(block);
221  unsigned int lastStateInBlock = firstStateInBlock + numStatesPerBlock.at(block)-1;
222  offset += (state - firstStateInBlock)*numStatesPerBlock.at(block);
223  if(linearIndex >= firstStateInBlock && linearIndex <= lastStateInBlock)
224  return eigenVectors[offset + (linearIndex - firstStateInBlock)];
225  else
226  return 0;
227 }
228 
229 inline const std::complex<double> BlockDiagonalizer::getAmplitude(
230  const Index &blockIndex,
231  int state,
232  const Index &intraBlockIndex
233 ){
234  int firstStateInBlock = getModel().getHoppingAmplitudeSet().getFirstIndexInBlock(
235  blockIndex
236  );
237  unsigned int block = stateToBlockMap.at(firstStateInBlock);
238  TBTKAssert(
239  state >= 0 && state < (int)numStatesPerBlock.at(block),
240  "BlockDiagonalizer::getAmplitude()",
241  "Out of bound error. The block with block Index "
242  << blockIndex.toString() << " has "
243  << numStatesPerBlock.at(block) << " states, but state "
244  << state << " was requested.",
245  ""
246  );
247  unsigned int offset = eigenVectorOffsets.at(block) + state*numStatesPerBlock.at(block);
248  unsigned int linearIndex = getModel().getBasisIndex(
249  Index(blockIndex, intraBlockIndex)
250  );
251 // Streams::out << linearIndex << "\t" << Index(blockIndex, intraBlockIndex).toString() << "\n";
252 
253  return eigenVectors[offset + (linearIndex - firstStateInBlock)];
254 }
255 
256 inline const double BlockDiagonalizer::getEigenValue(int state){
257  return eigenValues[state];
258 }
259 
261  const Index &blockIndex,
262  int state
263 ){
265  blockIndex
266  );
267 
268  return eigenValues[offset + state];
269 }
270 
272  const Index &index
273 ) const{
274  unsigned int linearIndex = getModel().getBasisIndex(index);
275  unsigned int block = stateToBlockMap.at(linearIndex);
276 
277  return blockToStateMap.at(block);
278 }
279 
281  const Index &index
282 ) const{
283  unsigned int linearIndex = getModel().getBasisIndex(index);
284  unsigned int block = stateToBlockMap.at(linearIndex);
285 
286  return getFirstStateInBlock(index) + numStatesPerBlock.at(block)-1;
287 }
288 
290  bool parallelExecution
291 ){
292  this->parallelExecution = parallelExecution;
293 }
294 
295 }; //End of namespace Solver
296 }; //End of namespace TBTK
297 
298 #endif
Base class for Solvers.
Definition: Solver.h:32
Model & getModel()
Definition: Solver.h:63
void setParallelExecution(bool parallelExecution)
Definition: BlockDiagonalizer.h:289
const std::complex< double > getAmplitude(int state, const Index &index)
Definition: BlockDiagonalizer.h:212
const double getEigenValue(int state)
Definition: BlockDiagonalizer.h:256
Container of Model related information.
Timer class providing stacked tick and tock functions.
void setMaxIterations(int maxIterations)
Definition: BlockDiagonalizer.h:208
void setSelfConsistencyCallback(bool(*selfConsistencyCallback)(BlockDiagonalizer &blockDiagonalizer))
Definition: BlockDiagonalizer.h:200
const HoppingAmplitudeSet & getHoppingAmplitudeSet() const
Definition: Model.h:344
int getBasisIndex(const Index &index) const
Definition: Model.h:296
unsigned int getLastStateInBlock(const Index &index) const
Definition: BlockDiagonalizer.h:280
Base class for Solvers.
Flexible physical index.
Definition: Index.h:69
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:271
Definition: Communicator.h:28
Container of Model related information.
Definition: Model.h:52
std::string toString() const
Definition: Index.h:282
int getFirstIndexInBlock(const Index &blockIndex) const
Definition: HoppingAmplitudeSet.h:457