bivariate interpolation

[Charles Banas]

I'm not sure if this is the right place to ask about this, but I've
seen several posts in the past regarding Akima's Bivariate
Interpolations routines, and i'm wondering if someone can give me some
ideas or PD code I can put to use right away.

At the moment, I'm maintaining a contour calculation and plotting
program for radio wave propagation studies. This program uses a
number of routines written by and for the FCC in Fortran in 1976 for
their Univac machine. It was more recently used on VAX machines, and
is now used on Sun Sparc stations. The code I'm using is directly
from FCC/OCE report RS76-01, January 1976. It has numerical tables
calculated from an unspecified bivariate polynomial, which are
interpolated using Akima's Bivariate Interpolation for Smooth Surface
Fitting algorithm.

Unfortunately, due to the nature of the code, I've been forced to
rewrite the code form its original Fortran V to much cleaner C.
(Actually, I've been using ISO C++, but that's beside the point.) My
problem is this: Akima's routine is extremely cryptic and
uncommented, impossible to compile with the tools I have available
(mostly stemming from the complexity of the EQUIVALENCE statements he
used), and so, is very hard to rewrite into C.

I tried just diving into a rewrite, but the biggest issue I ran across
was the way the DIMENSIONed arrays are accessed:

The routine has numerous 2-dimensional arrays that are being accessed
as one-dimensional arrays and then aliased with simple variables.
Since I'm very unfamiliar with Fortran's underpinnings and quirks, I'm
having a hard time with the code. And, unfortunately, because of the
turn this project took, I need it rewritten in C.

I have tried using f2c to translate the routine directly into C, but
I'm quite uncomfortable with how messy the code is, and I honestly
don't trust f2c all that much. (Especially since f2c returned with
well over 40 warnings.)

The name of the routine is ITPLBV, and this version has only 10
parameters. I'm wondering if anyone here knows of an existing
(non-f2c) C version of the routine - possibly one written by Akima
himself? Or maybe someone here is familiar enough with Fortran to
give me some insight on how to treat these 2-dimensional arrays in C?

I'd appreciate anything anyone has to say.

-- Charles Banas

 

[Malcolm]

I'm not sure if this is the right place to ask about this

Try comp.programming, however your post does take a C language turn.

Akima's routine is extremely cryptic and uncommented, impossible to
compile with the tools I have available (mostly stemming from the
complexity of the EQUIVALENCE statements he
used), and so, is very hard to rewrite into C.

I don't know much Fortran. One problem you may have is that in C 2d arrays
are of fixed dimensions. This can be fixed by using malloc and accessing by
array[y * width + x];

Another problem is that C is call by value. This can be fixed by liberal use
of pointers, though it will make the code very hard to read.

The thing to do is to use these techniques to write Fortran-like C. It will
be horrible but it should compile (BTW did you test the routine spat out by
your Fortran-to-C converter?).
Then rewrite it as more idiomatic C, keeping a working version for test
purposes.

 

[Sheldon Simms]

I'm not sure if this is the right place to ask about this,

It's not.

My problem is this: Akima's routine is extremely cryptic and
uncommented, impossible to compile with the tools I have available
(mostly stemming from the complexity of the EQUIVALENCE statements he
used), and so, is very hard to rewrite into C.

The equivalence statements are old style, they subscript multi-dimension
arrays with a single subscript. This was required by very old Fortrans
but g77 spews errors. After rewriting all of the subscripts in the
equivalence statements, g77 compiled the code without errors. Whether
it works is another matter.

This can be rewritten in C but given that it's easy to make it
compilable in Fortran, why not just compile it in fortran and
link it into your C++ progam?

-Sheldon

followups set to comp.lang.fortran

 

[Peter Pichler]

[lenghty and irrelevant background about converting a Fortran program
to C snipped.]

