File: solvelin.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 (148 lines) | stat: -rw-r--r-- 6,353 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
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
(*  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

(* solve a linear system Ax = b, with input nxn matrix A and output nx1
 * matrix b *)
let solve_linear (stack : rpc_stack) (evaln : int -> unit) =
   evaln 2;
   let gen_el2 = stack#pop () in
   let gen_el1 = stack#pop () in
   match gen_el1 with
   |RpcFloatMatrixUnit (el1, u1) ->
      begin match gen_el2 with
      |RpcFloatMatrixUnit (el2, u2) ->
         let n1, m1 = Gsl_matrix.dims el1 in
         if n1 <> m1 then
            (stack#push gen_el2;
            stack#push gen_el1;
            raise (Invalid_argument "multiplier matrix must be square"))
         else
            let n2, m2 = Gsl_matrix.dims el2 in
            if m2 <> 1 then begin
               stack#push gen_el1;
               stack#push gen_el2;
               raise (Invalid_argument "resultant matrix must be a column")
            end else if n2 <> m1 then begin
               stack#push gen_el1;
               stack#push gen_el2;
               raise (Invalid_argument ("dimensions of multiplier and " ^ 
               "resultant matrices do not match"))
            end else begin
               let b = Gsl_matrix.to_array el2 in
               let x = Gsl_linalg.solve_LU (`M el1) (`A b) in
               let x_mat = Gsl_matrix.of_array x m1 1 in
               stack#push (RpcFloatMatrixUnit (x_mat, Units.div u2 u1))
            end
      |RpcComplexMatrixUnit (el2, u2) ->
         let n1, m1 = Gsl_matrix.dims el1 in
         if n1 <> m1 then
            (stack#push gen_el2;
            stack#push gen_el1;
            raise (Invalid_argument "multiplier matrix must be square"))
         else
            let n2, m2 = Gsl_matrix_complex.dims el2 in
            if m2 <> 1 then begin
               stack#push gen_el1;
               stack#push gen_el2;
               raise (Invalid_argument "resultant matrix must be a column")
            end else if n2 <> m1 then begin
               stack#push gen_el1;
               stack#push gen_el2;
               raise (Invalid_argument ("dimensions of multiplier and" ^ 
               "resultant matrices do not match"))
            end else begin
               let a_cpx = Gsl_assist.cmat_of_fmat el1 in
               let b_arr = Gsl_matrix_complex.to_array el2 in
               let b_vec = Gsl_vector_complex.of_array b_arr in
               let x = Gsl_assist.solve_complex_LU (`CM a_cpx) b_vec in
               let x_mat = Gsl_matrix_complex.of_complex_array x m1 1 in
               stack#push (RpcComplexMatrixUnit (x_mat, Units.div u2 u1))
            end
      |_ ->
         (stack#push gen_el1;
         stack#push gen_el2;
         raise (Invalid_argument "both arguments of solve_linear must be matrices"))
      end
   |RpcComplexMatrixUnit (el1, u1) ->
      begin match gen_el2 with
      |RpcFloatMatrixUnit (el2, u2) ->
         let n1, m1 = Gsl_matrix_complex.dims el1 in
         if n1 <> m1 then begin
            stack#push gen_el2;
            stack#push gen_el1;
            raise (Invalid_argument "multiplier matrix must be square")
         end else
            let n2, m2 = Gsl_matrix.dims el2 in
            if m2 <> 1 then begin
               stack#push gen_el1;
               stack#push gen_el2;
               raise (Invalid_argument "resultant matrix must be a column")
            end else if n2 <> m1 then begin
               stack#push gen_el1;
               stack#push gen_el2;
               raise (Invalid_argument ("dimensions of multiplier and" ^ 
               "resultant matrices do not match"))
            end else begin
               let b_cpx = Gsl_assist.cmat_of_fmat el2 in
               let b_arr = Gsl_matrix_complex.to_array b_cpx in
               let b_vec = Gsl_vector_complex.of_array b_arr in
               let x = Gsl_assist.solve_complex_LU (`CM el1) b_vec in
               let x_mat = Gsl_matrix_complex.of_complex_array x m1 1 in
               stack#push (RpcComplexMatrixUnit (x_mat, Units.div u2 u1))
            end
      |RpcComplexMatrixUnit (el2, u2) ->
         let n1, m1 = Gsl_matrix_complex.dims el1 in
         if n1 <> m1 then begin
            stack#push gen_el2;
            stack#push gen_el1;
            raise (Invalid_argument "multiplier matrix must be square")
         end else
            let n2, m2 = Gsl_matrix_complex.dims el2 in
            if m2 <> 1 then begin
               stack#push gen_el1;
               stack#push gen_el2;
               raise (Invalid_argument "resultant matrix must be a column")
            end else if n2 <> m1 then begin
               stack#push gen_el1;
               stack#push gen_el2;
               raise (Invalid_argument ("dimensions of multiplier and" ^ 
               "resultant matrices do not match"))
            end else begin
               let b_arr = Gsl_matrix_complex.to_array el2 in
               let b_vec = Gsl_vector_complex.of_array b_arr in
               let x = Gsl_assist.solve_complex_LU (`CM el1) b_vec in
               let x_mat = Gsl_matrix_complex.of_complex_array x m1 1 in
               stack#push (RpcComplexMatrixUnit (x_mat, Units.div u2 u1))
            end
      |_ ->
         (stack#push gen_el1;
         stack#push gen_el2;
         raise (Invalid_argument "both arguments of solve_linear must be matrices"))
      end
   |_ ->
      (stack#push gen_el1;
      stack#push gen_el2;
      raise (Invalid_argument "both arguments of solve_linear must be matrices"))



(* arch-tag: DO_NOT_CHANGE_c11268a8-d37a-4573-98db-e985cc338e38 *)