utilities.c

Go to the documentation of this file.

/*! \file utilities.c
    \brief See documentation of header.h - this file contains source code only, with function prototypes in header.h.
 
 *  Numerical State Tomography
 *  Using GSL -> www.gnu.org/software/gsl
 *  Copyright (C) 2007  Max S Kaznady
 *
 *  This program 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;
 *
 *  AUTHOR
 *     Max S Kaznady,  [EMAIL PROTECTED]
 *     max.kaznady@gmail.com
 *     Summer, 2007
 *
 *
 ******************************************************************************/


#include <stdio.h>
#include <math.h>
#include <stdlib.h> //for writing files
#include <string.h>
#include <time.h>
#include <sys/time.h>

#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_math.h>

#include "header.h"

void
kron(gsl_matrix_complex *a, gsl_matrix_complex *b, gsl_matrix_complex *c) {

    //gsl_matrix_complex *c =
    //gsl_matrix_complex_alloc((a->size1)*(b->size1), (a->size1)*(b->size1));
    if ((c->size1)!=(a->size1)*(b->size1) || (c->size2)!=(a->size2)*(b->size2)) {
        fprintf(stderr, "kron: invalid size of output matrix\n");
        fprintf(stderr, "expected: size %d by %d\n", (a->size1)*(b->size1), (a->size1)*(b->size1));
        fprintf(stderr, "got: size %d by %d\n", c->size1, c->size2 );
        exit(1);
    }
    gsl_complex x, y;
    int i, j, k, l;
    for (i = 0; i<(a->size1); i++) {
        for (j = 0; j<(a->size2); j++) {
            for (k = 0; k<(b->size1); k++) {
                for (l = 0; l<(b->size2); l++) {
                    x = gsl_matrix_complex_get(a, i, j);
                    y = gsl_matrix_complex_get(b, k, l);
                    //printf("at pos %d %d\n", i*(b->size1)+k, j*(b->size2)+l);
                    gsl_matrix_complex_set(c, i*(b->size1)+k, j*(b->size2)+l,
                                           gsl_complex_mul(x,y));

                }
            }
        }
    }

    //printf("exiting kron\n");
    //return c;
}

void
outer_product(gsl_vector_complex *a, gsl_matrix_complex *c) {
    //check the dimensions
    if ((c->size1)!=(a->size) || (c->size2)!=(a->size)) {
        fprintf(stderr, "outer product: size mismatch\n");
        exit(1);
    }
    //gsl_matrix_complex *c = gsl_matrix_complex_alloc(a->size, a->size);
    gsl_complex x, y;
    int i, j;
    for (i=0; i<(a->size); i++)
        for (j=0; j<(a->size); j++) {
            x = gsl_vector_complex_get(a, i);
            y = gsl_vector_complex_get(a, j);
            gsl_matrix_complex_set(c, i, j, gsl_complex_mul(x, gsl_complex_conjugate(y)));
        }

}

gsl_complex
trace(gsl_matrix_complex *m) {
    if ((m->size1)!=m->size2) {
        fprintf(stderr, "Matrix must be square for trace to work\n");
        exit(1);
    }
    gsl_complex tr;
    gsl_complex temp;
    int i;
    GSL_SET_COMPLEX(&tr, 0, 0);
    for (i = 0; i<(m->size1); i++) {
        //have to make an awkward switch below - just in case
        temp=gsl_complex_add(tr, gsl_matrix_complex_get(m, i, i));
        tr=temp;
        //printf("%e\n", GSL_REAL(tr));
    }
    return tr;
}


