updated libquantum 1.1.1 source files

CHANGES

@@ -1,3 +1,12 @@
+libquantum 1.1.1:
+	- Added support for ground state calculations
+	- Added example program simulating the tranverse Ising chain
+	- Optimized memory layout for quantum registers (thanks to
+          Acumem, breaks backward compatiblity)
+	- Added OpenMP support
+	- Improved C99 compatibility for complex numbers
+	- Improved support of double precision arithmetic
+
 libquantum 1.1.0:
 	- Added exact diagonlization based on LAPACK
 	- Added flag in quantum_rk4 to preserve qureg (breaks backward 

Makefile.in

@@ -48,8 +48,8 @@ LIBTOOL=@LIBTOOL@
 
 # Flags passed to C compiler
 
-CFLAGS=@CFLAGS@
+CFLAGS=@CFLAGS@ @OPENMP_CFLAGS@
-LDFLAGS=-rpath $(LIBDIR) -version-info 7:0:0
+LDFLAGS=-rpath $(LIBDIR) -version-info 8:0:0
 
 # Dependencies
 
@@ -57,11 +57,11 @@ all:	libquantum.la
 
 libquantum.la: complex.lo measure.lo matrix.lo gates.lo qft.lo classic.lo \
 	qureg.lo decoherence.lo oaddn.lo omuln.lo expn.lo qec.lo version.lo \
-	objcode.lo density.lo error.lo qtime.lo lapack.lo Makefile
+	objcode.lo density.lo error.lo qtime.lo lapack.lo energy.lo Makefile
 	$(LIBTOOL) --mode=link $(CC) $(LDFLAGS) -o libquantum.la complex.lo \
 	measure.lo matrix.lo gates.lo oaddn.lo omuln.lo expn.lo qft.lo \
 	classic.lo qureg.lo decoherence.lo qec.lo version.lo objcode.lo \
-	density.lo error.lo qtime.lo lapack.lo @LIBS@
+	density.lo error.lo qtime.lo lapack.lo energy.lo @LIBS@
 
 complex.lo: complex.c complex.h config.h Makefile
 	$(LIBTOOL) --mode=compile $(CC) $(CFLAGS) -c complex.c
@@ -123,6 +123,9 @@ qtime.lo: qtime.c qtime.h qureg.h Makefile
 lapack.lo: lapack.c lapack.h matrix.h qureg.h config.h error.h Makefile
 	$(LIBTOOL) --mode=compile $(CC) $(CFLAGS) -c lapack.c
 
+energy.lo: energy.c energy.h qureg.h config.h error.h Makefile
+	$(LIBTOOL) --mode=compile $(CC) $(CFLAGS) -c energy.c
+
 # Autoconf stuff
 
 Makefile: config.status Makefile.in aclocal.m4 config.h.in types.h.in \
@@ -140,7 +143,7 @@ config.status: configure
 
 # Build demos of Shor's and Grover's algorithms
 
-demos: shor grover
+demos: shor grover ising
 
 shor: libquantum.la shor.c Makefile
 	$(LIBTOOL) --mode=link $(CC) $(CFLAGS) -o shor shor.c -I./ -lquantum \
@@ -150,6 +153,10 @@ grover: libquantum.la grover.c Makefile
 	$(LIBTOOL) --mode=link $(CC) $(CFLAGS) -o grover grover.c -I./ \
 	-lquantum -static -lm
 
+ising:	libquantum.la ising.c Makefile
+	$(LIBTOOL) --mode=link $(CC) $(CFLAGS) -o ising ising.c -I./ -lquantum \
+	-static -lm
+
 # Quantum object code tools
 
 quobtools: quobprint quobdump

complex.c

@@ -22,26 +22,27 @@
 */
 
 #include <math.h>
+#include <complex.h>
 
 #include "complex.h"
 #include "config.h"
 
 /* Return the complex conjugate of a complex number */
 
-COMPLEX_FLOAT
+/*COMPLEX_FLOAT
 quantum_conj(COMPLEX_FLOAT a)
 {
-  float r, i;
+  REAL_FLOAT r, i;
 
   r = quantum_real(a);
   i = quantum_imag(a);
 
   return r - IMAGINARY * i;
-}
+  }*/
 
 /* Calculate the square of a complex number (i.e. the probability) */
 
-float
+double
 quantum_prob(COMPLEX_FLOAT a)
 {
   return quantum_prob_inline(a);
@@ -49,7 +50,7 @@ quantum_prob(COMPLEX_FLOAT a)
 
 /* Calculate e^(i * phi) */
 
-COMPLEX_FLOAT quantum_cexp(float phi)
+COMPLEX_FLOAT quantum_cexp(REAL_FLOAT phi)
 {
   return cos(phi) + IMAGINARY * sin(phi);
 }

complex.h

@@ -25,37 +25,22 @@
 
 #define __COMPLEX_H
 
+#include <complex.h>
 #include "config.h"
 
-extern COMPLEX_FLOAT quantum_conj(COMPLEX_FLOAT a);
+#define quantum_conj(z) (conj(z))
+#define quantum_real(z) (creal(z))
+#define quantum_imag(z) (cimag(z))
 
-extern float quantum_prob (COMPLEX_FLOAT a);
+extern double quantum_prob (COMPLEX_FLOAT a);
-extern COMPLEX_FLOAT quantum_cexp(float phi);
+extern COMPLEX_FLOAT quantum_cexp(REAL_FLOAT phi);
-
-/* Return the real part of a complex number */
-
-static inline float
-quantum_real(COMPLEX_FLOAT a)
-{
-  float *p = (float *) &a;
-  return p[0];
-}
-
-/* Return the imaginary part of a complex number */
-
-static inline float
-quantum_imag(COMPLEX_FLOAT a)
-{
-  float *p = (float *) &a;
-  return p[1];
-}
 
 /* Calculate the square of a complex number (i.e. the probability) */
 
-static inline float 
+static inline double
 quantum_prob_inline(COMPLEX_FLOAT a)
 {
-  float r, i;
+  REAL_FLOAT r, i;
 
   r = quantum_real(a);
   i = quantum_imag(a);

config.h.in

@@ -94,4 +94,7 @@
 /* Define to 1 if you have LAPACK */
 #undef HAVE_LIBLAPACK
 
+/* Define to 1 if using double precision */
+#undef USE_DOUBLE
+
 #include "types.h"

configure.in

@@ -1,7 +1,7 @@
 # configure.in: Process this file with autoconf to produce a configure
 #               script.
 # 
-# Copyright 2003-2005 Bjoern Butscher, Hendrik Weimer
+# Copyright 2003-2013 Bjoern Butscher, Hendrik Weimer
 #
 # This file is part of libquantum
 #
@@ -20,7 +20,7 @@
 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 # MA 02110-1301, USA
 
-AC_INIT([libquantum], [1.1.0], [libquantum@libquantum.de])
+AC_INIT([libquantum], [1.1.1], [libquantum@libquantum.de])
 AC_CONFIG_SRCDIR([classic.c])
 AC_CONFIG_HEADER([config.h])
 
@@ -35,7 +35,7 @@ AC_PROG_LIBTOOL
 # Checks for libraries.
 AC_CHECK_LIB([m], [sqrt])
 AC_ARG_WITH(lapack,
-	[  --with-lapack        LAPACK support [default=yes]],
+	[  --with-lapack           LAPACK support [[default=yes]]],
 	[if test $withval = "yes"
 	then 
 	     AC_CHECK_LIB([lapack], [cheev_])
@@ -75,11 +75,16 @@ fi
 # Check for complex number type
 AC_ARG_WITH([complex-type],
 	[  --with-complex-type=ARG type for complex numbers], 
-	[CF_TYPE=$withval
+	[CF_TYPE=$withval], [CF_TYPE="none"
-], [CF_TYPE="none"
+	AC_CHECK_TYPE([float complex], [AC_DEFINE([COMPLEX_FLOAT], 
+					[float complex])
+			  		CF_TYPE="float complex"])])
+if test "$CF_TYPE" = "none"
+then
 	AC_CHECK_TYPE([float _Complex], [AC_DEFINE([COMPLEX_FLOAT],
-					[float _Complex])
+			     		[float _Complex])
-			     CF_TYPE="float _Complex"])])
+			  		CF_TYPE="float _Complex"])
+fi
 if test "$CF_TYPE" = "none"
 then
 	AC_CHECK_TYPE([__complex__ float], [AC_DEFINE([COMPLEX_FLOAT], 
@@ -91,6 +96,14 @@ then
 	AC_MSG_ERROR([No complex number type!])
 fi
 
+AC_MSG_CHECKING([for corresponding real data type])
+AC_RUN_IFELSE(
+	[AC_LANG_PROGRAM([], [return sizeof($CF_TYPE) == 2*sizeof(double)])],
+	[RF_TYPE="float"],
+	[RF_TYPE="double"; AC_DEFINE(USE_DOUBLE)], [float])
+AC_MSG_RESULT($RF_TYPE)
+
+
 # Check for the imaginary unit
 AC_MSG_CHECKING([for the imaginary unit])
 AC_ARG_WITH([imaginary],
@@ -108,9 +121,13 @@ then
 fi
 AC_MSG_RESULT($I)
 
+# Check for OpenMP support
+AC_OPENMP
+
 # Substitute fields in quantum.h.in and types.h
 AC_SUBST(MU_TYPE)
 AC_SUBST(CF_TYPE)
+AC_SUBST(RF_TYPE)
 AC_SUBST(I)
 
 # Profiling check
@@ -120,9 +137,6 @@ AC_ARG_ENABLE(profiling,
 	    then CFLAGS="$CFLAGS -pg -fprofile-arcs -ftest-coverage"
 	fi], [])
 
-# Disable LAPACK check
-
-
 # Enable -Wall for gcc
 if test $CC = "gcc"
 then

decoherence.c

@@ -114,13 +114,13 @@ quantum_decohere(quantum_reg *reg)
 
 	  for(j=0; j<reg->width; j++)
 	    {
-	      if(reg->node[i].state & ((MAX_UNSIGNED) 1 << j))
+	      if(reg->state[i] & ((MAX_UNSIGNED) 1 << j))
 		angle += nrands[j];
 	      else
 		angle -= nrands[j];
 	    }
 
-	  reg->node[i].amplitude *= quantum_cexp(angle);
+	  reg->amplitude[i] *= quantum_cexp(angle);
 	  
 	}
       free(nrands);

density.c

@@ -66,7 +66,8 @@ quantum_new_density_op(int num, float *prob, quantum_reg *reg)
 
   reg[0].size = 0;
   reg[0].width = 0;
-  reg[0].node = 0;
+  reg[0].state = 0;
+  reg[0].amplitude = 0;
   reg[0].hash = 0;
 
   for(i=1; i<num; i++)
@@ -78,7 +79,8 @@ quantum_new_density_op(int num, float *prob, quantum_reg *reg)
 
       reg[i].size = 0;
       reg[i].width = 0;
-      reg[i].node = 0;
+      reg[i].state = 0;
+      reg[i].amplitude = 0;
       reg[i].hash = 0;
     }
 
@@ -133,8 +135,8 @@ quantum_reduced_density_op(int pos, quantum_density_op *rho)
   
       for(j=0; j<rho->reg[i].size; j++)
 	{
-	  if(!(rho->reg[i].node[j].state & pos2))
+	  if(!(rho->reg[i].state[j] & pos2))
-	    p0 += quantum_prob_inline(rho->reg[i].node[j].amplitude);
+	    p0 += quantum_prob_inline(rho->reg[i].amplitude[j]);
 	}
 
       rho->prob[i] = ptmp * p0;
@@ -176,8 +178,8 @@ quantum_density_matrix(quantum_density_op *rho)
 	      l1 = quantum_get_state(i, rho->reg[k]);
 	      l2 = quantum_get_state(j, rho->reg[k]);
 	      if((l1 > -1) && (l2 > -1))
-		M(m, i, j) += rho->prob[k] * rho->reg[k].node[l2].amplitude
+		M(m, i, j) += rho->prob[k] * rho->reg[k].amplitude[l2]
-		  * quantum_conj(rho->reg[k].node[l1].amplitude);
+		  * quantum_conj(rho->reg[k].amplitude[l1]);
 	    }
 	}
     }
@@ -246,14 +248,14 @@ quantum_purity(quantum_density_op *rho)
 	      /* quantum_dot_product makes sure that rho->reg[j] has a
 		 correct hash table */
 
-	      l = quantum_get_state(rho->reg[i].node[k].state, rho->reg[j]);
+	      l = quantum_get_state(rho->reg[i].state[k], rho->reg[j]);
 
 	      /* Compute p_i p_j <k|\psi_iX\psi_i|\psi_jX\psi_j|k> */
 	      
 	      if(l > -1)
 		g = rho->prob[i] * rho->prob[j] * dp 
-		  * rho->reg[i].node[k].amplitude
+		  * rho->reg[i].amplitude[k]
-		  * quantum_conj(rho->reg[j].node[l].amplitude);
+		  * quantum_conj(rho->reg[j].amplitude[l]);
 	      else
 		g = 0;
 

energy.c

