Markus Mottl avatar Markus Mottl committed 0700f53

Added tpmv and tpsv

Comments (0)

Files changed (9)

+2008-11-02:  Added new BLAS functions:
+
+               * tpmv
+               * tpsv
+
 2008-11-01:  Added new BLAS function:
 
                * syr2k
-release-4-6-4
+release-4-6-5
 name="lacaml"
-version="4.6.4"
+version="4.6.5"
 description="LACAML - BLAS/LAPACK-interface for OCaml"
 
 requires="lacaml.core"
   direct_trsv ar ac a n uplo_char trans_char diag_char ofsx incx x
 
 
+(* TPMV *)
+
+external direct_tpmv :
+  int -> (* OFSAP *)
+  vec -> (* AP *)
+  int -> (* N *)
+  char -> (* UPLO *)
+  char -> (* TRANS *)
+  char -> (* DIAG *)
+  int -> (* OFSX *)
+  int -> (* INCX *)
+  vec (* X *)
+  -> unit = "lacaml_NPRECtpmv_stub_bc" "lacaml_NPRECtpmv_stub"
+
+let tpmv ?n ?(trans = `N) ?(diag = `N) ?(up = true) ?ofsap ap ?ofsx ?incx x =
+  let loc = "Lacaml.Impl.NPREC.tpmv" in
+  let n, ofsap, ofsx, incx, uplo_char, trans_char, diag_char =
+    tpXv_get_params loc ofsap ap ?n ofsx incx x up trans diag
+  in
+  direct_tpmv ofsap ap n uplo_char trans_char diag_char ofsx incx x
+
+
+(* TPSV *)
+
+external direct_tpsv :
+  int -> (* OFSAP *)
+  vec -> (* AP *)
+  int -> (* N *)
+  char -> (* UPLO *)
+  char -> (* TRANS *)
+  char -> (* DIAG *)
+  int -> (* OFSX *)
+  int -> (* INCX *)
+  vec (* X *)
+  -> unit = "lacaml_NPRECtpsv_stub_bc" "lacaml_NPRECtpsv_stub"
+
+let tpsv ?n ?(trans = `N) ?(diag = `N) ?(up = true) ?ofsap ap ?ofsx ?incx x =
+  let loc = "Lacaml.Impl.NPREC.tpsv" in
+  let n, ofsap, ofsx, incx, uplo_char, trans_char, diag_char =
+    tpXv_get_params loc ofsap ap ?n ofsx incx x up trans diag
+  in
+  direct_tpsv ofsap ap n uplo_char trans_char diag_char ofsx incx x
+
+
 (* BLAS-3 *)
 
 (* GEMM *)

lib/impl_SDCZ.mli

     @param ofsx default = 1
     @param incx default = 1 *)
 
+val tpmv :
+  ?n : int ->
+  ?trans : trans3 ->
+  ?diag : diag ->
+  ?up : bool ->
+  ?ofsap : int ->
+  vec ->
+  ?ofsx : int ->
+  ?incx : int ->
+  vec
+  -> unit
+(** [tpmv ?n ?trans ?diag ?up ?ofsap ap ?ofsx ?incx x]
+    see BLAS documentation!
+    @param n default = dimension of packed triangular matrix [ap]
+    @param trans default = `N
+    @param diag default = false (not a unit triangular matrix)
+    @param up default = true (upper triangular portion of [ap] is accessed)
+    @param ofsap default = 1
+    @param ofsx default = 1
+    @param incx default = 1 *)
+
+val tpsv :
+  ?n : int ->
+  ?trans : trans3 ->
+  ?diag : diag ->
+  ?up : bool ->
+  ?ofsap : int ->
+  vec ->
+  ?ofsx : int ->
+  ?incx : int ->
+  vec
+  -> unit
+(** [tpsv ?n ?trans ?diag ?up ?ofsap ap ?ofsx ?incx x]
+    see BLAS documentation!
+    @param n default = dimension of packed triangular matrix [ap]
+    @param trans default = `N
+    @param diag default = false (not a unit triangular matrix)
+    @param up default = true (upper triangular portion of [ap] is accessed)
+    @param ofsap default = 1
+    @param ofsx default = 1
+    @param incx default = 1 *)
+
 
 (** {6 BLAS-3 interface} *)
 

lib/impl_SDCZ_c.c

 }
 
 
