Re: need performance tips

----- Original Message ----- 
From: "Richard Fateman" <[email protected]>
To: "'Andrew K. Wolven'" <[email protected]>; <[email protected]>
Sent: Saturday, September 23, 2006 5:43 PM
Subject: RE: need performance tips


> You seem to be writing a fortran program even though you are using lisp. 
> You
> could just transform the fortran into lisp by using f2cl.

The computer would then know all about programming spline curves, and I 
would only know how to type in some fortran code and call f2cl on it which 
is completely useless as preparation for the exam.

>
> Since your programs have no comments, it is hard to see what you are 
> trying
> to do, but why is "derivative" so long??

First of all, it is a misomer.  It is a vector containing the second 
derivatives of a peicewise polynomial at it's knots.  As you can see, the 
vector T2 is set to zero at it's end points (natural spline).  Then to kick 
off the loading of the vector, one starts by working backwards by first 
computing the second to last second derivative.  Then the other second 
derivatives are computed by using the coefficients of the polynomial at that 
point and the value of the second derivative previously computed (the i+1th: 
we are going backwards).  The coefficient vectors are computed in the 
previous two loops.  The code for coming up with these coefficients comes 
from the algebraic solution of solving a system of equations which describe 
the spline, which in turn was written and optimized in fortran, and then 
manually translated into CL.  The fortran name for the subroutine was 
"derivative".  Unfortunately I got an A in the course and never got back to 
completely dissecting, lispifiying and reoptimizing the algorithm for 
computing the coefficients.  I just know it computed the right answers (once 
I found the typo in my test dataset).

I'm going to remimplement the program from scratch from the algebra and then 
optimize it.  I will probably do the parametric definition of instead of 
this one which is a function of x.  The parametric definition does not 
prohibit unique y values for a particular x value, so you can have loops and 
things, while this implementation is designed for curve fitting an (x_i, 
y_i) dataset.

It's getting late.  Why is a mathematician asking me about splines.

AKW



>
>
>> -----Original Message-----
>> From: Andrew K. Wolven [mailto:[email protected]]
>> Sent: Saturday, September 23, 2006 1:23 PM
>> To: [email protected]
>> Subject: need performance tips
>>
>> Hi all,
>> I am preparing to do a mini-lisp primer for a group of CAD students in
>> preparation for them to see a GDL presentation by David Cooper.
>>
>> Here is the low down:
>> I have chosen to do the cubic spline as an example.  It's is going to be 
>> a
>> comparison against Fortran which most of them have taken.  Therefore I
>> will
>> be pitting something like gcc against ACL.  *You know* somebody is going
>> to
>> ask which language is faster in execution.  This is going to be tough.
>>
>> I will do a "get it done correctly" example first, and then show a
>> pre-prepared optimized version.
>>
>> What can anyone tell me (or point me in the right direction) about 
>> getting
>> optimal numerical performance where I will be using arrays and iteration
>> to
>> do this.
>>
>> Here is some code that worked last time I checked: It's basically a
>> hand-translation of a Fortran program, so I know it might be a little
>> ugly.
>> I basically just wanted to get it done to turn in and get my grade and
>> move
>> on.  I do know somethings, like obviously I can have those arrays in
>> derivative function pre-made and reused, so long as I am careful about 
>> it.
>> Oh yeah, it's missing a macro vref which is something like (defmacro vref
>> (array index) `(aref ,array (1- ,index)))  --> which I probably should 
>> rid
>> myself of for the optimized version and perhaps use svref?
>> What about declarations?  If anyone happens to want to share a better
>> algorithm, I would also be very grateful.
>>
>> Thanks in advance all,
>> (I hope this is correct code!)
>> AKW
>>
>> (in-package :me428)
>>
>> (defun V-example-6.4 ()
>>   (make-array
>>    15
>>    :initial-contents
>>    (list 0.397 0.798 1.203 1.611 2.022 2.436 2.85 3.266 3.681 4.095 4.508
>> 4.919 5.327 5.733 6.137)))
>>
>> (defun Temp-example-6.4 ()
>>   (make-array
>>    15
>>    :initial-contents
>>    (list 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150)))
>>
>> (defun derivative (m V Temp)
>>   (let ((T2 (make-array m))
>>  (A (make-array m))
>>  (B (make-array m))
>>  (C (make-array m))
>>  (D (make-array m)))
>>     (setf (vref C 1) (- (vref V 2) (vref V 1)))
>>     (loop for i from 2 to (1- m)
>>  do (setf (vref A i) (- (vref V i) (vref V (1- i)))
>>    (vref B i) (* 2.0 (- (vref V (1+ i)) (vref V (1- i))))
>>    (vref C i) (- (vref V (1+ i)) (vref V i))
>>    (vref D i) (* 6.0 (- (/ (- (vref Temp (1+ i)) (vref Temp i)) (vref C
>> i))
>>           (/ (- (vref Temp i) (vref Temp (1- i))) (vref A i))))))
>>     (loop for i from 3 to (1- m)
>>  do
>>    (setf (vref B i) (- (vref B i)
>>          (/ (* (vref A i)
>>         (vref C (1- i)))
>>      (vref B (1- i))))
>>   (vref D i) (- (vref D i)
>>          (/ (* (vref A i)
>>         (vref D (1- i)))
>>      (vref B (1- i))))))
>>     (setf (vref T2 1) 0.0
>>    (vref T2 m) 0.0
>>    (vref T2 (1- m)) (/ (vref D (1- m)) (vref B (1- m))))
>>     (let (IN)
>>       (loop for i from 2 to (- m 2)
>>    do (setf IN (- m i)
>>      (vref T2 IN) (/ (- (vref D IN)
>>           (* (vref C IN)
>>       (vref T2 (1+ IN))))
>>        (vref B IN)))))
>>     T2))
>>
>> (defun spline-interpolate (V Temp vp)
>>   (let* ((m (length V))
>>   (T2 (derivative m V Temp))
>>   TP)
>>     (loop for i from 1 to (1- m)
>>  do (when (<= vp (vref V (1+ i)))
>>       (let ((S1 (- (vref V (1+ i)) (vref V i)))
>>      (S2 (- vp (vref V i)))
>>      (S3 (- (vref V (1+ i)) vp)))
>>         (setf TP (+ (/ (* (vref T2 i) S3 (- (/ (expt S3 2) S1) S1)) 6.0)
>>       (/ (* (vref T2 (1+ i)) S2 (- (/ (expt S2 2) S1) S1)) 6.0)
>>       (/ (* (vref Temp i) S3) S1)
>>       (/ (* (vref Temp (1+ i)) S2) S1)))
>>         (return))))
>>     TP))
>>
>
>