common.cpp

Go to the documentation of this file.

/*
Created on Fri Jun 26 14:13:26 2020
Copyright (C) 2020 Peter Rakyta, Ph.D.

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, either version 3 of the License, or
(at your option) any later version.

This program 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 this program.  If not, see http://www.gnu.org/licenses/.

@author: Peter Rakyta, Ph.D.
*/
//
// @brief A base class responsible for constructing matrices of C-NOT, U3
// gates acting on the N-qubit space

#include "qgd/common.h"



void* qgd_calloc( size_t element_num, size_t size, size_t alignment ) {

    void* ret = scalable_aligned_malloc( size*element_num, CACHELINE);
    memset(ret, 0, element_num*size);
    return ret;
}

void* qgd_realloc(void* aligned_ptr, size_t element_num, size_t size, size_t alignment ) {

    void* ret = scalable_aligned_realloc(aligned_ptr, size*element_num, CACHELINE);
    return ret;

}



void qgd_free( void* ptr ) {

    // preventing double free corruption
    if (ptr != NULL ) {
        scalable_aligned_free(ptr);
    }
    ptr = NULL;
}


int Power_of_2(int n) {
  if (n == 0) return 1;
  if (n == 1) return 2;

  return 2 * Power_of_2(n-1);
}



void add_unique_elelement( std::vector<int> &involved_qbits, int qbit ) {

    if ( involved_qbits.size() == 0 ) {
        involved_qbits.push_back( qbit );
    }

    for(std::vector<int>::iterator it = involved_qbits.begin(); it != involved_qbits.end(); ++it) {

        int current_val = *it;

        if (current_val == qbit) {
            return;
        }
        else if (current_val > qbit) {
            involved_qbits.insert( it, qbit );
            return;
        }

    }

    // add the qbit to the end if neither of the conditions were satisfied
    involved_qbits.push_back( qbit );

    return;

}


Matrix create_identity( int matrix_size ) {

    Matrix mtx = Matrix(matrix_size, matrix_size);
    memset(mtx.get_data(), 0, mtx.size()*sizeof(QGD_Complex16) );

    // setting the diagonal elelments to identity
    for(int idx = 0; idx < matrix_size; ++idx)
    {
        int element_index = idx*matrix_size + idx;
            mtx[element_index].real = 1.0;
    }

    return mtx;

}



int create_identity( QGD_Complex16* matrix, int matrix_size ) {

    memset( matrix, 0, matrix_size*matrix_size*sizeof(QGD_Complex16) );

    // setting the diagonal elelments to identity
    for(int idx = 0; idx < matrix_size; ++idx)
    {
        int element_index = idx*matrix_size + idx;
            matrix[element_index].real = 1.0;
    }

    return 0;

}


QGD_Complex16 scalar_product( QGD_Complex16* A, QGD_Complex16* B, int vector_size) {

    // parameters alpha and beta for the cblas_zgemm3m function
    QGD_Complex16 alpha;
    alpha.real = 1.0;
    alpha.imag = 0.0;
    QGD_Complex16 beta;
    beta.real = 0.0;
    beta.imag = 0.0;

    // preallocate array for the result
    QGD_Complex16 C;

    // calculate the product of A and B
#ifdef CBLAS
    cblas_zgemm3m (CblasRowMajor, CblasNoTrans, CblasConjTrans, 1, 1, vector_size, &alpha, A, vector_size, B, vector_size, &beta, &C, 1);
#else
    cblas_zgemm (CblasRowMajor, CblasNoTrans, CblasConjTrans, 1, 1, vector_size, &alpha, A, vector_size, B, vector_size, &beta, &C, 1);
#endif

    return C;
}



int zgemm3m_wrapper_adj( QGD_Complex16* A, QGD_Complex16* B, QGD_Complex16* C, int matrix_size) {

    // parameters alpha and beta for the cblas_zgemm3m function
    QGD_Complex16 alpha;
    alpha.real = 1.0;
    alpha.imag = 0.0;
    QGD_Complex16 beta;
    beta.real = 0.0;
    beta.imag = 0.0;

    // calculate the product of A and B
#ifdef CBLAS
    cblas_zgemm3m (CblasRowMajor, CblasNoTrans, CblasConjTrans, matrix_size, matrix_size, matrix_size, &alpha, A, matrix_size, B, matrix_size, &beta, C, matrix_size);
#else
    cblas_zgemm (CblasRowMajor, CblasNoTrans, CblasConjTrans, matrix_size, matrix_size, matrix_size, &alpha, A, matrix_size, B, matrix_size, &beta, C, matrix_size);
#endif

    return 0;
}



