mat_ops.hpp

Go to the documentation of this file.

#ifndef QNLP_MAT_OPS
#define QNLP_MAT_OPS

#include <complex>
#include <vector>
#include <iostream>

namespace QNLP{
    template <class Mat2x2Type>            
    const Mat2x2Type matrixSqrt(const Mat2x2Type& U){
        Mat2x2Type V(U);
        std::complex<double> delta = U(0,0)*U(1,1) - U(0,1)*U(1,0);
        std::complex<double> tau = U(0,0) + U(1,1);
        std::complex<double> s = sqrt(delta);
        std::complex<double> t = sqrt(tau + 2.0*s);

        //must be a way to vectorise these; TinyMatrix have a scale/shift option?
        V(0,0) += s;
        V(1,1) += s;
        std::complex<double> scale_factor(1.,0.);
        scale_factor /= t;
        V(0,0) *= scale_factor; //(std::complex<double>(1.,0.)/t);
        V(0,1) *= scale_factor; //(1/t);
        V(1,0) *= scale_factor; //(1/t);
        V(1,1) *= scale_factor; //(1/t);

        return V;
    }

    template <class Mat2x2Type>            
    Mat2x2Type adjointMatrix(const Mat2x2Type& U){
        Mat2x2Type Uadjoint(U);
        std::complex<double> tmp;
        tmp = Uadjoint(0,1);
        Uadjoint(0,1) = Uadjoint(1,0);
        Uadjoint(1,0) = tmp;
        Uadjoint(0,0) = std::conj(Uadjoint(0,0));
        Uadjoint(0,1) = std::conj(Uadjoint(0,1));
        Uadjoint(1,0) = std::conj(Uadjoint(1,0));
        Uadjoint(1,1) = std::conj(Uadjoint(1,1));
        return Uadjoint;
    }
}
#endif