init

LICENSE

@@ -0,0 +1 @@
+Copyright Aaron M. Schinder, 2023

amsculib3.code-workspace

@@ -0,0 +1,21 @@
+{
+	"folders": [
+		{
+			"path": "."
+		}
+	],
+	//"editor.quickSuggestions": false
+	//https://code.visualstudio.com/docs/editing/intellisense#_customizing-intellisense	
+	// Controls if quick suggestions should show up while typing
+	"settings": {
+		"editor.quickSuggestions": {
+			"other": false,
+			"comments": false,
+			"strings": false
+		},
+		"editor.suggestOnTriggerCharacters": false,
+		"editor.tabCompletion": "off",
+		"editor.parameterHints.enabled": false,
+		"editor.acceptSuggestionOnCommitCharacter": false
+	}
+}

backup.sh

@@ -0,0 +1,6 @@
+#!/usr/bin/bash 
+
+tar -czvf ../amsculib3.tar.gz ./* 
+
+scp ../amsculib3.tar.gz aschinder@amssolarempire.com:~/workspace/projects 
+

build/amsbuildlib4.py

@@ -0,0 +1,813 @@
+#!/usr/bin/python3
+
+import os,sys,math
+import subprocess
+
+"""
+Copyright Aaron M. Schinder, 2011 - MIT/BSD License
+
+This script contains a bunch of helper functions for generating simple, imperative, hopefully
+transparent build scripts using the python language (and nothing else). 
+
+I just want the script to do the compiling and linking operations I want it to do in the order
+I want it to do it in, finding every relevant source file. 
+
+That's it. That's what I want in a build system.
+"""
+
+def flist(pth,**kwargs): 
+    """
+    flist - list all files in a given directory pth
+    optional arguments:
+       recurse - (T/F): Whether to recursively search for files in directory tree
+       exts - (list): A list of file extensions to search for, otherwise all files
+       normpath (T/F): whether to normalize path variables after
+    filelist = flist(pth,**kwargs): 
+    """
+    flst = []
+    if(not('recurse' in kwargs)):
+        recurse_ = False
+    else:
+        recurse_ = kwargs['recurse']
+    if(not('exts' in kwargs)):
+        filterexts_ = False
+    else:
+        filterexts_ = True
+        exts = kwargs['exts']
+    if(not('normpath' in kwargs)):
+        normpath_ = True
+    else:
+        normpath_ = kwargs['normpath']
+    if(not('linuxpath' in kwargs)):
+        linuxpath_ = False
+    else:
+        linuxpath_ = kwargs['linuxpath']
+    if(not('followlinks' in kwargs)):
+        followlinks_ = False
+    else:
+        followlinks_ = kwargs['followlinks']
+
+    dirlist = []
+    rawlist = os.listdir(pth)
+
+    for F in rawlist:
+        F2 = os.path.join(pth,F)
+        if(os.path.isdir(F2)):
+            b = (followlinks_) or ((not followlinks_) and not(os.path.islink(F2)))
+            if(b):
+                if((F2!=".")&(F2!="..")):
+                    dirlist.append(F2)
+        elif(os.path.isfile(F2)):
+            flst.append(F2)
+
+    #Recurse through directories
+    if(recurse_):
+        for D in dirlist:
+            lst = flist(D,**kwargs)
+            for L in lst:
+                flst.append(L)
+
+    #Postprocess:
+    #Filter out all extensions except the selected ext list
+    if(filterexts_):
+        flst = filterexts(flst,exts)
+
+    #Normalize filename path according to os
+    if(normpath_):
+        flst2 = list(flst)
+        for I in range(0,len(flst2)):
+            flst[I] = os.path.normpath(flst2[I])
+
+    #If linuxpath, convert all \\ to /
+    #if(linuxpath_):
+    #    flst2 = list(flst)
+    #    for I in range(0,len(flst2)):
+    #        flst[I] = linuxpath(flst2[I])
+    
+    return flst
+
+
+def filterexts(flst,exts):
+    """
+    Filters by extensions in a list of files
+    flst = def filterexts(flst,exts):
+    """
+    flst2 = []
+    if(isinstance(exts,str)):
+        exts = list([exts])
+    for F in flst:
+        b = False
+        for ext in exts:
+            if(ext[0]!='.'):
+                ext = '.'+ext
+            F2 = os.path.splitext(F)
+            if(len(F2)>=2):
+                ex = F2[1]
+                if(len(ex)>0):            
+                    if(ex[0]!='.'):
+                        ex = '.'+ex
+                    if(ex==ext):
+                        b = True
+        if(b):
+            flst2.append(F)
+              
+    return flst2
+
+#Find a file fname, starting in pth and recursing
+#Used for finding library files to link
+def findfile(fname,pth,**kwargs):
+    fullfname = ""
+    flst = flist(pth,recurse=True)
+    for F in flst:
+        F2 = os.path.split(F)[1]
+        if(F2 == fname):
+            fullfname = F
+    
+    return fullfname
+
+def replaceext(fname,ext):
+    fname2 = ""
+    if(len(ext)>0):
+        if(ext[0]!='.'):
+            ext = '.'+ext
+        fname2 = os.path.splitext(fname)[0]+ext
+    else:
+        fname2 = os.path.splitext(fname)[0]
+    return fname2
+
+def replaceexts(fnamelist,ext):
+    """Takes a list of filenames and returns a list with the extensions replaced by ext """
+    fname2list = []
+    for F in fnamelist:
+        F2 = replaceext(F,ext)
+        fname2list.append(F2)
+    return fname2list
+
+def except_contains(lst1,exc):
+    """
+    Takes a list of file names lst1, and removes filenams that match the
+    list of exceptions exc. Returns a list without the exceptions.
+    """
+    lst2 = []
+    for item in lst1:
+        b = 1
+        for item2 in exc:
+            fsplit = os.path.split(item)
+            fn = fsplit[len(fsplit)-1]
+            if(fn==item2):
+                b = 0
+                break
+        if(b==1):
+            lst2.append(item)
+    return lst2
+
+
+def list_to_sss(lst):
+    """List of strings to space-seperated-string"""
+    lout = ""
+    for I in range(0,len(lst)-1):
+        lout = lout + lst[I] + " "
+    if(len(lst)>0):
+        lout = lout + lst[len(lst)-1]
+    return lout
+
+##########################
+##System Call Procedures##
+##########################
+
+def callproc(cmd, **kwargs):
+    if(not('logfile' in kwargs)):
+        use_lf = False
+    else:
+        logfile = kwargs['logfile']
+        if(logfile!=""):
+            fp = open(kwargs['logfile'],'a+')
+            use_lf = True
+        else:
+            use_lf = False
+        
+    if(not('echo' in kwargs)):
+        echo = True
+    else:
+        echo = kwargs['echo']
+
+    if(echo):
+        print(cmd)
+
+    #encoding/deconding to/from bytes is necessary to use the subprocess command
+    #in python3.7
+    #However, only do this in linux
+    if(sys.platform!='win32'):
+        cmd2 = cmd.encode(encoding='utf-8')
+    else:
+        cmd2 = cmd
+    proc = subprocess.Popen(cmd2,stderr = subprocess.STDOUT, stdout=subprocess.PIPE, shell=True)
+    (out, err) = proc.communicate()
+
+    out = out.decode(encoding='utf-8')
+    
+    if(echo):
+        print(out)
+        #print(err);
+    if(use_lf):
+        fp.writelines(cmd+'\n')
+        fp.writelines(out+'\n')
+
+    if(use_lf):
+        fp.close()
+
+
+#############################################
+## Compiler, Archive, and Linker Functions ##
+#############################################
+
+#MSVC compiler wrapper
+def msvc_compile(compilername, srcfile, **kwargs):
+
+    if(not('include' in kwargs)):
+        include = ''
+    else:
+        include = kwargs['include']
+        if(isinstance(include,list)):
+            include  = list_to_sss(include)
+    
+    if(not('flags' in kwargs)):
+        flags = ''
+    else:
+        flags = kwargs['flags']
+        if(isinstance(flags,list)):
+            flags = list_to_sss(flags)
+
+    if(not('objext' in kwargs)):
+        objext = '.obj'
+    else:
+        objext = kwargs['objext']
+    
+    if(not('srcfileflag' in kwargs)):
+        srcfileflag = '/c'
+    else:
+        srcfileflag = kwargs['srcfileflag']
+
+    if(not('outfileflag' in kwargs)):
+        outfileflag = '/Fo:'
+    else:
+        outfileflag = kwargs['outfileflag']
+    if(not('logfile' in kwargs)):
+        logfile = ""
+    else:
+        logfile = kwargs['logfile']
+
+    outfile = replaceext(srcfile,objext)
+    ln = compilername+" "+flags+" "+" "+srcfileflag+" "+srcfile+" "+outfileflag+'"'+outfile+'"'
+    ln = ln + " " + include
+
+    callproc(ln,echo=True,logfile=logfile)
+
+    return
+
+#MSVC compiler wrapper
+def msvc_compile_list(compiler,srclist,**kwargs):
+    for S in srclist:
+        msvc_compile(compiler,S,**kwargs)
+    return
+
+#gnu-style compiler compile: Should work with gcc, g++, gfortran
+def gs_compile(compiler,srcfile,**kwargs):
+    if(not('include' in kwargs)):
+        include = ''
+    else:
+        include = kwargs['include']
+        if(isinstance(include,list)):
+            include  = list_to_sss(include)
+    
+    if(not('flags' in kwargs)):
+        flags = ''
+    else:
+        flags = kwargs['flags']
+        if(isinstance(flags,list)):
+            flags = list_to_sss(flags)
+
+    if(not('objext' in kwargs)):
+        objext = '.o'
+    else:
+        objext = kwargs['objext']
+
+    if(not('srcfileflag' in kwargs)):
+        srcfileflag = '-c'
+    else:
+        srcfileflag = kwargs['srcfileflag']
+
+    if(not('outfileflag' in kwargs)):
+        outfileflag = '-o'
+    else:
+        outfileflag = kwargs['outfileflag']
+
+    if(not('logfile' in kwargs)):
+        logfile = ""
+    else:
+        logfile = kwargs['logfile']
+
+    if(not('smartcompile' in kwargs)):
+        _smartcompile = True
+    else:
+        _smartcompile = kwargs['smartcompile']
+    
+    #Do I want to make this thing this general?
+
+    # if(not(_smartcompile) or smartcompile(srcfile,objext)):
+    #     outfile = replaceext(srcfile,objext)
+    #     ln = compiler+" "+flags+" " + outfileflag+" "+outfile+" "+srcfileflag+" "+srcfile
+    #     ln = ln + " " + include
+
+    #     callproc(ln,echo=True,logfile=logfile)
+
+    outfile = replaceext(srcfile,objext)
+    ln = compiler+" "+flags+" " + outfileflag+" "+outfile+" "+srcfileflag+" "+srcfile
+    ln = ln + " " + include
+
+    callproc(ln,echo=True,logfile=logfile)
+
+    return
+
+def gs_compile_list(compiler,srclist,**kwargs):
+    for S in srclist:
+        gs_compile(compiler,S,**kwargs)
+    return
+
+def gs_compile_all(compiler,srcdir,srcexts,**kwargs):
+    if(not('recurse' in kwargs)):
+        recurse = True
+    else:
+        recurse = kwargs['recurse']
+
+    srcfils = flist(srcdir,exts=srcexts,recurse=recurse)
+    
+    for S in srcfils:
+        gs_compile(compiler,S,**kwargs)
+
+    return
+
+def gs_link_all(linker,srcpath,target,**kwargs):
+
+    if(not('objext' in kwargs)):
+        objext = '.o'
+    else:
+        objext = kwargs['objext']
+
+    if(not('recurse' in kwargs)):
+        recurse = True
+    else:
+        recurse = kwargs['recurse']
+
+
+    objfils = flist(srcpath,exts=objext,recurse=recurse)
+    oflst = list_to_sss(objfils)
+    
+    gs_link_list(linker,oflst,target,**kwargs)
+    
+    return
+
+def gs_link_list(linker,objlist,target,**kwargs):
+
+    if(not('objext' in kwargs)):
+        objext = '.o'
+    else:
+        objext = kwargs['objext']
+
+    if(not('libdir' in kwargs)):
+        libdir = ''
+    else:
+        libdir = kwargs['libdir']
+
+    if(not('staticlibs' in kwargs)):
+        staticlibs = ''
+    else:
+        staticlibs = kwargs['staticlibs']
+
+    if(not('libflags' in kwargs)):
+        libflags = ''
+    else:
+        libflags = kwargs['libflags']
+
+    if(not('linkerflags' in kwargs)):
+        linkerflags = ''
+    else:
+        linkerflags = kwargs['linkerflags']
+
+    if(not('recurse' in kwargs)):
+        recurse = True
+    else:
+        recurse = kwargs['recurse']
+
+    if(not('logfile' in kwargs)):
+        logfile = ''
+    else:
+        logfile = kwargs['logfile']
+
+    ln = linker+" -o "+target+" "+libdir
+    ln = ln+" "+objlist+" "+staticlibs+" "+libflags+" "+linkerflags
+    
+    callproc(ln,logfile=logfile)
+    return
+
+def msvc_link_list(objlist,target,**kwargs):
+
+    linker = 'link'
+
+    if(not('objext' in kwargs)):
+        objext = '.obj'
+    else:
+        objext = kwargs['objext']
+
+    if(not('libdir' in kwargs)):
+        libdir = ''
+    else:
+        libdir = kwargs['libdir']
+
+    if(not('staticlibs' in kwargs)):
+        staticlibs = ''
+    else:
+        staticlibs = kwargs['staticlibs']
+
+    if(not('libflags' in kwargs)):
+        libflags = ''
+    else:
+        libflags = kwargs['libflags']
+
+    if(not('linkerflags' in kwargs)):
+        linkerflags = ''
+    else:
+        linkerflags = kwargs['linkerflags']
+
+    if(not('recurse' in kwargs)):
+        recurse = True
+    else:
+        recurse = kwargs['recurse']
+
+    if(not('logfile' in kwargs)):
+        logfile = ''
+    else:
+        logfile = kwargs['logfile']
+
+    ln = linker+" "+libdir
+    ln = ln+" "+objlist+" "+staticlibs+" "+linkerflags
+    ln = ln+" /out:"+target+" "+libflags
+    
+    callproc(ln,logfile=logfile)
+
+    return
+
+def ar_all(srcpath,arname,**kwargs):
+    if(not('recurse' in kwargs)):
+        recurse = True
+    else:
+        recurse = kwargs['recurse']
+    if(not('objext' in kwargs)):
+        objext = '.o'
+    else:
+        objext = kwargs['objext']
+    
+    objlist = flist(srcpath,exts=objext,recurse=recurse)
+    ar_list(objlist,arname,**kwargs)
+    
+    return
+
+def msvc_lib_list(objlist,arname,**kwargs):
+    objlist2 = list_to_sss(objlist)
+    
+    ln = "lib "+objlist2+" /out:"+arname
+    callproc(ln)
+
+    return
+
+def ar_list(objlist,arname,**kwargs):
+    objlist2 = list_to_sss(objlist)
+
+    ln = "ar cr "+ arname+" "+objlist2
+    callproc(ln)
+    
+    return
+
+def ar_add_list(objlist,arname,**kwargs):
+    objlist2 = list_to_sss(objlist)
+
+    ln = "ar t "+arname+" "+objlist2
+    callproc(ln)
+    return
+
+#####################################
+## Incremental Compilation Library ##
+#####################################
+
+#silently read lines from a text file if exists
+def readtextlines(fname):
+    txtlns = []
+
+    if(not os.path.isfile(fname)):
+        return txtlns
+
+    try:
+        fp = open(fname,"r")
+    except:
+        return txtlns
+
+    ln = " "
+    while(ln!=""):
+        ln = fp.readline()
+        txtlns.append(ln)
+
+    fp.close()
+
+    return txtlns
+
+def getincludefnfrage(includeline):
+
+    fnfrag = ""
+    I1 = -1
+    I2 = -1
+
+    for I in range(0,len(includeline)):
+        if(I1<0 and (includeline[I]=='<' or includeline[I]=='"')):
+            I1 = I
+        if(I1>=0 and (includeline[I]=='>' or includeline[I]=='"')):
+            I2 = I
+            break
+    if(I1>=0 and I2>=0):
+        fnfrag = includeline[I1+1:I2]
+
+    return fnfrag
+
+#Returns the name of the source file fname (if it exists)
+#and all included filenames
+def getsrcandincludes(fname, incdirs):
+
+    flist = []
+    if(os.path.isfile(fname)):
+        flist.append(fname)
+
+        Ilist = 0
+        while(Ilist<len(flist)):
+            #recurse through files
+            f1 = flist[Ilist]
+            lns = readtextlines(f1)
+            for J in range(0,len(lns)):
+                if(lns[J].find("#include")>=0):
+                    fnfrag = getincludefnfrage(lns[J])
+                    for K in range(0,len(incdirs)):
+                        tfn = os.path.join(incdirs[K],fnfrag)
+                        if(os.path.isfile(tfn)):
+                            flist.append(tfn)
+                            break
+
+            Ilist = Ilist + 1
+
+    return flist
+
+#Returns the name of the object file associated with the source file
+#within the object store folder (if it exists)
+def getobjfile(fname,objstore,objext = ".o"):
+
+    fret = ""
+    f1 = os.path.split(fname)[1]
+    f2 = f1
+    while(os.path.splitext(f2)[1]!=""):
+        f2 = os.path.splitext(f2)[0]
+    objext = objext.strip('.')
+    f3 = os.path.join(objstore,"{}.{}".format(f2,objext))
+    if(os.path.exists(f3)):
+        fret = f3
+
+    return fret
+
+def getsrctimes(fname, incdirs):
+
+    ftimes = []
+    flst = getsrcandincludes(fname, incdirs)
+    for I in range(0,len(flst)):
+        f = flst[I]
+        mt = os.path.getmtime(f)
+        ftimes.append(mt)
+    
+    return ftimes
+
+def getobjtime(fname,objstore,objext=".o"):
+    ret = -1
+    fret = getobjfile(fname,objstore,objext)
+    if(fret!=""):
+        ret = os.path.getmtime(fret)
+
+    return ret
+
+#Decide whether or not to compile source file
+def decidecompile(fname,**kwargs):
+    ret = True
+
+    if(not os.path.isfile(fname)):
+        ret = False
+        return ret
+    
+    ##unpack kwargs
+    if("searchincdirs" in kwargs):
+        incdirs = kwargs["searchincdirs"]
+    else:
+        incdirs = ["./include"]
+    
+    if("objext" in kwargs):
+        objext = kwargs["objext"]
+    else:
+        objext = ".o"
+    if("objstore" in kwargs):
+        objstore = kwargs["objstore"]
+    else:
+        objstore = "./objstore"
+    
+    
+    srclist = getsrcandincludes(fname,incdirs)
+    srctlist = getsrctimes(fname,incdirs)
+    obj = getobjfile(fname,objstore,objext)
+    objt = getobjtime(fname,objstore,objext)
+    
+    if(obj!=""):
+        ret = False
+        for I in range(0,len(srctlist)):
+            if(srctlist[I]>objt):
+                ret = True
+                break
+
+    return ret
+
+def gs_incremental_compile(compiler,srcfile,**kwargs):
+
+    if(not('include' in kwargs)):
+        include = ''
+    else:
+        include = kwargs['include']
+        if(isinstance(include,list)):
+            include  = list_to_sss(include)
+    if(not('flags' in kwargs)):
+        flags = ''
+    else:
+        flags = kwargs['flags']
+        if(isinstance(flags,list)):
+            flags = list_to_sss(flags)
+    if(not('objext' in kwargs)):
+        objext = '.o'
+    else:
+        objext = kwargs['objext']
+    if(not('srcfileflag' in kwargs)):
+        srcfileflag = '-c'
+    else:
+        srcfileflag = kwargs['srcfileflag']
+    if(not('outfileflag' in kwargs)):
+        outfileflag = '-o'
+    else:
+        outfileflag = kwargs['outfileflag']
+
+    if(not('logfile' in kwargs)):
+        logfile = ""
+    else:
+        logfile = kwargs['logfile']
+    if(not('smartcompile' in kwargs)):
+        _smartcompile = True
+    else:
+        _smartcompile = kwargs['smartcompile']
+
+    #incrementalcompile
+    if("searchincdirs" in kwargs):
+        incdirs = kwargs["searchincdirs"]
+    else:
+        incdirs = ["./include"]
+    
+    if("objstore" in kwargs):
+        objstore = kwargs["objstore"]
+    else:
+        objstore = "./objstore"
+    
+    #Do I want to make this thing this general?
+
+    docompile = decidecompile(srcfile,**kwargs)
+
+    if(docompile):
+        f1 = os.path.split(srcfile)[1]
+        f2 = f1
+        while(os.path.splitext(f2)[1]!=""):
+            f2 = os.path.splitext(f2)[0]
+        outfile = os.path.join(objstore,"{}{}".format(f2,objext))
+
+        ln = compiler+" "+flags+" " + outfileflag+" "+outfile+" "+srcfileflag+" "+srcfile
+        ln = ln + " " + include
+
+        callproc(ln,echo=True,logfile=logfile)
+
+    return
+
+def msvc_incremental_compile(compiler,srcfile,**kwargs):
+
+    if(not('include' in kwargs)):
+        include = ''
+    else:
+        include = kwargs['include']
+        if(isinstance(include,list)):
+            include  = list_to_sss(include)
+    if(not('flags' in kwargs)):
+        flags = ''
+    else:
+        flags = kwargs['flags']
+        if(isinstance(flags,list)):
+            flags = list_to_sss(flags)
+    if(not('objext' in kwargs)):
+        objext = '.obj'
+    else:
+        objext = kwargs['objext']
+    if(not('srcfileflag' in kwargs)):
+        srcfileflag = '/c'
+    else:
+        srcfileflag = kwargs['srcfileflag']
+    if(not('outfileflag' in kwargs)):
+        outfileflag = '/Fo'
+    else:
+        outfileflag = kwargs['outfileflag']
+
+    if(not('logfile' in kwargs)):
+        logfile = ""
+    else:
+        logfile = kwargs['logfile']
+
+    #incrementalcompile
+    if("searchincdirs" in kwargs):
+        incdirs = kwargs["searchincdirs"]
+    else:
+        incdirs = ["./include"]
+    
+    if("objstore" in kwargs):
+        objstore = kwargs["objstore"]
+    else:
+        objstore = "./objstore"
+    
+    #Do I want to make this thing this general?
+
+    docompile = decidecompile(srcfile,**kwargs)
+
+    if(docompile):
+        f1 = os.path.split(srcfile)[1]
+        f2 = f1
+        while(os.path.splitext(f2)[1]!=""):
+            f2 = os.path.splitext(f2)[0]
+        outfile = os.path.join(objstore,"{}{}".format(f2,objext))
+        outfile = os.path.normpath(outfile)
+
+        ln = compiler+" "+flags+" "+srcfileflag+" "+srcfile+" "+ outfileflag+'"'+outfile+'"'
+        ln = ln + " " + include
+
+        callproc(ln,echo=True,logfile=logfile)
+
+    return
+
+
+def gs_incremental_compile_list(compiler,srclist,**kwargs):
+
+    for s in srclist:
+        gs_incremental_compile(compiler,s,**kwargs)
+
+    return
+
+def msvc_incremental_compile_list(compiler,srclist,**kwargs):
+
+    for s in srclist:
+        msvc_incremental_compile(compiler,s,**kwargs)
+
+    return
+
+#######################
+## Main Script Tests ##
+#######################
+
+def testtimes(args):
+    if(len(args)>=2):
+        flist = getsrcandincludes(args[1],["./include"])
+        ftlist = getsrctimes(args[1],["./include"])
+        for I in range(0,len(flist)):
+            print("{}\t\t{}".format(flist[I],ftlist[I]))
+        
+        print("associated obj file:")
+        fobj = getobjfile(args[1],"./objstore")
+        ftobj = getobjtime(args[1],"./objstore")
+        if(fobj!=""):
+            print("{}\t\t{}".format(fobj,ftobj))
+        else:
+            print("none found")
+
+        cflag = decidecompile(args[1])
+        print("compile? : {}".format(cflag))
+
+    
+    return
+
+# if(__name__ == "__main__"):
+
+#     args = sys.argv
+#     testtimes(args)
+    
+        
+
+        

build/make.linux64.lib.py

@@ -0,0 +1,58 @@
+#!/usr/bin/python3
+
+import os,sys,math
+import subprocess
+import shutil
+from shutil import copytree
+
+from amsbuildlib4 import *
+
+libname = "amsculib3.linux64" #static library name to generate
+binname = "test" #create this executable when compiling main.c or main.cpp
+commondir = "../../linux64" #common directory to pul libraries and includes from
+depdir = "./dependencies/linux64" #local pre-compiled dependency libraries and their includes
+installdir = "../../linux64" #directory to install to when finished
+builddir = "./build_linux64"
+
+doinstall = True   #copies the build_output to the install dir when finished
+cc = "nvcc" #compiler
+cflags = "-dc --compiler-options '-fPIC -O3'"
+libraries = "-l{}".format(libname)
+libdirs = "-L{} -L{}/lib -L{}/lib".format(builddir,commondir,depdir)
+linkerflags = " -Xlinker=-rpath,."
+srcexts = [".c",".cpp",".cu"]
+binsrc = ["main.c","main.cpp", "main.cu"] #ignore these files when compiling the static library
+
+#keyword list to control the compilers/linkers
+kwargs = dict()
+include = "-I./include -I{}/include -I{}/include".format(commondir, depdir)
+kwargs["include"] = include
+kwargs["flags"] = cflags
+kwargs["libdir"] = libdirs
+kwargs["libflags"] = libraries
+kwargs["linkerflags"] = linkerflags
+kwargs["recurse"] = True
+kwargs["objstore"] = "{}/objstore".format(builddir)
+kwargs["searchincdirs"] = "./include"
+
+#Find all source files, except the main project files
+srcfiles = flist('./src',exts = srcexts, recurse=True)
+srcfiles = except_contains(srcfiles,binsrc)
+
+#compile all the source files in the list
+#gs_compile_list(cc,files,**kwargs)
+gs_incremental_compile_list(cc,srcfiles,**kwargs)
+
+#archive all the source files into a static library
+objlist = flist(kwargs['objstore'],exts='.o',recurse=True)
+ar_list(objlist,'{}/lib{}.a'.format(builddir,libname))
+
+if(doinstall):
+    #Push any libraries to the common lib folder
+    shutil.copy(
+        '{}/lib{}.a'.format(builddir,libname),
+        "{}/lib".format(installdir)
+    )
+
+    #Copy include files to the common include folder
+    copytree('./include/',installdir+'/include/',dirs_exist_ok=True)

build/make.linux64.test.py

@@ -0,0 +1,49 @@
+#!/usr/bin/python3
+
+import os,sys,math
+import subprocess
+import shutil
+from shutil import copytree
+
+from amsbuildlib4 import *
+
+libname = "amsculib3.linux64" #static library name to generate
+binname = "test" #create this executable when compiling main.c or main.cpp
+commondir = "../../linux64" #common directory to pul libraries and includes from
+depdir = "./dependencies/linux64" #local pre-compiled dependency libraries and their includes
+installdir = "../../linux64" #directory to install to when finished
+builddir = "./build_linux64"
+
+doinstall = True   #copies the build_output to the install dir when finished
+cc = "nvcc" #compiler
+cflags = "-dc --compiler-options '-fPIC -O3'"
+libraries = "-l{}".format(libname)
+libdirs = "-L{} -L{}/lib -L{}/lib".format(builddir,commondir,depdir)
+linkerflags = " -Xlinker=-rpath,."
+srcexts = [".c",".cpp",".cu"]
+binsrc = ["main.c","main.cpp", "main.cu"] #ignore these files when compiling the static library
+
+#keyword list to control the compilers/linkers
+kwargs = dict()
+include = "-I./include -I{}/include -I{}/include".format(commondir, depdir)
+kwargs["include"] = include
+kwargs["flags"] = cflags
+kwargs["libdir"] = libdirs
+kwargs["libflags"] = libraries
+kwargs["linkerflags"] = linkerflags
+kwargs["recurse"] = True
+kwargs["objstore"] = "{}/objstore".format(builddir)
+kwargs["searchincdirs"] = "./include"
+
+#Pull required binary dynamic libraries to the bin folder
+#shutil.copy('{}/lib/libcamsimg3.linux64.so'.format(commondir),builddir);
+#shutil.copy('{}/lib/libamsimg.dll'.format(commondir),builddir);
+#shutil.copy('{}/lib/glew32.dll','./bin_winx64');
+
+#Designate source files for main test program
+fsrc = ['./src/main.cu']
+fobj = replaceexts(fsrc,'.o')
+
+#Compile test programs
+gs_compile_list(cc,fsrc,**kwargs)
+gs_link_list(cc,list_to_sss(fobj),'{}/{}'.format(builddir,binname),**kwargs)

build/make.msvc64.lib.py

@@ -0,0 +1,61 @@
+#!/usr/bin/python3
+
+import os,sys,math
+import subprocess
+import shutil
+from shutil import copytree
+
+from amsbuildlib4 import *
+
+libname = "amsculib3.msvc64" #static library name to generate
+binname = "test" #create this executable when compiling main.c or main.cpp
+commondir = "../../winx64" #common directory to pul libraries and includes from
+depdir = "./dependencies/winx64" #local pre-compiled dependency libraries and their includes
+installdir = "../../winx64" #directory to install to when finished
+builddir = "./build_msvc64"
+
+doinstall = True   #copies the build_output to the install dir when finished
+cc = "nvcc" #compiler
+cflags = "-dc --compiler-options '-fPIC -O3'"
+libraries = "-l{}".format(libname)
+libdirs = "-L{} -L{}/lib -L{}/lib".format(builddir,commondir,depdir)
+linkerflags = " -Xlinker=-rpath,."
+srcexts = [".c",".cpp",".cu"]
+binsrc = ["main.c","main.cpp","main.cu"] #ignore these files when compiling the static library
+
+
+#keyword list to control the compilers/linkers
+kwargs = dict()
+include = "-I./include -I{}/include -I{}/include".format(commondir, depdir)
+kwargs["include"] = include
+kwargs["flags"] = cflags
+kwargs["libdir"] = libdirs
+kwargs["libflags"] = libraries
+kwargs["linkerflags"] = linkerflags
+kwargs["recurse"] = True
+kwargs["objstore"] = "{}/objstore".format(builddir)
+kwargs["searchincdirs"] = "./include"
+kwargs["objext"] = ".obj"
+
+
+#Find all source files, except the main project files
+srcfiles = flist('./src',exts = srcexts, recurse=True)
+srcfiles = except_contains(srcfiles,binsrc)
+
+#compile all the source files in the list
+#gs_compile_list(cc,files,**kwargs)
+msvc_incremental_compile_list(cc,srcfiles,**kwargs)
+
+#archive all the source files into a static library
+objlist = flist(kwargs['objstore'],exts='.obj',recurse=True)
+msvc_lib_list(objlist,'{}/lib{}.lib'.format(builddir,libname))
+
+if(doinstall):
+    #Push any libraries to the common lib folder
+    shutil.copy(
+        '{}/lib{}.lib'.format(builddir,libname),
+        "{}/lib".format(installdir)
+    )
+
+    #Copy include files to the common include folder
+    copytree('./include/',installdir+'/include/',dirs_exist_ok=True)

build/make.msvc64.test.py

@@ -0,0 +1,49 @@
+#!/usr/bin/python3
+
+import os,sys,math
+import subprocess
+import shutil
+from shutil import copytree
+
+from amsbuildlib4 import *
+
+libname = "amsculib3.msvc64" #static library name to generate
+binname = "test.exe" #create this executable when compiling main.c or main.cpp
+commondir = "../../winx64" #common directory to pul libraries and includes from
+depdir = "./dependencies/winx64" #local pre-compiled dependency libraries and their includes
+installdir = "../../winx64" #directory to install to when finished
+builddir = "./build_msvc64"
+
+doinstall = False   #copies the build_output to the install dir when finished
+cc = "nvcc" #compiler
+cflags = "-dc --compiler-options '-fPIC -O3'"
+libraries = "-l{}".format(libname)
+libdirs = "-L{} -L{}/lib -L{}/lib".format(builddir,commondir,depdir)
+linkerflags = " -Xlinker=-rpath,."
+srcexts = [".c",".cpp",".cu"]
+binsrc = ["main.c","main.cpp","main.cu"] #ignore these files when compiling the static library
+
+#keyword list to control the compilers/linkers
+kwargs = dict()
+include = "-I./include -I{}/include -I{}/include".format(commondir, depdir)
+kwargs["include"] = include
+kwargs["flags"] = cflags
+kwargs["libdir"] = libdirs
+kwargs["libflags"] = libraries
+kwargs["linkerflags"] = linkerflags
+kwargs["recurse"] = True
+kwargs["objstore"] = "{}/objstore".format(builddir)
+kwargs["searchincdirs"] = "./include"
+
+#Pull required binary dynamic libraries to the bin folder
+#shutil.copy('{}/lib/libcamsimg3.linux64.so'.format(commondir),builddir);
+#shutil.copy('{}/lib/libamsimg.dll'.format(commondir),builddir);
+#shutil.copy('{}/lib/glew32.dll','./bin_winx64');
+
+#Designate source files for main test program
+fsrc = ['./src/main.cpp']
+fobj = replaceexts(fsrc,'.obj')
+
+#Compile test programs
+msvc_compile_list(cc,fsrc,**kwargs)
+msvc_link_list(list_to_sss(fobj),'{}/{}'.format(builddir,binname),**kwargs)

include/amsculib3/amscu_cudafunctions.hpp

@@ -0,0 +1,40 @@
+#ifndef __AMSCU_CUDAFUNCTIONS_HPP__
+#define __AMSCU_CUDAFUNCTIONS_HPP__
+
+
+namespace amscuda
+{
+    // device memory operations
+    // I'm trying to avoid some of the boilerplate mental overhead involved
+    // in calling cuda functions and handling errors
+
+    //frees devbuffer if it is not already NULL, and sets devbuffer to NULL
+    //wrapper to cudaFree
+    template<typename T> int cuda_free(T **devptr);
+
+    //copies hostbuffer to devbuffer
+    //initializes devbuffer from NULL if devbuffer is NULL
+    //if overwrite is true, deletes and reallocates devbuffer on device (for resizing)
+    template<typename T> int buffer_copytodevice(T *hostbuffer, T **devbuffer, long N, bool overwrite);
+
+    //copies info from devbuffer to hostbuffer
+    //initialzies hostbuffer from NULL if NULL
+    //if overwrite is true, deletes and reallocates hostbuffer on host with new[] (for resizing)
+    template<typename T> int buffer_copyfromdevice(T *devbuffer, T **hostbuffer, long N, bool overwrite);
+
+    //wrapper for cudaMemcpy - copies an item or struct (count 1) to the device
+    //initializes devptr from NULL if not already initialized
+    template<typename T> int cuda_copytodevice(T *hostptr, T **devptr);
+
+    //wrapper for cudaMemcpy - copies an item or struct (count 1) from device
+    //initializes hostptr from NULL with new if not already initialized
+    template<typename T> int cuda_copyfromdevice(T *devptr, T **hostptr);
+    
+    int cuda_errortrap(const char *msgheader);
+
+};
+
+#include <amsculib2/amscu_cudafunctions_impl.hpp>
+
+#endif
+

include/amsculib3/amscu_cudafunctions_impl.hpp

@@ -0,0 +1,228 @@
+#ifndef __AMSCU_CUDAFUNCTIONS_IMPL_HPP__
+#define __AMSCU_CUDAFUNCTIONS_IMPL_HPP__
+
+namespace amscuda
+{
+
+//frees devbuffer if it is not already NULL, and sets devbuffer to NULL
+//wrapper to cudaFree
+template<typename T> int cuda_free(T **devptr)
+{
+    int ret = 0;
+    cudaError_t err = cudaSuccess;
+
+    if(*devptr==NULL)
+    {
+        return ret; //devbuffer is already NULL/freed
+    }
+
+    err = cudaFree(*devptr);
+    if(err!=cudaSuccess)
+    {
+        ret = -1; //failed to free device pointer
+        *devptr = NULL; // - ? should only happen if I'm trying to double-free something
+    }
+    else
+    {
+        ret = 1;
+        *devptr = NULL;
+    }
+
+    return ret;
+}
+
+//copies hostbuffer to devbuffer
+//initializes devbuffer from NULL if devbuffer is NULL
+//if overwrite is true, deletes and reallocates devbuffer on device (for resizing)
+template<typename T> int buffer_copytodevice(T *hostbuffer, T **devbuffer, long N, bool overwrite)
+{
+    int ret = 0;
+    cudaError_t err = cudaSuccess;
+
+    if(N<=0)
+    {
+        ret = 0;
+        return ret;
+    }
+
+    if(hostbuffer==NULL)
+    {
+        ret = -2; //host buffer is NULL
+        return ret;
+    }
+
+    if(overwrite==1)
+    {
+        if(*devbuffer !=NULL)
+        {
+            cuda_free(devbuffer);
+        }
+    }
+
+    if(*devbuffer==NULL)
+    {
+        err = cudaMalloc(devbuffer,sizeof(T)*N);
+        if(err!=cudaSuccess)
+        {
+            ret = -3; //failed to allocate
+            *devbuffer = NULL;
+            return ret;
+        }
+    }
+
+    err = cudaMemcpy(*devbuffer,hostbuffer,sizeof(T)*N,cudaMemcpyHostToDevice);
+    if(err!=cudaSuccess)
+    {
+        ret = -4; //failed to copy
+    }
+    else
+    {
+        ret = 1;
+    }
+
+
+    return ret;
+}
+
+//copies info from devbuffer to hostbuffer
+//initialzies hostbuffer from NULL if NULL
+//if overwrite is true, deletes and reallocates hostbuffer on host (for resizing)
+template<typename T> int buffer_copyfromdevice(T *devbuffer, T **hostbuffer, long N, bool overwrite)
+{
+    int ret = 0;
+    cudaError_t err = cudaSuccess;
+
+    if(N<=0)
+    {
+        ret = 0;
+        return ret;
+    }
+
+    if(devbuffer==NULL)
+    {
+        ret = -5; //null dev buffer
+        return ret;
+    }
+
+    if(overwrite==1 && *hostbuffer!=NULL)
+    {
+        delete[] (*hostbuffer); hostbuffer = NULL;
+    }
+
+    if(*hostbuffer==NULL)
+    {
+        *hostbuffer = new(std::nothrow) T[N];
+        if(*hostbuffer==NULL)
+        {
+            ret = -6; //failed to allocate host buffer
+            return ret;
+        }
+    }
+
+    err = cudaMemcpy(*hostbuffer, devbuffer, sizeof(T)*N, cudaMemcpyDeviceToHost);
+    if(err!=cudaSuccess)
+    {
+        ret = -7; //failed to copy
+    }
+    else
+    {
+        ret = 1;
+    }
+    
+    return ret;
+}
+
+//wrapper for cudaMemcpy - copies an item or struct (count 1) to the device
+//initializes devptr from NULL if not already initialized
+template<typename T> int cuda_copytodevice(T *hostptr, T **devptr)
+{
+    int ret = 0;
+    cudaError_t err = cudaSuccess;
+    bool overwrite = 1;
+
+    if(hostptr==NULL)
+    {
+        ret = -2; //host buffer is NULL
+        return ret;
+    }
+
+    if(overwrite==1)
+    {
+        if(*devptr !=NULL)
+        {
+            cuda_free(devptr);
+        }
+    }
+
+    if(*devptr==NULL)
+    {
+        err = cudaMalloc(devptr,sizeof(T));
+        if(err!=cudaSuccess)
+        {
+            ret = -3; //failed to allocate
+            *devptr = NULL;
+            return ret;
+        }
+    }
+
+    err = cudaMemcpy(*devptr,hostptr,sizeof(T),cudaMemcpyHostToDevice);
+    if(err!=cudaSuccess)
+    {
+        ret = -4; //failed to copy
+    }
+    else
+    {
+        ret = 1;
+    }
+
+
+    return ret;
+}
+
+//wrapper for cudaMemcpy - copies an item or struct (count 1) from device
+//initializes hostptr from NULL with new if not already initialized
+template<typename T> int cuda_copyfromdevice(T *devptr, T **hostptr)
+{
+    int ret = 0;
+    cudaError_t err = cudaSuccess;
+    bool overwrite = 1;
+
+    if(devptr==NULL)
+    {
+        ret = -5; //null dev buffer
+        return ret;
+    }
+
+    if(overwrite==1 && *hostptr!=NULL)
+    {
+        delete (*hostptr); hostptr = NULL;
+    }
+
+    if(*hostptr==NULL)
+    {
+        *hostptr = new(std::nothrow) T;
+        if(*hostptr==NULL)
+        {
+            ret = -6; //failed to allocate host buffer
+            return ret;
+        }
+    }
+
+    err = cudaMemcpy(*hostptr, devptr, sizeof(T), cudaMemcpyDeviceToHost);
+    if(err!=cudaSuccess)
+    {
+        ret = -7; //failed to copy
+    }
+    else
+    {
+        ret = 1;
+    }
+    
+    return ret;
+}
+
+
+};
+
+#endif
+

include/amsculib3/amscu_random.hpp

@@ -0,0 +1,55 @@
+#ifndef __AMSCU_RANDOM_HPP__
+#define __AMSCU_RANDOM_HPP__
+
+namespace amscuda
+{
+
+// Random Number Gerneators
+
+
+// faster floating point hash function used in fractal generators
+__device__ __host__ float fhash1d_su(float x);
+
+__device__ __host__ float fhash3d_su(float x, float y, float z);
+
+__device__ __host__ float fhash4d_su(float x, float y, float z, float w);
+
+
+//////////////////////////////////////////////////
+// Deterministic Pseudorandom int32_t Generator //
+//////////////////////////////////////////////////
+
+//Next seed in simple 32 bit integer deterministic psuedo-rand generator
+__host__ __device__ void dpr32_nextseed(int32_t *rseed_inout);
+
+//Simple 32 bit integer deterministic pseudo-random generator
+// *not* for cryptography
+// Frequency of generated floats should be uniform [0,1)
+__host__ __device__ float dpr32_randf(int32_t *rseed_inout);
+
+//box muller standard normal pseudorandom variable
+__host__ __device__ float dpr32_randnf(int32_t *rseed_inout);
+
+//////////////////////////////////////////////////
+// Deterministic Pseudorandom int64_t Generator //
+//////////////////////////////////////////////////
+
+//operates without side-effects on explicit seed for threaded use
+//deterministic pseudorandom number generator - takes seed and returns next seed
+__host__ __device__ void dpr64_nextseed(int64_t *seedinout);
+
+//deterministic pseudorandom number generator - takes seed and returns next seed
+//returns uniformly distributed double
+__host__ __device__ double dpr64_randd(int64_t *seedinout);
+
+__host__ __device__ float dpr64_randf(int64_t *seedinout);
+
+
+void test_dprg64();
+void test_dprg32();
+
+
+}; //end namespace amscuda
+
+#endif
+

include/amsculib3/amscuarray.hpp

@@ -0,0 +1,47 @@
+#ifndef __CUARRAY_HPP__
+#define __CUARRAY_HPP__
+
+namespace amscuda
+{
+
+template<typename T> class cuarray
+{
+public:
+    int length;
+    T* data;
+
+    __device__ __host__ cuarray();
+    __device__ __host__ ~cuarray();
+
+    //Only call this on the device for thread/block local 
+    // dynamic arrays
+    __device__ __host__ int resize(const int _length);
+
+    __device__ __host__ int size() const;
+    __device__ __host__ T& at(const int I);
+    __device__ __host__ const T& at(const int I) const;
+
+    __device__ __host__ T& operator[](const int I);
+    __device__ __host__ const T& operator[](const int I) const;
+    
+    
+
+    __host__ int device_send(cuarray<T> **dptr);
+    __host__ int _device_send_overwrite(cuarray<T> **dptr);
+    __host__ int _device_send_copy(cuarray<T> *dptr);
+    
+    __host__ int device_pull(cuarray<T> *dptr);
+    __host__ int device_free(cuarray<T> **dptr);
+
+    __host__ int device_length(cuarray<T> *dptr);
+    __host__ T* device_data_ptr(cuarray<T> *dptr);
+    
+};
+
+void test_cuarray();
+
+};
+
+#include <amsculib2/amscuarray_impl.hpp>
+
+#endif

include/amsculib3/amscuarray_dops.hpp

@@ -0,0 +1,76 @@
+#ifndef __AMSCUARRAY_DOPS_HPP__
+#define __AMSCUARRAY_DOPS_HPP__
+
+//Device Operations on Arrays
+//
+
+//Device Operations on Device Buffers
+// dodb
+
+namespace amscuda
+{
+
+
+    //sum
+    template<typename T> T devcuarray_sum(cuarray<T> *devptr);
+
+    template<typename T> T dbuff_sum(T *devbuffer, int N);
+
+
+    struct dbuff_statstruct
+    {
+        public:
+        float min;
+        float max;
+        float mean;
+        float stdev;
+        float sum;
+    };
+    
+    //stats (min,max,mean,stdev)
+
+    template<typename T> void dbuff_minmax(T *devbuffer, int N, T *min, T *max);
+
+    template<typename T> dbuff_statstruct dbuff_stats(T *devbuffer, int N); //
+
+    //sets all elements to setto
+    template<typename T> void dbuff_setall(T *devbuffer, int N, T setto, int nblocks, int nthreads);
+
+    //random device buffer functions
+    void dbuff_rand_dpr32(float *devbuffer, int N, int32_t *rseedinout,  int nblocks, int nthreads); //
+    void dbuff_rand_dpr32n(float *devbuffer, int N, int32_t *rseedinout,  int nblocks, int nthreads); //
+    
+
+    void dbuff_rand_dpr64(float *devbuffer, int N, int64_t *rseedinout,  int nblocks, int nthreads); //
+
+    //Elementwise device-buffer vector binary operation
+    //takes two input arrays ( , ) --> one output array
+    template<typename T1, typename T2, typename T3> void dbuff_vectorbinop(T1 *dbuf_a, T2 *dbuf_b, T3 *dbuf_out, int N, T3 (*fpnt)(T1,T2), int nblocks, int nthreads);
+
+    //Elementwise device-buffer vector two-parameter operation
+    //takes one input array, and a constant paramter ( ) ---> one output array
+    template<typename T1, typename T2, typename T3> void dbuff_vectorbinop(T1 *dbuf_a, T2 par_b, T3 *dbuf_out, int N, T3 (*fpnt)(T1,T2), int nblocks, int nthreads);
+
+
+    //vector_add
+    template<typename T> void dbuff_add(T *dbuff_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+    template<typename T> void dbuff_add(T *dbuff_a, T par_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+    template<typename T> void dbuff_sub(T *dbuff_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+    template<typename T> void dbuff_sub(T *dbuff_a, T par_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+    template<typename T> void dbuff_mult(T *dbuff_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+    template<typename T> void dbuff_mult(T *dbuff_a, T par_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+    template<typename T> void dbuff_div(T *dbuff_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+    template<typename T> void dbuff_div(T *dbuff_a, T par_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+    template<typename T> void dbuff_div(T par_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+
+
+    // Tests //
+
+    void test_dbuff_rand_dpr32();
+
+};
+
+#include <amsculib2/amscuarray_dops_impl.hpp>
+
+#endif
+

include/amsculib3/amscuarray_dops_impl.hpp

@@ -0,0 +1,404 @@
+#ifndef __AMSCUARRAY_DOPS_IMPL_HPP__
+#define __AMSCUARRAY_DOPS_IMPL_HPP__
+
+namespace amscuda
+{
+
+template<typename T> __global__ void dbuff_sum_kf(T *devbuffer, int N, T *rets)
+{
+    int I0 = threadIdx.x + blockIdx.x*blockDim.x;
+    int Is = blockDim.x*gridDim.x;
+    int I;
+
+    T ret = (T) 0;
+    for(I=I0;I<N;I=I+Is)
+    {
+        ret = ret + devbuffer[I];
+    }
+    rets[I0] = ret;
+}
+
+template<typename T> T devcuarray_sum(cuarray<T> *devptr)
+{
+    T ret = T();
+    cudaError_t err = cudaSuccess;
+
+    cuarray<T> ldptr;
+    
+    cudaMemcpy(&ldptr,devptr,sizeof(cuarray<T>),cudaMemcpyDeviceToHost);
+    
+    ret = devbuffer_sum(ldptr.data,ldptr.length);
+
+    ldptr.data = NULL;
+    ldptr.length=0;
+
+    return ret;
+}
+
+template<typename T> T dbuff_sum(T *dbuff, int N)
+{
+    int I;
+    T ret = T();
+    cudaError_t err = cudaSuccess;
+
+    int nblocks;
+    int nthreads;
+
+    if(dbuff==NULL || N<=0)
+    {
+        return ret;
+    }
+
+    if(N>100)
+    {
+        nblocks = 10;
+        nthreads = (int)sqrt((float) (N/nblocks));
+        if(nthreads<=0) nthreads=1;
+        if(nthreads>512) nthreads=512;
+    }
+    else
+    {
+        nblocks = 1;
+        nthreads = 1;
+    }
+
+    T *rets = NULL;
+    T *devrets = NULL;
+
+    rets = new T[nblocks*nthreads];
+    cudaMalloc(&devrets,sizeof(T)*nblocks*nthreads);
+
+    dbuff_sum_kf<<<nblocks,nthreads>>>(dbuff,N,devrets);
+    cudaDeviceSynchronize();
+    err = cudaGetLastError();
+    if(err!=cudaSuccess)
+    {
+        printf("amscu::dbuff_sum error: %s\n",cudaGetErrorString(err));
+    }
+
+    cudaMemcpy(rets,devrets,sizeof(T)*nblocks*nthreads,cudaMemcpyDeviceToHost);
+
+    ret = (T)0;
+    for(I=0;I<nblocks*nthreads;I++)
+    {
+        ret = ret + rets[I];
+    }
+
+    cudaFree(devrets); devrets = NULL;
+    delete[] rets;
+
+    return ret;
+}
+
+
+template<typename T> __global__ void dbuff_minmax_kf(T *devbuffer, int N, T *maxs, T *mins)
+{
+    int I0 = threadIdx.x + blockIdx.x*blockDim.x;
+    int Is = blockDim.x*gridDim.x;
+    int I;
+
+    for(I=I0;I<N;I=I+Is)
+    {
+        if(I==I0)
+        {
+            maxs[I0] = devbuffer[I];
+            mins[I0] = devbuffer[I];
+        }
+        else
+        {
+            if(devbuffer[I]>maxs[I0])
+            {
+                maxs[I0] = devbuffer[I];
+            }
+            if(devbuffer[I]<mins[I0])
+            {
+                mins[I0] = devbuffer[I];
+            }
+        }
+    }
+
+    return;
+}
+
+template<typename T> void dbuff_minmax(T *devbuffer, int N, T *min, T *max)
+{
+    cudaError_t err = cudaSuccess;
+    int nblocks;
+    int nthreads;
+    int I;
+
+    T *maxs = NULL;
+    T *dev_maxs = NULL;
+    T *mins = NULL;
+    T *dev_mins = NULL;
+
+    T localmax = T(0);
+    T localmin = T(0);
+
+    if(devbuffer==NULL || N<=0)
+    {
+        if(min!=NULL) *min = T(0);
+        if(max!=NULL) *max = T(0);
+        return;
+    }
+
+    if(N>25)
+    {
+        nblocks = 25;
+        nthreads = (int) sqrt((float)(N/nblocks));
+        if(nthreads<1) nthreads = 1;
+        if(nthreads>512) nthreads = 512;
+    }
+    else
+    {
+        nblocks = 1;
+        nthreads = 1;
+    }
+
+    maxs = new T[nblocks*nthreads];
+    mins = new T[nblocks*nthreads];
+    cudaMalloc(&dev_maxs,nblocks*nthreads);
+    cudaMalloc(&dev_mins,nblocks*nthreads);
+    
+    dbuff_minmax_kf<<<nblocks,nthreads>>>(devbuffer,N,dev_maxs,dev_mins);
+    cudaDeviceSynchronize();
+    err = cudaGetLastError();
+    if(err!=cudaSuccess)
+    {
+        printf("amscu::dbuff_minmax error: %s\n",cudaGetErrorString(err));
+    }
+
+    cudaMemcpy(maxs,dev_maxs,sizeof(T)*nblocks*nthreads,cudaMemcpyDeviceToHost);
+    cudaMemcpy(mins,dev_mins,sizeof(T)*nblocks*nthreads,cudaMemcpyDeviceToHost);
+    
+
+    for(I=0;I<nblocks*nthreads;I++)
+    {
+        if(I==0)
+        {
+            localmax = maxs[0];
+            localmin = mins[0];
+        }
+        else
+        {
+            if(maxs[I]>localmax) localmax = maxs[I];
+            if(mins[I]<localmin) localmin = mins[I];
+        }
+    }
+
+    if(max!=NULL) *max = localmax;
+    if(min!=NULL) *min = localmin;
+
+    cudaFree(dev_maxs); dev_maxs = NULL;
+    cudaFree(dev_mins); dev_mins = NULL;
+    delete[] maxs; maxs = NULL;
+    delete[] mins; mins = NULL;
+
+    return;
+}
+
+template<typename T> __global__ void dbuff_setall_kf(T *devbuffer, int N, T setto)
+{
+    int I0 = threadIdx.x + blockIdx.x*blockDim.x;
+    int Is = blockDim.x*gridDim.x;
+    int I;
+
+    for(I=I0;I<N;I=I+Is)
+    {
+        devbuffer[I] = setto;
+    }
+    return;
+}
+
+template<typename T> void dbuff_setall(T *devbuffer, int N, T setto, int nblocks, int nthreads)
+{
+    cudaError_t err = cudaSuccess;
+
+    if(devbuffer==NULL || N<=0)
+    {
+        return;
+    }
+
+    dbuff_setall_kf<<<nblocks,nthreads>>>(devbuffer,N,setto);
+    cudaDeviceSynchronize();
+    err = cudaGetLastError();
+    if(err!=cudaSuccess)
+    {
+        printf("amscu::dbuff_setall error: %s\n",cudaGetErrorString(err));
+    }
+
+    return;
+}
+
+template<typename T1, typename T2, typename T3> __global__ void dbuff_vectorbinop_kf1(T1 *dbuf_a, T2 *dbuf_b, T3 *dbuf_out, int N, T3 (*fpnt)(T1,T2))
+{
+    int I0 = threadIdx.x + blockIdx.x*blockDim.x;
+    int Is = blockDim.x*gridDim.x;
+    int I;
+
+    T1 a;
+    T2 b;
+    T3 c;
+
+    for(I=I0;I<N;I=I+Is)
+    {
+        a = dbuf_a[I];
+        b = dbuf_b[I];
+        c = fpnt(a,b);
+        dbuf_out[I] = c;
+    }
+
+    return;
+}
+
+template<typename T1, typename T2, typename T3> __global__ void dbuff_vectorbinop_kf2(T1 *dbuf_a, T2 par_b, T3 *dbuf_out, int N, T3 (*fpnt)(T1,T2))
+{
+    int I0 = threadIdx.x + blockIdx.x*blockDim.x;
+    int Is = blockDim.x*gridDim.x;
+    int I;
+
+    T1 a;
+    T2 b;
+    T3 c;
+
+    for(I=I0;I<N;I=I+Is)
+    {
+        a = dbuf_a[I];
+        b = par_b;
+        c = fpnt(a,b);
+        dbuf_out[I] = c;
+    }
+
+    return;
+}
+
+
+//Elementwise device-buffer vector binary operation
+//takes two input arrays ( , ) --> one output array
+template<typename T1, typename T2, typename T3> void dbuff_vectorbinop(T1 *dbuf_a, T2 *dbuf_b, T3 *dbuf_out, int N, T3 (*fpnt)(T1,T2), int nblocks, int nthreads)
+{
+    cudaError_t err = cudaSuccess;
+
+    if(dbuf_a == NULL || dbuf_b == NULL || dbuf_out == NULL || N<=0)
+    {
+        return;
+    }
+
+    dbuff_vectorbinop_kf1<<<nblocks,nthreads>>>(dbuf_a,dbuf_b,dbuf_out,N);
+    cudaDeviceSynchronize();
+    err = cudaGetLastError();
+    if(err!=cudaSuccess)
+    {
+        printf("amscu::devbuffer_vectorbinop error: %s\n",cudaGetErrorString(err));
+    }
+
+    return;
+}
+
+//Elementwise device-buffer vector two-parameter operation
+//takes one input array, and a constant paramter ( ) ---> one output array
+template<typename T1, typename T2, typename T3> void dbuff_vectorbinop(T1 *dbuf_a, T2 par_b, T3 *dbuf_out, int N, T3 (*fpnt)(T1,T2), int nblocks, int nthreads)
+{
+    cudaError_t err = cudaSuccess;
+
+    if(dbuf_a == NULL || dbuf_out == NULL || N<=0)
+    {
+        return;
+    }
+
+    dbuff_vectorbinop_kf2<<<nblocks,nthreads>>>(dbuf_a,par_b,dbuf_out,N);
+    cudaDeviceSynchronize();
+    err = cudaGetLastError();
+    if(err!=cudaSuccess)
+    {
+        printf("amscu::devbuffer_vectorbinop error: %s\n",cudaGetErrorString(err));
+    }
+
+    return;
+}
+
+template<typename T> T dbuff_add_fn(T a, T b)
+{
+    return a+b;
+}
+
+template<typename T> void dbuff_add(T *dbuff_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_a,dbuff_b,dbuff_out,N,&dbuff_add_fn,nblocks,nthreads);
+    return;
+}
+
+template<typename T> void dbuff_add(T *dbuff_a, T par_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_a,par_b,dbuff_out,N,&dbuff_add_fn,nblocks,nthreads);
+    return;
+}
+
+template<typename T> T dbuff_sub_fn(T a, T b)
+{
+    return a-b;
+}
+
+template<typename T> void dbuff_sub(T *dbuff_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_a,dbuff_b,dbuff_out,N,&dbuff_sub_fn,nblocks,nthreads);
+    return;
+}
+
+template<typename T> void dbuff_sub(T *dbuff_a, T par_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_a,par_b,dbuff_out,N,&dbuff_sub_fn,nblocks,nthreads);
+    return;
+}
+
+template<typename T> T dbuff_mult_fn(T a, T b)
+{
+    return a*b;
+}
+
+template<typename T> void dbuff_mult(T *dbuff_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_a,dbuff_b,dbuff_out,N,&dbuff_mult_fn,nblocks,nthreads);
+    return;
+}
+
+template<typename T> void dbuff_mult(T *dbuff_a, T par_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_a,par_b,dbuff_out,N,&dbuff_mult_fn,nblocks,nthreads);
+    return;
+}
+
+template<typename T> T dbuff_div_fn(T a, T b)
+{
+    return a/b;
+}
+
+template<typename T> void dbuff_div(T *dbuff_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_a,dbuff_b,dbuff_out,N,&dbuff_div_fn,nblocks,nthreads);
+    return;
+}
+
+template<typename T> void dbuff_div(T *dbuff_a, T par_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_a,par_b,dbuff_out,N,&dbuff_div_fn,nblocks,nthreads);
+    return;
+}
+
+template<typename T> T dbuff_ldiv_fn(T a, T b)
+{
+    return b/a;
+}
+
+
+template<typename T> void dbuff_div(T par_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_b,par_a,dbuff_out,N,&dbuff_ldiv_fn,nblocks,nthreads);
+    return;
+}
+
+
+};
+
+#endif
+

include/amsculib3/amscuarray_impl.hpp

@@ -0,0 +1,323 @@
+#ifndef __CUARRAY_IMPL_HPP__
+#define __CUARRAY_IMPL_HPP__
+
+namespace amscuda
+{
+
+// New Version cuarray<T>
+// simpler, less crap going on
+
+template<typename T> __device__ __host__ cuarray<T>::cuarray()
+{
+    length = 0;
+    data = NULL;
+}
+
+template<typename T> __device__ __host__ cuarray<T>::~cuarray()
+{
+    if(data!=NULL)
+    {
+        delete[] data; data = NULL;
+    }
+    length = 0;
+}
+
+template<typename T> __device__ __host__ int cuarray<T>::resize(const int _length)
+{
+    int ret = 0;
+
+    T *newbuffer = NULL;
+
+    if(length==_length)
+    {
+        //do nothing
+        ret = 1;
+        return ret;
+    }
+    if(_length<=0)
+    {
+        if(data!=NULL)
+        {
+            delete[] data;
+            data = NULL;
+        }
+        length = 0;
+        ret = 1;
+    }
+
+    newbuffer = new T[_length];
+    if(newbuffer==NULL)
+    {
+        ret = -1; //failed to allocate memory
+        return ret;
+    }
+
+    int I;
+    T def;
+
+    if(data!=NULL)
+    {
+        for(I=0;I<length&&I<_length;I++)
+        {
+            newbuffer[I] = data[I];
+        }
+        for(I=length;I<_length;I++)
+        {
+            newbuffer[I] = def;
+        }
+        delete[] data; data=NULL;
+    }
+    else
+    {
+        for(I=0;I<_length;I++)
+        {
+            newbuffer[I] = def;
+        }
+    }
+
+    data = newbuffer;
+    length = _length;
+    ret = 1;
+
+    return ret;
+}
+
+
+template<typename T> __host__ int cuarray<T>::device_send(cuarray<T> **dptr)
+{
+    int ret = 0;
+    int dlength;
+
+    if(*dptr==NULL)
+    {
+        ret = _device_send_overwrite(dptr);
+    }
+    else
+    {
+        dlength = device_length(*dptr);
+        if(dlength=length)
+        {
+            ret = _device_send_copy(*dptr);
+        }
+        else
+        {
+            ret = _device_send_overwrite(dptr);
+        }
+    }
+
+    return ret;
+}
+
+template<typename T> __host__ int cuarray<T>::_device_send_overwrite(cuarray<T> **dptr)
+{
+    int ret = 0;
+    cuarray<T> dlocal;
+    cudaError_t err = cudaSuccess;
+    device_free(dptr);
+    
+    if(length>=0 && data!=NULL)
+    {
+        err = cudaMalloc(dptr,sizeof(cuarray<T>));
+        if(err==cudaSuccess)
+        {
+            err = cudaMalloc(&(dlocal.data),sizeof(T)*length);
+            dlocal.length = length;
+
+            if(err==cudaSuccess)
+            {
+                cudaMemcpy(*dptr,&dlocal,sizeof(cuarray<T>),cudaMemcpyHostToDevice);
+                if(data!=NULL)
+                    err = cudaMemcpy(dlocal.data,data,sizeof(T)*length,cudaMemcpyHostToDevice);
+                else
+                    err = cudaSuccess;
+                if(err==cudaSuccess)
+                {
+                    ret = 1;
+                }
+                else
+                {
+                    ret = -3;
+                }
+            }
+            else
+            {
+                ret = -2;
+            }
+        }
+        else
+        {
+            ret = -1;
+        }
+    }
+    else
+    {
+        dlocal.data = NULL;
+        dlocal.length = 0;
+        err = cudaMalloc(dptr,sizeof(cuarray<T>));
+        if(err==cudaSuccess)
+        {
+            cudaMemcpy(*dptr,&dlocal,sizeof(cuarray<T>),cudaMemcpyHostToDevice);
+            ret = 1;
+        }
+        else
+        {
+            ret = -4;
+        }
+    }
+
+
+    dlocal.data = NULL;
+    dlocal.length = -1;
+
+    return ret;
+}
+
+template<typename T> __host__ int cuarray<T>::_device_send_copy(cuarray<T> *dptr)
+{
+    int ret = 0;
+    cudaError_t err = cudaSuccess;
+    T* ddata = NULL;
+    ddata = device_data_ptr(dptr);
+
+    err = cudaMemcpy(ddata,data,sizeof(T)*length,cudaMemcpyHostToDevice);
+    if(err==cudaSuccess)
+    {
+        ret = 1;
+    }
+    else
+    {
+        ret = -1;
+    }
+
+    return ret;
+}
+
+template<typename T> __host__ int cuarray<T>::device_pull(cuarray<T> *dptr)
+{
+    int ret = 0;
+    int dlength;
+    T* ddata;
+    cudaError_t err;
+
+    if(dptr==NULL) 
+    {
+        ret = -1; // null d pointer
+        return ret;
+    }
+
+    dlength = device_length(dptr);
+    if(dlength!=length)
+    {
+        this->resize(dlength);
+    }
+
+    ddata = device_data_ptr(dptr);
+
+    if(length>0 && data!=NULL && ddata!=NULL)
+    {
+        err = cudaMemcpy(data,dptr,length*sizeof(T),cudaMemcpyDeviceToHost);
+        if(err==cudaSuccess)
+        {
+            ret = 1;
+        }
+        else
+        {
+            ret = -2;
+        }
+    }
+
+    return ret;
+}
+
+template<typename T> __host__ int cuarray<T>::device_free(cuarray<T> **dptr)
+{
+    int ret = 0;
+    cuarray<T> dlocal;
+
+    if(*dptr!=NULL)
+    {
+        cudaMemcpy(&dlocal,dptr,sizeof(cuarray<T>),cudaMemcpyDeviceToHost);
+        if(dlocal.data!=NULL)
+        {
+            cudaFree(dlocal.data);
+            dlocal.data = NULL;
+        }
+
+        cudaFree(*dptr);
+        *dptr = NULL;
+        ret = 1;
+    }
+
+    dlocal.data = NULL;
+    dlocal.length = -1;
+
+    return ret;
+}
+
+template<typename T> __host__ int cuarray<T>::device_length(cuarray<T> *dptr)
+{
+    int ret = -1;
+    cuarray<T> dlocal;
+
+    if(dptr==NULL)
+    {
+        return ret;
+    }
+
+    cudaMemcpy(&dlocal,dptr,sizeof(cuarray<T>),cudaMemcpyDeviceToHost);
+    ret = dlocal.length;
+
+    dlocal.data = NULL;
+    dlocal.length = -1;
+
+    return ret;
+}
+
+template<typename T> __host__ T* cuarray<T>::device_data_ptr(cuarray<T> *dptr)
+{
+    T* ret = NULL;
+    cuarray<T> dlocal;
+
+    if(dptr==NULL)
+    {
+        return ret;
+    }
+
+    cudaMemcpy(&dlocal,dptr,sizeof(cuarray<T>),cudaMemcpyDeviceToHost);
+    ret = dlocal.data;
+
+    dlocal.data = NULL;
+    dlocal.length = -1;
+
+    return ret;
+}
+
+template<typename T>  __device__ __host__ int cuarray<T>::size() const
+{
+    return this->length;
+}
+
+template<typename T>  __device__ __host__ T& cuarray<T>::at(const int I)
+{
+    return this->data[I];
+}
+
+template<typename T>  __device__ __host__ const T& cuarray<T>::at(const int I) const
+{
+    return this->data[I];
+}
+
+template<typename T> __device__ __host__ T& cuarray<T>::operator[](const int I)
+{
+    return this->data[I];
+}
+
+template<typename T> __device__ __host__ const T& cuarray<T>::operator[](const int I) const
+{
+    return this->data[I];
+}
+
+};
+
+
+#endif

include/amsculib3/amscuda_binarrrw.hpp

@@ -0,0 +1,19 @@
+#ifndef __AMSCUDA_BINARRRW_HPP__
+#define __AMSCUDA_BINARRRW_HPP__
+
+namespace amscuda
+{
+
+template<typename T> int fread_ndarray(FILE *fp, cuarray<int32_t> *shape, cuarray<T> *buffer);
+template<typename T> int fwrite_ndarray(FILE *fp, const cuarray<int32_t> *shape, const cuarray<T> *buffer);
+
+template<typename T> int fwrite_buffer(FILE *fp, const int N, const T *buffer);
+template<typename T> int fread_buffer(FILE *fp, const int Nmax, const T *buffer);
+
+
+}; //end namespace amscuda
+
+#include <amsculib2/amscuda_binarrrw_impl.hpp>
+
+#endif
+

include/amsculib3/amscuda_binarrrw_impl.hpp

@@ -0,0 +1,194 @@
+#ifndef __AMSCUDA_BINARRRW_IMPL_HPP__
+#define __AMSCUDA_BINARRRW_IMPL_HPP__
+
+namespace amscuda
+{
+
+template<typename T> int fread_ndarray(FILE *fp, cuarray<int32_t> *shape, cuarray<T> *buffer)
+{
+    int ret = 1;
+    int I;
+    long piprod;
+    int32_t q;
+    int cnt;
+
+    int32_t Nd;
+
+    if(fp!=NULL)
+    {
+        if(!feof(fp))
+        {
+            cnt = fread(&Nd,sizeof(int32_t),1,fp);
+            if(Nd>0 && cnt>0)
+            {
+                shape->resize(Nd);
+                piprod = 1;
+                for(I=0;I<Nd;I++)
+                {
+                    cnt = fread(&q,sizeof(int32_t),1,fp);
+                    shape->at(I) = q;
+                    if(q>0)
+                    {
+                        piprod = piprod*q;
+                    }
+                    else
+                    {
+                        piprod = 0;
+                    }
+                }
+
+                buffer->resize(piprod);
+                if(piprod>0)
+                {
+                    cnt = fread((buffer->data),sizeof(T),piprod,fp);
+                    if(piprod==cnt)
+                    {
+                        ret = 1;
+                    }
+                    else
+                    {
+                        printf("fread_ndarray, read %d values, expecting %ld\n",cnt,piprod);
+                        ret = 0;
+                    }
+                }
+            }
+            else
+            {
+                printf("fread_ndarray: Read a number of dimensions<=0.\n");
+                Nd = 0;
+                shape->resize(0);
+                buffer->resize(0);
+            }
+        }
+        else
+        {
+            printf("fread_ndarray: fp=NULL.\n");
+            ret = 0;
+        }
+    }
+    else
+    {
+        ret = 0;
+    }
+
+    return ret;
+}
+
+template<typename T> int fwrite_ndarray(FILE *fp, const cuarray<int32_t> *shape, const cuarray<T> *buffer)
+{
+    int ret = 1;
+    long piprod;
+    int I;
+    int32_t Nd;
+
+    if(fp==NULL)
+    {
+        ret = 0;
+        printf("fwrite_ndarray: fp=NULL\n");
+        return ret;
+    }
+
+    piprod = 1;
+    for(I=0;I<shape->size();I++)
+    {
+        if(shape->at(I)>0)
+        {
+            piprod = piprod*shape->at(I);
+        }
+        else
+        {
+            piprod = 0;
+        }
+    }
+
+    Nd = (int32_t) shape->size();
+
+    if(piprod!=buffer->size())
+    {
+        ret = 0;
+        printf("fwrite_ndarray: buffer is size %ld, while shape is size %ld\n",(long)buffer->size(),(long)piprod);
+        return ret;
+    }
+
+    fwrite(&Nd,sizeof(int32_t),1,fp);
+    if(Nd>0)
+    {
+        fwrite(shape->data,sizeof(int32_t),Nd,fp);
+        if(piprod>0)
+        {
+            fwrite(buffer->data,sizeof(T),buffer->size(),fp);
+        }
+    }
+
+    return ret;
+}
+
+template<typename T> int fwrite_buffer(FILE *fp, const int N, const T *buffer)
+{
+    int ret = 0;
+    int Nd = 1;
+    
+    if(fp==NULL)
+    {
+        ret = 0;
+        printf("fwrite_buffer: fp=NULL\n");
+        return ret;
+    }
+
+    fwrite(&Nd,sizeof(int32_t),1,fp);
+    fwrite(&N,sizeof(int32_t),1,fp);
+    fwrite(buffer,sizeof(T),N,fp);
+    
+    return ret;
+}
+
+template<typename T> int fread_buffer(FILE *fp, const int Nmax, const T *buffer)
+{
+    int ret = 0;
+    
+    int cnt;
+    int32_t Nd;
+    int32_t *dims = NULL;
+    int piprod;
+    int32_t q;
+    int I;
+
+    int Nr;
+
+
+    if(fp==NULL) {ret = -1; return ret;}
+    if(feof(fp)) {ret = -2; return ret;}
+
+    cnt = fread(&Nd,sizeof(int32_t),1,fp);
+    if(Nd>0 && cnt>0)
+    {
+        piprod = 1;
+        dims = new(std::nothrow) int32_t[Nd];
+        for(I=0;I<Nd;I++)
+        {
+            cnt = fread(&q,sizeof(int32_t),1,fp);
+            dims[I] = q;
+            piprod = piprod*dims[I];
+            if(piprod==cnt)
+            {
+                ret = 1;
+            }
+            else
+            {
+                printf("fwrite_buffer, read %d values, expecting %d\n",cnt,piprod);
+            }
+        }
+
+        Nr = amscuda::min<int32_t>(Nmax,piprod);
+        cnt = fread(buffer,sizeof(T),Nr,fp);
+    }
+
+    if(dims!=NULL) {delete[] dims; dims=NULL;}
+
+    return ret;
+}
+
+}; //end namespace amscuda
+
+#endif
+

include/amsculib3/amscugeom.hpp

@@ -0,0 +1,11 @@
+#ifndef __AMSCUGEOM_HPP__
+#define __AMSCUGEOM_HPP__
+
+namespace amscuda
+{
+
+
+}; //end namespace amscuda
+
+#endif
+

include/amsculib3/amsculib3.hpp

@@ -0,0 +1,63 @@
+#ifndef __AMSCULIB3_HPP__
+#define __AMSCULIB3_HPP__
+
+//Std Lib Includes
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <stdint.h>
+#include <time.h>
+#include <new>
+
+#include <cuda_runtime_api.h> //where all the cuda functions live
+#include <cuda_runtime.h>
+#include <cuda.h>
+
+//Dependencies
+
+//Predeclarations
+class cuvect2;
+class cuvect3;
+class cuvect4;
+class cuvect2f;
+class cuvect3f;
+class cuvect4f;
+
+//Need a way to define the same symbols using both host and device code
+//A solution was found here: https://stackoverflow.com/questions/9457572/cuda-host-and-device-using-same-constant-memory
+#ifdef __CUDA_ARCH__
+#define AMSCU_CONST __constant__
+#else
+#define AMSCU_CONST
+#endif
+
+namespace amscuda
+{
+
+    //default thread and block execution
+    AMSCU_CONST static const int amscu_defnblocks = 256;
+    AMSCU_CONST static const int amscu_defnthreads = 512;
+
+    //default numthreads to execute on cpu
+    AMSCU_CONST static const int amscu_defcputhreads = 8;
+
+}; //end namespace amscuda
+
+//Components
+#include <amsculib3/amscu_cudafunctions.hpp>
+#include <amsculib3/math/amscumath.hpp>
+
+#include <amsculib3/amscugeom.hpp>
+#include <amsculib3/amscuarray.hpp>
+#include <amsculib3/amscuda_binarrrw.hpp>
+#include <amsculib3/amscu_random.hpp>
+
+#include <amsculib3/amscuarray_dops.hpp>
+
+#include <amsculib3/amscurarray.cuh>
+
+
+
+
+#endif
+

include/amsculib3/amscurarray.cuh

@@ -0,0 +1,66 @@
+#ifndef __AMSCURARRAY_HPP__
+#define __AMSCURARRAY_HPP__
+
+namespace amscuda
+{
+
+//Cuda ragged array class
+template<typename T> class curarray
+{
+public:
+    int device;
+    curarray* devptr; //pointer to mirror class on the device
+
+    int Narrays; //number of arrays
+    
+    int *N; //dimension of each array
+    T** hostarrayptrs; //pointers to each array on the host - null on the device
+    T** devarrayptrs; //pointers to each array on the device
+        //the double pointer is a host pointer to device pointers on the host class
+        //for the device class, only the second set of arrays is in use
+    
+    //the constructor and destructor set all pointers to NULL, they
+    // do *not* manage memory. This is done with curarray_new and curarray_delete
+    __device__ __host__ curarray();
+    __device__ __host__ ~curarray();
+    
+    __host__ int push();
+    __host__ int pull();
+    //__device__ int dev_resizearray(int arraynum, int arraysize);
+    __host__ int resizearray(int arraynum, int arraysize);
+    // I may want a way to resize arrays on the device without pushing/pulling all the array contents
+    
+
+};
+
+template<typename T> int curarray_new(curarray<T>** ptr, int Narrays);
+
+template<typename T> int curarray_delete(curarray<T>** ptr);
+
+template<typename T> int curarray_device_new(curarray<T> *hostptr);
+
+template<typename T> int curarray_device_delete(curarray<T> *hostptr);
+
+template<typename T> int curarray_push(curarray<T> *hostptr);
+
+template<typename T> int curarray_pull(curarray<T> *hostptr);
+
+
+//template<typename T> int curarray_host_fillall(curarray<T> *hostptr, const T &val);
+//template<typename T> int curarray_device_fillall(curarray<T> *hostptr, const T &val);
+
+
+//template<typename T> __host__ int curarray_deletearray(curarray<T> *hostptr, int arrayindex);
+//template<typename T> __device__ int curarray_dev_deletearray(curarray<T> *devptr, int arrayindex);
+
+//template<typename T> __host__ int curarray_allocarray(curarray<T> *hostptr, int arrayindex, int size);
+//template<typename T> __device__ int curarray_dev_allocarray(curarray<T> *devptr, int arrayindex, int size);
+
+
+void test_amscurarray1();
+
+};
+
+#include <amsculib2/amscurarray_impl.cuh>
+
+#endif

include/amsculib3/amscurarray_impl.cuh

@@ -0,0 +1,529 @@
+#ifndef __AMSCURARRAY_IMPL_HPP__
+#define __AMSCURARRAY_IMPL_HPP__
+
+namespace amscuda
+{
+
+template<typename T> curarray<T>::curarray()
+{
+    device = -1;
+    devptr = NULL;
+    Narrays = 0;
+    N = NULL;
+    hostarrayptrs = NULL;
+    devarrayptrs = NULL;
+
+}
+
+template<typename T> curarray<T>::~curarray()
+{
+    device = -1;
+    devptr = NULL;
+    Narrays = 0;
+    N = NULL;
+    hostarrayptrs = NULL;
+    devarrayptrs = NULL;
+    
+}
+
+template<typename T> int curarray_new(curarray<T>** ptr, int Narrays)
+{
+    int ret = 0;
+    int device;
+    curarray<T> *lhptr = *ptr;
+
+    cudaGetDevice(&device);
+
+    if(lhptr!=NULL)
+    {
+        curarray_delete(ptr);
+    }
+
+    *ptr = new(std::nothrow) curarray<T>();
+    lhptr = *ptr;
+
+    int I;
+
+    if(Narrays<0) Narrays=0;
+
+    lhptr->Narrays = Narrays;
+    lhptr->device = device;
+    lhptr->N = new(std::nothrow) int[Narrays];
+    lhptr->hostarrayptrs = new(std::nothrow) T*[Narrays];
+    lhptr->devarrayptrs = new(std::nothrow) T*[Narrays];
+    
+    for(I=0;I<Narrays;I++)
+    {
+        lhptr->N[I] = 0;
+        lhptr->hostarrayptrs[I] = NULL;
+        lhptr->devarrayptrs[I] = NULL;
+    }
+
+    curarray_device_new(lhptr);
+
+    return ret;
+}
+
+template<typename T> int curarray_delete(curarray<T>** ptr)
+{
+    int ret = 0;
+    curarray<T> *lptr = NULL;
+    int olddev;
+
+    int I;
+
+    if(*ptr==NULL)
+    {
+        return 0;
+    }
+    lptr = *ptr;
+
+    cudaGetDevice(&olddev);
+    cudaSetDevice(lptr->device);
+
+    
+    if(lptr->devptr!=NULL)
+    {
+        curarray_device_delete(lptr);
+    }
+
+    lptr->device = -1;
+
+    for(I=0;I<lptr->Narrays;I++)
+    {
+        if(lptr->hostarrayptrs!=NULL)
+        {
+            if(lptr->hostarrayptrs[I]!=NULL)
+            {
+                delete[] lptr->hostarrayptrs[I];
+                lptr->hostarrayptrs[I] = NULL;
+            }
+        }
+
+        if(lptr->devarrayptrs!=NULL)
+        {
+            if(lptr->devarrayptrs[I]!=NULL)
+            {
+                //erasing device memory should have been handled in curarray_device_delete
+                lptr->devarrayptrs[I] = NULL;
+            }
+        }
+
+        lptr->N[I] = 0;
+    }
+
+    if(lptr->N != NULL) {delete[] lptr->N; lptr->N = NULL;}
+    if(lptr->hostarrayptrs!=NULL) {delete[] lptr->hostarrayptrs; lptr->hostarrayptrs=NULL;}
+    if(lptr->devarrayptrs!=NULL) {delete[] lptr->devarrayptrs; lptr->devarrayptrs=NULL;}
+
+    if(*ptr!=NULL) {delete *ptr; *ptr = NULL;}
+
+    cudaSetDevice(olddev);
+
+    return ret;
+}
+
+template<typename T> int curarray_device_new(curarray<T> *hostptr)
+{
+    int ret = 0;
+    curarray<T> ldevdata;
+
+    if(hostptr==NULL) return -1;
+    if(hostptr->devptr!=NULL)
+    {
+        curarray_device_delete(hostptr);
+    }
+
+    int I;
+
+    cudaGetDevice(&(hostptr->device));
+
+    ldevdata.device = hostptr->device;
+    ldevdata.Narrays = hostptr->Narrays;
+    int Narrays = hostptr->Narrays;
+
+    cudaMalloc(&(ldevdata.N),sizeof(int)*Narrays);
+    cudaMemcpy(ldevdata.N,hostptr->N,sizeof(int)*Narrays,cudaMemcpyHostToDevice);
+
+    ldevdata.hostarrayptrs = NULL;
+
+    for(I=0;I<Narrays;I++)
+    {
+        if(hostptr->N[I]>0)
+        {
+            if(hostptr->devarrayptrs[I]!=NULL)
+            {
+                cudaFree(hostptr->devarrayptrs[I]);
+                hostptr->devarrayptrs[I] = NULL;
+            }
+
+            cudaMalloc(&(hostptr->devarrayptrs[I]),sizeof(T)*hostptr->N[I]);
+            cudaMemcpy(hostptr->devarrayptrs[I],hostptr->hostarrayptrs[I],sizeof(T)*hostptr->N[I],cudaMemcpyHostToDevice);
+        }
+        else
+        {
+            if(hostptr->devarrayptrs[I]!=NULL)
+            {
+                cudaFree(hostptr->devarrayptrs[I]);
+                hostptr->devarrayptrs[I] = NULL;
+            }
+        }
+    }
+
+    cudaMalloc(&(ldevdata.devarrayptrs),sizeof(T*)*Narrays);
+    cudaMemcpy(ldevdata.devarrayptrs,hostptr->devarrayptrs,sizeof(T*)*Narrays,cudaMemcpyHostToDevice);
+
+    cudaMalloc(&(hostptr->devptr),sizeof(curarray<T>));
+    cudaMemcpy(hostptr->devptr,&ldevdata,sizeof(curarray<T>),cudaMemcpyHostToDevice);    
+    
+    ret = 1;
+
+    return ret;
+}
+
+template<typename T> int curarray_device_delete(curarray<T> *hostptr)
+{
+    int ret = 0;
+
+    curarray<T> ldevdata;
+    int olddev;
+
+    if(hostptr->devptr==NULL)
+    {
+        return 0;
+    }
+
+    cudaGetDevice(&olddev);
+    cudaSetDevice(hostptr->device);
+
+    cudaMemcpy(&ldevdata,hostptr->devptr,sizeof(curarray<T>),cudaMemcpyDeviceToHost);
+
+    int I;
+    for(I=0;I<hostptr->Narrays;I++)
+    {
+        if(hostptr->devarrayptrs[I]!=NULL)
+        {
+            cudaFree(hostptr->devarrayptrs[I]);
+            hostptr->devarrayptrs[I] = NULL;
+        }
+    }
+
+    cudaFree(ldevdata.devarrayptrs);
+    cudaFree(ldevdata.N);
+
+    cudaFree(hostptr->devptr);
+    hostptr->devptr = NULL;
+    hostptr->device = -1;
+
+    cudaSetDevice(olddev);
+
+    ret = 1;
+    return ret;
+}
+
+template<typename T> int curarray_push(curarray<T> *hostptr)
+{
+    int ret = 0;
+
+    int olddev;
+
+    curarray<T> ldevdata;
+    T** ldevarrayptrs = NULL;
+    int *devN = NULL;
+
+    if(hostptr==NULL) return -1;
+
+    cudaGetDevice(&olddev);
+    cudaSetDevice(hostptr->device);
+
+    int Narrays = hostptr->Narrays;
+
+    cudaMemcpy(&ldevdata,hostptr->devptr,sizeof(curarray<T>),cudaMemcpyDeviceToHost);
+    ldevarrayptrs = new(std::nothrow) T*[Narrays];
+    devN = new(std::nothrow) int[Narrays];
+
+    cudaMemcpy(ldevarrayptrs,ldevdata.devarrayptrs,sizeof(T*)*Narrays,cudaMemcpyDeviceToHost);
+    cudaMemcpy(devN,ldevdata.N,sizeof(int)*Narrays,cudaMemcpyDeviceToHost);
+
+
+    int I;
+
+    for(I=0;I<Narrays;I++)
+    {
+        //check to see that host size is the same as device size, and that 
+        //the host device pointer is the same as the device device pointer
+        if( (hostptr->N[I]!=devN[I]) || 
+            (hostptr->devarrayptrs[I] != ldevarrayptrs[I])
+        )
+        {
+            cudaFree(ldevarrayptrs[I]);
+            ldevarrayptrs[I] = NULL;
+            hostptr->devarrayptrs[I] = NULL;
+
+            if(hostptr->N[I]>0)
+            {
+                cudaMalloc(&(hostptr->devarrayptrs[I]),sizeof(T)*hostptr->N[I]);
+                ldevarrayptrs[I] = hostptr->devarrayptrs[I];
+                devN[I] = hostptr->N[I];
+            }
+            else
+            {
+                devN[I] = 0;
+            }
+        }
+        if(hostptr->N[I]>0)
+        {
+            //copy host data to device
+            cudaMemcpy(hostptr->devarrayptrs[I],hostptr->hostarrayptrs[I],sizeof(T)*hostptr->N[I],cudaMemcpyHostToDevice);
+        }
+    } //for each array
+
+    //rectify and copy device data structure to device
+    ldevdata.device = hostptr->device;
+    ldevdata.devptr = NULL;
+    ldevdata.Narrays = hostptr->Narrays; //later - logic for dealing with when this is not true
+    ldevdata.hostarrayptrs = NULL;
+
+
+    cudaMemcpy(ldevdata.N,hostptr->N,sizeof(int)*Narrays,cudaMemcpyHostToDevice);
+    cudaMemcpy(ldevdata.devarrayptrs,hostptr->devarrayptrs,sizeof(T*)*Narrays,cudaMemcpyHostToDevice);
+
+    cudaMemcpy(hostptr->devptr,&ldevdata,sizeof(curarray<T>),cudaMemcpyHostToDevice);
+
+    cuda_errortrap("curarray_push cuda error:");
+
+    cudaSetDevice(olddev);
+
+    delete[] ldevarrayptrs;
+    delete[] devN;
+
+    return ret;
+}
+
+template<typename T> int curarray_pull(curarray<T> *hostptr)
+{
+    int ret = 0;
+
+    int olddev;
+
+    curarray<T> ldevdata;
+    T** ldevarrayptrs = NULL;
+    int *devN = NULL;
+
+    if(hostptr==NULL) return -1;
+
+    cudaGetDevice(&olddev);
+    cudaSetDevice(hostptr->device);
+
+    cuda_errortrap("dbg1");
+
+    int Narrays = hostptr->Narrays;
+
+    cudaMemcpy(&ldevdata,hostptr->devptr,sizeof(curarray<T>),cudaMemcpyDeviceToHost);
+    ldevarrayptrs = new(std::nothrow) T*[Narrays];
+    devN = new(std::nothrow) int[Narrays];
+
+    cuda_errortrap("dbg2");
+
+    cudaMemcpy(ldevarrayptrs,ldevdata.devarrayptrs,sizeof(T*)*Narrays,cudaMemcpyDeviceToHost);
+    cudaMemcpy(devN,ldevdata.N,sizeof(int)*Narrays,cudaMemcpyDeviceToHost);
+
+    cuda_errortrap("dbg3");
+    char dbgjnk[50];
+
+    int I;
+    for(I=0;I<Narrays;I++)
+    {
+        //check to see that host size is the same as device size, and that 
+        //the host device pointer is the same as the device device pointer
+        if(hostptr->devarrayptrs[I] != ldevarrayptrs[I])
+        {
+            hostptr->devarrayptrs[I] = ldevarrayptrs[I];
+        }
+
+        if(hostptr->N[I]!=devN[I])
+        {
+            if(hostptr->hostarrayptrs[I]!=NULL) 
+            {
+                delete[] hostptr->hostarrayptrs[I];
+                hostptr->hostarrayptrs[I] = NULL;
+            }
+
+            if(devN[I]>0)
+            {
+                hostptr->hostarrayptrs[I] = new(std::nothrow) T[devN[I]];
+                hostptr->N[I] = devN[I];
+            }
+            else
+            {
+                hostptr->N[I] = 0;
+            }
+        }
+
+        if(hostptr->hostarrayptrs[I]!=NULL && hostptr->devarrayptrs[I]!=NULL)
+        {
+            cudaMemcpy(hostptr->hostarrayptrs[I],hostptr->devarrayptrs[I],sizeof(T)*hostptr->N[I],cudaMemcpyDeviceToHost);
+            sprintf(dbgjnk,"%d dbg %d",I,hostptr->N[I]);
+            cuda_errortrap(dbgjnk);
+        }
+    } //for each array
+
+    //for the pull operation, I don't think any update of the device data structure is necessary
+
+    cudaSetDevice(olddev);
+
+    delete[] ldevarrayptrs;
+    delete[] devN;
+
+    return ret;
+}
+
+template<typename T> __host__ int curarray<T>::push()
+{
+    return curarray_push(this);
+}
+template<typename T> __host__ int curarray<T>::pull()
+{
+    return curarray_pull(this);
+}
+
+
+
+/*
+https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#memory-allocation-and-lifetime%5B/url%5D
+
+cudaMalloc() and cudaFree() have distinct semantics between the host and 
+device environments. When invoked from the host, cudaMalloc() allocates a 
+new region from unused device memory. When invoked from the device runtime 
+these functions map to device-side malloc() and free(). This implies that 
+within the device environment the total allocatable memory is limited to the 
+device malloc() heap size, which may be smaller than the available unused 
+device memory. Also, it is an error to invoke cudaFree() from the host 
+program on a pointer which was allocated by cudaMalloc() on the device 
+or vice-versa.
+
+So, basically this entire function is not going to work. I'll be unable to resize within 
+a kernel.
+*/
+
+/*
+
+template<typename T> __device__ int curarray<T>::dev_resizearray(int arraynum, int arraysize)
+{
+    int ret = 0;
+    T* newptr = NULL;
+    int I;
+    T def;
+
+    if(arraynum>=0 && arraynum<Narrays)
+    {
+        if(N[arraynum]!=arraysize)
+        {
+            if(arraysize<=0)
+            {
+                if(devarrayptrs[arraynum]!=NULL) cudaFree(devarrayptrs[arraynum]);
+                devarrayptrs[arraynum] = NULL;
+                N[arraynum] = 0;
+                ret = 1;
+                return ret;
+            }
+
+            cudaMalloc(&newptr,arraysize*sizeof(T));
+            if(newptr!=NULL)
+            {
+                //do I want to assume there is a copy operator? (operator=)
+                //for now, yes - write a more restrictive class later if I don't want it
+                if(devarrayptrs[arraynum]!=NULL)
+                {
+                    for(I=0;I<N[arraynum]&&I<arraysize;I++)
+                    {
+                        newptr[I] = devarrayptrs[arraynum][I];
+                    }
+                }
+                for(I=N[arraynum];I<arraysize;I++)
+                {
+                    newptr[I] = def;
+                }
+
+                if(devarrayptrs[arraynum]!=NULL) cudaFree(devarrayptrs[arraynum]);
+                devarrayptrs[arraynum] = newptr;
+                N[arraynum] = arraysize;
+                ret = 1;
+            }
+            else
+            {
+                ret = -1;
+            }
+        }
+        else
+        {
+            ret = 1;
+        }
+    }
+
+    return ret;
+}
+
+*/
+
+template<typename T> __host__ int curarray<T>::resizearray(int arraynum, int arraysize)
+{
+    int ret = 0;
+    T* newptr = NULL;
+    int I;
+    T def;
+
+    if(arraynum>=0 && arraynum<Narrays)
+    {
+        if(N[arraynum]!=arraysize)
+        {
+            if(arraysize<=0)
+            {
+                delete[] hostarrayptrs[arraynum];
+                hostarrayptrs[arraynum] = NULL;
+                N[arraynum] = 0;
+                ret = 1;
+                return ret;
+            }
+
+            newptr = new(std::nothrow) T[arraysize];
+            if(newptr!=NULL)
+            {
+                //do I want to assume there is a copy operator? (operator=)
+                //for now, yes - write a more restrictive class later if I don't want it
+                
+                if(hostarrayptrs[arraynum]!=NULL)
+                {
+                    for(I=0;I<N[arraynum]&&I<arraysize;I++)
+                    {
+                        newptr[I] = hostarrayptrs[arraynum][I];
+                    }
+                }
+                for(I=N[arraynum];I<arraysize;I++)
+                {
+                    newptr[I] = def;
+                }
+
+                //cudaFree(hostarrayptrs[arraynum]);
+                delete[] hostarrayptrs[arraynum];
+                hostarrayptrs[arraynum] = newptr;
+                N[arraynum] = arraysize;
+                ret = 1;
+            }
+            else
+            {
+                ret = -1;
+            }
+        }
+        else
+        {
+            ret = 1;
+        }
+    }
+
+    return ret;
+}
+
+};
+
+#endif

