← Back to branch summary

~ubuntu-branches/ubuntu/raring/orpie/raring

« back to all changes in this revision

Viewing changes to inv.ml

  • Committer: Bazaar Package Importer
  • Author(s): Uwe Steinmann
  • Date: 2004-09-20 14:18:45 UTC
  • Revision ID: james.westby@ubuntu.com-20040920141845-j092sbrg4hd0nfsf
Tags: upstream-1.4.1
Import upstream version 1.4.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
inv.ml
 
1
(*  Orpie -- a stack-based RPN calculator for the console
 
2
 *  Copyright (C) 2003-2004  Paul Pelzl
 
3
 *
 
4
 *  This program is free software; you can redistribute it and/or modify
 
5
 *  it under the terms of the GNU General Public License as published by
 
6
 *  the Free Software Foundation; either version 2 of the License, or
 
7
 *  (at your option) any later version.
 
8
 *
 
9
 *  This program is distributed in the hope that it will be useful,
 
10
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
 *  GNU General Public License for more details.
 
13
 *
 
14
 *  You should have received a copy of the GNU General Public License
 
15
 *  along with this program; if not, write to the Free Software
 
16
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
17
 *
 
18
 *  Please send bug reports, patches, etc. to Paul Pelzl at 
 
19
 *  <pelzlpj@eecs.umich.edu>.
 
20
 *)
 
21
 
 
22
open Rpc_stack
 
23
open Gsl_error
 
24
open Gsl_assist
 
25
 
 
26
 
 
27
let inv (stack : rpc_stack) (evaln : int -> unit) =
 
28
   evaln 1;
 
29
   let gen_el = stack#pop () in
 
30
   match gen_el with
 
31
   |RpcFloatUnit el ->
 
32
      stack#push (RpcFloatUnit (Units.div (funit_of_float 1.0) el))
 
33
   |RpcComplexUnit el ->
 
34
      stack#push (RpcComplexUnit (Units.div 
 
35
      (cunit_of_cpx Complex.one) el))
 
36
   |RpcFloatMatrixUnit (el, uu) ->
 
37
      let new_unit = Units.pow uu (-1.0) in
 
38
      let n, m = (Gsl_matrix.dims el) in
 
39
      if n = m then
 
40
         let copy_el1 = Gsl_matrix.copy el
 
41
         and copy_el2 = Gsl_matrix.copy el
 
42
         and perm     = Gsl_permut.create m
 
43
         and inv      = Gsl_matrix.create m m
 
44
         and vv       = Gsl_matrix.create m m
 
45
         and ss       = Gsl_vector.create m
 
46
         and work     = Gsl_vector.create m in
 
47
         begin
 
48
            (* test for singular matrix *)
 
49
            (* first factor out matrix norm, since the GSL SVD algorithm has
 
50
             * issues with large magnitude matrices *)
 
51
            let norm = Gsl_assist.one_norm copy_el1 in
 
52
            Gsl_matrix.scale copy_el1 (1.0 /. norm);
 
53
            (* now compute condition number as largest singular value
 
54
             * divided by smallest singular value *)
 
55
            Gsl_linalg._SV_decomp (`M copy_el1) (`M vv) (`V ss) (`V work);
 
56
            let condition_number = 
 
57
               (Gsl_vector.get ss 0) /. (Gsl_vector.get ss (pred m)) 
 
58
            in
 
59
            (* if the condition number is too large for machine precision,
 
60
             * then abort with error *)
 
61
            if condition_number > 1e14 then
 
62
               (stack#push gen_el;
 
63
               raise (Invalid_argument "cannot invert ill-conditioned matrix"))
 
64
            else
 
65
               let sign = Gsl_linalg._LU_decomp (`M copy_el2) perm in
 
66
               (Gsl_linalg._LU_invert (`M copy_el2) perm (`M inv);
 
67
               stack#push (RpcFloatMatrixUnit (inv, new_unit)))
 
68
         end
 
69
      else
 
70
         (stack#push gen_el;
 
71
         raise (Invalid_argument "cannot invert non-square matrix"))
 
72
   |RpcComplexMatrixUnit (el, uu) ->
 
73
      let new_unit = Units.pow uu (-1.0) in
 
74
      let n, m = (Gsl_matrix_complex.dims el) in
 
75
      if n = m then
 
76
         let copy_el = Gsl_vectmat.cmat_convert ~protect:true (`CM el) and
 
77
         perm = Gsl_permut.create m and
 
78
         inv = Gsl_matrix_complex.create m m in
 
79
         try
 
80
            let sign = Gsl_linalg.complex_LU_decomp copy_el perm in
 
81
            Gsl_linalg.complex_LU_invert copy_el perm (`CM inv);
 
82
            stack#push (RpcComplexMatrixUnit (inv, new_unit))
 
83
         with Gsl_exn _ ->
 
84
            (stack#push gen_el;
 
85
            raise (Invalid_argument "cannot invert singular matrix"))
 
86
      else
 
87
         (stack#push gen_el;
 
88
         raise (Invalid_argument "cannot invert non-square matrix"))
 
89
   |_ ->
 
90
      (stack#push gen_el;
 
91
      raise (Invalid_argument "inversion is undefined for this data
 
92
      type"))
 
93
 
 
94
 
 
95
 
 
96
(* arch-tag: DO_NOT_CHANGE_d8ce074c-3d77-4448-b3c6-9e239b853aad *)