lapack.c 3.73 KB
Newer Older
1 2
/* lapack.c: LAPACK interface

3
   Copyright 2008-2013 Hendrik Weimer
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

   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 <stdlib.h>
25
#include <math.h>
26 27 28 29 30 31 32 33 34 35 36 37

#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);

38 39 40 41
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);

42
void 
43
quantum_diag_time(double t, quantum_reg *reg0, quantum_reg *regt, 
44
		  quantum_reg *tmp1, quantum_reg *tmp2, quantum_matrix H, 
45
		  REAL_FLOAT **w)
46 47 48 49 50 51 52
{
#ifdef HAVE_LIBLAPACK
  char jobz = 'V';
  char uplo = 'U';
  int dim = H.cols;
  COMPLEX_FLOAT *work;
  int lwork = -1;
53
  REAL_FLOAT rwork[3*dim-2];
54
  int info;
55
  int i, j;
56 57 58 59 60 61
  void *p;
  
  if(tmp2->size != reg0->size)
    {
      /* perform diagonalization */

62 63 64 65 66 67 68 69 70 71 72
      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;
73
      *regt = *reg0;
74
      regt->amplitude = realloc(p, regt->size*sizeof(COMPLEX_FLOAT));
75
      
76
      p = tmp1->amplitude;
77
      *tmp1 = *reg0;
78
      tmp1->amplitude = realloc(p, regt->size*sizeof(COMPLEX_FLOAT));
79

80
      p = tmp2->amplitude;
81
      *tmp2 = *reg0;
82 83 84 85
      tmp2->amplitude = realloc(p, regt->size*sizeof(COMPLEX_FLOAT));

      if(!(regt->amplitude && tmp1->amplitude && tmp2->amplitude))
	quantum_error(QUANTUM_ENOMEM);
86 87 88 89 90 91 92 93 94 95 96

      *w = malloc(dim*sizeof(float));

      if(!*w)
	quantum_error(QUANTUM_ENOMEM);

      work = malloc(sizeof(COMPLEX_FLOAT));

      if(!work)
	quantum_error(QUANTUM_ENOMEM);

97 98
      QUANTUM_LAPACK_SOLVER(&jobz, &uplo, &dim, H.t, &dim, *w, work, &lwork, 
			    rwork, &info);
99 100 101 102 103

      if(info < 0)
	quantum_error(QUANTUM_ELAPACKARG);

      else if(info > 0)
104
	quantum_error(QUANTUM_ELAPACKCONV);
105 106 107 108 109 110 111
      
      lwork = (int) work[0];
      work = realloc(work, lwork*sizeof(COMPLEX_FLOAT));

      if(!work)
	quantum_error(QUANTUM_ENOMEM);

112 113
      QUANTUM_LAPACK_SOLVER(&jobz, &uplo, &dim, H.t, &dim, *w, work, &lwork, 
			    rwork, &info);
114 115 116 117 118

      if(info < 0)
	quantum_error(QUANTUM_ELAPACKARG);

      else if(info > 0)
119
	quantum_error(QUANTUM_ELAPACKCONV);
120 121 122 123 124 125 126 127 128 129
      
      free(work);

      quantum_mvmult(tmp1, H, reg0);

      quantum_adjoint(&H);
    }

  if(tmp1->size != reg0->size)
    {
130
      p = regt->amplitude;
131
      *regt = *reg0;
132
      regt->amplitude = realloc(p, regt->size*sizeof(COMPLEX_FLOAT));
133

134
      p = tmp1->amplitude;
135
      *tmp1 = *reg0;
136
      tmp1->amplitude = realloc(p, regt->size*sizeof(COMPLEX_FLOAT));
137 138 139 140 141 142 143 144 145

      quantum_adjoint(&H);
      
      quantum_mvmult(tmp1, H, reg0);

      quantum_adjoint(&H);
    }

  for(i=0; i<dim; i++)
146
    tmp2->amplitude[i] = quantum_cexp(-(*w)[i]*t)*tmp1->amplitude[i];
147 148 149 150 151 152 153 154 155 156 157

  quantum_mvmult(regt, H, tmp2);

#else
  quantum_error(QUANTUM_ENOLAPACK);

#endif /* HAVE_LIBLAPACK */
  
}