@@ -0,0 +1,320 @@
+/* energy.c: Compute energetic properties of quantum systems
+
+   Copyright 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 <float.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "energy.h"
+#include "qureg.h"
+#include "qtime.h"
+#include "complex.h"
+
+extern void dstevd_(char *jobz, int *n, double *d, double *e, double *z, 
+		    int *ldz, double *work, int *lwork, int *iwork, int *liwork,
+		    int *info);
+
+
+/* Modified Lanczos algorithm that iterates over a series of 2x2
+   matrix diagonalzations [E. Dagotto & A. Moreo, Phys. Rev. D 31, 865
+   (1985)] */
+
+double 
+quantum_lanczos_modified(quantum_reg H(MAX_UNSIGNED, double), double epsilon, 
+			 quantum_reg *reg)
+{
+  double E0=DBL_MAX, Eold=DBL_MAX, E1, E2, t;
+  quantum_reg tmp, tmp2;
+  int i;
+  COMPLEX_FLOAT h01;
+  double h00, h11;
+
+  for(i=0; i<reg->size; i++)
+    {
+      quantum_normalize(reg);
+
+      tmp = quantum_matrix_qureg(H, 0, reg, QUANTUM_RK4_NODELETE);
+
+      h00 = quantum_real(quantum_dot_product(&tmp, reg));
+
+      E0 = h00;
+
+      if(fabs(E0-Eold)<epsilon)
+	return E0;
+
+      Eold = E0;
+
+      quantum_copy_qureg(reg, &tmp2);
+      quantum_scalar_qureg(-h00, &tmp2);
+
+      quantum_vectoradd_inplace(&tmp, &tmp2);
+
+      quantum_normalize(&tmp);
+
+      quantum_delete_qureg(&tmp2);
+
+      tmp2 = quantum_matrix_qureg(H, 0, &tmp, QUANTUM_RK4_NODELETE);
+
+      h11 = quantum_real(quantum_dot_product(&tmp2, &tmp));
+      h01 = quantum_dot_product(&tmp2, reg);
+
+      t = sqrt(h11*h11-2*h00*h11+4*h01*quantum_conj(h01)+h00*h00);
+
+      E1 = -(t-h11-h00)/2.;
+      E2 = (t+h11+h00)/2.;
+
+      if(E1<E2)
+	{
+	  quantum_scalar_qureg(-(t-h11+h00)/2./h01, &tmp);
+	  quantum_vectoradd_inplace(reg, &tmp);
+	}
+      else
+	{
+	  quantum_scalar_qureg((t+h11-h00)/2./h01, &tmp);
+	  quantum_vectoradd_inplace(reg, &tmp);
+	}
+      
+      quantum_delete_qureg(&tmp);
+      quantum_delete_qureg(&tmp2);
+	
+    
+    }
+
+  quantum_error(QUANTUM_ENOCONVERGE);  
+  return nan("0");
+  
+}
+
+/* Standard Lanczos algorithm without reorthogonalization (see, e.g.,
+   [E. Dagotto, Rev. Mod. Phys. 66, 763 (1994)]. */
+
+double 
+quantum_lanczos(quantum_reg H(MAX_UNSIGNED, double), double epsilon, 
+		quantum_reg *reg)
+{
+#ifdef HAVE_LIBLAPACK
+  double E0=DBL_MAX, Eold=DBL_MAX, *a, *b, *d, *e, norm, *eig, *work;
+  quantum_reg *phi, tmp;
+  int n, i, j;
+  char jobz = 'V';
+  int lwork, *iwork, liwork, info;
+
+  phi = calloc(2, sizeof(quantum_reg));
+  a = calloc(2, sizeof(double));
+  b = calloc(2, sizeof(double));
+
+  work = malloc(sizeof(double));
+  iwork = malloc(sizeof(int));
+
+  eig = malloc(sizeof(double));
+  d = malloc(sizeof(double));
+  e = malloc(sizeof(double));
+
+  if(!(phi && a && b && work && iwork && eig && d && e))
+    quantum_error(QUANTUM_ENOMEM);
+
+  quantum_memman(2*sizeof(quantum_reg)+4*sizeof(double));
+
+  quantum_copy_qureg(reg, &phi[0]);
+  quantum_normalize(&phi[0]);
+
+  tmp = quantum_matrix_qureg(H, 0, &phi[0], QUANTUM_RK4_NODELETE);
+
+  a[0] = quantum_dot_product(&tmp, &phi[0]);
+  
+  quantum_copy_qureg(&phi[0], &phi[1]);
+  quantum_scalar_qureg(-a[0], &phi[1]);
+  quantum_vectoradd_inplace(&phi[1], &tmp);
+
+  quantum_delete_qureg(&tmp);
+
+  tmp = quantum_matrix_qureg(H, 0, &phi[1], QUANTUM_RK4_NODELETE);
+  
+  norm = quantum_dot_product(&phi[1], &phi[1]);
+
+  a[1] = quantum_dot_product(&tmp, &phi[1]) / norm;
+  b[0] = norm / quantum_dot_product(&phi[0], &phi[0]);
+
+  for(n=2; n<reg->size; n++)
+    {
+      lwork = n*n+4*n+1;
+      work = realloc(work, lwork*sizeof(double));
+
+      liwork = 5*n+3;
+      iwork = realloc(iwork, lwork*sizeof(int));
+
+      eig = realloc(eig, n*n*sizeof(double));
+      d = realloc(d, n*sizeof(double));
+      e = realloc(e, n*sizeof(double));
+
+      if(!(work && iwork && eig && d && e))
+	quantum_error(QUANTUM_ENOMEM);
+
+      memcpy(d, a, n*sizeof(double));
+
+      for(i=0; i<n; i++)
+	e[i] = sqrt(b[i]);
+
+      dstevd_(&jobz, &n, d, e, eig, &n, work, &lwork, iwork, &liwork, &info);
+
+      if(info < 0)
+	quantum_error(QUANTUM_ELAPACKARG);
+
+      else if(info > 0)
+	quantum_error(QUANTUM_ELAPACKCONV);
+
+      E0 = d[0];
+
+      if(fabs(E0-Eold) < epsilon)
+	break;
+
+      Eold = E0;
+
+      phi = realloc(phi, (n+1)*sizeof(quantum_reg));
+      a = realloc(a, (n+1)*sizeof(double));
+      b = realloc(b, (n+1)*sizeof(double));
+
+      if(!(phi && a && b))
+	quantum_error(QUANTUM_ENOMEM);
+
+      quantum_memman(sizeof(quantum_reg)+2*sizeof(double));
+
+      quantum_copy_qureg(&phi[n-1], &phi[n]);
+      quantum_scalar_qureg(-a[n-1], &phi[n]);
+      quantum_vectoradd_inplace(&phi[n], &tmp);
+
+      quantum_delete_qureg(&tmp);
+
+      quantum_copy_qureg(&phi[n-2], &tmp);
+      quantum_scalar_qureg(-b[n-2], &tmp);
+      quantum_vectoradd_inplace(&phi[n], &tmp);
+
+      /*      printf("%i %f\n", n, quantum_prob(quantum_dot_product(&phi[n], 
+	      &phi[0]))); */
+
+      quantum_delete_qureg(&tmp);
+
+      tmp = quantum_matrix_qureg(H, 0, &phi[n], QUANTUM_RK4_NODELETE);  
+      
+      norm = quantum_dot_product(&phi[n], &phi[n]);
+
+      a[n] = quantum_dot_product(&tmp, &phi[n]) / norm;
+      b[n-1] = norm / quantum_dot_product(&phi[n-1], &phi[n-1]);
+
+    }
+
+  if(n == reg->size)
+    {
+      quantum_error(QUANTUM_ENOCONVERGE);  
+      return nan("0");
+    }
+
+  for(i=0; i<n; i++)
+    quantum_normalize(&phi[i]);
+
+  for(i=0; i<reg->size; i++)
+    {
+      reg->amplitude[i] = 0;
+      for(j=0; j<n; j++)
+	reg->amplitude[i] += eig[j]*phi[j].amplitude[i];
+    }
+
+  quantum_delete_qureg(&tmp);
+
+  for(i=0; i<n; i++)
+    quantum_delete_qureg(&phi[i]);
+
+  free(phi);
+  free(a);
+  free(b);
+  free(d);
+  free(e);
+  free(eig);
+  free(work);
+  free(iwork);
+
+  return E0;
+
+
+#else
+  quantum_error(QUANTUM_ENOLAPACK);
+
+#endif /* HAVE_LIBLAPACK */
+
+}
+
+/* Imaginary time evolution algorithm */
+
+double 
+quantum_imaginary_time(quantum_reg H(MAX_UNSIGNED, double),  double epsilon, 
+		       double dt, quantum_reg *reg)
+{
+  double E0=DBL_MAX, Eold=DBL_MAX;
+  quantum_reg reg2;
+  int i;
+
+  for(i=0; i<reg->size; i++)
+    {
+      quantum_rk4(reg, 0, dt, H, QUANTUM_RK4_IMAGINARY | QUANTUM_RK4_NODELETE);
+      reg2 = quantum_matrix_qureg(H, 0, reg, QUANTUM_RK4_NODELETE);
+
+      E0 =  quantum_real(quantum_dot_product(&reg2, reg));
+
+      quantum_delete_qureg(&reg2);
+
+      if(fabs(Eold-E0)<epsilon)
+      	break;
+      Eold = E0;
+    }
+
+  if(i == reg->size)
+    {
+      quantum_error(QUANTUM_ENOCONVERGE);
+      return nan("0");
+    }
+  
+  else
+    return E0;
+}
+
+/* Wrapper around the various solver functions */
+
+double 
+quantum_groundstate(quantum_reg *reg, double epsilon, 
+		    quantum_reg H(MAX_UNSIGNED, double), int solver,
+		    double stepsize)
+{
+  switch(solver)
+    {
+    case QUANTUM_SOLVER_LANCZOS:
+      return quantum_lanczos(H, epsilon, reg);
+    case QUANTUM_SOLVER_LANCZOS_MODIFIED:
+      return quantum_lanczos_modified(H, epsilon, reg);
+    case QUANTUM_SOLVER_IMAGINARY_TIME:
+      return quantum_imaginary_time(H, epsilon, stepsize, reg);
+    default:
+      quantum_error(QUANTUM_ENOSOLVER);
+      return nan("0");
+    }
+}

energy.h

@@ -0,0 +1,41 @@
+/* energy.h: Declarations for energy.c
+
+   Copyright 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
+
+*/
+
+#ifndef __ENERGY_H
+
+#define __ENERGY_H
+
+#include "config.h"
+#include "qureg.h"
+
+enum {
+  QUANTUM_SOLVER_LANCZOS,
+  QUANTUM_SOLVER_LANCZOS_MODIFIED,
+  QUANTUM_SOLVER_IMAGINARY_TIME
+};
+
+extern double quantum_groundstate(quantum_reg *reg, double epsilon, 
+				  quantum_reg H(MAX_UNSIGNED, double), 
+				  int solver, double stepsize);
+
+#endif

error.c

@@ -54,12 +54,16 @@ quantum_strerr(int errno)
       return "wrong matrix size";
     case QUANTUM_EHASHFULL:
       return "hash table full";
+    case QUANTUM_EHERMITIAN:
+      return "matrix not Hermitian";
+    case QUANTUM_ENOCONVERGE:
+      return "method failed to converge";
     case QUANTUM_ENOLAPACK:
       return "LAPACK support not compiled in";
     case QUANTUM_ELAPACKARG:
       return "wrong arguments supplied to LAPACK";
-    case QUANTUM_ELAPACKCHEEV:
+    case QUANTUM_ELAPACKCONV:
-      return "LAPACK's CHEEV failed to converge";
+      return "LAPACK failed to converge";
     case QUANTUM_EMCMATRIX:
       return "single-column matrix expected";
     case QUANTUM_EOPCODE:

error.h

@@ -32,9 +32,12 @@ enum {
   QUANTUM_EMLARGE      = 3,
   QUANTUM_EMSIZE       = 4,
   QUANTUM_EHASHFULL    = 5,
+  QUANTUM_EHERMITIAN   = 6, 
+  QUANTUM_ENOCONVERGE  = 7,
+  QUANTUM_ENOSOLVER    = 8,
   QUANTUM_ENOLAPACK    = 32768, /* LAPACK errors start at 32768 */
   QUANTUM_ELAPACKARG   = 32769,
-  QUANTUM_ELAPACKCHEEV = 32770,
+  QUANTUM_ELAPACKCONV  = 32770,
   QUANTUM_EMCMATRIX    = 65536, /* internal errors start at 65536 */
   QUANTUM_EOPCODE      = 65537
 };

gates.c

@@ -52,12 +52,15 @@ quantum_cnot(int control, int target, quantum_reg *reg)
       if(quantum_objcode_put(CNOT, control, target))
 	return;
 
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif      
       for(i=0; i<reg->size; i++)
 	{
 	  /* Flip the target bit of a basis state if the control bit is set */
       
-	  if((reg->node[i].state & ((MAX_UNSIGNED) 1 << control)))
+	  if((reg->state[i] & ((MAX_UNSIGNED) 1 << control)))
-	    reg->node[i].state ^= ((MAX_UNSIGNED) 1 << target);
+	    reg->state[i] ^= ((MAX_UNSIGNED) 1 << target);
 	}
       quantum_decohere(reg);
     }
