Complex : mixing C and fortran
relaxmike
I am working on an interface between a C code and a Lapack routine,
which computes the generalized eigenvalues of a couple of complex
matrices (zggev).
What we have here is a C internal representation of complex arrays,
where we have a block of double values for the real
part and a block of double values for the imaginary part.
Was this a good choice or not, I don't know, but this is the current
algorithm.
I see that there is a process to convert from/to the internal by-
block
representation used in the C code and from/to the representation used
in
the fortran 77 routine. The Lapack routine is based on the COMPLEX*16
statement.
The map from C to fortran is based on a C struct :
struct complex {
double x, y;
};
I understand from the code I have under the eyes that the fortran 77
array of complex stores the values with an alternation, that is,
R1, I1, R2, I2, etc... where R stands for the real part and I for the
imaginary
part.
The algorithm written in C then rearrange the data to reorder
the real and imaginary parts by block.
What I would like to ask to fortran experts is the following :
does the C struct above will allways map to the fortran COMPLEX*16
data type, or is this compiler-dependent ?
What would happen if I had a fortran 90 compiler instead ?
Is this written in the fortran standard ?
Regards,
Michaël
Jugoslav Dujic
<snip>
> The map from C to fortran is based on a C struct :
>
> struct complex {
> double x, y;
> };
>
>
> What I would like to ask to fortran experts is the following :
> does the C struct above will allways map to the fortran COMPLEX*16
> data type, or is this compiler-dependent ?
> What would happen if I had a fortran 90 compiler instead ?
> Is this written in the fortran standard ?
The briefest answer would be, IMHO:
This is compiler-dependent in theory. In practice, the above
struct will always map to the Fortran complex*16 for every
plausible combination of compilers.
--
Jugoslav
www.xeffort.com
Please reply to the newsgroup.
You can find my real e-mail on my home page above.
Tobias Burnus
> relaxmike wrote:
> > The map from C to fortran is based on a C struct :
> >
> > struct complex {
> > double x, y;
> > };
> >
> > does the C struct above will allways map to the fortran COMPLEX*16
> > data type, or is this compiler-dependent?
No, it will definitely not map to always Fortran's COMPLEX*16 data
type. In additional, whether it works is also platform dependent.
> > What would happen if I had a fortran 90 compiler instead?
> > Is this written in the fortran standard ?
Well, on the C side, one can use the C99 datatype:
double _Complex z;
and on the Fortran side using a Fortran 2003 feature:
USE ISO_C_Binding
complex(c_double_complex) :: z
> The briefest answer would be, IMHO:
>
> This is compiler-dependent in theory. In practice, the above
> struct will always map to the Fortran complex*16 for every
> plausible combination of compilers.
Using a struct definitely breaks with gfortran on some systems. We had
such a construct in our testsuite and it failed on some systems. I
think among the systems was an IBM s390; on x86/x86-64 the test case
worked, though. Note: For gfortran, using "complex*16" is equivalent
to "complex(c_double_complex)" thus the problem is (for gfortran) only
on the C side and has something to do with alignment/padding. [Maybe
the problem was only for "float" and not for "double", I forgot.]
I think the f2c compiler and maybe also g77 were using (internally)
such a struct, i.e. in that case "double _Complex" might not work
properly.
Tobias
JB
> I am working on an interface between a C code and a Lapack routine,
Why not use one of the available C/C++ Lapack wrappers rather than
reinventing the wheel?
> What we have here is a C internal representation of complex arrays,
> where we have a block of double values for the real
> part and a block of double values for the imaginary part.
> Was this a good choice or not, I don't know, but this is the current
> algorithm.
A priori, it doesn't seem to be a particularly good choice, no. Not
only is it incompatible with Fortran, it's incompatible with C99, C++
and how C89 libraries typically implement complex, which is similar to
the struct you show below. So if your code makes use of external
libraries, you're going to be spending time converting between the
"normal" representation and yours.
I also suspect that in most cases the "normal" representation will
give better performance due to keeping the real and complex parts of
each number close to each other, thus improving cache efficiency.
> Is this written in the fortran standard ?
The standard says that a complex variable of kind X will take up twice
the storage space as a real variable of kind X, which in practice
means that a Fortran complex will have the real and complex parts
stored after each other. But as other have said, there might be
alignment issues and additionally on the C side there might be struct
packing issues. Best to use the ISO_C_BINDING approach as explained by
Tobias Burnus.
--
JB
relaxmike
(I reply to all message at once.)
> Jugoslav Dujic wrote :
> This is compiler-dependent in theory. In practice, the above
> struct will always map to the Fortran complex*16 for every
> plausible combination of compilers.
> JB wrote :
> Why not use one of the available C/C++ Lapack wrappers rather than
> reinventing the wheel?
The code we have was written long before Lapack was available
and has been maintained "as is" since.
But there is a massive work of the team to move to a new architecture,
based on
C/C++, so this may be the opportunity to use a new connection to
Lapack.
The constraint is that the source code must be available with
a LGPL-like (MIT, BSD, etc...) license. If not, the source code cannot
be
integrated in our software.
After some fast Google searches, I found these projects :
http://www.netlib.org/clapack/
http://math.nist.gov/lapack++/
Are there other C or C++ wrappers for Lapack ?
> JB wrote :
> A priori, it doesn't seem to be a particularly good choice, no. [...]
That's clear now. May be it was not so obvious when the software
was designed, 25 years ago. There may also be other constraints.
The current storage allows a fast extraction of all real or imaginary
parts of an array, which may be efficient in some situations : but
I admit that I do not have any evidence of practical situations where
this
advantage may be of practical value.
> Tobias Burnus wrote :
> Well, on the C side, one can use the C99 datatype:
> double _Complex z;
> and on the Fortran side using a Fortran 2003 feature:
> USE ISO_C_Binding
> complex(c_double_complex) :: z
Using "double _Complex" on the C side makes sense.
But I cannot modify F77's Lapack source code.
> Tobias Burnus wrote :
> I think the f2c compiler and maybe also g77 were using (internally)
> such a struct, i.e. in that case "double _Complex" might not work
> properly.
This explains why our current code compiles (and works !) both with
f2c and
gfortran, based on our handcrafted C struct, on windows and linux.
Regards,
Michaël
relaxmike
C/C++ interfaces to Lapack :
Lapack ++
v1.1, Roldan Pozo, http://math.nist.gov/lapack++/
v1.9, Christian Stimming, http://lapackpp.sourceforge.net/
CLAPACK : http://www.netlib.org/clapack/
LPP : http://sourceforge.net/projects/lpp/
CVMlib : http://www.cvmlib.com/
Regards,
Michaël
George
Though LAPACK is a bunch of fortran routines, I found that the literature
that you find through agoogle search of keywords is most often how to call
them using c, c++.
Is your platform windows? What were you able to accomplish? (If your
answer is "nothing but complete failure," it was the answer I came up with
the first couple timesI tried to use these routines calling from fortran on
windows.
Elliot Chandler and user1 are chipping away at this using g95, so you might
want to focus on their posts and or replies in particular.
Fortran has all this incredible machinery to deal with arrays, but you
can't really do much of anything if you can't calculate a determinant, for
example.
The task that remains after you create the library is to wade through the
documentation, the lion's share of which is the comments for the calling
the routine.
As far as mixing C with fortran with complextypes, there's a lot of
guarantees about width and alignmnet in both C99 and F03, so it would be a
good way to get your feet wet with interop.
--
George
It's going to be the year of the sharp elbow and the quick tongue.
George W. Bush
Picture of the Day http://apod.nasa.gov/apod/