archive/src/fortdef.lsp

   1 (in-package "MAXIMA")
   2 ;; see example of specification at end of this file.
   3
   4 ;;for maxima  g(t,u) == (($g simp) $t $u)
   5 ;; (coerce-call-lisp '(($g simp) $t $u)) --> ($g $t $u)
   6 (defun coerce-call-lisp (fcall)
   7   (cond ((atom fcall) fcall)
   8         ((and (consp fcall)
   9               (consp (car fcall)))
  10          (cond ((eq (caar fcall) 'mlist)
  11                 (mapcar 'coerce-call-lisp (cdr fcall)))
  12                (t (cons (caar fcall) (cdr fcall)))))
  13         (t fcall)))
  14
  15 (defvar *fortran-types*)
  16 (setq  *fortran-types* '(($integer . "int")
  17                           ($double . "double")
  18                           ($dimension . "dimension")
  19                           ($external . "void *")
  20
  21                           ))
  22
  23 (defun get-fortran-type (type)
  24   (or (cdr (assoc type *fortran-types*)) (error "unrecognized type ~a " type)))
  25
  26 ;; if true then
  27 (defvar $fortran_force_new_compile nil)
  28 (defvar $fortran_path '((mlist) "./" "/usr/local/lib/linpack/"))
  29
  30
  31
  32 (defmacro $make_fortran (fun declarations &rest l &aux path tem locals)
  33   (let ((keys (parse_keys l '($defaults $result $assertions $locals $requires_subroutines
  34                                         $libs)
  35                           nil nil)))
  36     (setq fun (coerce-call-lisp fun))
  37     (setq declarations  (coerce-call-lisp declarations))
  38
  39     (let ((result  (coerce-call-lisp (cadr (assoc '$result keys))))
  40           (defaults  (cadr (assoc '$defaults keys)))
  41           (assertions  (cadr (assoc '$assertions keys)))
  42           (locals  (cadr (assoc '$locals keys)))
  43           (name  (with-output-to-string  (*fortran-out*)(wt (car fun))))
  44           (requires_subroutines (cadr (assoc '$requires_subroutines)))
  45           (libs (if (assoc '$libs)
  46                              (coerce (mstring(cadr (assoc '$libs))) 'string)
  47                            "-llinpack -lblass")
  48
  49                            )
  50           )
  51
  52       (when
  53        (or $fortran_force_new_compile
  54            (not
  55             (dolist (v (cdr  $fortran_path))
  56                     (if (probe-file (setq tem (si::string-concatenate v name
  57                                                                       ".exe")))
  58                         (return (setq path (namestring (truename tem))))))))
  59
  60        (make-fortran-fun fun declarations result name locals)
  61                                         ;
  62        (compiler::safe-system
  63         (si::string-concatenate
  64          "cc -c -O " name".c ; "
  65          "cc  -L/usr/local/lib " libs
  66          " -L/usr/lang/SC1.0 -lF77 -lm fony.o -o "    name".exe " name ".o"
  67                                         ;";rm -f "  name ".o " name ".c "
  68          ))
  69        (setq path (namestring
  70                    (truename (si::string-concatenate name ".exe")))))
  71       (define-maxima-function fun declarations defaults result path assertions)
  72       (list 'quote (car fun)))))
  73
  74 (defstruct fortran args requireds decls defaults results exec-path assertions
  75   subroutines
  76   libs
  77   locals)
  78
  79 (defun default-type-value (type) (or (cdr (assoc type '(($integer . 0)
  80                                                   ($double . 0.0d0)
  81                                                   ($float . 0.0s0))))
  82                                (error "no default for ~a" type)))
  83
  84 (defun $default_array (type dim1 &optional dim2 &aux res)
  85   (or (and (typep dim1 'fixnum)
  86            (or (null dim2) (typep dim2 'fixnum)))
  87       (error "dimension not fixnum"))
  88   (cond ((null dim2) (cons '(mlist)(make-list dim1 :initial-element
  89                                               (default-type-value type))))
  90         (t `(($matrix),@
  91              (dotimes (i dim1 res)
  92                       (push ($default_array type dim2) res))))))
  93
  94 (defun maxima-defaults (deflt )
  95       (declare (special *decls*))
  96    (cond ((and (consp deflt)
  97                (consp (car deflt))
  98                (eq (caar deflt) 'mequal))
  99           deflt)
 100          ((and (symbolp deflt)
 101                (let* ((type (cdr (get-arg-decl deflt *decls*)))
 102                       (tem (if (atom type) (default-type-value type)
 103                              `(($default_array) ,@ type))))
 104                  (and tem (list '(mequal) deflt tem)))))
 105          (t
 106           (error "don't know how to give default to ~a" deflt))))
 107
 108 (defun define-maxima-function (fun declarations defaults results path
 109                                    assertions)
 110   (let (args (*decls* declarations))
 111     (declare (special *decls*))
 112     (setq defaults (mapcar 'maxima-defaults (cdr defaults)))
 113     (dolist (v (cdr fun) (setq args (nreverse args)))
 114             (or  (find v defaults :key 'cadr) (push v args)))
 115     (setf (get (car fun) 'fortran) (make-fortran :args (cdr fun)
 116                                            :requireds args
 117                                            :decls declarations
 118                                            :defaults
 119                                            defaults
 120                                            :results (or results (cdr fun))
 121                                            :exec-path path
 122                                            :assertions assertions
 123                                            :locals locals
 124                                            ))
 125     (setf (get (car fun) 'mfexpr*)
 126               #'(LAMBDA (FORTRAN-ARG) (MAXIMA-FORTRAN-INVOKE FORTRAN-ARG)))
 127     ))
 128
 129 (defvar $tmp_prefix nil)
 130 (defun get-temp-path (prefix suffix &aux tem)
 131   (dotimes (i 10000)
 132            (unless (probe-file (setq tem (format nil "~a~a" prefix i suffix)))
 133                    (return tem))))
 134
 135
 136 (defun link-with-subroutines (f arg-values &aux fun name)
 137   (let ((subrs (fortran-subroutines f)))
 138     (sloop for v in (cdr (fortran-subroutines f))
 139        do (or ($listp v) (error "subr decl should be list ~M" v))
 140        (setq fun (nth 1 v))
 141        (or (and (consp fun)
 142                 (symbolp (caar fun)))
 143            (merror "bad declaration ~M" v))
 144        (setq name (caar fun))
 145        (let* ((keys (parse_keys (cddr v) '($decls) nil nil))
 146               (decls (nth 2 (assoc '$decls keys)))
 147               (val (nth (position name (fortran-args f)) arg-values)))
 148          (cond ((symbolp f)
 149                 (setq def (mget '$ff 'mexpr)))
 150
 151                )))))
 152
 153
 154
 155
 156
 157
 158 (defun maxima-fortran-invoke (arg  &aux arg-values tem io-path exec-path)
 159   (or $tmp_prefix (setq $tmp_prefix (format nil "/tmp/~a_fort"
 160                                             (si::getenv "USER"))))
 161   (let* ((f (get (caar arg) 'fortran))
 162          (arg-values-c)
 163          (default-pos
 164            (+ (length (fortran-requireds f)) 1))
 165          exec-path
 166
 167          (defaults
 168            (cond ((eql (length arg) default-pos)
 169                   (fortran-defaults f))
 170                  ((eql (length arg) (+ 1 default-pos))
 171                   (or ($listp (nth default-pos arg))
 172                       (error "defaults not a list"))
 173                   (prog1
 174                       (append (mapcar 'maxima-defaults
 175                                       (cdr(nth default-pos arg)))
 176                               (fortran-defaults f))
 177                     (setq arg (butlast arg))))
 178                  (t (error "~a needs ~a or ~a args" (caar arg)(- default-pos 1)
 179                            (- default-pos 2))))))
 180     (progv (fortran-requireds f)
 181            (mapcar 'meval* (cdr arg))
 182            ;;eval defaults
 183            (setq defaults (parse_keys defaults (Fortran-args f)
 184                                        t nil))
 185            (dolist (v (fortran-args f))
 186                    (if (setq tem (assoc v defaults))
 187                        (push (list v (second tem)) arg-values)
 188                      (push (list v (symbol-value v)) arg-values))))
 189     (setq arg-values (nreverse     arg-values))
 190     (setq exec-path
 191           (if (fortran-subroutines f)
 192               (link-with-subroutines f arg-values )
 193             (fortran-exec-path f)))
 194
 195     (with-open-file (st (setq io-path (get-temp-path $tmp_prefix "Fw"))
 196                         :direction :output)
 197       (let ((xdr (si::xdr-open st)))
 198         (dolist (v arg-values)
 199                 (push (maxima-coerce-to-xdr
 200                        (second v)
 201                        (cdr (get-arg-decl (car v)
 202                                           (fortran-decls f))))
 203                       arg-values-c)
 204                 (si::xdr-write xdr (car arg-values-c)))))
 205     (setq arg-values-c (nreverse arg-values-c))
 206     (system
 207      (format nil " ~a <  ~a > ~ao ; rm -f ~a"
 208              exec-path
 209              io-path io-path io-path
 210              ))
 211    (unwind-protect
 212     (with-open-file (st (setq io-path (concatenate 'string io-path "o")))
 213       (let ((xdr (si::xdr-open st))
 214             (pairs (pairlis (fortran-args f) arg-values-c))
 215             result)
 216         (dolist (v (fortran-results f))
 217                 (push (si::xdr-read xdr (cdr (assoc v pairs))) result))
 218         (setq result (nreverse result))
 219         (cons '(mlist)
 220               (sloop::sloop for v in result
 221                             for na in (fortran-results f)
 222                             collect
 223                             (coerce-maxima-from-xdr v
 224                                                     (nth
 225                                                      (position na
 226                                                              (fortran-args f))
 227                                                          arg-values))))
 228
 229         ))
 230     (delete-file io-path)
 231
 232     )))
 233
 234 #+when-run-process-deallocates
 235 (defun maxima-fortran-invoke (arg &aux arg-values tem)
 236   (let* ((f (get (caar arg) 'fortran))
 237          (arg-values-c)
 238          (default-pos
 239            (+ (length (fortran-requireds f)) 1))
 240          (defaults
 241            (cond ((eql (length arg) default-pos)
 242                   (fortran-defaults f))
 243                  ((eql (length arg) (+ 1 default-pos))
 244                   (or ($listp (nth default-pos arg))
 245                       (error "defaults not a list"))
 246                   (prog1
 247                       (append (mapcar 'maxima-defaults
 248                                       (cdr(nth default-pos arg)))
 249                               (fortran-defaults f))
 250                     (setq arg (butlast arg))))
 251                  (t (error "~a needs ~a or ~a args" (caar arg)(- default-pos 1)
 252                            (- default-pos 2))))))
 253     (progv (fortran-requireds f)
 254            (mapcar 'meval* (cdr arg))
 255            ;;eval defaults
 256            (setq defaults (parse_keys defaults (Fortran-args f)
 257                                       t nil))
 258            (dolist (v (fortran-args f))
 259                    (if (setq tem (assoc v defaults))
 260                        (push (list v (second tem)) arg-values)
 261                      (push (list v (symbol-value v)) arg-values))))
 262     (setq arg-values (nreverse     arg-values))
 263     (let* ((tem (si::run-process (fortran-exec-path f) nil))
 264            (out (si::fp-output-stream tem))
 265            (in (si::fp-input-stream tem)))
 266       (let ((xdr (si::xdr-open out)))
 267         (dolist (v arg-values)
 268                 (push (maxima-coerce-to-xdr (second v)
 269                                             (cdr (get-arg-decl (car v)
 270                                                                (fortran-decls f))))
 271                       arg-values-c)
 272
 273                 (si::xdr-write xdr (car arg-values-c))))
 274     (setq arg-values-c (nreverse arg-values-c))
 275     ;; invoked
 276     (force-output out)
 277     (let ((xdr (si::xdr-open in))
 278           (pairs (pairlis (fortran-args f) arg-values-c))
 279           result)
 280       (dolist (v (fortran-results f))
 281               (push (si::xdr-read xdr (cdr (assoc v pairs))) result))
 282       (setq result (nreverse result))
 283       (cons '(mlist)
 284             (sloop::sloop for v in result
 285                           for na in (fortran-results f)
 286                           collect
 287                           (coerce-maxima-from-xdr v
 288                                                   (nth (position na
 289                                                                  (fortran-args f))
 290                                                        arg-values))))))))
 291
 292 (defun meval-* (lis &aux tem)
 293   (cond ((null lis) lis)
 294         (t (progv (list (caar lis))
 295                   (list (setq tem (meval* (second (car lis)))))
 296                   (cons (list (caar lis) tem)
 297                         (meval-* (cdr lis)))))))
 298
 299
 300 ;; return an list of  (list keys  value)
 301 ;; where each value is meval* in with previous keys bound to
 302 ;; previous values.
 303 (defun parse_keys (list keys meval done &aux sym tem)
 304   (cond ((null list) done)
 305         ((and (consp list)
 306               (consp (setq tem (car list)))
 307               (consp (car tem))
 308               (eq (caar tem) 'mequal)
 309               (or (eq keys '$allow_any_key)
 310                   (and (consp keys)
 311                        (member (setq sym (second tem)) keys)))
 312               (symbolp sym))
 313          (cond ((assoc sym done)
 314                 (parse_keys (cdr list) keys meval done))
 315                (t
 316                 (progv (list sym)
 317                        (list
 318                         (setq tem (if meval (meval* (third (car list)))
 319                                     (third (car list)))))
 320                        (setq done (nconc done (list (list sym tem))))
 321                        (parse_keys (cdr list) keys meval done)))))
 322         ((error "unrecognized key = ~a" (or  sym tem)))))
 323
 324
 325 (defvar *lisp-types* '(($double . long-float) ($integer . fixnum)
 326                        ($single . short-float)))
 327
 328 (defun maxima-coerce-to-xdr (v decl)
 329   (case decl
 330     ($integer v)
 331     ($double (coerce v 'long-float))
 332     ($single (coerce v 'short-float))
 333     ($external (coerce v 'fixnum))
 334     (t (cond  ((consp decl)
 335                (let* ((dim (cond ((atom (second v)) ($length v))
 336                                  (t (* ($length v) ($length (second v))))))
 337                       (ar (lisp::make-array dim :static t
 338                               :element-type (cdr (assoc (car decl)
 339                                                         *lisp-types*)))))
 340                  (or (numberp (aref ar 0)) (error "bad array type"))
 341                  (maxima-flatten ($transpose v) ar 0)
 342                  ar))
 343               (t (error "unknown type"))))))
 344
 345 (defun maxima-flatten (v ar i)
 346   (declare (fixnum i))
 347   (cond ((mbagp v)
 348          (cond ((mbagp (second v))
 349                 (dolist (w (cdr v))
 350                         (setq i (maxima-flatten w ar i))))
 351                (t (dolist (u (cdr v))
 352                           (setf (aref ar i) u)
 353                           (incf i))
 354                   i)))
 355         (t (error "bad elt"))))
 356
 357 (defun coerce-maxima-from-xdr (v elt)
 358   (setq elt (second elt))
 359   (cond ((numberp v)
 360          v)
 361         ((arrayp v)
 362          (cond ((and ($listp elt) (numberp (second elt)))
 363                 (cons '(mlist) (coerce v 'list)))
 364                (($matrixp elt)
 365                 (let ((n (length (cdr elt)))
 366                          ans)
 367                   (setq ans
 368                         (cons '($matrix)
 369                         (sloop for i below (floor (length v) n)
 370                          with ind = 0
 371                          declare (fixnum ind)
 372                          collect
 373                          (cons '(mlist)
 374                                (sloop for j below n
 375                                       collect (aref v ind)
 376                                       do (incf ind))))))
 377                   (setq ans ($transpose ans))))
 378                (t (error "unkown type"))))))
 379
 380
 381
 382
 383 ;; variables to forget from names
 384 (defvar *forget* '(#\$))
 385 (defvar *fortran-out* nil)
 386 (defun wt1 (x &aux (ch #\a))
 387   (declare (character ch))
 388   (cond ((symbolp x)
 389          (dotimes (i (length (symbol-name x)))
 390             (declare (fixnum i))
 391             (setq ch (schar (symbol-name x) i))
 392             (if (member ch *forget*) nil
 393               (princ (char-downcase ch) *fortran-out*))))
 394         ((stringp x)
 395          (princ x *fortran-out*))
 396         (t (princ x *fortran-out*))))
 397
 398 (defmacro wt (&rest l &aux lis)
 399   `(progn ,@ (dolist (v l (nreverse lis)) (push `(wt1 ,v) lis))))
 400
 401 (defmacro wt-nl (&rest l)
 402   `(progn (terpri *fortran-out*)(wt1 "      ") (wt ,@l)))
 403
 404
 405 (defun get-arg-decl (arg decls &aux size tem)
 406   (cond ((setq tem (assoc '$dimension decls))
 407          (dolist (v (cdr tem))
 408            (when (and (consp v) (eq (car v) arg))
 409              (setq size (cdr v))
 410              (return nil)))))
 411   (dolist (v decls)
 412     (dolist (u (cdr v))
 413          (cond ((atom u)
 414                 (if (eq u arg) (return-from get-arg-decl
 415                                             (cons arg
 416                                                   (if size
 417                                                       (cons (car v) size)
 418                                                     (car v))))))
 419                ((eq (car  u) arg)
 420                 (return-from  get-arg-decl(cons arg (cons (car v) (cdr u))))))))
 421   (error "undeclared arg ~a" arg))
 422
 423 (defun xdr-fun (type)
 424   (if (consp type) "array"
 425     (or (cdr (assoc type '(($integer . "int")
 426                            ($double . "double")
 427                            ($complex . "double")
 428                            ($single . "float")
 429                            ($external . "int")
 430                            )))
 431         (error "unknown type ~a" type))))
 432
 433 (defun wt-xdr (v stream)
 434   (wt-nl "CHECK(xdr_"  (xdr-fun (cdr v))"(" stream ",")
 435   (wt "&" (car v))
 436   (cond ((consp (cdr v))
 437          (wt ",&" (car v)"_length, MAX_ARRAY("(car v)"_length),")
 438          (wt-nl "  sizeof(" (get-fortran-type (second v)) "),xdr_"
 439                 (xdr-fun (cadr v)))))
 440   (wt "));"))
 441
 442
 443 (defun make-fortran-fun (fun+args decls  results name locals)
 444   (let ((file (concatenate 'string name ".c"))
 445         (args (cdr fun+args)))
 446     (or results (setq results (cdr fun+args)))
 447     (with-open-file (*fortran-out* file :direction :output)
 448       (wt "
 449 #include <stdio.h>
 450 #ifdef AIX3
 451 #include <sys/select.h>
 452 #endif
 453 #include <rpc/rpc.h>
 454 #define MAX_ARRAY(x) (x ? x :  20000)
 455 #define CHECK(x) if (!x) {fprintf(stderr,\"xdr failed\"); exit(1);}
 456 ")
 457
 458       (wt-nl)
 459       (wt "main()
 460  {XDR xdrs;
 461   int invoked=0;
 462   xdrstdio_create(&xdrs, stdin, XDR_DECODE);
 463   {")
 464       (let (res)(dolist (v args (setq args (nreverse res)))
 465                         (push (get-arg-decl v decls) res)))
 466       (dolist (v args)
 467               (cond ((atom (cdr v)) (wt-nl  (get-fortran-type(cdr v)) " " (car v)  ";"))
 468                     (t (wt-nl (get-fortran-type(cadr v)) " *" (car v) "= 0 ;")
 469                        (wt-nl "u_int " (car v) "_length = 0;"))))
 470
 471       (wt-nl "xdrstdio_create(&xdrs, stdin, XDR_DECODE);")
 472       (terpri *fortran-out*)
 473       (wt "DO_ARGS:")
 474       (dolist (v args)
 475               (wt-xdr v "&xdrs")
 476               (unless (atom (cdr v))
 477                       (if (and (numberp (third v))
 478                                (null (cdddr v)))
 479                           (wt-nl "if ("(car v)"_length != " (third v)")"
 480                                  "fprintf(stderr,\"Wrong length for " (car v)" \");")
 481                         )))
 482       (if (eq results (cdr fun+args))
 483           (wt-nl "if (invoked) exit(0);"))
 484       (dolist (v (cdr (assoc '$external decls)))
 485         (wt-nl "{extern int "v "_ext_();" v "="v "_ext_;}"))
 486
 487       ;; call the routine:
 488       (wt-nl "/* invoke the function */
 489  ")
 490       (wt-nl  (car fun+args) "_" "(")
 491       (do ((v args (cdr v)))
 492           ((null v))
 493           (or (consp (cdar v))
 494               (wt "&"))
 495           (wt (caar v))
 496           (if (cdr v) (wt ",")))
 497       (wt ");
 498    ")
 499       (wt-nl "/* write the results out */")
 500       (wt-nl "xdrstdio_create(&xdrs, stdout, XDR_ENCODE);")
 501       (cond ((equal results (cdr fun+args))
 502              (wt-nl "invoked=1 ; goto DO_ARGS;}}"))
 503             (t
 504              (dolist (v results)
 505              (wt-xdr (get-arg-decl v decls) "&xdrs"))
 506             (wt-nl "exit(0);}}"))))))
 507 #+debug
 508 (defun try-xdr (lis)
 509   (with-open-file (st "joe" :direction :output)
 510     (let ((xdrs (si::xdr-open st)))
 511       (dolist (v lis) (si::xdr-write xdrs v))))
 512   (with-open-file (st "joe" )
 513     (let ((xdrs (si::xdr-open st)))
 514       (sloop::sloop for w in lis collect (si::xdr-read xdrs w)))))
 515
 516 #|
 517 /* example of specification of routine */
 518 make_fortran(dtrdi(t,ldt,n,det,job,info),
 519                    [[ integer ,ldt,n,job,info],
 520                      [double, t(ldt,1),det(2)]],
 521                     defaults=
 522                        [ldt=length(t[1]),
 523                         det, info,
 524                         n = ldt,
 525                         job = 100 + if (t[1,2]=0) then 10  else 11],
 526                    result=[t,det,info]);
 527
 528 make_fortran(dgedi(a,lda,n,ipvt,det,work,job),
 529   [ [integer, lda,n,ipvt(n),job],
 530   [ double , a(lda,1),det(2),work(n)]],
 531   defaults=[lda=length(a),n=length(a[1]),ipvt=dgeco(a)[2],det,job=11,work],
 532   result=[a,det,ipvt])$
 533
 534 make_fortran(dgeco(a,lda,n,ipvt,rcond,z),
 535   [[  integer, lda,n,ipvt(n)],
 536    [   double , a(lda,1),z(n)],
 537     [  double, rcond]],
 538   defaults=[lda=length(a),n=length(a[1]),ipvt,rcond=0.0,z],
 539   result=[a,ipvt,rcond,z])
 540
 541
 542 >(run)
 543 Maxima 4.130 Sun Nov 12 15:59:51 CST 1989 (with enhancements by W. Schelter).
 544 (C1) make_fortran(dtrdi(t,ldt,n,det,job,info),
 545                    [[ integer ,ldt,n,job,info],
 546                      [double, t(ldt,1),det(2)]],
 547                     defaults=
 548                        [ldt=length(t[1]),
 549                         det, info,
 550                         n = ldt,
 551                         job = 100 + if (t[1,2]=0) then 10  else 11],
 552                    result=[t,det,info]);
 553
 554 (D1)                                 DTRDI
 555 (C2) dtrdi(matrix([11.00,12.0],[0.0,13.0]));
 556
 557       [ 0.090909090909090912  - 0.083916083916083919 ]
 558 (D2) [[                                              ],
 559       [         0.0            0.076923076923076927  ]
 560
 561                                                   [1.4300000000000002, 2.0], 0]
 562
 563
 564 |#
 565
 566
 567
 568
 569