Friday, October 25, 2013

Passing Assumed Shape Arrays between Fortran and C++


Fortran 90 introduced assumed shape arrays (colloquially referred to as "Fortran 90 Arrays"), a big step forward.  Fortran 2003 introduced interoperability with C/C++, another step forward.  Unfortunately, there is no well-known way to pass assumed shape arrays between Fortran and C++.  This is a serious problem: assumed shape arrays are widely used in modern Fortran code, and not being able to interoperate easily with C/C++ in those cases can make the C interop features nearly useless.

Here, we present a standards-compliant way to allocate assumed shape arrays in Fortran and pass them to code written in C++.  This method is fully compliant with the Fortran 2003 and C++98/11 language standards.


We begin by using a C++ library that provides functionality equivalent to Fortran's assumed shape arrays.  We used Blitz++; however, others have suggested Armadillo and Eigen.  Although we have not investigated, we believe the techniques outlined here would work for these other packages as well.

These techniques will also require a Fortran 2003 compiler that supports the ISO standard Fortran/C interop.


With appropriate information on memory layout, Blitz++ is able to create array structures that reference already-allocated areas of memory.  The required information is stored in the Fortran "dope vector."  Unfortunately, dope vectors are implement-dependent, and cannot be read directly.

We proceed by reconstructing the Fortran dope vector using a series of standard Fortran 2003 calls.  Consider the following structure, which can store information on a 2-dimensional dope vector:
type, bind(c) :: arr_spec_2    type(c_ptr) :: base    type(c_ptr) :: deltas(2    integer(c_int) :: lbounds(2)    integer(c_int) :: ubounds(2)end type
We then fill it in.  Here is how that can be done for 2-D assumed shape arrays of type real*8:
    function c_loc_double(x)
        use, intrinsic :: iso_c_binding
        real*8, target :: x
        type(c_ptr) :: c_loc_double
        c_loc_double = c_loc(x)
    end function
    ! ================ Type (real*8, double), Rank 2
    subroutine get_spec_double_2(arr, low1,low2, spec)
    implicit none
    real*8, dimension(:,:), target :: arr
    integer :: low1,low2
    type(arr_spec_2) :: spec
        spec%base = c_loc_double( arr(lbound(arr,1),lbound(arr,2)) )
        ! ------- Dimension 1
        spec%lbounds(1) = low1
        spec%ubounds(1) = low1 + ubound(arr,1) - 1
        if (spec%lbounds(1) < spec%ubounds(1)) then
            spec%deltas(1) = c_loc_double(arr(lbound(arr,1)+1,lbound(arr,2)))
            spec%deltas(1) = spec%base
        end if
        ! ------- Dimension 2
        spec%lbounds(2) = low2
        spec%ubounds(2) = low2 + ubound(arr,2) - 1
        if (spec%lbounds(2) < spec%ubounds(2)) then
            spec%deltas(2) = c_loc_double(arr(lbound(arr,1),lbound(arr,2)+1))
            spec%deltas(2) = spec%base
        end if
    end subroutine
A few things to note:

  1. The function c_loc_double(x) is required because c_loc() does not seem to work directly on elements of arrays (at least not in GNU's gfortran).
  2. The deltas() array contains the address of the second element in each dimension.  True skip offsets for each dimension can be computed by spec%deltas(i) - spec%base.  This is not done because Fortran does not have pointer arithmetic.

Fortran 2003 does not preserve the lower bound of an array when you pass it.  Therefore, the lower bounds of each dimension have to be passed into the subroutine get_spec_double_2().  In order to make life more convenient, this can be hidden under the hood with the following macro (if one is using Fortran with the C pre-processor):
#define GET_SPEC_DOUBLE_2(arr, spec) \    get_spec_double_2(arr, lbound(arr,1),lbound(arr,2), spec)

With the dope vector stored in a type(arr_spec_2) object, that information may now be passed to a C++ subroutine via a matching C++ structure --- and then reconstituted into a Blitz++ array.  We can use C++ templates to make a version that will work for all types and ranks:
/** Used to accept Fortran 90 arrays and re-constitute them as blitz++ arrays.
@see */
template<class ArrT, int rank>
struct F90Array {
    ArrT *base;
    ArrT *deltas[rank];
    int lbounds[rank];
    int ubounds[rank];
    /** Extract F90Array info from existing blitz::Array.  Used to
    write C++ test code for Fortrn APIs. */
    F90Array(blitz::Array<ArrT,rank> &arr) {
        this->base =;
        blitz::TinyVector<int, rank> idx(0);
        for (int i=0; i<rank; ++i) {
            this->deltas[i] = this->base + arr.stride(i);
            this->lbounds[i] = arr.lbound(i);
            this->ubounds[i] = arr.ubound(i);

    blitz::Array<ArrT,rank> to_blitz()
        blitz::TinyVector<int, rank> shape, stride;
        blitz::GeneralArrayStorage<rank> stor;
        for (int i=0; i<rank; ++i) {
            shape[i] = ubounds[i] - lbounds[i] + 1;
            stride[i] = deltas[i] - base;
            stor.base()[i] = lbounds[i];
            // Ordering is not needed because we're using stride
            // stor.ordering()[i] = i;      // Fortran ordering, blitz++ manual p31
        return blitz::Array<ArrT,rank>(base, shape, stride,
            blitz::neverDeleteData, stor);


The code may now be used as follows:

interface    subroutine c_function(arr_spec) bind(c)        type(arr_spec_2) :: arr_spec    end subroutine c_functionend interface
real*8, dimension(:,:) :: arrtype(arr_spec_2) :: arr_spec
    call GET_SPEC_DOUBLE_2(arr, arr_spec)    call c_function(arr_spec)
-------------------------------extern "C" void c_function(F90Array<double, 2> &arr_spec){    auto arr(arr_spec.to_blitz());    printf("Here is a value: %f\n", arr(17,34));}


The above code works only for rank-2 arrays of real*8.  The same principles can be applied to extract and use dope vectors for any size and kind of array; but because Fortran lacks templates, writing out the required code would be tedious.  For this reason, we provide a Python script, to generate the required Fortran code.  This script really needs only be run once.  Its usage is:
python <output.f90> <output.h>
The variables fctypes and ranks must be modified inside the script, to control which data types and ranks of arrays get generated.  This can be customized as needed for any particular program.

The required code is only two files: a templated C++ header file, and a Python script to generate the required Fortran files.  That code may be downloaded here:

Tuesday, August 27, 2013

Can Fracking Save Us from Climate Catastrophe?

It has been noted that US carbon emissions have gone down more than that of any other industrialized nation in the past few years.  This has been used as evidence to support the American laissez-faire approach to the entire issue.

The decrease is due mainly to fracking and the shale gas boom.  There is so much cheap natural gas out there that utilities are switching power generation from coal to gas.  It's a win for utilities, win for consumers, and win for the environment!  If only it were that simple.  For one thing, any change in the relative price of coal vs. natural gas can reverse the "progress" on emissions:

But there is a more serious problem.  Let's suppose for now that fracking releases so much gas that the world stops burning coal, switching to gas instead.  Heck, let's suppose that we stop burning oil as well, and convert all our fossil fuel use to gas!  And let's imagine that this conversion will reduce the world's carbon footprint by 50%.  Where does that get us?

Well... if the world was hurtling toward a cliff at 100mph, now we're speeding toward that very same cliff at 50mph.  It means that all the climate-related problems we will be facing will take twice as long to materialize.  That would certainly be good news.

But it would not really solve the problem.  Humanity is currently releasing huge amounts of CO2 into the atmosphere.  When similar amounts of CO2 were released "suddenly" in the past, bad things happened... 10 degrees (F) of warming, mass extinctions, long-term changes in the course of the history of life.  Look, for example, at the end-Permian and end-Triassic extinctions.  The extinction issues will likely persist long after all the CO2 we've released into the atmosphere gets mopped up into the rocks due to the long-term carbon cycle.  

The problem here is that both of the above extinction events were the result of "sudden" releases of CO2.  How sudden?  Paleontologists think about a million years.  Which is 10,000 times slower than we're releasing CO2 right now.  It doesn't matter if we release that much CO2 over 100 years or 200 years, the end result will be the same.  All that would have to change is Jim Hansen would have to change the title of his book from "Storms of my Grandchildren" to "Storms of my Great-Great-Grandchildren."

So let's quit the charade of how natural gas and fracking will save us.  It will only ever get us halfway there.  Instead, let's focus on REAL solutions to the climate problem, which require zero-carbon technologies.  The best that natural gas will ever buy us is a little bit of time.

[By the same token, one could consider one of my favorite technologies, the electric bicycle.  Suppose we stop driving cars and drive electric bikes instead.  These use dramatically less fuel than automobiles, about 10% that of a Nissan Leaf.  Heck, let's suppose we cut ALL our fossil fuel use by 90%.  Would we still need to worry about climate change?  Unfortunately, yes.  Now the catastrophe that would have taken 100 years will take 1000 years to unfold.  Maybe that's OK for some people, although the mass extinctions would still happen similarly in either case.]

Saturday, July 27, 2013

Is Solar Net Metering REALLY a Threat to Utilities?

Intereting article in the NY Times Blog:

For years, power companies have watched warily as solar panels have sprouted across the nation’s rooftops. Now, in almost panicked tones, they are fighting hard to slow the spread...        

The real problem seems to be that solar customers pay nothing under net metering, even though they use the grid.  A few thoughts:

1. Develop a rate for reverse electric flow, as well as forward.  Then solar customers will get charged for their use of the grid, whether it was forward or reverse.  Meanwhile, solar customers should be paid the prevailing wholesale rate for the electricity they provide AT THE TIME THEY PROVIDE IT.  That is, they should get the high daytime wholesale rate if they provide electricity during a hot summer day with all the A/C going.  This seems to be the most economically straightforward way to price things.

2. Stop subsidizing solar, and put in place a sensible carbon tax (aka Jim Hansen).  With such a tax in place, suddenly solar won't need any subsidies.

3. If utilities don't like distributed generation, regulators can suggest a willingness to change policy --- but ONLY IF the utilities meet certain benchmark targets for renewable content in the electricity they provide.