//returns a complex random matrix, sampled from Uniform(-1, 1)
void
matrix_complex_random_init(gsl_matrix_complex *random) {
    gsl_complex rnd;
    gsl_rng * r;
    if ((r= gsl_rng_alloc (gsl_rng_taus)) == NULL) {
        perror("gsl_rng_alloc");
        exit(1);
    }

    //set seed - random - use local time in microseconds
    struct timeval tv;
    struct timezone tz;
    struct tm *tm;
    gettimeofday(&tv, &tz);
    tm = localtime(&tv.tv_sec);
    if (VERBOSE)
        printf(" %d:%02d:%02d %d \n", tm->tm_hour, tm->tm_min, tm->tm_sec, (int)tv.tv_usec);
    //gsl_rng_set (r, time(NULL));
    gsl_rng_set (r, tv.tv_usec);

    int i, j;
    for (i=0; i<(random->size1); i++) {
        for (j=0; j<(random->size2); j++) {
            //printf("%e\n", 2*gsl_rng_uniform_pos(r)-1);
            //set seed
            //gsl_rng_set (r, time(NULL));
            GSL_SET_COMPLEX(&rnd, 2*gsl_rng_uniform_pos(r)-1,2*gsl_rng_uniform_pos(r)-1);
            gsl_matrix_complex_set(random, i, j, rnd);
        }
    }
    gsl_rng_free(r);
}

//sets a random density matrix
void
random_density_matrix(gsl_matrix_complex *density_random) {
    gsl_matrix_complex *random = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    matrix_complex_random_init(random);
    gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one, random, random, zero, density_random);
    //renormalize trace
    gsl_matrix_complex_scale(density_random, gsl_complex_rect(1/GSL_REAL(trace(density_random)), 0));
    gsl_matrix_complex_free(random);
}

//increment counter base-4 rep from nu-1 to nu
void
base_4_rep_inc(int counter[], int n) {
    int c = n - 1;
    int truth = 1, subtruth;
    int current = NQ-1;
    while (truth) {
        if (c==n) //reached the increment, terminate
            truth=0;
        else { //increment the counter
            if (counter[current] == 3) { //last digit is 3, need to shift
                subtruth = 1;
                while (subtruth) {
                    counter[current]=0; //count up current value
                    current--; //go to next number
                    if (counter[current]!=3) { //if can increment, do so
                        counter[current]++;
                        subtruth = 0;
                    } //otherwise zero out again, approach next available number
                }
            } else { //simple case - here we can just increment the number
                counter[current]++;
            }
        }
        c++;
    }
}

void
mu_tensor(gsl_matrix_complex *mu_cell[], gsl_matrix_complex *mu_array[], int mu_counter[], int counter[]) {
    //gsl_matrix_complex *mu_nu = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    int i, j;
    for (i=0; i<NQ; i++) {
        if (counter[i]!=mu_counter[i]) { //locate place from which to update
            if (i==0) { //base case
                //mu_cell[i] = mu_array[0];
                //mu_cell[0] = gsl_matrix_complex_alloc(2, 2);
                gsl_matrix_complex_memcpy(mu_cell[0], mu_array[counter[i]]);
                i++;
            }
            for (j=i; j<NQ; j++) { //update the rest of elements
                //gsl_matrix_complex_free(mu_cell[i]);
                //printf("mu_tensor, j=%d and i=%d, counter[j]=%d \n", j, i, counter[j]);
                //mu_cell[j] = gsl_matrix_complex_alloc((mu_cell[j-1]->size1)*2, (mu_cell[j-1]->size2)*2);
                kron(mu_cell[j-1], mu_array[counter[j]], mu_cell[j]);
            }
            break; // done updating
        }
    }
}

void
sig_tensor(gsl_matrix_complex *sigma_cell[], gsl_matrix_complex *sig_array[], int mu_counter[], int counter[]) {
    //gsl_matrix_complex *mu_nu = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    int i, j;
    for (i=0; i<NQ; i++) {
        if (counter[i]!=mu_counter[i]) { //locate place from which to update
            if (i==0) { //base case
                //mu_cell[i] = mu_array[0];
                //mu_cell[0] = gsl_matrix_complex_alloc(2, 2);
                gsl_matrix_complex_memcpy(sigma_cell[0], sig_array[counter[i]]);
                i++;
            }
            for (j=i; j<NQ; j++) { //update the rest of elements
                //gsl_matrix_complex_free(mu_cell[i]);
                //printf("mu_tensor, j=%d and i=%d, counter[j]=%d \n", j, i, counter[j]);
                //sigma_cell[j] = gsl_matrix_complex_alloc((mu_cell[j-1]->size1)*2, (mu_cell[j-1]->size2)*2);
                kron(sigma_cell[j-1], sig_array[counter[j]], sigma_cell[j]);
            }
            break; // done updating
        }
    }
}