include/amsculib3/math/amscu_comp128.hpp

@@ -0,0 +1,89 @@
+#ifndef __AMSCU_COMP128_HPP__
+#define __AMSCU_COMP128_HPP__
+
+namespace amscuda
+{
+namespace cmp
+{
+
+    class cucomp128
+    {
+        public:
+        double real;
+        double imag;
+
+        __host__ __device__ cucomp128();
+        __host__ __device__ ~cucomp128();
+        __host__ __device__ cucomp128(const cucomp128 &other);
+        __host__ __device__ cucomp128(const double &other);
+
+        __host__ __device__ cucomp128& operator=(cucomp128& other);
+        __host__ __device__ const cucomp128& operator=(const cucomp128& other);
+        __host__ __device__ cucomp128& operator=(double& other);
+        __host__ __device__ const cucomp128& operator=(const double& other);
+        
+        __host__ __device__ double& operator[](int& ind);
+        __host__ __device__ const double& operator[](const int& ind) const;
+
+        __host__ __device__ cucomp128 operator+(const cucomp128& z);
+        __host__ __device__ cucomp128 operator-(const cucomp128& z);
+        __host__ __device__ cucomp128 operator*(const cucomp128& z);
+        __host__ __device__ cucomp128 operator/(const cucomp128& z);
+
+        __host__ __device__ cucomp128 operator+(const double& z);
+        __host__ __device__ cucomp128 operator-(const double& z);
+        __host__ __device__ cucomp128 operator*(const double& z);
+        __host__ __device__ cucomp128 operator/(const double& z);
+
+        __host__ __device__ friend cucomp128 operator-(const cucomp128& z); //negation sign
+
+        //comparison operators
+        __host__ __device__ bool operator==(const cucomp128& z) const;
+        __host__ __device__ bool operator!=(const cucomp128& z) const;
+        __host__ __device__ bool operator>(const cucomp128& z) const;
+        __host__ __device__ bool operator<(const cucomp128& z) const;
+        __host__ __device__ bool operator>=(const cucomp128& z) const;
+        __host__ __device__ bool operator<=(const cucomp128& z) const;
+
+        __host__ __device__ bool isnan() const;
+        __host__ __device__ bool isinf() const;
+
+        __host__ __device__ bool isreal() const;
+        __host__ __device__ bool isimag() const;
+        __host__ __device__ bool iszero() const;
+        __host__ __device__ double arg() const;
+        __host__ __device__ double mag() const;
+        __host__ __device__ cucomp128 conj() const;
+    };
+
+    __host__ __device__ double arg(cucomp128 z);
+
+    __host__ __device__ cucomp128 dtocomp(double _r, double _i);
+    __host__ __device__ double real(cucomp128 z);
+    __host__ __device__ double imag(cucomp128 z);
+    __host__ __device__ cucomp128 sin(cucomp128 z);
+    __host__ __device__ cucomp128 cos(cucomp128 z);
+    __host__ __device__ cucomp128 tan(cucomp128 z);
+    __host__ __device__ cucomp128 exp(cucomp128 z);
+    __host__ __device__ cucomp128 log(cucomp128 z);
+    __host__ __device__ double abs(cucomp128 z);
+    __host__ __device__ cucomp128 conj(cucomp128 z);
+
+    // //need hyperbolic trig Functions
+    __host__ __device__ cucomp128 cosh(cucomp128 z);
+    __host__ __device__ cucomp128 sinh(cucomp128 z);
+    __host__ __device__ cucomp128 tanh(cucomp128 z);
+
+    __host__ __device__ cucomp128 pow(cucomp128 z1, cucomp128 z2);
+
+    // //returns "complex sign" of complex number - 0, or a unit number with same argument
+    __host__ __device__ cucomp128 csgn(cucomp128 z);
+
+void test_cucomp128_1();
+
+
+}; //end namespace cmp
+}; //end namespace amscuda
+
+#endif
+

