File: inv.ml

package info (click to toggle)
orpie 1.5.2-2
  • links: PTS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 1,924 kB
  • ctags: 2,720
  • sloc: ml: 13,872; ansic: 3,754; makefile: 310; sh: 11; python: 11
file content (95 lines) | stat: -rw-r--r-- 3,898 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
(*  Orpie -- a fullscreen RPN calculator for the console
 *  Copyright (C) 2003-2004, 2005, 2006-2007, 2010 Paul Pelzl
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License, Version 2,
 *  as published by the Free Software Foundation.
 *
 *  This program 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 General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 *  Please send bug reports, patches, etc. to Paul Pelzl at 
 *  <pelzlpj@gmail.com>.
 *)

open Rpc_stack
open Gsl_error
open Gsl_assist


let inv (stack : rpc_stack) (evaln : int -> unit) =
   evaln 1;
   let gen_el = stack#pop () in
   match gen_el with
   |RpcFloatUnit (el, uu) ->
      stack#push (RpcFloatUnit (1.0 /. el, Units.div Units.empty_unit uu))
   |RpcComplexUnit (el, uu) ->
      stack#push (RpcComplexUnit 
         (Complex.inv el, Units.div Units.empty_unit uu))
   |RpcFloatMatrixUnit (el, uu) ->
      let new_unit = Units.pow uu (~-. 1.0) in
      let n, m = (Gsl_matrix.dims el) in
      if n = m then
         let copy_el1 = Gsl_matrix.copy el
         and copy_el2 = Gsl_matrix.copy el
         and perm     = Gsl_permut.create m
         and inv      = Gsl_matrix.create m m
         and vv       = Gsl_matrix.create m m
         and ss       = Gsl_vector.create m
         and work     = Gsl_vector.create m in
         begin
            (* test for singular matrix *)
            (* first factor out matrix norm, since the GSL SVD algorithm has
             * issues with large magnitude matrices *)
            let norm = Gsl_assist.one_norm copy_el1 in
            Gsl_matrix.scale copy_el1 (1.0 /. norm);
            (* now compute condition number as largest singular value
             * divided by smallest singular value *)
            Gsl_linalg._SV_decomp (`M copy_el1) (`M vv) (`V ss) (`V work);
            let condition_number = 
               (Gsl_vector.get ss 0) /. (Gsl_vector.get ss (pred m)) 
            in
            (* if the condition number is too large for machine precision,
             * then abort with error *)
            if condition_number > 1e14 then
               (stack#push gen_el;
               raise (Invalid_argument "cannot invert ill-conditioned matrix"))
            else
               let _ = Gsl_linalg._LU_decomp (`M copy_el2) perm in
               (Gsl_linalg._LU_invert (`M copy_el2) perm (`M inv);
               stack#push (RpcFloatMatrixUnit (inv, new_unit)))
         end
      else
         (stack#push gen_el;
         raise (Invalid_argument "cannot invert non-square matrix"))
   |RpcComplexMatrixUnit (el, uu) ->
      let new_unit = Units.pow uu (~-. 1.0) in
      let n, m = (Gsl_matrix_complex.dims el) in
      if n = m then
         let copy_el = Gsl_vectmat.cmat_convert ~protect:true (`CM el) and
         perm = Gsl_permut.create m and
         inv = Gsl_matrix_complex.create m m in
         try
            let _ = Gsl_linalg.complex_LU_decomp copy_el perm in
            Gsl_linalg.complex_LU_invert copy_el perm (`CM inv);
            stack#push (RpcComplexMatrixUnit (inv, new_unit))
         with Gsl_exn _ ->
            (stack#push gen_el;
            raise (Invalid_argument "cannot invert singular matrix"))
      else
         (stack#push gen_el;
         raise (Invalid_argument "cannot invert non-square matrix"))
   |_ ->
      (stack#push gen_el;
      raise (Invalid_argument "inversion is undefined for this data
      type"))



(* arch-tag: DO_NOT_CHANGE_d8ce074c-3d77-4448-b3c6-9e239b853aad *)