Quandary
Loading...
Searching...
No Matches
mastereq.hpp
Go to the documentation of this file.
1#include "defs.hpp"
2#include "oscillator.hpp"
3#include "util.hpp"
4#include <optional>
5#include <petscts.h>
6#include <vector>
7#include <assert.h>
8#include <iostream>
9#include "gate.hpp"
10#include "hamiltonianfilereader.hpp"
11
12#pragma once
13
14
21typedef struct {
22 PetscInt dim;
23 std::vector<size_t> nlevels;
24 IS *isu, *isv;
25 Oscillator** oscil_vec;
26 std::vector<double> crosskerr;
27 std::vector<double> Jkl;
28 std::vector<double> eta;
29 DecoherenceType decoherence_type;
30 bool addT1, addT2;
31 std::vector<double> control_Re;
32 std::vector<double> control_Im;
33 std::vector<Mat> Ac_vec;
34 std::vector<Mat> Bc_vec;
35 Mat *Ad;
36 Mat *Bd;
37 std::vector<Mat> Ad_vec;
38 std::vector<Mat> Bd_vec;
39 std::vector<double> Bd_coeffs; //< Time-dependent coefficients for dipole-dipole coupling matrices: cos(eta_k*t)
40 std::vector<double> Ad_coeffs; //< Time-dependent coefficients for dipole-dipole coupling matrices: sin(eta_k*t)
41 Vec *aux;
42 double time;
43} MatShellCtx;
44
45
54int applyRHS_matfree_1Osc(Mat RHS, Vec x, Vec y);
55int applyRHS_matfree_transpose_1Osc(Mat RHS, Vec x, Vec y);
56int applyRHS_matfree_2Osc(Mat RHS, Vec x, Vec y);
57int applyRHS_matfree_transpose_2Osc(Mat RHS, Vec x, Vec y);
58int applyRHS_matfree_3Osc(Mat RHS, Vec x, Vec y);
59int applyRHS_matfree_transpose_3Osc(Mat RHS, Vec x, Vec y);
60int applyRHS_matfree_4Osc(Mat RHS, Vec x, Vec y);
61int applyRHS_matfree_transpose_4Osc(Mat RHS, Vec x, Vec y);
62int applyRHS_matfree_5Osc(Mat RHS, Vec x, Vec y);
63int applyRHS_matfree_transpose_5Osc(Mat RHS, Vec x, Vec y);
64int applyRHS_sparsemat(Mat RHS, Vec x, Vec y);
65int applyRHS_sparsemat_transpose(Mat RHS, Vec x, Vec y);
66
67
92class MasterEq{
93
94 protected:
95 PetscInt dim;
96 PetscInt dim_rho;
97 PetscInt dim_ess;
98 int noscillators;
99 Oscillator** oscil_vec;
100
101 Mat RHS;
102 MatShellCtx RHSctx;
103
104 std::vector<Mat> Ac_vec; // Vector of constant mats for time-varying control term (real). One for each oscillators.
105 std::vector<Mat> Bc_vec; // Vector of constant mats for time-varying control term (imag). One for each oscillators.
106 Mat Ad, Bd; // Real and imaginary part of constant system matrix
107 std::vector<Mat> Ad_vec; // Vector of constant mats for Dipole-Dipole coupling term in drift Hamiltonian (real)
108 std::vector<Mat> Bd_vec; // Vector of constant mats for Dipole-Dipole coupling term in drift Hamiltonian (imag)
109
110 std::vector<double> crosskerr;
111 std::vector<double> Jkl;
112 std::vector<double> eta;
113 bool addT1, addT2;
114
115 int mpirank_world;
116 int mpirank_petsc;
117 int mpisize_petsc;
118 PetscInt localsize_u;
119 PetscInt ilow;
120 PetscInt iupp;
121
122 IS isu, isv;
123 Vec aux;
124 bool quietmode;
125 std::optional<std::string> hamiltonian_file_Hsys;
126 std::optional<std::string> hamiltonian_file_Hc;
127
128 public:
129 std::vector<size_t> nlevels;
130 std::vector<size_t> nessential;
131 bool usematfree;
132 DecoherenceType decoherence_type;
133
134 public:
135 MasterEq();
136
144 MasterEq(const Config& config, Oscillator** oscil_vec_, bool quietmode=false);
145
146 ~MasterEq();
147
154 void set_RHS_MatMult_operation();
155
159 void initSparseMatSolver();
160
167 Oscillator* getOscillator(const size_t i) { return oscil_vec[i]; }
168
174 size_t getNOscillators() { return noscillators; }
175
181 PetscInt getDim(){ return dim; }
182
188 PetscInt getDimEss(){ return dim_ess; }
189
195 PetscInt getDimRho(){ return dim_rho; }
196
206 int assemble_RHS(const double t);
207
213 Mat getRHS();
214
227 void compute_dRHS_dParams(const double t,const Vec x,const Vec x_bar, const double alpha, Vec grad);
228
234 void setControlAmplitudes(const Vec x);
235
242 double expectedEnergy(const Vec x);
243
250 void population(const Vec x, std::vector<double> &population_com);
251
252 // /* Compute reduced density operator for a sub-system defined by IDs in the oscilIDs vector */
253 // void createReducedDensity(const Vec rho, Vec *reduced, const std::vector<int>& oscilIDs);
254 // /* Derivative of reduced density computation */
255 // void createReducedDensity_diff(Vec rhobar, const Vec reducedbar, const std::vector<int>& oscilIDs);
256};
257
258
278void compute_dRHS_dParams_sparsemat(const double t,const Vec x,const Vec x_bar, const double alpha, Vec grad, std::vector<size_t>& nlevels, IS isu, IS isv, std::vector<Mat>& Ac_vec, std::vector<Mat>& Bc_vec, Vec aux, Oscillator** oscil_vec);
279
296void compute_dRHS_dParams_matfree(const PetscInt dim, const double t,const Vec x,const Vec x_bar, const double alpha, Vec grad, std::vector<size_t>& nlevels, DecoherenceType decoherence_type, Oscillator** oscil_vec);
297
298
299
300// Inline functions for the Matrix-free RHS application
301
309inline double H_detune(const double detuning0, const int a) {
310 return detuning0*a;
311};
312
320inline double H_selfkerr(const double xi0, const int a) {
321 return - xi0 / 2.0 * a * (a-1);
322};
323
332inline double L2(const double dephase0, const int i0, const int i0p){
333 return dephase0 * ( i0*i0p - 1./2. * (i0*i0 + i0p*i0p) );
334};
335
344inline double L1diag(const double decay0, const int i0, const int i0p){
345 return - decay0 / 2.0 * ( i0 + i0p ) ;
346};
347
356inline int TensorGetIndex(const int nlevels0, const int i0, const int i0p){
357 return i0 + nlevels0 * i0p;
358};
359
369inline double H_detune(const double detuning0, const double detuning1, const int a, const int b) {
370 return detuning0*a + detuning1*b;
371};
372
382inline double H_selfkerr(const double xi0, const double xi1, const int a, const int b) {
383 return - xi0 / 2.0 * a * (a-1) - xi1 / 2.0 * b * (b-1);
384};
385
394inline double H_crosskerr(const double xi01, const int a, const int b) {
395 return - xi01 * a * b;
396};
397
409inline double L2(const double dephase0, const double dephase1, const int i0, const int i1, const int i0p, const int i1p){
410 return dephase0 * ( i0*i0p - 1./2. * (i0*i0 + i0p*i0p) ) + dephase1 * ( i1*i1p - 1./2. * (i1*i1 + i1p*i1p) );
411};
412
424inline double L1diag(const double decay0, const double decay1, const int i0, const int i1, const int i0p, const int i1p){
425 return - decay0 / 2.0 * ( i0 + i0p ) - decay1 / 2.0 * ( i1 + i1p );
426};
427
439inline int TensorGetIndex(const int nlevels0, const int nlevels1,const int i0, const int i1, const int i0p, const int i1p){
440 return i0*nlevels1 + i1 + (nlevels0 * nlevels1) * ( i0p * nlevels1 + i1p);
441};
442
443
444// Inlines for Matrix-free RHS for 3 oscillators
445// See documentation for 1-2 oscillator functions above for parameter descriptions
446inline double H_detune(const double detuning0, const double detuning1, const double detuning2, const int i0, const int i1, const int i2) {
447 return detuning0*i0 + detuning1*i1 + detuning2*i2;
448};
449inline double H_selfkerr(const double xi0, const double xi1, const double xi2, const int i0, const int i1, const int i2) {
450 return - xi0 / 2.0 * i0 * (i0-1) - xi1 / 2.0 * i1 * (i1-1) - xi2 / 2.0 * i2 * (i2-1);
451};
452inline double H_crosskerr(const double xi01, const double xi02, const double xi12, const int i0, const int i1, const int i2) {
453 return - xi01 * i0 * i1 - xi02 * i0 * i2 - xi12 * i1 * i2;
454};
455inline double L2(const double dephase0, const double dephase1, const double dephase2, const int i0, const int i1, const int i2, const int i0p, const int i1p, const int i2p){
456 return dephase0 * ( i0*i0p - 1./2. * (i0*i0 + i0p*i0p) )
457 + dephase1 * ( i1*i1p - 1./2. * (i1*i1 + i1p*i1p) )
458 + dephase2 * ( i2*i2p - 1./2. * (i2*i2 + i2p*i2p) );
459};
460inline double L1diag(const double decay0, const double decay1, const double decay2, const int i0, const int i1, const int i2, const int i0p, const int i1p, const int i2p){
461 return - decay0 / 2.0 * ( i0 + i0p )
462 - decay1 / 2.0 * ( i1 + i1p )
463 - decay2 / 2.0 * ( i2 + i2p );
464};
465inline int TensorGetIndex(const int nlevels0, const int nlevels1, const int nlevels2, const int i0, const int i1, const int i2, const int i0p, const int i1p, const int i2p){
466 return i0*nlevels1*nlevels2 + i1*nlevels2 + i2 + (nlevels0 * nlevels1 * nlevels2) * ( i0p * nlevels1*nlevels2 + i1p*nlevels2 + i2p);
467};
468
469
470// Matrix-free solver inlines for 4 oscillators
471// See documentation for 1-2 oscillator functions above for parameter descriptions
472inline double H_detune(const double detuning0, const double detuning1, const double detuning2, const double detuning3, const int i0, const int i1, const int i2, const int i3) {
473 return detuning0*i0 + detuning1*i1 + detuning2*i2 + detuning3*i3;
474};
475inline double H_selfkerr(const double xi0, const double xi1, const double xi2, const double xi3, const int i0, const int i1, const int i2, const int i3) {
476 return - xi0 / 2.0 * i0 * (i0-1) - xi1 / 2.0 * i1 * (i1-1) - xi2 / 2.0 * i2 * (i2-1) - xi3/2.0 * i3 * (i3-1);
477};
478inline double H_crosskerr(const double xi01, const double xi02, const double xi03, const double xi12, const double xi13, const double xi23, const int i0, const int i1, const int i2, const int i3) {
479 return - xi01 * i0 * i1 - xi02 * i0 * i2 - xi03*i0*i3 - xi12 * i1 * i2 - xi13*i1*i3 - xi23*i2*i3;
480};
481inline double L2(const double dephase0, const double dephase1, const double dephase2, const double dephase3, const int i0, const int i1, const int i2, const int i3, const int i0p, const int i1p, const int i2p, const int i3p){
482 return dephase0 * ( i0*i0p - 1./2. * (i0*i0 + i0p*i0p) )
483 + dephase1 * ( i1*i1p - 1./2. * (i1*i1 + i1p*i1p) )
484 + dephase2 * ( i2*i2p - 1./2. * (i2*i2 + i2p*i2p) )
485 + dephase3 * ( i3*i3p - 1./2. * (i3*i3 + i3p*i3p) );
486};
487inline double L1diag(const double decay0, const double decay1, const double decay2, const double decay3, const int i0, const int i1, const int i2, const int i3, const int i0p, const int i1p, const int i2p, const int i3p){
488 return - decay0 / 2.0 * ( i0 + i0p )
489 - decay1 / 2.0 * ( i1 + i1p )
490 - decay2 / 2.0 * ( i2 + i2p )
491 - decay3 / 2.0 * ( i3 + i3p );
492};
493inline int TensorGetIndex(const int nlevels0, const int nlevels1, const int nlevels2, const int nlevels3, const int i0, const int i1, const int i2, const int i3, const int i0p, const int i1p, const int i2p, const int i3p){
494 return i0*nlevels1*nlevels2*nlevels3 + i1*nlevels2*nlevels3 + i2*nlevels3 + i3 + (nlevels0 * nlevels1 * nlevels2 * nlevels3) * ( i0p * nlevels1*nlevels2*nlevels3 + i1p*nlevels2*nlevels3 + i2p*nlevels3 + i3p);
495}
496
497// Matrix-free solver inlines for 5 oscillators
498// See documentation for 1-2 oscillator functions above for parameter descriptions
499inline double H_detune(const double detuning0, const double detuning1, const double detuning2, const double detuning3, const double detuning4, const int i0, const int i1, const int i2, const int i3, const int i4) {
500 return detuning0*i0 + detuning1*i1 + detuning2*i2 + detuning3*i3 + detuning4*i4;
501};
502inline double H_selfkerr(const double xi0, const double xi1, const double xi2, const double xi3, const double xi4, const int i0, const int i1, const int i2, const int i3, const int i4) {
503 return - xi0 / 2.0 * i0 * (i0-1) - xi1 / 2.0 * i1 * (i1-1) - xi2 / 2.0 * i2 * (i2-1) - xi3/2.0 * i3 * (i3-1) - xi4/2.0 * i4 * (i4-1);
504};
505inline double H_crosskerr(const double xi01, const double xi02, const double xi03, const double xi04, const double xi12, const double xi13, const double xi14, const double xi23, const double xi24, const double xi34, const int i0, const int i1, const int i2, const int i3, const int i4) {
506 return - xi01 * i0 * i1 - xi02 * i0 * i2 - xi03*i0*i3 - xi04*i0*i4 - xi12 * i1 * i2 - xi13*i1*i3 - xi14*i1*i4 - xi23*i2*i3 - xi24*i2*i4 - xi34*i3*i4;
507};
508inline double L2(const double dephase0, const double dephase1, const double dephase2, const double dephase3, const double dephase4, const int i0, const int i1, const int i2, const int i3, const int i4, const int i0p, const int i1p, const int i2p, const int i3p, const int i4p){
509 return dephase0 * ( i0*i0p - 1./2. * (i0*i0 + i0p*i0p) )
510 + dephase1 * ( i1*i1p - 1./2. * (i1*i1 + i1p*i1p) )
511 + dephase2 * ( i2*i2p - 1./2. * (i2*i2 + i2p*i2p) )
512 + dephase3 * ( i3*i3p - 1./2. * (i3*i3 + i3p*i3p) )
513 + dephase4 * ( i4*i4p - 1./2. * (i4*i4 + i4p*i4p) );
514};
515inline double L1diag(const double decay0, const double decay1, const double decay2, const double decay3, const double decay4, const int i0, const int i1, const int i2, const int i3, const int i4, const int i0p, const int i1p, const int i2p, const int i3p, const int i4p){
516 return - decay0 / 2.0 * ( i0 + i0p )
517 - decay1 / 2.0 * ( i1 + i1p )
518 - decay2 / 2.0 * ( i2 + i2p )
519 - decay3 / 2.0 * ( i3 + i3p )
520 - decay4 / 2.0 * ( i4 + i4p );
521};
522inline int TensorGetIndex(const int nlevels0, const int nlevels1, const int nlevels2, const int nlevels3, const int nlevels4, const int i0, const int i1, const int i2, const int i3, const int i4, const int i0p, const int i1p, const int i2p, const int i3p, const int i4p){
523 return i0*nlevels1*nlevels2*nlevels3*nlevels4 + i1*nlevels2*nlevels3*nlevels4 + i2*nlevels3*nlevels4 + i3*nlevels4 + i4 + (nlevels0 * nlevels1 * nlevels2 * nlevels3*nlevels4) * ( i0p * nlevels1*nlevels2*nlevels3*nlevels4 + i1p*nlevels2*nlevels3*nlevels4 + i2p*nlevels3*nlevels4+ i3p*nlevels4 + i4p);
524}
525
526
546inline void dRHSdp_getcoeffs(const PetscInt dim, const int it, const int n, const int np, const int i, const int ip, const int stridei, const int strideip, const double* xptr, double* res_p_re, double* res_p_im, double* res_q_re, double* res_q_im) {
547
548 *res_p_re = 0.0;
549 *res_p_im = 0.0;
550 *res_q_re = 0.0;
551 *res_q_im = 0.0;
552
553 /* ik+1..,ik'.. term */
554 if (i < n-1) {
555 int itx = it + stridei;
556 double xre = xptr[itx];
557 double xim = xptr[itx+dim];
558 double sq = sqrt(i + 1);
559 *res_p_re += sq * xim;
560 *res_p_im += - sq * xre;
561 *res_q_re += sq * xre;
562 *res_q_im += sq * xim;
563 }
564 /* \rho(ik..,ik'+1..) */
565 if (ip < np-1) {
566 int itx = it + strideip;
567 double xre = xptr[itx];
568 double xim = xptr[itx+dim];
569 double sq = sqrt(ip + 1);
570 *res_p_re += - sq * xim;
571 *res_p_im += + sq * xre;
572 *res_q_re += sq * xre;
573 *res_q_im += sq * xim;
574 }
575 /* \rho(ik-1..,ik'..) */
576 if (i > 0) {
577 int itx = it - stridei;
578 double xre = xptr[itx];
579 double xim = xptr[itx+dim];
580 double sq = sqrt(i);
581 *res_p_re += + sq * xim;
582 *res_p_im += - sq * xre;
583 *res_q_re += - sq * xre;
584 *res_q_im += - sq * xim;
585 }
586 /* \rho(ik..,ik'-1..) */
587 if (ip > 0) {
588 int itx = it - strideip;
589 double xre = xptr[itx];
590 double xim = xptr[itx+dim];
591 double sq = sqrt(ip);
592 *res_p_re += - sq * xim;
593 *res_p_im += + sq * xre;
594 *res_q_re += - sq * xre;
595 *res_q_im += - sq * xim;
596 }
597}
598
625inline void Jkl_coupling(const PetscInt dim, const int it, const int ni, const int nj, const int nip, const int njp, const int i, const int ip, const int j, const int jp, const int stridei, const int strideip, const int stridej, const int stridejp, const double* xptr, const double Jij, const double cosij, const double sinij, double* yre, double* yim) {
626 if (fabs(Jij)>1e-10) {
627 // 1) J_kl (-icos + sin) * ρ_{E−k+l i, i′}
628 if (i > 0 && j < nj-1) {
629 int itx = it - stridei + stridej;
630 double xre = xptr[itx];
631 double xim = xptr[itx+dim];
632 double sq = sqrt(i * (j + 1));
633 // sin u + cos v + i ( -cos u + sin v)
634 *yre += Jij * sq * ( cosij * xim + sinij * xre);
635 *yim += Jij * sq * ( - cosij * xre + sinij * xim);
636 }
637 // 2) J_kl (−icos − sin)sqrt(il*(ik +1)) ρ_{E+k−li,i′}
638 if (i < ni-1 && j > 0) {
639 int itx = it + stridei - stridej; // E+k-l i, i'
640 double xre = xptr[itx];
641 double xim = xptr[itx+dim];
642 double sq = sqrt(j * (i + 1)); // sqrt( il*(ik+1))
643 // -sin u + cos v + i (-cos u - sin v)
644 *yre += Jij * sq * ( cosij * xim - sinij * xre);
645 *yim += Jij * sq * ( - cosij * xre - sinij * xim);
646 }
647 // 3) J_kl ( icos + sin)sqrt(ik'*(il' +1)) ρ_{i,E-k+li'}
648 if (ip > 0 && jp < njp-1) {
649 int itx = it - strideip + stridejp; // i, E-k+l i'
650 double xre = xptr[itx];
651 double xim = xptr[itx+dim];
652 double sq = sqrt(ip * (jp + 1)); // sqrt( ik'*(il'+1))
653 // sin u - cos v + i ( cos u + sin v)
654 *yre += Jij * sq * ( - cosij * xim + sinij * xre);
655 *yim += Jij * sq * ( cosij * xre + sinij * xim);
656 }
657 // 4) J_kl ( icos - sin)sqrt(il'*(ik' +1)) ρ_{i,E+k-li'}
658 if (ip < nip-1 && jp > 0) {
659 int itx = it + strideip - stridejp; // i, E+k-l i'
660 double xre = xptr[itx];
661 double xim = xptr[itx+dim];
662 double sq = sqrt(jp * (ip + 1)); // sqrt( il'*(ik'+1))
663 // - sin u - cos v + i ( cos u - sin v)
664 *yre += Jij * sq * ( - cosij * xim - sinij * xre);
665 *yim += Jij * sq * ( cosij * xre - sinij * xim);
666 }
667 }
668}
669
694inline void Jkl_coupling_T(const PetscInt dim, const int it, const int ni, const int nj, const int nip, const int njp, const int i, const int ip, const int j, const int jp, const int stridei, const int strideip, const int stridej, const int stridejp, const double* xptr, const double Jij, const double cosij, const double sinij, double* yre, double* yim) {
695 if (fabs(Jij)>1e-10) {
696 // 1) [...] * \bar y_{E+k-l i, i′}
697 if (i < ni-1 && j > 0) {
698 int itx = it + stridei - stridej;
699 double xre = xptr[itx];
700 double xim = xptr[itx+dim];
701 double sq = sqrt(j * (i + 1));
702 *yre += Jij * sq * ( - cosij * xim + sinij * xre);
703 *yim += Jij * sq * ( + cosij * xre + sinij * xim);
704 }
705 // 2) J_kl (−icos − sin)sqrt(ik*(il +1)) \bar y_{E-k+li,i′}
706 if (i > 0 && j < nj-1) {
707 int itx = it - stridei + stridej; // E-k+l i, i'
708 double xre = xptr[itx];
709 double xim = xptr[itx+dim];
710 double sq = sqrt(i * (j + 1)); // sqrt( ik*(il+1))
711 *yre += Jij * sq * ( - cosij * xim - sinij * xre);
712 *yim += Jij * sq * ( + cosij * xre - sinij * xim);
713 }
714 // 3) J_kl ( icos + sin)sqrt(il'*(ik' +1)) \bar y_{i,E+k-li'}
715 if (ip < nip-1 && jp > 0) {
716 int itx = it + strideip - stridejp; // i, E+k-l i'
717 double xre = xptr[itx];
718 double xim = xptr[itx+dim];
719 double sq = sqrt(jp * (ip + 1)); // sqrt( il'*(ik'+1))
720 *yre += Jij * sq * ( cosij * xim + sinij * xre);
721 *yim += Jij * sq * ( - cosij * xre + sinij * xim);
722 }
723 // 4) J_kl ( icos - sin)sqrt(ik'*(il' +1)) \bar y_{i,E-k+li'}
724 if (ip > 0 && jp < njp-1) {
725 int itx = it - strideip + stridejp; // i, E-k+l i'
726 double xre = xptr[itx];
727 double xim = xptr[itx+dim];
728 double sq = sqrt(ip * (jp + 1)); // sqrt( ik'*(il'+1))
729 *yre += Jij * sq * ( cosij * xim - sinij * xre);
730 *yim += Jij * sq * ( - cosij * xre - sinij * xim);
731 }
732 }
733
734}
735
751inline void L1decay(const PetscInt dim, const int it, const int n, const int i, const int ip, const int stridei, const int strideip, const double* xptr, const double decayi, double* yre, double* yim){
752 if (fabs(decayi) > 1e-12) {
753 if (i < n-1 && ip < n-1) {
754 double l1off = decayi * sqrt((i+1)*(ip+1));
755 int itx = it + stridei + strideip;
756 double xre = xptr[itx];
757 double xim = xptr[itx+dim];
758 *yre += l1off * xre;
759 *yim += l1off * xim;
760 }
761 }
762}
763
764
779inline void L1decay_T(const PetscInt dim, const int it, const int i, const int ip, const int stridei, const int strideip, const double* xptr, const double decayi, double* yre, double* yim){
780 if (fabs(decayi) > 1e-12) {
781 if (i > 0 && ip > 0) {
782 double l1off = decayi * sqrt(i*ip);
783 int itx = it - stridei - strideip;
784 double xre = xptr[itx];
785 double xim = xptr[itx+dim];
786 *yre += l1off * xre;
787 *yim += l1off * xim;
788 }
789 }
790}
791
811inline void control(const PetscInt dim, const int it, const int n, const int i, const int np, const int ip, const int stridei, const int strideip, const double* xptr, const double pt, const double qt, double* yre, double* yim){
812 /* \rho(ik+1..,ik'..) term */
813 if (i < n-1) {
814 int itx = it + stridei;
815 double xre = xptr[itx];
816 double xim = xptr[itx+dim];
817 double sq = sqrt(i + 1);
818 *yre += sq * ( pt * xim + qt * xre);
819 *yim += sq * ( - pt * xre + qt * xim);
820 }
821 /* \rho(ik..,ik'+1..) */
822 if (ip < np-1) {
823 int itx = it + strideip;
824 double xre = xptr[itx];
825 double xim = xptr[itx+dim];
826 double sq = sqrt(ip + 1);
827 *yre += sq * ( -pt * xim + qt * xre);
828 *yim += sq * ( pt * xre + qt * xim);
829 }
830 /* \rho(ik-1..,ik'..) */
831 if (i > 0) {
832 int itx = it - stridei;
833 double xre = xptr[itx];
834 double xim = xptr[itx+dim];
835 double sq = sqrt(i);
836 *yre += sq * ( pt * xim - qt * xre);
837 *yim += sq * (- pt * xre - qt * xim);
838 }
839 /* \rho(ik..,ik'-1..) */
840 if (ip > 0) {
841 int itx = it - strideip;
842 double xre = xptr[itx];
843 double xim = xptr[itx+dim];
844 double sq = sqrt(ip);
845 *yre += sq * (- pt * xim - qt * xre);
846 *yim += sq * ( pt * xre - qt * xim);
847 }
848}
849
850
868inline void control_T(const PetscInt dim, const int it, const int n, const int i, const int np, const int ip, const int stridei, const int strideip, const double* xptr, const double pt, const double qt, double* yre, double* yim){
869 /* \rho(ik+1..,ik'..) term */
870 if (i > 0) {
871 int itx = it - stridei;
872 double xre = xptr[itx];
873 double xim = xptr[itx+dim];
874 double sq = sqrt(i);
875 *yre += sq * ( - pt * xim + qt * xre);
876 *yim += sq * ( pt * xre + qt * xim);
877 }
878 /* \rho(ik..,ik'+1..) */
879 if (ip > 0) {
880 int itx = it - strideip;
881 double xre = xptr[itx];
882 double xim = xptr[itx+dim];
883 double sq = sqrt(ip);
884 *yre += sq * ( pt * xim + qt * xre);
885 *yim += sq * ( -pt * xre + qt * xim);
886 }
887 /* \rho(ik-1..,ik'..) */
888 if (i < n-1) {
889 int itx = it + stridei;
890 double xre = xptr[itx];
891 double xim = xptr[itx+dim];
892 double sq = sqrt(i+1);
893 *yre += sq * (- pt * xim - qt * xre);
894 *yim += sq * ( pt * xre - qt * xim);
895 }
896 /* \rho(ik..,ik'-1..) */
897 if (ip < np-1) {
898 int itx = it + strideip;
899 double xre = xptr[itx];
900 double xim = xptr[itx+dim];
901 double sq = sqrt(ip+1);
902 *yre += sq * (+ pt * xim - qt * xre);
903 *yim += sq * (- pt * xre - qt * xim);
904 }
905}
Config
Configuration class containing all validated settings.
Definition config.hpp:57
MasterEq
Implementation of the real-valued right-hand-side (RHS) system matrix of the quantum dynamical equati...
Definition mastereq.hpp:92
MasterEq::addT1
bool addT1
Definition mastereq.hpp:113
MasterEq::addT2
bool addT2
Flags for including T1 decay and T2 dephasing Lindblad operators.
Definition mastereq.hpp:113
MasterEq::Ac_vec
std::vector< Mat > Ac_vec
Definition mastereq.hpp:104
MasterEq::decoherence_type
DecoherenceType decoherence_type
Type of Lindblad operators to include (NONE means Schroedinger equation)
Definition mastereq.hpp:132
MasterEq::getRHS
Mat getRHS()
Retrieves the right-hand-side system matrix.
Definition mastereq.cpp:694
MasterEq::nlevels
std::vector< size_t > nlevels
Number of levels per oscillator.
Definition mastereq.hpp:129
MasterEq::mpirank_petsc
int mpirank_petsc
Rank of PETSc's communicator.
Definition mastereq.hpp:116
MasterEq::MasterEq
MasterEq()
Definition mastereq.cpp:3
MasterEq::Bd
Mat Bd
Definition mastereq.hpp:106
MasterEq::usematfree
bool usematfree
Flag for using matrix-free solver.
Definition mastereq.hpp:131
MasterEq::dim_rho
PetscInt dim_rho
Dimension of Hilbert space = N.
Definition mastereq.hpp:96
MasterEq::quietmode
bool quietmode
Flag for quiet mode operation.
Definition mastereq.hpp:124
MasterEq::nessential
std::vector< size_t > nessential
Number of essential levels per oscillator.
Definition mastereq.hpp:130
MasterEq::crosskerr
std::vector< double > crosskerr
Cross-Kerr coefficients (rad/time) for ZZ-coupling .
Definition mastereq.hpp:110
MasterEq::isu
IS isu
Definition mastereq.hpp:122
MasterEq::hamiltonian_file_Hc
std::optional< std::string > hamiltonian_file_Hc
Filename if a custom control Hamiltonians are read from file.
Definition mastereq.hpp:126
MasterEq::isv
IS isv
Vector strides for accessing real and imaginary parts u=Re(x), v=Im(x)
Definition mastereq.hpp:122
MasterEq::localsize_u
PetscInt localsize_u
Size of local sub vector u or v in state x=[u,v].
Definition mastereq.hpp:118
MasterEq::getDim
PetscInt getDim()
Retrieves the dimension of the vectorized system.
Definition mastereq.hpp:181
MasterEq::getNOscillators
size_t getNOscillators()
Retrieves the number of oscillators in the system.
Definition mastereq.hpp:174
MasterEq::getDimRho
PetscInt getDimRho()
Retrieves the dimension of the density matrix.
Definition mastereq.hpp:195
MasterEq::RHSctx
MatShellCtx RHSctx
Context containing data for applying the RHS to a state.
Definition mastereq.hpp:102
MasterEq::oscil_vec
Oscillator ** oscil_vec
Array of pointers to oscillator objects.
Definition mastereq.hpp:99
MasterEq::Bd_vec
std::vector< Mat > Bd_vec
Definition mastereq.hpp:108
MasterEq::hamiltonian_file_Hsys
std::optional< std::string > hamiltonian_file_Hsys
Filename if a custom system Hamiltonian is read from file.
Definition mastereq.hpp:125
MasterEq::Ad
Mat Ad
Definition mastereq.hpp:106
MasterEq::getOscillator
Oscillator * getOscillator(const size_t i)
Retrieves the i-th oscillator.
Definition mastereq.hpp:167
MasterEq::eta
std::vector< double > eta
Frequency differences in rotating frame (rad/time) for dipole-dipole coupling.
Definition mastereq.hpp:112
MasterEq::population
void population(const Vec x, std::vector< double > &population_com)
Computes population of the full composite system.
Definition mastereq.cpp:2950
MasterEq::set_RHS_MatMult_operation
void set_RHS_MatMult_operation()
Pass MatMult operations for applying the RHS to PETSc.
Definition mastereq.cpp:724
MasterEq::ilow
PetscInt ilow
First index of the local sub vector u,v.
Definition mastereq.hpp:119
MasterEq::noscillators
int noscillators
Number of oscillators in the system.
Definition mastereq.hpp:98
MasterEq::initSparseMatSolver
void initSparseMatSolver()
Initializes matrices needed for the sparse matrix solver.
Definition mastereq.cpp:206
MasterEq::assemble_RHS
int assemble_RHS(const double t)
Assembles the real-valued system matrix (RHS) at time t.
Definition mastereq.cpp:671
MasterEq::expectedEnergy
double expectedEnergy(const Vec x)
Computes expected energy of the full composite system.
Definition mastereq.cpp:2911
MasterEq::Ad_vec
std::vector< Mat > Ad_vec
Definition mastereq.hpp:107
MasterEq::RHS
Mat RHS
MatShell for real-valued, vectorized system matrix (size 2N^2 x 2N^2)
Definition mastereq.hpp:101
MasterEq::mpisize_petsc
int mpisize_petsc
Size of PETSc's communicator.
Definition mastereq.hpp:117
MasterEq::Jkl
std::vector< double > Jkl
Dipole-dipole coupling coefficients (rad/time), multiplies .
Definition mastereq.hpp:111
MasterEq::mpirank_world
int mpirank_world
Rank of global MPI communicator.
Definition mastereq.hpp:115
MasterEq::aux
Vec aux
Auxiliary vector for computations.
Definition mastereq.hpp:123
MasterEq::dim
PetscInt dim
Dimension of full vectorized system: N^2 if Lindblad, N if Schroedinger.
Definition mastereq.hpp:95
MasterEq::getDimEss
PetscInt getDimEss()
Retrieves the dimension of the essential level system.
Definition mastereq.hpp:188
MasterEq::dim_ess
PetscInt dim_ess
Dimension of essential levels = N_e.
Definition mastereq.hpp:97
MasterEq::iupp
PetscInt iupp
Last index (+1) of the local sub vector u,v.
Definition mastereq.hpp:120
MasterEq::Bc_vec
std::vector< Mat > Bc_vec
Definition mastereq.hpp:105
MasterEq::setControlAmplitudes
void setControlAmplitudes(const Vec x)
Pass control parameters from global design vector to each oscillator.
Definition mastereq.cpp:707
MasterEq::~MasterEq
~MasterEq()
Definition mastereq.cpp:172
MasterEq::compute_dRHS_dParams
void compute_dRHS_dParams(const double t, const Vec x, const Vec x_bar, const double alpha, Vec grad)
Computes gradient of RHS with respect to control parameters.
Definition mastereq.cpp:697
Oscillator
Quantum oscillator (multi-level qubit) with control capabilities.
Definition oscillator.hpp:37
defs.hpp
Core type definitions and enumerations for Quandary quantum optimal control.
DecoherenceType
DecoherenceType
Available Lindblad operator types for open quantum systems, or NONE for closed quantum systems.
Definition defs.hpp:31
compute_dRHS_dParams_sparsemat
void compute_dRHS_dParams_sparsemat(const double t, const Vec x, const Vec x_bar, const double alpha, Vec grad, std::vector< size_t > &nlevels, IS isu, IS isv, std::vector< Mat > &Ac_vec, std::vector< Mat > &Bc_vec, Vec aux, Oscillator **oscil_vec)
: Sparse-matrix version to compute gradient of RHS with respect to parameters
Definition mastereq.cpp:939
Jkl_coupling_T
void Jkl_coupling_T(const PetscInt dim, const int it, const int ni, const int nj, const int nip, const int njp, const int i, const int ip, const int j, const int jp, const int stridei, const int strideip, const int stridej, const int stridejp, const double *xptr, const double Jij, const double cosij, const double sinij, double *yre, double *yim)
Transpose of dipole-dipole coupling for adjoint computations.
Definition mastereq.hpp:694
applyRHS_sparsemat
int applyRHS_sparsemat(Mat RHS, Vec x, Vec y)
Sparse matrix MatMult.
Definition mastereq.cpp:757
compute_dRHS_dParams_matfree
void compute_dRHS_dParams_matfree(const PetscInt dim, const double t, const Vec x, const Vec x_bar, const double alpha, Vec grad, std::vector< size_t > &nlevels, DecoherenceType decoherence_type, Oscillator **oscil_vec)
: Matrix free version to compute gradient of RHS with respect to parameters
Definition mastereq.cpp:984
applyRHS_matfree_1Osc
int applyRHS_matfree_1Osc(Mat RHS, Vec x, Vec y)
Matrix-vector products to apply the RHS system matrix to a state vector.
Definition mastereq.cpp:2991
L2
double L2(const double dephase0, const int i0, const int i0p)
Inline for Matrix-free RHS to compute L2 dephasing coefficient for 1 oscillator.
Definition mastereq.hpp:332
applyRHS_matfree_transpose_1Osc
int applyRHS_matfree_transpose_1Osc(Mat RHS, Vec x, Vec y)
Transpose matrix-free MatMult for 1 oscillator.
Definition mastereq.cpp:3020
applyRHS_matfree_transpose_5Osc
int applyRHS_matfree_transpose_5Osc(Mat RHS, Vec x, Vec y)
Transpose matrix-free MatMult for 5 oscillators.
Definition mastereq.cpp:3232
applyRHS_matfree_transpose_2Osc
int applyRHS_matfree_transpose_2Osc(Mat RHS, Vec x, Vec y)
Transpose matrix-free MatMult for 2 oscillators.
Definition mastereq.cpp:3077
applyRHS_matfree_4Osc
int applyRHS_matfree_4Osc(Mat RHS, Vec x, Vec y)
Matrix-free MatMult for 4 oscillators.
Definition mastereq.cpp:3155
applyRHS_matfree_3Osc
int applyRHS_matfree_3Osc(Mat RHS, Vec x, Vec y)
Matrix-free MatMult for 3 oscillators.
Definition mastereq.cpp:3105
dRHSdp_getcoeffs
void dRHSdp_getcoeffs(const PetscInt dim, const int it, const int n, const int np, const int i, const int ip, const int stridei, const int strideip, const double *xptr, double *res_p_re, double *res_p_im, double *res_q_re, double *res_q_im)
Inline for Matrix-free RHS for gradient updates.
Definition mastereq.hpp:546
Jkl_coupling
void Jkl_coupling(const PetscInt dim, const int it, const int ni, const int nj, const int nip, const int njp, const int i, const int ip, const int j, const int jp, const int stridei, const int strideip, const int stridej, const int stridejp, const double *xptr, const double Jij, const double cosij, const double sinij, double *yre, double *yim)
Matrix-free solver inline for dipole-dipole coupling between oscillators.
Definition mastereq.hpp:625
L1decay
void L1decay(const PetscInt dim, const int it, const int n, const int i, const int ip, const int stridei, const int strideip, const double *xptr, const double decayi, double *yre, double *yim)
Matrix-free solver inline for off-diagonal L1 decay term.
Definition mastereq.hpp:751
H_crosskerr
double H_crosskerr(const double xi01, const int a, const int b)
Inline for Matrix-free RHS, Computes cross-Kerr coupling between 2 oscillators.
Definition mastereq.hpp:394
TensorGetIndex
int TensorGetIndex(const int nlevels0, const int i0, const int i0p)
Inline for Matrix-free RHS to compute tensor product index for 1 oscillator system.
Definition mastereq.hpp:356
control
void control(const PetscInt dim, const int it, const int n, const int i, const int np, const int ip, const int stridei, const int strideip, const double *xptr, const double pt, const double qt, double *yre, double *yim)
Matrix-free solver inline for control terms.
Definition mastereq.hpp:811
applyRHS_matfree_5Osc
int applyRHS_matfree_5Osc(Mat RHS, Vec x, Vec y)
Matrix-free MatMult for 5 oscillators.
Definition mastereq.cpp:3206
H_selfkerr
double H_selfkerr(const double xi0, const int a)
Inline for Matrix-free RHS to Compute self-Kerr coefficient for 1 oscillator.
Definition mastereq.hpp:320
applyRHS_matfree_2Osc
int applyRHS_matfree_2Osc(Mat RHS, Vec x, Vec y)
Matrix-free MatMult for 2 oscillators.
Definition mastereq.cpp:3047
applyRHS_matfree_transpose_4Osc
int applyRHS_matfree_transpose_4Osc(Mat RHS, Vec x, Vec y)
Transpose matrix-free MatMult for 4 oscillators.
Definition mastereq.cpp:3181
control_T
void control_T(const PetscInt dim, const int it, const int n, const int i, const int np, const int ip, const int stridei, const int strideip, const double *xptr, const double pt, const double qt, double *yre, double *yim)
Matrix-free Transpose of control terms for adjoint computations.
Definition mastereq.hpp:868
L1diag
double L1diag(const double decay0, const int i0, const int i0p)
Inline for Matrix-free RHS to compute L1 decay coefficient for 1 oscillator.
Definition mastereq.hpp:344
applyRHS_matfree_transpose_3Osc
int applyRHS_matfree_transpose_3Osc(Mat RHS, Vec x, Vec y)
Transpose matrix-free MatMult for 3 oscillators.
Definition mastereq.cpp:3131
applyRHS_sparsemat_transpose
int applyRHS_sparsemat_transpose(Mat RHS, Vec x, Vec y)
Transpose sparse matrix MatMult.
Definition mastereq.cpp:848
H_detune
double H_detune(const double detuning0, const int a)
Inline for Matrix-free RHS to Compute detuning for 1 oscillator.
Definition mastereq.hpp:309
L1decay_T
void L1decay_T(const PetscInt dim, const int it, const int i, const int ip, const int stridei, const int strideip, const double *xptr, const double decayi, double *yre, double *yim)
Matrix-free inline Transpose of off-diagonal L1 decay for adjoint computations.
Definition mastereq.hpp:779
MatShellCtx
Matrix shell context containing data needed for applying the right-hand-side (RHS) system matrix to a...
Definition mastereq.hpp:21
MatShellCtx::Ad_vec
std::vector< Mat > Ad_vec
Vector of real parts of dipole-dipole coupling system matrices.
Definition mastereq.hpp:37
MatShellCtx::oscil_vec
Oscillator ** oscil_vec
Array of pointers to the oscillators.
Definition mastereq.hpp:25
MatShellCtx::aux
Vec * aux
Auxiliary vector for computations.
Definition mastereq.hpp:41
MatShellCtx::Jkl
std::vector< double > Jkl
Dipole-dipole coupling strength.
Definition mastereq.hpp:27
MatShellCtx::eta
std::vector< double > eta
Frequency differences of the rotating frames.
Definition mastereq.hpp:28
MatShellCtx::nlevels
std::vector< size_t > nlevels
Number of levels per oscillator.
Definition mastereq.hpp:23
MatShellCtx::crosskerr
std::vector< double > crosskerr
Cross-Kerr coupling coefficients.
Definition mastereq.hpp:26
MatShellCtx::time
double time
Current time.
Definition mastereq.hpp:42
MatShellCtx::Ad
Mat * Ad
Real parts of time-independent system matrix.
Definition mastereq.hpp:35
MatShellCtx::control_Im
std::vector< double > control_Im
Imaginary parts of control pulse .
Definition mastereq.hpp:32
MatShellCtx::dim
PetscInt dim
Dimension of full vectorized system: N^2 if Lindblad, N if Schroedinger.
Definition mastereq.hpp:22
MatShellCtx::control_Re
std::vector< double > control_Re
Real parts of control pulse .
Definition mastereq.hpp:31
MatShellCtx::Bd_coeffs
std::vector< double > Bd_coeffs
Definition mastereq.hpp:39
MatShellCtx::Ac_vec
std::vector< Mat > Ac_vec
Vector of real parts of control matrices per oscillator.
Definition mastereq.hpp:33
MatShellCtx::isu
IS * isu
Definition mastereq.hpp:24
MatShellCtx::Bd
Mat * Bd
Imaginary parts of time-independent system matrix.
Definition mastereq.hpp:36
MatShellCtx::Ad_coeffs
std::vector< double > Ad_coeffs
Definition mastereq.hpp:40
MatShellCtx::decoherence_type
DecoherenceType decoherence_type
Type of Lindblad operators to include.
Definition mastereq.hpp:29
MatShellCtx::addT1
bool addT1
Definition mastereq.hpp:30
MatShellCtx::Bd_vec
std::vector< Mat > Bd_vec
Vector of imaginary parts of dipole-dipole coupling system matrices.
Definition mastereq.hpp:38
MatShellCtx::Bc_vec
std::vector< Mat > Bc_vec
Vector of imaginary parts of control matrices per oscillator.
Definition mastereq.hpp:34
Generated by doxygen 1.9.8