Markus Mottl avatar Markus Mottl committed 8d004ac

Added Mat.transpose_copy

Comments (0)

Files changed (6)

              Some matrix and vector types were polymorphic, thus
              leading to much worse element access speed.
 
+	     Improved compiler selection in Makefile.
+
+	     Improved some matrix operations.
+
+	     Added new function:
+
+               * Mat.transpose_copy
+
 2009-06-13:  Added new function:
 
                * Mat.as_vec
 name="lacaml"
-version="5.4.2"
+version="5.4.3"
 description="LACAML - BLAS/LAPACK-interface for OCaml"
 
 requires="lacaml.core"
 install:	libinstall
 uninstall:	libuninstall
 
+CC = $(shell ocamlc -config | grep native_c_compiler | cut -d' ' -f2)
+
 # Generation rules for precision-dependent C-code
 %2_S_c.o:	%_SD_c.c
 	$(CC) -c $(CFLAGS) $(CINCFLAGS) -I'$(OCAMLLIBPATH)' $< -o $@
   let gen = genarray_of_array2 mat in
   reshape_1 gen (dim1 mat * dim2 mat)
 
-let transpose ?m ?n ?(ar = 1) ?(ac = 1) (a : mat) =
+external direct_transpose_copy :
+  m : int ->
+  n : int ->
+  ar : int ->
+  ac : int ->
+  a : mat ->
+  br : int ->
+  bc : int ->
+  b : mat ->
+  unit = "lacaml_NPRECtranspose_copy_stub_bc" "lacaml_NPRECtranspose_copy_stub"
+
+let transpose_copy ?m ?n ?(ar = 1) ?(ac = 1) a ?(br = 1) ?(bc = 1) b =
+  let loc = "Lacaml.Impl.NPREC.Mat.transpose_copy" in
+  let m = get_dim1_mat loc a_str a ar m_str m in
+  let n = get_dim2_mat loc a_str a ac n_str n in
+  check_dim_mat loc b_str br bc b n m;
+  direct_transpose_copy ~m ~n ~ar ~ac ~a ~br ~bc ~b
+
+let transpose ?m ?n ?(ar = 1) ?(ac = 1) a =
   let loc = "Lacaml.Impl.NPREC.Mat.transpose" in
   let m = get_dim1_mat loc a_str a ar m_str m in
   let n = get_dim2_mat loc a_str a ac n_str n in
-  let res = create n m in
-  let ar_1 = ar - 1 in
-  let ac_1 = ac - 1 in
-  for c = 1 to n do
-    let ac_col_1 = ac_1 + c in
-    for r = 1 to m do res.{c, r} <- a.{ar_1 + r, ac_col_1} done
-  done;
-  res
+  let b = create n m in
+  direct_transpose_copy ~m ~n ~ar ~ac ~a ~br:1 ~bc:1 ~b;
+  b
 
 let detri ?(up = true) ?n ?(ar = 1) ?(ac = 1) (a : mat) =
   let loc = "Lacaml.Impl.NPREC.Mat.detri" in
     for c = 1 to n - 1 do
       let ar_c = ar + c in
       let ac_c = ac + c in
-      for r = 0 to c - 1 do
-        a.{ar_c, ac + r} <- a.{ar + r, ac_c}
-      done
+      for r = 0 to c - 1 do a.{ar_c, ac + r} <- a.{ar + r, ac_c} done
     done
   else
     for c = 1 to n - 1 do
       let ar_c = ar + c in
       let ac_c = ac + c in
-      for r = 0 to c - 1 do
-        a.{ar + r, ac_c} <- a.{ar_c, ac + r}
-      done
+      for r = 0 to c - 1 do a.{ar + r, ac_c} <- a.{ar_c, ac + r} done
     done
 
 let packed ?(up = true) ?n ?(ar = 1) ?(ac = 1) (a : mat) =
 
 (** {6 Matrix transformations} *)
 
+val transpose_copy :
+  ?m : int -> ?n : int ->
+  ?ar : int -> ?ac : int -> mat ->
+  ?br : int -> ?bc : int -> mat ->
+  unit
+(** [transpose_copy ?m ?n ?ar ?ac a ?br ?bc b] copy the transpose
+    of (sub-)matrix [a] into (sub-)matrix [b].
+
+    @param m default = [Mat.dim1 a]
+    @param n default = [Mat.dim2 a]
+    @param ar default = [1]
+    @param ac default = [1]
+    @param br default = [1]
+    @param bc default = [1]
+*)
+
+
 val transpose : ?m : int -> ?n : int -> ?ar : int -> ?ac : int -> mat -> mat