include/amsculib3/math/amscu_comp64.hpp

@@ -0,0 +1,88 @@
+#ifndef __AMSCU_COMP64_HPP__
+#define __AMSCU_COMP64_HPP__
+
+namespace amscuda
+{
+namespace cmp
+{
+
+    class cucomp64
+    {
+        public:
+        float real;
+        float imag;
+
+        __host__ __device__ cucomp64();
+        __host__ __device__ ~cucomp64();
+        __host__ __device__ cucomp64(const cucomp64 &other);
+        __host__ __device__ cucomp64(const float &other);
+
+        __host__ __device__ cucomp64& operator=(cucomp64& other);
+        __host__ __device__ const cucomp64& operator=(const cucomp64& other);
+        __host__ __device__ cucomp64& operator=(float& other);
+        __host__ __device__ const cucomp64& operator=(const float& other);
+        
+        __host__ __device__ float& operator[](int& ind);
+        __host__ __device__ const float& operator[](const int& ind) const;
+
+        __host__ __device__ cucomp64 operator+(const cucomp64& z);
+        __host__ __device__ cucomp64 operator-(const cucomp64& z);
+        __host__ __device__ cucomp64 operator*(const cucomp64& z);
+        __host__ __device__ cucomp64 operator/(const cucomp64& z);
+
+        __host__ __device__ cucomp64 operator+(const float& z);
+        __host__ __device__ cucomp64 operator-(const float& z);
+        __host__ __device__ cucomp64 operator*(const float& z);
+        __host__ __device__ cucomp64 operator/(const float& z);
+
+        __host__ __device__ friend cucomp64 operator-(const cucomp64& z); //negation sign
+
+        //comparison operators
+        __host__ __device__ bool operator==(const cucomp64& z) const;
+        __host__ __device__ bool operator!=(const cucomp64& z) const;
+        __host__ __device__ bool operator>(const cucomp64& z) const;
+        __host__ __device__ bool operator<(const cucomp64& z) const;
+        __host__ __device__ bool operator>=(const cucomp64& z) const;
+        __host__ __device__ bool operator<=(const cucomp64& z) const;
+
+        __host__ __device__ bool isnan() const;
+        __host__ __device__ bool isinf() const;
+
+        __host__ __device__ bool isreal() const;
+        __host__ __device__ bool isimag() const;
+        __host__ __device__ bool iszero() const;
+        __host__ __device__ float arg() const;
+        __host__ __device__ float mag() const;
+        __host__ __device__ cucomp64 conj() const;
+    };
+
+    __host__ __device__ float arg(cucomp64 z);
+
+    __host__ __device__ cucomp64 dtocomp64(float _r, float _i);
+    __host__ __device__ float real(cucomp64 z);
+    __host__ __device__ float imag(cucomp64 z);
+    __host__ __device__ cucomp64 sin(cucomp64 z);
+    __host__ __device__ cucomp64 cos(cucomp64 z);
+    __host__ __device__ cucomp64 tan(cucomp64 z);
+    __host__ __device__ cucomp64 exp(cucomp64 z);
+    __host__ __device__ cucomp64 log(cucomp64 z);
+    __host__ __device__ float abs(cucomp64 z);
+    __host__ __device__ cucomp64 conj(cucomp64 z);
+
+    // //need hyperbolic trig Functions
+    __host__ __device__ cucomp64 cosh(cucomp64 z);
+    __host__ __device__ cucomp64 sinh(cucomp64 z);
+    __host__ __device__ cucomp64 tanh(cucomp64 z);
+
+    __host__ __device__ cucomp64 pow(cucomp64 z1, cucomp64 z2);
+
+    // //returns "complex sign" of complex number - 0, or a unit number with same argument
+    __host__ __device__ cucomp64 csgn(cucomp64 z);
+
+void test_cucomp64_1();
+
+}; //end namespace cmp
+}; //end namespace amscuda
+
+#endif
+

