gpuRcuda: CUDA GPU functions for R Objects

# Currently separate methods until integer methods
# written in CUDA, then will consolidate
# This separate allows for using the parent gpuR methods
# Note - when consolidating will need to make sure to import typeof from gpuR




#' Single Precision CUDA Matrix Addition/Subtraction
#' @param e1 An fcudaMatrix object
#' @param e2 An fcudaMatrix object
#' @export
setMethod("Arith", c(e1="fcudaMatrix", e2="fcudaMatrix"),
					function(e1, e2)
					{						
						Z <- cudaMatrix(nrow=nrow(e1), ncol=ncol(e2), type="float")
						if(length(e2[]) != length(e1[])) stop("Lengths of matrices must match")
						Z@address <- e1@address
						
						op = .Generic[[1]]
						
						switch(op,
									 `+` = cpp_vienna_cudaMatrix_saxpy(1, 
									 																	Z@address, 
									 																	e2@address),
									 `-` = cpp_vienna_cudaMatrix_saxpy(-1, 
									 																	Z@address, 
									 																	e2@address),
									{
										stop("undefined operation")	
									}
						)
						
						return(Z)
					},
valueClass = "fcudaMatrix"
)

#' Double Precision CUDA Matrix Addition/Subtraction
#' @param e1 An dcudaMatrix object
#' @param e2 An dcudaMatrix object
#' @export
setMethod("Arith", c(e1="dcudaMatrix", e2="dcudaMatrix"),
					function(e1, e2)
					{
						if(!deviceHasDouble()){
							stop("Selected GPU does not support double precision")
						}else{
							op = .Generic[[1]]
							
							Z <- cudaMatrix(nrow=nrow(e1), ncol=ncol(e2), type="double")
							
							if(length(e2[]) != length(e1[])) stop("Lengths of matrices must match")
							Z@address <- e1@address
							
							op = .Generic[[1]]
							
							switch(op,
										 `+` = cpp_vienna_cudaMatrix_daxpy(1, 
										 																	Z@address, 
										 																	e2@address),
										 `-` = cpp_vienna_cudaMatrix_daxpy(-1, 
										 																	Z@address, 
										 																	e2@address),
										{
											stop("undefined operation")	
										}
							)
							return(Z)
						}
					},
valueClass = "dcudaMatrix"
)


#' CUDA Matrix Multiplication
#' @param x An fcudaMatrix object
#' @param y An fcudaMatrix object
#' @export
setMethod("%*%", signature(x="cudaMatrix", y = "cudaMatrix"),
					function(x,y)
					{
						if( dim(x)[2] != dim(y)[1]){
							stop("Non-conformant matrices")
						}
						return(gpu_Mat_mult(x, y))
					},
					valueClass = "cudaMatrix"
)

#' @title The Number of Rows/Columns of a gpuRcudaMatrix
#' @description \code{nrow} and \code{ncol} return the number of rows
#' or columns present in \code{x}.
#' @param x A gpuRcudaMatrix object
#' @return An integer of length 1
#' @docType methods
#' @rdname nrow-gpuRcudaMatrix
#' @aliases nrow,cudaMatrix
#' @aliases nrow,nvMatrix
#' @aliases ncol,cudaMatrix
#' @aliases ncol,nvMatrix
#' @author Charles Determan Jr.
#' @export
setMethod('nrow', signature(x="gpuRcudaMatrix"), 
					function(x) {
						switch(typeof(x),
									 "integer" = return(cpp_nrow(x@address, 4L)),
									 "float" = return(cpp_nrow(x@address, 6L)),
									 "double" = return(cpp_nrow(x@address, 8L))
						)
					}
)

#' @rdname nrow-gpuRcudaMatrix
#' @export
setMethod('ncol', signature(x="gpuRcudaMatrix"),
					function(x) {
						switch(typeof(x),
									 "integer" = return(cpp_ncol(x@address, 4L)),
									 "float" = return(cpp_ncol(x@address, 6L)),
									 "double" = return(cpp_ncol(x@address, 8L))
						)
					}
)


#' @title gpuRcuda dim method
#' @param x A gpuRcuda matrix object
#' @return A length 2 vector of the number of rows and columns respectively.
#' @docType methods
#' @rdname dim-methods
#' @author Charles Determan Jr.
#' @aliases dim,cudaMatrix
#' @aliases dim,nvMatrix
#' @export
setMethod('dim', signature(x="gpuRcudaMatrix"),
					function(x) return(c(nrow(x), ncol(x))))

