diff --git a/CHANGES b/CHANGES index edbb7e7..d6896c6 100644 --- a/CHANGES +++ b/CHANGES @@ -1,4 +1,7 @@ -libquantum 1.0.0: +libquantum 1.1.0: + - Added exact diagonlization based on LAPACK + - Added flag in quantum_rk4 to preserve qureg (breaks backward + compatibility) - Fixed quantum_gate1 to work properly with sorted regs - Fixed several bugs in quantum_rk4a diff --git a/Makefile.in b/Makefile.in index 1004f7e..6280197 100644 --- a/Makefile.in +++ b/Makefile.in @@ -1,6 +1,6 @@ # Makefile: Build libquantum # -# Copyright 2003-2005 Bjoern Butscher, Hendrik Weimer +# Copyright 2003-2008 Bjoern Butscher, Hendrik Weimer # # This file is part of libquantum # @@ -49,7 +49,7 @@ LIBTOOL=@LIBTOOL@ # Flags passed to C compiler CFLAGS=@CFLAGS@ -LDFLAGS=-rpath $(LIBDIR) -version-info 6:1:3 +LDFLAGS=-rpath $(LIBDIR) -version-info 7: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 Makefile + objcode.lo density.lo error.lo qtime.lo lapack.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 -lm + density.lo error.lo qtime.lo lapack.lo @LIBS@ complex.lo: complex.c complex.h config.h Makefile $(LIBTOOL) --mode=compile $(CC) $(CFLAGS) -c complex.c @@ -120,6 +120,9 @@ error.lo: error.c error.h Makefile qtime.lo: qtime.c qtime.h qureg.h Makefile $(LIBTOOL) --mode=compile $(CC) $(CFLAGS) -c qtime.c +lapack.lo: lapack.c lapack.h matrix.h qureg.h config.h error.h Makefile + $(LIBTOOL) --mode=compile $(CC) $(CFLAGS) -c lapack.c + # Autoconf stuff Makefile: config.status Makefile.in aclocal.m4 config.h.in types.h.in \ diff --git a/config.h.in b/config.h.in index b673056..e2fd729 100644 --- a/config.h.in +++ b/config.h.in @@ -1,6 +1,6 @@ /* config.h: configure settings - Copyright 2003-2005 Bjoern Butscher, Hendrik Weimer + Copyright 2003-2008 Bjoern Butscher, Hendrik Weimer This file is part of libquantum @@ -91,4 +91,7 @@ /* Define to 1 if you have GCC */ #undef HAVE_GCC +/* Define to 1 if you have LAPACK */ +#undef HAVE_LIBLAPACK + #include "types.h" diff --git a/configure b/configure index adfa4d4..1800d21 100755 --- a/configure +++ b/configure @@ -1,6 +1,6 @@ #! /bin/sh # Guess values for system-dependent variables and create Makefiles. -# Generated by GNU Autoconf 2.61 for libquantum 1.0.0. +# Generated by GNU Autoconf 2.61 for libquantum 1.1.0. # # Report bugs to . # @@ -724,8 +724,8 @@ SHELL=${CONFIG_SHELL-/bin/sh} # Identity of this package. PACKAGE_NAME='libquantum' PACKAGE_TARNAME='libquantum' -PACKAGE_VERSION='1.0.0' -PACKAGE_STRING='libquantum 1.0.0' +PACKAGE_VERSION='1.1.0' +PACKAGE_STRING='libquantum 1.1.0' PACKAGE_BUGREPORT='libquantum@libquantum.de' ac_unique_file="classic.c" @@ -1346,7 +1346,7 @@ if test "$ac_init_help" = "long"; then # Omit some internal or obsolete options to make the list less imposing. # This message is too long to be a string in the A/UX 3.1 sh. cat <<_ACEOF -\`configure' configures libquantum 1.0.0 to adapt to many kinds of systems. +\`configure' configures libquantum 1.1.0 to adapt to many kinds of systems. Usage: $0 [OPTION]... [VAR=VALUE]... @@ -1411,7 +1411,7 @@ fi if test -n "$ac_init_help"; then case $ac_init_help in - short | recursive ) echo "Configuration of libquantum 1.0.0:";; + short | recursive ) echo "Configuration of libquantum 1.1.0:";; esac cat <<\_ACEOF @@ -1432,9 +1432,11 @@ Optional Packages: --with-pic try to use only PIC/non-PIC objects [default=use both] --with-tags[=TAGS] include additional configurations [automatic] - --with-max-unsigned-type=ARG integer type for quantum registers - --with-complex-type=ARG type for complex numbers - --with-imaginary=ARG name of the imaginary unit + --with-lapack LAPACK support default=yes + --with-max-unsigned-type=ARG + integer type for quantum registers + --with-complex-type=ARG type for complex numbers + --with-imaginary=ARG name of the imaginary unit Some influential environment variables: CC C compiler command @@ -1510,7 +1512,7 @@ fi test -n "$ac_init_help" && exit $ac_status if $ac_init_version; then cat <<\_ACEOF -libquantum configure 1.0.0 +libquantum configure 1.1.0 generated by GNU Autoconf 2.61 Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001, @@ -1524,7 +1526,7 @@ cat >config.log <<_ACEOF This file contains any messages produced by compilers while running configure, to aid debugging if configure makes a mistake. -It was created by libquantum $as_me 1.0.0, which was +It was created by libquantum $as_me 1.1.0, which was generated by GNU Autoconf 2.61. Invocation command line was $ $0 $@ @@ -3695,7 +3697,7 @@ ia64-*-hpux*) ;; *-*-irix6*) # Find out which ABI we are using. - echo '#line 3698 "configure"' > conftest.$ac_ext + echo '#line 3700 "configure"' > conftest.$ac_ext if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5 (eval $ac_compile) 2>&5 ac_status=$? @@ -5414,11 +5416,11 @@ else -e 's:.*FLAGS}? :&$lt_compiler_flag :; t' \ -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ -e 's:$: $lt_compiler_flag:'` - (eval echo "\"\$as_me:5417: $lt_compile\"" >&5) + (eval echo "\"\$as_me:5419: $lt_compile\"" >&5) (eval "$lt_compile" 2>conftest.err) ac_status=$? cat conftest.err >&5 - echo "$as_me:5421: \$? = $ac_status" >&5 + echo "$as_me:5423: \$? = $ac_status" >&5 if (exit $ac_status) && test -s "$ac_outfile"; then # The compiler can only warn and ignore the option if not recognized # So say no if there are warnings @@ -5647,11 +5649,11 @@ else -e 's:.*FLAGS}? :&$lt_compiler_flag :; t' \ -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ -e 's:$: $lt_compiler_flag:'` - (eval echo "\"\$as_me:5650: $lt_compile\"" >&5) + (eval echo "\"\$as_me:5652: $lt_compile\"" >&5) (eval "$lt_compile" 2>conftest.err) ac_status=$? cat conftest.err >&5 - echo "$as_me:5654: \$? = $ac_status" >&5 + echo "$as_me:5656: \$? = $ac_status" >&5 if (exit $ac_status) && test -s "$ac_outfile"; then # The compiler can only warn and ignore the option if not recognized # So say no if there are warnings @@ -5707,11 +5709,11 @@ else -e 's:.*FLAGS}? :&$lt_compiler_flag :; t' \ -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ -e 's:$: $lt_compiler_flag:'` - (eval echo "\"\$as_me:5710: $lt_compile\"" >&5) + (eval echo "\"\$as_me:5712: $lt_compile\"" >&5) (eval "$lt_compile" 2>out/conftest.err) ac_status=$? cat out/conftest.err >&5 - echo "$as_me:5714: \$? = $ac_status" >&5 + echo "$as_me:5716: \$? = $ac_status" >&5 if (exit $ac_status) && test -s out/conftest2.$ac_objext then # The compiler can only warn and ignore the option if not recognized @@ -7867,7 +7869,7 @@ else lt_dlunknown=0; lt_dlno_uscore=1; lt_dlneed_uscore=2 lt_status=$lt_dlunknown cat > conftest.$ac_ext < conftest.$ac_ext <&5 +echo $ECHO_N "checking for cheev_ in -llapack... $ECHO_C" >&6; } +if test "${ac_cv_lib_lapack_cheev_+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + ac_check_lib_save_LIBS=$LIBS +LIBS="-llapack $LIBS" +cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ + +/* Override any GCC internal prototype to avoid an error. + Use char because int might match the return type of a GCC + builtin and then its argument prototype would still apply. */ +#ifdef __cplusplus +extern "C" +#endif +char cheev_ (); +int +main () +{ +return cheev_ (); + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext conftest$ac_exeext +if { (ac_try="$ac_link" +case "(($ac_try" in + *\"* | *\`* | *\\*) ac_try_echo=\$ac_try;; + *) ac_try_echo=$ac_try;; +esac +eval "echo \"\$as_me:$LINENO: $ac_try_echo\"") >&5 + (eval "$ac_link") 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && { + test -z "$ac_c_werror_flag" || + test ! -s conftest.err + } && test -s conftest$ac_exeext && + $as_test_x conftest$ac_exeext; then + ac_cv_lib_lapack_cheev_=yes +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + + ac_cv_lib_lapack_cheev_=no +fi + +rm -f core conftest.err conftest.$ac_objext conftest_ipa8_conftest.oo \ + conftest$ac_exeext conftest.$ac_ext +LIBS=$ac_check_lib_save_LIBS +fi +{ echo "$as_me:$LINENO: result: $ac_cv_lib_lapack_cheev_" >&5 +echo "${ECHO_T}$ac_cv_lib_lapack_cheev_" >&6; } +if test $ac_cv_lib_lapack_cheev_ = yes; then + cat >>confdefs.h <<_ACEOF +#define HAVE_LIBLAPACK 1 +_ACEOF + + LIBS="-llapack $LIBS" + +fi + + fi +else + +{ echo "$as_me:$LINENO: checking for cheev_ in -llapack" >&5 +echo $ECHO_N "checking for cheev_ in -llapack... $ECHO_C" >&6; } +if test "${ac_cv_lib_lapack_cheev_+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + ac_check_lib_save_LIBS=$LIBS +LIBS="-llapack $LIBS" +cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ + +/* Override any GCC internal prototype to avoid an error. + Use char because int might match the return type of a GCC + builtin and then its argument prototype would still apply. */ +#ifdef __cplusplus +extern "C" +#endif +char cheev_ (); +int +main () +{ +return cheev_ (); + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext conftest$ac_exeext +if { (ac_try="$ac_link" +case "(($ac_try" in + *\"* | *\`* | *\\*) ac_try_echo=\$ac_try;; + *) ac_try_echo=$ac_try;; +esac +eval "echo \"\$as_me:$LINENO: $ac_try_echo\"") >&5 + (eval "$ac_link") 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && { + test -z "$ac_c_werror_flag" || + test ! -s conftest.err + } && test -s conftest$ac_exeext && + $as_test_x conftest$ac_exeext; then + ac_cv_lib_lapack_cheev_=yes +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + + ac_cv_lib_lapack_cheev_=no +fi + +rm -f core conftest.err conftest.$ac_objext conftest_ipa8_conftest.oo \ + conftest$ac_exeext conftest.$ac_ext +LIBS=$ac_check_lib_save_LIBS +fi +{ echo "$as_me:$LINENO: result: $ac_cv_lib_lapack_cheev_" >&5 +echo "${ECHO_T}$ac_cv_lib_lapack_cheev_" >&6; } +if test $ac_cv_lib_lapack_cheev_ = yes; then + cat >>confdefs.h <<_ACEOF +#define HAVE_LIBLAPACK 1 +_ACEOF + + LIBS="-llapack $LIBS" + +fi + +fi + + + # Checks for header files. { echo "$as_me:$LINENO: checking for ANSI C header files" >&5 echo $ECHO_N "checking for ANSI C header files... $ECHO_C" >&6; } @@ -9700,6 +9854,9 @@ if test "${enable_profiling+set}" = set; then fi +# Disable LAPACK check + + # Enable -Wall for gcc if test $CC = "gcc" then @@ -10109,7 +10266,7 @@ exec 6>&1 # report actual input values of CONFIG_FILES etc. instead of their # values after options handling. ac_log=" -This file was extended by libquantum $as_me 1.0.0, which was +This file was extended by libquantum $as_me 1.1.0, which was generated by GNU Autoconf 2.61. Invocation command line was CONFIG_FILES = $CONFIG_FILES @@ -10158,7 +10315,7 @@ Report bugs to ." _ACEOF cat >>$CONFIG_STATUS <<_ACEOF ac_cs_version="\\ -libquantum config.status 1.0.0 +libquantum config.status 1.1.0 configured by $0, generated by GNU Autoconf 2.61, with options \\"`echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`\\" diff --git a/configure.in b/configure.in index 3750490..9be62f8 100644 --- a/configure.in +++ b/configure.in @@ -20,7 +20,7 @@ # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, # MA 02110-1301, USA -AC_INIT([libquantum], [1.0.0], [libquantum@libquantum.de]) +AC_INIT([libquantum], [1.1.0], [libquantum@libquantum.de]) AC_CONFIG_SRCDIR([classic.c]) AC_CONFIG_HEADER([config.h]) @@ -34,6 +34,13 @@ AC_PROG_LIBTOOL # Checks for libraries. AC_CHECK_LIB([m], [sqrt]) +AC_ARG_WITH(lapack, + [ --with-lapack LAPACK support [default=yes]], + [if test $withval = "yes" + then + AC_CHECK_LIB([lapack], [cheev_]) + fi] , [AC_CHECK_LIB([lapack], [cheev_])]) + # Checks for header files. AC_HEADER_STDC @@ -44,7 +51,8 @@ AC_C_INLINE # Check for 64-bit integer AC_ARG_WITH([max-unsigned-type], - [ --with-max-unsigned-type=ARG integer type for quantum registers], + [ --with-max-unsigned-type=ARG + integer type for quantum registers], [MU_TYPE=$withval], [MU_TYPE="none" AC_CHECK_TYPE([uint_64t], [AC_DEFINE([MAX_UNSIGNED], [uint_64t]) MU_TYPE="uint_64t"])]) @@ -66,7 +74,7 @@ fi # Check for complex number type AC_ARG_WITH([complex-type], - [ --with-complex-type=ARG type for complex numbers], + [ --with-complex-type=ARG type for complex numbers], [CF_TYPE=$withval ], [CF_TYPE="none" AC_CHECK_TYPE([float _Complex], [AC_DEFINE([COMPLEX_FLOAT], @@ -86,7 +94,7 @@ fi # Check for the imaginary unit AC_MSG_CHECKING([for the imaginary unit]) AC_ARG_WITH([imaginary], - [ --with-imaginary=ARG name of the imaginary unit], + [ --with-imaginary=ARG name of the imaginary unit], [I=$withval], [I="none" AC_COMPILE_IFELSE([AC_LANG_PROGRAM([$CF_TYPE z;], [z=I;])], [AC_DEFINE([IMAGINARY], [I], [Imaginary unit]) I="I"]) @@ -112,6 +120,9 @@ 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 diff --git a/error.c b/error.c index d868dbd..18c284c 100644 --- a/error.c +++ b/error.c @@ -54,6 +54,12 @@ quantum_strerr(int errno) return "wrong matrix size"; case QUANTUM_EHASHFULL: return "hash table full"; + case QUANTUM_ENOLAPACK: + return "LAPACK support not compiled in"; + case QUANTUM_ELAPACKARG: + return "wrong arguments supplied to LAPACK"; + case QUANTUM_ELAPACKCHEEV: + return "LAPACK's CHEEV failed to converge"; case QUANTUM_EMCMATRIX: return "single-column matrix expected"; case QUANTUM_EOPCODE: diff --git a/error.h b/error.h index 4f699dd..34b4f9a 100644 --- a/error.h +++ b/error.h @@ -26,14 +26,17 @@ #define __ERROR_H enum { - QUANTUM_SUCCESS = 0, - QUANTUM_FAILURE = 1, - QUANTUM_ENOMEM = 2, - QUANTUM_EMLARGE = 3, - QUANTUM_EMSIZE = 4, - QUANTUM_EHASHFULL = 5, - QUANTUM_EMCMATRIX = 65536, /* internal errors start at 65536 */ - QUANTUM_EOPCODE = 65537 + QUANTUM_SUCCESS = 0, + QUANTUM_FAILURE = 1, + QUANTUM_ENOMEM = 2, + QUANTUM_EMLARGE = 3, + QUANTUM_EMSIZE = 4, + QUANTUM_EHASHFULL = 5, + QUANTUM_ENOLAPACK = 32768, /* LAPACK errors start at 32768 */ + QUANTUM_ELAPACKARG = 32769, + QUANTUM_ELAPACKCHEEV = 32770, + QUANTUM_EMCMATRIX = 65536, /* internal errors start at 65536 */ + QUANTUM_EOPCODE = 65537 }; extern void *quantum_error_handler(void *f(int)); diff --git a/lapack.c b/lapack.c new file mode 100644 index 0000000..76ff101 --- /dev/null +++ b/lapack.c @@ -0,0 +1,147 @@ +/* lapack.c: LAPACK interface + + Copyright 2008 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 + +#include "lapack.h" +#include "matrix.h" +#include "complex.h" +#include "qureg.h" +#include "error.h" +#include "config.h" + +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); + +void +quantum_diag_time(float t, quantum_reg *reg0, quantum_reg *regt, + quantum_reg *tmp1, quantum_reg *tmp2, quantum_matrix H, + float **w) +{ +#ifdef HAVE_LIBLAPACK + char jobz = 'V'; + char uplo = 'U'; + int dim = H.cols; + COMPLEX_FLOAT *work; + int lwork = -1; + float rwork[3*dim-2]; + int info; + int i; + void *p; + + if(tmp2->size != reg0->size) + { + /* perform diagonalization */ + + p = regt->node; + *regt = *reg0; + regt->node = realloc(p, regt->size*sizeof(quantum_reg_node)); + for(i=0; isize; i++) + regt->node[i].state = i; + + p = tmp1->node; + *tmp1 = *reg0; + tmp1->node = realloc(p, regt->size*sizeof(quantum_reg_node)); + for(i=0; isize; i++) + tmp1->node[i].state = i; + + p = tmp2->node; + *tmp2 = *reg0; + tmp2->node = realloc(p, regt->size*sizeof(quantum_reg_node)); + for(i=0; isize; i++) + tmp2->node[i].state = i; + + *w = malloc(dim*sizeof(float)); + + if(!*w) + quantum_error(QUANTUM_ENOMEM); + + work = malloc(sizeof(COMPLEX_FLOAT)); + + if(!work) + quantum_error(QUANTUM_ENOMEM); + + cheev_(&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); + + lwork = (int) work[0]; + work = realloc(work, lwork*sizeof(COMPLEX_FLOAT)); + + if(!work) + quantum_error(QUANTUM_ENOMEM); + + cheev_(&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); + + free(work); + + quantum_mvmult(tmp1, H, reg0); + + quantum_adjoint(&H); + } + + if(tmp1->size != reg0->size) + { + p = regt->node; + *regt = *reg0; + regt->node = realloc(p, regt->size*sizeof(quantum_reg_node)); + for(i=0; isize; i++) + regt->node[i].state = i; + + p = tmp1->node; + *tmp1 = *reg0; + tmp1->node = realloc(p, regt->size*sizeof(quantum_reg_node)); + for(i=0; isize; i++) + tmp1->node[i].state = i; + + quantum_adjoint(&H); + + quantum_mvmult(tmp1, H, reg0); + + quantum_adjoint(&H); + } + + for(i=0; inode[i].amplitude = quantum_cexp(-(*w)[i]*t)*tmp1->node[i].amplitude; + + quantum_mvmult(regt, H, tmp2); + +#else + quantum_error(QUANTUM_ENOLAPACK); + +#endif /* HAVE_LIBLAPACK */ + +} + + diff --git a/lapack.h b/lapack.h new file mode 100644 index 0000000..ec21050 --- /dev/null +++ b/lapack.h @@ -0,0 +1,36 @@ +/* time.c: Declarations for lapack.h + + Copyright 2006 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 + +*/ + +#ifndef __LAPACK_H + +#define __LAPACK_H + +#include "config.h" +#include "matrix.h" +#include "qureg.h" + +extern void quantum_diag_time(float t, quantum_reg *reg0, quantum_reg *regt, + quantum_reg *tmp1, quantum_reg *tmp2, + quantum_matrix H, float **w); + +#endif diff --git a/matrix.c b/matrix.c index e039753..5485e57 100644 --- a/matrix.c +++ b/matrix.c @@ -134,4 +134,24 @@ quantum_matrix quantum_mmult(quantum_matrix A, quantum_matrix B) return C; } - + +/* Compute the adjoint of a matrix */ + +void +quantum_adjoint(quantum_matrix *m) +{ + int i, j; + COMPLEX_FLOAT tmp; + quantum_matrix A = *m; + + for(i=0; icols; i++) + { + for(j=0;jhashw = 0; /* k1 */ - k = quantum_matrix_qureg(H, t, reg); + k = quantum_matrix_qureg(H, t, reg, flags); quantum_scalar_qureg(-IMAGINARY*dt/2.0, &k); tmp = quantum_vectoradd(reg, &k); quantum_scalar_qureg(1.0/3.0, &k); @@ -57,7 +57,7 @@ quantum_rk4(quantum_reg *reg, double t, double dt, quantum_delete_qureg(&k); /* k2 */ - k = quantum_matrix_qureg(H, t+dt/2.0, &tmp); + k = quantum_matrix_qureg(H, t+dt/2.0, &tmp, flags); quantum_delete_qureg(&tmp); quantum_scalar_qureg(-IMAGINARY*dt/2.0, &k); tmp = quantum_vectoradd(reg, &k); @@ -66,7 +66,7 @@ quantum_rk4(quantum_reg *reg, double t, double dt, quantum_delete_qureg(&k); /* k3 */ - k = quantum_matrix_qureg(H, t+dt/2.0, &tmp); + k = quantum_matrix_qureg(H, t+dt/2.0, &tmp, flags); quantum_delete_qureg(&tmp); quantum_scalar_qureg(-IMAGINARY*dt, &k); tmp = quantum_vectoradd(reg, &k); @@ -75,7 +75,7 @@ quantum_rk4(quantum_reg *reg, double t, double dt, quantum_delete_qureg(&k); /* k4 */ - k = quantum_matrix_qureg(H, t+dt, &tmp); + k = quantum_matrix_qureg(H, t+dt, &tmp, flags); quantum_delete_qureg(&tmp); quantum_scalar_qureg(-IMAGINARY*dt/6.0, &k); quantum_vectoradd_inplace(&out, &k); @@ -103,7 +103,7 @@ quantum_rk4(quantum_reg *reg, double t, double dt, double quantum_rk4a(quantum_reg *reg, double t, double *dt, double epsilon, - quantum_reg H(MAX_UNSIGNED, double)) + quantum_reg H(MAX_UNSIGNED, double), int flags) { quantum_reg reg2, old; double delta, r, dtused; @@ -122,9 +122,9 @@ quantum_rk4a(quantum_reg *reg, double t, double *dt, double epsilon, do { - quantum_rk4(reg, t, *dt, H); - quantum_rk4(®2, t, *dt/2.0, H); - quantum_rk4(®2, t, *dt/2.0, H); + quantum_rk4(reg, t, *dt, H, flags); + quantum_rk4(®2, t, *dt/2.0, H, flags); + quantum_rk4(®2, t+*dt/2.0, *dt/2.0, H, flags); delta = 0; diff --git a/qtime.h b/qtime.h index b90244b..59a8256 100644 --- a/qtime.h +++ b/qtime.h @@ -29,9 +29,9 @@ #include "config.h" extern void quantum_rk4(quantum_reg *reg, double t, double dt, - quantum_reg H(MAX_UNSIGNED, double)); + quantum_reg H(MAX_UNSIGNED, double), int flags); extern double quantum_rk4a(quantum_reg *reg, double t, double *dt, double epsilon, - quantum_reg H(MAX_UNSIGNED, double)); + quantum_reg H(MAX_UNSIGNED, double), int flags); #endif diff --git a/quantum.h.in b/quantum.h.in index c125499..e020828 100644 --- a/quantum.h.in +++ b/quantum.h.in @@ -1,6 +1,6 @@ /* quantum.h: Header file for libquantum - Copyright 2003-2007 Bjoern Butscher, Hendrik Weimer + Copyright 2003-2008 Bjoern Butscher, Hendrik Weimer This file is part of libquantum @@ -138,7 +138,7 @@ extern COMPLEX_FLOAT quantum_dot_product(quantum_reg *reg1, quantum_reg *reg2); extern quantum_reg quantum_vectoradd(quantum_reg *reg1, quantum_reg *reg2); extern void quantum_vectoradd_inplace(quantum_reg *reg1, quantum_reg *reg2); extern quantum_reg quantum_matrix_qureg(quantum_reg A(MAX_UNSIGNED, double), - double t, quantum_reg *reg); + double t, quantum_reg *reg, int flags); extern void quantum_scalar_qureg(COMPLEX_FLOAT r, quantum_reg *reg); extern void quantum_print_timeop(int width, void f(quantum_reg *)); @@ -166,9 +166,13 @@ extern const char *quantum_strerr(int errno); extern void quantum_error(int errno); extern void quantum_rk4(quantum_reg *reg, double t, double dt, - quantum_reg H(MAX_UNSIGNED, double)); + quantum_reg H(MAX_UNSIGNED, double), int flags); extern double quantum_rk4a(quantum_reg *reg, double t, double *dt, double epsilon, - quantum_reg H(MAX_UNSIGNED, double)); + quantum_reg H(MAX_UNSIGNED, double), int flags); + +extern void quantum_diag_time(float t, quantum_reg *reg0, quantum_reg *regt, + quantum_reg *tmp1, quantum_reg *tmp2, + quantum_matrix H, float **w); #endif diff --git a/qureg.c b/qureg.c index b0a3fe8..56f5dc9 100644 --- a/qureg.c +++ b/qureg.c @@ -333,7 +333,7 @@ quantum_kronecker(quantum_reg *reg1, quantum_reg *reg2) reg.width = reg1->width+reg2->width; reg.size = reg1->size*reg2->size; - reg.hashw = reg.width+2; + reg.hashw = reg.width + 2; /* allocate memory for the new basis states */ @@ -605,7 +605,7 @@ quantum_vectoradd_inplace(quantum_reg *reg1, quantum_reg *reg2) quantum_reg quantum_matrix_qureg(quantum_reg A(MAX_UNSIGNED, double), double t, - quantum_reg *reg) + quantum_reg *reg, int flags) { MAX_UNSIGNED i; quantum_reg reg2; @@ -627,7 +627,8 @@ quantum_matrix_qureg(quantum_reg A(MAX_UNSIGNED, double), double t, reg2.node[i].state = i; tmp = A(i, t); reg2.node[i].amplitude = quantum_dot_product_noconj(&tmp, reg); - quantum_delete_qureg(&tmp); + if(!(flags & 1)) + quantum_delete_qureg(&tmp); } @@ -635,6 +636,22 @@ quantum_matrix_qureg(quantum_reg A(MAX_UNSIGNED, double), double t, } +/* Matrix-vector multiplication using a quantum_matrix */ + +void +quantum_mvmult(quantum_reg *y, quantum_matrix A, quantum_reg *x) +{ + int i, j; + + for(i=0; inode[i].amplitude = 0; + for(j=0; jnode[i].amplitude += M(A, j, i)*x->node[j].amplitude; + } +} + + /* Scalar multiplication of a quantum register. This is a purely mathematical operation without any physical meaning, so only use it if you know what you are doing. */ diff --git a/qureg.h b/qureg.h index 2ee5871..1ab3eda 100644 --- a/qureg.h +++ b/qureg.h @@ -79,8 +79,9 @@ extern COMPLEX_FLOAT quantum_dot_product(quantum_reg *reg1, quantum_reg *reg2); extern quantum_reg quantum_vectoradd(quantum_reg *reg1, quantum_reg *reg2); extern void quantum_vectoradd_inplace(quantum_reg *reg1, quantum_reg *reg2); extern quantum_reg quantum_matrix_qureg(quantum_reg A(MAX_UNSIGNED, double), - double t, quantum_reg *reg); + double t, quantum_reg *reg, int flags); 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 *));