Markus Mottl avatar Markus Mottl committed a012a57

Added packed/unpacked

Comments (0)

Files changed (5)

+2008-10-26:  Added two new matrix functions:
+
+               * packed
+               * unpacked
+
 2008-10-15:  Added new LAPACK function:
 
                * geqrf
-release-4-6-0
+release-4-6-1
 name="lacaml"
-version="4.6.0"
+version="4.6.1"
 description="LACAML - BLAS/LAPACK-interface for OCaml"
 
 requires="lacaml.core"
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *)
 
+open Printf
 open Bigarray
 open Numberxx
 open Utils
   done;
   res
 
-let detri ?(up = true) ?(ar = 1) ?(ac = 1) ?n a =
+let detri ?(up = true) ?n ?(ar = 1) ?(ac = 1) a =
   let loc = "Lacaml.Impl.NPREC.Mat.detri" in
   let n = get_n_of_square "a" loc ar ac a n in
   if up then
       done
     done
 
+let packed ?(up = true) ?n ?(ar = 1) ?(ac = 1) a =
+  let loc = "Lacaml.Impl.NPREC.Mat.packed" in
+  let n = get_n_of_square "a" loc ar ac a n in
+  let dst = Array1.create prec fortran_layout n in
+  let pos_ref = ref 1 in
+  if up then
+    for c = 1 to n do
+      for r = 1 to c do
+        dst.{!pos_ref} <- a.{r, c};
+        incr pos_ref;
+      done
+    done
+  else
+    for c = 1 to n do
+      for r = c to n do
+        dst.{!pos_ref} <- a.{r, c};
+        incr pos_ref;
+      done
+    done;
+  dst
+
+let unpacked ?(up = true) src =
+  let n_vec = Array1.dim src in
+  let n = truncate (sqrt (float (8 * n_vec + 1)) /. 2.) in
+  if n * n + n <> n_vec then
+    let loc = "Lacaml.Impl.NPREC.Mat.unpacked" in
+    failwith (sprintf "%s: illegal vector length: %d" loc n_vec)
+  else
+    let a = make0 n n in
+    let pos_ref = ref 1 in
+    if up then
+      for c = 1 to n do
+        for r = 1 to c do
+          a.{r, c} <- src.{!pos_ref};
+          incr pos_ref;
+        done
+      done
+    else
+      for c = 1 to n do
+        for r = c to n do
+          a.{r, c} <- src.{!pos_ref};
+          incr pos_ref;
+        done
+      done;
+    a
+
 external direct_scal_mat :
   int -> (* M *)
   int -> (* N *)
 (** {6 Matrix transformations} *)
 
 val transpose : ?m : int -> ?n : int -> ?ar : int -> ?ac : int -> mat -> mat
-(** [transpose a] @return the transpose of matrix [a]. *)
+(** [transpose a] @return the transpose of (sub-)matrix [a].
 
-val detri : ?up : bool -> ?ar : int -> ?ac : int -> ?n : int -> mat -> unit
-(** [detri ?up ?ar ?ac ?n a] takes a triangular (sub-)matrix [a], i.e. one
+    @param m default = [Mat.dim1 a]
+    @param n default = [Mat.dim2 a]
+    @param ar default = [1]
+    @param ac default = [1]
+*)
+
+val detri : ?up : bool -> ?n : int -> ?ar : int -> ?ac : int -> mat -> unit
+(** [detri ?up ?n ?ar ?ac a] takes a triangular (sub-)matrix [a], i.e. one
     where only the upper (iff [up] is true) or lower triangle is defined,
     and makes it a symmetric matrix by mirroring the defined triangle
     along the diagonal.
 
     @param up default = [true]
+    @param n default = [Mat.dim1 a]
     @param ar default = [1]
     @param ac default = [1]
-    @param n default = [Mat.dim1 a]
+*)
+
+val packed : ?up : bool -> ?n : int -> ?ar : int -> ?ac : int -> mat -> vec
+(** [packed ?up ?n ?ar ?ac a] @return (sub-)matrix [a] in packed
+    storage format.
+
+    @param up default = [true]
+    @param n default = [Mat.dim2 a]
+    @param ar default = [1]
+    @param ac default = [1]
+*)
+
+val unpacked : ?up : bool -> vec -> mat
+(** [packed ?up x] @return an upper or lower (depending on [up])
+    triangular matrix from packed representation [vec].  The other
+    triangle of the matrix will be filled with zeros.
+
+    @param up default = [true]
 *)
 
 
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.