Markus Mottl avatar Markus Mottl committed 47ad11d

Added orgqr + QR example

Comments (0)

Files changed (11)

 examples/lin_eq_comp/lin_eq_comp$
 examples/lin_reg/lin_reg$
 examples/svd/svd$
+examples/qr/qr$
+2009-09-16:  Added new functions:
+
+               * sygv
+
+             Thanks to Christophe Troestler <Christophe.Troestler@umh.ac.be>
+             for contributing sygv!
+
+               * orgqr
+
+             Added example for QR-factorization.
+
 2009-09-09:  Fixed a few performance bugs.
 
              Some matrix and vector types were polymorphic, thus
              leading to much worse element access speed.
 
-	     Improved compiler selection in Makefile.
+             Improved compiler selection in Makefile.
 
-	     Improved some matrix operations.
+             Improved some matrix operations.
 
-	     Added new function:
+             Added new function:
 
                * Mat.transpose_copy
 
 -include Makefile.conf
 
-EXAMPLES = $(filter-out examples/OMakefile examples/CVS examples/Makefile.examples, $(wildcard examples/*))
+EXAMPLES = $(filter-out examples/OMakefile examples/Makefile.examples, $(wildcard examples/*))
 
 .PHONY: all
 all:

examples/qr/Makefile

+SOURCES = qr.ml
+RESULT = qr
+
+-include ../Makefile.examples

examples/qr/OMakefile

+OCamlMakeProjDefaults(qr.exe)
+
+InstantiateOCamlEnv()

examples/qr/qr.ml

+(* File: qr.ml
+
+   Copyright (C) 2009-
+
+     Markus Mottl
+     email: markus.mottl@gmail.com
+     WWW: http://www.ocaml.info
+
+   This library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2 of the License, or (at your option) any later version.
+
+   This library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with this library; if not, write to the Free Software
+   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+*)
+
+open Format
+open Bigarray
+
+open Lacaml.Impl.D
+open Lacaml.Io
+
+let () = Random.self_init ()
+
+let () =
+  let m = 10 in
+  let n = 5 in
+  let a = Mat.random m n in
+  printf "@[<2>A =@\n@\n@[%a@]@]@\n@\n" pp_mat a;
+  let tau = geqrf a in
+  let r = Mat.make0 n n in
+  let r = lacpy ~m:n ~uplo:`U a ~b:r in
+  orgqr ~tau a;
+  printf "@[<2>Q =@\n@\n@[%a@]@]@\n@\n" pp_mat a;
+  printf "@[<2>R =@\n@\n@[%a@]@]@\n@\n" pp_mat r;
+  printf "@[<2>QR = A =@\n@\n@[%a@]@]@\n" pp_mat (gemm a r)
 name="lacaml"
-version="5.4.3"
+version="5.4.4"
 description="LACAML - BLAS/LAPACK-interface for OCaml"
 
 requires="lacaml.core"
   let n = get_n_of_a loc ar ac a n in
   let uplo = get_uplo_char up in
   let min_work = lansy_min_lwork n norm in
-  let work, _lwork = get_work loc Vec.create work min_work min_work "lwork" in
+  let work, _lwork = get_work loc Vec.create work min_work min_work lwork_str in
   let norm = get_norm_char norm in
   direct_lansy ~norm ~uplo ~n ~ar ~ac ~a ~work
 
 
 (* Linear equations (computational routines) *)
 
+(* ORGQR *)
+
+external direct_orgqr :
+  m : int ->
+  n : int ->
+  k : int ->
+  work : vec ->
+  lwork : int ->
+  tau : vec ->
+  ar : int ->
+  ac : int ->
+  a : mat ->
+  int = "lacaml_FPRECorgqr_stub_bc" "lacaml_FPRECorgqr_stub"
+
+let orgqr_min_lwork ~n = max 1 n
+
+let orgqr_get_opt_lwork loc ~m ~n ~k ~tau ~ar ~ac ~a =
+  let work = Vec.create 1 in
+  let info = direct_orgqr ~m ~n ~k ~work ~lwork:~-1 ~tau ~ar ~ac ~a in
+  if info = 0 then int_of_float work.{1}
+  else orgqr_err ~loc ~m ~n ~k ~work ~a ~err:info
+
+let orgqr_opt_lwork ?m ?n ?k ~tau ?(ar = 1) ?(ac = 1) a =
+  let loc = "Lacaml.Impl.FPREC.orgqr_opt_lwork" in
+  let m, n, k = orgqr_get_params loc ?m ?n ?k ~tau ?ar ?ac a in
+  orgqr_get_opt_lwork loc ~m ~n ~k ~tau ~ar ~ac ~a
+
+let orgqr ?m ?n ?k ?work ~tau ?(ar = 1) ?(ac = 1) a =
+  let loc = "Lacaml.Impl.FPREC.orgqr" in
+  let m, n, k = orgqr_get_params loc ?m ?n ?k ~tau ?ar ?ac a in
+  let work, lwork =
+    let min_lwork = orgqr_min_lwork ~n in
+    get_work loc Vec.create work min_lwork min_lwork lwork_str
+  in
+  let info = direct_orgqr ~m ~n ~k ~work ~lwork ~tau ~ar ~ac ~a in
+  if info = 0 then ()
+  else orgqr_err ~loc ~m ~n ~k ~work ~a ~err:info
+
+
 (* GECON *)
 
 external direct_gecon :
   let loc = "Lacaml.Impl.FPREC.gecon" in
   let n = get_n_of_a loc ar ac a n in
   let work, _lwork =
-    get_work
-      loc Vec.create work (gecon_min_lwork n) (gecon_min_lwork n) "lwork" in
+    let min_lwork = gecon_min_lwork n in
+    get_work loc Vec.create work min_lwork min_lwork lwork_str
+  in
   let iwork, _liwork =
     get_work
       loc Common.create_int_vec iwork
-      (gecon_min_liwork n) (gecon_min_liwork n) "liwork" in
+      (gecon_min_liwork n) (gecon_min_liwork n) liwork_str in
   let anorm =
     match anorm with
     | None -> lange ~norm:(norm :> norm4) ~m:n ~n a
   let uplo = get_uplo_char up in
   let work, _lwork =
     get_work
-      loc Vec.create work (sycon_min_lwork n) (sycon_min_lwork n) "lwork" in
+      loc Vec.create work (sycon_min_lwork n) (sycon_min_lwork n) lwork_str in
   let iwork, _liwork =
     get_work
       loc Common.create_int_vec iwork
-      (sycon_min_liwork n) (sycon_min_liwork n) "liwork" in
+      (sycon_min_liwork n) (sycon_min_liwork n) liwork_str in
   let ipiv =
     if ipiv = None then sytrf ~n ~up ~work ~ar ~ac a
     else sytrf_get_ipiv loc ipiv n in
   let uplo = get_uplo_char up in
   let work, _lwork =
     get_work
-      loc Vec.create work (pocon_min_lwork n) (pocon_min_lwork n) "lwork" in
+      loc Vec.create work (pocon_min_lwork n) (pocon_min_lwork n) lwork_str in
   let iwork, _liwork =
     get_work
       loc Common.create_int_vec iwork
-      (pocon_min_liwork n) (pocon_min_liwork n) "liwork" in
+      (pocon_min_liwork n) (pocon_min_liwork n) liwork_str in
   let anorm =
     match anorm with
     | None -> lange ~m:n ~n ~work ~ar ~ac a
 
 (** Linear equations (computational routines) *)
 
+(* ORGQR *)
+
+val orgqr_min_lwork : n : int -> int
+(** [orgqr_min_lwork ~n] @return the minimum length of the
+    work-array used by the [orgqr]-function if the matrix has [n]
+    columns. *)
+
+val orgqr_opt_lwork :
+  ?m : int ->
+  ?n : int ->
+  ?k : int ->
+  tau : vec ->
+  ?ar : int ->
+  ?ac : int ->
+  mat ->
+  int
+(** [orgqr_opt_lwork ?m ?n ?k ~tau ?ar ?ac a] @return the optimum
+    length of the work-array used by the [orgqr_opt_lwork]-function
+    given matrix [a], optionally its logical dimensions [m] and
+    [n], and the number of reflectors [k].
+
+    @param m default = available number of rows in matrix [a]
+    @param n default = available number of columns in matrix [a]
+    @param k default = available number of elements in vector [tau]
+*)
+
+val orgqr :
+  ?m : int ->
+  ?n : int ->
+  ?k : int ->
+  ?work : vec ->
+  tau : vec ->
+  ?ar : int ->
+  ?ac : int ->
+  mat ->
+  unit
+(** [orgqr ?m ?n ?k ?work ~tau ?ar ?ac a] see LAPACK documentation!
+
+    @param m default = available number of rows in matrix [a]
+    @param n default = available number of columns in matrix [a]
+    @param k default = available number of elements in vector [tau]
+*)
+
+
 (* GECON *)
 
 val gecon_min_lwork : int -> int
 }
 
 
+/** ORGQR */
+
+extern void FUN(orgqr)(
+  integer *M,
+  integer *N,
+  integer *K,
+  REAL *A, integer *LDA,
+  REAL *TAU,
+  REAL *WORK,
+  integer *LWORK,
+  integer *INFO);
+
+CAMLprim value LFUN(orgqr_stub)(
+  value vM,
+  value vN,
+  value vK,
+  value vWORK,
+  value vLWORK,
+  value vTAU,
+  value vAR,
+  value vAC,
+  value vA)
+{
+  CAMLparam2(vTAU, vA);
+
+  integer GET_INT(M), GET_INT(N), GET_INT(K), GET_INT(LWORK), INFO;
+
+  VEC_PARAMS1(WORK);
+  VEC_PARAMS1(TAU);
+  MAT_PARAMS(A);
+
+  caml_enter_blocking_section();  /* Allow other threads */
+  FUN(orgqr)(
+    &M, &N, &K,
+    A_data, &rows_A,
+    TAU_data,
+    WORK_data, &LWORK, &INFO);
+  caml_leave_blocking_section(); /* Disallow other threads */
+
+  CAMLreturn(Val_int(INFO));
+}
+
+CAMLprim value LFUN(orgqr_stub_bc)(value *argv, int argn)
+{
+  return
+    LFUN(orgqr_stub)(
+      argv[0], argv[1], argv[2], argv[3], argv[4],
+      argv[5], argv[6], argv[7], argv[8]);
+}
+
+
 /** GECON */
 
 extern void FUN(gecon)(
 let ipiv_str = "ipiv"
 let iseed_str = "iseed"
 let k_str = "k"
+let lwork_str = "lwork"
+let liwork_str = "liwork"
 let kd_str = "kd"
 let kl_str = "kl"
 let ku_str = "ku"
 let n_str = "n"
 let nrhs_str = "nrhs"
 let s_str = "s"
+let tau_str = "tau"
 let u_str = "u"
 let um_str = "um"
 let un_str = "un"
 
 (* ??MV auxiliary functions *)
 
-(* [get_dim_vec loc ofsvec incvec vec n_name n] if the dimension [n]
-   is given, check that the vector [vec] is big enough, otherwise
-   return the maximal [n] for the given vector [vec]. *)
+(* [get_dim_vec loc vec_name ofsvec incvec vec n_name n] if the
+   dimension [n] is given, check that the vector [vec] is big enough,
+   otherwise return the maximal [n] for the given vector [vec]. *)
 let get_dim_vec loc vec_name ofsvec incvec vec n_name = function
   | Some n ->
       if n < 0 then invalid_arg (sprintf "%s: %s < 0" loc n_name);
   nrhs
 
 
+(* ORGQR - Auxiliary Functions *)
+
+let orgqr_err ~loc ~m ~n ~k ~work ~a ~err =
+  let msg =
+    match err with
+    | -1 -> sprintf "m: valid=[0..[ got=%d" m
+    | -2 -> sprintf "n: valid=[0..%d] got=%d" m n
+    | -3 -> sprintf "k: valid=[0..%d] got=%d" n k
+    | -5 -> sprintf "dim2(a): valid=[%d..[ got=%d" n (Array2.dim2 a)
+    | -8 ->
+        sprintf "dim1(work): valid=[%d..[ got=%d" (max 1 n) (Array1.dim work)
+    | n -> raise (InternalError (sprintf "%s: error code %d" loc n))
+  in
+  invalid_arg (sprintf "%s: %s" loc msg)
+
+let orgqr_get_params loc ?m ?n ?k ~tau ~ar ~ac a =
+  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 k = get_dim_vec loc tau_str 1 1 tau k_str k in
+  m, n, k
+
+
 (* GELS? - Auxiliary Functions *)
 
 let gelsX_err loc gelsX_min_work ar a m n lwork nrhs br b err =
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.