The routine has numerous 2-dimensional arrays that are being accessed
as one-dimensional arrays and then aliased with simple variables.
Since I'm very unfamiliar with Fortran's underpinnings and quirks, I'm
having a hard time with the code. And, unfortunately, because of the
turn this project took, I need it rewritten in C.

I have tried using f2c to translate the routine directly into C, but
I'm quite uncomfortable with how messy the code is, and I honestly
don't trust f2c all that much. (Especially since f2c returned with
well over 40 warnings.)

Charles,
I am not at all surprised about the warnings. Given the way the program
treats 2D arrays, there are almost certainly other "interesting", ehm,
"features"

The name of the routine is ITPLBV, and this version has only 10
parameters. I'm wondering if anyone here knows of an existing
(non-f2c) C version of the routine - possibly one written by Akima
himself? Or maybe someone here is familiar enough with Fortran to
give me some insight on how to treat these 2-dimensional arrays in C?

I'd appreciate anything anyone has to say.

I'm afraid I know very little about Fortran, but in C multidimensional
arrays are emulated by arrays of arrays. Let's say we have a declaration

double xy[4][3];

This is an array of 4 elements, each of them an array of 3 doubles. It
could be rewritten as

mytype y[4];

where y is an array of 4 elements of type mytype. You can draw it like
this:

+------+------+------+------+
| y[0] | y[1] | y[2] | y[3] | (you know that C arrays are 0-based, right?)
+------+------+------+------+

Because of C pointer arithmetic, arrays must be contiguous. If you have a
pointer to mytype and point it to y[0], incrementing that pointer must
yield a pointer to y[1].

Now if you remember, each of these y is an array itself. Let's draw another
picture:

+--------------+--------------+--------------+--------------+
| y[0] | y[1] | y[2] | y[3] |
+----+----+----+----+----+----+----+----+----+----+----+----+
|x[0]|x[1]|x[2]|x[0]|x[1]|x[2]|x[0]|x[1]|x[2]|x[0]|x[1]|x[2]|
+----+----+----+----+----+----+----+----+----+----+----+----+

As you can see, inside each y there are 3 x, because each y is an array of
3 x.

Now because there is no checking for array bounds in C, it is very easy to
skip the fence and peek into the neighbour's garden by using wrong a index.
For example, y[1]x[1] (just an illustration, NOT a valid C syntax!) could
be accessed as y[0]x[4] (mind you, the valid C syntax would be xy[1][1] and
xy[0][4] respectively, given my first declaration).

This is not as unsafe as it may seem, if you know what you are doing, but
it does lead to a very n#unreadable code.

It appears that your Fortran code uses a similar technique, so I guess it
must work in Fortran too. However, I heard that the indexes (indices?) are
reversed in Fortran, i.e. the first one is the one changing the fastest,
not the ast one as in C.

If I were to rewrite a piece of code from one language into another, I
would first try to understand what the code does and how it does it. It
seems from your post that you may have understand what, but have problems
understanding how. With that, we cannot help you, try a Fortran newsgroup
or Google for bivariate interpolation.

Or hire an experienced C programmer ;-)

 

[CBFalconer]

I'm not sure if this is the right place to ask about this, but
I've seen several posts in the past regarding Akima's Bivariate
Interpolations routines, and i'm wondering if someone can give
me some ideas or PD code I can put to use right away.
.... snip ...

Unfortunately, due to the nature of the code, I've been forced
to rewrite the code form its original Fortran V to much cleaner
C. (Actually, I've been using ISO C++, but that's beside the
point.) My problem is this: Akima's routine is extremely
cryptic and uncommented, impossible to compile with the tools I
have available (mostly stemming from the complexity of the
EQUIVALENCE statements he used), and so, is very hard to
rewrite into C.
.... snip ...

I'd appreciate anything anyone has to say.

This is neither fish nor fowl. It is not as far off-topic as it
first appears, because the fundamental question is how to
translate bad Fortran into C. Is is not really an algorithmic
question, which would belong on comp.programming.

