Anonymous avatar Anonymous committed 83cbb21

add more functions + more tests

Comments (0)

Files changed (3)

 // simplistic wrapper around the C-interface to the BLAS
 package cblas
 
-//#include "cblas.h"
-//#include <stdlib.h>
-//#include <string.h>
+/*
+ #include "cblas.h"
+*/
 import "C"
 
-func init() {
-	
-}
+/*
+ * Enumerated and derived types
+ */
+
+type Index int
+
+type Order C.enum_CBLAS_ORDER
+const (
+	RowMajor = Order(101)
+	ColMajor = Order(102)
+)
+
+type Transpose C.enum_CBLAS_TRANSPOSE
+const (
+	NoTrans   = Transpose(111)
+	Trans     = Transpose(112)
+	ConjTrans = Transpose(113)
+)
+
+type UpLo C.enum_CBLAS_UPLO
+const (
+	Upper = UpLo(121)
+	Lower = UpLo(122)
+)
+
+type Diag C.enum_CBLAS_DIAG
+const (
+	NonUnit = Diag(131)
+	Unit    = Diag(132)
+	)
+
+type Side C.enum_CBLAS_SIDE
+const (
+	Left = Side(141)
+	Right= Side(142)
+)
 
 // EOF

pkg/cblas_lvl1.go

 // simplistic wrapper around the C-interface to the BLAS
 package cblas
 
-//#include "cblas.h"
-//#include <stdlib.h>
-//#include <string.h>
+/*
+ #include <complex.h>
+ #undef I
+ #include "cblas.h"
+*/
 import "C"
 import "unsafe"
 
-func Sdsdot(N int, alpha float, 
-	X []float, incX int,
-	Y []float, incY int) float {
-	c_X := (*C.float)(unsafe.Pointer(&X[0]))
-	c_Y := (*C.float)(unsafe.Pointer(&Y[0]))
+/*
+ * ===========================================================================
+ * Prototypes for level 1 BLAS functions (complex are recast as routines)
+ * ===========================================================================
+ */
+
+/*
+ float  cblas_sdsdot(const int N, const float alpha, 
+ const float *X, const int incX, 
+ const float *Y, const int incY);
+ */
+func Sdsdot(alpha float, x, y []float) float {
+	if len(x) != len(y) {
+		panic("slices' size differ")
+	}
+
+	c_N := C.int(len(x))
+	c_alpha := C.float(alpha)
+
+	c_X    := (*C.float)(unsafe.Pointer(&x[0]))
+	c_incX := C.int(1)
+
+	c_Y    := (*C.float)(unsafe.Pointer(&y[0]))
+	c_incY := C.int(-1)
 
 	return float(
-		C.cblas_sdsdot(C.int(N), C.float(alpha),
-		 c_X, C.int(incX),
-		 c_Y, C.int(incY)))
+		C.cblas_sdsdot(c_N, c_alpha, 
+		c_X, c_incX,
+		c_Y, c_incY))
 }
 