int
write_data(char *name, gsl_matrix_complex *data) {
    FILE *fp_real, *fp_imag;

    char file_real[MAXCHAR], file_imag[MAXCHAR];

    if ((strncpy(file_real, DATA, MAXCHAR)) == NULL) {
        perror("strncpy");
        exit(1);
    }
    if ((strncat(file_real, name, MAXCHAR)) == NULL) {
        perror("strncat");
        exit(1);
    }
    if ((strncat(file_real, "_real.txt", MAXCHAR)) == NULL) {
        perror("strncat");
        exit(1);
    }
    if ((strncpy(file_imag, DATA, MAXCHAR)) == NULL) {
        perror("strncpy");
        exit(1);
    }
    if ((strncat(file_imag, name, MAXCHAR)) == NULL) {
        perror("strncat");
        exit(1);
    }
    if ((strncat(file_imag, "_imag.txt", MAXCHAR)) == NULL) {
        perror("strncat");
        exit(1);
    }

    //printf("directory %s \n", file_real);
    //printf("directory %s \n", file_imag);
    //open a file for write operation
    if ((fp_real = fopen(file_real, "w")) == NULL) {
        perror("fopen");
        exit(1);
    }
    if ((fp_imag = fopen(file_imag, "w")) == NULL) {
        perror("fopen");
        exit(1);
    }
    int i, j;
    double number;
    char tab='\t';
    char newline='\n';
    char *write = malloc(MAXCHAR);
    char buffer[MAXCHAR];
    gsl_complex temp;
    for (i=0; i<(data->size1); i++) {
        for (j=0; j<(data->size2); j++) {

            temp = gsl_matrix_complex_get(data, i, j);

            number = GSL_REAL(temp);

            //printf("%d\n", sizeof(write));
            snprintf(write, MAXCHAR, "%e", number);
            strncpy(buffer, write, MAXCHAR);

            Fwrite(buffer, sizeof(buffer), 1, fp_real);
            Fwrite(&tab, sizeof(tab), 1, fp_real);

            //problem here
            number = GSL_IMAG(temp);
            snprintf(write, MAXCHAR, "%e", number);
            strncpy(buffer, write, MAXCHAR);
            //printf("%s\n", write);
            Fwrite(buffer, sizeof(buffer), 1, fp_imag);
            Fwrite(&tab, sizeof(tab), 1, fp_imag);
        }
        Fwrite(&newline, sizeof(newline), 1, fp_real);
        Fwrite(&newline, sizeof(newline), 1, fp_imag);
    }

    if ((fclose(fp_real)) != 0) {
        perror("fclose");
        exit(1);
    }

    if ((fclose(fp_imag)) != 0) {
        perror("fclose");
        exit(1);
    }


    free(write);
    return 0;
}

//wrapper function
int
Fwrite(const void *ptr, size_t size, size_t nmemb, FILE *stream) {
    int n = fwrite(ptr, size, nmemb, stream);
    if(n < nmemb) {
        perror("fwrite");
        exit(1);
    }
    return(n);
}

void
display_matrix(gsl_matrix_complex *data) {
    int i, j;
    for (i=0; i<(data->size1); i++) {
        for (j=0; j<(data->size2); j++) {
            printf("%e+i%e\t", GSL_REAL(gsl_matrix_complex_get(data, i, j)),
                   GSL_IMAG(gsl_matrix_complex_get(data, i, j)));
        }
        printf("\n");
    }
    //printf("size2: %d\n", data->size2);
}