Given that you are translating shaky concepts I think you should
use a language with very strong typing, either Pascal or Ada.
Don't worry about efficiency, just get the algorithms right.

 

[Julian V. Noble]

I'm not sure if this is the right place to ask about this, but I've
seen several posts in the past regarding Akima's Bivariate
Interpolations routines, and i'm wondering if someone can give me some
ideas or PD code I can put to use right away.

At the moment, I'm maintaining a contour calculation and plotting
program for radio wave propagation studies. This program uses a
number of routines written by and for the FCC in Fortran in 1976 for
their Univac machine. It was more recently used on VAX machines, and
is now used on Sun Sparc stations. The code I'm using is directly
from FCC/OCE report RS76-01, January 1976. It has numerical tables
calculated from an unspecified bivariate polynomial, which are
interpolated using Akima's Bivariate Interpolation for Smooth Surface
Fitting algorithm.

Unfortunately, due to the nature of the code, I've been forced to
rewrite the code form its original Fortran V to much cleaner C.
(Actually, I've been using ISO C++, but that's beside the point.) My
problem is this: Akima's routine is extremely cryptic and
uncommented, impossible to compile with the tools I have available
(mostly stemming from the complexity of the EQUIVALENCE statements he
used), and so, is very hard to rewrite into C.

I tried just diving into a rewrite, but the biggest issue I ran across
was the way the DIMENSIONed arrays are accessed:

The routine has numerous 2-dimensional arrays that are being accessed
as one-dimensional arrays and then aliased with simple variables.
Since I'm very unfamiliar with Fortran's underpinnings and quirks, I'm
having a hard time with the code. And, unfortunately, because of the
turn this project took, I need it rewritten in C.

I have tried using f2c to translate the routine directly into C, but
I'm quite uncomfortable with how messy the code is, and I honestly
don't trust f2c all that much. (Especially since f2c returned with
well over 40 warnings.)

The name of the routine is ITPLBV, and this version has only 10
parameters. I'm wondering if anyone here knows of an existing
(non-f2c) C version of the routine - possibly one written by Akima
himself? Or maybe someone here is familiar enough with Fortran to
give me some insight on how to treat these 2-dimensional arrays in C?

I'd appreciate anything anyone has to say.

-- Charles Banas

Let me take a swing at this for you. First of all, Fortran stores
2-dimensional arrays by columns, not by rows. This is a frequent
problem when trying to call a compiled Fortran routine from C or
other language that stores row-wise.

The point of EQUIVALENCE in Fortran was to tell the compiler to
put two arrays (or chunks of arrays) in the same memory locations:
thus you might have

REAL MAT(10,10)
REAL A(100)
EQUIVALENCE M(1,1), A(1)

Old Fortran did not have dynamically dimensioned arrays. That is,
you could not tell it how much space to allocate at run time. I
think more modern versions let you do it.


Finally, as far as I know, bivariate interpolation is no harder than
univariate interpolation (at least in concept). Let's say you have
a table of 121 points, 11 in the x-direction and 11 in the y-direction.

And suppose, to make things simple, x runs from 0.0 to 10.0 in steps of
1.0, and the same for y. Then MAT(1,1) is the point x=0.0, y=0.0
and MAT(11,11) is the point x=10.0, y=10.0 . Similarly the value at
x=3.0, y=7.0 would be stored in MAT(4,8).

(I am doing it this way because old Fortran used 1-based indexing rather
than 0-based as in C, to dereference arrays. If it were a C array the
indices would run from 0 to 10 rather than 1 to 11.)

Now, to interpolate linearly in Fortran you could write

K = INT(X) *** integer part of x ***
L = INT(Y) *** etc ***
U = MAT(K,L) + (X-K)*( MAT(K+1,L) - MAT(K,L) )
1 + (Y-L)*(MAT(K,L+1)-MAT(K,L))