@@ -80,16 +83,19 @@ quantum_toffoli(int control1, int control2, int target, quantum_reg *reg)
       if(quantum_objcode_put(TOFFOLI, control1, control2, target))
 	return;
 
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
       for(i=0; i<reg->size; i++)
 	{
 	  /* Flip the target bit of a basis state if both control bits are
 	     set */
 
-	  if(reg->node[i].state & ((MAX_UNSIGNED) 1 << control1))
+	  if(reg->state[i] & ((MAX_UNSIGNED) 1 << control1))
 	    {
-	      if(reg->node[i].state & ((MAX_UNSIGNED) 1 << control2))
+	      if(reg->state[i] & ((MAX_UNSIGNED) 1 << control2))
 		{
-		  reg->node[i].state ^= ((MAX_UNSIGNED) 1 << target);
+		  reg->state[i] ^= ((MAX_UNSIGNED) 1 << target);
 		}
 	    }
 	}
@@ -126,13 +132,16 @@ quantum_unbounded_toffoli(int controlling, quantum_reg *reg, ...)
 
   va_end(bits);
 
+#ifdef _OPENMP
+#pragma omp parallel for private (j)
+#endif      
   for(i=0; i<reg->size; i++)
     {
       for(j=0; (j < controlling) && 
-	    (reg->node[i].state & (MAX_UNSIGNED) 1 << controls[j]); j++);
+	    (reg->state[i] & (MAX_UNSIGNED) 1 << controls[j]); j++);
       
       if(j == controlling) /* all control bits are set */
-	reg->node[i].state ^= ((MAX_UNSIGNED) 1 << target);
+	reg->state[i] ^= ((MAX_UNSIGNED) 1 << target);
     }
 
   free(controls);
@@ -160,11 +169,14 @@ quantum_sigma_x(int target, quantum_reg *reg)
       if(quantum_objcode_put(SIGMA_X, target))
 	return;
 
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif      
       for(i=0; i<reg->size; i++)
 	{
 	  /* Flip the target bit of each basis state */
 
-	  reg->node[i].state ^= ((MAX_UNSIGNED) 1 << target);
+	  reg->state[i] ^= ((MAX_UNSIGNED) 1 << target);
 	} 
       quantum_decohere(reg);
     }
@@ -180,17 +192,20 @@ quantum_sigma_y(int target, quantum_reg *reg)
   if(quantum_objcode_put(SIGMA_Y, target))
     return;
 
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif        
   for(i=0; i<reg->size;i++)
     {
       /* Flip the target bit of each basis state and multiply with 
 	 +/- i */
 
-      reg->node[i].state ^= ((MAX_UNSIGNED) 1 << target);
+      reg->state[i] ^= ((MAX_UNSIGNED) 1 << target);
       
-      if(reg->node[i].state & ((MAX_UNSIGNED) 1 << target))
+      if(reg->state[i] & ((MAX_UNSIGNED) 1 << target))
-	reg->node[i].amplitude *= IMAGINARY;
+	reg->amplitude[i] *= IMAGINARY;
       else
-	reg->node[i].amplitude *= -IMAGINARY;
+	reg->amplitude[i] *= -IMAGINARY;
     }
 
   quantum_decohere(reg);
@@ -206,12 +221,15 @@ quantum_sigma_z(int target, quantum_reg *reg)
   if(quantum_objcode_put(SIGMA_Z, target))
     return;
 
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif      
   for(i=0; i<reg->size; i++)
     {
       /* Multiply with -1 if the target bit is set */
 
-      if(reg->node[i].state & ((MAX_UNSIGNED) 1 << target))
+      if(reg->state[i] & ((MAX_UNSIGNED) 1 << target))
-	reg->node[i].amplitude *= -1;
+	reg->amplitude[i] *= -1;
     }
   quantum_decohere(reg);
 }
@@ -248,21 +266,21 @@ quantum_swaptheleads(int width, quantum_reg *reg)
 
 	  /* calculate left bit pattern */
 	  
-	  pat1 = reg->node[i].state % ((MAX_UNSIGNED) 1 << width);
+	  pat1 = reg->state[i] % ((MAX_UNSIGNED) 1 << width);
 	  
 	  /*calculate right but pattern */
 	  
 	  pat2 = 0;
 
 	  for(j=0; j<width; j++)
-	    pat2 += reg->node[i].state & ((MAX_UNSIGNED) 1 << (width + j));
+	    pat2 += reg->state[i] & ((MAX_UNSIGNED) 1 << (width + j));
 	  
 	  /* construct the new basis state */
 	  
-	  l = reg->node[i].state - (pat1 + pat2);
+	  l = reg->state[i] - (pat1 + pat2);
 	  l += (pat1 << width);
 	  l += (pat2 >> width);
-	  reg->node[i].state = l;
+	  reg->state[i] = l;
 	}
     }
 }
