share/colnew/lisp/approx.lisp

   1 ;;; Compiled by f2cl version:
   2 ;;; ("f2cl1.l,v 1.221 2010/05/26 19:25:52 rtoy Exp $"
   3 ;;;  "f2cl2.l,v 1.37 2008/02/22 22:19:33 rtoy Exp $"
   4 ;;;  "f2cl3.l,v 1.6 2008/02/22 22:19:33 rtoy Exp $"
   5 ;;;  "f2cl4.l,v 1.7 2008/02/22 22:19:34 rtoy Exp $"
   6 ;;;  "f2cl5.l,v 1.204 2010/02/23 05:21:30 rtoy Exp $"
   7 ;;;  "f2cl6.l,v 1.48 2008/08/24 00:56:27 rtoy Exp $"
   8 ;;;  "macros.l,v 1.114 2010/05/17 01:42:14 rtoy Exp $")
   9
  10 ;;; Using Lisp CMU Common Lisp CVS Head 2010-05-25 18:21:07 (20A Unicode)
  11 ;;;
  12 ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  13 ;;;           (:coerce-assigns :as-needed) (:array-type ':array)
  14 ;;;           (:array-slicing t) (:declare-common nil)
  15 ;;;           (:float-format double-float))
  16
  17 (in-package :colnew)
  18
  19
  20 (defun approx (i x zval a coef xi n z dmz k ncomp mmax m mstar mode dmval modm)
  21   (declare (type (array f2cl-lib:integer4 (*)) m)
  22            (type (array double-float (*)) a)
  23            (type (array double-float (*)) dmval dmz z xi coef zval)
  24            (type double-float x)
  25            (type (f2cl-lib:integer4) modm mode mstar mmax ncomp k n i))
  26   (let ()
  27     (symbol-macrolet ((precis (aref (colout-part-0 *colout-common-block*) 0))
  28                       (iout (aref (colout-part-1 *colout-common-block*) 0))
  29                       (iprint (aref (colout-part-1 *colout-common-block*) 1)))
  30       (f2cl-lib:with-multi-array-data
  31           ((zval double-float zval-%data% zval-%offset%)
  32            (coef double-float coef-%data% coef-%offset%)
  33            (xi double-float xi-%data% xi-%offset%)
  34            (z double-float z-%data% z-%offset%)
  35            (dmz double-float dmz-%data% dmz-%offset%)
  36            (dmval double-float dmval-%data% dmval-%offset%)
  37            (a double-float a-%data% a-%offset%)
  38            (m f2cl-lib:integer4 m-%data% m-%offset%))
  39         (prog ((fact 0.0) (lb 0) (ll 0) (zsum 0.0) (ind 0) (mj 0) (jcomp 0)
  40                (idmz 0) (ir 0) (s 0.0) (iright 0) (l 0) (ileft 0) (j 0) (iz 0)
  41                (dm (make-array 7 :element-type 'double-float))
  42                (bm (make-array 4 :element-type 'double-float)))
  43           (declare (type (array double-float (4)) bm)
  44                    (type (array double-float (7)) dm)
  45                    (type (f2cl-lib:integer4) iz j ileft l iright ir idmz jcomp
  46                                              mj ind ll lb)
  47                    (type double-float s zsum fact))
  48           (f2cl-lib:computed-goto (label10 label30 label80 label90) mode)
  49          label10
  50           (setf x (f2cl-lib:fref xi-%data% (i) ((1 1)) xi-%offset%))
  51           (setf iz (f2cl-lib:int-mul (f2cl-lib:int-sub i 1) mstar))
  52           (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
  53                         ((> j mstar) nil)
  54             (tagbody
  55               (setf iz (f2cl-lib:int-add iz 1))
  56               (setf (f2cl-lib:fref zval-%data% (j) ((1 1)) zval-%offset%)
  57                       (f2cl-lib:fref z-%data% (iz) ((1 1)) z-%offset%))
  58              label20))
  59           (go end_label)
  60          label30
  61           (if
  62            (and
  63             (>= x (- (f2cl-lib:fref xi-%data% (1) ((1 1)) xi-%offset%) precis))
  64             (<= x
  65                 (+
  66                  (f2cl-lib:fref xi-%data%
  67                                 ((f2cl-lib:int-add n 1))
  68                                 ((1 1))
  69                                 xi-%offset%)
  70                  precis)))
  71            (go label40))
  72           (if (< iprint 1)
  73               (f2cl-lib:fformat iout
  74                                 (" ****** DOMAIN ERROR IN APPROX ******" "~%"
  75                                  " X =" 1 (("~20,10,2,0,'*,,'DE")) "   ALEFT ="
  76                                  1 (("~20,10,2,0,'*,,'DE")) "   ARIGHT =" 1
  77                                  (("~20,10,2,0,'*,,'DE")) "~%")
  78                                 x
  79                                 (f2cl-lib:fref xi-%data%
  80                                                (1)
  81                                                ((1 1))
  82                                                xi-%offset%)
  83                                 (f2cl-lib:fref xi-%data%
  84                                                ((f2cl-lib:int-add n 1))
  85                                                ((1 1))
  86                                                xi-%offset%)))
  87           (if (< x (f2cl-lib:fref xi-%data% (1) ((1 1)) xi-%offset%))
  88               (setf x (f2cl-lib:fref xi-%data% (1) ((1 1)) xi-%offset%)))
  89           (if
  90            (> x
  91               (f2cl-lib:fref xi-%data%
  92                              ((f2cl-lib:int-add n 1))
  93                              ((1 1))
  94                              xi-%offset%))
  95            (setf x
  96                    (f2cl-lib:fref xi-%data%
  97                                   ((f2cl-lib:int-add n 1))
  98                                   ((1 1))
  99                                   xi-%offset%)))
 100          label40
 101           (if (or (> i n) (< i 1))
 102               (setf i (the f2cl-lib:integer4 (truncate (+ n 1) 2))))
 103           (setf ileft i)
 104           (if (< x (f2cl-lib:fref xi-%data% (ileft) ((1 1)) xi-%offset%))
 105               (go label60))
 106           (f2cl-lib:fdo (l ileft (f2cl-lib:int-add l 1))
 107                         ((> l n) nil)
 108             (tagbody
 109               (setf i l)
 110               (if
 111                (< x
 112                   (f2cl-lib:fref xi-%data%
 113                                  ((f2cl-lib:int-add l 1))
 114                                  ((1 1))
 115                                  xi-%offset%))
 116                (go label80))
 117              label50))
 118           (go label80)
 119          label60
 120           (setf iright (f2cl-lib:int-sub ileft 1))
 121           (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 122                         ((> l iright) nil)
 123             (tagbody
 124               (setf i (f2cl-lib:int-sub (f2cl-lib:int-add iright 1) l))
 125               (if (>= x (f2cl-lib:fref xi-%data% (i) ((1 1)) xi-%offset%))
 126                   (go label80))
 127              label70))
 128          label80
 129           (setf s
 130                   (/ (- x (f2cl-lib:fref xi-%data% (i) ((1 1)) xi-%offset%))
 131                      (-
 132                       (f2cl-lib:fref xi-%data%
 133                                      ((f2cl-lib:int-add i 1))
 134                                      ((1 1))
 135                                      xi-%offset%)
 136                       (f2cl-lib:fref xi-%data% (i) ((1 1)) xi-%offset%))))
 137           (rkbas s coef k mmax a dm modm)
 138          label90
 139           (setf (f2cl-lib:fref bm (1) ((1 4)))
 140                   (- x (f2cl-lib:fref xi-%data% (i) ((1 1)) xi-%offset%)))
 141           (f2cl-lib:fdo (l 2 (f2cl-lib:int-add l 1))
 142                         ((> l mmax) nil)
 143             (tagbody
 144               (setf (f2cl-lib:fref bm (l) ((1 4)))
 145                       (/ (f2cl-lib:fref bm (1) ((1 4))) (f2cl-lib:dfloat l)))
 146              label95))
 147          label100
 148           (setf ir 1)
 149           (setf iz
 150                   (f2cl-lib:int-add
 151                    (f2cl-lib:int-mul (f2cl-lib:int-sub i 1) mstar)
 152                    1))
 153           (setf idmz (f2cl-lib:int-mul (f2cl-lib:int-sub i 1) k ncomp))
 154           (f2cl-lib:fdo (jcomp 1 (f2cl-lib:int-add jcomp 1))
 155                         ((> jcomp ncomp) nil)
 156             (tagbody
 157               (setf mj (f2cl-lib:fref m-%data% (jcomp) ((1 1)) m-%offset%))
 158               (setf ir (f2cl-lib:int-add ir mj))
 159               (setf iz (f2cl-lib:int-add iz mj))
 160               (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 161                             ((> l mj) nil)
 162                 (tagbody
 163                   (setf ind (f2cl-lib:int-add idmz jcomp))
 164                   (setf zsum 0.0)
 165                   (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 166                                 ((> j k) nil)
 167                     (tagbody
 168                       (setf zsum
 169                               (+ zsum
 170                                  (*
 171                                   (f2cl-lib:fref a-%data%
 172                                                  (j l)
 173                                                  ((1 7) (1 1))
 174                                                  a-%offset%)
 175                                   (f2cl-lib:fref dmz-%data%
 176                                                  (ind)
 177                                                  ((1 1))
 178                                                  dmz-%offset%))))
 179                      label110
 180                       (setf ind (f2cl-lib:int-add ind ncomp))))
 181                   (f2cl-lib:fdo (ll 1 (f2cl-lib:int-add ll 1))
 182                                 ((> ll l) nil)
 183                     (tagbody
 184                       (setf lb (f2cl-lib:int-sub (f2cl-lib:int-add l 1) ll))
 185                      label120
 186                       (setf zsum
 187                               (+ (* zsum (f2cl-lib:fref bm (lb) ((1 4))))
 188                                  (f2cl-lib:fref z-%data%
 189                                                 ((f2cl-lib:int-sub iz ll))
 190                                                 ((1 1))
 191                                                 z-%offset%)))))
 192                  label130
 193                   (setf (f2cl-lib:fref zval-%data%
 194                                        ((f2cl-lib:int-sub ir l))
 195                                        ((1 1))
 196                                        zval-%offset%)
 197                           zsum)))
 198              label140))
 199           (if (= modm 0) (go end_label))
 200           (f2cl-lib:fdo (jcomp 1 (f2cl-lib:int-add jcomp 1))
 201                         ((> jcomp ncomp) nil)
 202             (tagbody
 203              label150
 204               (setf (f2cl-lib:fref dmval-%data% (jcomp) ((1 1)) dmval-%offset%)
 205                       0.0)))
 206           (setf idmz (f2cl-lib:int-add idmz 1))
 207           (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 208                         ((> j k) nil)
 209             (tagbody
 210               (setf fact (f2cl-lib:fref dm (j) ((1 7))))
 211               (f2cl-lib:fdo (jcomp 1 (f2cl-lib:int-add jcomp 1))
 212                             ((> jcomp ncomp) nil)
 213                 (tagbody
 214                   (setf (f2cl-lib:fref dmval-%data%
 215                                        (jcomp)
 216                                        ((1 1))
 217                                        dmval-%offset%)
 218                           (+
 219                            (f2cl-lib:fref dmval-%data%
 220                                           (jcomp)
 221                                           ((1 1))
 222                                           dmval-%offset%)
 223                            (* fact
 224                               (f2cl-lib:fref dmz-%data%
 225                                              (idmz)
 226                                              ((1 1))
 227                                              dmz-%offset%))))
 228                   (setf idmz (f2cl-lib:int-add idmz 1))
 229                  label160))
 230              label170))
 231           (go end_label)
 232          end_label
 233           (return
 234            (values i
 235                    x
 236                    nil
 237                    nil
 238                    nil
 239                    nil
 240                    nil
 241                    nil
 242                    nil
 243                    nil
 244                    nil
 245                    nil
 246                    nil
 247                    nil
 248                    nil
 249                    nil
 250                    nil)))))))
 251
 252 (in-package #-gcl #:cl-user #+gcl "CL-USER")
 253 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 254 (eval-when (:load-toplevel :compile-toplevel :execute)
 255   (setf (gethash 'fortran-to-lisp::approx
 256                  fortran-to-lisp::*f2cl-function-info*)
 257           (fortran-to-lisp::make-f2cl-finfo
 258            :arg-types '((fortran-to-lisp::integer4) double-float
 259                         (array double-float (1)) (array double-float (7))
 260                         (array double-float (1)) (array double-float (1))
 261                         (fortran-to-lisp::integer4) (array double-float (1))
 262                         (array double-float (1)) (fortran-to-lisp::integer4)
 263                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 264                         (array fortran-to-lisp::integer4 (1))
 265                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 266                         (array double-float (1)) (fortran-to-lisp::integer4))
 267            :return-values '(fortran-to-lisp::i fortran-to-lisp::x nil nil nil
 268                             nil nil nil nil nil nil nil nil nil nil nil nil)
 269            :calls '(fortran-to-lisp::rkbas))))
 270