Passing two-dimensional dynamic array from C to Fortran

[Intel_C_Intel]

Hello

I have some trouble. I want to pass two-dimensional dynamic array from VC++ to IVCs DLL. My DLL subroutine:

Code:

Subroutine DecompS( N, A, cond, IPvt, MatrixType )
	!DEC$ ATTRIBUTES DLLEXPORT, C, DECORATE, ALIAS : 'DecompS' :: DecompS
	!DEC$ ATTRIBUTES REFERENCE :: cond

	integer, parameter :: R_ = Single_
	
	integer(I_), intent(in) :: N, MatrixType
	real(R_), intent(inout) :: A(N, N)
	integer(I_), intent(out) :: IPvt(N)
	real(R_), intent(out) :: cond

	print *, A
end subroutine DecompS

In VC++ I declare:

Code:

extern "C" __declspec(dllimport) void DecompS( int N, float **A, float *cond,
	int *IPvt, int MatrixType );

int _tmain(int argc, _TCHAR* argv[])
{
	int N;
	int i, j;
	int *IndexPivot;
	float **A;
	float cond;

	N = 3;
	IndexPivot = new int;
	A = new float*;
	for( i = 0; i < N; i ++ ) {
		A = new float;
	}
	A[0][0] = 10;
	A[0][1] = -7;
	A[0][2] = 0;
	A[1][0] = -3;
	A[1][1] = 2;
	A[1][2] = 6;
	A[2][0] = 5;
	A[2][1] = -1;
	A[2][2] = 5;
	for( i = 0; i < N; i ++ ) {
		for( j = 0; j < N; j ++ )
            printf( "%5.0f", A );
		printf( "
" );
	}

	DecompS( N, A, &cond, IndexPivot, 1 );
	
}

DecompS prints another numbers. If I write in VC++:

Code:

extern "C" __declspec(dllimport) void DecompS( int N, float *A, float *cond,
	int *IPvt, int MatrixType );
DecompS( N, *A, &cond, IndexPivot, 1 );

DecompS prints first string Ok, thats because C++ is row-ordered, but next numbers from nowhere, they are absolutely another. In memory they must be together. Whats going on there?

[Jugoslav_Dujic]

Well, you described it yourself in a previous thread -- C (generally) doesn't have dynamic 2D arrays, but only arrays of pointers to 1D arrays. And that's not the same as Fortran 2D array.

So, you either have to

- have a 1D array on C side, which you have to manipulate with address mathematics rather than [][], or

- write a more appropriate Fortran structure, which is a pain. That should be like:

  type t_row
    integer(int_ptr_kind(prow)):: pRow
  end type t_row
  ...
  Subroutine DecompS( N, A, ...)
  type(t_row):: A(*)
  
  real(R_)::    row(N); pointer(pRow, row) !A "cray" pointer
  ...
  do i=1, N
    pRow = A(i)%pRow
    do j=1,N
      row(N) = whatever
    end do
  end do

(Replace "row" with "column" depending on which is easier to read for you).

Neither solution is particularly appealing though. To me, it looks best to define a class/struct in C++ holding a 2-D array in 1-D form with dimension information, then define operator [] for it. In that way, you will get both C++ and Fortran code "looking like" you work with a 2-D array all the time.

Jugoslav

[Intel_C_Intel]

Thank you, Jugoslav. I preferred first method. Second requires serious changes in Fortran program thats real pain, for it. But the last your approach is the most professional. But it seems to have some strange indexation: Matrix[Index(i, j)]. Index is a structure {int i, int j}. We know that impossible overload double index like in sharp. So this is what I get:

Code:

extern "C" __declspec(dllimport) void DecompS( int N, float *A, float *cond,
	int IPvt[], int MatrixType );

int _tmain(int argc, _TCHAR* argv[]) {
	int N;
	int i, j;
	int *IndexPivot;
	float *A;
	float cond;

	N = 3;
	IndexPivot = new int;

	A = new float[N * N];
	A[0*N + 0] = 10;
	A[0*N + 1] = -7;
	A[0*N + 2] = 0;
	A[1*N + 0] = -3;
	A[1*N + 1] = 2;
	A[1*N + 2] = 6;
	A[2*N + 0] = 5;
	A[2*N + 1] = -1;
	A[2*N + 2] = 5;
	for( i = 0; i < N; i ++ ) {
		for( j = 0; j < N; j ++ )
            printf( "%5.0f", A[i*N + j] );
		printf( "
" );
	}

	DecompS( N, A, &cond, IndexPivot, 1 );
	return 0;
}

It works but looks a little bit baldy for C++. Next probably stepshould beoperator [] and class Matrix: A[Index(i, j)]. But I wouldnt do this too say and probably slower. "That's C++", what can I say more? Anyway, I get what I want. Thanks again.

[Jugoslav_Dujic]

Shouldn't be too difficult. My C++ is rusty, but since I'm curious too, here's a quick sample:

#include "stdio.h"

struct Matrix
{
  float* pf;
  int nRows;
  int nCols;
  float* operator [] (int iRow)
      {return (pf + iRow*nCols);};
  Matrix(int iRows, int iCols)
      {
         nRows = iRows;
         nCols = iCols;
         pf = new float[nRows*nCols];
      }
  ~Matrix()
    {delete[] pf;};
};

int main(void)
{
   Matrix A(6,5);
   for (int i=0;i<6;i++)
      for (int j=0;j<5;j++)
        A = float(100*i + j+1);
   
   for (i=0;i<6;i++)
   {
      for (int j=0;j<5;j++) printf("(%10.1f  )", A);
      printf ("
");
   }

   return 1;
}

Jugoslav

Message Edited by JugoslavDujic on 08-24-2005 09:46 AM

[Jugoslav_Dujic]

P.S. The code above works (and does not produce any memory problems), but please carefully check what I wrote regarding rows vs. columns -- it's quite possible that I mixed them up so the poor reader gets confused. I can correctly process that kind of stuff only after very careful thinking, but I don't feel like doing it right now :-).

Jugoslav

[Intel_C_Intel]

Thanks fro approach, Jugoslav.

Cool. I would never guessed. Actually thats because I do not understand why does double operator A work. Theoretically its a (pf + iRows*nCols)+jRows*nCols but works like a (pf + iRows*nCols)+jRows. It seems that second operator means just + for address. This useful and means that with tree-dimensional array problems may happen.

But it isnt important for me now. Ive some trouble with this code:

Code:

#include "stdafx.h"

struct Matrix {
	float* pf;
	int nRows;
	int nCols;
	
	float* operator []( int iRow ) {
		return( pf + iRow*nCols );
	}
	
	Matrix() {
	}
	
	Matrix( int iRows, int iCols ) {
		nRows = iRows;
		nCols = iCols;
		pf = new float[nRows * nCols];
	}

	~Matrix() {
		delete [] pf;
	}
};

int _tmain(int argc, _TCHAR* argv[])
{
	Matrix A;
	// Matrix A = Matrix(1, 2);
	A = Matrix(6, 5);
	
	for( int i = 0; i < 6; i ++ )
		for( int j = 0; j < 5; j ++ )
			A = float( ( i + 1)*100 + j + 1 );

	for( i = 0; i < 6; i ++ ) {
		for( int j = 0; j < 5; j ++ )
			printf( "(%10.1f  )", A );
		printf( "
" );
	}

	getchar();
	return 0;
}

I added default constructor and default Matrix A declaration. Program throws exception. Somethings going on with memory. The same will be if uncom-ment string and comment Matrix A. Do you have any ideas?

[Jugoslav_Dujic]

As I said, I'm not a C++ expert :-). I played a bit with the code and concluded that you mustn't do:

A = Matrix(6, 5);

When I step into the constructor in the debugger, it says:

this	CXX0069: Error: variable needs stack frame

i.e. the constructor is called for a non-existing object. Instead, you have to give it a frame somehow:

int main(void)
{
   Matrix *A;
   A = new Matrix(6,5);
   for (int i=0;i<6;i++)
      for (int j=0;j<5;j++)
        (*A) = float(100*i + j+1);
   
   for (i=0;i<6;i++)
   {
      for (int j=0;j<5;j++) printf("(%10.1f  )", (*A));
      printf ("
");
   }

   return 0;
}

In this way, the object is constructed on the heap.

Jugoslav

[Intel_C_Intel]

Right yours. I forgot that struct is not in a heap. Now its looking and working good. Thanks for help, Jugoslav.