/pliant/math/matrix.pli
1	# Copyright  Hubert Tonneau  hubert.tonneau@pliant.cx
2	#
3	# This program is free software; you can redistribute it and/or
4	# modify it under the terms of the GNU General Public License version 2
5	# as published by the Free Software Foundation.
6	#
7	# This program is distributed in the hope that it will be useful,
8	# but WITHOUT ANY WARRANTY; without even the implied warranty of
9	# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
10	# GNU General Public License for more details.
11	#
12	# You should have received a copy of the GNU General Public License
13	# version 2 along with this program; if not, write to the Free Software
14	# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
15
16	module "/pliant/language/unsafe.pli"
17	module "/pliant/language/context.pli"
18
19	constant T Float
20
21	type Matrix
22	field Array:T coefs
23	field Int ny nx
24
25	function build  m
26	arg_w Matrix m
27	m ny := 0 ; m nx := 0
28
29	method m nb_lines -> l
30	arg Matrix m ; arg Int l
31	l := m ny
32
33	method m nb_columns -> c
34	arg Matrix m ; arg Int c
35	c := m nx
36
37	method m resize nb_lines nb_columns
38	arg_rw Matrix m ; arg Int nb_lines nb_columns
39	m ny := nb_lines ; m nx := nb_columns
40	m:coefs size := nb_lines*nb_columns
41
42	method m '' line column -> coef
43	arg Matrix m ; arg Int line column ; arg_C T coef
44	check line>=0 and line<m:nb_lines
45	check column>=0 and column<m:nb_columns
46	coef :> m:coefs line*m:nx+column
47
48	function '+' a b -> r
49	arg Matrix a b r
50	check a:nb_lines=b:nb_lines and a:nb_columns=b:nb_columns
51	r resize a:nb_lines a:nb_columns
52	for (var Int y) 0 r:nb_lines-1
53	for (var Int x) 0 r:nb_columns-1
54	r y x := (a y x)+(b y x)
55
56	function '-' a b -> r
57	arg Matrix a b r
58	check a:nb_lines=b:nb_lines and a:nb_columns=b:nb_columns
59	r resize a:nb_lines a:nb_columns
60	for (var Int y) 0 r:nb_lines-1
61	for (var Int x) 0 r:nb_columns-1
62	r y x := (a y x)-(b y x)
63
64	function '*' a b -> r
65	arg Matrix a b r
66	check a:nb_columns=b:nb_lines
67	r resize a:nb_lines b:nb_columns
68	for (var Int y) 0 r:nb_lines-1
69	for (var Int x) 0 r:nb_columns-1
70	var T t := 0
71	for (var Int k) 0 a:nb_columns-1
72	t += (a y k)*(b k x)
73	r y x := t
74
75	function '*' c a -> r
76	arg T c ; arg Matrix a r
77	r resize a:nb_lines a:nb_columns
78	for (var Int y) 0 r:nb_lines-1
79	for (var Int x) 0 r:nb_columns-1
80	r y x := c*(a y x)
81
82	function transp a -> r
83	arg Matrix a r
84	r resize a:nb_columns a:nb_lines
85	for (var Int y) 0 r:nb_lines-1
86	for (var Int x) 0 r:nb_columns-1
87	r y x := a x y
88
89	method m deti py px dim -> d
90	arg Matrix m ; arg_rw Array:Int py px ; arg Int dim ; arg T d
91	if dim=1
92	return (m py:0 px:0)
93	d := 0
94	var Int extra
95	for (var Int x) dim-1 0 step -1
96	var Int t := px x ; px x := extra ; extra := t
97	var T e := (m py:(dim-1) extra)*(m deti py px dim-1)
98	if (x+(dim-1))%2<>0
99	d -= e
100	else
101	d += e
102	for x dim-1 1 step -1
103	px x := px x-1
104	px 0 := extra
105
106	function det m -> d
107	arg Matrix m ; arg T d
108	check m:nb_lines=m:nb_columns
109	var Int dim := m nb_lines ; check dim<>0
110	(var Array:Int py) size := dim
111	for (var Int i) 0 dim-1
112	py i := i
113	var Array:Int px := py
114	d := m deti py px dim
115	if pliant_debugging_level>=2
116	for (var Int i) 0 dim-1
117	check py:i=i
118	check px:i=i
119
120
121	function invert m -> r
122	arg Matrix m r
123	check m:nb_lines=m:nb_columns
124	var Int dim := m nb_lines ; check dim<>0
125	r resize dim dim
126	if dim=1
127	r 0 0 := 1/(m 0 0)
128	return
129	(var Array:Int py) size := dim-1
130	(var Array:Int px) size := dim-1
131	var T d := 0
132	for (var Int y) 0 dim-1
133	for (var Int z) 0 y-1 { py z := z }
134	for (var Int z) y+1 dim-1 { py z-1 := z }
135	for (var Int x) 0 dim-1
136	for (var Int z) 0 x-1 { px z := z }
137	for (var Int z) x+1 dim-1 { px z-1 := z }
138	var T e := m deti py px dim-1
139	if (x+y)%2<>0
140	e := -e
141	r x y := e
142	if y=0
143	d += (m y x)*e
144	if d=0
145	r resize 0 0
146	return
147	for y 0 dim-1
148	for x 0 dim-1
149	r y x /= d
150
151	function '^' m i -> r
152	arg Matrix m r ; arg Int i
153	check i=(-1)
154	r := invert m
155
156	function 'cast Status' m -> s
157	arg Matrix m ; arg Status s
158	explicit
159	s := shunt m:ny<>0 defined undefined
160
161	export Matrix '. resize' '' '. nb_lines' '. nb_columns'
162	export '+' '-' '*' '^' 'cast Status'