void
display_matrix_real(gsl_matrix *data) {
    int i, j;
    for (i=0; i<(data->size1); i++) {
        for (j=0; j<(data->size2); j++) {
            printf("%e\t", gsl_matrix_get(data, i, j));
        }
        printf("\n");
    }
    //printf("size2: %d\n", data->size2);
}

void
make_T(const gsl_vector *t, gsl_matrix_complex *T) {
    gsl_matrix_complex_set_zero(T);
    int limiter, nu;
    int counter = 0;
    for (limiter = 0; limiter<pow(2, NQ); limiter++) {
        for (nu=limiter; nu<pow(2, NQ); nu++) {
            if (counter<pow(2,NQ)) {
                gsl_matrix_complex_set(T, nu, nu, gsl_complex_rect(gsl_vector_get(t, counter),0));
            } else {
                gsl_matrix_complex_set(T, nu, nu-limiter, gsl_complex_rect(gsl_vector_get(t,counter),gsl_vector_get(t,counter+1)));
                //nu++;
                counter++;
            }
            counter++;
        }
    }

}

int
write_vector_real(char *name, gsl_vector *data) {
    FILE *fp;

    char file[MAXCHAR];

    if ((strncpy(file, DATA, MAXCHAR)) == NULL) {
        perror("strncpy");
        exit(1);
    }
    if ((strncat(file, name, MAXCHAR)) == NULL) {
        perror("strncat");
        exit(1);
    }
    if ((strncat(file, ".txt", MAXCHAR)) == NULL) {
        perror("strncat");
        exit(1);
    }

    //printf("directory %s \n", file_real);
    //printf("directory %s \n", file_imag);
    //open a file for write operation
    if ((fp = fopen(file, "w")) == NULL) {
        perror("fopen");
        exit(1);
    }

    int i;
    double number;
    char tab='\t';
    char newline='\n';
    char *write = malloc(MAXCHAR);
    char buffer[MAXCHAR];

    for (i=0; i<(data->size); i++) {

        number = gsl_vector_get(data, i);

        snprintf(write, MAXCHAR, "%e", number);
        strncpy(buffer, write, MAXCHAR);
        Fwrite(buffer, sizeof(buffer), 1, fp);
        Fwrite(&tab, sizeof(tab), 1, fp);

    }
    //newline termination
    Fwrite(&newline, sizeof(newline), 1, fp);

    if ((fclose(fp)) != 0) {
        perror("fclose");
        exit(1);
    }

    free(write);
    return 0;
}

int
write_matrix_real(char *name, gsl_matrix *data) {
    FILE *fp;

    char file[MAXCHAR];

    if ((strncpy(file, DATA, MAXCHAR)) == NULL) {
        perror("strncpy");
        exit(1);
    }
    if ((strncat(file, name, MAXCHAR)) == NULL) {
        perror("strncat");
        exit(1);
    }
    if ((strncat(file, ".txt", MAXCHAR)) == NULL) {
        perror("strncat");
        exit(1);
    }

    //printf("directory %s \n", file_real);
    //printf("directory %s \n", file_imag);
    //open a file for write operation
    if ((fp = fopen(file, "w")) == NULL) {
        perror("fopen");
        exit(1);
    }

    int i, j;
    double number;
    char tab='\t';
    char newline='\n';
    char *write = malloc(MAXCHAR);
    char buffer[MAXCHAR];

    for (i=0; i<(data->size1); i++) {
        for (j=0; j<(data->size2); j++) {

            number = gsl_matrix_get(data, i, j);

            snprintf(write, MAXCHAR, "%e", number);
            strncpy(buffer, write, MAXCHAR);
            Fwrite(buffer, sizeof(buffer), 1, fp);
            Fwrite(&tab, sizeof(tab), 1, fp);
        }
        Fwrite(&newline, sizeof(newline), 1, fp);

    }
    //newline termination
    //Fwrite(&newline, sizeof(newline), 1, fp);

    if ((fclose(fp)) != 0) {
        perror("fclose");
        exit(1);
    }

    free(write);
    return 0;
}