include/amsculib3/math/amscumath.hpp

@@ -0,0 +1,70 @@
+#ifndef __AMSCUMATH_HPP__
+#define __AMSCUMATH_HPP__
+
+namespace amscuda
+{
+
+
+
+    //Problem: These are not in the namespace
+    //#define nan NAN
+    //#define fnan (float) NAN
+    //#define inf INFINITY
+    //#define finf (float) INFINITY
+    //#define pi 3.1415926535897936
+
+    //These need to be the same symbol for both host and device code
+    AMSCU_CONST static const double nan = NAN;
+    AMSCU_CONST static const float fnan = (float) NAN;
+    AMSCU_CONST static const double inf = INFINITY;
+    AMSCU_CONST static const float finf = (float) INFINITY;
+    AMSCU_CONST static const double pi = 3.1415926535897936;
+    AMSCU_CONST static const float pif = 3.1415926535897936;
+
+    __host__ __device__ double dabs(const double &x);
+    __host__ __device__ float fabs(const float &x);
+    
+    template<typename T> __host__ __device__ T abs(const T in)
+    {
+        T ret = in;
+        if(in<0) ret = -in;
+        return ret;
+    }
+
+    __host__ __device__ double mod(const double &a, const double &md);
+    __host__ __device__ float mod(const float &a, const float &md);
+    __host__ __device__ int mod(const int &x, const int &n);
+    __host__ __device__ long mod(const long &x, const long &n);
+
+    __host__ __device__ int truediv(const int &x, const int &y);
+    __host__ __device__ long truediv(const long &x, const long &y);
+
+    template<typename T> __host__ __device__ T min(const T &a, const T &b);
+    template<typename T> __host__ __device__ T max(const T &a, const T &b);
+
+    __device__ __host__ double arg(const double &x, const double &y);
+    __device__ __host__ void get_azel(const double &x, const double &y, const double &z, double *az, double *el);
+    __device__ __host__ void get_azel(const double &x, const double &y, const double &z, double &az_out, double &el_out);
+
+    __device__ __host__ double argf(const float &x, const float &y);
+    __device__ __host__ void get_azelf(const float &x, const float &y, const float &z, float &az_out, float &el_out);
+
+    void test_amscumath1();
+
+
+}; //end namespace amscuda
+
+#include <amsculib3/math/amscumath_impl.hpp>
+
+//component headers
+#include <amsculib3/math/amscu_comp64.hpp>
+#include <amsculib3/math/amscu_comp128.hpp>
+#include <amsculib3/math/cuvect2.hpp>
+#include <amsculib3/math/cuvect3.hpp>
+#include <amsculib3/math/cuvect4.hpp>
+#include <amsculib3/math/cuvect2f.hpp>
+#include <amsculib3/math/cuvect3f.hpp>
+#include <amsculib3/math/cuvect4f.hpp>
+
+#endif
+

include/amsculib3/math/amscumath_impl.hpp

@@ -0,0 +1,42 @@
+#ifndef __AMSCUMATH_IMPL_HPP__
+#define __AMSCUMATH_IMPL_HPP__
+
+namespace amscuda
+{
+
+template<typename T> __host__ __device__ T min(const T &a, const T &b)
+{
+    if(a>b)
+    {
+        return b;
+    }
+    else
+    {
+        return a;
+    }
+    return a;
+}
+
+template<typename T> __host__ __device__ T max(const T &a, const T &b)
+{
+    if(a>b)
+    {
+        return a;
+    }
+    else
+    {
+        return b;
+    }
+    return a;
+}
+
+template<> __host__ __device__ double min(const double &a, const double &b);
+template<> __host__ __device__ float min(const float &a, const float &b);
+template<> __host__ __device__ double max(const double &a, const double &b);
+template<> __host__ __device__ float max(const float &a, const float &b);
+
+
+}; //end namespace amscuda
+
+#endif
+

include/amsculib3/math/cuvect2.hpp

@@ -0,0 +1,87 @@
+#ifndef __CUVECT2_HPP__
+#define __CUVECT2_HPP__
+
+namespace amscuda
+{
+
+    class cuvect2
+    {
+        public:
+        double x;
+        double y;
+
+
+        __host__ __device__ cuvect2();
+        __host__ __device__ ~cuvect2();
+        __host__ __device__ cuvect2(double _x, double _y);
+
+        __host__ __device__ double& operator[](const int I);
+        __host__ __device__ const double& operator[](const int I) const;
+
+        __host__ __device__ cuvect2 operator+(cuvect2 lhs);
+        __host__ __device__ cuvect2 operator-(cuvect2 lhs);
+        __host__ __device__ cuvect2 operator*(double lhs);
+        __host__ __device__ cuvect2 operator/(double lhs);
+        __host__ __device__ friend cuvect2 operator-(cuvect2 rhs);
+    };
+
+    class cumat2
+    {
+        public:
+        double dat[4];
+
+        __host__ __device__ cumat2();
+        __host__ __device__ ~cumat2();
+        __host__ __device__ double& operator[](const int I);
+        __host__ __device__ double& operator()(const int I, const int J);
+        __host__ __device__ double& at(const int I, const int J);
+
+        __host__ __device__ cumat2 operator+(cumat2 lhs);
+        __host__ __device__ cumat2 operator-(cumat2 lhs);
+        __host__ __device__ cumat2 operator*(double lhs);
+        __host__ __device__ cumat2 operator/(double lhs);
+        __host__ __device__ cuvect2 operator*(cuvect2 lhs);
+        __host__ __device__ cumat2 operator*(cumat2 lhs);
+        __host__ __device__ friend cumat2 operator-(cumat2 rhs);
+
+        __host__ __device__ double det();
+        __host__ __device__ cumat2 transpose();
+        __host__ __device__ cumat2 inverse();
+    };
+
+    __host__ __device__ double cuvect2_dot(cuvect2 a, cuvect2 b);
+    __host__ __device__ double cuvect2_cross(cuvect2 a, cuvect2 b);
+    __host__ __device__ double cuvect2_norm(cuvect2 a);
+    __host__ __device__ cuvect2 cuvect2_normalize(cuvect2 a);
+    __host__ __device__ cuvect2 cuvect2_proj(cuvect2 a, cuvect2 b);
+
+    //2x2 matrix operations
+    //matrix order is assumed to be mat[I,J] = mat[I+3*J]
+
+    //transpose a 2x2 matrix in place
+    __host__ __device__ void mat2_transpose(double *mat2inout);
+    
+    //copies src to dest
+    __host__ __device__ void mat2_copy(double *mat2_dest, const double *mat2_src);
+    
+    //inverts mat?inout[4]
+    __host__ __device__ void mat2_inverse(double *mat2inout);
+    
+    //rotatin matrix from angle
+    __host__ __device__ void mat2_rot_from_angle(double angle, double *mat2);
+
+    //multiplies c = a*b
+    __host__ __device__ void mat2_mult(double *mat2a, double *mat2b, double *mat2c);
+
+    // ret = a*b
+    __host__ __device__ cuvect2 mat2_mult(double *mat2a, cuvect2 b);
+
+
+
+    void test_cuvect2_1();
+    
+
+}; //end namespace amscuda
+
+#endif
+

include/amsculib3/math/cuvect2f.hpp

@@ -0,0 +1,128 @@
+#ifndef __CUVECT2F_HPP__
+#define __CUVECT2F_HPP__
+
+namespace amscuda
+{
+
+    class cuvect2f
+    {
+        public:
+        float x;
+        float y;
+
+        __host__ __device__ cuvect2f();
+        __host__ __device__ ~cuvect2f();
+        __host__ __device__ cuvect2f(const float &_x, const float &_y);
+        
+        
+        __host__ __device__ float& operator[](const int &I);
+        __host__ __device__ const float& operator[](const int &I) const;
+
+        __host__ __device__ cuvect2f operator+(const cuvect2f &rhs);
+        __host__ __device__ cuvect2f operator-(const cuvect2f &rhs);
+        __host__ __device__ cuvect2f operator*(const float &rhs);
+        __host__ __device__ cuvect2f operator/(const float &rhs);
+        __host__ __device__ friend cuvect2f operator-(const cuvect2f &rhs);
+
+        __host__ __device__ cuvect2f& operator+=(const cuvect2f &rhs);
+        __host__ __device__ cuvect2f& operator-=(const cuvect2f &rhs);
+        __host__ __device__ cuvect2f& operator/=(const float &rhs);
+        __host__ __device__ cuvect2f& operator*=(const float &rhs);
+
+    };
+
+    class cumat2f
+    {
+        public:
+        float m00,m10; //named references to force register use?
+        float m01,m11; //switched to column-major-order to match GLSL/lapack
+
+        __host__ __device__ cumat2f();
+        __host__ __device__ ~cumat2f();
+
+        __host__ __device__ cumat2f(
+            const float & _m00, const float & _m01,
+            const float & _m10, const float & _m11
+        );
+
+        __host__ __device__ explicit cumat2f(const float* data2x2);
+
+        __host__ __device__ float& operator[](const int &I);
+        __host__ __device__ float& operator()(const int &I, const int &J);
+        __host__ __device__ float& at(const int &I, const int &J);
+
+        __host__ __device__ const float& operator[](const int &I) const;
+        __host__ __device__ const float& operator()(const int &I, const int &J) const;
+        __host__ __device__ const float& at(const int &I, const int &J) const;
+
+        __host__ __device__ cumat2f operator+(const cumat2f &rhs);
+        __host__ __device__ cumat2f operator-(const cumat2f &rhs);
+        __host__ __device__ cumat2f operator*(const float &rhs);
+        __host__ __device__ cumat2f operator/(const float &rhs);
+        __host__ __device__ cuvect2f operator*(const cuvect2f &rhs);
+        __host__ __device__ cumat2f operator*(const cumat2f &rhs);
+        __host__ __device__ friend cumat2f operator-(const cumat2f &rhs);
+        
+        __host__ __device__ float det();
+        __host__ __device__ cumat2f transpose();
+        __host__ __device__ cumat2f inverse();
+
+        __host__ __device__ float* data(); //pointer to float[4] representation of matrix
+        __host__ __device__ const float* data() const; //pointer to float[4] representation of matrix
+
+        //In place operations (to save GPU register use)
+        __host__ __device__ cumat2f& operator+=(const cumat2f &rhs);
+        __host__ __device__ cumat2f& operator-=(const cumat2f &rhs);
+        __host__ __device__ cumat2f& operator/=(const float &rhs);
+        __host__ __device__ cumat2f& operator*=(const float &rhs);
+        __host__ __device__ cumat2f& operator*=(const cumat2f &rhs);
+    };
+
+    __host__ __device__ float cuvect2f_dot(const cuvect2f &a, const cuvect2f &b);
+    __host__ __device__ float cuvect2f_cross(const cuvect2f &a, const cuvect2f &b);
+    __host__ __device__ float cuvect2f_norm(const cuvect2f &a);
+    __host__ __device__ cuvect2f cuvect2f_normalize(const cuvect2f &a);
+    __host__ __device__ cuvect2f cuvect2f_proj(const cuvect2f &a, const cuvect2f &b);
+
+    __host__ __device__ cumat2f cumat2f_rot_from_angle(const float &angle);
+
+
+    ///////////
+    // Tests //
+    ///////////
+    void test_cuvect2f_1();
+
+
+    ///////////////////////////
+    //legacy array operations//
+    ///////////////////////////
+
+    //2x2 matrix operations
+    //matrix order is assumed to be mat[I,J] = mat[I+2*J]
+
+    //transpose a 2x2 matrix in place
+    __host__ __device__ void mat2f_transpose(float *mat2inout);
+    
+    //copies src to dest
+    __host__ __device__ void mat2f_copy(float *mat2f_dest, const float *mat2f_src);
+    
+    //inverts mat?inout[4]
+    __host__ __device__ void mat2f_inverse(float *mat2inout);
+    
+    //rotatin matrix from angle
+    __host__ __device__ void mat2f_rot_from_angle(float angle, float *mat2);
+
+    //multiplies c = a*b
+    __host__ __device__ void mat2f_mult(float *mat2a, float *mat2b, float *mat2c);
+
+    // ret = a*b
+    __host__ __device__ cuvect2f mat2f_mult(float *mat2a, const cuvect2f &b);
+
+
+
+
+
+};
+
+#endif
+

include/amsculib3/math/cuvect3.hpp

@@ -0,0 +1,122 @@
+#ifndef __CUVECT3_HPP__
+#define __CUVECT3_HPP__
+
+namespace amscuda
+{
+
+    class cuvect3
+    {
+        public:
+        double x;
+        double y;
+        double z;
+
+        __host__ __device__ cuvect3();
+        __host__ __device__ ~cuvect3();
+        __host__ __device__ cuvect3(const double &_x, const double &_y, const double &_z);
+        
+        
+        __host__ __device__ double& operator[](const int &I);
+        __host__ __device__ const double& operator[](const int &I) const;
+
+        __host__ __device__ cuvect3 operator+(const cuvect3 &rhs);
+        __host__ __device__ cuvect3 operator-(const cuvect3 &rhs);
+        __host__ __device__ cuvect3 operator*(const double &rhs);
+        __host__ __device__ cuvect3 operator/(const double &rhs);
+        __host__ __device__ friend cuvect3 operator-(const cuvect3 &rhs);
+
+        __host__ __device__ cuvect3& operator+=(const cuvect3 &rhs);
+        __host__ __device__ cuvect3& operator-=(const cuvect3 &rhs);
+        __host__ __device__ cuvect3& operator/=(const double &rhs);
+        __host__ __device__ cuvect3& operator*=(const double &rhs);
+
+    };
+
+    class cumat3
+    {
+        public:
+        double m00,m10,m20; //named references to force register use?
+        double m01,m11,m21; //switched to column-major-order to match GLSL/lapack
+        double m02,m12,m22;
+
+        __host__ __device__ cumat3();
+        __host__ __device__ ~cumat3();
+        __host__ __device__ cumat3(
+            const double & _m00, const double & _m01, const double & _m02,
+            const double & _m10, const double & _m11, const double & _m12,
+            const double & _m20, const double & _m21, const double & _m22
+        );
+
+        __host__ __device__ double& operator[](const int &I);
+        __host__ __device__ double& operator()(const int &I, const int &J);
+        __host__ __device__ double& at(const int &I, const int &J);
+
+        __host__ __device__ const double& operator[](const int &I) const;
+        __host__ __device__ const double& operator()(const int &I, const int &J) const;
+        __host__ __device__ const double& at(const int &I, const int &J) const;
+
+        __host__ __device__ cumat3 operator+(const cumat3 &rhs);
+        __host__ __device__ cumat3 operator-(const cumat3 &rhs);
+        __host__ __device__ cumat3 operator*(const double &rhs);
+        __host__ __device__ cumat3 operator/(const double &rhs);
+        __host__ __device__ cuvect3 operator*(const cuvect3 &rhs);
+        __host__ __device__ cumat3 operator*(const cumat3 &rhs);
+        __host__ __device__ friend cumat3 operator-(const cumat3 &rhs);
+        
+        __host__ __device__ double det();
+        __host__ __device__ cumat3 transpose();
+        __host__ __device__ cumat3 inverse();
+
+        __host__ __device__ double* data(); //pointer to double[9] representation of matrix
+        __host__ __device__ const double* data() const; //pointer to double[9] representation of matrix
+
+        //In place operations (to save GPU register use)
+        __host__ __device__ cumat3& operator+=(const cumat3 &rhs);
+        __host__ __device__ cumat3& operator-=(const cumat3 &rhs);
+        __host__ __device__ cumat3& operator/=(const double &rhs);
+        __host__ __device__ cumat3& operator*=(const double &rhs);
+        __host__ __device__ cumat3& operator*=(const cumat3 &rhs);
+        
+    };
+
+    __host__ __device__ double cuvect3_dot(const cuvect3 &a,const cuvect3 &b);
+    __host__ __device__ cuvect3 cuvect3_cross(const cuvect3 &a, const cuvect3 &b);
+    __host__ __device__ double cuvect3_norm(const cuvect3 &a);
+    __host__ __device__ cuvect3 cuvect3_normalize(const cuvect3 &a);
+    __host__ __device__ cuvect3 cuvect3_proj(const cuvect3 &a, const cuvect3 &b);
+
+    __host__ __device__ cumat3 hodge_dual(const cuvect3 &vin);
+    __host__ __device__ cuvect3 hodge_dual(const cumat3 &min);
+    __host__ __device__ cumat3 rotmat_from_axisangle(const cuvect3 &axis, const double &angle);
+
+
+    //3x3 matrix operations
+    //matrix order is assumed to be mat[I,J] = mat[I+3*J]
+    
+    //transposes a 3x3 (9 element) matrix
+    __host__ __device__ void mat3f_transpose(double *mat3inout);
+
+    //copies src to dest
+    __host__ __device__ void mat3f_copy(double *mat3f_dest, const double *mat3f_src);
+
+    //returns determinant of 3x3 matrix
+    __host__ __device__ double mat3f_det(double *mat3in);
+
+    //inverts a 3x3 (9 element) matrix
+    __host__ __device__ void mat3f_inverse(double *mat3inout);
+
+    __host__ __device__ cuvect3 mat3f_mult(double *mat3in, const cuvect3 &cvin);
+    __host__ __device__ void mat3f_mult(double *matina, double *matinb, double *matout);
+
+    __host__ __device__ void mat3f_hodgedual(const cuvect3 &vecin, double *matout);
+    __host__ __device__ void mat3f_hodgedual(double *matin, cuvect3 &vecout);
+
+    //returns direction cosine rotation matrix from axis and angle
+    __host__ __device__ void mat3f_rot_from_axisangle(cuvect3 axis, double angle, double *matout);
+
+    __host__ void test_cudavect_logic1();
+
+}; //end namespace amscuda
+
+#endif
+

include/amsculib3/math/cuvect3f.hpp

@@ -0,0 +1,102 @@
+#ifndef __CUVECT3F_HPP__
+#define __CUVECT3F_HPP__
+
+namespace amscuda
+{
+
+    class cuvect3f
+    {
+        public:
+        float x;
+        float y;
+        float z;
+
+        __host__ __device__ cuvect3f();
+        __host__ __device__ ~cuvect3f();
+        __host__ __device__ cuvect3f(const float &_x, const float &_y, const float &_z);
+        
+        
+        __host__ __device__ float& operator[](const int &I);
+        __host__ __device__ const float& operator[](const int &I) const;
+
+        __host__ __device__ cuvect3f operator+(const cuvect3f &rhs);
+        __host__ __device__ cuvect3f operator-(const cuvect3f &rhs);
+        __host__ __device__ cuvect3f operator*(const float &rhs);
+        __host__ __device__ cuvect3f operator/(const float &rhs);
+        __host__ __device__ friend cuvect3f operator-(const cuvect3f &rhs);
+
+        __host__ __device__ cuvect3f& operator+=(const cuvect3f &rhs);
+        __host__ __device__ cuvect3f& operator-=(const cuvect3f &rhs);
+        __host__ __device__ cuvect3f& operator/=(const float &rhs);
+        __host__ __device__ cuvect3f& operator*=(const float &rhs);
+
+    };
+
+    class cumat3f
+    {
+        public:
+        float m00,m10,m20; //named references to force register use?
+        float m01,m11,m21; //switched to column-major-order to match GLSL/lapack
+        float m02,m12,m22;
+
+        __host__ __device__ cumat3f();
+        __host__ __device__ ~cumat3f();
+        __host__ __device__ cumat3f(
+            const float & _m00, const float & _m01, const float & _m02,
+            const float & _m10, const float & _m11, const float & _m12,
+            const float & _m20, const float & _m21, const float & _m22
+        );
+
+        __host__ __device__ explicit cumat3f(const float *data9);
+        
+
+        __host__ __device__ float& operator[](const int &I);
+        __host__ __device__ float& operator()(const int &I, const int &J);
+        __host__ __device__ float& at(const int &I, const int &J);
+
+        __host__ __device__ const float& operator[](const int &I) const;
+        __host__ __device__ const float& operator()(const int &I, const int &J) const;
+        __host__ __device__ const float& at(const int &I, const int &J) const;
+
+        __host__ __device__ cumat3f operator+(const cumat3f &rhs);
+        __host__ __device__ cumat3f operator-(const cumat3f &rhs);
+        __host__ __device__ cumat3f operator*(const float &rhs);
+        __host__ __device__ cumat3f operator/(const float &rhs);
+        __host__ __device__ cuvect3f operator*(const cuvect3f &rhs);
+        __host__ __device__ cumat3f operator*(const cumat3f &rhs);
+        __host__ __device__ friend cumat3f operator-(const cumat3f &rhs);
+        
+        __host__ __device__ float det();
+        __host__ __device__ cumat3f transpose();
+        __host__ __device__ cumat3f inverse();
+
+        __host__ __device__ float* data(); //pointer to float[9] representation of matrix
+        __host__ __device__ const float* data() const; //pointer to float[9] representation of matrix
+
+        //In place operations (to save GPU register use)
+        __host__ __device__ cumat3f& operator+=(const cumat3f &rhs);
+        __host__ __device__ cumat3f& operator-=(const cumat3f &rhs);
+        __host__ __device__ cumat3f& operator/=(const float &rhs);
+        __host__ __device__ cumat3f& operator*=(const float &rhs);
+        __host__ __device__ cumat3f& operator*=(const cumat3f &rhs);
+        
+    };
+
+    __host__ __device__ float cuvect3f_dot(const cuvect3f &a,const cuvect3f &b);
+    __host__ __device__ cuvect3f cuvect3f_cross(const cuvect3f &a, const cuvect3f &b);
+    __host__ __device__ float cuvect3f_norm(const cuvect3f &a);
+    __host__ __device__ cuvect3f cuvect3f_normalize(const cuvect3f &a);
+    __host__ __device__ cuvect3f cuvect3f_proj(const cuvect3f &a, const cuvect3f &b);
+
+    __host__ __device__ cumat3f hodge_dual(const cuvect3f &vin);
+    __host__ __device__ cuvect3f hodge_dual(const cumat3f &min);
+    __host__ __device__ cumat3f rotmat_from_axisangle(const cuvect3f &axis, const float &angle);
+
+    __host__ void test_cudavectf_logic1();
+
+};
+
+#include <amsculib2/cuvect3f_legacy.hpp>
+
+#endif
+

include/amsculib3/math/cuvect4.hpp

