F2CL ODEPACK problems.

Hi


I think I got another problem with odepack-lsode.

I wrote the following routine, which, AFAIU, is a translation of the
Fortran example given in the main opkdmain.f file from netlib (sorry if
it is long.)


> (defun kinetic-example ()
>   (flet ((fex (neq time y ydot)
>            (declare (ignore neq time))
>            (setf (aref ydot 0) #I(-0.04d0 * y[0] + 1.0D4 * y[1] * y[2])
>                  (aref ydot 1) #I(3.0D7 * y[1] * y[1]))
>            (setf (aref ydot 2) #I(-ydot[0] - ydot[1]))
>            )
> 	 (jex (neq time y ml mu pd nrpd)
>            (declare (ignore neq time ml mu nrpd))
>            (setf (aref pd 0 0) -0.04D0)
>            (setf (aref pd 0 1) #I(1.0D4 * y[2]))
>            (setf (aref pd 0 2) #I(1.0D4 * y[1]))
>            (setf (aref pd 1 0) 0.04D0)
>            (setf (aref pd 1 2) (- (aref pd 0 2)))
>            (setf (aref pd 2 1) #I(6.0D7 * y[1]))
>            (setf (aref pd 1 1) #I(-pd[0, 1] - pd[2, 1]))
>            )
> 	 )
>     (let ((neq (make-array 1 :initial-element 3 :element-type 'f2cl-lib:integer4))
>           ;; (neq 3) ; Fortran is fun!
> 	  (time 0.D0)
> 	  (tout 0.4D0)
> 	  (itol 2)
>           (rtol (make-array 1 :initial-element 1.0D-4 :element-type 'double-float))
> 	  ;; (rtol 1.0D-4) ; Fortran is fun!
> 	  (itask 1)
> 	  (istate 1)
> 	  (iopt 1)
> 	  (lrw 58)
> 	  (liw 23)
> 	  (mf 21)
> 	  ;;
> 	  (y (make-array 3
> 			 :element-type 'double-float
> 			 :initial-contents '(1.0D0 0.0D0 0.0D0)))
> 	  (atol (make-array 3
> 			    :element-type 'double-float
> 			    :initial-contents '(1.0D-6 1.0D-10 1.0D-6)))
> 
> 	  (iwork (make-array 23
> 			     :initial-element 0
> 			     :element-type 'fixnum))
> 	  (rwork (make-array 58
> 			     :initial-element 0.0D0
> 			     :element-type 'double-float))
> 	  )
>       (loop for iout from 1 upto 12
> 	    for t-out = tout then (* tout 10.0D0)
> 	    do (multiple-value-bind (fex
> 			             neq
> 			             y
> 			             time
> 			             tout
> 			             itol
> 			             rtol
> 			             atol
> 			             itask
> 			             istate
> 			             iopt
> 			             rwork
> 			             lrw
> 			             iwork
> 			             liw
> 			             jex
> 			             mf)
> 		   (dlsode #'fex
> 			   neq
> 			   y
> 			   time
> 			   t-out
> 			   itol
> 			   rtol
> 			   atol
> 			   itask
> 			   istate
> 			   iopt
> 			   rwork
> 			   lrw
> 			   iwork
> 			   liw
> 			   #'jex
> 			   mf)
>                  (declare (ignore fex
>                                   neq
>                                   ;; y
>                                   ;; time
>                                   tout
>                                   itol
>                                   rtol
>                                   atol
>                                   itask
>                                   ;; istate
>                                   iopt
>                                   rwork
>                                   lrw
>                                   iwork
>                                   liw
>                                   jex
>                                   mf))
> 		 (format t "At t = ~F   Y = ~S.~%" time y)
> 		 (when (minusp istate)
> 		   (error "Error halt. ISTATE = ~D." istate)))
> 	    finally (format t "~%N. steps = ~D, N. f-s = ~D, N. J-s = ~D.~%"
> 			    (aref iwork 11)
> 			    (aref iwork 12)
> 			    (aref iwork 13))
> 	    ))
>     ))

When I try to execute the above in LW (I'd bet the same would happen in
other lisps) I get the following error


Error: The array #(0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 9900.990099009901
1.0E10 1000000.0 -0.039995298965865594 0.0 0.039995298965865594 0.0 0.0
0.0) requires 1 subscript but 2 were given.
   1 (abort) Return to level 0.
   2 Return to top loop level 0.

Type :b for backtrace, :c <option number> to proceed,  or :? for other
options


I tracked down this error to the DPREPJ routine in the call to the JEX
local function, which expects a 2D array, but it seems to get a 1D array
slice.

Is this the correct diagnosis?  Is this a Fortran trick to pass around
things that can be seen of different shapes at different times?

thanks
--
Marco




-------------------------------------------------------
SF.Net email is sponsored by Shop4tech.com-Lowest price on Blank Media
100pk Sonic DVD-R 4x for only $29 -100pk Sonic DVD+R for only $33
Save 50% off Retail on Ink & Toner - Free Shipping and Free Gift.
http://www.shop4tech.com/z/Inkjet_Cartridges/9_108_r285