(note the 1 in column 6 of the punched card to indicate continuation ;-)

To translate all this into C isn't that hard. C has automatic type
conversion so you only have to say

int k,l;
double x,y;

k = x;
l = y;

to get the integer parts.

Don't forget that in general you will need a scaling factor to map the
range of x or y into an integer range. That is, if you have 101 points
from x=0.0 to 1.0 at intervals of 0.01, you will need to multiply by
100.

It is possible to use higher-order interpolation schemes--see Abramowitz
& Stegun (available as a Dover PB --you should own it, it's very useful).

I never heard of Akima or his bivariate interpolation. But it must be some-
thing like the above.

Have a look at the lecture notes available at the site

http://www.phys.virginia.edu/classes/551.jvn.fall01/

if you need further interpolation ideas.


--
Julian V. Noble
Professor Emeritus of Physics
(e-mail address removed)
^^^^^^^^^^^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/

"Science knows only one commandment: contribute to science."
-- Bertolt Brecht, "Galileo".

 

[nobody]

[lenghty and irrelevant background about converting a Fortran program
to C snipped.]
[more snip]

Because of C pointer arithmetic, arrays must be contiguous. If you have a

Actually, I believe that it is the other way around. Arrays are (*must be*)
contiguous by definition /cause/, and we have pointer arithmetic because
of this /effect/.

<quote>
6.2.5 Types
....
19 Any number of derived types can be constructed from the object, function,
and
incomplete types, as follows:
- An array type describes a contiguously allocated nonempty set of objects
with a
particular member object type, called the element type.33) ...
</quote>

[snip]

 

[Peter Pichler]

a

Actually, I believe that it is the other way around. Arrays are (*must be*)
contiguous by definition /cause/, and we have pointer arithmetic because
of this /effect/.

<Quote from the Standard snipped>

Err, yes, indeed. I was a bit ambiguous, I guess. What I had in mind (and
failed to express clearly) was that in theory, a programming language may
implement an array in such a way that the memory it occupies does not have
to be contiguous. Indexing would do the magic behind the programmer's back.
But this would be more difficult, runtime consuming and require array bounds
checking. I have no idea whether any implementation of any language bothers
doing something like that, but C does not. And yes, that makes pointer
arithmetics easier.

 

[Charles Banas]

It's not.


The equivalence statements are old style, they subscript multi-dimension
arrays with a single subscript. This was required by very old Fortrans
but g77 spews errors. After rewriting all of the subscripts in the
equivalence statements, g77 compiled the code without errors. Whether
it works is another matter.

Which is the core of the issue for me. I'm basically moving from a MS
VC++ / Intel Fortran environment to an all-GCC environment, mainly
because I've been unable to get the other tools I need.

This can be rewritten in C but given that it's easy to make it
compilable in Fortran, why not just compile it in fortran and
link it into your C++ progam?

I considered that. But because I'm so unfamiliar with Fortran arrays,
I decided it would probably be easier to do a rewrite. Even after
reading through my Fortran books, I came away with a less complete
understanding of the way array elements are ordered than before.

Besides, I am a huge fan of clean code, and I wanted code that
wouldn't be grossly over-populated with pointers and ugly constructs.
I rewrote the original calling function (a freespace calculator) using
very clean-looking C arrays, and that is clearly imcompatible with
Fortran arrays.

Hence my problem.

-Sheldon

followups set to comp.lang.fortran

Sorry, the reason I posted to comp.lang.c in the first place was
because I was hoping someone knew of a C version of the routine
already available.

I apologize in advance for making this a cross-post.

-- Charles Banas

 

[Charles Banas]

[lenghty and irrelevant background about converting a Fortran program
to C snipped.]

The routine has numerous 2-dimensional arrays that are being accessed
as one-dimensional arrays and then aliased with simple variables.
Since I'm very unfamiliar with Fortran's underpinnings and quirks, I'm
having a hard time with the code. And, unfortunately, because of the
turn this project took, I need it rewritten in C.

