#!/usr/bin/python3

import os,sys,math
import numpy as np

cudadec = "__host__ __device__"

def typechar(dtype):
    
    sshort = ""
    #return short typechar given typestring
    match(dtype):
        case "double":
            sshort = ""
        case "float":
            sshort = "f"
        case "int":
            sshort = "i"

    return sshort

def vtypename(dim,dtype):
    name = "cuvec{}{}".format(dim,typechar(dtype))
    return name

def mtypename(dim,dtype):
    name = "cumat{}{}".format(dim,typechar(dtype))
    return name

def vdc(dim):
    dcomp = "x"
    match(dim):
        case 0:
            dcomp = "x"
        case 1:
            dcomp = "y"
        case 2:
            dcomp = "z"
        case 3:
            dcomp = "w"

    return dcomp


def genvecheaderv1(dim,dtype):

    lns = ""
    name = vtypename(dim,dtype)
    lns += "class {}\n".format(name)
    lns += "{\n\tpublic:\n"

    for I in range(0,dim):
        lns += "\t{} {};\n".format(dtype,vdc(I))

    lns += "\n"
    lns += "\t{} {}();\n".format(cudadec,name)
    lns += "\t{} ~{}();\n".format(cudadec,name)
    lns += "\t{} {}(".format(cudadec,name)
    for I in range(0,dim):
        lns += "const {} &_{}".format(dtype,vdc(I))
        if(I<dim-1):
            lns +=", "
    lns +=");\n\n"


    

    lns += "};"

    return lns

def fwrite(fname,lns):
    fp = open(fname,"w+")
    fp.writelines(lns)
    fp.close()

def vconstructor1(dim,dtype):
    lnsh = ""
    lnss = ""

    vname = vtypename(dim,dtype)

    lnsh += "\t{} {}(); \n".format(cudadec,vname)

    lnss += "\t{} {}::{}() \n".format(cudadec,vname,vname)
    lnss += "\t{\n"

    lnss += "\t\t"
    for I in range(0,dim):
        lnss += "{} = 0; ".format(vdc(I))
    lnss += "\n"
    

    lnss += "\t\treturn;\n"
    lnss += "\t}\n\n"
    
    return [lnsh,lnss]

def vdestructor1(dim,dtype):
    lnsh = ""
    lnss = ""

    vname = vtypename(dim,dtype)

    lnsh += "\t{} ~{}(); \n".format(cudadec,vname)

    lnss += "\t{} {}::~{}() \n".format(cudadec,vname,vname)
    lnss += "\t{\n"

    lnss += "\t\t"
    for I in range(0,dim):
        lnss += "{} = 0; ".format(vdc(I))
    lnss += "\n"
    

    lnss += "\t\treturn;\n"
    lnss += "\t}\n\n"
    
    return [lnsh,lnss]

def vconstructor2(dim,dtype):
    lnsh = ""
    lnss = ""

    vname = vtypename(dim,dtype)

    lnsh += "\t{} {}(".format(cudadec,vname)
    for I in range(0,dim):
        lnsh += "const {} &_{}".format(dtype,vdc(I))
        if(I<dim-1):
            lnsh += ", "
    lnsh += "); \n"

    lnss += "\t{} {}::{}(".format(cudadec,vname,vname)
    for I in range(0,dim):
        lnss += "const {} &_{}".format(dtype,vdc(I))
        if(I<dim-1):
            lnss += ", "
    lnss += ") \n"
    lnss += "\t{\n"

    lnss += "\t\t"
    for I in range(0,dim):
        lnss += "{} = _{}; ".format(vdc(I),vdc(I))
    lnss += "\n"
    

    lnss += "\t\treturn;\n"
    lnss += "\t}\n\n"
    
    return [lnsh,lnss]