#' @title Extract gpuRcuda elements
#' @param x A gpuRcuda object
#' @param i missing
#' @param j missing
#' @param drop missing
#' @docType methods
#' @rdname extract-methods
#' @author Charles Determan Jr.
#' @export
setMethod("[",
					signature(x = "cudaMatrix", i = "missing", j = "missing", drop = "missing"),
					function(x, i, j, drop) {
						
						init <- ifelse(typeof(x) == "double", 0, 0L)
						out <- matrix(init, nrow = nrow(x), ncol = ncol(x))
						
						switch(typeof(x),
									 "integer" = {
									 		cudaMatToSEXP(x@address, out, 4L)
									 		return(out)
									 },
									 "float" = {
									 		cudaMatToSEXP(x@address, out, 6L)
									 		return(out)
									 },
									 "double" = {
									 		cudaMatToSEXP(x@address, out, 8L)
									 		return(out)
									 }
						)
					})

#' @rdname extract-methods
#' @export
setMethod("[",
					signature(x = "nvMatrix", i = "missing", j = "missing", drop = "missing"),
					function(x, i, j, drop) {
						
						init <- ifelse(typeof(x) == "double", 0, 0L)
						out <- matrix(init, nrow = nrow(x), ncol = ncol(x))
						
						switch(typeof(x),
									 "integer" = {
									 		nvMatToSEXP(x@address, out, 4L)
										 	return(out)
									 	},
									 "float" = {
									 		nvMatToSEXP(x@address, out, 6L)
										 	return(out)
									 	},
									 "double" = {
									 		nvMatToSEXP(x@address, out, 8L)
									 		return(out)
									 }
						)
					})

#' Single Precision CUDA Matrix Multiplication
#' @param x An fcudaMatrix object
#' @param y An fcudaMatrix object
#' @export
setMethod("%*%", signature(x="fcudaMatrix", y = "fcudaMatrix"),
					function(x,y)
					{
						if( dim(x)[2] != dim(y)[1]){
							stop("Non-conformant matrices")
						}
						
						C <- cudaMatrix(nrow=nrow(x), ncol=ncol(y), type="float")
						
						cpp_vienna_cudaMatrix_sgemm(x@address,
																				y@address,
																				C@address)
						
						return(C)
					},
					valueClass = "fcudaMatrix"
)

#' Double Precision CUDA Matrix Multiplication
#' @param x An cudaMatrix object
#' @param y An cudaMatrix object
#' @export
setMethod("%*%", signature(x="cudaMatrix", y = "cudaMatrix"),
					function(x,y)
					{
						if( dim(x)[2] != dim(y)[1]){
							stop("Non-conformant matrices")
						}
						
						if(!deviceHasDouble()){
							stop("Selected GPU does not support double precision")
						}else{
							C <- cudaMatrix(nrow=nrow(x), ncol=ncol(y), type="double")
							
							cpp_vienna_cudaMatrix_dgemm(x@address,
																					y@address,
																					C@address)
						}
						
						return(C)
					},
					valueClass = "dcudaMatrix"
)
          
# #' Single Precision CUDA Matrix Multiplication
# #' @param x An fcudaMatrix object
# #' @param y An fcudaMatrix object
# #' @export
 setMethod("%*%", signature(x="fcudaMatrix", y = "fcudaMatrix"),
 					function(x,y)
 					{
 						if( dim(x)[2] != dim(y)[1]){
 							stop("Non-conformant matrices")
 						}
 						
 						C <- cudaMatrix(nrow=nrow(x), ncol=ncol(y), type="float")
 						
 						cpp_vienna_cudaMatrix_sgemm(x@address,
 																				y@address,
 																				C@address)
 						
 						return(C)
 					},
 					valueClass = "fcudaMatrix"
 )
 
# #' Double Precision CUDA Matrix Multiplication
# #' @param x An fcudaMatrix object
# #' @param y An fcudaMatrix object
# #' @export
 setMethod("%*%", signature(x="cudaMatrix", y = "cudaMatrix"),
 					function(x,y)
 					{
 						if( dim(x)[2] != dim(y)[1]){
 							stop("Non-conformant matrices")
 						}
 						
 						if(!deviceHasDouble()){
 							stop("Selected GPU does not support double precision")
 						}else{
 							C <- cudaMatrix(nrow=nrow(x), ncol=ncol(y), type="double")
 							
 							cpp_vienna_cudaMatrix_dgemm(x@address,
 																					y@address,
 																					C@address)
 						}
 						
 						return(C)
 					},
 					valueClass = "dcudaMatrix"
 )