@@ -0,0 +1,61 @@
+#ifndef __CUVECT4_HPP__
+#define __CUVECT4_HPP__
+
+namespace amscuda
+{
+
+    class cuvect4
+    {
+        public:
+        double x;
+        double y;
+        double z;
+        double w;
+
+        __host__ __device__ cuvect4();
+        __host__ __device__ ~cuvect4();
+        __host__ __device__ cuvect4(double _x, double _y, double _z, double _w);
+                
+        __host__ __device__ double& operator[](const int I);
+        __host__ __device__ const double& operator[](const int I) const;
+
+        __host__ __device__ cuvect4 operator+(cuvect4 lhs);
+        __host__ __device__ cuvect4 operator-(cuvect4 lhs);
+        __host__ __device__ cuvect4 operator*(double lhs);
+        __host__ __device__ cuvect4 operator/(double lhs);
+        __host__ __device__ friend cuvect4 operator-(cuvect4 rhs);
+    };
+
+    class cumat4
+    {
+        public:
+        double dat[16];
+
+        __host__ __device__ cumat4();
+        __host__ __device__ ~cumat4();
+        __host__ __device__ double& operator[](const int I);
+        __host__ __device__ double& operator()(const int I, const int J);
+        __host__ __device__ double& at(const int I, const int J);
+
+        __host__ __device__ cumat4 operator+(cumat4 lhs);
+        __host__ __device__ cumat4 operator-(cumat4 lhs);
+        __host__ __device__ cumat4 operator*(double lhs);
+        __host__ __device__ cumat4 operator/(double lhs);
+        __host__ __device__ cuvect4 operator*(cuvect4 lhs);
+        __host__ __device__ cumat4 operator*(cumat4 lhs);
+        __host__ __device__ friend cumat4 operator-(cumat4 rhs);
+
+        __host__ __device__ double det();
+        __host__ __device__ cumat4 transpose();
+        __host__ __device__ cumat4 inverse();
+    };
+
+    __host__ __device__ double cuvect4_dot(cuvect4 a, cuvect4 b);
+    __host__ __device__ double cuvect4_norm(cuvect4 a);
+    __host__ __device__ cuvect4 cuvect4_normalize(cuvect4 a);
+    __host__ __device__ cuvect4 cuvect4_proj(cuvect4 a, cuvect4 b);
+
+}; //end namespace amscuda
+
+#endif
+

include/amsculib3/math/cuvect4f.hpp

@@ -0,0 +1,97 @@
+#ifndef __CUVECT4F_HPP__
+#define __CUVECT4F_HPP__
+
+namespace amscuda
+{
+
+    class cuvect4f
+    {
+        public:
+        float x;
+        float y;
+        float z;
+        float w;
+
+        __host__ __device__ cuvect4f();
+        __host__ __device__ ~cuvect4f();
+        __host__ __device__ cuvect4f(const float &_x, const float &_y, const float &_z, const float &_w);
+        
+        
+        __host__ __device__ float& operator[](const int &I);
+        __host__ __device__ const float& operator[](const int &I) const;
+
+        __host__ __device__ cuvect4f operator+(const cuvect4f &rhs);
+        __host__ __device__ cuvect4f operator-(const cuvect4f &rhs);
+        __host__ __device__ cuvect4f operator*(const float &rhs);
+        __host__ __device__ cuvect4f operator/(const float &rhs);
+        __host__ __device__ friend cuvect4f operator-(const cuvect4f &rhs);
+
+        __host__ __device__ cuvect4f& operator+=(const cuvect4f &rhs);
+        __host__ __device__ cuvect4f& operator-=(const cuvect4f &rhs);
+        __host__ __device__ cuvect4f& operator/=(const float &rhs);
+        __host__ __device__ cuvect4f& operator*=(const float &rhs);
+    };
+
+    class cumat4f
+    {
+        public:
+        //float dat[16];
+
+        float m00,m10,m20,m30; //named references to force register use?
+        float m01,m11,m21,m31; //switched to column-major-order to match GLSL/lapack
+        float m02,m12,m22,m32;
+        float m03,m13,m23,m33;
+
+        __host__ __device__ cumat4f();
+        __host__ __device__ ~cumat4f();
+        __host__ __device__ cumat4f(
+            const float & _m00, const float & _m01, const float & _m02, const float & _m03,
+            const float & _m10, const float & _m11, const float & _m12, const float & _m13,
+            const float & _m20, const float & _m21, const float & _m22, const float & _m23,
+            const float & _m30, const float & _m31, const float & _m32, const float & _m33
+        );
+
+        __host__ __device__ explicit cumat4f(const float *data16);
+        
+        //__forceinline__
+        __host__ __device__ float& operator[](const int &I);
+        __host__ __device__ float& operator()(const int &I, const int &J);
+        __host__ __device__ float& at(const int &I, const int &J);
+
+        __host__ __device__ const float& operator[](const int &I) const;
+        __host__ __device__ const float& operator()(const int &I, const int &J) const;
+        __host__ __device__ const float& at(const int &I, const int &J) const;
+
+        __host__ __device__ cumat4f operator+(const cumat4f &rhs);
+        __host__ __device__ cumat4f operator-(const cumat4f &rhs);
+        __host__ __device__ cumat4f operator*(const float &rhs);
+        __host__ __device__ cumat4f operator/(const float &rhs);
+        __host__ __device__ cuvect4f operator*(const cuvect4f &rhs);
+        __host__ __device__ cumat4f operator*(const cumat4f &rhs);
+        __host__ __device__ friend cumat4f operator-(const cumat4f &rhs);
+        
+        __host__ __device__ float det();
+        __host__ __device__ cumat4f transpose();
+        __host__ __device__ cumat4f inverse();
+
+        __host__ __device__ float* data(); //pointer to float[9] representation of matrix
+        __host__ __device__ const float* data() const; //pointer to float[9] representation of matrix
+
+        //In place operations (to save GPU register use)
+        __host__ __device__ cumat4f& operator+=(const cumat4f &rhs);
+        __host__ __device__ cumat4f& operator-=(const cumat4f &rhs);
+        __host__ __device__ cumat4f& operator/=(const float &rhs);
+        __host__ __device__ cumat4f& operator*=(const float &rhs);
+        __host__ __device__ cumat4f& operator*=(const cumat4f &rhs);
+    };
+
+    __host__ __device__ float cuvect4f_dot(cuvect4f &a, cuvect4f &b);
+    __host__ __device__ float cuvect4f_norm(cuvect4f &a);
+    __host__ __device__ cuvect4f cuvect4f_normalize(cuvect4f &a);
+    __host__ __device__ cuvect4f cuvect4f_proj(cuvect4f &a, cuvect4f &b);
+
+
+};
+
+#endif
+

make_linux.py

@@ -0,0 +1,22 @@
+#!/usr/bin/python3 
+
+import os,sys,math
+from build.amsbuildlib4 import *
+
+if(len(sys.argv)>=2):
+    if(sys.argv[1]=="clean"):
+        obj_list = flist('./build_linux64',recurse=True,exts=['.o'])
+        for o in obj_list:
+            os.remove('{}'.format(o))
+    exit()
+
+os.system('python3 ./build/make.linux64.lib.py')
+os.system('python3 ./build/make.linux64.test.py')
+
+obj_list = flist('./src',recurse=True,exts=['.o'])
+for o in obj_list:
+    os.remove('{}'.format(o))
+
+os.chdir('./build_linux64')
+callproc('./test')
+os.chdir('..')

make_mingw.py

@@ -0,0 +1,28 @@
+#!/usr/bin/python3
+#!/usr/bin/python3 
+
+import os,sys,math
+from build.amsbuildlib4 import *
+
+if(len(sys.argv)>=2):
+    if(sys.argv[1]=="clean"):
+        obj_list = flist('./build_mingw64',recurse=True,exts=['.o'])
+        for o in obj_list:
+            os.remove('{}'.format(o))
+    exit()
+
+os.system('python3 ./build/make.mingw64.lib.py')
+os.system('python3 ./build/make.mingw64.test.py')
+
+obj_list = flist('./src',recurse=True,exts=['.o','.obj'])
+for o in obj_list:
+    os.remove('{}'.format(o))
+
+if(sys.platform!="win32"):
+    os.chdir('./build_mingw64')
+    callproc('wine ./test.exe')
+    os.chdir('..')
+else:
+    os.chdir('./build_mingw64')
+    callproc('test.exe')
+    os.chdir('..')

old/9apr26_prerefactor/include/amsculib2/amscu_comp128.hpp

@@ -0,0 +1,89 @@
+#ifndef __AMSCU_COMP128_HPP__
+#define __AMSCU_COMP128_HPP__
+
+namespace amscuda
+{
+namespace cmp
+{
+
+    class cucomp128
+    {
+        public:
+        double real;
+        double imag;
+
+        __host__ __device__ cucomp128();
+        __host__ __device__ ~cucomp128();
+        __host__ __device__ cucomp128(const cucomp128 &other);
+        __host__ __device__ cucomp128(const double &other);
+
+        __host__ __device__ cucomp128& operator=(cucomp128& other);
+        __host__ __device__ const cucomp128& operator=(const cucomp128& other);
+        __host__ __device__ cucomp128& operator=(double& other);
+        __host__ __device__ const cucomp128& operator=(const double& other);
+        
+        __host__ __device__ double& operator[](int& ind);
+        __host__ __device__ const double& operator[](const int& ind) const;
+
+        __host__ __device__ cucomp128 operator+(const cucomp128& z);
+        __host__ __device__ cucomp128 operator-(const cucomp128& z);
+        __host__ __device__ cucomp128 operator*(const cucomp128& z);
+        __host__ __device__ cucomp128 operator/(const cucomp128& z);
+
+        __host__ __device__ cucomp128 operator+(const double& z);
+        __host__ __device__ cucomp128 operator-(const double& z);
+        __host__ __device__ cucomp128 operator*(const double& z);
+        __host__ __device__ cucomp128 operator/(const double& z);
+
+        __host__ __device__ friend cucomp128 operator-(const cucomp128& z); //negation sign
+
+        //comparison operators
+        __host__ __device__ bool operator==(const cucomp128& z) const;
+        __host__ __device__ bool operator!=(const cucomp128& z) const;
+        __host__ __device__ bool operator>(const cucomp128& z) const;
+        __host__ __device__ bool operator<(const cucomp128& z) const;
+        __host__ __device__ bool operator>=(const cucomp128& z) const;
+        __host__ __device__ bool operator<=(const cucomp128& z) const;
+
+        __host__ __device__ bool isnan() const;
+        __host__ __device__ bool isinf() const;
+
+        __host__ __device__ bool isreal() const;
+        __host__ __device__ bool isimag() const;
+        __host__ __device__ bool iszero() const;
+        __host__ __device__ double arg() const;
+        __host__ __device__ double mag() const;
+        __host__ __device__ cucomp128 conj() const;
+    };
+
+    __host__ __device__ double arg(cucomp128 z);
+
+    __host__ __device__ cucomp128 dtocomp(double _r, double _i);
+    __host__ __device__ double real(cucomp128 z);
+    __host__ __device__ double imag(cucomp128 z);
+    __host__ __device__ cucomp128 sin(cucomp128 z);
+    __host__ __device__ cucomp128 cos(cucomp128 z);
+    __host__ __device__ cucomp128 tan(cucomp128 z);
+    __host__ __device__ cucomp128 exp(cucomp128 z);
+    __host__ __device__ cucomp128 log(cucomp128 z);
+    __host__ __device__ double abs(cucomp128 z);
+    __host__ __device__ cucomp128 conj(cucomp128 z);
+
+    // //need hyperbolic trig Functions
+    __host__ __device__ cucomp128 cosh(cucomp128 z);
+    __host__ __device__ cucomp128 sinh(cucomp128 z);
+    __host__ __device__ cucomp128 tanh(cucomp128 z);
+
+    __host__ __device__ cucomp128 pow(cucomp128 z1, cucomp128 z2);
+
+    // //returns "complex sign" of complex number - 0, or a unit number with same argument
+    __host__ __device__ cucomp128 csgn(cucomp128 z);
+
+void test_cucomp128_1();
+
+
+}; //end namespace cmp
+}; //end namespace amscuda
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/amscu_comp64.hpp

@@ -0,0 +1,88 @@
+#ifndef __AMSCU_COMP64_HPP__
+#define __AMSCU_COMP64_HPP__
+
+namespace amscuda
+{
+namespace cmp
+{
+
+    class cucomp64
+    {
+        public:
+        float real;
+        float imag;
+
+        __host__ __device__ cucomp64();
+        __host__ __device__ ~cucomp64();
+        __host__ __device__ cucomp64(const cucomp64 &other);
+        __host__ __device__ cucomp64(const float &other);
+
+        __host__ __device__ cucomp64& operator=(cucomp64& other);
+        __host__ __device__ const cucomp64& operator=(const cucomp64& other);
+        __host__ __device__ cucomp64& operator=(float& other);
+        __host__ __device__ const cucomp64& operator=(const float& other);
+        
+        __host__ __device__ float& operator[](int& ind);
+        __host__ __device__ const float& operator[](const int& ind) const;
+
+        __host__ __device__ cucomp64 operator+(const cucomp64& z);
+        __host__ __device__ cucomp64 operator-(const cucomp64& z);
+        __host__ __device__ cucomp64 operator*(const cucomp64& z);
+        __host__ __device__ cucomp64 operator/(const cucomp64& z);
+
+        __host__ __device__ cucomp64 operator+(const float& z);
+        __host__ __device__ cucomp64 operator-(const float& z);
+        __host__ __device__ cucomp64 operator*(const float& z);
+        __host__ __device__ cucomp64 operator/(const float& z);
+
+        __host__ __device__ friend cucomp64 operator-(const cucomp64& z); //negation sign
+
+        //comparison operators
+        __host__ __device__ bool operator==(const cucomp64& z) const;
+        __host__ __device__ bool operator!=(const cucomp64& z) const;
+        __host__ __device__ bool operator>(const cucomp64& z) const;
+        __host__ __device__ bool operator<(const cucomp64& z) const;
+        __host__ __device__ bool operator>=(const cucomp64& z) const;
+        __host__ __device__ bool operator<=(const cucomp64& z) const;
+
+        __host__ __device__ bool isnan() const;
+        __host__ __device__ bool isinf() const;
+
+        __host__ __device__ bool isreal() const;
+        __host__ __device__ bool isimag() const;
+        __host__ __device__ bool iszero() const;
+        __host__ __device__ float arg() const;
+        __host__ __device__ float mag() const;
+        __host__ __device__ cucomp64 conj() const;
+    };
+
+    __host__ __device__ float arg(cucomp64 z);
+
+    __host__ __device__ cucomp64 dtocomp64(float _r, float _i);
+    __host__ __device__ float real(cucomp64 z);
+    __host__ __device__ float imag(cucomp64 z);
+    __host__ __device__ cucomp64 sin(cucomp64 z);
+    __host__ __device__ cucomp64 cos(cucomp64 z);
+    __host__ __device__ cucomp64 tan(cucomp64 z);
+    __host__ __device__ cucomp64 exp(cucomp64 z);
+    __host__ __device__ cucomp64 log(cucomp64 z);
+    __host__ __device__ float abs(cucomp64 z);
+    __host__ __device__ cucomp64 conj(cucomp64 z);
+
+    // //need hyperbolic trig Functions
+    __host__ __device__ cucomp64 cosh(cucomp64 z);
+    __host__ __device__ cucomp64 sinh(cucomp64 z);
+    __host__ __device__ cucomp64 tanh(cucomp64 z);
+
+    __host__ __device__ cucomp64 pow(cucomp64 z1, cucomp64 z2);
+
+    // //returns "complex sign" of complex number - 0, or a unit number with same argument
+    __host__ __device__ cucomp64 csgn(cucomp64 z);
+
+void test_cucomp64_1();
+
+}; //end namespace cmp
+}; //end namespace amscuda
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/amscu_cudafunctions.hpp

@@ -0,0 +1,40 @@
+#ifndef __AMSCU_CUDAFUNCTIONS_HPP__
+#define __AMSCU_CUDAFUNCTIONS_HPP__
+
+
+namespace amscuda
+{
+    // device memory operations
+    // I'm trying to avoid some of the boilerplate mental overhead involved
+    // in calling cuda functions and handling errors
+
+    //frees devbuffer if it is not already NULL, and sets devbuffer to NULL
+    //wrapper to cudaFree
+    template<typename T> int cuda_free(T **devptr);
+
+    //copies hostbuffer to devbuffer
+    //initializes devbuffer from NULL if devbuffer is NULL
+    //if overwrite is true, deletes and reallocates devbuffer on device (for resizing)
+    template<typename T> int buffer_copytodevice(T *hostbuffer, T **devbuffer, long N, bool overwrite);
+
+    //copies info from devbuffer to hostbuffer
+    //initialzies hostbuffer from NULL if NULL
+    //if overwrite is true, deletes and reallocates hostbuffer on host with new[] (for resizing)
+    template<typename T> int buffer_copyfromdevice(T *devbuffer, T **hostbuffer, long N, bool overwrite);
+
+    //wrapper for cudaMemcpy - copies an item or struct (count 1) to the device
+    //initializes devptr from NULL if not already initialized
+    template<typename T> int cuda_copytodevice(T *hostptr, T **devptr);
+
+    //wrapper for cudaMemcpy - copies an item or struct (count 1) from device
+    //initializes hostptr from NULL with new if not already initialized
+    template<typename T> int cuda_copyfromdevice(T *devptr, T **hostptr);
+    
+    int cuda_errortrap(const char *msgheader);
+
+};
+
+#include <amsculib2/amscu_cudafunctions_impl.hpp>
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/amscu_cudafunctions_impl.hpp

@@ -0,0 +1,228 @@
+#ifndef __AMSCU_CUDAFUNCTIONS_IMPL_HPP__
+#define __AMSCU_CUDAFUNCTIONS_IMPL_HPP__
+
+namespace amscuda
+{
+
+//frees devbuffer if it is not already NULL, and sets devbuffer to NULL
+//wrapper to cudaFree
+template<typename T> int cuda_free(T **devptr)
+{
+    int ret = 0;
+    cudaError_t err = cudaSuccess;
+
+    if(*devptr==NULL)
+    {
+        return ret; //devbuffer is already NULL/freed
+    }
+
+    err = cudaFree(*devptr);
+    if(err!=cudaSuccess)
+    {
+        ret = -1; //failed to free device pointer
+        *devptr = NULL; // - ? should only happen if I'm trying to double-free something
+    }
+    else
+    {
+        ret = 1;
+        *devptr = NULL;
+    }
+
+    return ret;
+}
+
+//copies hostbuffer to devbuffer
+//initializes devbuffer from NULL if devbuffer is NULL
+//if overwrite is true, deletes and reallocates devbuffer on device (for resizing)
+template<typename T> int buffer_copytodevice(T *hostbuffer, T **devbuffer, long N, bool overwrite)
+{
+    int ret = 0;
+    cudaError_t err = cudaSuccess;
+
+    if(N<=0)
+    {
+        ret = 0;
+        return ret;
+    }
+
+    if(hostbuffer==NULL)
+    {
+        ret = -2; //host buffer is NULL
+        return ret;
+    }
+
+    if(overwrite==1)
+    {
+        if(*devbuffer !=NULL)
+        {
+            cuda_free(devbuffer);
+        }
+    }
+
+    if(*devbuffer==NULL)
+    {
+        err = cudaMalloc(devbuffer,sizeof(T)*N);
+        if(err!=cudaSuccess)
+        {
+            ret = -3; //failed to allocate
+            *devbuffer = NULL;
+            return ret;
+        }
+    }
+
+    err = cudaMemcpy(*devbuffer,hostbuffer,sizeof(T)*N,cudaMemcpyHostToDevice);
+    if(err!=cudaSuccess)
+    {
+        ret = -4; //failed to copy
+    }
+    else
+    {
+        ret = 1;
+    }
+
+
+    return ret;
+}
+
+//copies info from devbuffer to hostbuffer
+//initialzies hostbuffer from NULL if NULL
+//if overwrite is true, deletes and reallocates hostbuffer on host (for resizing)
+template<typename T> int buffer_copyfromdevice(T *devbuffer, T **hostbuffer, long N, bool overwrite)
+{
+    int ret = 0;
+    cudaError_t err = cudaSuccess;
+
+    if(N<=0)
+    {
+        ret = 0;
+        return ret;
+    }
+
+    if(devbuffer==NULL)
+    {
+        ret = -5; //null dev buffer
+        return ret;
+    }
+
+    if(overwrite==1 && *hostbuffer!=NULL)
+    {
+        delete[] (*hostbuffer); hostbuffer = NULL;
+    }
+
+    if(*hostbuffer==NULL)
+    {
+        *hostbuffer = new(std::nothrow) T[N];
+        if(*hostbuffer==NULL)
+        {
+            ret = -6; //failed to allocate host buffer
+            return ret;
+        }
+    }
+
+    err = cudaMemcpy(*hostbuffer, devbuffer, sizeof(T)*N, cudaMemcpyDeviceToHost);
+    if(err!=cudaSuccess)
+    {
+        ret = -7; //failed to copy
+    }
+    else
+    {
+        ret = 1;
+    }
+    
+    return ret;
+}
+
+//wrapper for cudaMemcpy - copies an item or struct (count 1) to the device
+//initializes devptr from NULL if not already initialized
+template<typename T> int cuda_copytodevice(T *hostptr, T **devptr)
+{
+    int ret = 0;
+    cudaError_t err = cudaSuccess;
+    bool overwrite = 1;
+
+    if(hostptr==NULL)
+    {
+        ret = -2; //host buffer is NULL
+        return ret;
+    }
+
+    if(overwrite==1)
+    {
+        if(*devptr !=NULL)
+        {
+            cuda_free(devptr);
+        }
+    }
+
+    if(*devptr==NULL)
+    {
+        err = cudaMalloc(devptr,sizeof(T));
+        if(err!=cudaSuccess)
+        {
+            ret = -3; //failed to allocate
+            *devptr = NULL;
+            return ret;
+        }
+    }
+
+    err = cudaMemcpy(*devptr,hostptr,sizeof(T),cudaMemcpyHostToDevice);
+    if(err!=cudaSuccess)
+    {
+        ret = -4; //failed to copy
+    }
+    else
+    {
+        ret = 1;
+    }
+
+
+    return ret;
+}
+
+//wrapper for cudaMemcpy - copies an item or struct (count 1) from device
+//initializes hostptr from NULL with new if not already initialized
+template<typename T> int cuda_copyfromdevice(T *devptr, T **hostptr)
+{
+    int ret = 0;
+    cudaError_t err = cudaSuccess;
+    bool overwrite = 1;
+
+    if(devptr==NULL)
+    {
+        ret = -5; //null dev buffer
+        return ret;
+    }
+
+    if(overwrite==1 && *hostptr!=NULL)
+    {
+        delete (*hostptr); hostptr = NULL;
+    }
+
+    if(*hostptr==NULL)
+    {
+        *hostptr = new(std::nothrow) T;
+        if(*hostptr==NULL)
+        {
+            ret = -6; //failed to allocate host buffer
+            return ret;
+        }
+    }
+
+    err = cudaMemcpy(*hostptr, devptr, sizeof(T), cudaMemcpyDeviceToHost);
+    if(err!=cudaSuccess)
+    {
+        ret = -7; //failed to copy
+    }
+    else
+    {
+        ret = 1;
+    }
+    
+    return ret;
+}
+
+
+};
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/amscu_random.hpp

@@ -0,0 +1,55 @@
+#ifndef __AMSCU_RANDOM_HPP__
+#define __AMSCU_RANDOM_HPP__
+
+namespace amscuda
+{
+
+// Random Number Gerneators
+
+
+// faster floating point hash function used in fractal generators
+__device__ __host__ float fhash1d_su(float x);
+
+__device__ __host__ float fhash3d_su(float x, float y, float z);
+
+__device__ __host__ float fhash4d_su(float x, float y, float z, float w);
+
+
+//////////////////////////////////////////////////
+// Deterministic Pseudorandom int32_t Generator //
+//////////////////////////////////////////////////
+
+//Next seed in simple 32 bit integer deterministic psuedo-rand generator
+__host__ __device__ void dpr32_nextseed(int32_t *rseed_inout);
+
+//Simple 32 bit integer deterministic pseudo-random generator
+// *not* for cryptography
+// Frequency of generated floats should be uniform [0,1)
+__host__ __device__ float dpr32_randf(int32_t *rseed_inout);
+
+//box muller standard normal pseudorandom variable
+__host__ __device__ float dpr32_randnf(int32_t *rseed_inout);
+
+//////////////////////////////////////////////////
+// Deterministic Pseudorandom int64_t Generator //
+//////////////////////////////////////////////////
+
+//operates without side-effects on explicit seed for threaded use
+//deterministic pseudorandom number generator - takes seed and returns next seed
+__host__ __device__ void dpr64_nextseed(int64_t *seedinout);
+
+//deterministic pseudorandom number generator - takes seed and returns next seed
+//returns uniformly distributed double
+__host__ __device__ double dpr64_randd(int64_t *seedinout);
+
+__host__ __device__ float dpr64_randf(int64_t *seedinout);
+
+
+void test_dprg64();
+void test_dprg32();
+
+
+}; //end namespace amscuda
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/amscuarray.hpp

@@ -0,0 +1,47 @@
+#ifndef __CUARRAY_HPP__
+#define __CUARRAY_HPP__
+
+namespace amscuda
+{
+
+template<typename T> class cuarray
+{
+public:
+    int length;
+    T* data;
+
+    __device__ __host__ cuarray();
+    __device__ __host__ ~cuarray();
+
+    //Only call this on the device for thread/block local 
+    // dynamic arrays
+    __device__ __host__ int resize(const int _length);
+
+    __device__ __host__ int size() const;
+    __device__ __host__ T& at(const int I);
+    __device__ __host__ const T& at(const int I) const;
+
+    __device__ __host__ T& operator[](const int I);
+    __device__ __host__ const T& operator[](const int I) const;
+    
+    
+
+    __host__ int device_send(cuarray<T> **dptr);
+    __host__ int _device_send_overwrite(cuarray<T> **dptr);
+    __host__ int _device_send_copy(cuarray<T> *dptr);
+    
+    __host__ int device_pull(cuarray<T> *dptr);
+    __host__ int device_free(cuarray<T> **dptr);
+
+    __host__ int device_length(cuarray<T> *dptr);
+    __host__ T* device_data_ptr(cuarray<T> *dptr);
+    
+};
+
+void test_cuarray();
+
+};
+
+#include <amsculib2/amscuarray_impl.hpp>
+
+#endif

old/9apr26_prerefactor/include/amsculib2/amscuarray_dops.hpp

