Velocity Reviews - Computer Hardware Reviews

Velocity Reviews > Newsgroups > Programming > C Programming > bivariate interpolation

Reply
Thread Tools

bivariate interpolation

 
 
Charles Banas
Guest
Posts: n/a
 
      11-03-2003
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
 
Reply With Quote
 
 
 
 
Malcolm
Guest
Posts: n/a
 
      11-03-2003

"Charles Banas" <(E-Mail Removed)> wrote in message
>
> 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.


 
Reply With Quote
 
 
 
 
Sheldon Simms
Guest
Posts: n/a
 
      11-03-2003
On Mon, 03 Nov 2003 13:04:09 -0800, Charles Banas wrote:

> 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


 
Reply With Quote
 
Peter Pichler
Guest
Posts: n/a
 
      11-03-2003
"Charles Banas" <(E-Mail Removed)> wrote in message
news:(E-Mail Removed) m...

[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


 
Reply With Quote
 
CBFalconer
Guest
Posts: n/a
 
      11-03-2003
Charles Banas wrote:
>
> 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.

--
Chuck F ((E-Mail Removed)) ((E-Mail Removed))
Available for consulting/temporary embedded and systems.
<http://cbfalconer.home.att.net> USE worldnet address!


 
Reply With Quote
 
Julian V. Noble
Guest
Posts: n/a
 
      11-04-2003
Charles Banas wrote:
>
> 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,.

(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
http://www.velocityreviews.com/forums/(E-Mail Removed)
^^^^^^^^^^^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/

"Science knows only one commandment: contribute to science."
-- Bertolt Brecht, "Galileo".
 
Reply With Quote
 
nobody
Guest
Posts: n/a
 
      11-04-2003
"Peter Pichler" <(E-Mail Removed)> wrote in message
news:haBpb.1303$(E-Mail Removed)...
> "Charles Banas" <(E-Mail Removed)> wrote in message
> news:(E-Mail Removed) m...
>
> [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]


 
Reply With Quote
 
Peter Pichler
Guest
Posts: n/a
 
      11-04-2003
"nobody" <(E-Mail Removed)> wrote in message
news:MvGpb.176735$(E-Mail Removed). cable.rogers.com...
> "Peter Pichler" <(E-Mail Removed)> wrote in message
> news:haBpb.1303$(E-Mail Removed)...
> >
> > 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 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.


 
Reply With Quote
 
Charles Banas
Guest
Posts: n/a
 
      11-04-2003
Sheldon Simms <(E-Mail Removed)> wrote in message news:<(E-Mail Removed)>...
> On Mon, 03 Nov 2003 13:04:09 -0800, Charles Banas wrote:
>
> > 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.
>

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
 
Reply With Quote
 
Charles Banas
Guest
Posts: n/a
 
      11-04-2003
"Peter Pichler" <(E-Mail Removed)> wrote in message news:<haBpb.1303$(E-Mail Removed)>...
> "Charles Banas" <(E-Mail Removed)> wrote in message
> news:(E-Mail Removed) m...
>
> [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.

> > 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
>

<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
 
Reply With Quote
 
 
 
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are Off


Similar Threads
Thread Thread Starter Forum Replies Last Post
Linear interpolation in vhdl df84077@gmail.com VHDL 8 10-07-2011 08:38 AM
2D - bivariate random number generatin John Java 10 06-15-2006 10:22 AM
How to stop interpolation of Escape character in variable? \Rob\ ASP .Net 8 02-24-2006 08:52 AM
Opening files with interpolation Mat W Perl 3 05-21-2004 08:40 AM
Your votes please ( drawImage bad interpolation bug id 4950176) Michele Puccini Java 0 11-09-2003 11:06 AM



Advertisments