+/*
+ double cblas_dsdot(const int N, const float *X, const int incX, const float *Y,
+ const int incY);
+ */
+func Dsdot(x,y []float) float64 {
+	if len(x) != len(y) {
+		panic("slices' size differ")
+	}
+
+	c_N := C.int(len(x))
+	
+	c_X    := (*C.float)(unsafe.Pointer(&x[0]))
+	c_incX := C.int(1)
+
+	c_Y    := (*C.float)(unsafe.Pointer(&y[0]))
+	c_incY := C.int(-1)
+
+	return float64(
+		C.cblas_dsdot(c_N,
+		c_X, c_incX,
+		c_Y, c_incY))
+}
+
+/*
+ float  cblas_sdot(const int N, const float  *X, const int incX,
+ const float  *Y, const int incY);
+ */
+func Sdot(x,y []float) float {
+	if len(x) != len(y) {
+		panic("slices' size differ")
+	}
+
+	c_N := C.int(len(x))
+	
+	c_X    := (*C.float)(unsafe.Pointer(&x[0]))
+	c_incX := C.int(1)
+
+	c_Y    := (*C.float)(unsafe.Pointer(&y[0]))
+	c_incY := C.int(-1)
+
+	return float(
+		C.cblas_sdot(c_N,
+		c_X, c_incX,
+		c_Y, c_incY))
+	
+}
+
+/*
+ double cblas_ddot(const int N, const double *X, const int incX,
+ const double *Y, const int incY);
+ */
+func Ddot(x,y []float64) float64 {
+	if len(x) != len(y) {
+		panic("slices' size differ")
+	}
+
+	c_N := C.int(len(x))
+	
+	c_X    := (*C.double)(unsafe.Pointer(&x[0]))
+	c_incX := C.int(1)
+
+	c_Y    := (*C.double)(unsafe.Pointer(&y[0]))
+	c_incY := C.int(-1)
+
+	return float64(
+		C.cblas_ddot(c_N,
+		c_X, c_incX,
+		c_Y, c_incY))
+}
+
+/*
+ void   cblas_cdotu_sub(const int N, const void *X, const int incX,
+ const void *Y, const int incY, void *dotu);
+ */
+func Cdotu(x,y []complex64) complex64 {
+	
+	if len(x) != len(y) {
+		panic("slices' size differ")
+	}
+
+	var ret C.complex64
+	c_ret := unsafe.Pointer(&ret)
+
+	c_N := C.int(len(x))
+	
+	c_X    := unsafe.Pointer(&x[0])
+	c_incX := C.int(1)
+
+	c_Y    := unsafe.Pointer(&y[0])
+	c_incY := C.int(-1)
+
+	C.cblas_cdotu_sub(c_N,
+		c_X, c_incX,
+		c_Y, c_incY,
+		c_ret)
+
+	return cmplx(
+		float32(C.crealf(ret)), 
+		float32(C.cimagf(ret)))
+}
+
+/*
+ void   cblas_cdotc_sub(const int N, const void *X, const int incX,
+ const void *Y, const int incY, void *dotc);
+ */
+func Cdotc(x,y []complex64) complex64 {
+	
+	if len(x) != len(y) {
+		panic("slices' size differ")
+	}
+
+	var ret C.complex64
+	c_ret := unsafe.Pointer(&ret)
+
+	c_N := C.int(len(x))
+	
+	c_X    := unsafe.Pointer(&x[0])
+	c_incX := C.int(1)
+
+	c_Y    := unsafe.Pointer(&y[0])
+	c_incY := C.int(-1)
+
+	C.cblas_cdotc_sub(c_N,
+		c_X, c_incX,
+		c_Y, c_incY,
+		c_ret)
+
+	return cmplx(
+		float32(C.crealf(ret)), 
+		float32(C.cimagf(ret)))
+}
+
+/*
+ void   cblas_zdotu_sub(const int N, const void *X, const int incX,
+                        const void *Y, const int incY, void *dotu);
+ */
+func Zdotu(x,y []complex128) complex128 {
+	
+	if len(x) != len(y) {
+		panic("slices' size differ")
+	}
+
+	var ret C.complex128
+
+	c_ret := unsafe.Pointer(&ret)
+
+	c_N := C.int(len(x))
+	
+	c_X    := unsafe.Pointer(&x[0])
+	c_incX := C.int(1)
+
+	c_Y    := unsafe.Pointer(&y[0])
+	c_incY := C.int(-1)
+
+	C.cblas_zdotu_sub(c_N,
+		c_X, c_incX,
+		c_Y, c_incY,
+		c_ret)
+
+	return cmplx(
+		float64(C.creal(ret)), 
+		float64(C.cimag(ret)))
+}
+
+/*
+ void   cblas_zdotc_sub(const int N, const void *X, const int incX,
+ const void *Y, const int incY, void *dotc);
+ func Zdotc(x,y []complex64) complex64 {
+ 
+ if len(x) != len(y) {
+ panic("slices' size differ")
+ }
+
+ var ret C.complex64_t
+ c_ret := (*C.char)(unsafe.Pointer(ret))
+
+ c_N := C.int(len(x))
+ 
+ c_X    := (*C.complex64_t)(unsafe.Pointer(&x[0]))
+ c_incX := C.int(1)
+
+ c_Y    := (*C.complex64_t)(unsafe.Pointer(&y[0]))
+ c_incY := C.int(-1)
+
+ C.cblas_cdotu_sub(c_N,
+ c_X, c_incX,
+ c_Y, c_incY,
+ c_ret)
+
+ return cmplx(
+ float64(ret.real), 
+ float64(ret.imag))
+ }
+ */
+
+/*
+ float  cblas_snrm2(const int N, const float *X, const int incX);
+ */
+func Snrm2(x []float) float {
+	
+	c_N := C.int(len(x))
+	
+	c_X    := (*C.float)(unsafe.Pointer(&x[0]))
+	c_incX := C.int(1)
+
+	return float(C.cblas_snrm2(c_N, c_X, c_incX))
+}
+
+/*
+ float  cblas_sasum(const int N, const float *X, const int incX);
+ */
+func Sasum(x []float) float {
+	
+	c_N := C.int(len(x))
+	
+	c_X    := (*C.float)(unsafe.Pointer(&x[0]))
+	c_incX := C.int(1)
+
+	return float(C.cblas_sasum(c_N, c_X, c_incX))
+}
+
+/*
+ double cblas_dnrm2(const int N, const double *X, const int incX);
+ */
+func Dnrm2(x []float64) float64 {
+	c_N := C.int(len(x))
+	
+	c_X    := (*C.double)(unsafe.Pointer(&x[0]))
+	c_incX := C.int(1)
+
+	return float64(C.cblas_dnrm2(c_N, c_X, c_incX))
+}
+
+/*
+ double cblas_dasum(const int N, const double *X, const int incX);
+ */
+func Dasum(x []float64) float64 {
+	c_N := C.int(len(x))
+	
+	c_X    := (*C.double)(unsafe.Pointer(&x[0]))
+	c_incX := C.int(1)
+
+	return float64(C.cblas_dasum(c_N, c_X, c_incX))
+}
+
+/*
+ float  cblas_scnrm2(const int N, const void *X, const int incX);
+ func Snrm2(x []complex) float {
+ 
+ c_N := C.int(len(x))
+ 
+ c_X    := (*C.float)(unsafe.Pointer(&x[0]))
+ c_incX := C.int(1)
+
+ return float(C.cblas_scnrm2(c_N, c_X, c_incX))
+ }
+ */
+
+/*
+ float  cblas_scasum(const int N, const void *X, const int incX);
+ */
+func Scasum(x []complex) float {
+	
+	c_N := C.int(len(x))
+	
+	c_X    := unsafe.Pointer(&x[0])
+	c_incX := C.int(1)
+
+	return float(C.cblas_scasum(c_N, c_X, c_incX))
+}
+
+/*
+ double  cblas_dznrm2(const int N, const void *X, const int incX);
+ func Dznrm2(x []complex64) float64 {
+ 
+ c_N := C.int(len(x))
+ 
+ c_X    := (*C.double)(unsafe.Pointer(&x[0]))
+ c_incX := C.int(1)
+
+ return float64(C.cblas_dznrm2(c_N, c_X, c_incX))
+ }
+ */
+
+/*
+ double  cblas_dzasum(const int N, const void *X, const int incX);
+ */
+func Dzasum(x []complex128) float64 {
+	
+	c_N := C.int(len(x))
+	
+	c_X    := unsafe.Pointer(&x[0])
+	c_incX := C.int(1)
+
+	return float64(C.cblas_dzasum(c_N, c_X, c_incX))
+}
+
+
 // functions having standard 4 prefixes (S D C Z)
 