@@ -0,0 +1,76 @@
+#ifndef __AMSCUARRAY_DOPS_HPP__
+#define __AMSCUARRAY_DOPS_HPP__
+
+//Device Operations on Arrays
+//
+
+//Device Operations on Device Buffers
+// dodb
+
+namespace amscuda
+{
+
+
+    //sum
+    template<typename T> T devcuarray_sum(cuarray<T> *devptr);
+
+    template<typename T> T dbuff_sum(T *devbuffer, int N);
+
+
+    struct dbuff_statstruct
+    {
+        public:
+        float min;
+        float max;
+        float mean;
+        float stdev;
+        float sum;
+    };
+    
+    //stats (min,max,mean,stdev)
+
+    template<typename T> void dbuff_minmax(T *devbuffer, int N, T *min, T *max);
+
+    template<typename T> dbuff_statstruct dbuff_stats(T *devbuffer, int N); //
+
+    //sets all elements to setto
+    template<typename T> void dbuff_setall(T *devbuffer, int N, T setto, int nblocks, int nthreads);
+
+    //random device buffer functions
+    void dbuff_rand_dpr32(float *devbuffer, int N, int32_t *rseedinout,  int nblocks, int nthreads); //
+    void dbuff_rand_dpr32n(float *devbuffer, int N, int32_t *rseedinout,  int nblocks, int nthreads); //
+    
+
+    void dbuff_rand_dpr64(float *devbuffer, int N, int64_t *rseedinout,  int nblocks, int nthreads); //
+
+    //Elementwise device-buffer vector binary operation
+    //takes two input arrays ( , ) --> one output array
+    template<typename T1, typename T2, typename T3> void dbuff_vectorbinop(T1 *dbuf_a, T2 *dbuf_b, T3 *dbuf_out, int N, T3 (*fpnt)(T1,T2), int nblocks, int nthreads);
+
+    //Elementwise device-buffer vector two-parameter operation
+    //takes one input array, and a constant paramter ( ) ---> one output array
+    template<typename T1, typename T2, typename T3> void dbuff_vectorbinop(T1 *dbuf_a, T2 par_b, T3 *dbuf_out, int N, T3 (*fpnt)(T1,T2), int nblocks, int nthreads);
+
+
+    //vector_add
+    template<typename T> void dbuff_add(T *dbuff_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+    template<typename T> void dbuff_add(T *dbuff_a, T par_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+    template<typename T> void dbuff_sub(T *dbuff_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+    template<typename T> void dbuff_sub(T *dbuff_a, T par_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+    template<typename T> void dbuff_mult(T *dbuff_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+    template<typename T> void dbuff_mult(T *dbuff_a, T par_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+    template<typename T> void dbuff_div(T *dbuff_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+    template<typename T> void dbuff_div(T *dbuff_a, T par_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+    template<typename T> void dbuff_div(T par_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads); 
+
+
+    // Tests //
+
+    void test_dbuff_rand_dpr32();
+
+};
+
+#include <amsculib2/amscuarray_dops_impl.hpp>
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/amscuarray_dops_impl.hpp

@@ -0,0 +1,404 @@
+#ifndef __AMSCUARRAY_DOPS_IMPL_HPP__
+#define __AMSCUARRAY_DOPS_IMPL_HPP__
+
+namespace amscuda
+{
+
+template<typename T> __global__ void dbuff_sum_kf(T *devbuffer, int N, T *rets)
+{
+    int I0 = threadIdx.x + blockIdx.x*blockDim.x;
+    int Is = blockDim.x*gridDim.x;
+    int I;
+
+    T ret = (T) 0;
+    for(I=I0;I<N;I=I+Is)
+    {
+        ret = ret + devbuffer[I];
+    }
+    rets[I0] = ret;
+}
+
+template<typename T> T devcuarray_sum(cuarray<T> *devptr)
+{
+    T ret = T();
+    cudaError_t err = cudaSuccess;
+
+    cuarray<T> ldptr;
+    
+    cudaMemcpy(&ldptr,devptr,sizeof(cuarray<T>),cudaMemcpyDeviceToHost);
+    
+    ret = devbuffer_sum(ldptr.data,ldptr.length);
+
+    ldptr.data = NULL;
+    ldptr.length=0;
+
+    return ret;
+}
+
+template<typename T> T dbuff_sum(T *dbuff, int N)
+{
+    int I;
+    T ret = T();
+    cudaError_t err = cudaSuccess;
+
+    int nblocks;
+    int nthreads;
+
+    if(dbuff==NULL || N<=0)
+    {
+        return ret;
+    }
+
+    if(N>100)
+    {
+        nblocks = 10;
+        nthreads = (int)sqrt((float) (N/nblocks));
+        if(nthreads<=0) nthreads=1;
+        if(nthreads>512) nthreads=512;
+    }
+    else
+    {
+        nblocks = 1;
+        nthreads = 1;
+    }
+
+    T *rets = NULL;
+    T *devrets = NULL;
+
+    rets = new T[nblocks*nthreads];
+    cudaMalloc(&devrets,sizeof(T)*nblocks*nthreads);
+
+    dbuff_sum_kf<<<nblocks,nthreads>>>(dbuff,N,devrets);
+    cudaDeviceSynchronize();
+    err = cudaGetLastError();
+    if(err!=cudaSuccess)
+    {
+        printf("amscu::dbuff_sum error: %s\n",cudaGetErrorString(err));
+    }
+
+    cudaMemcpy(rets,devrets,sizeof(T)*nblocks*nthreads,cudaMemcpyDeviceToHost);
+
+    ret = (T)0;
+    for(I=0;I<nblocks*nthreads;I++)
+    {
+        ret = ret + rets[I];
+    }
+
+    cudaFree(devrets); devrets = NULL;
+    delete[] rets;
+
+    return ret;
+}
+
+
+template<typename T> __global__ void dbuff_minmax_kf(T *devbuffer, int N, T *maxs, T *mins)
+{
+    int I0 = threadIdx.x + blockIdx.x*blockDim.x;
+    int Is = blockDim.x*gridDim.x;
+    int I;
+
+    for(I=I0;I<N;I=I+Is)
+    {
+        if(I==I0)
+        {
+            maxs[I0] = devbuffer[I];
+            mins[I0] = devbuffer[I];
+        }
+        else
+        {
+            if(devbuffer[I]>maxs[I0])
+            {
+                maxs[I0] = devbuffer[I];
+            }
+            if(devbuffer[I]<mins[I0])
+            {
+                mins[I0] = devbuffer[I];
+            }
+        }
+    }
+
+    return;
+}
+
+template<typename T> void dbuff_minmax(T *devbuffer, int N, T *min, T *max)
+{
+    cudaError_t err = cudaSuccess;
+    int nblocks;
+    int nthreads;
+    int I;
+
+    T *maxs = NULL;
+    T *dev_maxs = NULL;
+    T *mins = NULL;
+    T *dev_mins = NULL;
+
+    T localmax = T(0);
+    T localmin = T(0);
+
+    if(devbuffer==NULL || N<=0)
+    {
+        if(min!=NULL) *min = T(0);
+        if(max!=NULL) *max = T(0);
+        return;
+    }
+
+    if(N>25)
+    {
+        nblocks = 25;
+        nthreads = (int) sqrt((float)(N/nblocks));
+        if(nthreads<1) nthreads = 1;
+        if(nthreads>512) nthreads = 512;
+    }
+    else
+    {
+        nblocks = 1;
+        nthreads = 1;
+    }
+
+    maxs = new T[nblocks*nthreads];
+    mins = new T[nblocks*nthreads];
+    cudaMalloc(&dev_maxs,nblocks*nthreads);
+    cudaMalloc(&dev_mins,nblocks*nthreads);
+    
+    dbuff_minmax_kf<<<nblocks,nthreads>>>(devbuffer,N,dev_maxs,dev_mins);
+    cudaDeviceSynchronize();
+    err = cudaGetLastError();
+    if(err!=cudaSuccess)
+    {
+        printf("amscu::dbuff_minmax error: %s\n",cudaGetErrorString(err));
+    }
+
+    cudaMemcpy(maxs,dev_maxs,sizeof(T)*nblocks*nthreads,cudaMemcpyDeviceToHost);
+    cudaMemcpy(mins,dev_mins,sizeof(T)*nblocks*nthreads,cudaMemcpyDeviceToHost);
+    
+
+    for(I=0;I<nblocks*nthreads;I++)
+    {
+        if(I==0)
+        {
+            localmax = maxs[0];
+            localmin = mins[0];
+        }
+        else
+        {
+            if(maxs[I]>localmax) localmax = maxs[I];
+            if(mins[I]<localmin) localmin = mins[I];
+        }
+    }
+
+    if(max!=NULL) *max = localmax;
+    if(min!=NULL) *min = localmin;
+
+    cudaFree(dev_maxs); dev_maxs = NULL;
+    cudaFree(dev_mins); dev_mins = NULL;
+    delete[] maxs; maxs = NULL;
+    delete[] mins; mins = NULL;
+
+    return;
+}
+
+template<typename T> __global__ void dbuff_setall_kf(T *devbuffer, int N, T setto)
+{
+    int I0 = threadIdx.x + blockIdx.x*blockDim.x;
+    int Is = blockDim.x*gridDim.x;
+    int I;
+
+    for(I=I0;I<N;I=I+Is)
+    {
+        devbuffer[I] = setto;
+    }
+    return;
+}
+
+template<typename T> void dbuff_setall(T *devbuffer, int N, T setto, int nblocks, int nthreads)
+{
+    cudaError_t err = cudaSuccess;
+
+    if(devbuffer==NULL || N<=0)
+    {
+        return;
+    }
+
+    dbuff_setall_kf<<<nblocks,nthreads>>>(devbuffer,N,setto);
+    cudaDeviceSynchronize();
+    err = cudaGetLastError();
+    if(err!=cudaSuccess)
+    {
+        printf("amscu::dbuff_setall error: %s\n",cudaGetErrorString(err));
+    }
+
+    return;
+}
+
+template<typename T1, typename T2, typename T3> __global__ void dbuff_vectorbinop_kf1(T1 *dbuf_a, T2 *dbuf_b, T3 *dbuf_out, int N, T3 (*fpnt)(T1,T2))
+{
+    int I0 = threadIdx.x + blockIdx.x*blockDim.x;
+    int Is = blockDim.x*gridDim.x;
+    int I;
+
+    T1 a;
+    T2 b;
+    T3 c;
+
+    for(I=I0;I<N;I=I+Is)
+    {
+        a = dbuf_a[I];
+        b = dbuf_b[I];
+        c = fpnt(a,b);
+        dbuf_out[I] = c;
+    }
+
+    return;
+}
+
+template<typename T1, typename T2, typename T3> __global__ void dbuff_vectorbinop_kf2(T1 *dbuf_a, T2 par_b, T3 *dbuf_out, int N, T3 (*fpnt)(T1,T2))
+{
+    int I0 = threadIdx.x + blockIdx.x*blockDim.x;
+    int Is = blockDim.x*gridDim.x;
+    int I;
+
+    T1 a;
+    T2 b;
+    T3 c;
+
+    for(I=I0;I<N;I=I+Is)
+    {
+        a = dbuf_a[I];
+        b = par_b;
+        c = fpnt(a,b);
+        dbuf_out[I] = c;
+    }
+
+    return;
+}
+
+
+//Elementwise device-buffer vector binary operation
+//takes two input arrays ( , ) --> one output array
+template<typename T1, typename T2, typename T3> void dbuff_vectorbinop(T1 *dbuf_a, T2 *dbuf_b, T3 *dbuf_out, int N, T3 (*fpnt)(T1,T2), int nblocks, int nthreads)
+{
+    cudaError_t err = cudaSuccess;
+
+    if(dbuf_a == NULL || dbuf_b == NULL || dbuf_out == NULL || N<=0)
+    {
+        return;
+    }
+
+    dbuff_vectorbinop_kf1<<<nblocks,nthreads>>>(dbuf_a,dbuf_b,dbuf_out,N);
+    cudaDeviceSynchronize();
+    err = cudaGetLastError();
+    if(err!=cudaSuccess)
+    {
+        printf("amscu::devbuffer_vectorbinop error: %s\n",cudaGetErrorString(err));
+    }
+
+    return;
+}
+
+//Elementwise device-buffer vector two-parameter operation
+//takes one input array, and a constant paramter ( ) ---> one output array
+template<typename T1, typename T2, typename T3> void dbuff_vectorbinop(T1 *dbuf_a, T2 par_b, T3 *dbuf_out, int N, T3 (*fpnt)(T1,T2), int nblocks, int nthreads)
+{
+    cudaError_t err = cudaSuccess;
+
+    if(dbuf_a == NULL || dbuf_out == NULL || N<=0)
+    {
+        return;
+    }
+
+    dbuff_vectorbinop_kf2<<<nblocks,nthreads>>>(dbuf_a,par_b,dbuf_out,N);
+    cudaDeviceSynchronize();
+    err = cudaGetLastError();
+    if(err!=cudaSuccess)
+    {
+        printf("amscu::devbuffer_vectorbinop error: %s\n",cudaGetErrorString(err));
+    }
+
+    return;
+}
+
+template<typename T> T dbuff_add_fn(T a, T b)
+{
+    return a+b;
+}
+
+template<typename T> void dbuff_add(T *dbuff_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_a,dbuff_b,dbuff_out,N,&dbuff_add_fn,nblocks,nthreads);
+    return;
+}
+
+template<typename T> void dbuff_add(T *dbuff_a, T par_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_a,par_b,dbuff_out,N,&dbuff_add_fn,nblocks,nthreads);
+    return;
+}
+
+template<typename T> T dbuff_sub_fn(T a, T b)
+{
+    return a-b;
+}
+
+template<typename T> void dbuff_sub(T *dbuff_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_a,dbuff_b,dbuff_out,N,&dbuff_sub_fn,nblocks,nthreads);
+    return;
+}
+
+template<typename T> void dbuff_sub(T *dbuff_a, T par_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_a,par_b,dbuff_out,N,&dbuff_sub_fn,nblocks,nthreads);
+    return;
+}
+
+template<typename T> T dbuff_mult_fn(T a, T b)
+{
+    return a*b;
+}
+
+template<typename T> void dbuff_mult(T *dbuff_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_a,dbuff_b,dbuff_out,N,&dbuff_mult_fn,nblocks,nthreads);
+    return;
+}
+
+template<typename T> void dbuff_mult(T *dbuff_a, T par_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_a,par_b,dbuff_out,N,&dbuff_mult_fn,nblocks,nthreads);
+    return;
+}
+
+template<typename T> T dbuff_div_fn(T a, T b)
+{
+    return a/b;
+}
+
+template<typename T> void dbuff_div(T *dbuff_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_a,dbuff_b,dbuff_out,N,&dbuff_div_fn,nblocks,nthreads);
+    return;
+}
+
+template<typename T> void dbuff_div(T *dbuff_a, T par_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_a,par_b,dbuff_out,N,&dbuff_div_fn,nblocks,nthreads);
+    return;
+}
+
+template<typename T> T dbuff_ldiv_fn(T a, T b)
+{
+    return b/a;
+}
+
+
+template<typename T> void dbuff_div(T par_a, T *dbuff_b, T *dbuff_out, int N, int nblocks, int nthreads)
+{
+    dbuff_vectorbinop(dbuff_b,par_a,dbuff_out,N,&dbuff_ldiv_fn,nblocks,nthreads);
+    return;
+}
+
+
+};
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/amscuarray_impl.hpp

@@ -0,0 +1,323 @@
+#ifndef __CUARRAY_IMPL_HPP__
+#define __CUARRAY_IMPL_HPP__
+
+namespace amscuda
+{
+
+// New Version cuarray<T>
+// simpler, less crap going on
+
+template<typename T> __device__ __host__ cuarray<T>::cuarray()
+{
+    length = 0;
+    data = NULL;
+}
+
+template<typename T> __device__ __host__ cuarray<T>::~cuarray()
+{
+    if(data!=NULL)
+    {
+        delete[] data; data = NULL;
+    }
+    length = 0;
+}
+
+template<typename T> __device__ __host__ int cuarray<T>::resize(const int _length)
+{
+    int ret = 0;
+
+    T *newbuffer = NULL;
+
+    if(length==_length)
+    {
+        //do nothing
+        ret = 1;
+        return ret;
+    }
+    if(_length<=0)
+    {
+        if(data!=NULL)
+        {
+            delete[] data;
+            data = NULL;
+        }
+        length = 0;
+        ret = 1;
+    }
+
+    newbuffer = new T[_length];
+    if(newbuffer==NULL)
+    {
+        ret = -1; //failed to allocate memory
+        return ret;
+    }
+
+    int I;
+    T def;
+
+    if(data!=NULL)
+    {
+        for(I=0;I<length&&I<_length;I++)
+        {
+            newbuffer[I] = data[I];
+        }
+        for(I=length;I<_length;I++)
+        {
+            newbuffer[I] = def;
+        }
+        delete[] data; data=NULL;
+    }
+    else
+    {
+        for(I=0;I<_length;I++)
+        {
+            newbuffer[I] = def;
+        }
+    }
+
+    data = newbuffer;
+    length = _length;
+    ret = 1;
+
+    return ret;
+}
+
+
+template<typename T> __host__ int cuarray<T>::device_send(cuarray<T> **dptr)
+{
+    int ret = 0;
+    int dlength;
+
+    if(*dptr==NULL)
+    {
+        ret = _device_send_overwrite(dptr);
+    }
+    else
+    {
+        dlength = device_length(*dptr);
+        if(dlength=length)
+        {
+            ret = _device_send_copy(*dptr);
+        }
+        else
+        {
+            ret = _device_send_overwrite(dptr);
+        }
+    }
+
+    return ret;
+}
+
+template<typename T> __host__ int cuarray<T>::_device_send_overwrite(cuarray<T> **dptr)
+{
+    int ret = 0;
+    cuarray<T> dlocal;
+    cudaError_t err = cudaSuccess;
+    device_free(dptr);
+    
+    if(length>=0 && data!=NULL)
+    {
+        err = cudaMalloc(dptr,sizeof(cuarray<T>));
+        if(err==cudaSuccess)
+        {
+            err = cudaMalloc(&(dlocal.data),sizeof(T)*length);
+            dlocal.length = length;
+
+            if(err==cudaSuccess)
+            {
+                cudaMemcpy(*dptr,&dlocal,sizeof(cuarray<T>),cudaMemcpyHostToDevice);
+                if(data!=NULL)
+                    err = cudaMemcpy(dlocal.data,data,sizeof(T)*length,cudaMemcpyHostToDevice);
+                else
+                    err = cudaSuccess;
+                if(err==cudaSuccess)
+                {
+                    ret = 1;
+                }
+                else
+                {
+                    ret = -3;
+                }
+            }
+            else
+            {
+                ret = -2;
+            }
+        }
+        else
+        {
+            ret = -1;
+        }
+    }
+    else
+    {
+        dlocal.data = NULL;
+        dlocal.length = 0;
+        err = cudaMalloc(dptr,sizeof(cuarray<T>));
+        if(err==cudaSuccess)
+        {
+            cudaMemcpy(*dptr,&dlocal,sizeof(cuarray<T>),cudaMemcpyHostToDevice);
+            ret = 1;
+        }
+        else
+        {
+            ret = -4;
+        }
+    }
+
+
+    dlocal.data = NULL;
+    dlocal.length = -1;
+
+    return ret;
+}
+
+template<typename T> __host__ int cuarray<T>::_device_send_copy(cuarray<T> *dptr)
+{
+    int ret = 0;
+    cudaError_t err = cudaSuccess;
+    T* ddata = NULL;
+    ddata = device_data_ptr(dptr);
+
+    err = cudaMemcpy(ddata,data,sizeof(T)*length,cudaMemcpyHostToDevice);
+    if(err==cudaSuccess)
+    {
+        ret = 1;
+    }
+    else
+    {
+        ret = -1;
+    }
+
+    return ret;
+}
+
+template<typename T> __host__ int cuarray<T>::device_pull(cuarray<T> *dptr)
+{
+    int ret = 0;
+    int dlength;
+    T* ddata;
+    cudaError_t err;
+
+    if(dptr==NULL) 
+    {
+        ret = -1; // null d pointer
+        return ret;
+    }
+
+    dlength = device_length(dptr);
+    if(dlength!=length)
+    {
+        this->resize(dlength);
+    }
+
+    ddata = device_data_ptr(dptr);
+
+    if(length>0 && data!=NULL && ddata!=NULL)
+    {
+        err = cudaMemcpy(data,dptr,length*sizeof(T),cudaMemcpyDeviceToHost);
+        if(err==cudaSuccess)
+        {
+            ret = 1;
+        }
+        else
+        {
+            ret = -2;
+        }
+    }
+
+    return ret;
+}
+
+template<typename T> __host__ int cuarray<T>::device_free(cuarray<T> **dptr)
+{
+    int ret = 0;
+    cuarray<T> dlocal;
+
+    if(*dptr!=NULL)
+    {
+        cudaMemcpy(&dlocal,dptr,sizeof(cuarray<T>),cudaMemcpyDeviceToHost);
+        if(dlocal.data!=NULL)
+        {
+            cudaFree(dlocal.data);
+            dlocal.data = NULL;
+        }
+
+        cudaFree(*dptr);
+        *dptr = NULL;
+        ret = 1;
+    }
+
+    dlocal.data = NULL;
+    dlocal.length = -1;
+
+    return ret;
+}
+
+template<typename T> __host__ int cuarray<T>::device_length(cuarray<T> *dptr)
+{
+    int ret = -1;
+    cuarray<T> dlocal;
+
+    if(dptr==NULL)
+    {
+        return ret;
+    }
+
+    cudaMemcpy(&dlocal,dptr,sizeof(cuarray<T>),cudaMemcpyDeviceToHost);
+    ret = dlocal.length;
+
+    dlocal.data = NULL;
+    dlocal.length = -1;
+
+    return ret;
+}
+
+template<typename T> __host__ T* cuarray<T>::device_data_ptr(cuarray<T> *dptr)
+{
+    T* ret = NULL;
+    cuarray<T> dlocal;
+
+    if(dptr==NULL)
+    {
+        return ret;
+    }
+
+    cudaMemcpy(&dlocal,dptr,sizeof(cuarray<T>),cudaMemcpyDeviceToHost);
+    ret = dlocal.data;
+
+    dlocal.data = NULL;
+    dlocal.length = -1;
+
+    return ret;
+}
+
+template<typename T>  __device__ __host__ int cuarray<T>::size() const
+{
+    return this->length;
+}
+
+template<typename T>  __device__ __host__ T& cuarray<T>::at(const int I)
+{
+    return this->data[I];
+}
+
+template<typename T>  __device__ __host__ const T& cuarray<T>::at(const int I) const
+{
+    return this->data[I];
+}
+
+template<typename T> __device__ __host__ T& cuarray<T>::operator[](const int I)
+{
+    return this->data[I];
+}
+
+template<typename T> __device__ __host__ const T& cuarray<T>::operator[](const int I) const
+{
+    return this->data[I];
+}
+
+};
+
+
+#endif

old/9apr26_prerefactor/include/amsculib2/amscuda_binarrrw.hpp

@@ -0,0 +1,19 @@
+#ifndef __AMSCUDA_BINARRRW_HPP__
+#define __AMSCUDA_BINARRRW_HPP__
+
+namespace amscuda
+{
+
+template<typename T> int fread_ndarray(FILE *fp, cuarray<int32_t> *shape, cuarray<T> *buffer);
+template<typename T> int fwrite_ndarray(FILE *fp, const cuarray<int32_t> *shape, const cuarray<T> *buffer);
+
+template<typename T> int fwrite_buffer(FILE *fp, const int N, const T *buffer);
+template<typename T> int fread_buffer(FILE *fp, const int Nmax, const T *buffer);
+
+
+}; //end namespace amscuda
+
+#include <amsculib2/amscuda_binarrrw_impl.hpp>
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/amscuda_binarrrw_impl.hpp

@@ -0,0 +1,194 @@
+#ifndef __AMSCUDA_BINARRRW_IMPL_HPP__
+#define __AMSCUDA_BINARRRW_IMPL_HPP__
+
+namespace amscuda
+{
+
+template<typename T> int fread_ndarray(FILE *fp, cuarray<int32_t> *shape, cuarray<T> *buffer)
+{
+    int ret = 1;
+    int I;
+    long piprod;
+    int32_t q;
+    int cnt;
+
+    int32_t Nd;
+
+    if(fp!=NULL)
+    {
+        if(!feof(fp))
+        {
+            cnt = fread(&Nd,sizeof(int32_t),1,fp);
+            if(Nd>0 && cnt>0)
+            {
+                shape->resize(Nd);
+                piprod = 1;
+                for(I=0;I<Nd;I++)
+                {
+                    cnt = fread(&q,sizeof(int32_t),1,fp);
+                    shape->at(I) = q;
+                    if(q>0)
+                    {
+                        piprod = piprod*q;
+                    }
+                    else
+                    {
+                        piprod = 0;
+                    }
+                }
+
+                buffer->resize(piprod);
+                if(piprod>0)
+                {
+                    cnt = fread((buffer->data),sizeof(T),piprod,fp);
+                    if(piprod==cnt)
+                    {
+                        ret = 1;
+                    }
+                    else
+                    {
+                        printf("fread_ndarray, read %d values, expecting %ld\n",cnt,piprod);
+                        ret = 0;
+                    }
+                }
+            }
+            else
+            {
+                printf("fread_ndarray: Read a number of dimensions<=0.\n");
+                Nd = 0;
+                shape->resize(0);
+                buffer->resize(0);
+            }
+        }
+        else
+        {
+            printf("fread_ndarray: fp=NULL.\n");
+            ret = 0;
+        }
+    }
+    else
+    {
+        ret = 0;
+    }
+
+    return ret;
+}
+
+template<typename T> int fwrite_ndarray(FILE *fp, const cuarray<int32_t> *shape, const cuarray<T> *buffer)
+{
+    int ret = 1;
+    long piprod;
+    int I;
+    int32_t Nd;
+
+    if(fp==NULL)
+    {
+        ret = 0;
+        printf("fwrite_ndarray: fp=NULL\n");
+        return ret;
+    }
+
+    piprod = 1;
+    for(I=0;I<shape->size();I++)
+    {
+        if(shape->at(I)>0)
+        {
+            piprod = piprod*shape->at(I);
+        }
+        else
+        {
+            piprod = 0;
+        }
+    }
+
+    Nd = (int32_t) shape->size();
+
+    if(piprod!=buffer->size())
+    {
+        ret = 0;
+        printf("fwrite_ndarray: buffer is size %ld, while shape is size %ld\n",(long)buffer->size(),(long)piprod);
+        return ret;
+    }
+
+    fwrite(&Nd,sizeof(int32_t),1,fp);
+    if(Nd>0)
+    {
+        fwrite(shape->data,sizeof(int32_t),Nd,fp);
+        if(piprod>0)
+        {
+            fwrite(buffer->data,sizeof(T),buffer->size(),fp);
+        }
+    }
+
+    return ret;
+}
+
+template<typename T> int fwrite_buffer(FILE *fp, const int N, const T *buffer)
+{
+    int ret = 0;
+    int Nd = 1;
+    
+    if(fp==NULL)
+    {
+        ret = 0;
+        printf("fwrite_buffer: fp=NULL\n");
+        return ret;
+    }
+
+    fwrite(&Nd,sizeof(int32_t),1,fp);
+    fwrite(&N,sizeof(int32_t),1,fp);
+    fwrite(buffer,sizeof(T),N,fp);
+    
+    return ret;
+}
+
+template<typename T> int fread_buffer(FILE *fp, const int Nmax, const T *buffer)
+{
+    int ret = 0;
+    
+    int cnt;
+    int32_t Nd;
+    int32_t *dims = NULL;
+    int piprod;
+    int32_t q;
+    int I;
+
+    int Nr;
+
+
+    if(fp==NULL) {ret = -1; return ret;}
+    if(feof(fp)) {ret = -2; return ret;}
+
+    cnt = fread(&Nd,sizeof(int32_t),1,fp);
+    if(Nd>0 && cnt>0)
+    {
+        piprod = 1;
+        dims = new(std::nothrow) int32_t[Nd];
+        for(I=0;I<Nd;I++)
+        {
+            cnt = fread(&q,sizeof(int32_t),1,fp);
+            dims[I] = q;
+            piprod = piprod*dims[I];
+            if(piprod==cnt)
+            {
+                ret = 1;
+            }
+            else
+            {
+                printf("fwrite_buffer, read %d values, expecting %d\n",cnt,piprod);
+            }
+        }
+
+        Nr = amscuda::min<int32_t>(Nmax,piprod);
+        cnt = fread(buffer,sizeof(T),Nr,fp);
+    }
+
+    if(dims!=NULL) {delete[] dims; dims=NULL;}
+
+    return ret;
+}
+
+}; //end namespace amscuda
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/amscugeom.hpp

@@ -0,0 +1,11 @@
+#ifndef __AMSCUGEOM_HPP__
+#define __AMSCUGEOM_HPP__
+
+namespace amscuda
+{
+
+
+}; //end namespace amscuda
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/amsculib2.hpp

@@ -0,0 +1,70 @@
+#ifndef __AMSCULIB2_HPP__
+#define __AMSCULIB2_HPP__
+
+//Std Lib Includes
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <stdint.h>
+#include <time.h>
+#include <new>
+
+#include <cuda_runtime_api.h> //where all the cuda functions live
+#include <cuda_runtime.h>
+#include <cuda.h>
+
+//Dependencies
+
+//Predeclarations
+class cuvect2;
+class cuvect3;
+class cuvect4;
+class cuvect2f;
+class cuvect3f;
+class cuvect4f;
+
+//Need a way to define the same symbols using both host and device code
+//A solution was found here: https://stackoverflow.com/questions/9457572/cuda-host-and-device-using-same-constant-memory
+#ifdef __CUDA_ARCH__
+#define AMSCU_CONST __constant__
+#else
+#define AMSCU_CONST
+#endif
+
+namespace amscuda
+{
+
+    //default thread and block execution
+    AMSCU_CONST static const int amscu_defnblocks = 256;
+    AMSCU_CONST static const int amscu_defnthreads = 512;
+
+    //default numthreads to execute on cpu
+    AMSCU_CONST static const int amscu_defcputhreads = 8;
+
+}; //end namespace amscuda
+
+//Components
+#include <amsculib2/amscu_cudafunctions.hpp>
+#include <amsculib2/amscumath.hpp>
+#include <amsculib2/amscu_comp64.hpp>
+#include <amsculib2/amscu_comp128.hpp>
+#include <amsculib2/cuvect2.hpp>
+#include <amsculib2/cuvect3.hpp>
+#include <amsculib2/cuvect4.hpp>
+#include <amsculib2/cuvect2f.hpp>
+#include <amsculib2/cuvect3f.hpp>
+#include <amsculib2/cuvect4f.hpp>
+#include <amsculib2/amscugeom.hpp>
+#include <amsculib2/amscuarray.hpp>
+#include <amsculib2/amscuda_binarrrw.hpp>
+#include <amsculib2/amscu_random.hpp>
+
+#include <amsculib2/amscuarray_dops.hpp>
+
+#include <amsculib2/amscurarray.cuh>
+
+
+
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/amscumath.hpp

@@ -0,0 +1,56 @@
+#ifndef __AMSCUMATH_HPP__
+#define __AMSCUMATH_HPP__
+
+namespace amscuda
+{
+
+
+
+    //Problem: These are not in the namespace
+    //#define nan NAN
+    //#define fnan (float) NAN
+    //#define inf INFINITY
+    //#define finf (float) INFINITY
+    //#define pi 3.1415926535897936
+
+    //These need to be the same symbol for both host and device code
+    AMSCU_CONST static const double nan = NAN;
+    AMSCU_CONST static const float fnan = (float) NAN;
+    AMSCU_CONST static const double inf = INFINITY;
+    AMSCU_CONST static const float finf = (float) INFINITY;
+    AMSCU_CONST static const double pi = 3.1415926535897936;
+    AMSCU_CONST static const float pif = 3.1415926535897936;
+
+    __host__ __device__ double dabs(double x);
+    __host__ __device__ float fabs(float x);
+    
+    template<typename T> __host__ __device__ T abs(const T in)
+    {
+        T ret = in;
+        if(in<0) ret = -in;
+        return ret;
+    }
+
+    __host__ __device__ double mod(double a, double md);
+    __host__ __device__ float mod(float a, float md);
+    __host__ __device__ int mod(int x, int n);
+    __host__ __device__ long mod(long x, long n);
+
+    __host__ __device__ int truediv(int x, int y);
+    __host__ __device__ long truediv(long x, long y);
+
+    template<typename T> __host__ __device__ T min(T a, T b);
+    template<typename T> __host__ __device__ T max(T a, T b);
+
+    __device__ __host__ double arg(double x, double y);
+    __device__ __host__ void get_azel(double x, double y, double z, double *az, double *el);
+
+    void test_amscumath1();
+
+
+}; //end namespace amscuda
+
+#include <amsculib2/amscumath_impl.hpp>
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/amscumath_impl.hpp

@@ -0,0 +1,42 @@
+#ifndef __AMSCUMATH_IMPL_HPP__
+#define __AMSCUMATH_IMPL_HPP__
+
+namespace amscuda
+{
+
+template<typename T> __host__ __device__ T min(T a, T b)
+{
+    if(a>b)
+    {
+        return b;
+    }
+    else
+    {
+        return a;
+    }
+    return a;
+}
+
+template<typename T> __host__ __device__ T max(T a, T b)
+{
+    if(a>b)
+    {
+        return a;
+    }
+    else
+    {
+        return b;
+    }
+    return a;
+}
+
+template<> __host__ __device__ double min(double a, double b);
+template<> __host__ __device__ float min(float a, float b);
+template<> __host__ __device__ double max(double a, double b);
+template<> __host__ __device__ float max(float a, float b);
+
+
+}; //end namespace amscuda
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/amscurarray.cuh

@@ -0,0 +1,66 @@
+#ifndef __AMSCURARRAY_HPP__
+#define __AMSCURARRAY_HPP__
+
+namespace amscuda
+{
+
+//Cuda ragged array class
+template<typename T> class curarray
+{
+public:
+    int device;
+    curarray* devptr; //pointer to mirror class on the device
+
+    int Narrays; //number of arrays
+    
+    int *N; //dimension of each array
+    T** hostarrayptrs; //pointers to each array on the host - null on the device
+    T** devarrayptrs; //pointers to each array on the device
+        //the double pointer is a host pointer to device pointers on the host class
+        //for the device class, only the second set of arrays is in use
+    
+    //the constructor and destructor set all pointers to NULL, they
+    // do *not* manage memory. This is done with curarray_new and curarray_delete
+    __device__ __host__ curarray();
+    __device__ __host__ ~curarray();
+    
+    __host__ int push();
+    __host__ int pull();
+    //__device__ int dev_resizearray(int arraynum, int arraysize);
+    __host__ int resizearray(int arraynum, int arraysize);
+    // I may want a way to resize arrays on the device without pushing/pulling all the array contents
+    
+
+};
+
+template<typename T> int curarray_new(curarray<T>** ptr, int Narrays);
+
+template<typename T> int curarray_delete(curarray<T>** ptr);
+
+template<typename T> int curarray_device_new(curarray<T> *hostptr);
+
+template<typename T> int curarray_device_delete(curarray<T> *hostptr);
+
+template<typename T> int curarray_push(curarray<T> *hostptr);
+
+template<typename T> int curarray_pull(curarray<T> *hostptr);
+
+
+//template<typename T> int curarray_host_fillall(curarray<T> *hostptr, const T &val);
+//template<typename T> int curarray_device_fillall(curarray<T> *hostptr, const T &val);
+
+
+//template<typename T> __host__ int curarray_deletearray(curarray<T> *hostptr, int arrayindex);
+//template<typename T> __device__ int curarray_dev_deletearray(curarray<T> *devptr, int arrayindex);
+
+//template<typename T> __host__ int curarray_allocarray(curarray<T> *hostptr, int arrayindex, int size);
+//template<typename T> __device__ int curarray_dev_allocarray(curarray<T> *devptr, int arrayindex, int size);
+
+
+void test_amscurarray1();
+
+};
+
+#include <amsculib2/amscurarray_impl.cuh>
+
+#endif

