share/colnew/colnew-if.lisp

   1 ;;; -*- Mode: lisp; Package: CL-USER -*-
   2
   3 (in-package :maxima)
   4
   5 ;; Raw interface to colnew where essentially all of the parameters for
   6 ;; colnew are exposed.  The f, df, g, dg, and init-guess parameters
   7 ;; must be functions.
   8 (defun $colnew_expert (ncomp m aleft aright zeta ipar ltol tol fixpnt ispace fspace
   9                        iflag f df g dg init-guess)
  10   (flet ((convert-to-array (mlist atype)
  11            ;; Convert the Maxima lists to arrays
  12            (make-array (length (cdr mlist)) :element-type atype
  13                        :initial-contents (cdr mlist)))
  14          (convert-to-mlist (v)
  15            (list* '(mlist) (coerce v 'list))))
  16     (let ((a-m (convert-to-array m 'f2cl-lib:integer4))
  17           (a-zeta (convert-to-array zeta 'double-float))
  18           (a-ipar (convert-to-array ipar 'f2cl-lib:integer4))
  19           (a-ltol (convert-to-array ltol 'f2cl-lib:integer4))
  20           (a-tol (convert-to-array tol 'double-float))
  21           (a-fixpnt (convert-to-array fixpnt 'double-float))
  22           (a-ispace (convert-to-array ispace 'f2cl-lib:integer4))
  23           (a-fspace (convert-to-array fspace 'double-float))
  24           (f (coerce-float-fun f))
  25           (df (coerce-float-fun df))
  26           (g (coerce-float-fun g))
  27           (dg (coerce-float-fun dg))
  28           (nz (reduce #'+ (cdr m))))
  29       (flet ((fsub (x z f-array)
  30                (declare (type (cl:array double-float (*)) z f-array)
  31                         (double-float x))
  32                (let ((res (apply f x (coerce (subseq z 0 nz) 'list))))
  33                  (when (listp res)
  34                    (loop for k from 0
  35                       for ff in (cdr res)
  36                       do
  37                       (setf (aref f-array k) ($float ff)))))
  38                (values nil nil f-array))
  39              (dfsub (x z df-array)
  40                (declare (type (cl:array double-float (*)) z df-array)
  41                         (double-float x))
  42                (let ((res (apply df x (coerce (subseq z 0 nz) 'list))))
  43                  ;; res is a Maxima matrix.
  44                  (when (listp res)
  45                    (loop for row from 1
  46                       for row-list in (cdr res)
  47                       do
  48                       (loop for col from 1
  49                          for element in (cdr row-list)
  50                          do
  51                          (setf (f2cl-lib::fref df-array (row col) ((1 ncomp) (1 *))) ($float element))))))
  52                (values nil nil df-array))
  53              (gsub (i z dummy)
  54                (declare (type (cl:array double-float (*)) z)
  55                         (type f2cl-lib:integer4 i))
  56                (declare (ignore dummy))
  57                (let ((res (apply g i (coerce (subseq z 0 nz) 'list))))
  58                  (values nil nil res)))
  59              (dgsub (i z dg-array)
  60                (declare (type (cl:array double-float (*)) z dg-array)
  61                         (type f2cl-lib:integer4 i))
  62                (let ((res (apply dg i (coerce (subseq z 0 nz) 'list))))
  63                  (loop for k from 0
  64                     for ff in (cdr res)
  65                     do
  66                     (setf (aref dg-array k) ($float ff))))
  67                (values nil nil dg-array))
  68              (guess (x z dmval)
  69                (declare (double-float x)
  70                         (type (cl:array double-float (*)) z dmval))
  71                (let* ((res ($float (mcall init-guess x)))
  72                       (new-z (second res))
  73                       (new-d (third res)))
  74                  (loop for k from 0
  75                     for ff in (cdr new-d)
  76                     do (setf (aref dmval k) ($float ff)))
  77                  (loop for k from 0
  78                     for ff in (cdr new-z)
  79                     do (setf (aref z k) ($float ff))))
  80                (values nil z dmval)))
  81         (multiple-value-bind (z-ncomp z-m z-aleft z-aright z-zeta z-ipar z-ltol z-tol
  82                                       z-fixpnt z-ispace z-fspace o-iflag)
  83             (colnew:colnew ncomp
  84                            a-m
  85                            ($float aleft)
  86                            ($float aright)
  87                            a-zeta
  88                            a-ipar
  89                            a-ltol
  90                            a-tol
  91                            a-fixpnt
  92                            a-ispace
  93                            a-fspace
  94                            iflag
  95                            #'fsub
  96                            #'dfsub
  97                            #'gsub
  98                            #'dgsub
  99                            #'guess)
 100           (declare (ignore z-ncomp z-m z-aleft z-aright z-zeta z-ipar z-ltol z-tol
 101                            z-fixpnt z-ispace z-fspace))
 102           (list '(mlist)
 103                 o-iflag
 104                 (convert-to-mlist a-fspace)
 105                 (convert-to-mlist a-ispace)))))))
 106
 107 ;; This is a slight extension of the actual APPSLN function.  The X
 108 ;; must be a list.  The output is then a list of the z values for each
 109 ;; x.  The z values themselves are lists.
 110 (defun $colnew_appsln (x zlen fspace ispace)
 111   (flet ((convert-to-array (mlist atype)
 112            (make-array (length (cdr mlist)) :element-type atype
 113                        :initial-contents (cdr mlist)))
 114          (convert-to-mlist (v)
 115            (list* '(mlist) (coerce v 'list))))
 116     (let ((a-fspace (convert-to-array fspace 'double-float))
 117           (a-ispace (convert-to-array ispace 'f2cl-lib:integer4))
 118           (z (make-array zlen :element-type 'double-float))
 119           (result nil))
 120       (dolist (pnt (cdr x))
 121         (colnew:appsln ($float pnt) z a-fspace a-ispace)
 122         (push (convert-to-mlist z) result))
 123       (list* '(mlist) (nreverse result)))))