@@ -307,26 +325,29 @@ quantum_gate1(int target, quantum_matrix m, quantum_reg *reg)
 	{
 	  /* determine whether XORed basis state already exists */
 
-	  if(quantum_get_state(reg->node[i].state 
+	  if(quantum_get_state(reg->state[i] 
 			       ^ ((MAX_UNSIGNED) 1 << target), *reg) == -1)
 	    addsize++;
 	}
       
       /* allocate memory for the new basis states */
   
-      reg->node = realloc(reg->node, 
+      reg->state = realloc(reg->state, 
-			  (reg->size + addsize) * sizeof(quantum_reg_node));
+			       (reg->size + addsize) * sizeof(MAX_UNSIGNED));
+      reg->amplitude = realloc(reg->amplitude, 
+			       (reg->size + addsize) * sizeof(COMPLEX_FLOAT));
       
-      if(!reg->node) 
+      if(reg->size && !(reg->state && reg->amplitude)) 
 	quantum_error(QUANTUM_ENOMEM);
 
-      quantum_memman(addsize*sizeof(quantum_reg_node));
+      quantum_memman(addsize*(sizeof(COMPLEX_FLOAT) + sizeof(MAX_UNSIGNED)));
       
       for(i=0; i<addsize; i++)
 	{
-	  reg->node[i+reg->size].state = 0;
+	  reg->state[i+reg->size] = 0;
-	  reg->node[i+reg->size].amplitude = 0;
+	  reg->amplitude[i+reg->size] = 0;
 	}
+       
     }
 
   done = calloc(reg->size + addsize, sizeof(char));
@@ -348,29 +369,29 @@ quantum_gate1(int target, quantum_matrix m, quantum_reg *reg)
 	{
 	  /* determine if the target of the basis state is set */
 	  
-	  iset = reg->node[i].state & ((MAX_UNSIGNED) 1 << target);
+	  iset = reg->state[i] & ((MAX_UNSIGNED) 1 << target);
 
 	  tnot = 0;
-	  j = quantum_get_state(reg->node[i].state 
+	  j = quantum_get_state(reg->state[i] 
 				^ ((MAX_UNSIGNED) 1<<target), *reg);
-	  t = reg->node[i].amplitude;
-
 	  if(j >= 0)
-	    tnot = reg->node[j].amplitude;
+	    tnot = reg->amplitude[j];
+
+	  t = reg->amplitude[i];
 
 	  if(iset)
-	    reg->node[i].amplitude = m.t[2] * tnot + m.t[3] * t;
+	    reg->amplitude[i] = m.t[2] * tnot + m.t[3] * t;
 
 	  else
-	    reg->node[i].amplitude = m.t[0] * t + m.t[1] * tnot;
+	    reg->amplitude[i] = m.t[0] * t + m.t[1] * tnot;
 
 	  if(j >= 0)
 	    {
 	      if(iset)
-		reg->node[j].amplitude = m.t[0] * tnot + m.t[1] * t;
+		reg->amplitude[j] = m.t[0] * tnot + m.t[1] * t;
 
 	      else
-		reg->node[j].amplitude = m.t[2] * t + m.t[3] * tnot;
+		reg->amplitude[j] = m.t[2] * t + m.t[3] * tnot;
 	    }
 
 	  
@@ -382,14 +403,14 @@ quantum_gate1(int target, quantum_matrix m, quantum_reg *reg)
 	      if((m.t[2] == 0) && !(iset))
 		 break; 
 
-	      reg->node[k].state = reg->node[i].state 
+	      reg->state[k] = reg->state[i] 
 		^ ((MAX_UNSIGNED) 1 << target);
 
 	      if(iset)
-		reg->node[k].amplitude = m.t[1] * t;
+		reg->amplitude[k] = m.t[1] * t;
 
 	      else
-		reg->node[k].amplitude = m.t[2] * t;
+		reg->amplitude[k] = m.t[2] * t;
 
 	      k++;
 	    }
@@ -411,7 +432,7 @@ quantum_gate1(int target, quantum_matrix m, quantum_reg *reg)
     {
       for(i=0, j=0; i<reg->size; i++)
 	{
-	  if(quantum_prob_inline(reg->node[i].amplitude) < limit)
+	  if(quantum_prob_inline(reg->amplitude[i]) < limit)
 	    {
 	      j++;
 	      decsize++;
@@ -419,23 +440,33 @@ quantum_gate1(int target, quantum_matrix m, quantum_reg *reg)
 	  
 	  else if(j)
 	    {
-	      reg->node[i-j].state = reg->node[i].state;
+	      reg->state[i-j] = reg->state[i];
-	      reg->node[i-j].amplitude = reg->node[i].amplitude;
+	      reg->amplitude[i-j] = reg->amplitude[i];
 	    }
 	}
     
       if(decsize)
 	{
 	  reg->size -= decsize;
-	  reg->node = realloc(reg->node, reg->size * sizeof(quantum_reg_node));
+	  reg->amplitude = realloc(reg->amplitude, 
+				   reg->size * sizeof(COMPLEX_FLOAT));
+	  reg->state = realloc(reg->state, 
+			       reg->size * sizeof(MAX_UNSIGNED));
 	  
-	  if(!reg->node) 
+	  
+	  if(reg->size && !(reg->state && reg->amplitude)) 
 	    quantum_error(QUANTUM_ENOMEM);
 
-	  quantum_memman(-decsize * sizeof(quantum_reg_node));
+	  quantum_memman(-decsize * (sizeof(MAX_UNSIGNED) 
+				     + sizeof(COMPLEX_FLOAT)));
 	}
     }
 
+  if(reg->size > (1 << (reg->hashw-1)))
+    fprintf(stderr, "Warning: inefficient hash table (size %i vs hash %i)\n", 
+	    reg->size, 1<<reg->hashw);
+
+
   quantum_decohere(reg);
 }
 
@@ -464,34 +495,36 @@ quantum_gate2(int target1, int target2, quantum_matrix m, quantum_reg *reg)
     reg->hash[i] = 0;
       
   for(i=0; i<reg->size; i++)
-    quantum_add_hash(reg->node[i].state, i, reg);
+    quantum_add_hash(reg->state[i], i, reg);
 
   /* calculate the number of basis states to be added */
 
   for(i=0; i<reg->size; i++)
     {
-      if(quantum_get_state(reg->node[i].state ^ ((MAX_UNSIGNED) 1 << target1),
+      if(quantum_get_state(reg->state[i] ^ ((MAX_UNSIGNED) 1 << target1),
 			   *reg) == -1)
 	addsize++;
-      if(quantum_get_state(reg->node[i].state ^ ((MAX_UNSIGNED) 1 << target2),
+      if(quantum_get_state(reg->state[i] ^ ((MAX_UNSIGNED) 1 << target2),
 			   *reg) == -1)
 	addsize++;
     }
 
   /* allocate memory for the new basis states */
 
-  reg->node = realloc(reg->node, 
+  reg->state = realloc(reg->state, 
-		      (reg->size + addsize) * sizeof(quantum_reg_node));
+		       (reg->size + addsize) * sizeof(MAX_UNSIGNED));
+  reg->amplitude = realloc(reg->amplitude, 
+			   (reg->size + addsize) * sizeof(COMPLEX_FLOAT));
       
-  if(!reg->node) 
+  if(reg->size && !(reg->state && reg->amplitude)) 
-    quantum_error(QUANTUM_EMSIZE);
+    quantum_error(QUANTUM_ENOMEM);
 
-  quantum_memman(addsize*sizeof(quantum_reg_node));
+  quantum_memman(addsize*(sizeof(COMPLEX_FLOAT) + sizeof(MAX_UNSIGNED)));
 
   for(i=0; i<addsize; i++)
     {
-      reg->node[i+reg->size].state = 0;
+      reg->state[i+reg->size] = 0;
-      reg->node[i+reg->size].amplitude = 0;
+      reg->amplitude[i+reg->size] = 0;
     }
 
   done = calloc(reg->size + addsize, sizeof(char));
@@ -514,15 +547,15 @@ quantum_gate2(int target1, int target2, quantum_matrix m, quantum_reg *reg)
     {
       if(!done[i])
 	{
-	  j = quantum_bitmask(reg->node[i].state, 2, bits);
+	  j = quantum_bitmask(reg->state[i], 2, bits);
 	  base[j] = i;
-	  base[j ^ 1] = quantum_get_state(reg->node[i].state 
+	  base[j ^ 1] = quantum_get_state(reg->state[i] 
 					  ^ ((MAX_UNSIGNED) 1 << target2),
 					  *reg);
-	  base[j ^ 2] = quantum_get_state(reg->node[i].state
+	  base[j ^ 2] = quantum_get_state(reg->state[i]
 					  ^ ((MAX_UNSIGNED) 1 << target1), 
 					  *reg);
-	  base[j ^ 3] = quantum_get_state(reg->node[i].state
+	  base[j ^ 3] = quantum_get_state(reg->state[i]
 					  ^ ((MAX_UNSIGNED) 1 << target1)
 					  ^ ((MAX_UNSIGNED) 1 << target2),
 					  *reg);
@@ -532,17 +565,17 @@ quantum_gate2(int target1, int target2, quantum_matrix m, quantum_reg *reg)
 	      if(base[j] == -1)
 		{
 		  base[j] = l;
-		  //		  reg->node[l].state = reg->node[i].state
+		  //		  reg->node[l].state = reg->state[i]
 		  l++;
 		}
-	      psi_sub[j] = reg->node[base[j]].amplitude;
+	      psi_sub[j] = reg->amplitude[base[j]];
 	    }
 
 	  for(j=0; j<4; j++)
 	    {
-	      reg->node[base[j]].amplitude = 0;
+	      reg->amplitude[base[j]] = 0;
 	      for(k=0; k<4; k++)
-		reg->node[base[j]].amplitude += M(m, k, j) * psi_sub[k];
+		reg->amplitude[base[j]] += M(m, k, j) * psi_sub[k];
 
 	      done[base[j]] = 1;
 	    }
@@ -560,7 +593,7 @@ quantum_gate2(int target1, int target2, quantum_matrix m, quantum_reg *reg)
 
   for(i=0, j=0; i<reg->size; i++)
     {
-      if(quantum_prob_inline(reg->node[i].amplitude) < limit)
+      if(quantum_prob_inline(reg->amplitude[i]) < limit)
 	{
 	  j++;
 	  decsize++;
@@ -568,20 +601,26 @@ quantum_gate2(int target1, int target2, quantum_matrix m, quantum_reg *reg)
       
       else if(j)
 	{
-	  reg->node[i-j].state = reg->node[i].state;
+	  reg->state[i-j] = reg->state[i];
-	  reg->node[i-j].amplitude = reg->node[i].amplitude;
+	  reg->amplitude[i-j] = reg->amplitude[i];
 	}
     }
 
   if(decsize)
     {
       reg->size -= decsize;
-      reg->node = realloc(reg->node, reg->size * sizeof(quantum_reg_node));
+      reg->amplitude = realloc(reg->amplitude, 
+			       reg->size * sizeof(COMPLEX_FLOAT));
+      reg->state = realloc(reg->state, 
+			   reg->size * sizeof(MAX_UNSIGNED));
 	  
-      if(!reg->node) 
+	  
+      if(reg->size && !(reg->state && reg->amplitude)) 
 	quantum_error(QUANTUM_ENOMEM);
 
-      quantum_memman(-decsize * sizeof(quantum_reg_node));
+      quantum_memman(-decsize * (sizeof(MAX_UNSIGNED) 
+				 + sizeof(COMPLEX_FLOAT)));
+      
     }
 
   quantum_decohere(reg);
@@ -677,10 +716,10 @@ quantum_r_z(int target, float gamma, quantum_reg *reg)
   
   for(i=0; i<reg->size; i++)
     {
-      if(reg->node[i].state & ((MAX_UNSIGNED) 1 << target))
+      if(reg->state[i] & ((MAX_UNSIGNED) 1 << target))
-	reg->node[i].amplitude *= z;
+	reg->amplitude[i] *= z;
       else
-	reg->node[i].amplitude /= z;
+	reg->amplitude[i] /= z;
     }
 
   quantum_decohere(reg);
@@ -699,9 +738,12 @@ quantum_phase_scale(int target, float gamma, quantum_reg *reg)
 
   z = quantum_cexp(gamma);
 
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif        
   for(i=0; i<reg->size; i++)
     {
-      reg->node[i].amplitude *= z;
+      reg->amplitude[i] *= z;
     }
 
   quantum_decohere(reg);
@@ -721,10 +763,13 @@ quantum_phase_kick(int target, float gamma, quantum_reg *reg)
 
   z = quantum_cexp(gamma);
 
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif        
   for(i=0; i<reg->size; i++)
     {
-      if(reg->node[i].state & ((MAX_UNSIGNED) 1 << target))
+      if(reg->state[i] & ((MAX_UNSIGNED) 1 << target))
-	reg->node[i].amplitude *= z;
+	reg->amplitude[i] *= z;
     }
 
   quantum_decohere(reg);
@@ -743,12 +788,15 @@ quantum_cond_phase(int control, int target, quantum_reg *reg)
 
   z = quantum_cexp(pi / ((MAX_UNSIGNED) 1 << (control - target)));
 
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif      
   for(i=0; i<reg->size; i++)
     {
-      if(reg->node[i].state & ((MAX_UNSIGNED) 1 << control))
+      if(reg->state[i] & ((MAX_UNSIGNED) 1 << control))
 	{
-	  if(reg->node[i].state & ((MAX_UNSIGNED) 1 << target))
+	  if(reg->state[i] & ((MAX_UNSIGNED) 1 << target))
-	    reg->node[i].amplitude *= z;
+	    reg->amplitude[i] *= z;
 	}
     }
 
@@ -764,12 +812,15 @@ quantum_cond_phase_inv(int control, int target, quantum_reg *reg)
 
   z = quantum_cexp(-pi / ((MAX_UNSIGNED) 1 << (control - target)));
 
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif      
   for(i=0; i<reg->size; i++)
     {
-      if(reg->node[i].state & ((MAX_UNSIGNED) 1 << control))
+      if(reg->state[i] & ((MAX_UNSIGNED) 1 << control))
 	{
-	  if(reg->node[i].state & ((MAX_UNSIGNED) 1 << target))
+	  if(reg->state[i] & ((MAX_UNSIGNED) 1 << target))
-	    reg->node[i].amplitude *= z;
+	    reg->amplitude[i] *= z;
 	}
     }
 
@@ -788,17 +839,49 @@ quantum_cond_phase_kick(int control, int target, float gamma, quantum_reg *reg)
 
   z = quantum_cexp(gamma);
 
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif      
   for(i=0; i<reg->size; i++)
     {
-      if(reg->node[i].state & ((MAX_UNSIGNED) 1 << control))
+      if(reg->state[i] & ((MAX_UNSIGNED) 1 << control))
 	{
-	  if(reg->node[i].state & ((MAX_UNSIGNED) 1 << target))
+	  if(reg->state[i] & ((MAX_UNSIGNED) 1 << target))
-	    reg->node[i].amplitude *= z;
+	    reg->amplitude[i] *= z;
 	}
      }
   quantum_decohere(reg);
 }
 
+void
+quantum_cond_phase_shift(int control, int target, float gamma, quantum_reg *reg)
+{
+  int i;
+  COMPLEX_FLOAT z;
+
+  if(quantum_objcode_put(COND_PHASE, control, target, (double) gamma))
+    return;  
+
+  z = quantum_cexp(gamma/2);
+
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif      
+  for(i=0; i<reg->size; i++)
+    {
+      if(reg->state[i] & ((MAX_UNSIGNED) 1 << control))
+	{
+	  if(reg->state[i] & ((MAX_UNSIGNED) 1 << target))
+	    reg->amplitude[i] *= z;
+	  else
+	    reg->amplitude[i] /= z;
+	}
+     }
+  quantum_decohere(reg);
+}
+
+
+
 
 /* Increase the gate counter by INC steps or reset it if INC < 0. The
    current value of the counter is returned. */

grover.c

@@ -145,6 +145,7 @@ int main(int argc, char **argv)
     width = quantum_getwidth(N+1);
 
   reg = quantum_new_qureg(0, width);
+  //  reg.width--;
 
   quantum_sigma_x(reg.width, &reg);
 
@@ -169,9 +170,9 @@ int main(int argc, char **argv)
   
   for(i=0; i<reg.size; i++)
     {
-      if(reg.node[i].state == N)
+      if(reg.state[i] == N)
 	printf("\nFound %i with a probability of %f\n\n", N, 
-	       quantum_prob(reg.node[i].amplitude));
+	       quantum_prob(reg.amplitude[i]));
     }
 
   quantum_delete_qureg(&reg);

ising.c

@@ -0,0 +1,182 @@
+/* ising.c: Calculate the ground state of the transverse field Ising model
+
+   Copyright 2013 Bjoern Butscher, 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 <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+#include <math.h>
+
+#include <quantum.h>
+
+quantum_reg *hreg;
+int N;
+double g;
+int *V;
+
+quantum_reg H(MAX_UNSIGNED i, double t)
+{
+  quantum_reg reg;
+  int j;
+
+  reg = quantum_new_qureg_sparse(N+1, N);
+
+  /* Transverse field part */
+
+  for(j=0; j<N; j++)
+    {
+      reg.state[j] = i^(1 << j);
+      reg.amplitude[j] = g;
+    }
+
+  reg.state[N] = i;
+
+  /* Interaction part */
+
+  reg.amplitude[N] = V[i];
+
+  return reg;
+}
+
+quantum_reg H2(MAX_UNSIGNED i, double t)
+{
+  return hreg[i];
+}
+
+int main()
+{
+  quantum_reg reg;
+  int i, j, k;
+  double E0, m, m2;
+
+  printf("# Ground state properties of the transverse Ising chain\n");
+  printf("# g: Transverse field in units of the Ising interaction\n");
+  printf("# N: Number of spins\n");
+  printf("# E_0: Ground state energy\n");
+  printf("# m: Spontaneous magnetization\n");
+  printf("# x: Spin susceptibility\n");
+  printf("# g\t\tN\tE_0\t\tm\t\tx\n");
+
+  for(N=8; N<=18; N+=2)
+    {
+
+      /* Precompute interaction energies */
+
+      V = calloc(1<<N, sizeof(int));
+
+      assert(V);
+
+      for(i=0; i<(1<<N); i++)
+	{
+	  k = 0;
+	  for(j=0; j<N-1; j++)
+	    {
+	      if(i & (1<<j))
+		{
+		  if(i & (1<<(j+1)))
+		    k--;
+		  else
+		    k++;
+		}
+	      else
+		{
+		  if(i & (1<<(j+1)))
+		    k++;
+		  else
+		    k--;
+		}
+	    }
+
+	  /* Periodic boundary conditions */
+
+	  if(i & (1<<(N-1)))
+	    {
+	      if(i & 1)
+		k--;
+		  else
+		    k++;
+	    }
+	    else
+	      {
+		if(i & 1)
+		  k++;
+		else
+		  k--;
+	      }
+	    
+	  V[i] = k;
+	}
+
+      for(g=0.9; g<1.1; g+=0.01)
+	{
+
+	  reg = quantum_new_qureg_size(1<<N, N);
+
+	  for(i=0; i<(1<<N); i++)
+	    reg.amplitude[i] = rand();
+
+	  hreg = calloc(1<<N, sizeof(quantum_reg));
+
+	  assert(hreg);
+
+	  for(i=0; i<(1<<N); i++)
+	      hreg[i] = H(i, 0);
+
+	  E0 = quantum_groundstate(&reg, 1e-12, H2, QUANTUM_SOLVER_LANCZOS, 0);
+
+	  m = 0;
+	  m2 = 0;
+
+	  for(i=0; i<(1<<N); i++)
+	    {
+	      k = 0;
+	      for(j=0; j<N; j++)
+		{
+		  if(i & (1<<j))
+		    k--;
+		  else
+		    k++;
+		}
+	      m += quantum_prob(reg.amplitude[i])*abs(k);
+	      m2 += quantum_prob(reg.amplitude[i])*k*k;
+	    }
+
+	  m /= N;
+	  m2 /= N;
+
+	  printf("%f\t%i\t%f\t%f\t%f\n", g, N, E0, m, m2-m*m);
+
+	  quantum_delete_qureg(&reg);
+
+	  for(i=0; i<(1<<N); i++)
+	    quantum_delete_qureg(&hreg[i]);
+
+	  free(hreg);
+	}
+
+      free(V);
+    
+    }
+
+  return 0;
+
+}

lapack.c

@@ -1,6 +1,6 @@
 /* lapack.c: LAPACK interface
 
-   Copyright 2008 Bjoern Butscher, Hendrik Weimer
+   Copyright 2008-2013 Hendrik Weimer
 
    This file is part of libquantum
 
@@ -22,6 +22,7 @@
 */
 
 #include <stdlib.h>
+#include <math.h>
 
 #include "lapack.h"
 #include "matrix.h"
@@ -34,10 +35,14 @@ extern void cheev_(char *jobz, char *uplo, int *n, float _Complex *A, int *lda,
 		   float *w, float _Complex *work, int *lwork, float *rwork,
 		   int *info);
 
+extern void zheev_(char *jobz, char *uplo, int *n, double _Complex *A, int *lda,
+		   double *w, double _Complex *work, int *lwork, double *rwork,
+		   int *info);
+
 void 
-quantum_diag_time(float t, quantum_reg *reg0, quantum_reg *regt, 
+quantum_diag_time(double t, quantum_reg *reg0, quantum_reg *regt, 
 		  quantum_reg *tmp1, quantum_reg *tmp2, quantum_matrix H, 
-		  float **w)
+		  REAL_FLOAT **w)
 {
 #ifdef HAVE_LIBLAPACK
   char jobz = 'V';
@@ -45,32 +50,39 @@ quantum_diag_time(float t, quantum_reg *reg0, quantum_reg *regt,
   int dim = H.cols;
   COMPLEX_FLOAT *work;
   int lwork = -1;
-  float rwork[3*dim-2];
+  REAL_FLOAT rwork[3*dim-2];
   int info;
-  int i;
+  int i, j;
   void *p;
   
   if(tmp2->size != reg0->size)
     {
       /* perform diagonalization */
 
-      p = regt->node;
+      for(i=0; i<dim; i++)
+	{
+	  for(j=0; j<dim; j++)
+	    {
+	      if(sqrt(quantum_prob(M(H, i, j) - quantum_conj(M(H, j, i)))) 
+		 > 1e-6)
+		quantum_error(QUANTUM_EHERMITIAN);
+	    }
+	}
+
+      p = regt->amplitude;
       *regt = *reg0;
-      regt->node = realloc(p, regt->size*sizeof(quantum_reg_node));
+      regt->amplitude = realloc(p, regt->size*sizeof(COMPLEX_FLOAT));
-      for(i=0; i<reg0->size; i++)
-	regt->node[i].state = i;
       
-      p = tmp1->node;
+      p = tmp1->amplitude;
       *tmp1 = *reg0;
-      tmp1->node = realloc(p, regt->size*sizeof(quantum_reg_node));
+      tmp1->amplitude = realloc(p, regt->size*sizeof(COMPLEX_FLOAT));
-      for(i=0; i<reg0->size; i++)
-	tmp1->node[i].state = i;
 
-      p = tmp2->node;
+      p = tmp2->amplitude;
       *tmp2 = *reg0;
-      tmp2->node = realloc(p, regt->size*sizeof(quantum_reg_node));
+      tmp2->amplitude = realloc(p, regt->size*sizeof(COMPLEX_FLOAT));
-      for(i=0; i<reg0->size; i++)
+
-	tmp2->node[i].state = i;
+      if(!(regt->amplitude && tmp1->amplitude && tmp2->amplitude))
+	quantum_error(QUANTUM_ENOMEM);
 
       *w = malloc(dim*sizeof(float));
 
@@ -82,13 +94,14 @@ quantum_diag_time(float t, quantum_reg *reg0, quantum_reg *regt,
       if(!work)
 	quantum_error(QUANTUM_ENOMEM);
 
-      cheev_(&jobz, &uplo, &dim, H.t, &dim, *w, work, &lwork, rwork, &info);
+      QUANTUM_LAPACK_SOLVER(&jobz, &uplo, &dim, H.t, &dim, *w, work, &lwork, 
+			    rwork, &info);
 
       if(info < 0)
 	quantum_error(QUANTUM_ELAPACKARG);
 
       else if(info > 0)
-	quantum_error(QUANTUM_ELAPACKCHEEV);
+	quantum_error(QUANTUM_ELAPACKCONV);
       
       lwork = (int) work[0];
       work = realloc(work, lwork*sizeof(COMPLEX_FLOAT));
@@ -96,13 +109,14 @@ quantum_diag_time(float t, quantum_reg *reg0, quantum_reg *regt,
       if(!work)
 	quantum_error(QUANTUM_ENOMEM);
 
-      cheev_(&jobz, &uplo, &dim, H.t, &dim, *w, work, &lwork, rwork, &info);
+      QUANTUM_LAPACK_SOLVER(&jobz, &uplo, &dim, H.t, &dim, *w, work, &lwork, 
+			    rwork, &info);
 
       if(info < 0)
 	quantum_error(QUANTUM_ELAPACKARG);
 
       else if(info > 0)
-	quantum_error(QUANTUM_ELAPACKCHEEV);
+	quantum_error(QUANTUM_ELAPACKCONV);
       
       free(work);
 
@@ -113,17 +127,13 @@ quantum_diag_time(float t, quantum_reg *reg0, quantum_reg *regt,
 
   if(tmp1->size != reg0->size)
     {
-      p = regt->node;
+      p = regt->amplitude;
       *regt = *reg0;
-      regt->node = realloc(p, regt->size*sizeof(quantum_reg_node));
+      regt->amplitude = realloc(p, regt->size*sizeof(COMPLEX_FLOAT));
-      for(i=0; i<reg0->size; i++)
-	regt->node[i].state = i;
 
-      p = tmp1->node;
+      p = tmp1->amplitude;
       *tmp1 = *reg0;
-      tmp1->node = realloc(p, regt->size*sizeof(quantum_reg_node));
+      tmp1->amplitude = realloc(p, regt->size*sizeof(COMPLEX_FLOAT));
-      for(i=0; i<reg0->size; i++)
-	tmp1->node[i].state = i;
 
       quantum_adjoint(&H);
       
@@ -133,7 +143,7 @@ quantum_diag_time(float t, quantum_reg *reg0, quantum_reg *regt,
     }
 
   for(i=0; i<dim; i++)
-    tmp2->node[i].amplitude = quantum_cexp(-(*w)[i]*t)*tmp1->node[i].amplitude;
+    tmp2->amplitude[i] = quantum_cexp(-(*w)[i]*t)*tmp1->amplitude[i];
 
   quantum_mvmult(regt, H, tmp2);
 

lapack.h

@@ -1,6 +1,6 @@
-/* time.c: Declarations for lapack.h
+/* lapack.h: Declarations for lapack.c
 
-   Copyright 2006 Bjoern Butscher, Hendrik Weimer
+   Copyright 2006-2013 Hendrik Weimer
 
    This file is part of libquantum
 
@@ -29,8 +29,14 @@
 #include "matrix.h"
 #include "qureg.h"
 
-extern void quantum_diag_time(float t, quantum_reg *reg0, quantum_reg *regt, 
+#ifdef USE_DOUBLE
+#define QUANTUM_LAPACK_SOLVER zheev_
+#else
+#define QUANTUM_LAPACK_SOLVER cheev_
+#endif
+
+extern void quantum_diag_time(double t, quantum_reg *reg0, quantum_reg *regt, 
 			      quantum_reg *tmp1, quantum_reg *tmp2, 
-			      quantum_matrix H, float **w);
+			      quantum_matrix H, REAL_FLOAT **w);
 
 #endif

matrix.c

@@ -89,8 +89,18 @@ void
 quantum_print_matrix(quantum_matrix m) 
 {
   int i, j, z=0;
+  int print_imag = 0;
   /* int l; */
 
+  for(i=0; i<m.rows; i++) 
+    {
+      for(j=0; j<m.cols; j++)
+	{
+	  if(quantum_imag(M(m, j, i))/quantum_real(M(m, j, i)) > 1e-3)
+	    print_imag = 1;
+	}
+    }
+
   while ((1 << z++) < m.rows);
   z--;
 
@@ -104,8 +114,15 @@ quantum_print_matrix(quantum_matrix m)
 	  } */
 
       for(j=0; j<m.cols; j++)
-	printf("%g %+gi ", quantum_real(M(m, j, i)), 
+	{
-	       quantum_imag(M(m, j, i)));
+	  if(print_imag)
+	    printf("%3.3f%+.3fi ", quantum_real(M(m, j, i)), 
+		   quantum_imag(M(m, j, i)));
+	  else
+	    //	    printf("%3.3f ", quantum_real(M(m, j, i)));
+	    printf("%+.1f ", quantum_real(M(m, j, i)));
+	}
+	    
       printf("\n");
     }
   printf("\n");

measure.c

@@ -64,9 +64,9 @@ quantum_measure(quantum_reg reg)
 	 given base state - r, return the base state as the
 	 result. Otherwise, continue with the next base state. */
 
-      r -= quantum_prob_inline(reg.node[i].amplitude);
+      r -= quantum_prob_inline(reg.amplitude[i]);
       if(0 >= r)
-	return reg.node[i].state;
+	return reg.state[i];
     }
 
   /* The sum of all probabilities is less than 1. Usually, the cause
@@ -100,8 +100,8 @@ quantum_bmeasure(int pos, quantum_reg *reg)
 
   for(i=0; i<reg->size; i++)
     {
-      if(!(reg->node[i].state & pos2))
+      if(!(reg->state[i] & pos2))
-	pa += quantum_prob_inline(reg->node[i].amplitude);
+	pa += quantum_prob_inline(reg->amplitude[i]);
     }
 
   /* Compare the probability for 0 with a random number and determine
@@ -141,8 +141,8 @@ quantum_bmeasure_bitpreserve(int pos, quantum_reg *reg)
 
   for(i=0; i<reg->size; i++)
     {
-      if(!(reg->node[i].state & pos2))
+      if(!(reg->state[i] & pos2))
-	pa += quantum_prob_inline(reg->node[i].amplitude);
+	pa += quantum_prob_inline(reg->amplitude[i]);
     }
 
   /* Compare the probability for 0 with a random number and determine
@@ -158,23 +158,23 @@ quantum_bmeasure_bitpreserve(int pos, quantum_reg *reg)
 
   for(i=0;i<reg->size;i++)
     {
-      if(reg->node[i].state & pos2)
+      if(reg->state[i] & pos2)
 	{
 	  if(!result)
-	    reg->node[i].amplitude = 0;
+	    reg->amplitude[i] = 0;
 	  else
 	    {
-	      d += quantum_prob_inline(reg->node[i].amplitude);
+	      d += quantum_prob_inline(reg->amplitude[i]);
 	      size++;
 	    }
 	}
       else
 	{
 	  if(result)
-	    reg->node[i].amplitude = 0;
+	    reg->amplitude[i] = 0;
 	  else
 	    {
-	      d += quantum_prob_inline(reg->node[i].amplitude);
+	      d += quantum_prob_inline(reg->amplitude[i]);
 	      size++;
 	    }
 	}
@@ -183,12 +183,14 @@ quantum_bmeasure_bitpreserve(int pos, quantum_reg *reg)
   /* Build the new quantum register */
 
   out.size = size;
-  out.node = calloc(size, sizeof(quantum_reg_node));
+  out.state = calloc(size, sizeof(MAX_UNSIGNED));
+  out.amplitude = calloc(size, sizeof(COMPLEX_FLOAT));
 
-  if(!out.node)
+  if(!(out.state && out.amplitude))
     quantum_error(QUANTUM_ENOMEM);
 
-  quantum_memman(size * sizeof(quantum_reg_node));
+  quantum_memman(size * (sizeof(MAX_UNSIGNED) + sizeof(COMPLEX_FLOAT)));
+
   out.hashw = reg->hashw;
   out.hash = reg->hash;
   out.width = reg->width;
@@ -198,10 +200,10 @@ quantum_bmeasure_bitpreserve(int pos, quantum_reg *reg)
   
   for(i=0, j=0; i<reg->size; i++)
     {
-      if(reg->node[i].amplitude)
+      if(reg->amplitude[i])
 	{
-	  out.node[j].state = reg->node[i].state;
+	  out.state[j] = reg->state[i];
-	  out.node[j].amplitude = reg->node[i].amplitude * 1 / (float) sqrt(d);
+	  out.amplitude[j] = reg->amplitude[i] * 1 / (float) sqrt(d);
 	
 	  j++;
 	}

objcode.c

@@ -336,16 +336,30 @@ quantum_objcode_run(char *file, quantum_reg *reg)
       switch(operation)
 	{
 	case INIT:
-	  fread(buf, sizeof(MAX_UNSIGNED), 1, fhd);
+	  if(!fread(buf, sizeof(MAX_UNSIGNED), 1, fhd))
+	    {
+	      quantum_error(QUANTUM_FAILURE);
+	      break;
+	    }
 	  mu = quantum_char2mu(buf);
 	  *reg = quantum_new_qureg(mu, 12);
+
 	  break;
 
 	case CNOT:
 	case COND_PHASE:
-	  fread(buf, sizeof(int), 1, fhd);
+	  if(!fread(buf, sizeof(int), 1, fhd))
+	    {
+	      quantum_error(QUANTUM_FAILURE);
+	      break;
+	    }
 	  j = quantum_char2int(buf);
-	  fread(buf, sizeof(int), 1, fhd);
+	  if(!fread(buf, sizeof(int), 1, fhd))
+	    {
+	      quantum_error(QUANTUM_FAILURE);
+	      break;
+	    }
+	  
 	  k = quantum_char2int(buf);
 	  switch(operation)
 	    {
@@ -357,11 +371,23 @@ quantum_objcode_run(char *file, quantum_reg *reg)
 	  break;
 
 	case TOFFOLI:
-	  fread(buf, sizeof(int), 1, fhd);
+	  if(!fread(buf, sizeof(int), 1, fhd))
+	    {
+	      quantum_error(QUANTUM_FAILURE);
+	      break;
+	    }
 	  j = quantum_char2int(buf);
-	  fread(buf, sizeof(int), 1, fhd);
+	  if(!fread(buf, sizeof(int), 1, fhd))
+	    {
+	      quantum_error(QUANTUM_FAILURE);
+	      break;
+	    }
 	  k = quantum_char2int(buf);
-	  fread(buf, sizeof(int), 1, fhd);
+	  if(!fread(buf, sizeof(int), 1, fhd))
+	    {
+	      quantum_error(QUANTUM_FAILURE);
+	      break;
+	    }
 	  l = quantum_char2int(buf);
 	  quantum_toffoli(j, k, l, reg);
 	  break;
@@ -373,7 +399,11 @@ quantum_objcode_run(char *file, quantum_reg *reg)
 	case BMEASURE:
 	case BMEASURE_P:
 	case SWAPLEADS:
-	  fread(buf, sizeof(int), 1, fhd);
+	  if(!fread(buf, sizeof(int), 1, fhd))
+	    {
+	      quantum_error(QUANTUM_FAILURE);
+	      break;
+	    }
 	  j = quantum_char2int(buf);
 	  switch(operation)
 	    {
@@ -399,9 +429,17 @@ quantum_objcode_run(char *file, quantum_reg *reg)
 	case ROT_Z:
 	case PHASE_KICK:
 	case PHASE_SCALE:
-	  fread(buf, sizeof(int), 1, fhd);
+	  if(!fread(buf, sizeof(int), 1, fhd))
+	    {
+	      quantum_error(QUANTUM_FAILURE);
+	      break;
+	    }
 	  j = quantum_char2int(buf);
-	  fread(buf, sizeof(double), 1, fhd);
+	  if(!fread(buf, sizeof(double), 1, fhd))
+	    {
+	      quantum_error(QUANTUM_FAILURE);
+	      break;
+	    }
 	  d = quantum_char2double(buf);
 	  switch(operation)
 	    {
@@ -419,11 +457,23 @@ quantum_objcode_run(char *file, quantum_reg *reg)
 	  break;
 
 	case CPHASE_KICK:
-	  fread(buf, sizeof(int), 1, fhd);
+	  if(!fread(buf, sizeof(int), 1, fhd))	
+	    {
+	      quantum_error(QUANTUM_FAILURE);
+	      break;
+	    }
 	  j = quantum_char2int(buf);
-	  fread(buf, sizeof(int), 1, fhd);
+	  if(!fread(buf, sizeof(int), 1, fhd))
+	    {
+	      quantum_error(QUANTUM_FAILURE);
+	      break;
+	    }
 	  k = quantum_char2int(buf);
-	  fread(buf, sizeof(double), 1, fhd);
+	  if(!fread(buf, sizeof(double), 1, fhd))
+	    {
+	      quantum_error(QUANTUM_FAILURE);
+	      break;
+	    }
 	  d = quantum_char2double(buf);
 	  quantum_cond_phase_kick(j, k, d, reg);
 	  break;

qec.c

@@ -243,34 +243,34 @@ quantum_toffoli_ft(int control1, int control2, int target, quantum_reg *reg)
       c1 = 0;
       c2 = 0;
 
-      if(reg->node[i].state & ((MAX_UNSIGNED) 1 << control1))
+      if(reg->state[i] & ((MAX_UNSIGNED) 1 << control1))
 	c1 = 1;
-      if(reg->node[i].state 
+      if(reg->state[i] 
 	 & ((MAX_UNSIGNED) 1 << (control1+width)))
 	{
 	  c1 ^= 1;
 	}
-      if(reg->node[i].state 
+      if(reg->state[i] 
 	 & ((MAX_UNSIGNED) 1 << (control1+2*width)))
 	{
 	  c1 ^= 1;
 	}
 
-      if(reg->node[i].state & ((MAX_UNSIGNED) 1 << control2))
+      if(reg->state[i] & ((MAX_UNSIGNED) 1 << control2))
 	c2 = 1;
-      if(reg->node[i].state 
+      if(reg->state[i] 
 	 & ((MAX_UNSIGNED) 1 << (control2+width)))
 	{
 	  c2 ^= 1;
 	}
-      if(reg->node[i].state 
+      if(reg->state[i] 
 	 & ((MAX_UNSIGNED) 1 << (control2+2*width)))
 	{
 	  c2 ^= 1;
 	}
 
       if(c1 == 1 && c2 == 1)
-	reg->node[i].state = reg->node[i].state ^ mask;
+	reg->state[i] = reg->state[i] ^ mask;
 
     }
 

qtime.c

@@ -1,6 +1,6 @@
 /* qtime.c: Time evolution of a quantum system
 
-   Copyright 2006,2007 Bjoern Butscher, Hendrik Weimer
+   Copyright 2006-2013 Hendrik Weimer
 
    This file is part of libquantum
 
@@ -30,7 +30,10 @@
 #include "complex.h"
 #include "config.h"
 
-/* Forth-order Runge-Kutta */
+/* 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, 
@@ -41,6 +44,7 @@ quantum_rk4(quantum_reg *reg, double t, double dt,
   int i;
   void *hash;
   int hashw;
+  COMPLEX_FLOAT step = dt;
 
   hash = reg->hash;
   reg->hash = 0;
@@ -48,36 +52,39 @@ quantum_rk4(quantum_reg *reg, double t, double dt,
   hashw = reg->hashw;
   reg->hashw = 0;
 
+  if(!(flags & QUANTUM_RK4_IMAGINARY))
+    step *= IMAGINARY;
+
   /* k1 */
-  k = quantum_matrix_qureg(H, t, reg, flags);
+  k = quantum_matrix_qureg(H, t, reg, flags & QUANTUM_RK4_NODELETE);
-  quantum_scalar_qureg(-IMAGINARY*dt/2.0, &k);
+  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);
+  k = quantum_matrix_qureg(H, t+dt/2.0, &tmp, flags & QUANTUM_RK4_NODELETE);
   quantum_delete_qureg(&tmp);
-  quantum_scalar_qureg(-IMAGINARY*dt/2.0, &k);
+  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);
+  k = quantum_matrix_qureg(H, t+dt/2.0, &tmp, flags & QUANTUM_RK4_NODELETE);
   quantum_delete_qureg(&tmp);
-  quantum_scalar_qureg(-IMAGINARY*dt, &k);
+  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);
+  k = quantum_matrix_qureg(H, t+dt, &tmp, flags & QUANTUM_RK4_NODELETE);
   quantum_delete_qureg(&tmp);
-  quantum_scalar_qureg(-IMAGINARY*dt/6.0, &k);
+  quantum_scalar_qureg(-step/6.0, &k);
   quantum_vectoradd_inplace(&out, &k);
   
   quantum_delete_qureg(&k);
@@ -85,10 +92,14 @@ quantum_rk4(quantum_reg *reg, double t, double dt,
 
   /* Normalize quantum register */
 
-  for(i=0; i<out.size; i++)
+  if(flags & QUANTUM_RK4_IMAGINARY)
-    r += quantum_prob(out.node[i].amplitude);
+    {
 
-  //  quantum_scalar_qureg(sqrt(1.0/r), &out);
+      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;
@@ -130,10 +141,8 @@ quantum_rk4a(quantum_reg *reg, double t, double *dt, double epsilon,
 
       for(i=0;i<reg->size;i++)
 	{
-	  r = 2*sqrt(quantum_prob(reg->node[i].amplitude 
+	  r = 2*sqrt(quantum_prob(reg->amplitude[i] - reg2.amplitude[i])/
-				  - reg2.node[i].amplitude)/
+		     quantum_prob(reg->amplitude[i] + reg2.amplitude[i]));
-		     quantum_prob(reg->node[i].amplitude 
-				  + reg2.node[i].amplitude));
 	  
 	  if(r > delta)
 	    delta = r;
@@ -155,8 +164,14 @@ quantum_rk4a(quantum_reg *reg, double t, double *dt, double epsilon,
 
       if(delta > epsilon)
 	{
-	  memcpy(reg->node, old.node, reg->size*sizeof(quantum_reg_node));
+	  memcpy(reg->amplitude, old.amplitude, 
-	  memcpy(reg2.node, old.node, reg->size*sizeof(quantum_reg_node));
+		 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);

qtime.h

@@ -1,6 +1,6 @@
 /* qtime.h: Declarations for qtime.c
 
-   Copyright 2006,2007 Bjoern Butscher, Hendrik Weimer
+   Copyright 2006-2013 Hendrik Weimer
 
    This file is part of libquantum
 
@@ -28,6 +28,11 @@
 #include "qureg.h"
 #include "config.h"
 
+enum {
+  QUANTUM_RK4_NODELETE  = 1,
+  QUANTUM_RK4_IMAGINARY = 2
+};
+
 extern void quantum_rk4(quantum_reg *reg, double t, double dt, 
 			quantum_reg H(MAX_UNSIGNED, double), int flags);
 extern double quantum_rk4a(quantum_reg *reg, double t, double *dt, 

quantum.h.in

@@ -1,6 +1,6 @@
 /* quantum.h: Header file for libquantum
 
-   Copyright 2003-2008 Bjoern Butscher, Hendrik Weimer
+   Copyright 2003-2013 Bjoern Butscher, Hendrik Weimer
 
    This file is part of libquantum
 
@@ -45,14 +45,6 @@ struct quantum_matrix_struct {
 
 typedef struct quantum_matrix_struct quantum_matrix;
 
-struct quantum_reg_node_struct
-{
-  COMPLEX_FLOAT amplitude; /* alpha_j */
-  MAX_UNSIGNED state;      /* j */
-};
-
-typedef struct quantum_reg_node_struct quantum_reg_node;
-
 /* The quantum register */
 
 struct quantum_reg_struct
@@ -60,7 +52,8 @@ struct quantum_reg_struct
   int width;    /* number of qubits in the qureg */
   int size;     /* number of non-zero vectors */
   int hashw;    /* width of the hash array */
-  quantum_reg_node *node;
+  COMPLEX_FLOAT *amplitude;
+  MAX_UNSIGNED *state;
   int *hash;
 };
 
@@ -75,8 +68,15 @@ struct quantum_density_op_struct
 
 typedef struct quantum_density_op_struct quantum_density_op;
 
+enum {
+  QUANTUM_SOLVER_LANCZOS,
+  QUANTUM_SOLVER_LANCZOS_MODIFIED,
+  QUANTUM_SOLVER_IMAGINARY_TIME
+};
+
 extern quantum_reg quantum_new_qureg(MAX_UNSIGNED initval, int width);
 extern quantum_reg quantum_new_qureg_size(int n, int width);
+extern quantum_reg quantum_new_qureg_sparse(int n, int width);
 extern void quantum_delete_qureg(quantum_reg *reg);
 extern void quantum_print_qureg(quantum_reg reg);
 extern void quantum_addscratch(int bits, quantum_reg *reg);
@@ -103,6 +103,8 @@ extern void quantum_cond_phase(int control, int target, quantum_reg *reg);
 extern void quantum_cond_phase_inv(int control, int target, quantum_reg *reg);
 extern void quantum_cond_phase_kick(int control, int target, float gamma, 
 				    quantum_reg *reg);
+extern void quantum_cond_phase_shift(int control, int target, float gamma, 
+				    quantum_reg *reg);
 extern int quantum_gate_counter(int inc);
 
 extern void quantum_qft(int width, quantum_reg *reg);
@@ -125,7 +127,7 @@ extern void quantum_cancel(int *a, int *b);
 extern void quantum_frac_approx(int *a, int *b, int width);
 extern int quantum_getwidth(int n);
 
-extern float quantum_prob(COMPLEX_FLOAT a);
+extern double quantum_prob(COMPLEX_FLOAT a);
 
 extern float quantum_get_decoherence();
 extern void quantum_set_decoherence(float lambda);
@@ -171,8 +173,13 @@ extern double quantum_rk4a(quantum_reg *reg, double t, double *dt,
 			   double epsilon, 
 			   quantum_reg H(MAX_UNSIGNED, double), int flags);
 
-extern void quantum_diag_time(float t, quantum_reg *reg0, quantum_reg *regt, 
+extern void quantum_diag_time(double t, quantum_reg *reg0, quantum_reg *regt, 
 			      quantum_reg *tmp1, quantum_reg *tmp2, 
-			      quantum_matrix H, float **w);
+			      quantum_matrix H, @RF_TYPE@ **w);
+
+
+extern double quantum_groundstate(quantum_reg *reg, double epsilon, 
+				  quantum_reg H(MAX_UNSIGNED, double), 
+				  int solver, double stepsize);
 
 #endif

qureg.c

@@ -1,6 +1,6 @@
 /* qureg.c: Quantum register management
 
-   Copyright 2003, 2004, 2006 Bjoern Butscher, Hendrik Weimer
+   Copyright 2003-2013 Bjoern Butscher, Hendrik Weimer
 
    This file is part of libquantum
 
@@ -59,12 +59,13 @@ quantum_matrix2qureg(quantum_matrix *m, int width)
   reg.size = size;
   reg.hashw = width + 2;
 
-  reg.node = calloc(size, sizeof(quantum_reg_node));
+  reg.amplitude = calloc(size, sizeof(COMPLEX_FLOAT));
+  reg.state = calloc(size, sizeof(MAX_UNSIGNED));
 
-  if(!reg.node)
+  if(!(reg.state && reg.amplitude))
     quantum_error(QUANTUM_ENOMEM);
 
-  quantum_memman(size * sizeof(quantum_reg_node));
+  quantum_memman(size * (sizeof(COMPLEX_FLOAT) + sizeof(MAX_UNSIGNED)));
 
   /* Allocate the hash table */
 
@@ -82,8 +83,8 @@ quantum_matrix2qureg(quantum_matrix *m, int width)
     {
       if(m->t[i])
 	{
-	  reg.node[j].state = i;
+	  reg.state[j] = i;
-	  reg.node[j].amplitude = m->t[i];
+	  reg.amplitude[j] = m->t[i];
 	  j++;
 	}
     }
@@ -105,12 +106,13 @@ quantum_new_qureg(MAX_UNSIGNED initval, int width)
 
   /* Allocate memory for 1 base state */
 
-  reg.node = calloc(1, sizeof(quantum_reg_node));
+  reg.state = calloc(1, sizeof(MAX_UNSIGNED));
+  reg.amplitude = calloc(1, sizeof(COMPLEX_FLOAT));
 
-  if(!reg.node)
+  if(!(reg.state && reg.amplitude))
     quantum_error(QUANTUM_ENOMEM);
 
-  quantum_memman(sizeof(quantum_reg_node));
+  quantum_memman(sizeof(MAX_UNSIGNED) + sizeof(COMPLEX_FLOAT));
 
   /* Allocate the hash table */
 
@@ -123,8 +125,8 @@ quantum_new_qureg(MAX_UNSIGNED initval, int width)
 
   /* Initialize the quantum register */
   
-  reg.node[0].state = initval;
+  reg.state[0] = initval;
-  reg.node[0].amplitude = 1;
+  reg.amplitude[0] = 1;
 
   /* Initialize the PRNG */
 
@@ -158,12 +160,38 @@ quantum_new_qureg_size(int n, int width)
 
   /* Allocate memory for n basis states */
 
-  reg.node = calloc(n, sizeof(quantum_reg_node));
+  reg.amplitude = calloc(n, sizeof(COMPLEX_FLOAT));
+  reg.state = 0;
 
-  if(!reg.node)
+  if(!reg.amplitude)
     quantum_error(QUANTUM_ENOMEM);
 
-  quantum_memman(n*sizeof(quantum_reg_node));
+  quantum_memman(n*sizeof(COMPLEX_FLOAT));
+
+  return reg;
+}
+
+/* Returns an empty sparse quantum register of size N */
+
+quantum_reg
+quantum_new_qureg_sparse(int n, int width)
+{
+  quantum_reg reg;
+
+  reg.width = width;
+  reg.size = n;
+  reg.hashw = 0;
+  reg.hash = 0;
+
+  /* Allocate memory for n basis states */
+
+  reg.amplitude = calloc(n, sizeof(COMPLEX_FLOAT));
+  reg.state = calloc(n, sizeof(MAX_UNSIGNED));
+
+  if(!(reg.amplitude && reg.state))
+    quantum_error(QUANTUM_ENOMEM);
+
+  quantum_memman(n*(sizeof(COMPLEX_FLOAT)+sizeof(MAX_UNSIGNED)));
 
   return reg;
 }
@@ -179,7 +207,7 @@ quantum_qureg2matrix(quantum_reg reg)
   m = quantum_new_matrix(1, 1 << reg.width);
   
   for(i=0; i<reg.size; i++)
-    m.t[reg.node[i].state] = reg.node[i].amplitude;
+    m.t[reg.state[i]] = reg.amplitude[i];
 
   return m;
 }
@@ -201,9 +229,18 @@ quantum_delete_qureg(quantum_reg *reg)
 {
   if(reg->hashw && reg->hash)
     quantum_destroy_hash(reg);
-  free(reg->node);
+
-  quantum_memman(-reg->size * sizeof(quantum_reg_node));
+  free(reg->amplitude);
-  reg->node = 0;
+  quantum_memman(-reg->size * sizeof(COMPLEX_FLOAT));
+  reg->amplitude = 0;
+
+  if(reg->state)
+    {
+      free(reg->state);
+      quantum_memman(-reg->size * sizeof(MAX_UNSIGNED));
+      reg->state = 0;
+    }
+
 }
 
 /* Delete a quantum register but leave the hash table alive */
@@ -211,9 +248,16 @@ quantum_delete_qureg(quantum_reg *reg)
 void
 quantum_delete_qureg_hashpreserve(quantum_reg *reg)
 {
-  free(reg->node);
+  free(reg->amplitude);
-  quantum_memman(-reg->size * sizeof(quantum_reg_node));
+  quantum_memman(-reg->size * sizeof(COMPLEX_FLOAT));
-  reg->node = 0;
+  reg->amplitude = 0;
+
+  if(reg->state)
+    {
+      free(reg->state);
+      quantum_memman(-reg->size * sizeof(MAX_UNSIGNED));
+      reg->state = 0;
+    }
 }
 
 /* Copy the contents of src to dst */
@@ -225,12 +269,27 @@ quantum_copy_qureg(quantum_reg *src, quantum_reg *dst)
   
   /* Allocate memory for basis states */
 
-  dst->node = calloc(dst->size, sizeof(quantum_reg_node));
+  dst->amplitude = calloc(dst->size, sizeof(COMPLEX_FLOAT));
 
-  if(!dst->node)
+  if(!dst->amplitude)
     quantum_error(QUANTUM_ENOMEM);
 
-  quantum_memman(dst->size*sizeof(quantum_reg_node));
+  quantum_memman(dst->size*sizeof(COMPLEX_FLOAT));
+
+  memcpy(dst->amplitude, src->amplitude, src->size*sizeof(COMPLEX_FLOAT));
+
+  if(src->state)
+    {
+      dst->state = calloc(dst->size, sizeof(MAX_UNSIGNED));
+
+      if(!dst->state)
+	quantum_error(QUANTUM_ENOMEM);
+      
+      quantum_memman(dst->size*sizeof(MAX_UNSIGNED));
+
+      memcpy(dst->state, src->state, src->size*sizeof(MAX_UNSIGNED));
+
+    }
 
   /* Allocate the hash table */
 
@@ -244,8 +303,6 @@ quantum_copy_qureg(quantum_reg *src, quantum_reg *dst)
       quantum_memman((1 << dst->hashw) * sizeof(int));
     }
 
-  memcpy(dst->node, src->node, src->size*sizeof(quantum_reg_node));
-
 }
 
 /* Print the contents of a quantum register to stdout */
@@ -257,14 +314,14 @@ quantum_print_qureg(quantum_reg reg)
   
   for(i=0; i<reg.size; i++)
     {
-      printf("% f %+fi|%lli> (%e) (|", quantum_real(reg.node[i].amplitude),
+      printf("% f %+fi|%lli> (%e) (|", quantum_real(reg.amplitude[i]),
-	     quantum_imag(reg.node[i].amplitude), reg.node[i].state, 
+	     quantum_imag(reg.amplitude[i]), reg.state[i], 
-	     quantum_prob_inline(reg.node[i].amplitude));
+	     quantum_prob_inline(reg.amplitude[i]));
       for(j=reg.width-1;j>=0;j--)
 	{
 	  if(j % 4 == 3)
 	    printf(" ");
-	  printf("%i", ((((MAX_UNSIGNED) 1 << j) & reg.node[i].state) > 0));
+	  printf("%i", ((((MAX_UNSIGNED) 1 << j) & reg.state[i]) > 0));
 	}
 
       printf(">)\n");
@@ -282,7 +339,7 @@ quantum_print_expn(quantum_reg reg)
   
   for(i=0; i<reg.size; i++)
     {
-      printf("%i: %lli\n", i, reg.node[i].state - i * (1 << (reg.width / 2)));
+      printf("%i: %lli\n", i, reg.state[i] - i * (1 << (reg.width / 2)));
     }
 }
 
@@ -292,17 +349,15 @@ quantum_print_expn(quantum_reg reg)
 void
 quantum_addscratch(int bits, quantum_reg *reg)
 {
-  int i, oldwidth;
+  int i;
   MAX_UNSIGNED l;
   
-  oldwidth = reg->width;
-
   reg->width += bits;
 
   for(i=0; i<reg->size; i++)
     {
-      l = reg->node[i].state << bits;
+      l = reg->state[i] << bits;
-      reg->node[i].state = l;
+      reg->state[i] = l;
     }
 }
 
@@ -318,7 +373,7 @@ quantum_print_hash(quantum_reg reg)
     {
       if(i)
 	printf("%i: %i %llu\n", i, reg.hash[i]-1, 
-	       reg.node[reg.hash[i]-1].state);
+	       reg.state[reg.hash[i]-1]);
     }
 
 }
@@ -335,15 +390,15 @@ quantum_kronecker(quantum_reg *reg1, quantum_reg *reg2)
   reg.size = reg1->size*reg2->size;
   reg.hashw = reg.width + 2;
 
-
   /* allocate memory for the new basis states */
 
-  reg.node = calloc(reg.size, sizeof(quantum_reg_node));
+  reg.amplitude = calloc(reg.size, sizeof(COMPLEX_FLOAT));
-  if(!reg.node) 
+  reg.state = calloc(reg.size, sizeof(MAX_UNSIGNED));
+
+  if(!(reg.state && reg.amplitude))
     quantum_error(QUANTUM_ENOMEM);
 
-  quantum_memman((reg.size)*sizeof(quantum_reg_node));
+  quantum_memman(reg.size * (sizeof(COMPLEX_FLOAT) + sizeof(MAX_UNSIGNED)));
-
 
   /* Allocate the hash table */
 
@@ -356,14 +411,13 @@ quantum_kronecker(quantum_reg *reg1, quantum_reg *reg2)
   for(i=0; i<reg1->size; i++)
     for(j=0; j<reg2->size; j++)
     {
-      /* printf("processing |%lli> x |%lli>\n", reg1->node[i].state, 
+      /* printf("processing |%lli> x |%lli>\n", reg1->state[i], 
-	     reg2->node[j].state);
+	     reg2->state[j]);
-         printf("%lli\n", (reg1->node[i].state) << reg2->width); */
+         printf("%lli\n", (reg1->state[i]) << reg2->width); */
 
-      reg.node[i*reg2->size+j].state = 	((reg1->node[i].state) << reg2->width)
+      reg.state[i*reg2->size+j] = ((reg1->state[i]) << reg2->width) 
-	| reg2->node[j].state;
+	| reg2->state[j];
-      reg.node[i*reg2->size+j].amplitude = 
+      reg.amplitude[i*reg2->size+j] = reg1->amplitude[i] * reg2->amplitude[j];
-	reg1->node[i].amplitude * reg2->node[j].amplitude;
     }
 
   return reg;
@@ -387,10 +441,10 @@ quantum_state_collapse(int pos, int value, quantum_reg reg)
   
   for(i=0;i<reg.size;i++)
     {
-      if(((reg.node[i].state & pos2) && value) 
+      if(((reg.state[i] & pos2) && value) 
-	 || (!(reg.node[i].state & pos2) && !value))
+	 || (!(reg.state[i] & pos2) && !value))
 	{
-	  d += quantum_prob_inline(reg.node[i].amplitude);
+	  d += quantum_prob_inline(reg.amplitude[i]);
 	  size++;
 	}
     }
@@ -399,12 +453,13 @@ quantum_state_collapse(int pos, int value, quantum_reg reg)
 
   out.width = reg.width-1;
   out.size = size;
-  out.node = calloc(size, sizeof(quantum_reg_node));
+  out.amplitude = calloc(size, sizeof(COMPLEX_FLOAT));
+  out.state = calloc(size, sizeof(MAX_UNSIGNED));
 
-  if(!out.node)
+  if(!(out.state && out.amplitude))
     quantum_error(QUANTUM_ENOMEM);
 
-  quantum_memman(size * sizeof(quantum_reg_node));
+  quantum_memman(size * (sizeof(COMPLEX_FLOAT) + sizeof(MAX_UNSIGNED)));
   out.hashw = reg.hashw;
   out.hash = reg.hash;
 
@@ -413,21 +468,21 @@ quantum_state_collapse(int pos, int value, quantum_reg reg)
 
   for(i=0, j=0; i<reg.size; i++)
     {
-      if(((reg.node[i].state & pos2) && value) 
+      if(((reg.state[i] & pos2) && value) 
-	 || (!(reg.node[i].state & pos2) && !value))
+	 || (!(reg.state[i] & pos2) && !value))
 	{
 	  for(k=0, rpat=0; k<pos; k++)
 	    rpat += (MAX_UNSIGNED) 1 << k;
 
-	  rpat &= reg.node[i].state;
+	  rpat &= reg.state[i];
 
 	  for(k=sizeof(MAX_UNSIGNED)*8-1, lpat=0; k>pos; k--)
 	    lpat += (MAX_UNSIGNED) 1 << k;
 
-	  lpat &= reg.node[i].state;
+	  lpat &= reg.state[i];
 
-	  out.node[j].state = (lpat >> 1) | rpat;
+	  out.state[j] = (lpat >> 1) | rpat;
-	  out.node[j].amplitude = reg.node[i].amplitude * 1 / (float) sqrt(d);
+	  out.amplitude[j] = reg.amplitude[i] * 1 / (float) sqrt(d);
 	
 	  j++;
 	}
@@ -450,12 +505,26 @@ quantum_dot_product(quantum_reg *reg1, quantum_reg *reg2)
   if(reg2->hashw)
     quantum_reconstruct_hash(reg2);
 
-  for(i=0; i<reg1->size; i++)
+  if(reg1->state)
     {
-      j = quantum_get_state(reg1->node[i].state, *reg2);
+      for(i=0; i<reg1->size; i++)
+	{
+	  j = quantum_get_state(reg1->state[i], *reg2);
 
-      if(j > -1) /* state exists in reg2 */
+	  if(j > -1) /* state exists in reg2 */
-	f += quantum_conj(reg1->node[i].amplitude) * reg2->node[j].amplitude;
+	    f += quantum_conj(reg1->amplitude[i]) * reg2->amplitude[j];
+	}
+    }
+
+  else
+    {
+      for(i=0; i<reg1->size; i++)
+	{
+	  j = quantum_get_state(i, *reg2);
+
+	  if(j > -1) /* state exists in reg2 */
+	    f += quantum_conj(reg1->amplitude[i]) * reg2->amplitude[j];
+	}
     }
       
   return f;
@@ -475,13 +544,21 @@ quantum_dot_product_noconj(quantum_reg *reg1, quantum_reg *reg2)
   if(reg2->hashw)
     quantum_reconstruct_hash(reg2);
 
-  for(i=0; i<reg1->size; i++)
+  if(!reg2->state)
     {
-      j = quantum_get_state(reg1->node[i].state, *reg2);
+      for(i=0; i<reg1->size; i++)
+	f += reg1->amplitude[i] * reg2->amplitude[reg1->state[i]];
+    }
 
-      if(j > -1) /* state exists in reg2 */
+  else
-	f += reg1->node[i].amplitude * reg2->node[j].amplitude;
+    {
+      for(i=0; i<reg1->size; i++)
+	{
+	  j = quantum_get_state(reg1->state[i], *reg2);
 
+	  if(j > -1) /* state exists in reg2 */
+	    f += reg1->amplitude[i] * reg2->amplitude[j];
+	}
     }
 
   return f;
@@ -510,32 +587,47 @@ quantum_vectoradd(quantum_reg *reg1, quantum_reg *reg2)
 
       for(i=0; i<reg2->size; i++)
 	{
-	  if(quantum_get_state(reg2->node[i].state, *reg1) == -1)
+	  if(quantum_get_state(reg2->state[i], *reg1) == -1)
 	    addsize++;
 	}
     }
 
-  reg.size += addsize;
+  if(addsize)
-  reg.node = realloc(reg.node, (reg.size)*sizeof(quantum_reg_node));
+    {
-  if(!reg.node)
+      reg.size += addsize;
-    quantum_error(QUANTUM_ENOMEM);
 
-  quantum_memman(addsize*sizeof(quantum_reg_node));
+      reg.amplitude = realloc(reg.amplitude, reg.size*sizeof(COMPLEX_FLOAT));
+      reg.state = realloc(reg.state, reg.size*sizeof(MAX_UNSIGNED));
+
+      if(!(reg.state && reg.amplitude))
+	quantum_error(QUANTUM_ENOMEM);
+
+      quantum_memman(addsize * (sizeof(COMPLEX_FLOAT) + sizeof(MAX_UNSIGNED)));
+    }
 
   k = reg1->size;
 
-  for(i=0; i<reg2->size; i++)
+  if(!reg2->state)
     {
-      j = quantum_get_state(reg2->node[i].state, *reg1);
+      for(i=0; i<reg2->size; i++)
+	reg.amplitude[i] += reg2->amplitude[i];
+    }
 
-      if(j >= 0)
+  else
-	reg.node[j].amplitude += reg2->node[i].amplitude;
+    {
-
+      for(i=0; i<reg2->size; i++)
-      else
 	{
-	  reg.node[k].state = reg2->node[i].state;
+	  j = quantum_get_state(reg2->state[i], *reg1);
-	  reg.node[k].amplitude = reg2->node[i].amplitude;
+	  
-	  k++;
+	  if(j >= 0)
+	    reg.amplitude[j] += reg2->amplitude[i];
+
+	  else
+	    {
+	      reg.state[k] = reg2->state[i];
+	      reg.amplitude[k] = reg2->amplitude[i];
+	      k++;
+	    }
 	}
     }
   
@@ -559,41 +651,55 @@ quantum_vectoradd_inplace(quantum_reg *reg1, quantum_reg *reg2)
 
       for(i=0; i<reg2->size; i++)
 	{
-	  if(quantum_get_state(reg2->node[i].state, *reg1) == -1)
+	  if(quantum_get_state(reg2->state[i], *reg1) == -1)
 	    addsize++;
 	}
     }
 
-  /* Allocate memory for basis states */
+  if(addsize)
+    {
 
-  reg1->node = realloc(reg1->node, (reg1->size+addsize) 
+      /* Allocate memory for basis states */
-		       * sizeof(quantum_reg_node));
 
-  if(!reg1->node)
+      reg1->amplitude = realloc(reg1->amplitude, 
-    quantum_error(QUANTUM_ENOMEM);
+				(reg1->size+addsize)*sizeof(COMPLEX_FLOAT));
+      reg1->state = realloc(reg1->state, (reg1->size+addsize)
+			    *sizeof(MAX_UNSIGNED));
 
-  quantum_memman(addsize*sizeof(quantum_reg_node));
+      if(!(reg1->state && reg1->amplitude))
+	quantum_error(QUANTUM_ENOMEM);
 
-  /* Allocate the hash table */
+      quantum_memman(addsize * (sizeof(COMPLEX_FLOAT) + sizeof(MAX_UNSIGNED)));
+
+    }
 
   k = reg1->size;
 
-  for(i=0; i<reg2->size; i++)
+  if(!reg2->state)
     {
-      j = quantum_get_state(reg2->node[i].state, *reg1);
+      for(i=0; i<reg2->size; i++)
-
+	reg1->amplitude[i] += reg2->amplitude[i];
-      if(j >= 0)
-	reg1->node[j].amplitude += reg2->node[i].amplitude;
-
-      else
-	{
-	  reg1->node[k].state = reg2->node[i].state;
-	  reg1->node[k].amplitude = reg2->node[i].amplitude;
-	  k++;
-	}
     }
 
-  reg1->size += addsize;
+  else
+    {
+      for(i=0; i<reg2->size; i++)
+	{
+	  j = quantum_get_state(reg2->state[i], *reg1);
+
+	  if(j >= 0)
+	    reg1->amplitude[j] += reg2->amplitude[i];
+
+	  else
+	    {
+	      reg1->state[k] = reg2->state[i];
+	      reg1->amplitude[k] = reg2->amplitude[i];
+	      k++;
+	    }
+	}
+
+      reg1->size += addsize;
+    }
       
 }
 
@@ -607,31 +713,46 @@ quantum_reg
 quantum_matrix_qureg(quantum_reg A(MAX_UNSIGNED, double), double t,
 		     quantum_reg *reg, int flags)
 {
-  MAX_UNSIGNED i;
+  int i;
   quantum_reg reg2;
   quantum_reg tmp;
 
   reg2.width = reg->width;
-  reg2.size = 1 << reg2.width;
+  reg2.size = reg->size;
   reg2.hashw = 0;
   reg2.hash = 0;
 
-  reg2.node = calloc(reg2.size, sizeof(quantum_reg_node));
+  reg2.amplitude = calloc(reg2.size, sizeof(COMPLEX_FLOAT));
-  if(!reg2.node)
+  reg2.state = 0;
+
+  if(!reg2.amplitude)
     quantum_error(QUANTUM_ENOMEM);
 
-  quantum_memman(reg2.size*sizeof(quantum_reg_node));
+  quantum_memman(reg2.size * sizeof(COMPLEX_FLOAT));
 
-  for(i=0; i<(1<<reg->width); i++)
+  if(reg->state)
     {
-      reg2.node[i].state = i;
+      reg2.state = calloc(reg2.size, sizeof(MAX_UNSIGNED));
+
+      if(!reg2.state)
+	quantum_error(QUANTUM_ENOMEM);
+
+      quantum_memman(reg2.size * sizeof(MAX_UNSIGNED));
+    }
+
+#ifdef _OPENMP
+  #pragma omp parallel for private (tmp)
+#endif
+  for(i=0; i<reg->size; i++)
+    {
+      if(reg2.state)
+	reg2.state[i] = i;
       tmp = A(i, t);
-      reg2.node[i].amplitude = quantum_dot_product_noconj(&tmp, reg);
+      reg2.amplitude[i] = quantum_dot_product_noconj(&tmp, reg);
       if(!(flags & 1))
 	quantum_delete_qureg(&tmp);
     }
  
- 
   return reg2;
 
 }
@@ -645,9 +766,9 @@ quantum_mvmult(quantum_reg *y, quantum_matrix A, quantum_reg *x)
 
   for(i=0; i<A.cols; i++)
     {
-      y->node[i].amplitude = 0;
+      y->amplitude[i] = 0;
       for(j=0; j<A.cols; j++)
-	y->node[i].amplitude += M(A, j, i)*x->node[j].amplitude;
+	y->amplitude[i] += M(A, j, i)*x->amplitude[j];
     }
 } 
 
@@ -662,7 +783,7 @@ quantum_scalar_qureg(COMPLEX_FLOAT r, quantum_reg *reg)
   int i;
   
   for(i=0; i<reg->size; i++)
-      reg->node[i].amplitude *= r;
+      reg->amplitude[i] *= r;
 }
 
 /* Print the time evolution matrix for a series of gates */
@@ -681,7 +802,7 @@ quantum_print_timeop(int width, void f(quantum_reg *))
       tmp = quantum_new_qureg(i, width);
       f(&tmp);
       for(j=0; j<tmp.size; j++)
-	M(m, tmp.node[j].state, i) = tmp.node[j].amplitude;
+	M(m, tmp.state[j], i) = tmp.amplitude[j];
 
       quantum_delete_qureg(&tmp);
 	  
@@ -692,3 +813,18 @@ quantum_print_timeop(int width, void f(quantum_reg *))
   quantum_delete_matrix(&m);
   
 }
+
+/* Normalize a quantum register */
+
+void
+quantum_normalize(quantum_reg *reg)
+{
+  int i;
+  double r = 0;
+
+  for(i=0; i<reg->size; i++)
+    r += quantum_prob(reg->amplitude[i]);
+
+  quantum_scalar_qureg(1./sqrt(r), reg);
+
+}

qureg.h

@@ -1,6 +1,6 @@
 /* qureg.h: Declarations for qureg.c and inline hashing functions
 
-   Copyright 2003, 2004 Bjoern Butscher, Hendrik Weimer
+   Copyright 2003-2013 Bjoern Butscher, Hendrik Weimer
 
    This file is part of libquantum
 
@@ -31,16 +31,6 @@
 #include "matrix.h"
 #include "error.h"
 
-/* Representation of a base state of a quantum register: alpha_j |j> */
-
-struct quantum_reg_node_struct
-{
-  COMPLEX_FLOAT amplitude; /* alpha_j */
-  MAX_UNSIGNED state;      /* j */
-};
-
-typedef struct quantum_reg_node_struct quantum_reg_node;
-
 /* The quantum register */
 
 struct quantum_reg_struct
@@ -48,7 +38,8 @@ struct quantum_reg_struct
   int width;    /* number of qubits in the qureg */
   int size;     /* number of non-zero vectors */
   int hashw;    /* width of the hash array */
-  quantum_reg_node *node;
+  COMPLEX_FLOAT *amplitude;
+  MAX_UNSIGNED *state;
   int *hash;
 };
 
@@ -57,6 +48,7 @@ typedef struct quantum_reg_struct quantum_reg;
 extern quantum_reg quantum_matrix2qureg(quantum_matrix *m, int width);
 extern quantum_reg quantum_new_qureg(MAX_UNSIGNED initval, int width);
 extern quantum_reg quantum_new_qureg_size(int n, int width);
+extern quantum_reg quantum_new_qureg_sparse(int n, int width);
 extern quantum_matrix quantum_qureg2matrix(quantum_reg reg);
 extern void quantum_destroy_hash(quantum_reg *reg);
 extern void quantum_delete_qureg(quantum_reg *reg);
@@ -84,6 +76,7 @@ extern void quantum_scalar_qureg(COMPLEX_FLOAT r, quantum_reg *reg);
 extern void quantum_mvmult(quantum_reg *y, quantum_matrix A, quantum_reg *x);
 
 extern void quantum_print_timeop(int width, void f(quantum_reg *));
+extern void quantum_normalize(quantum_reg *reg);
 
 /* Our 64-bit multiplicative hash function */
 
@@ -114,7 +107,7 @@ quantum_get_state(MAX_UNSIGNED a, quantum_reg reg)
 
   while(reg.hash[i])
     {
-      if(reg.node[reg.hash[i]-1].state == a)
+      if(reg.state[reg.hash[i]-1] == a)
 	return reg.hash[i]-1;
       i++;
       if(i == (1 << reg.hashw))
@@ -168,7 +161,7 @@ quantum_reconstruct_hash(quantum_reg *reg)
   for(i=0; i<(1 << reg->hashw); i++)
     reg->hash[i] = 0;
   for(i=0; i<reg->size; i++)
-    quantum_add_hash(reg->node[i].state, i, reg);
+    quantum_add_hash(reg->state[i], i, reg);
 }
       
 /* Return the reduced bitmask of a basis state */

types.h.in

@@ -29,6 +29,10 @@
   #define COMPLEX_FLOAT @CF_TYPE@
 #endif
 
+#ifndef REAL_FLOAT
+  #define REAL_FLOAT @RF_TYPE@
+#endif
+
 #ifndef MAX_UNSIGNED
   #define MAX_UNSIGNED @MU_TYPE@
 #endif