old/9apr26_prerefactor/include/amsculib2/amscurarray_impl.cuh

@@ -0,0 +1,529 @@
+#ifndef __AMSCURARRAY_IMPL_HPP__
+#define __AMSCURARRAY_IMPL_HPP__
+
+namespace amscuda
+{
+
+template<typename T> curarray<T>::curarray()
+{
+    device = -1;
+    devptr = NULL;
+    Narrays = 0;
+    N = NULL;
+    hostarrayptrs = NULL;
+    devarrayptrs = NULL;
+
+}
+
+template<typename T> curarray<T>::~curarray()
+{
+    device = -1;
+    devptr = NULL;
+    Narrays = 0;
+    N = NULL;
+    hostarrayptrs = NULL;
+    devarrayptrs = NULL;
+    
+}
+
+template<typename T> int curarray_new(curarray<T>** ptr, int Narrays)
+{
+    int ret = 0;
+    int device;
+    curarray<T> *lhptr = *ptr;
+
+    cudaGetDevice(&device);
+
+    if(lhptr!=NULL)
+    {
+        curarray_delete(ptr);
+    }
+
+    *ptr = new(std::nothrow) curarray<T>();
+    lhptr = *ptr;
+
+    int I;
+
+    if(Narrays<0) Narrays=0;
+
+    lhptr->Narrays = Narrays;
+    lhptr->device = device;
+    lhptr->N = new(std::nothrow) int[Narrays];
+    lhptr->hostarrayptrs = new(std::nothrow) T*[Narrays];
+    lhptr->devarrayptrs = new(std::nothrow) T*[Narrays];
+    
+    for(I=0;I<Narrays;I++)
+    {
+        lhptr->N[I] = 0;
+        lhptr->hostarrayptrs[I] = NULL;
+        lhptr->devarrayptrs[I] = NULL;
+    }
+
+    curarray_device_new(lhptr);
+
+    return ret;
+}
+
+template<typename T> int curarray_delete(curarray<T>** ptr)
+{
+    int ret = 0;
+    curarray<T> *lptr = NULL;
+    int olddev;
+
+    int I;
+
+    if(*ptr==NULL)
+    {
+        return 0;
+    }
+    lptr = *ptr;
+
+    cudaGetDevice(&olddev);
+    cudaSetDevice(lptr->device);
+
+    
+    if(lptr->devptr!=NULL)
+    {
+        curarray_device_delete(lptr);
+    }
+
+    lptr->device = -1;
+
+    for(I=0;I<lptr->Narrays;I++)
+    {
+        if(lptr->hostarrayptrs!=NULL)
+        {
+            if(lptr->hostarrayptrs[I]!=NULL)
+            {
+                delete[] lptr->hostarrayptrs[I];
+                lptr->hostarrayptrs[I] = NULL;
+            }
+        }
+
+        if(lptr->devarrayptrs!=NULL)
+        {
+            if(lptr->devarrayptrs[I]!=NULL)
+            {
+                //erasing device memory should have been handled in curarray_device_delete
+                lptr->devarrayptrs[I] = NULL;
+            }
+        }
+
+        lptr->N[I] = 0;
+    }
+
+    if(lptr->N != NULL) {delete[] lptr->N; lptr->N = NULL;}
+    if(lptr->hostarrayptrs!=NULL) {delete[] lptr->hostarrayptrs; lptr->hostarrayptrs=NULL;}
+    if(lptr->devarrayptrs!=NULL) {delete[] lptr->devarrayptrs; lptr->devarrayptrs=NULL;}
+
+    if(*ptr!=NULL) {delete *ptr; *ptr = NULL;}
+
+    cudaSetDevice(olddev);
+
+    return ret;
+}
+
+template<typename T> int curarray_device_new(curarray<T> *hostptr)
+{
+    int ret = 0;
+    curarray<T> ldevdata;
+
+    if(hostptr==NULL) return -1;
+    if(hostptr->devptr!=NULL)
+    {
+        curarray_device_delete(hostptr);
+    }
+
+    int I;
+
+    cudaGetDevice(&(hostptr->device));
+
+    ldevdata.device = hostptr->device;
+    ldevdata.Narrays = hostptr->Narrays;
+    int Narrays = hostptr->Narrays;
+
+    cudaMalloc(&(ldevdata.N),sizeof(int)*Narrays);
+    cudaMemcpy(ldevdata.N,hostptr->N,sizeof(int)*Narrays,cudaMemcpyHostToDevice);
+
+    ldevdata.hostarrayptrs = NULL;
+
+    for(I=0;I<Narrays;I++)
+    {
+        if(hostptr->N[I]>0)
+        {
+            if(hostptr->devarrayptrs[I]!=NULL)
+            {
+                cudaFree(hostptr->devarrayptrs[I]);
+                hostptr->devarrayptrs[I] = NULL;
+            }
+
+            cudaMalloc(&(hostptr->devarrayptrs[I]),sizeof(T)*hostptr->N[I]);
+            cudaMemcpy(hostptr->devarrayptrs[I],hostptr->hostarrayptrs[I],sizeof(T)*hostptr->N[I],cudaMemcpyHostToDevice);
+        }
+        else
+        {
+            if(hostptr->devarrayptrs[I]!=NULL)
+            {
+                cudaFree(hostptr->devarrayptrs[I]);
+                hostptr->devarrayptrs[I] = NULL;
+            }
+        }
+    }
+
+    cudaMalloc(&(ldevdata.devarrayptrs),sizeof(T*)*Narrays);
+    cudaMemcpy(ldevdata.devarrayptrs,hostptr->devarrayptrs,sizeof(T*)*Narrays,cudaMemcpyHostToDevice);
+
+    cudaMalloc(&(hostptr->devptr),sizeof(curarray<T>));
+    cudaMemcpy(hostptr->devptr,&ldevdata,sizeof(curarray<T>),cudaMemcpyHostToDevice);    
+    
+    ret = 1;
+
+    return ret;
+}
+
+template<typename T> int curarray_device_delete(curarray<T> *hostptr)
+{
+    int ret = 0;
+
+    curarray<T> ldevdata;
+    int olddev;
+
+    if(hostptr->devptr==NULL)
+    {
+        return 0;
+    }
+
+    cudaGetDevice(&olddev);
+    cudaSetDevice(hostptr->device);
+
+    cudaMemcpy(&ldevdata,hostptr->devptr,sizeof(curarray<T>),cudaMemcpyDeviceToHost);
+
+    int I;
+    for(I=0;I<hostptr->Narrays;I++)
+    {
+        if(hostptr->devarrayptrs[I]!=NULL)
+        {
+            cudaFree(hostptr->devarrayptrs[I]);
+            hostptr->devarrayptrs[I] = NULL;
+        }
+    }
+
+    cudaFree(ldevdata.devarrayptrs);
+    cudaFree(ldevdata.N);
+
+    cudaFree(hostptr->devptr);
+    hostptr->devptr = NULL;
+    hostptr->device = -1;
+
+    cudaSetDevice(olddev);
+
+    ret = 1;
+    return ret;
+}
+
+template<typename T> int curarray_push(curarray<T> *hostptr)
+{
+    int ret = 0;
+
+    int olddev;
+
+    curarray<T> ldevdata;
+    T** ldevarrayptrs = NULL;
+    int *devN = NULL;
+
+    if(hostptr==NULL) return -1;
+
+    cudaGetDevice(&olddev);
+    cudaSetDevice(hostptr->device);
+
+    int Narrays = hostptr->Narrays;
+
+    cudaMemcpy(&ldevdata,hostptr->devptr,sizeof(curarray<T>),cudaMemcpyDeviceToHost);
+    ldevarrayptrs = new(std::nothrow) T*[Narrays];
+    devN = new(std::nothrow) int[Narrays];
+
+    cudaMemcpy(ldevarrayptrs,ldevdata.devarrayptrs,sizeof(T*)*Narrays,cudaMemcpyDeviceToHost);
+    cudaMemcpy(devN,ldevdata.N,sizeof(int)*Narrays,cudaMemcpyDeviceToHost);
+
+
+    int I;
+
+    for(I=0;I<Narrays;I++)
+    {
+        //check to see that host size is the same as device size, and that 
+        //the host device pointer is the same as the device device pointer
+        if( (hostptr->N[I]!=devN[I]) || 
+            (hostptr->devarrayptrs[I] != ldevarrayptrs[I])
+        )
+        {
+            cudaFree(ldevarrayptrs[I]);
+            ldevarrayptrs[I] = NULL;
+            hostptr->devarrayptrs[I] = NULL;
+
+            if(hostptr->N[I]>0)
+            {
+                cudaMalloc(&(hostptr->devarrayptrs[I]),sizeof(T)*hostptr->N[I]);
+                ldevarrayptrs[I] = hostptr->devarrayptrs[I];
+                devN[I] = hostptr->N[I];
+            }
+            else
+            {
+                devN[I] = 0;
+            }
+        }
+        if(hostptr->N[I]>0)
+        {
+            //copy host data to device
+            cudaMemcpy(hostptr->devarrayptrs[I],hostptr->hostarrayptrs[I],sizeof(T)*hostptr->N[I],cudaMemcpyHostToDevice);
+        }
+    } //for each array
+
+    //rectify and copy device data structure to device
+    ldevdata.device = hostptr->device;
+    ldevdata.devptr = NULL;
+    ldevdata.Narrays = hostptr->Narrays; //later - logic for dealing with when this is not true
+    ldevdata.hostarrayptrs = NULL;
+
+
+    cudaMemcpy(ldevdata.N,hostptr->N,sizeof(int)*Narrays,cudaMemcpyHostToDevice);
+    cudaMemcpy(ldevdata.devarrayptrs,hostptr->devarrayptrs,sizeof(T*)*Narrays,cudaMemcpyHostToDevice);
+
+    cudaMemcpy(hostptr->devptr,&ldevdata,sizeof(curarray<T>),cudaMemcpyHostToDevice);
+
+    cuda_errortrap("curarray_push cuda error:");
+
+    cudaSetDevice(olddev);
+
+    delete[] ldevarrayptrs;
+    delete[] devN;
+
+    return ret;
+}
+
+template<typename T> int curarray_pull(curarray<T> *hostptr)
+{
+    int ret = 0;
+
+    int olddev;
+
+    curarray<T> ldevdata;
+    T** ldevarrayptrs = NULL;
+    int *devN = NULL;
+
+    if(hostptr==NULL) return -1;
+
+    cudaGetDevice(&olddev);
+    cudaSetDevice(hostptr->device);
+
+    cuda_errortrap("dbg1");
+
+    int Narrays = hostptr->Narrays;
+
+    cudaMemcpy(&ldevdata,hostptr->devptr,sizeof(curarray<T>),cudaMemcpyDeviceToHost);
+    ldevarrayptrs = new(std::nothrow) T*[Narrays];
+    devN = new(std::nothrow) int[Narrays];
+
+    cuda_errortrap("dbg2");
+
+    cudaMemcpy(ldevarrayptrs,ldevdata.devarrayptrs,sizeof(T*)*Narrays,cudaMemcpyDeviceToHost);
+    cudaMemcpy(devN,ldevdata.N,sizeof(int)*Narrays,cudaMemcpyDeviceToHost);
+
+    cuda_errortrap("dbg3");
+    char dbgjnk[50];
+
+    int I;
+    for(I=0;I<Narrays;I++)
+    {
+        //check to see that host size is the same as device size, and that 
+        //the host device pointer is the same as the device device pointer
+        if(hostptr->devarrayptrs[I] != ldevarrayptrs[I])
+        {
+            hostptr->devarrayptrs[I] = ldevarrayptrs[I];
+        }
+
+        if(hostptr->N[I]!=devN[I])
+        {
+            if(hostptr->hostarrayptrs[I]!=NULL) 
+            {
+                delete[] hostptr->hostarrayptrs[I];
+                hostptr->hostarrayptrs[I] = NULL;
+            }
+
+            if(devN[I]>0)
+            {
+                hostptr->hostarrayptrs[I] = new(std::nothrow) T[devN[I]];
+                hostptr->N[I] = devN[I];
+            }
+            else
+            {
+                hostptr->N[I] = 0;
+            }
+        }
+
+        if(hostptr->hostarrayptrs[I]!=NULL && hostptr->devarrayptrs[I]!=NULL)
+        {
+            cudaMemcpy(hostptr->hostarrayptrs[I],hostptr->devarrayptrs[I],sizeof(T)*hostptr->N[I],cudaMemcpyDeviceToHost);
+            sprintf(dbgjnk,"%d dbg %d",I,hostptr->N[I]);
+            cuda_errortrap(dbgjnk);
+        }
+    } //for each array
+
+    //for the pull operation, I don't think any update of the device data structure is necessary
+
+    cudaSetDevice(olddev);
+
+    delete[] ldevarrayptrs;
+    delete[] devN;
+
+    return ret;
+}
+
+template<typename T> __host__ int curarray<T>::push()
+{
+    return curarray_push(this);
+}
+template<typename T> __host__ int curarray<T>::pull()
+{
+    return curarray_pull(this);
+}
+
+
+
+/*
+https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#memory-allocation-and-lifetime%5B/url%5D
+
+cudaMalloc() and cudaFree() have distinct semantics between the host and 
+device environments. When invoked from the host, cudaMalloc() allocates a 
+new region from unused device memory. When invoked from the device runtime 
+these functions map to device-side malloc() and free(). This implies that 
+within the device environment the total allocatable memory is limited to the 
+device malloc() heap size, which may be smaller than the available unused 
+device memory. Also, it is an error to invoke cudaFree() from the host 
+program on a pointer which was allocated by cudaMalloc() on the device 
+or vice-versa.
+
+So, basically this entire function is not going to work. I'll be unable to resize within 
+a kernel.
+*/
+
+/*
+
+template<typename T> __device__ int curarray<T>::dev_resizearray(int arraynum, int arraysize)
+{
+    int ret = 0;
+    T* newptr = NULL;
+    int I;
+    T def;
+
+    if(arraynum>=0 && arraynum<Narrays)
+    {
+        if(N[arraynum]!=arraysize)
+        {
+            if(arraysize<=0)
+            {
+                if(devarrayptrs[arraynum]!=NULL) cudaFree(devarrayptrs[arraynum]);
+                devarrayptrs[arraynum] = NULL;
+                N[arraynum] = 0;
+                ret = 1;
+                return ret;
+            }
+
+            cudaMalloc(&newptr,arraysize*sizeof(T));
+            if(newptr!=NULL)
+            {
+                //do I want to assume there is a copy operator? (operator=)
+                //for now, yes - write a more restrictive class later if I don't want it
+                if(devarrayptrs[arraynum]!=NULL)
+                {
+                    for(I=0;I<N[arraynum]&&I<arraysize;I++)
+                    {
+                        newptr[I] = devarrayptrs[arraynum][I];
+                    }
+                }
+                for(I=N[arraynum];I<arraysize;I++)
+                {
+                    newptr[I] = def;
+                }
+
+                if(devarrayptrs[arraynum]!=NULL) cudaFree(devarrayptrs[arraynum]);
+                devarrayptrs[arraynum] = newptr;
+                N[arraynum] = arraysize;
+                ret = 1;
+            }
+            else
+            {
+                ret = -1;
+            }
+        }
+        else
+        {
+            ret = 1;
+        }
+    }
+
+    return ret;
+}
+
+*/
+
+template<typename T> __host__ int curarray<T>::resizearray(int arraynum, int arraysize)
+{
+    int ret = 0;
+    T* newptr = NULL;
+    int I;
+    T def;
+
+    if(arraynum>=0 && arraynum<Narrays)
+    {
+        if(N[arraynum]!=arraysize)
+        {
+            if(arraysize<=0)
+            {
+                delete[] hostarrayptrs[arraynum];
+                hostarrayptrs[arraynum] = NULL;
+                N[arraynum] = 0;
+                ret = 1;
+                return ret;
+            }
+
+            newptr = new(std::nothrow) T[arraysize];
+            if(newptr!=NULL)
+            {
+                //do I want to assume there is a copy operator? (operator=)
+                //for now, yes - write a more restrictive class later if I don't want it
+                
+                if(hostarrayptrs[arraynum]!=NULL)
+                {
+                    for(I=0;I<N[arraynum]&&I<arraysize;I++)
+                    {
+                        newptr[I] = hostarrayptrs[arraynum][I];
+                    }
+                }
+                for(I=N[arraynum];I<arraysize;I++)
+                {
+                    newptr[I] = def;
+                }
+
+                //cudaFree(hostarrayptrs[arraynum]);
+                delete[] hostarrayptrs[arraynum];
+                hostarrayptrs[arraynum] = newptr;
+                N[arraynum] = arraysize;
+                ret = 1;
+            }
+            else
+            {
+                ret = -1;
+            }
+        }
+        else
+        {
+            ret = 1;
+        }
+    }
+
+    return ret;
+}
+
+};
+
+#endif

old/9apr26_prerefactor/include/amsculib2/cuvect2.hpp

@@ -0,0 +1,87 @@
+#ifndef __CUVECT2_HPP__
+#define __CUVECT2_HPP__
+
+namespace amscuda
+{
+
+    class cuvect2
+    {
+        public:
+        double x;
+        double y;
+
+
+        __host__ __device__ cuvect2();
+        __host__ __device__ ~cuvect2();
+        __host__ __device__ cuvect2(double _x, double _y);
+
+        __host__ __device__ double& operator[](const int I);
+        __host__ __device__ const double& operator[](const int I) const;
+
+        __host__ __device__ cuvect2 operator+(cuvect2 lhs);
+        __host__ __device__ cuvect2 operator-(cuvect2 lhs);
+        __host__ __device__ cuvect2 operator*(double lhs);
+        __host__ __device__ cuvect2 operator/(double lhs);
+        __host__ __device__ friend cuvect2 operator-(cuvect2 rhs);
+    };
+
+    class cumat2
+    {
+        public:
+        double dat[4];
+
+        __host__ __device__ cumat2();
+        __host__ __device__ ~cumat2();
+        __host__ __device__ double& operator[](const int I);
+        __host__ __device__ double& operator()(const int I, const int J);
+        __host__ __device__ double& at(const int I, const int J);
+
+        __host__ __device__ cumat2 operator+(cumat2 lhs);
+        __host__ __device__ cumat2 operator-(cumat2 lhs);
+        __host__ __device__ cumat2 operator*(double lhs);
+        __host__ __device__ cumat2 operator/(double lhs);
+        __host__ __device__ cuvect2 operator*(cuvect2 lhs);
+        __host__ __device__ cumat2 operator*(cumat2 lhs);
+        __host__ __device__ friend cumat2 operator-(cumat2 rhs);
+
+        __host__ __device__ double det();
+        __host__ __device__ cumat2 transpose();
+        __host__ __device__ cumat2 inverse();
+    };
+
+    __host__ __device__ double cuvect2_dot(cuvect2 a, cuvect2 b);
+    __host__ __device__ double cuvect2_cross(cuvect2 a, cuvect2 b);
+    __host__ __device__ double cuvect2_norm(cuvect2 a);
+    __host__ __device__ cuvect2 cuvect2_normalize(cuvect2 a);
+    __host__ __device__ cuvect2 cuvect2_proj(cuvect2 a, cuvect2 b);
+
+    //2x2 matrix operations
+    //matrix order is assumed to be mat[I,J] = mat[I+3*J]
+
+    //transpose a 2x2 matrix in place
+    __host__ __device__ void mat2_transpose(double *mat2inout);
+    
+    //copies src to dest
+    __host__ __device__ void mat2_copy(double *mat2_dest, const double *mat2_src);
+    
+    //inverts mat?inout[4]
+    __host__ __device__ void mat2_inverse(double *mat2inout);
+    
+    //rotatin matrix from angle
+    __host__ __device__ void mat2_rot_from_angle(double angle, double *mat2);
+
+    //multiplies c = a*b
+    __host__ __device__ void mat2_mult(double *mat2a, double *mat2b, double *mat2c);
+
+    // ret = a*b
+    __host__ __device__ cuvect2 mat2_mult(double *mat2a, cuvect2 b);
+
+
+
+    void test_cuvect2_1();
+    
+
+}; //end namespace amscuda
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/cuvect2f.hpp

@@ -0,0 +1,86 @@
+#ifndef __CUVECT2F_HPP__
+#define __CUVECT2F_HPP__
+
+namespace amscuda
+{
+
+    class cuvect2f
+    {
+        public:
+        float x;
+        float y;
+
+
+        __host__ __device__ cuvect2f();
+        __host__ __device__ ~cuvect2f();
+        __host__ __device__ cuvect2f(const float &_x, const float &_y);
+
+        __host__ __device__ float& operator[](const int &I);
+        __host__ __device__ const float& operator[](const int &I) const;
+
+        __host__ __device__ cuvect2f operator+(const cuvect2f &lhs);
+        __host__ __device__ cuvect2f operator-(const cuvect2f &lhs);
+        __host__ __device__ cuvect2f operator*(const float &lhs);
+        __host__ __device__ cuvect2f operator/(const float &lhs);
+        __host__ __device__ friend cuvect2f operator-(const cuvect2f &rhs);
+    };
+
+    class cumat2f
+    {
+        public:
+        float dat[4];
+
+        __host__ __device__ cumat2f();
+        __host__ __device__ ~cumat2f();
+        __host__ __device__ float& operator[](const int I);
+        __host__ __device__ float& operator()(const int I, const int J);
+        __host__ __device__ float& at(const int I, const int J);
+
+        __host__ __device__ cumat2f operator+(cumat2f lhs);
+        __host__ __device__ cumat2f operator-(cumat2f lhs);
+        __host__ __device__ cumat2f operator*(float lhs);
+        __host__ __device__ cumat2f operator/(float lhs);
+        __host__ __device__ cuvect2f operator*(cuvect2f lhs);
+        __host__ __device__ cumat2f operator*(cumat2f lhs);
+        __host__ __device__ friend cumat2f operator-(cumat2f rhs);
+
+        __host__ __device__ float det();
+        __host__ __device__ cumat2f transpose();
+        __host__ __device__ cumat2f inverse();
+    };
+
+    __host__ __device__ float cuvect2f_dot(const cuvect2f &a, const cuvect2f &b);
+    __host__ __device__ float cuvect2f_cross(const cuvect2f &a, const cuvect2f &b);
+    __host__ __device__ float cuvect2f_norm(const cuvect2f &a);
+    __host__ __device__ cuvect2f cuvect2f_normalize(const cuvect2f &a);
+    __host__ __device__ cuvect2f cuvect2f_proj(const cuvect2f &a, const cuvect2f &b);
+
+    //2x2 matrix operations
+    //matrix order is assumed to be mat[I,J] = mat[I+3*J]
+
+    //transpose a 2x2 matrix in place
+    __host__ __device__ void mat2f_transpose(float *mat2inout);
+    
+    //copies src to dest
+    __host__ __device__ void mat2f_copy(float *mat2f_dest, const float *mat2f_src);
+    
+    //inverts mat?inout[4]
+    __host__ __device__ void mat2f_inverse(float *mat2inout);
+    
+    //rotatin matrix from angle
+    __host__ __device__ void mat2f_rot_from_angle(float angle, float *mat2);
+
+    //multiplies c = a*b
+    __host__ __device__ void mat2f_mult(float *mat2a, float *mat2b, float *mat2c);
+
+    // ret = a*b
+    __host__ __device__ cuvect2f mat2f_mult(float *mat2a, const cuvect2f &b);
+
+
+    void test_cuvect2f_1();
+
+
+};
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/cuvect3.hpp

@@ -0,0 +1,88 @@
+#ifndef __CUVECT3_HPP__
+#define __CUVECT3_HPP__
+
+namespace amscuda
+{
+
+    class cuvect3
+    {
+        public:
+        double x;
+        double y;
+        double z;
+
+        __host__ __device__ cuvect3();
+        __host__ __device__ ~cuvect3();
+        __host__ __device__ cuvect3(double _x, double _y, double _z);
+        
+        
+        __host__ __device__ double& operator[](const int I);
+        __host__ __device__ const double& operator[](const int I) const;
+
+        __host__ __device__ cuvect3 operator+(cuvect3 lhs);
+        __host__ __device__ cuvect3 operator-(cuvect3 lhs);
+        __host__ __device__ cuvect3 operator*(double lhs);
+        __host__ __device__ cuvect3 operator/(double lhs);
+        __host__ __device__ friend cuvect3 operator-(cuvect3 rhs);
+    };
+
+    class cumat3
+    {
+        public:
+        double dat[9];
+
+        __host__ __device__ cumat3();
+        __host__ __device__ ~cumat3();
+        __host__ __device__ double& operator[](const int I);
+        __host__ __device__ double& operator()(const int I, const int J);
+        __host__ __device__ double& at(const int I, const int J);
+
+        __host__ __device__ cumat3 operator+(cumat3 lhs);
+        __host__ __device__ cumat3 operator-(cumat3 lhs);
+        __host__ __device__ cumat3 operator*(double lhs);
+        __host__ __device__ cumat3 operator/(double lhs);
+        __host__ __device__ cuvect3 operator*(cuvect3 lhs);
+        __host__ __device__ cumat3 operator*(cumat3 lhs);
+        __host__ __device__ friend cumat3 operator-(cumat3 rhs);
+
+        __host__ __device__ double det();
+        __host__ __device__ cumat3 transpose();
+        __host__ __device__ cumat3 inverse();
+    };
+
+    __host__ __device__ double cuvect3_dot(cuvect3 a, cuvect3 b);
+    __host__ __device__ cuvect3 cuvect3_cross(cuvect3 a, cuvect3 b);
+    __host__ __device__ double cuvect3_norm(cuvect3 a);
+    __host__ __device__ cuvect3 cuvect3_normalize(cuvect3 a);
+    __host__ __device__ cuvect3 cuvect3_proj(cuvect3 a, cuvect3 b);
+
+    //3x3 matrix operations
+    //matrix order is assumed to be mat[I,J] = mat[I+3*J]
+    
+    //transposes a 3x3 (9 element) matrix
+    __host__ __device__ void mat3_transpose(double *mat3inout);
+
+    //copies src to dest
+    __host__ __device__ void mat3_copy(double *mat3_dest, const double *mat3_src);
+
+    //returns determinant of 3x3 matrix
+    __host__ __device__ double mat3_det(double *mat3in);
+
+    //inverts a 3x3 (9 element) matrix
+    __host__ __device__ void mat3_inverse(double *mat3inout);
+
+    __host__ __device__ cuvect3 mat3_mult(double *mat3in, cuvect3 cvin);
+    __host__ __device__ void mat3_mult(double *matina, double *matinb, double *matout);
+
+    __host__ __device__ void mat3_hodgedual(cuvect3 vecin, double *matout);
+    __host__ __device__ void mat3_hodgedual(double *matin, cuvect3 vecout);
+
+    //returns direction cosine rotation matrix from axis and angle
+    __host__ __device__ void mat3_rot_from_axisangle(cuvect3 axis, double angle, double *matout);
+
+    __host__ void test_cudavect_logic1();
+
+}; //end namespace amscuda
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/cuvect3f.hpp

@@ -0,0 +1,122 @@
+#ifndef __CUVECT3F_HPP__
+#define __CUVECT3F_HPP__
+
+namespace amscuda
+{
+
+    class cuvect3f
+    {
+        public:
+        float x;
+        float y;
+        float z;
+
+        __host__ __device__ cuvect3f();
+        __host__ __device__ ~cuvect3f();
+        __host__ __device__ cuvect3f(const float &_x, const float &_y, const float &_z);
+        
+        
+        __host__ __device__ float& operator[](const int &I);
+        __host__ __device__ const float& operator[](const int &I) const;
+
+        __host__ __device__ cuvect3f operator+(const cuvect3f &rhs);
+        __host__ __device__ cuvect3f operator-(const cuvect3f &rhs);
+        __host__ __device__ cuvect3f operator*(const float &rhs);
+        __host__ __device__ cuvect3f operator/(const float &rhs);
+        __host__ __device__ friend cuvect3f operator-(const cuvect3f &rhs);
+
+        __host__ __device__ cuvect3f& operator+=(const cuvect3f &rhs);
+        __host__ __device__ cuvect3f& operator-=(const cuvect3f &rhs);
+        __host__ __device__ cuvect3f& operator/=(const float &rhs);
+        __host__ __device__ cuvect3f& operator*=(const float &rhs);
+
+    };
+
+    class cumat3f
+    {
+        public:
+        float m00,m10,m20; //named references to force register use?
+        float m01,m11,m21; //switched to column-major-order to match GLSL/lapack
+        float m02,m12,m22;
+
+        __host__ __device__ cumat3f();
+        __host__ __device__ ~cumat3f();
+        __host__ __device__ cumat3f(
+            const float & _m00, const float & _m01, const float & _m02,
+            const float & _m10, const float & _m11, const float & _m12,
+            const float & _m20, const float & _m21, const float & _m22
+        );
+
+        __host__ __device__ float& operator[](const int &I);
+        __host__ __device__ float& operator()(const int &I, const int &J);
+        __host__ __device__ float& at(const int &I, const int &J);
+
+        __host__ __device__ const float& operator[](const int &I) const;
+        __host__ __device__ const float& operator()(const int &I, const int &J) const;
+        __host__ __device__ const float& at(const int &I, const int &J) const;
+
+        __host__ __device__ cumat3f operator+(const cumat3f &rhs);
+        __host__ __device__ cumat3f operator-(const cumat3f &rhs);
+        __host__ __device__ cumat3f operator*(const float &rhs);
+        __host__ __device__ cumat3f operator/(const float &rhs);
+        __host__ __device__ cuvect3f operator*(const cuvect3f &rhs);
+        __host__ __device__ cumat3f operator*(const cumat3f &rhs);
+        __host__ __device__ friend cumat3f operator-(const cumat3f &rhs);
+        
+        __host__ __device__ float det();
+        __host__ __device__ cumat3f transpose();
+        __host__ __device__ cumat3f inverse();
+
+        __host__ __device__ float* data(); //pointer to float[9] representation of matrix
+        __host__ __device__ const float* data() const; //pointer to float[9] representation of matrix
+
+        //In place operations (to save GPU register use)
+        __host__ __device__ cumat3f& operator+=(const cumat3f &rhs);
+        __host__ __device__ cumat3f& operator-=(const cumat3f &rhs);
+        __host__ __device__ cumat3f& operator/=(const float &rhs);
+        __host__ __device__ cumat3f& operator*=(const float &rhs);
+        __host__ __device__ cumat3f& operator*=(const cumat3f &rhs);
+        
+    };
+
+    __host__ __device__ float cuvect3f_dot(const cuvect3f &a,const cuvect3f &b);
+    __host__ __device__ cuvect3f cuvect3f_cross(const cuvect3f &a, const cuvect3f &b);
+    __host__ __device__ float cuvect3f_norm(const cuvect3f &a);
+    __host__ __device__ cuvect3f cuvect3f_normalize(const cuvect3f &a);
+    __host__ __device__ cuvect3f cuvect3f_proj(const cuvect3f &a, const cuvect3f &b);
+
+    __host__ __device__ cumat3f hodge_dual(const cuvect3f &vin);
+    __host__ __device__ cuvect3f hodge_dual(const cumat3f &min);
+    __host__ __device__ cumat3f rotmat_from_axisangle(const cuvect3f &axis, const float &angle);
+
+
+    //3x3 matrix operations
+    //matrix order is assumed to be mat[I,J] = mat[I+3*J]
+    
+    //transposes a 3x3 (9 element) matrix
+    __host__ __device__ void mat3f_transpose(float *mat3inout);
+
+    //copies src to dest
+    __host__ __device__ void mat3f_copy(float *mat3f_dest, const float *mat3f_src);
+
+    //returns determinant of 3x3 matrix
+    __host__ __device__ float mat3f_det(float *mat3in);
+
+    //inverts a 3x3 (9 element) matrix
+    __host__ __device__ void mat3f_inverse(float *mat3inout);
+
+    __host__ __device__ cuvect3f mat3f_mult(float *mat3in, const cuvect3f &cvin);
+    __host__ __device__ void mat3f_mult(float *matina, float *matinb, float *matout);
+
+    __host__ __device__ void mat3f_hodgedual(const cuvect3f &vecin, float *matout);
+    __host__ __device__ void mat3f_hodgedual(float *matin, cuvect3f &vecout);
+
+    //returns direction cosine rotation matrix from axis and angle
+    __host__ __device__ void mat3f_rot_from_axisangle(cuvect3f axis, float angle, float *matout);
+
+    __host__ void test_cudavectf_logic1();
+
+};
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/cuvect4.hpp