int
read_vector_real(char *fname, gsl_vector *data) {
    FILE *fname_fp;
    if((fname_fp = fopen(fname, "r")) == NULL) {
        perror("fopen");
        exit(1);
    }

    int i=0;
    char line[MAXCHAR];
    while((fgets(line, MAXCHAR, fname_fp)) != NULL) {
        if (i==data->size) {
            fprintf(stderr, "file read mismatch - too many elements to read, aborting\n");
            exit(1);
        }
        line[strlen(line)-1] = '\0';
        gsl_vector_set(data, i, atof(line));
        i++;
    }

    if (i!=data->size) {
        fprintf(stderr, "file read mismatch - too few elements to read, aborting\n");
        exit(1);
    }

    if((fclose(fname_fp))) {
        perror("fclose");
    }

    return 0;

}

void
display_vector_real(gsl_vector *data) {
    int i;
    for (i=0; i<data->size; i++)
        printf("%e\t", gsl_vector_get(data, i));

    printf("\n");
}

//modifies parameter matrix
void
sqrtm(gsl_matrix_complex *matrix) {
    int i;
    //allocate workspace
    gsl_eigen_hermv_workspace *w = gsl_eigen_hermv_alloc (pow(2,NQ));
    //allocate eigenvalue vector
    gsl_vector *eval = gsl_vector_alloc(pow(2,NQ));
    //allocate eigenvector matrix - mutually orthogonal, normalized to 1
    gsl_matrix_complex *evec = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));

    if ((gsl_eigen_hermv (matrix, eval, evec, w)) != 0) {
        //could not find eigenvalues
        fprintf(stderr, "Failed at finding eigenvalues and eigenvectors\n");
        fprintf(stderr, "SQRTM failed, exiting...\n");
        exit(1);
    } else {
        if (VERBOSE)
            printf("SQRTM - succeeded in finding eigenvalues and eigenvectors\n");
    }

    //take square root of the eigenvalues, put them into a matrix
    //   for (i=0; i<pow(2,NQ); i++)
    //     gsl_vector_set(eval, i, sqrt(gsl_vector_get(eval, i)));

    gsl_eigen_hermv_sort(eval, evec, GSL_EIGEN_SORT_ABS_ASC);

    //take square root of the eigenvalues, put them into a matrix
    //matrix with eigenvalues across the diagonal
    gsl_matrix_complex *eval_matrix = gsl_matrix_complex_alloc(pow(2,NQ), pow(2, NQ));
    gsl_matrix_complex_set_zero(eval_matrix);
    //display_vector_real(eval);
    for (i=0; i<pow(2,NQ); i++) {
        gsl_matrix_complex_set(eval_matrix, i, i, gsl_complex_sqrt_real(gsl_vector_get(eval, i)));
    }
    //display_matrix(eval_matrix);

    //perform V*D^1/2*inv(V)
    gsl_matrix_complex *tmp = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, evec, eval_matrix, zero, tmp);
    gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one, tmp, evec, zero, matrix);

    //re-normalize the trace
    //gsl_matrix_complex_scale(density_quick, gsl_complex_rect(1.0/GSL_REAL(trace(density_quick)), 0));

    //deallocate memory
    gsl_matrix_complex_free(eval_matrix);
    gsl_eigen_hermv_free(w);
    gsl_vector_free(eval);
    gsl_matrix_complex_free(evec);
    gsl_matrix_complex_free(tmp);
}

//does not modify 2 parameter matrices
double
fidelity(gsl_matrix_complex *m1, gsl_matrix_complex *m2) {

    if ((m1->size1 != m2->size1) || (m1->size2 != m2->size2)) {
        fprintf(stderr, "inconsistent matrix size on fidelity()\n");
        exit(1);
    }

    gsl_matrix_complex *rho1 = gsl_matrix_complex_alloc(m1->size1, m1->size2);
    gsl_matrix_complex *rho2 = gsl_matrix_complex_alloc(m2->size1, m2->size2);
    gsl_matrix_complex *tmp = gsl_matrix_complex_alloc(m2->size1, m2->size2);
    gsl_matrix_complex_memcpy(rho1, m1);
    gsl_matrix_complex_memcpy(rho2, m2);

    double f = 0.0;

    sqrtm(rho1);

    gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, rho1, rho2, zero, tmp);
    gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, tmp, rho1, zero, rho2);

    //rho2 contains the answer
    sqrtm(rho2);
    f = GSL_REAL(trace(rho2));
    f = pow(f, 2);

    gsl_matrix_complex_free(rho1);
    gsl_matrix_complex_free(rho2);
    gsl_matrix_complex_free(tmp);

    return f;
}