def v_operatorbk1(dim,dtype):
    lnsh = ""
    lnss = ""
    vname = vtypename(dim,dtype)

    lnsh += f"\t{cudadec} {dtype}& operator[](const int &I);\n"

    lnss += f"\t{cudadec} {dtype}& {vname}::operator[](const int &I)\n"
    lnss += "\t{\n"
    lnss += "\t\tswitch(I)\n\t\t{\n"
    for I in range(0,dim):
        lnss += "\t\t\tcase {}:\n".format(I)
        lnss += "\t\t\t\treturn {};\n".format(vdc(I))
    lnss += "\t\t}\n"
    lnss += "\n\t\treturn {};\n".format(vdc(0))
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def v_operatorbk2(dim,dtype):
    lnsh = ""
    lnss = ""
    vname = vtypename(dim,dtype)

    lnsh += f"\t{cudadec} const {dtype}& operator[](const int &I) const;\n"

    lnss += f"\t{cudadec} const {dtype}& {vname}::operator[](const int &I) const\n"
    lnss += "\t{\n"
    lnss += "\t\tswitch(I)\n\t\t{\n"
    for I in range(0,dim):
        lnss += "\t\t\tcase {}:\n".format(I)
        lnss += "\t\t\t\treturn {};\n".format(vdc(I))
    lnss += "\t\t}\n"
    lnss += "\n\t\treturn {};\n".format(vdc(0))
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def v_operatorplus1(dim,dtype):
    lnsh = ""
    lnss = ""
    vname = vtypename(dim,dtype)

    lnsh += f"\t{cudadec} {vname} operator+(const {vname}& rhs) const;\n"

    lnss += f"\t{cudadec} {vname} {vname}::operator+(const {vname}& rhs) const\n"
    lnss += "\t{\n"
    lnss += f"\t\t{vname} ret;\n"
    for I in range(0,dim):
        lnss += f"\t\tret.{vdc(I)} = {vdc(I)} + rhs.{vdc(I)};\n"
    lnss += "\t\treturn ret;\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def v_operatorminus1(dim,dtype):
    lnsh = ""
    lnss = ""
    vname = vtypename(dim,dtype)

    lnsh += f"\t{cudadec} {vname} operator-(const {vname}& rhs) const;\n"

    lnss += f"\t{cudadec} {vname} {vname}::operator-(const {vname}& rhs) const\n"
    lnss += "\t{\n"
    lnss += f"\t\t{vname} ret;\n"
    for I in range(0,dim):
        lnss += f"\t\tret.{vdc(I)} = {vdc(I)} - rhs.{vdc(I)};\n"
    lnss += "\t\treturn ret;\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def v_operatortimes3(dim,dtype):
    lnsh = ""
    lnss = ""
    vname = vtypename(dim,dtype)

    lnsh += f"\t{cudadec} {vname} operator*(const {vname}& rhs) const; \t //elementwise product\n"

    lnss += f"\t{cudadec} {vname} {vname}::operator*(const {vname}& rhs) const\n"
    lnss += "\t{\n"
    lnss += "\t\t//Elementwise product\n"
    lnss += f"\t\t{vname} ret;\n"
    for I in range(0,dim):
        lnss += f"\t\tret.{vdc(I)} = {vdc(I)} * rhs.{vdc(I)};\n"
    lnss += "\t\treturn ret;\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def v_operatordiv3(dim,dtype):
    lnsh = ""
    lnss = ""
    vname = vtypename(dim,dtype)

    lnsh += f"\t{cudadec} {vname} operator/(const {vname}& rhs) const; \t //elementwise division\n"

    lnss += f"\t{cudadec} {vname} {vname}::operator/(const {vname}& rhs) const\n"
    lnss += "\t{\n"
    lnss += "\t\t//Elementwise division\n"
    lnss += f"\t\t{vname} ret;\n"
    for I in range(0,dim):
        lnss += f"\t\tret.{vdc(I)} = {vdc(I)} / rhs.{vdc(I)};\n"
    lnss += "\t\treturn ret;\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def v_operatortimes1(dim,dtype):
    lnsh = ""
    lnss = ""
    vname = vtypename(dim,dtype)

    lnsh += f"\t{cudadec} friend {vname} operator*(const {vname}& lhs, const {dtype}& rhs);\n"

    lnss += f"\t{cudadec} {vname} operator*(const {vname}& lhs, const {dtype}& rhs)\n"
    lnss += "\t{\n"
    lnss += f"\t\t{vname} ret;\n"
    for I in range(0,dim):
        lnss += f"\t\tret.{vdc(I)} = lhs.{vdc(I)}*rhs;\n"
    lnss += "\t\treturn ret;\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def v_operatortimes2(dim,dtype):
    lnsh = ""
    lnss = ""
    vname = vtypename(dim,dtype)

    lnsh += f"\t{cudadec} friend {vname} operator*(const {dtype}& lhs, const {vname}& rhs);\n"

    lnss += f"\t{cudadec} {vname} operator*(const {dtype}& lhs, const {vname}& rhs)\n"
    lnss += "\t{\n"
    lnss += f"\t\t{vname} ret;\n"
    for I in range(0,dim):
        lnss += f"\t\tret.{vdc(I)} = lhs*rhs.{vdc(I)};\n"
    lnss += "\t\treturn ret;\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def v_operatordiv1(dim,dtype):
    lnsh = ""
    lnss = ""
    vname = vtypename(dim,dtype)

    lnsh += f"\t{cudadec} friend {vname} operator/(const {vname}& lhs, const {dtype}& rhs);\n"

    lnss += f"\t{cudadec} {vname} operator/(const {vname}& lhs, const {dtype}& rhs)\n"
    lnss += "\t{\n"
    lnss += f"\t\t{vname} ret;\n"
    for I in range(0,dim):
        lnss += f"\t\tret.{vdc(I)} = lhs.{vdc(I)}/rhs;\n"
    lnss += "\t\treturn ret;\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def v_operatordiv2(dim,dtype):
    lnsh = ""
    lnss = ""
    vname = vtypename(dim,dtype)

    lnsh += f"\t{cudadec} friend {vname} operator/(const {dtype}& lhs, const {vname}& rhs);\n"

    lnss += f"\t{cudadec} {vname} operator/(const {dtype}& lhs, const {vname}& rhs)\n"
    lnss += "\t{\n"
    lnss += f"\t\t{vname} ret;\n"
    for I in range(0,dim):
        lnss += f"\t\tret.{vdc(I)} = lhs/rhs.{vdc(I)};\n"
    lnss += "\t\treturn ret;\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def v_operatorneg(dim,dtype):
    lnsh = ""
    lnss = ""
    vname = vtypename(dim,dtype)

    lnsh += f"\t{cudadec} friend {vname} operator-(const {vname}& other);\n"

    lnss += f"\t{cudadec} {vname} operator-(const {vname}& other)\n"
    lnss += "\t{\n"
    lnss += f"\t\t{vname} ret;\n"
    for I in range(0,dim):
        lnss += f"\t\tret.{vdc(I)} = -other.{vdc(I)};\n"
    lnss += "\t\treturn ret;\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def v_operatorpluseq1(dim,dtype):
    lnsh = ""
    lnss = ""
    vname = vtypename(dim,dtype)

    lnsh += f"\t{cudadec} {vname}& operator+=(const {vname}& rhs);\n"

    lnss += f"\t{cudadec} {vname}& {vname}::operator+=(const {vname}& rhs)\n"
    lnss += "\t{\n"
    for I in range(0,dim):
        lnss += f"\t\t{vdc(I)} += rhs.{vdc(I)};\n"
    lnss += "\t\treturn *this;\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def v_operatorminuseq1(dim,dtype):
    lnsh = ""
    lnss = ""
    vname = vtypename(dim,dtype)

    lnsh += f"\t{cudadec} {vname}& operator-=(const {vname}& rhs);\n"

    lnss += f"\t{cudadec} {vname}& {vname}::operator-=(const {vname}& rhs)\n"
    lnss += "\t{\n"
    for I in range(0,dim):
        lnss += f"\t\t{vdc(I)} -= rhs.{vdc(I)};\n"
    lnss += "\t\treturn *this;\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def v_operatortimeseq1(dim,dtype):
    lnsh = ""
    lnss = ""
    vname = vtypename(dim,dtype)

    lnsh += f"\t{cudadec} {vname}& operator*=(const {dtype}& rhs);\n"

    lnss += f"\t{cudadec} {vname}& {vname}::operator*=(const {dtype}& rhs)\n"
    lnss += "\t{\n"
    for I in range(0,dim):
        lnss += f"\t\t{vdc(I)} *= rhs;\n"
    lnss += "\t\treturn *this;\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]