@@ -0,0 +1,61 @@
+#ifndef __CUVECT4_HPP__
+#define __CUVECT4_HPP__
+
+namespace amscuda
+{
+
+    class cuvect4
+    {
+        public:
+        double x;
+        double y;
+        double z;
+        double w;
+
+        __host__ __device__ cuvect4();
+        __host__ __device__ ~cuvect4();
+        __host__ __device__ cuvect4(double _x, double _y, double _z, double _w);
+                
+        __host__ __device__ double& operator[](const int I);
+        __host__ __device__ const double& operator[](const int I) const;
+
+        __host__ __device__ cuvect4 operator+(cuvect4 lhs);
+        __host__ __device__ cuvect4 operator-(cuvect4 lhs);
+        __host__ __device__ cuvect4 operator*(double lhs);
+        __host__ __device__ cuvect4 operator/(double lhs);
+        __host__ __device__ friend cuvect4 operator-(cuvect4 rhs);
+    };
+
+    class cumat4
+    {
+        public:
+        double dat[16];
+
+        __host__ __device__ cumat4();
+        __host__ __device__ ~cumat4();
+        __host__ __device__ double& operator[](const int I);
+        __host__ __device__ double& operator()(const int I, const int J);
+        __host__ __device__ double& at(const int I, const int J);
+
+        __host__ __device__ cumat4 operator+(cumat4 lhs);
+        __host__ __device__ cumat4 operator-(cumat4 lhs);
+        __host__ __device__ cumat4 operator*(double lhs);
+        __host__ __device__ cumat4 operator/(double lhs);
+        __host__ __device__ cuvect4 operator*(cuvect4 lhs);
+        __host__ __device__ cumat4 operator*(cumat4 lhs);
+        __host__ __device__ friend cumat4 operator-(cumat4 rhs);
+
+        __host__ __device__ double det();
+        __host__ __device__ cumat4 transpose();
+        __host__ __device__ cumat4 inverse();
+    };
+
+    __host__ __device__ double cuvect4_dot(cuvect4 a, cuvect4 b);
+    __host__ __device__ double cuvect4_norm(cuvect4 a);
+    __host__ __device__ cuvect4 cuvect4_normalize(cuvect4 a);
+    __host__ __device__ cuvect4 cuvect4_proj(cuvect4 a, cuvect4 b);
+
+}; //end namespace amscuda
+
+#endif
+

old/9apr26_prerefactor/include/amsculib2/cuvect4f.hpp

@@ -0,0 +1,62 @@
+#ifndef __CUVECT4F_HPP__
+#define __CUVECT4F_HPP__
+
+namespace amscuda
+{
+
+    class cuvect4f
+    {
+        public:
+        float x;
+        float y;
+        float z;
+        float w;
+
+        __host__ __device__ cuvect4f();
+        __host__ __device__ ~cuvect4f();
+        __host__ __device__ cuvect4f(float _x, float _y, float _z, float _w);
+                
+        __host__ __device__ float& operator[](const int I);
+        __host__ __device__ const float& operator[](const int I) const;
+
+        __host__ __device__ cuvect4f operator+(cuvect4f lhs);
+        __host__ __device__ cuvect4f operator-(cuvect4f lhs);
+        __host__ __device__ cuvect4f operator*(float lhs);
+        __host__ __device__ cuvect4f operator/(float lhs);
+        __host__ __device__ friend cuvect4f operator-(cuvect4f rhs);
+    };
+
+    class cumat4f
+    {
+        public:
+        float dat[16];
+
+        __host__ __device__ cumat4f();
+        __host__ __device__ ~cumat4f();
+        __host__ __device__ float& operator[](const int I);
+        __host__ __device__ float& operator()(const int I, const int J);
+        __host__ __device__ float& at(const int I, const int J);
+
+        __host__ __device__ cumat4f operator+(cumat4f lhs);
+        __host__ __device__ cumat4f operator-(cumat4f lhs);
+        __host__ __device__ cumat4f operator*(float lhs);
+        __host__ __device__ cumat4f operator/(float lhs);
+        __host__ __device__ cuvect4f operator*(cuvect4f lhs);
+        __host__ __device__ cumat4f operator*(cumat4f lhs);
+        __host__ __device__ friend cumat4f operator-(cumat4f rhs);
+
+        __host__ __device__ float det();
+        __host__ __device__ cumat4f transpose();
+        __host__ __device__ cumat4f inverse();
+    };
+
+    __host__ __device__ float cuvect4f_dot(cuvect4f a, cuvect4f b);
+    __host__ __device__ float cuvect4f_norm(cuvect4f a);
+    __host__ __device__ cuvect4f cuvect4f_normalize(cuvect4f a);
+    __host__ __device__ cuvect4f cuvect4f_proj(cuvect4f a, cuvect4f b);
+
+
+};
+
+#endif
+

old/9apr26_prerefactor/src/amsculib2/amscu_comp128.cu

@@ -0,0 +1,476 @@
+#include <amsculib2/amsculib2.hpp>
+
+namespace amscuda
+{
+namespace cmp
+{
+
+__host__ __device__ cucomp128::cucomp128()
+{
+    real = 0.0;
+    imag = 0.0;
+    return;
+}
+
+__host__ __device__ cucomp128::~cucomp128()
+{
+    real = 0.0;
+    imag = 0.0;
+    return;
+}
+
+__host__ __device__ cucomp128::cucomp128(const cucomp128 &other)
+{
+    real = other.real;
+    imag = other.imag;
+
+    return;
+}
+
+__host__ __device__ cucomp128::cucomp128(const double &other)
+{
+    real = other;
+    imag = 0.0;
+    return;
+}
+
+__host__ __device__ cucomp128& cucomp128::operator=(cucomp128& other)
+{
+    real = other.real;
+    imag = other.imag;
+    return *this;
+}
+
+__host__ __device__ const cucomp128& cucomp128::operator=(const cucomp128& other) 
+{
+    this->real = other.real;
+    this->imag = other.imag;
+    return *this;
+}
+
+__host__ __device__ cucomp128& cucomp128::operator=(double& other)
+{
+    real = other;
+    imag = 0.0;
+    return *this;
+}
+
+__host__ __device__ const cucomp128& cucomp128::operator=(const double& other) 
+{
+    this->real = other;
+    this->imag = 0.0;
+    return *this;
+}
+
+__host__ __device__ double& cucomp128::operator[](int& ind)
+{
+    if(ind==0)
+    {
+        return this->real;
+    }
+    else
+    {
+        return this->imag;
+    }
+}
+
+__host__ __device__ const double& cucomp128::operator[](const int& ind) const
+{
+    if(ind==0)
+    {
+        return this->real;
+    }
+    else
+    {
+        return this->imag;
+    }
+}
+
+__host__ __device__ cucomp128 cucomp128::operator+(const cucomp128& z)
+{
+  cucomp128 ret;
+  ret.real = real + z.real;
+  ret.imag = imag + z.imag;
+  return ret;  
+}
+
+__host__ __device__ cucomp128 cucomp128::operator-(const cucomp128& z)
+{
+  cucomp128 ret;
+  ret.real = real - z.real;
+  ret.imag = imag - z.imag;
+  return ret;  
+}
+
+__host__ __device__ cucomp128 cucomp128::operator*(const cucomp128& z)
+{
+  cucomp128 ret;
+  ret.real = (real*z.real - imag*z.imag);
+  ret.imag = (imag*z.real + real*z.imag);
+  return ret;  
+}
+
+__host__ __device__ cucomp128 cucomp128::operator/(const cucomp128& z)
+{
+  cucomp128 ret;
+  double zm2 = z.real*z.real+z.imag*z.imag;
+  
+  if(zm2>0.0)
+  {
+      ret.real = (this->real*z.real+this->imag*z.imag)/zm2;
+      ret.imag = (this->imag*z.real-this->real*z.imag)/zm2;
+  }
+  else
+  {
+      ret.real = (double) finf;
+      ret.imag = (double) finf;
+  }
+
+  return ret;  
+}
+
+__host__ __device__ cucomp128 cucomp128::operator+(const double& z)
+{
+    cucomp128 ret;
+    ret.real = this->real + z;
+    ret.imag = this->imag;
+    return ret;
+}
+__host__ __device__ cucomp128 cucomp128::operator-(const double& z)
+{
+    cucomp128 ret;
+    ret.real = real-z;
+    ret.imag = imag;
+    return ret;
+}
+__host__ __device__ cucomp128 cucomp128::operator*(const double& z)
+{
+    cucomp128 ret;
+    ret.real = real*z;
+    ret.imag = imag*z;
+    return ret;
+}
+__host__ __device__ cucomp128 cucomp128::operator/(const double& z)
+{
+    cucomp128 ret;
+    if(z!=0.0f)
+    {
+        ret.real = real/z;
+        ret.imag = imag/z;
+    }
+    else
+    {
+        ret.real = finf;
+        ret.imag = finf;
+    }
+
+    return ret;
+}
+
+__host__ __device__ bool cucomp128::operator==(const cucomp128& z) const
+{
+    bool ret = 0;
+    if(z.real == real && z.imag == imag)
+    {
+        ret = 1;
+    }
+    return ret;
+}
+
+__host__ __device__ bool cucomp128::operator!=(const cucomp128& z) const
+{
+    return !(*this==z);
+}
+
+//sort first by real value, then by imaginary value
+//this is done so that an ordering exists, as long as two values aren't equal
+__host__ __device__ bool cucomp128::operator>(const cucomp128& z) const
+{
+    bool ret = 0;
+    if(this->real>z.real)
+    {
+        ret = 1;
+    }
+    else if(this->real==z.real)
+    {
+        if(this->imag>z.imag)
+        {
+            ret = 1;
+        }
+        else
+        {
+            ret = 0;
+        }
+    }
+    else
+    {
+        ret = 0;
+    }
+    return ret;
+}
+
+__host__ __device__ bool cucomp128::operator<(const cucomp128& z) const
+{
+    bool ret = 0;
+    if(this->real<z.real)
+    {
+        ret = 1;
+    }
+    else if(this->real==z.real)
+    {
+        if(this->imag<z.imag)
+        {
+            ret = 1;
+        }
+        else
+        {
+            ret = 0;
+        }
+    }
+    else
+    {
+        ret = 0;
+    }
+    return ret;
+}
+
+__host__ __device__ cucomp128 operator-(const cucomp128 &z)
+{
+	cucomp128 ret;
+	ret.real = -z.real;
+	ret.imag = -z.imag;
+	return ret;
+}
+
+__host__ __device__ bool cucomp128::operator>=(const cucomp128& z) const
+{
+    bool ret = (*this==z || *this>z);
+    return ret;
+}
+
+__host__ __device__ bool cucomp128::operator<=(const cucomp128& z) const
+{
+    bool ret = (*this==z || *this<z);
+    return ret;
+}
+
+__host__ __device__ bool cucomp128::isnan() const
+{
+  bool ret = 0;
+  if(::isnan(this->real) || ::isnan(this->imag))
+  {
+    ret = 1;
+  }
+  return ret;
+}
+
+__host__ __device__ bool cucomp128::isinf() const
+{
+  bool ret = 0;
+  //calls math.h isinf()
+  if(::isinf(this->real) || ::isinf(this->imag))
+  {
+    ret = 1;
+  }
+  return ret;
+}
+
+__host__ __device__ bool cucomp128::isreal() const
+{
+    bool ret = 1;
+    if(imag!=0.0f)
+    {
+        ret = 0;
+    }
+    return ret;
+}
+
+__host__ __device__ bool cucomp128::isimag() const
+{
+    bool ret = 1;
+    if(real!=0.0f)
+    {
+        ret = 0;
+    }
+    return ret;
+}
+
+__host__ __device__ bool cucomp128::iszero() const
+{
+    bool ret = 1;
+    if(real!=0.0f || imag!=0.0f)
+    {
+        ret = 0;
+    }
+    return ret;
+}
+
+__host__ __device__ double cucomp128::arg() const
+{
+    double ret = 0.0;
+    ret = (double) amscuda::arg((double)real,(double)imag);
+    return ret;
+}
+
+__host__ __device__ double cucomp128::mag() const
+{
+    double ret = 0.0;
+    ret = ::sqrt(real*real+imag*imag);
+    return ret;
+}
+
+__host__ __device__ cucomp128 cucomp128::conj() const
+{
+    cucomp128 ret;
+    ret.real = real;
+    ret.imag = -imag;
+    return ret;
+}
+
+__host__ __device__ double arg(cucomp128 z)
+{
+    return z.arg();
+}
+
+__host__ __device__ double abs(cucomp128 z)
+{
+    return z.mag();
+}
+
+__host__ __device__ cucomp128 dtocomp(double _r, double _i)
+{
+    cucomp128 ret;
+    ret.real = _r;
+    ret.imag = _i;
+    return ret;
+}
+
+__host__ __device__ double real(cucomp128 z)
+{
+    return z.real;
+}
+
+__host__ __device__ double imag(cucomp128 z)
+{
+    return z.imag;
+}
+
+__host__ __device__ cucomp128 sin(cucomp128 z)
+{
+    cucomp128 ret;
+    cucomp128 im1 = dtocomp(0.0f,1.0f);
+    cucomp128 div = dtocomp(0.0f,2.0f);
+
+    ret = (exp(im1*z)-exp(-im1*z))/div;
+
+    return ret;
+}
+
+__host__ __device__ cucomp128 cos(cucomp128 z)
+{
+    cucomp128 ret;
+    cucomp128 im1 = dtocomp(0.0f,1.0f);
+    cucomp128 div = dtocomp(2.0f,0.0f);
+
+    ret = (exp(im1*z)+exp(-im1*z))/div;
+
+    return ret;
+}
+
+__host__ __device__ cucomp128 tan(cucomp128 z)
+{
+    return sin(z)/cos(z);
+}
+
+__host__ __device__ cucomp128 exp(cucomp128 z)
+{
+    cucomp128 ret;
+    ret.real = ::exp(z.real)*::cos(z.imag);
+    ret.imag = ::exp(z.real)*::sin(z.imag);
+    return ret;
+}
+
+__host__ __device__ cucomp128 log(cucomp128 z)
+{
+    cucomp128 ret;
+    ret.real = ::log(::sqrt(z.real*z.real+z.imag*z.imag));
+    ret.imag = amscuda::arg(z.real,z.imag);
+    return ret;
+}
+
+__host__ __device__ cucomp128 conj(cucomp128 z)
+{
+    return z.conj();
+}
+
+__host__ __device__ cucomp128 cosh(cucomp128 z)
+{
+    cucomp128 ret;
+    cucomp128 div = dtocomp(2.0f,0.0f);
+
+    ret = (exp(z)+exp(-z))/div;
+
+    return ret;
+}
+
+__host__ __device__ cucomp128 sinh(cucomp128 z)
+{
+    cucomp128 ret;
+    cucomp128 div = dtocomp(2.0f,0.0f);
+
+    ret = (exp(z)-exp(-z))/div;
+
+    return ret;
+}
+
+__host__ __device__ cucomp128 tanh(cucomp128 z)
+{
+    return sinh(z)/cosh(z);
+}
+
+__host__ __device__ cucomp128 pow(cucomp128 z1, cucomp128 z2)
+{
+	cucomp128 ret;
+	if(z1.mag()>0.0)
+		ret = exp(z2*log(z1));
+	else
+		ret = dtocomp(0.0f,0.0f);
+	return ret;
+}
+
+
+void test_cucomp128_1()
+{
+    cucomp128 z1;
+    cucomp128 a,b,c;
+    double d1;
+    double f1;
+
+    printf("sizeof double=%ld\n",(long)(8*sizeof(f1)));
+    printf("sizeof double=%ld\n",(long)(8*sizeof(d1)));
+    printf("sizeof complex=%ld\n",(long)(8*sizeof(z1)));
+    printf("sizeof cucomp128=%ld\n",(long)(8*sizeof(a)));
+
+    a = dtocomp(1.0,1.0);
+    b = dtocomp(1.0,-1.0);
+
+    printf("a=%1.4f + %1.4fi\n",a[0],a[1]);
+    printf("b=%1.4f + %1.4fi\n",b[0],b[1]);
+    c = a+b;
+    printf("c=a+b: c=%1.4f + %1.4fi\n",c[0],c[1]);
+    c = a-b;
+    printf("c=a-b: c=%1.4f + %1.4fi\n",c[0],c[1]);
+    c = a*b;
+    printf("c=a*b: c=%1.4f + %1.4fi\n",c[0],c[1]);
+    c = a/b;
+    printf("c=a/b: c=%1.4f + %1.4fi\n",c[0],c[1]);
+    f1 = abs(a);
+    printf("abs(a)=%1.4f\n",f1);
+    f1 = arg(a);
+    printf("abs(a)=%1.4f pi\n",f1/pi);
+
+
+}
+
+}; //end namespace cmp
+}; //end namespace amscuda

old/9apr26_prerefactor/src/amsculib2/amscu_comp64.cu

@@ -0,0 +1,476 @@
+#include <amsculib2/amsculib2.hpp>
+
+namespace amscuda
+{
+namespace cmp
+{
+
+__host__ __device__ cucomp64::cucomp64()
+{
+    real = 0.0;
+    imag = 0.0;
+    return;
+}
+
+__host__ __device__ cucomp64::~cucomp64()
+{
+    real = 0.0;
+    imag = 0.0;
+    return;
+}
+
+__host__ __device__ cucomp64::cucomp64(const cucomp64 &other)
+{
+    real = other.real;
+    imag = other.imag;
+
+    return;
+}
+
+__host__ __device__ cucomp64::cucomp64(const float &other)
+{
+    real = other;
+    imag = 0.0;
+    return;
+}
+
+__host__ __device__ cucomp64& cucomp64::operator=(cucomp64& other)
+{
+    real = other.real;
+    imag = other.imag;
+    return *this;
+}
+
+__host__ __device__ const cucomp64& cucomp64::operator=(const cucomp64& other) 
+{
+    this->real = other.real;
+    this->imag = other.imag;
+    return *this;
+}
+
+__host__ __device__ cucomp64& cucomp64::operator=(float& other)
+{
+    real = other;
+    imag = 0.0;
+    return *this;
+}
+
+__host__ __device__ const cucomp64& cucomp64::operator=(const float& other) 
+{
+    this->real = other;
+    this->imag = 0.0;
+    return *this;
+}
+
+__host__ __device__ float& cucomp64::operator[](int& ind)
+{
+    if(ind==0)
+    {
+        return this->real;
+    }
+    else
+    {
+        return this->imag;
+    }
+}
+
+__host__ __device__ const float& cucomp64::operator[](const int& ind) const
+{
+    if(ind==0)
+    {
+        return this->real;
+    }
+    else
+    {
+        return this->imag;
+    }
+}
+
+__host__ __device__ cucomp64 cucomp64::operator+(const cucomp64& z)
+{
+  cucomp64 ret;
+  ret.real = real + z.real;
+  ret.imag = imag + z.imag;
+  return ret;  
+}
+
+__host__ __device__ cucomp64 cucomp64::operator-(const cucomp64& z)
+{
+  cucomp64 ret;
+  ret.real = real - z.real;
+  ret.imag = imag - z.imag;
+  return ret;  
+}
+
+__host__ __device__ cucomp64 cucomp64::operator*(const cucomp64& z)
+{
+  cucomp64 ret;
+  ret.real = (real*z.real - imag*z.imag);
+  ret.imag = (imag*z.real + real*z.imag);
+  return ret;  
+}
+
+__host__ __device__ cucomp64 cucomp64::operator/(const cucomp64& z)
+{
+  cucomp64 ret;
+  float zm2 = z.real*z.real+z.imag*z.imag;
+  
+  if(zm2>0.0)
+  {
+      ret.real = (this->real*z.real+this->imag*z.imag)/zm2;
+      ret.imag = (this->imag*z.real-this->real*z.imag)/zm2;
+  }
+  else
+  {
+      ret.real = (float) finf;
+      ret.imag = (float) finf;
+  }
+
+  return ret;  
+}
+
+__host__ __device__ cucomp64 cucomp64::operator+(const float& z)
+{
+    cucomp64 ret;
+    ret.real = this->real + z;
+    ret.imag = this->imag;
+    return ret;
+}
+__host__ __device__ cucomp64 cucomp64::operator-(const float& z)
+{
+    cucomp64 ret;
+    ret.real = real-z;
+    ret.imag = imag;
+    return ret;
+}
+__host__ __device__ cucomp64 cucomp64::operator*(const float& z)
+{
+    cucomp64 ret;
+    ret.real = real*z;
+    ret.imag = imag*z;
+    return ret;
+}
+__host__ __device__ cucomp64 cucomp64::operator/(const float& z)
+{
+    cucomp64 ret;
+    if(z!=0.0f)
+    {
+        ret.real = real/z;
+        ret.imag = imag/z;
+    }
+    else
+    {
+        ret.real = finf;
+        ret.imag = finf;
+    }
+
+    return ret;
+}
+
+__host__ __device__ bool cucomp64::operator==(const cucomp64& z) const
+{
+    bool ret = 0;
+    if(z.real == real && z.imag == imag)
+    {
+        ret = 1;
+    }
+    return ret;
+}
+
+__host__ __device__ bool cucomp64::operator!=(const cucomp64& z) const
+{
+    return !(*this==z);
+}
+
+//sort first by real value, then by imaginary value
+//this is done so that an ordering exists, as long as two values aren't equal
+__host__ __device__ bool cucomp64::operator>(const cucomp64& z) const
+{
+    bool ret = 0;
+    if(this->real>z.real)
+    {
+        ret = 1;
+    }
+    else if(this->real==z.real)
+    {
+        if(this->imag>z.imag)
+        {
+            ret = 1;
+        }
+        else
+        {
+            ret = 0;
+        }
+    }
+    else
+    {
+        ret = 0;
+    }
+    return ret;
+}
+
+__host__ __device__ bool cucomp64::operator<(const cucomp64& z) const
+{
+    bool ret = 0;
+    if(this->real<z.real)
+    {
+        ret = 1;
+    }
+    else if(this->real==z.real)
+    {
+        if(this->imag<z.imag)
+        {
+            ret = 1;
+        }
+        else
+        {
+            ret = 0;
+        }
+    }
+    else
+    {
+        ret = 0;
+    }
+    return ret;
+}
+
+__host__ __device__ cucomp64 operator-(const cucomp64 &z)
+{
+	cucomp64 ret;
+	ret.real = -z.real;
+	ret.imag = -z.imag;
+	return ret;
+}
+
+__host__ __device__ bool cucomp64::operator>=(const cucomp64& z) const
+{
+    bool ret = (*this==z || *this>z);
+    return ret;
+}
+
+__host__ __device__ bool cucomp64::operator<=(const cucomp64& z) const
+{
+    bool ret = (*this==z || *this<z);
+    return ret;
+}
+
+__host__ __device__ bool cucomp64::isnan() const
+{
+  bool ret = 0;
+  if(::isnan(this->real) || ::isnan(this->imag))
+  {
+    ret = 1;
+  }
+  return ret;
+}
+
+__host__ __device__ bool cucomp64::isinf() const
+{
+  bool ret = 0;
+  //calls math.h isinf()
+  if(::isinf(this->real) || ::isinf(this->imag))
+  {
+    ret = 1;
+  }
+  return ret;
+}
+
+__host__ __device__ bool cucomp64::isreal() const
+{
+    bool ret = 1;
+    if(imag!=0.0f)
+    {
+        ret = 0;
+    }
+    return ret;
+}
+
+__host__ __device__ bool cucomp64::isimag() const
+{
+    bool ret = 1;
+    if(real!=0.0f)
+    {
+        ret = 0;
+    }
+    return ret;
+}
+
+__host__ __device__ bool cucomp64::iszero() const
+{
+    bool ret = 1;
+    if(real!=0.0f || imag!=0.0f)
+    {
+        ret = 0;
+    }
+    return ret;
+}
+
+__host__ __device__ float cucomp64::arg() const
+{
+    float ret = 0.0;
+    ret = (float) amscuda::arg((double)real,(double)imag);
+    return ret;
+}
+
+__host__ __device__ float cucomp64::mag() const
+{
+    float ret = 0.0;
+    ret = ::sqrt(real*real+imag*imag);
+    return ret;
+}
+
+__host__ __device__ cucomp64 cucomp64::conj() const
+{
+    cucomp64 ret;
+    ret.real = real;
+    ret.imag = -imag;
+    return ret;
+}
+
+__host__ __device__ float arg(cucomp64 z)
+{
+    return z.arg();
+}
+
+__host__ __device__ float abs(cucomp64 z)
+{
+    return z.mag();
+}
+
+__host__ __device__ cucomp64 dtocomp64(float _r, float _i)
+{
+    cucomp64 ret;
+    ret.real = _r;
+    ret.imag = _i;
+    return ret;
+}
+
+__host__ __device__ float real(cucomp64 z)
+{
+    return z.real;
+}
+
+__host__ __device__ float imag(cucomp64 z)
+{
+    return z.imag;
+}
+
+__host__ __device__ cucomp64 sin(cucomp64 z)
+{
+    cucomp64 ret;
+    cucomp64 im1 = dtocomp64(0.0f,1.0f);
+    cucomp64 div = dtocomp64(0.0f,2.0f);
+
+    ret = (exp(im1*z)-exp(-im1*z))/div;
+
+    return ret;
+}
+
+__host__ __device__ cucomp64 cos(cucomp64 z)
+{
+    cucomp64 ret;
+    cucomp64 im1 = dtocomp64(0.0f,1.0f);
+    cucomp64 div = dtocomp64(2.0f,0.0f);
+
+    ret = (exp(im1*z)+exp(-im1*z))/div;
+
+    return ret;
+}
+
+__host__ __device__ cucomp64 tan(cucomp64 z)
+{
+    return sin(z)/cos(z);
+}
+
+__host__ __device__ cucomp64 exp(cucomp64 z)
+{
+    cucomp64 ret;
+    ret.real = ::exp(z.real)*::cos(z.imag);
+    ret.imag = ::exp(z.real)*::sin(z.imag);
+    return ret;
+}
+
+__host__ __device__ cucomp64 log(cucomp64 z)
+{
+    cucomp64 ret;
+    ret.real = ::log(::sqrt(z.real*z.real+z.imag*z.imag));
+    ret.imag = amscuda::arg(z.real,z.imag);
+    return ret;
+}
+
+__host__ __device__ cucomp64 conj(cucomp64 z)
+{
+    return z.conj();
+}
+
+__host__ __device__ cucomp64 cosh(cucomp64 z)
+{
+    cucomp64 ret;
+    cucomp64 div = dtocomp64(2.0f,0.0f);
+
+    ret = (exp(z)+exp(-z))/div;
+
+    return ret;
+}
+
+__host__ __device__ cucomp64 sinh(cucomp64 z)
+{
+    cucomp64 ret;
+    cucomp64 div = dtocomp64(2.0f,0.0f);
+
+    ret = (exp(z)-exp(-z))/div;
+
+    return ret;
+}
+
+__host__ __device__ cucomp64 tanh(cucomp64 z)
+{
+    return sinh(z)/cosh(z);
+}
+
+__host__ __device__ cucomp64 pow(cucomp64 z1, cucomp64 z2)
+{
+	cucomp64 ret;
+	if(z1.mag()>0.0)
+		ret = exp(z2*log(z1));
+	else
+		ret = dtocomp64(0.0f,0.0f);
+	return ret;
+}
+
+
+void test_cucomp64_1()
+{
+    cucomp64 z1;
+    cucomp64 a,b,c;
+    double d1;
+    float f1;
+
+    printf("sizeof double=%ld\n",(long)(8*sizeof(f1)));
+    printf("sizeof double=%ld\n",(long)(8*sizeof(d1)));
+    printf("sizeof complex=%ld\n",(long)(8*sizeof(z1)));
+    printf("sizeof cucomp128=%ld\n",(long)(8*sizeof(a)));
+
+    a = dtocomp64(1.0,1.0);
+    b = dtocomp64(1.0,-1.0);
+
+    printf("a=%1.4f + %1.4fi\n",a[0],a[1]);
+    printf("b=%1.4f + %1.4fi\n",b[0],b[1]);
+    c = a+b;
+    printf("c=a+b: c=%1.4f + %1.4fi\n",c[0],c[1]);
+    c = a-b;
+    printf("c=a-b: c=%1.4f + %1.4fi\n",c[0],c[1]);
+    c = a*b;
+    printf("c=a*b: c=%1.4f + %1.4fi\n",c[0],c[1]);
+    c = a/b;
+    printf("c=a/b: c=%1.4f + %1.4fi\n",c[0],c[1]);
+    f1 = abs(a);
+    printf("abs(a)=%1.4f\n",f1);
+    f1 = arg(a);
+    printf("abs(a)=%1.4f pi\n",f1/pi);
+
+
+}
+
+}; //end namespace cmp
+}; //end namespace amscuda

old/9apr26_prerefactor/src/amsculib2/amscu_cudafunctions.cu

@@ -0,0 +1,21 @@
+#include <amsculib2/amsculib2.hpp>
+
+namespace amscuda
+{
+
+int cuda_errortrap(const char *msgheader)
+{
+    int ret = 0;
+    cudaError_t err = cudaSuccess;
+
+    err = cudaGetLastError();
+    if(err!=cudaSuccess)
+    {
+        printf("%s :%s\n",msgheader,cudaGetErrorString(err));
+        ret = 1;
+    }
+
+    return ret;
+}
+
+}; //end namespace amscuda