double
entropy_linear(gsl_matrix_complex *matrix) {

    double s = 0.0;

    gsl_matrix_complex *tmp = gsl_matrix_complex_alloc(matrix->size1, matrix->size2);
    gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, matrix, matrix, zero, tmp);

    s = (1.0-GSL_REAL(trace(tmp)))*pow(2, NQ)/(pow(2,NQ) - 1);

    gsl_matrix_complex_free(tmp);

    return s;
}

//lower bound on concurrence
//does not modify original matrix
double
concurrence(gsl_matrix_complex *matrix) {

    if ((matrix->size1) != (matrix->size2)) {
        fprintf(stderr, "CONCURRENCE - must be a square matrix\n");
        exit(1);
    }

    int i, j;
    //allocate workspace
    gsl_eigen_herm_workspace *w = gsl_eigen_herm_alloc (pow(2,NQ));
    //allocate eigenvalue vector
    gsl_vector *eval = gsl_vector_alloc(pow(2,NQ));

    //create a stunt double - leave "matrix" intact
    gsl_matrix_complex *tmp = gsl_matrix_complex_alloc(matrix->size1, matrix->size2);
    gsl_matrix_complex_memcpy(tmp, matrix);

    gsl_matrix_complex *before, *after;
    after = gsl_matrix_complex_alloc(2, 2);
    before = gsl_matrix_complex_alloc(2, 2);
    gsl_matrix_complex_memcpy(after, sig_array[2]);
    gsl_matrix_complex_memcpy(before, sig_array[2]);

    for (i=0; i<NQ-1; i++) {
        gsl_matrix_complex_free(after);
        after = gsl_matrix_complex_alloc((before->size1)*2, (before->size2)*2);
        kron(before, sig_array[2], after);
        gsl_matrix_complex_free(before);
        before = gsl_matrix_complex_alloc(after->size1, after->size2);
        gsl_matrix_complex_memcpy(before, after);
    }

    //display_matrix(after);

    gsl_matrix_complex *conj = gsl_matrix_complex_alloc(matrix->size1, matrix->size2);
    for (i=0; i<matrix->size1; i++)
        for (j=0; j<matrix->size1; j++)
            gsl_matrix_complex_set(conj, i, j, gsl_complex_conjugate(gsl_matrix_complex_get(matrix, i, j)));

    /////////////////////////
    // IMPORTANT:
    //   You're right, it's not Hermitian, but there is a fix:  use the eigenvalues of
    //   R =   \sqrt{\rho} \Lamba \rho* \Lamba \sqrt{\rho}
    //   which is Hermitian.
    /////////////////////////

    //after is the SPIN-FLIP matrix - reusing memory space of before and tmp
    gsl_matrix_complex *matrix_copy = gsl_matrix_complex_alloc(matrix->size1, matrix->size2);
    gsl_matrix_complex_memcpy(matrix_copy, matrix);
    sqrtm(matrix_copy);
    gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, matrix_copy, after, zero, tmp);
    gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, tmp, conj, zero, before);
    gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, before, after, zero, tmp);
    gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, tmp, matrix_copy, zero, before);
    gsl_matrix_complex_free(conj);
    gsl_matrix_complex_free(matrix_copy);

    //display_matrix(before);

    if ((gsl_eigen_herm(before, eval, w)) != 0) {
        //could not find eigenvalues
        fprintf(stderr, "Failed at finding eigenvalues\n");
        fprintf(stderr, "CONCURRENCE failed, exiting...\n");
        exit(1);
    } else {
        if (VERBOSE)
            printf("CONCURRENCE - succeeded in finding eigenvalues\n");
    }

    //display_matrix(tmp);
    //display_vector_real(eval);
    //take the square roots of eigenvalues
    for (i=0; i<matrix->size1; i++) {
        //take the square root
        gsl_vector_set(eval, i, GSL_REAL(gsl_complex_sqrt_real(gsl_vector_get(eval, i))));
    }

    //display_vector_real(eval);
    //sort in ascending order
    //gsl_eigen_herm_sort(eval, evec, GSL_EIGEN_SORT_ABS_ASC);
    gsl_sort_vector(eval);

    double sum = 0.0;

    for (i=0; i<(matrix->size1)-1; i++) {
        sum += gsl_vector_get(eval, i);
    }

    double largest_eval = gsl_vector_get(eval, (matrix->size1)-1);

    //deallocate memory
    gsl_eigen_herm_free(w);
    gsl_vector_free(eval);
    gsl_matrix_complex_free(tmp);
    gsl_matrix_complex_free(before);
    gsl_matrix_complex_free(after);

    //return max of zero and eigenvalue difference
    if ((largest_eval - sum) < 0.0) {
        return 0.0;
    } else {
        return largest_eval - sum;
    }

}