def v_operatordiveq1(dim,dtype):
    lnsh = ""
    lnss = ""
    vname = vtypename(dim,dtype)

    lnsh += f"\t{cudadec} {vname}& operator/=(const {dtype}& rhs);\n"

    lnss += f"\t{cudadec} {vname}& {vname}::operator/=(const {dtype}& rhs)\n"
    lnss += "\t{\n"
    for I in range(0,dim):
        lnss += f"\t\t{vdc(I)} /= rhs;\n"
    lnss += "\t\treturn *this;\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def genvectorstuff(dim,dtype):

    lnsh = ""
    lnss = ""

    [lnsh2,lnss2] = vconstructor1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = vdestructor1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = vconstructor2(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = v_operatorbk1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = v_operatorbk2(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = v_operatorplus1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = v_operatorminus1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = v_operatortimes3(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = v_operatordiv3(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = v_operatortimes1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = v_operatortimes2(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = v_operatordiv1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = v_operatordiv2(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = v_operatorneg(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = v_operatorpluseq1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = v_operatorminuseq1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = v_operatortimeseq1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2
    [lnsh2,lnss2] = v_operatordiveq1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2

    return [lnsh,lnss]

def matrix_data(dim,dtype):
    lnsh = ""
    lnss = ""

    
    for J in range(0,dim):
        lnsh += f"{dtype} "
        for I in range(0,dim):
            lnsh += f"m{I}{J}"
            if(I<dim-1):
                lnsh += ","
        lnsh += ";\n"

    lnsh+="\n"

    return [lnsh, lnss]


def mconstructor1(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"\t{cudadec} {mname}(); \n"

    lnss += f"\t{cudadec} {mname}::{mname}() \n"
    lnss += "\t{\n"

    for I in range(0,dim):
        for J in range(0,dim):
            lnss += f"\t\tm{I}{J} = 0;\n"
        lnss += "\n"
    
    lnss += "\t\treturn;\n"
    lnss += "\t}\n\n"
    
    return [lnsh,lnss]

def mdestructor1(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"\t{cudadec} ~{mname}(); \n"

    lnss += f"\t{cudadec} {mname}::~{mname}() \n"
    lnss += "\t{\n"

    for I in range(0,dim):
        for J in range(0,dim):
            lnss += f"\t\t//m{I}{J} = 0;\n"
        lnss += "\n"
    

    lnss += "\t\treturn;\n"
    lnss += "\t}\n\n"
    
    return [lnsh,lnss]

def mconstructor2(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"\t{cudadec} {mname}(\n"
    for I in range(0,dim):
        lnsh+="\t\t"
        for J in range(0,dim):
            lnsh += f"const {dtype}& _m{I}{J}"
            if(J+I*dim<dim*dim-1):
                lnsh+=", "
        lnsh+="\n"
    lnsh += "\t);\n"

    lnss += f"\t{cudadec} {mname}::{mname}(\n"
    for I in range(0,dim):
        lnss+="\t\t"
        for J in range(0,dim):
            lnss += f"const {dtype}& _m{I}{J}"
            if(J+I*dim<dim*dim-1):
                lnss+=", "
        lnss+="\n"
    lnss += "\t)\n"
    lnss += "\t{\n"

    for I in range(0,dim):
        for J in range(0,dim):
            lnss += f"\t\tm{J}{I} = _m{J}{I}; \n"
        lnss+="\n"
    lnss += "\n"
    

    lnss += "\t\treturn;\n"
    lnss += "\t}\n\n"
    
    return [lnsh,lnss]

def mconstructor3(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"\t{cudadec} {mname}(const {dtype}* data{dim*dim});\n"

    lnss += f"\t{cudadec} {mname}::{mname}(const {dtype}* data{dim*dim})\n"
    lnss += "\t{\n"

    for I in range(0,dim):
        for J in range(0,dim):
            lnss += f"\t\tm{J}{I} = data{dim*dim}[{J+dim*I}]; \n"
        lnss+="\n"
    lnss += "\n"
    

    lnss += "\t\treturn;\n"
    lnss += "\t}\n\n"
    
    return [lnsh,lnss]

def moperatorbk1(dim,dtype):
    lnsh = ""
    lnss = ""
    mname = mtypename(dim,dtype)

    lnsh += f"\t{cudadec} {dtype}& operator[](const int &I);\n"

    lnss += f"\t{cudadec} {dtype}& {mname}::operator[](const int &I)\n"
    lnss += "\t{\n"
    lnss += "\t\tswitch(I)\n\t\t{\n"
    for I in range(0,dim*dim):
        Ix = I%dim
        Iy = I//dim
        lnss += "\t\t\tcase {}:\n".format(I)
        lnss += f"\t\t\t\treturn m{Ix}{Iy};\n"
    lnss += "\t\t}\n"
    lnss += "\n\t\treturn m00;\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def moperatorbk2(dim,dtype):
    lnsh = ""
    lnss = ""
    mname = mtypename(dim,dtype)

    lnsh += f"\t{cudadec} const {dtype}& operator[](const int &I) const;\n"

    lnss += f"\t{cudadec} const {dtype}& {mname}::operator[](const int &I) const\n"
    lnss += "\t{\n"
    lnss += "\t\tswitch(I)\n\t\t{\n"
    for I in range(0,dim*dim):
        Ix = I%dim
        Iy = I//dim
        lnss += "\t\t\tcase {}:\n".format(I)
        lnss += f"\t\t\t\treturn m{Ix}{Iy};\n"
    lnss += "\t\t}\n"
    lnss += "\n\t\treturn m00;\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def moperatorbk3(dim,dtype):
    lnsh = ""
    lnss = ""
    mname = mtypename(dim,dtype)

    lnsh += f"\t{cudadec} {dtype}& operator()(const int &I, const int &J);\n"

    lnss += f"\t{cudadec} {dtype}& {mname}::operator()(const int &I, const int &J)\n"
    lnss += "\t{\n"
    lnss += f"\t\treturn (*this)[I+{dim}*J];\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def moperatorbk4(dim,dtype):
    lnsh = ""
    lnss = ""
    mname = mtypename(dim,dtype)

    lnsh += f"\t{cudadec} const {dtype}& operator()(const int &I, const int &J) const;\n"

    lnss += f"\t{cudadec} const {dtype}& {mname}::operator()(const int &I, const int &J) const\n"
    lnss += "\t{\n"
    lnss += f"\t\treturn (*this)[I+{dim}*J];\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def moperatorbk5(dim,dtype):
    lnsh = ""
    lnss = ""
    mname = mtypename(dim,dtype)

    lnsh += f"\t{cudadec} {dtype}& at(const int &I, const int &J);\n"

    lnss += f"\t{cudadec} {dtype}& {mname}::at(const int &I, const int &J)\n"
    lnss += "\t{\n"
    lnss += f"\t\treturn (*this)[I+{dim}*J];\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

def moperatorbk6(dim,dtype):
    lnsh = ""
    lnss = ""
    mname = mtypename(dim,dtype)

    lnsh += f"\t{cudadec} const {dtype}& at(const int &I, const int &J) const;\n"

    lnss += f"\t{cudadec} const {dtype}& {mname}::at(const int &I, const int &J) const\n"
    lnss += "\t{\n"
    lnss += f"\t\treturn (*this)[I+{dim}*J];\n"
    lnss += "\t}\n\n"

    return [lnsh,lnss]

#matrix arithmetic operators

def mdataoperator1(dim,dtype):
    lnsh = ""
    lnss = ""
    mname = mtypename(dim,dtype)

    lnsh += f"{cudadec} {dtype}* data(); //pointer to {dtype}{dim*dim} representation of matrix\n"

    lnss += f"{cudadec} {dtype}* {mname}::data()\n"
    lnss += "{\n"
    lnss += f"\treturn ({dtype}*)this;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]

def mdataoperator2(dim,dtype):
    lnsh = ""
    lnss = ""
    mname = mtypename(dim,dtype)

    lnsh += f"{cudadec} const {dtype}* data() const; //pointer to {dtype}{dim*dim} representation of matrix\n"

    lnss += f"{cudadec} const {dtype}* {mname}::data() const\n"
    lnss += "{\n"
    lnss += f"\treturn ({dtype}*)this;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]


def moperatorplus1(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"{cudadec} {mname} operator+(const {mname}& rhs) const;\n"
    lnss += f"{cudadec} {mname} {mname}::operator+(const {mname}& rhs) const\n"
    lnss += "{\n"
    lnss += f"\t{mname} ret;\n"
    for I in range(0,dim):
        for J in range(0,dim):
            lnss += f"\tret.m{J}{I} = m{J}{I} + rhs.m{J}{I};\n"
        lnss += "\n"
    lnss += "\treturn ret;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]

def moperatorminus1(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"{cudadec} {mname} operator-(const {mname}& rhs) const;\n"
    lnss += f"{cudadec} {mname} {mname}::operator-(const {mname}& rhs) const\n"
    lnss += "{\n"
    lnss += f"\t{mname} ret;\n"
    for I in range(0,dim):
        for J in range(0,dim):
            lnss += f"\tret.m{J}{I} = m{J}{I} - rhs.m{J}{I};\n"
        lnss += "\n"
    lnss += "\treturn ret;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]

def moperatortimes1(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"{cudadec} {mname} operator*(const {mname}& rhs) const;\n"
    lnss += f"{cudadec} {mname} {mname}::operator*(const {mname}& rhs) const\n"
    lnss += "{\n"
    lnss += f"\t{mname} ret; //should be zeroed in constructor\n"
    for I in range(0,dim):
        for K in range(0,dim):
            lnss += f"\tret.m{I}{K} = "
            for J in range(0,dim):
                lnss += f"m{I}{J}*rhs.m{J}{K}"
                if(J<dim-1):
                    lnss += " + "
                else:
                    lnss += ";\n"
        lnss += "\n"
    lnss += "\treturn ret;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]

def moperatortimes2(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"{cudadec} friend {mname} operator*(const {mname}& lhs, const {dtype}& rhs);\n"
    lnss += f"{cudadec} {mname} operator*(const {mname}& lhs, const {dtype}& rhs)\n"
    lnss += "{\n"
    lnss += f"\t{mname} ret;\n"
    for I in range(0,dim):
        for J in range(0,dim):
            lnss += f"\tret.m{J}{I}=lhs.m{J}{I}*rhs;\n"
        #lnss += "\n"
    lnss += "\treturn ret;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]

def moperatordiv2(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"{cudadec} friend {mname} operator/(const {mname}& lhs, const {dtype}& rhs);\n"
    lnss += f"{cudadec} {mname} operator/(const {mname}& lhs, const {dtype}& rhs)\n"
    lnss += "{\n"
    lnss += f"\t{mname} ret;\n"
    for I in range(0,dim):
        for J in range(0,dim):
            lnss += f"\tret.m{J}{I}=lhs.m{J}{I}/rhs;\n"
        #lnss += "\n"
    lnss += "\treturn ret;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]

def moperatortimes3(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"{cudadec} friend {mname} operator*(const {dtype}& lhs, const {mname}& rhs);\n"
    lnss += f"{cudadec} {mname} operator*(const {dtype}& lhs, const {mname}& rhs)\n"
    lnss += "{\n"
    lnss += f"\t{mname} ret;\n"
    for I in range(0,dim):
        for J in range(0,dim):
            lnss += f"\tret.m{J}{I}=lhs*rhs.m{J}{I};\n"
        #lnss += "\n"
    lnss += "\treturn ret;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]

def moperatortimes4(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)
    vname = vtypename(dim,dtype)

    lnsh += f"{cudadec} friend {vname} operator*(const {mname}& lhs, const {vname}& rhs);\n"
    lnss += f"{cudadec} {vname} operator*(const {mname}& lhs, const {vname}& rhs)\n"
    lnss += "{\n"
    lnss += f"\t{vname} ret;\n"
    for I in range(0,dim):
        lnss += f"\tret.{vdc(I)} = "
        for J in range(0,dim):
            lnss += f"lhs.m{I}{J}*rhs.{vdc(J)}"
            if(J<dim-1):
                lnss += " + "
            else:
                lnss += ";\n"
        #lnss += "\n"
    lnss += "\treturn ret;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]

def moperatortimes5(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)
    vname = vtypename(dim,dtype)

    lnsh += f"{cudadec} friend {vname} operator*(const {vname}& lhs, const {mname}& rhs);\n"
    lnss += f"{cudadec} {vname} operator*(const {vname}& lhs, const {mname}& rhs)\n"
    lnss += "{\n"
    lnss += f"\t{vname} ret;\n"
    for I in range(0,dim):
        lnss += f"\tret.{vdc(I)} = "
        for J in range(0,dim):
            #lnss += f"lhs.m{I}{J}*rhs.{vdc(J)}"
            lnss += f"lhs.{vdc(J)}*rhs.m{J}{I}"
            if(J<dim-1):
                lnss += " + "
            else:
                lnss += ";\n"
        #lnss += "\n"
    lnss += "\treturn ret;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]

def moperatorneg(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"{cudadec} friend {mname} operator-(const {mname}& rhs);\n"
    lnss += f"{cudadec} {mname} operator-(const {mname}& rhs)\n"
    lnss += "{\n"
    lnss += f"\t{mname} ret;\n"
    for I in range(0,dim):
        for J in range(0,dim):
            lnss += f"\tret.m{J}{I} = -rhs.m{J}{I};\n"
        #lnss += "\n"
    lnss += "\treturn ret;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]

def moperatorpluseq(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"{cudadec} {mname}& operator+=(const {mname}& rhs);\n"
    lnss += f"{cudadec} {mname}& {mname}::operator+=(const {mname}& rhs)\n"
    lnss += "{\n"
    for I in range(0,dim):
        for J in range(0,dim):
            lnss += f"\tm{J}{I} += rhs.m{J}{I};\n"
        #lnss += "\n"
    lnss += "\treturn *this;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]

def moperatorminuseq(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"{cudadec} {mname}& operator-=(const {mname}& rhs);\n"
    lnss += f"{cudadec} {mname}& {mname}::operator-=(const {mname}& rhs)\n"
    lnss += "{\n"
    for I in range(0,dim):
        for J in range(0,dim):
            lnss += f"\tm{J}{I} -= rhs.m{J}{I};\n"
        #lnss += "\n"
    lnss += "\treturn *this;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]

def moperatortimeseq1(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"{cudadec} {mname}& operator*=(const {dtype}& rhs);\n"
    lnss += f"{cudadec} {mname}& {mname}::operator*=(const {dtype}& rhs)\n"
    lnss += "{\n"
    for I in range(0,dim):
        for J in range(0,dim):
            lnss += f"\tm{J}{I} *= rhs;\n"
        #lnss += "\n"
    lnss += "\treturn *this;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]

def moperatordiveq1(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"{cudadec} {mname}& operator/=(const {dtype}& rhs);\n"
    lnss += f"{cudadec} {mname}& {mname}::operator/=(const {dtype}& rhs)\n"
    lnss += "{\n"
    for I in range(0,dim):
        for J in range(0,dim):
            lnss += f"\tm{J}{I} /= rhs;\n"
        #lnss += "\n"
    lnss += "\treturn *this;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]

def moperatortimeseq2(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"{cudadec} {mname}& operator*=(const {mname}& rhs);\n"
    lnss += f"{cudadec} {mname}& {mname}::operator*=(const {mname}& rhs)\n"
    lnss += "{\n"
    lnss += f"\t{mname} tmp = *this;\n"
    for I in range(0,dim):
        for K in range(0,dim):
            lnss += f"\tm{I}{K} = "
            for J in range(0,dim):
                lnss += f"tmp.m{I}{J}*rhs.m{J}{K}"
                if(J<dim-1):
                    lnss += " + "
                else:
                    lnss += ";\n"
        #lnss += "\n"
    lnss += "\treturn *this;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]

def mtranspose(dim,dtype):
    lnsh = ""
    lnss = ""

    mname = mtypename(dim,dtype)

    lnsh += f"{cudadec} {mname} transpose() const;\n"
    lnss += f"{cudadec} {mname} {mname}::transpose() const\n"
    lnss += "{\n"
    lnss += f"\t{mname} ret;\n"
    for I in range(0,dim):
        for J in range(0,dim):
            lnss += f"\tret.m{J}{I} = m{I}{J};\n"
        #lnss += "\n"
    lnss += "\treturn ret;\n"
    lnss += "}\n\n"

    return [lnsh,lnss]


    

def genmatrixstuff(dim,dtype):

    lnsh = ""
    lnss = ""

    lnsh += "//\n//Matrix Header Stuff\n//\n"
    lnss += "\n//\n//Matrix Header Stuff\n//\n\n"

    [lnsh2,lnss2] = matrix_data(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = mconstructor1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = mdestructor1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = mconstructor2(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = mconstructor3(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatorbk1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatorbk2(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatorbk3(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatorbk4(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatorbk5(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatorbk6(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = mdataoperator1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = mdataoperator2(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatorplus1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatorminus1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatortimes1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatortimes2(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatordiv2(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatortimes3(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatortimes4(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatortimes5(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatorneg(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatorpluseq(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatorminuseq(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatortimeseq1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatordiveq1(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = moperatortimeseq2(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;
    [lnsh2,lnss2] = mtranspose(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;

    return [lnsh,lnss]

def gencode(dim,dtype):
    genheader = "./test_scripts/cuvec_codegen/{}_codegen1.hpp".format(vtypename(dim,dtype))
    gensource = "./test_scripts/cuvec_codegen/{}_codegen1.cu".format(vtypename(dim,dtype))

    lnsh = ""
    lnss = ""

    [lnsh2,lnss2] = genvectorstuff(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;

    [lnsh2,lnss2] = genmatrixstuff(dim,dtype)
    lnsh+=lnsh2; lnss+=lnss2;

    fwrite(genheader,lnsh)
    fwrite(gensource,lnss)

if(__name__=="__main__"):

    dim = 4
    dtype = "int"
    
    for I in range(2,5):
        for d in ["float","double","int"]:
            gencode(I,d)