+/** TPMV */
+
+extern void FUN(tpmv)(
+  char *UPLO,
+  char *TRANS,
+  char *DIAG,
+  integer *N,
+  NUMBER *AP,
+  NUMBER *X, integer *INCX);
+
+CAMLprim value LFUN(tpmv_stub)(
+  value vOFSAP,
+  value vAP,
+  value vN,
+  value vUPLO,
+  value vTRANS,
+  value vDIAG,
+  value vOFSX, value vINCX, value vX)
+{
+  CAMLparam2(vAP, vX);
+
+  char GET_INT(UPLO),
+       GET_INT(TRANS),
+       GET_INT(DIAG);
+
+  integer GET_INT(N),
+          GET_INT(INCX);
+
+  VEC_PARAMS(AP);
+  VEC_PARAMS(X);
+
+  caml_enter_blocking_section();  /* Allow other threads */
+  FUN(tpmv)(
+    &UPLO,
+    &TRANS,
+    &DIAG,
+    &N,
+    AP_data,
+    X_data, &INCX);
+  caml_leave_blocking_section();  /* Disallow other threads */
+
+  CAMLreturn(Val_unit);
+}
+
+CAMLprim value LFUN(tpmv_stub_bc)(value *argv, int argn)
+{
+  return
+    LFUN(tpmv_stub)(
+      argv[0], argv[1], argv[2], argv[3], argv[4], argv[5], argv[6],
+      argv[7], argv[8]);
+}
+
+
+/** TPSV */
+
+extern void FUN(tpsv)(
+  char *UPLO,
+  char *TRANS,
+  char *DIAG,
+  integer *N,
+  NUMBER *AP,
+  NUMBER *X, integer *INCX);
+
+CAMLprim value LFUN(tpsv_stub)(
+  value vOFSAP,
+  value vAP,
+  value vN,
+  value vUPLO,
+  value vTRANS,
+  value vDIAG,
+  value vOFSX, value vINCX, value vX)
+{
+  CAMLparam2(vAP, vX);
+
+  char GET_INT(UPLO),
+       GET_INT(TRANS),
+       GET_INT(DIAG);
+
+  integer GET_INT(N),
+          GET_INT(INCX);
+
+  VEC_PARAMS(AP);
+  VEC_PARAMS(X);
+
+  caml_enter_blocking_section();  /* Allow other threads */
+  FUN(tpsv)(
+    &UPLO,
+    &TRANS,
+    &DIAG,
+    &N,
+    AP_data,
+    X_data, &INCX);
+  caml_leave_blocking_section();  /* Disallow other threads */
+
+  CAMLreturn(Val_unit);
+}
+
+CAMLprim value LFUN(tpsv_stub_bc)(value *argv, int argn)
+{
+  return
+    LFUN(tpsv_stub)(
+      argv[0], argv[1], argv[2], argv[3], argv[4], argv[5], argv[6],
+      argv[7], argv[8]);
+}
+
+
 /** TODO: SPMV */
 
 /** TODO: TBMV */
 
-/** TODO: TPMV */
-
 /** TODO: TBSV */
 
-/** TODO: TPSV */
-
 
 /*** BLAS-3 */
 
     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) / 2 <> n_vec then
-    let loc = "Lacaml.Impl.NPREC.Mat.unpacked" in
-    failwith (sprintf "%s: illegal vector length: %d" loc n_vec)
+let unpacked ?(up = true) ?n src =
+  let loc = "Lacaml.Impl.NPREC.Mat.unpacked" in
+  let n = get_unpacked_dim loc ?n (Array1.dim src) in
+  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
-    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
+    for c = 1 to n do
+      for r = c to n do
+        a.{r, c} <- src.{!pos_ref};
+        incr pos_ref;
       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
+    done;
+  a
 
 external direct_scal_mat :
   int -> (* M *)
     @param ac default = [1]
 *)
 
-val unpacked : ?up : bool -> vec -> mat
+val unpacked : ?up : bool -> ?n : int -> vec -> mat
 (** [unpacked ?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]
+    @param n default = [Vec.dim x]
 *)
 
 
 
 (* Name information *)
 let a_str = "a"
+let ap_str = "ap"
 let b_str = "b"
 let c_str = "c"
 let k_str = "k"
       else work, lwork
   | None -> vec_create opt_lwork, opt_lwork
 
+let calc_unpacked_dim loc n_vec =
+  let n = truncate (sqrt (float (8 * n_vec + 1)) /. 2.) in
+  if (n * n + n) / 2 <> n_vec then
+    failwith (sprintf "%s: illegal vector length: %d" loc n_vec)
+  else n
+
+(* Calculate the dimension of a packed square matrix given the vector length *)
+let get_unpacked_dim loc ?n n_vec =
+  match n with
+  | None -> calc_unpacked_dim loc n_vec
+  | Some n ->
+      let n_unpacked = calc_unpacked_dim loc n_vec in
+      if n < 0 || n > n_unpacked then
+        invalid_arg (sprintf "%s: n: valid=[0..%d] got=%d" loc n_unpacked n)
+      else n
+
 (* Fetches problem-dependent parameters for LAPACK-functions *)
 external ilaenv : int -> string -> string -> int -> int -> int -> int -> int
   = "lacaml_ilaenv_stub_bc" "lacaml_ilaenv_stub" "noalloc"
   check_vec loc x_str x (ofsx + (n - 1) * abs incx);
   n, ofsx, incx, get_uplo_char up, trans_char, diag_char
 
+(* tp?v -- auxiliary functions *)
+let tpXv_get_params loc ofsap ap ?n ofsx incx x up trans unit_triangular =
+  let ofsap = get_ofs loc ap_str ofsap in
+  let n = get_unpacked_dim loc ?n (Array1.dim ap - ofsap + 1) in
+  let trans_char = get_trans_char trans in
+  let diag_char = get_diag_char unit_triangular in
+  let ofsx, incx = get_vec_geom loc x_str ofsx incx in
+  check_vec loc x_str x (ofsx + (n - 1) * abs incx);
+  n, ofsap, ofsx, incx, get_uplo_char up, trans_char, diag_char
+
 (* gemm -- auxiliary functions *)
 
 let get_c loc mat_create cr cc c m n = get_mat loc c_str mat_create cr cc c m n
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.