void
gsl_vector_add_element(gsl_vector *vector, double element) {

    //append to existing elements
    int i, orig_size = vector->size;
    printf("Alloc tmp\n");
    gsl_vector *tmp = gsl_vector_alloc(vector->size);
    printf("Alloc OK\n");
    for (i=0; i<vector->size; i++) {
        printf("Set tmp\n");
        gsl_vector_set(tmp, i, gsl_vector_get(vector, i));
        printf("Set OK\n");
    }

    gsl_vector_free(vector);
    printf("Alloc vector\n");
    vector = gsl_vector_alloc(orig_size + 1);
    printf("Alloc OK\n");
    for (i=0; i<tmp->size; i++) {
        printf("Set vector\n");
        gsl_vector_set(vector, i, gsl_vector_get(tmp, i));
        printf("Set OK\n");
    }


    printf("Set vector 2\n");
    gsl_vector_set(vector, orig_size, element);
    printf("Set OK\n");

    gsl_vector_free(tmp);
}

void
PSI_matrix(gsl_matrix_complex *psi, double theta) {
    gsl_matrix_complex_set_zero(psi);
    double temp = 1-pow(theta, 2);
    gsl_matrix_complex_set(psi, 0, 0, gsl_complex_rect(temp, 0.0));
    gsl_matrix_complex_set(psi, 0, (psi->size2)-1, gsl_complex_rect(theta*sqrt(temp), 0.0));
    gsl_matrix_complex_set(psi, (psi->size1)-1, 0, gsl_complex_rect(theta*sqrt(temp), 0.0));
    gsl_matrix_complex_set(psi, (psi->size1)-1, (psi->size2)-1, gsl_complex_rect(pow(theta, 2), 0.0));
}

//returns absolute value of a vector
double
abs_val(gsl_vector *vec) {

    int i;
    double ret = 0.0;
    for (i=0; i < vec->size; i++) {
        ret += gsl_pow_2(gsl_vector_get(vec, i));
    }

    return sqrt(ret);
}

//MEM State - only defined for 2 qubits
void
rho_MEMS(gsl_matrix_complex *rho_mems, double gamma) {

    double g;

    if (gamma >= 2.0/3.0)
        g = (double) gamma/2.0;
    else
        g = (double) 1.0/3.0;

    gsl_matrix_complex_set_zero(rho_mems);
    gsl_matrix_complex_set(rho_mems, 0, 0, gsl_complex_rect(g, 0.0));
    gsl_matrix_complex_set(rho_mems, 0, 3, gsl_complex_rect((double) gamma/2.0, 0.0));
    gsl_matrix_complex_set(rho_mems, 1, 1, gsl_complex_rect((double) (1.0 - 2.0*g), 0.0));
    gsl_matrix_complex_set(rho_mems, 3, 0, gsl_complex_rect((double) gamma/2.0, 0.0));
    gsl_matrix_complex_set(rho_mems, 3, 3, gsl_complex_rect(g, 0.0));
}

