/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 Float 
 20   
 21  type Matrix 
 22    field Array:coefs 
 23    field Int ny nx 
 24   
 25  function build  m 
 26    arg_w Matrix m 
 27    ny := 0 ; nx := 0 
 28     
 29  method m nb_lines -> l 
 30    arg Matrix m ; arg Int l 
 31    := ny 
 32   
 33  method m nb_columns -> c 
 34    arg Matrix m ; arg Int c 
 35    := nx 
 36   
 37  method m resize nb_lines nb_columns 
 38    arg_rw Matrix m ; arg Int nb_lines nb_columns 
 39    ny := nb_lines ; 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 coef 
 44    check line>=and line<m:nb_lines 
 45    check column>=and column<m:nb_columns 
 46    coef :> m:coefs line*m:nx+column 
 47     
 48  function '+' a b -> r 
 49    arg Matrix r 
 50    check a:nb_lines=b:nb_lines and a:nb_columns=b:nb_columns 
 51    resize a:nb_lines a:nb_columns 
 52    for (var Int y) r:nb_lines-1 
 53      for (var Int x) r:nb_columns-1 
 54        := (x)+(x) 
 55     
 56  function '-' a b -> r 
 57    arg Matrix r 
 58    check a:nb_lines=b:nb_lines and a:nb_columns=b:nb_columns 
 59    resize a:nb_lines a:nb_columns 
 60    for (var Int y) r:nb_lines-1 
 61      for (var Int x) r:nb_columns-1 
 62        := (x)-(x) 
 63     
 64  function '*' a b -> r 
 65    arg Matrix r 
 66    check a:nb_columns=b:nb_lines 
 67    resize a:nb_lines b:nb_columns 
 68    for (var Int y) r:nb_lines-1 
 69      for (var Int x) r:nb_columns-1 
 70        var := 0 
 71        for (var Int k) a:nb_columns-1 
 72          += (k)*(x) 
 73        := t 
 74   
 75  function '*' c a -> r 
 76    arg c ; arg Matrix r 
 77    resize a:nb_lines a:nb_columns 
 78    for (var Int y) r:nb_lines-1 
 79      for (var Int x) r:nb_columns-1 
 80        := c*(x) 
 81     
 82  function transp a -> r 
 83    arg Matrix r 
 84    resize a:nb_columns a:nb_lines 
 85    for (var Int y) r:nb_lines-1 
 86      for (var Int x) r:nb_columns-1 
 87        := y 
 88   
 89  method m deti py px dim -> d 
 90    arg Matrix m ; arg_rw Array:Int py px ; arg Int dim ; arg  
 91    if dim=1 
 92      return (py:px:0) 
 93    := 0 
 94    var Int extra 
 95    for (var Int x) dim-1 0 step -1 
 96      var Int := px x ; px := extra ; extra := t 
 97      var := (py:(dim-1) extra)*(deti py px dim-1) 
 98      if (x+(dim-1))%2<>0 
 99        -= e 
 100      else 
 101        += e 
 102    for dim-1 1 step -1 
 103      px := px x-1 
 104    px := extra 
 105   
 106  function det m -> d 
 107    arg Matrix m ; arg d 
 108    check m:nb_lines=m:nb_columns 
 109    var Int dim := nb_lines ; check dim<>0 
 110    (var Array:Int py) size := dim 
 111    for (var Int i) dim-1 
 112      py := i 
 113    var Array:Int px := py 
 114    := 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 r 
 123    check m:nb_lines=m:nb_columns 
 124    var Int dim := nb_lines ; check dim<>0 
 125    resize dim dim 
 126    if dim=1 
 127      0 0 := 1/(0 0) 
 128      return 
 129    (var Array:Int py) size := dim-1 
 130    (var Array:Int px) size := dim-1 
 131    var := 0 
 132    for (var Int y) dim-1 
 133      for (var Int z) y-1 { py := z }   
 134      for (var Int z) y+dim-1 { py z-:= z } 
 135      for (var Int x) dim-1 
 136        for (var Int z) x-1 { px := z }   
 137        for (var Int z) x+dim-1 { px z-:= z } 
 138        var := deti py px dim-1 
 139        if (x+y)%2<>0 
 140          := -e 
 141        := e 
 142        if y=0 
 143          += (x)*e 
 144    if d=0 
 145      resize 0 0 
 146      return 
 147    for dim-1 
 148      for dim-1 
 149        /= d 
 150   
 151  function '^' m i -> r 
 152    arg Matrix r ; arg Int i 
 153    check i=(-1) 
 154    := invert m 
 155   
 156  function 'cast Status' m -> s 
 157    arg Matrix m ; arg Status s 
 158    explicit 
 159    := shunt m:ny<>0 defined undefined 
 160   
 161  export Matrix '. resize' '' '. nb_lines' '. nb_columns' 
 162  export '+' '-' '*' '^' 'cast Status'