QGD_Complex16* zgemm3m_wrapper( QGD_Complex16* A, QGD_Complex16* B, int matrix_size) {

    // parameters alpha and beta for the cblas_zgemm3m function
    QGD_Complex16 alpha;
    alpha.real = 1.0;
    alpha.imag = 0.0;
    QGD_Complex16 beta;
    beta.real = 0.0;
    beta.imag = 0.0;

    // preallocate array for the result
    QGD_Complex16* C = (QGD_Complex16*)qgd_calloc(matrix_size*matrix_size, sizeof(QGD_Complex16), 64);

    // calculate the product of A and B
#ifdef CBLAS
    cblas_zgemm3m (CblasRowMajor, CblasNoTrans, CblasNoTrans, matrix_size, matrix_size, matrix_size, &alpha, A, matrix_size, B, matrix_size, &beta, C, matrix_size);
#else
    cblas_zgemm (CblasRowMajor, CblasNoTrans, CblasNoTrans, matrix_size, matrix_size, matrix_size, &alpha, A, matrix_size, B, matrix_size, &beta, C, matrix_size);
#endif

    return C;
}


int zgemm3m_wrapper( QGD_Complex16* A, QGD_Complex16* B, QGD_Complex16* C, int matrix_size) {

    // parameters alpha and beta for the cblas_zgemm3m function
    QGD_Complex16 alpha;
    alpha.real = 1.0;
    alpha.imag = 0.0;
    QGD_Complex16 beta;
    beta.real = 0.0;
    beta.imag = 0.0;

    // calculate the product of A and B
#ifdef CBLAS
    cblas_zgemm3m (CblasRowMajor, CblasNoTrans, CblasNoTrans, matrix_size, matrix_size, matrix_size, &alpha, A, matrix_size, B, matrix_size, &beta, C, matrix_size);
#else
    cblas_zgemm (CblasRowMajor, CblasNoTrans, CblasNoTrans, matrix_size, matrix_size, matrix_size, &alpha, A, matrix_size, B, matrix_size, &beta, C, matrix_size);
#endif

    return 0;
}



Matrix
reduce_zgemm( std::vector<Matrix>& mtxs ) {



    if (mtxs.size() == 0 ) {
        return Matrix();
    }

    // pointers to matrices to be used in the multiplications
    Matrix A;
    Matrix B;
    Matrix C;

    // the iteration number
    int iteration = 0;


    A = *mtxs.begin();

    // calculate the product of complex matrices
    for(std::vector<Matrix>::iterator it=++mtxs.begin(); it != mtxs.end(); ++it) {

        iteration++;
        B = *it;

        if ( iteration>1 ) {
            A = C;
        }

        // calculate the product of A and B
        C = dot(A, B);
/*
A.print_matrix();
B.print_matrix();
C.print_matrix();
*/
    }

    //C.print_matrix();
    return C;

}


void subtract_diag( Matrix& mtx, QGD_Complex16 scalar ) {

    size_t matrix_size = mtx.rows;

    for(size_t idx = 0; idx < matrix_size; idx++)   {
        size_t element_idx = idx*matrix_size+idx;
        mtx[element_idx].real = mtx[element_idx].real - scalar.real;
        mtx[element_idx].imag = mtx[element_idx].imag - scalar.imag;
    }

}




QGD_Complex16 mult( QGD_Complex16 a, QGD_Complex16 b ) {

    QGD_Complex16 ret;
    ret.real = a.real*b.real - a.imag*b.imag;
    ret.imag = a.real*b.imag + a.imag*b.real;

    return ret;

}

QGD_Complex16 mult( double a, QGD_Complex16 b ) {

    QGD_Complex16 ret;
    ret.real = a*b.real;
    ret.imag = a*b.imag;

    return ret;

}



void mult( QGD_Complex16 a, Matrix& b ) {

    size_t element_num = b.size();

    for (size_t idx=0; idx<element_num; idx++) {
        QGD_Complex16 tmp = b[idx];
        b[idx].real = a.real*tmp.real - a.imag*tmp.imag;
        b[idx].imag = a.real*tmp.imag + a.imag*tmp.real;
    }

    return;

}