void
W_state(gsl_matrix_complex *density) {
    int size = density->size1;
    int i, j;
    int index_array[NQ];
    gsl_matrix_complex *a, *b, *c;
    gsl_vector_complex *psi = gsl_vector_complex_alloc(size);
    for (j=0; j<NQ; j++)
        index_array[j]=0;

    gsl_vector_complex_set_zero(psi);
    b = gsl_matrix_complex_alloc(2, 1);
    //temporary allocation for c
    c = gsl_matrix_complex_alloc(1, 1);
    //compiler complains, init a
    a = gsl_matrix_complex_alloc(1,1);
    gsl_matrix_complex_free(a);

    for (i=0; i<NQ; i++) {
        //reset the index array
        for (j=0; j<NQ; j++)
            index_array[j]=0;
        index_array[i] = 1;

        //create W state component vector
        for (j=1; j<NQ; j++) {
            if (j==1) {
                a = gsl_matrix_complex_alloc(2, 1);
                if (index_array[0] == 0) {
                    gsl_matrix_complex_set(a, 0, 0, gsl_complex_rect(1, 0));
                    gsl_matrix_complex_set(a, 1, 0, gsl_complex_rect(0, 0));
                } else {
                    gsl_matrix_complex_set(a, 0, 0, gsl_complex_rect(0, 0));
                    gsl_matrix_complex_set(a, 1, 0, gsl_complex_rect(1, 0));
                }
            } else {
                gsl_matrix_complex_free(a);
                a = gsl_matrix_complex_alloc(c->size1, 1);
                gsl_matrix_complex_memcpy(a, c);
            }

            if (index_array[j] == 0) {
                gsl_matrix_complex_set(b, 0, 0, gsl_complex_rect(1, 0));
                gsl_matrix_complex_set(b, 1, 0, gsl_complex_rect(0, 0));
            } else {
                gsl_matrix_complex_set(b, 0, 0, gsl_complex_rect(0, 0));
                gsl_matrix_complex_set(b, 1, 0, gsl_complex_rect(1, 0));
            }
            gsl_matrix_complex_free(c);
            c = gsl_matrix_complex_alloc(2*(a->size1), 1);
            kron(a, b, c);
        }

        //psi + W state component vector
        for (j=0; j<size; j++)
            gsl_vector_complex_set(psi, j,
                                   gsl_complex_add(gsl_vector_complex_get(psi,j), gsl_matrix_complex_get(c, j, 0)));


        gsl_matrix_complex_free(a);
    }

    outer_product(psi, density);

    gsl_matrix_complex_scale(density, gsl_complex_div(gsl_complex_rect(1,0), trace(density)));

    //gsl_matrix_complex_free(a);
    gsl_matrix_complex_free(b);
    gsl_matrix_complex_free(c);
    gsl_vector_complex_free(psi);

}

// int
// write_data_binary(char *name, gsl_matrix_complex *data) {
//   FILE *fp;
//
//   char file[MAXCHAR];
//
//   if ((strncpy(file, DATA, MAXCHAR)) == NULL) {
//     perror("strncpy");
//     exit(1);
//   }
//   if ((strncat(file, name, MAXCHAR)) == NULL) {
//     perror("strncat");
//     exit(1);
//   }
//   if ((strncat(file, ".bin", MAXCHAR)) == NULL) {
//     perror("strncat");
//     exit(1);
//   }
//
//   //open a file for write operation
//   if ((fp = fopen(file, "w")) == NULL) {
//     perror("fopen");
//     exit(1);
//   }
//
//   gsl_matrix_complex_fwrite(fp, data);
//
//   if ((fclose(fp)) != 0) {
//     perror("fclose");
//     exit(1);
//   }
//
//   return 0;
// }

---