Files
libquantum/qtime.c
2016-10-27 04:32:19 +09:00

187 lines
4.5 KiB
C

/* qtime.c: Time evolution of a quantum system
Copyright 2006-2013 Hendrik Weimer
This file is part of libquantum
libquantum is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published
by the Free Software Foundation; either version 3 of the License,
or (at your option) any later version.
libquantum is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with libquantum; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
MA 02110-1301, USA
*/
#include <math.h>
#include <string.h>
#include <stdio.h>
#include "qtime.h"
#include "qureg.h"
#include "complex.h"
#include "config.h"
/* Forth-order Runge-Kutta
Flags: QUANTUM_RK4_NODELETE: Do not delete quantum_reg returned by H
QUANTUM_RK4_IMAGINARY: Imaginary time evolution */
void
quantum_rk4(quantum_reg *reg, double t, double dt,
quantum_reg H(MAX_UNSIGNED, double), int flags)
{
quantum_reg k, out, tmp;
double r = 0;
int i;
void *hash;
int hashw;
COMPLEX_FLOAT step = dt;
hash = reg->hash;
reg->hash = 0;
hashw = reg->hashw;
reg->hashw = 0;
if(!(flags & QUANTUM_RK4_IMAGINARY))
step *= IMAGINARY;
/* k1 */
k = quantum_matrix_qureg(H, t, reg, flags & QUANTUM_RK4_NODELETE);
quantum_scalar_qureg(-step/2.0, &k);
tmp = quantum_vectoradd(reg, &k);
quantum_scalar_qureg(1.0/3.0, &k);
out = quantum_vectoradd(reg, &k);
quantum_delete_qureg(&k);
/* k2 */
k = quantum_matrix_qureg(H, t+dt/2.0, &tmp, flags & QUANTUM_RK4_NODELETE);
quantum_delete_qureg(&tmp);
quantum_scalar_qureg(-step/2.0, &k);
tmp = quantum_vectoradd(reg, &k);
quantum_scalar_qureg(2.0/3.0, &k);
quantum_vectoradd_inplace(&out, &k);
quantum_delete_qureg(&k);
/* k3 */
k = quantum_matrix_qureg(H, t+dt/2.0, &tmp, flags & QUANTUM_RK4_NODELETE);
quantum_delete_qureg(&tmp);
quantum_scalar_qureg(-step, &k);
tmp = quantum_vectoradd(reg, &k);
quantum_scalar_qureg(1.0/3.0, &k);
quantum_vectoradd_inplace(&out, &k);
quantum_delete_qureg(&k);
/* k4 */
k = quantum_matrix_qureg(H, t+dt, &tmp, flags & QUANTUM_RK4_NODELETE);
quantum_delete_qureg(&tmp);
quantum_scalar_qureg(-step/6.0, &k);
quantum_vectoradd_inplace(&out, &k);
quantum_delete_qureg(&k);
quantum_delete_qureg(reg);
/* Normalize quantum register */
if(flags & QUANTUM_RK4_IMAGINARY)
{
for(i=0; i<out.size; i++)
r += quantum_prob(out.amplitude[i]);
quantum_scalar_qureg(sqrt(1.0/r), &out);
}
out.hash = hash;
out.hashw = hashw;
*reg = out;
}
/* Adaptive Runge-Kutta. Stores the new stepsize in dt and returns the
stepsize actually used. For further details, see Press et al.,
Numerical Recipes in C (Second Edition, CUP, 1992), Sec. 16.3 */
double
quantum_rk4a(quantum_reg *reg, double t, double *dt, double epsilon,
quantum_reg H(MAX_UNSIGNED, double), int flags)
{
quantum_reg reg2, old;
double delta, r, dtused;
int i;
void *hash;
int hashw;
hash = reg->hash;
reg->hash = 0;
hashw = reg->hashw;
reg->hashw = 0;
quantum_copy_qureg(reg, &old);
quantum_copy_qureg(reg, &reg2);
do
{
quantum_rk4(reg, t, *dt, H, flags);
quantum_rk4(&reg2, t, *dt/2.0, H, flags);
quantum_rk4(&reg2, t+*dt/2.0, *dt/2.0, H, flags);
delta = 0;
for(i=0;i<reg->size;i++)
{
r = 2*sqrt(quantum_prob(reg->amplitude[i] - reg2.amplitude[i])/
quantum_prob(reg->amplitude[i] + reg2.amplitude[i]));
if(r > delta)
delta = r;
}
dtused = *dt;
if(delta < epsilon)
*dt *= 0.9*pow(epsilon/delta, 0.2);
else
*dt *= 0.9*pow(epsilon/delta, 0.25);
if(*dt > 4*dtused)
*dt = 4*dtused;
else if(*dt < 0.25*dtused)
*dt = 0.25*dtused;
if(delta > epsilon)
{
memcpy(reg->amplitude, old.amplitude,
reg->size*sizeof(COMPLEX_FLOAT));
memcpy(reg2.amplitude, old.amplitude,
reg->size*sizeof(COMPLEX_FLOAT));
if(reg->state && old.state)
memcpy(reg->state, old.state, reg->size*sizeof(MAX_UNSIGNED));
if(reg2.state && old.state)
memcpy(reg2.state, old.state, reg->size*sizeof(MAX_UNSIGNED));
}
} while(delta > epsilon);
reg->hash = hash;
reg->hashw = hashw;
quantum_delete_qureg(&old);
quantum_delete_qureg(&reg2);
return dtused;
}