-func Isamax(N int, X []float, incX int) int {
-	c_X := (*C.float)(unsafe.Pointer(&X[0]))
-	return int(C.cblas_isamax(C.int(N), c_X, C.int(incX)))
+/*
+ CBLAS_INDEX cblas_isamax(const int N, const float  *X, const int incX);
+ */
+func Isamax(x []float) Index {
+	c_N := C.int(len(x))
+	c_X := (*C.float)(unsafe.Pointer(&x[0]))
+	c_incX := C.int(1)
+	return Index(C.cblas_isamax(c_N, c_X, c_incX))
 }
+
+/*
+ CBLAS_INDEX cblas_idamax(const int N, const double *X, const int incX);
+ */
+func Idamax(x []float64) Index {
+	c_N := C.int(len(x))
+	c_X := (*C.double)(unsafe.Pointer(&x[0]))
+	c_incX := C.int(1)
+	return Index(C.cblas_idamax(c_N, c_X, c_incX))
+}
+
 package main
 
-import _ "cblas"
+import "fmt"
+import "cblas"
 
 func main() {
-	
+	{
+		x := []float{1., 2., 3.}
+		fmt.Printf("x=%v\n", x)
+		sumx := cblas.Sasum(x)
+		fmt.Printf("sasum(x)=%v\n", sumx)
+	}
+
+	{
+		x := []float64{1., 2., 3.}
+		fmt.Printf("x=%v\n", x)
+		sumx := cblas.Dasum(x)
+		fmt.Printf("dasum(x)=%v\n", sumx)
+	}
+
+	{
+		x := []complex{ cmplx(1., 1.),
+			cmplx(2., 2.),
+			cmplx(3., 3.)}
+		fmt.Printf("x=%v\n", x)
+		sumx := cblas.Scasum(x)
+		fmt.Printf("scasum(x)=%v\n", sumx)
+	}
+
+	{
+		x := []complex128{ cmplx(1., 1.),
+			cmplx(2., 2.),
+			cmplx(3., 3.)}
+		fmt.Printf("x=%v\n", x)
+		sumx := cblas.Dzasum(x)
+		fmt.Printf("dzasum(x)=%v\n", sumx)
+	}
+
+	{
+		x := []float{ 0.733 }
+		y := []float{ 0.825 }
+		alpha := float(0.)
+
+		exp := float(0.604725)
+		val := cblas.Sdsdot(alpha, x,y)
+		fmt.Printf("sdsdot: exp: %v\n", exp)
+		fmt.Printf("sdsdot: val: %v\n", val)
+	}
+	{
+		x := []complex128{ cmplx(-0.87, -0.631) }
+		y := []complex128{ cmplx(-0.7,  -0.224) }
+
+		exp := complex128(cmplx(0.467656, 0.63658))
+		val := cblas.Zdotu(x,y)
+
+		fmt.Printf("zdotu: exp: %v\n", exp)
+		fmt.Printf("zdotu: val: %v\n", val)
+	}
+	return
 }
 
 /*
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.