-(** [transpose a] @return the transpose of (sub-)matrix [a].
+(** [transpose ?m ?n ?ar ?ac aa] @return the transpose of (sub-)matrix [a].
 
     @param m default = [Mat.dim1 a]
     @param n default = [Mat.dim2 a]
 static NUMBER number_minus_one = NUMBER_MINUS_ONE;
 
 
+/* transpose_copy */
+
+extern void FUN(copy)(
+  integer *N,
+  NUMBER *X, integer *INCX,
+  NUMBER *Y, integer *INCY);
+
+CAMLprim value LFUN(transpose_copy_stub)(
+  value vM, value vN,
+  value vAR, value vAC, value vA,
+  value vBR, value vBC, value vB)
+{
+  CAMLparam2(vA, vB);
+  integer GET_INT(M), GET_INT(N);
+
+  if (M > 0 && N > 0) {
+    MAT_PARAMS(A);
+    MAT_PARAMS(B);
+    NUMBER *A_last = A_data + rows_A * N;
+    caml_enter_blocking_section();
+      do {
+        FUN(copy)(&M, A_data, &integer_one, B_data, &rows_B);
+        A_data += rows_A;
+        B_data++;
+      } while (A_data != A_last);
+    caml_leave_blocking_section();
+  }
+
+  CAMLreturn(Val_unit);
+}
+
+CAMLprim value LFUN(transpose_copy_stub_bc)(value *argv, int argn)
+{
+  return LFUN(transpose_copy_stub)(
+    argv[0], argv[1], argv[2], argv[3], argv[4], argv[5], argv[6], argv[7]);
+}
+
+
 /* scal */
 
 extern void FUN(scal)(
   value vAR, value vAC, value vA)
 {
   CAMLparam1(vA);
+  integer GET_INT(M), GET_INT(N);
 
-  integer GET_INT(M), GET_INT(N);
-  CREATE_NUMBER(ALPHA);
-
-  MAT_PARAMS(A);
-
-  INIT_NUMBER(ALPHA);
-
-  caml_enter_blocking_section();
-    if (rows_A == M) {
-      integer MN = M * N;
-      FUN(scal)(&MN, &ALPHA, A_data, &integer_one);
-    } else {
-      NUMBER *A_src = A_data + rows_A * (N - 1);
-      while (A_src >= A_data) {
-        FUN(scal)(&M, &ALPHA, A_src, &integer_one);
-        A_src -= rows_A;
+  if ( M > 0 && N > 0) {
+    CREATE_NUMBER(ALPHA);
+    MAT_PARAMS(A);
+    INIT_NUMBER(ALPHA);
+    caml_enter_blocking_section();
+      if (rows_A == M) {
+        integer MN = M * N;
+        FUN(scal)(&MN, &ALPHA, A_data, &integer_one);
+      } else {
+        NUMBER *A_last = A_data + rows_A * N;
+        do {
+          FUN(scal)(&M, &ALPHA, A_data, &integer_one);
+          A_data += rows_A;
+        } while (A_data != A_last);
       }
-    }
-  caml_leave_blocking_section();
+    caml_leave_blocking_section();
+  }
 
   CAMLreturn(Val_unit);
 }
   value vOFSALPHAs, value vALPHAs)
 {
   CAMLparam2(vALPHAs, vA);
-
   integer GET_INT(M), GET_INT(N);
 
-  VEC_PARAMS(ALPHAs);
-  MAT_PARAMS(A);
-
-  NUMBER *A_src = A_data + rows_A * (N - 1);
-  NUMBER *ALPHAs_src = ALPHAs_data + (N - 1);
-
-  caml_enter_blocking_section();
-    while (A_src >= A_data) {
-      FUN(scal)(&M, ALPHAs_src, A_src, &integer_one);
-      A_src -= rows_A;
-      ALPHAs_src--;
-    }
-  caml_leave_blocking_section();
+  if (M > 0 && N > 0) {
+    VEC_PARAMS(ALPHAs);
+    MAT_PARAMS(A);
+    NUMBER *A_last = A_data + rows_A * N;
+    caml_enter_blocking_section();
+      do {
+        FUN(scal)(&M, ALPHAs_data, A_data, &integer_one);
+        A_data += rows_A;
+        ALPHAs_data++;
+      } while (A_data != A_last);
+    caml_leave_blocking_section();
+  }
 
   CAMLreturn(Val_unit);
 }
   value vYR, value vYC, value vY)
 {
   CAMLparam2(vX, vY);
-
   integer GET_INT(M), GET_INT(N);
-  CREATE_NUMBER(ALPHA);
-
-  MAT_PARAMS(X);
-  MAT_PARAMS(Y);
-
-  INIT_NUMBER(ALPHA);
-
-  caml_enter_blocking_section();
-    if (rows_X == M && rows_Y == M) {
-      integer MN = M * N;
-      FUN(axpy)(&MN, &ALPHA, X_data, &integer_one, Y_data, &integer_one);
-    } else {
-      NUMBER *X_src = X_data + rows_X * (N - 1);
-      NUMBER *Y_dst = Y_data + rows_Y * (N - 1);
-      while (X_src >= X_data) {
-        FUN(axpy)(&M, &ALPHA, X_src, &integer_one, Y_dst, &integer_one);
-        X_src -= rows_X;
-        Y_dst -= rows_Y;
+  if (M > 0 && N > 0) {
+    CREATE_NUMBER(ALPHA);
+    MAT_PARAMS(X);
+    MAT_PARAMS(Y);
+    INIT_NUMBER(ALPHA);
+    caml_enter_blocking_section();
+      if (rows_X == M && rows_Y == M) {
+        integer MN = M * N;
+        FUN(axpy)(&MN, &ALPHA, X_data, &integer_one, Y_data, &integer_one);
+      } else {
+        NUMBER *X_last = X_data + rows_X * N;
+        do {
+          FUN(axpy)(&M, &ALPHA, X_data, &integer_one, Y_data, &integer_one);
+          X_data += rows_X;
+          Y_data += rows_Y;
+        } while (X_data != X_last);
       }
-    }
-  caml_leave_blocking_section();
+    caml_leave_blocking_section();
+  }
 
   CAMLreturn(Val_unit);
 }
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.