I have tried using f2c to translate the routine directly into C, but
I'm quite uncomfortable with how messy the code is, and I honestly
don't trust f2c all that much. (Especially since f2c returned with
well over 40 warnings.)

Charles,
I am not at all surprised about the warnings. Given the way the program
treats 2D arrays, there are almost certainly other "interesting", ehm,
"features"

Which scare the jeepers out of me. I'd rather rewrite something than
make it unmaintainable with a machine translator.

I'm afraid I know very little about Fortran, but in C multidimensional
arrays are emulated by arrays of arrays. Let's say we have a declaration

<snip arrays tutorial>

Sorry, I've already got that down.

It appears that your Fortran code uses a similar technique, so I guess it
must work in Fortran too. However, I heard that the indexes (indices?) are
reversed in Fortran, i.e. the first one is the one changing the fastest,
not the ast one as in C.

The difference is that (in my understanding) Fortran arrays are
linearly addressed, unlike the array-of-pointers in C.

If I were to rewrite a piece of code from one language into another, I
would first try to understand what the code does and how it does it. It
seems from your post that you may have understand what, but have problems
understanding how. With that, we cannot help you, try a Fortran newsgroup
or Google for bivariate interpolation.

I probably should have posted to comp.lang.fortran. Sorry for the
hassle.

And yes, I did google for Bivariate Interpolation, but the only
information I can find is the ACM article (which requires an expensive
membership).

Or hire an experienced C programmer ;-)

Thing is, that was /me/. Small company doing big things, and I'm one
of very few local programmers.

Thank you anyway.

-- Charles Banas

 

[Charles Banas]

This is neither fish nor fowl. It is not as far off-topic as it
first appears, because the fundamental question is how to
translate bad Fortran into C. Is is not really an algorithmic
question, which would belong on comp.programming.

I'm *really* glad you think so. That was my biggest concern about
posting.

Given that you are translating shaky concepts I think you should
use a language with very strong typing, either Pascal or Ada.
Don't worry about efficiency, just get the algorithms right.

Which is why I posted here. I'm using ISO C++ because of two main
facts: 1. it's rather strongly typed, and 2. I'm familiar/comfortable
with it.

 

[Charles Banas]

Let me take a swing at this for you. First of all, Fortran stores
2-dimensional arrays by columns, not by rows. This is a frequent
problem when trying to call a compiled Fortran routine from C or
other language that stores row-wise.

You just pinpointed my main point of confusion. I know this, yet I
have a good deal of trouble understanding.

The point of EQUIVALENCE in Fortran was to tell the compiler to
put two arrays (or chunks of arrays) in the same memory locations:
thus you might have

REAL MAT(10,10)
REAL A(100)
EQUIVALENCE M(1,1), A(1)

Ready for a scare? Akima's EQUIVALENCE statements are 13 cards long.
Like this:

DIMENSION ZA(5,2), ....
EQUIVALENCE ( Z3A1, ZA(1) ), ( Z3A2, ZA(2) ), ( Z3A3, ZA(3) ),
& ( Z3A4, ZA(4) ), ( Z3A5, ZA(5) ), ( Z4A1, ZA(6) ),
etc.

EQUIVALENCES like this are confounding to me.

Old Fortran did not have dynamically dimensioned arrays. That is,
you could not tell it how much space to allocate at run time. I
think more modern versions let you do it.


Finally, as far as I know, bivariate interpolation is no harder than
univariate interpolation (at least in concept). Let's say you have
a table of 121 points, 11 in the x-direction and 11 in the y-direction.

<snip useful information>

You wouldn't know that by looking at Akima's code. It's literally 20
pages of quite confusing algorithmic code. I haven't had time to
analyze it.

It is possible to use higher-order interpolation schemes--see Abramowitz
& Stegun (available as a Dover PB --you should own it, it's very useful).

Could you point me to a location I can acquire it? It does sound as
if it would be useful for other projects as well.

I never heard of Akima or his bivariate interpolation. But it must be some-
thing like the above.

I haven't been able to find much information, but it appears that
Akima used a triangle-based algorithm for the interpolation. It
doesn't look one bit like any interpolation I've seen or done before.

Have a look at the lecture notes available at the site

http://www.phys.virginia.edu/classes/551.jvn.fall01/

if you need further interpolation ideas.

I will file this away for future reference. Thank you for the info.

--
Julian V. Noble
Professor Emeritus of Physics
(e-mail address removed)
^^^^^^^^^^^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/

"Science knows only one commandment: contribute to science."
-- Bertolt Brecht, "Galileo".

-- Charles Banas

 

[Charles Banas]

Try comp.programming, however your post does take a C language turn.

That's why I'm here and not at c.p.

Akima's routine is extremely cryptic and uncommented, impossible to
compile with the tools I have available (mostly stemming from the
complexity of the EQUIVALENCE statements he
used), and so, is very hard to rewrite into C.

I don't know much Fortran. One problem you may have is that in C 2d arrays
are of fixed dimensions. This can be fixed by using malloc and accessing by
array[y * width + x];

Sorry, I'm having a little trouble following this. I'm not quite sure
what you mean.

Another problem is that C is call by value. This can be fixed by liberal use
of pointers, though it will make the code very hard to read.

Which is one of the reasons for my rewrite. I'd rather /not/ have to
muck things up with references and pointers. I have a hard enough
time following the thread about casting malloc() calls. (Fun
read, BTW. Almost a flame war, but not quite.

The thing to do is to use these techniques to write Fortran-like C. It will
be horrible but it should compile (BTW did you test the routine spat out by
your Fortran-to-C converter?).
Then rewrite it as more idiomatic C, keeping a working version for test
purposes.

I will try the f2c routine as I progress and use it for referencial
behavior. It's just butt-ugly.

As I mentioned, I have been doing this in ISO C++, mostly for the type
checking. But after reading the malloc() casting arguments, what
c.l.c considers "Idomatic C" honestly scares me. Such Holy Wars over
casting. Wow.

However, there is one thing I think I should have mentioned. I'm
considering using sed script to remove all the EQUIVALENCEs. I
wonder, is that safe to do, or am I playing with fire now (since my
sed knowledge is limited)?

-- Charles Banas

 

[Alan Balmer]

Which is why I posted here. I'm using ISO C++ because of two main
facts: 1. it's rather strongly typed, and 2. I'm familiar/comfortable
with it.

After seeing this several times, I went back to your original post,
and found: "(Actually, I've been using ISO C++, but that's beside the
point.)"

Perhaps you should consider that it might *not* be beside the point. C
and C++ are different languages, and their practitioners often have
different approaches to problems. In any case, it wouldn't hurt to
discuss this in comp.lang.c++.

I'm peculiar in that I rather enjoy deciphering and cleaning up other
people's code, and I have some pretty good tools for doing it, so if I
had your problem, I'd probably take a shot at making the f2c output
palatable.

 

[Charles Banas]

After seeing this several times, I went back to your original post,
and found: "(Actually, I've been using ISO C++, but that's beside the
point.)"

Perhaps you should consider that it might *not* be beside the point. C
and C++ are different languages, and their practitioners often have
different approaches to problems. In any case, it wouldn't hurt to
discuss this in comp.lang.c++.

I posted here because I felt a C discussion would be more relevant
than one in C++. In any case, you may be right.

I'm peculiar in that I rather enjoy deciphering and cleaning up other
people's code, and I have some pretty good tools for doing it, so if I
had your problem, I'd probably take a shot at making the f2c output
palatable.

Which is something I just may do. I'm sure you're well aware of how
ugly f2c's output is, though.

Thank you.

-- Charles Banas