On Sep 30, 2013, at 3:01 AM, Christophe Ortiz <[email protected]>
wrote:
> Hi Barry,
>
> It works ! Great job ! Here is the output:
>
>
> Components have the right values.
>
> Do you think it can be adapted to work with DMDAVecGetArrayF90(da,...) /
> DMDAVecRestoreArrayF90(da,...) in order to work with distributed arrays ?
What dimension do you need? For one dimension, which I think you mentioned
before it is almost identical to what I sent you
PETSC_EXTERN void PETSC_STDCALL dmdavecgetarrayf901_(DM *da,Vec *v,F90Array1d
*a,PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
PetscInt xs,ys,zs,xm,ym,zm,gxs,gys,gzs,gxm,gym,gzm,N,dim,dof;
PetscScalar *aa;
*ierr = DMDAGetCorners(*da,&xs,&ys,&zs,&xm,&ym,&zm);if (*ierr) return;
*ierr = DMDAGetGhostCorners(*da,&gxs,&gys,&gzs,&gxm,&gym,&gzm);if (*ierr)
return;
*ierr = DMDAGetInfo(*da,&dim,0,0,0,0,0,0,&dof,0,0,0,0,0);if (*ierr) return;
/* Handle case where user passes in global vector as opposed to local */
*ierr = VecGetLocalSize(*v,&N);if (*ierr) return;
if (N == xm*ym*zm*dof) {
gxm = xm;
gym = ym;
gzm = zm;
gxs = xs;
gys = ys;
gzs = zs;
} else if (N != gxm*gym*gzm*dof) {
*ierr = PETSC_ERR_ARG_INCOMP;
}
*ierr = VecGetArray(*v,&aa);if (*ierr) return;
*ierr = F90Array1dCreate(aa,PETSC_SCALAR,gxs,gxm,a
PETSC_F90_2PTR_PARAM(ptrd));if (*ierr) return;
}
except it (1) handles both the local vector or global vector version (the if
above) and (2) handles that start index at gxs instead of 0. All you should
have to do is rename the beast above, provide an FORTRAN 90 interface for it,
have it call the specialized version of F90Array1dCreate I already provided.
Two and three dimensions are a bit more involved since you need to copy over
and modify more functions but I don't think there are any tricky things.
Barry
>
> Best regards,
> Christophe
>
>
> 2013/9/29 Barry Smith <[email protected]>
>
>
> Sorry, that is what I get for having a ex1.F and ex1F90.F in the same
> directory.
>
> Barry
>
> On Sep 29, 2013, at 6:16 AM, Christophe Ortiz <[email protected]>
> wrote:
>
> > Hi Barry,
> >
> > I just looked at exf1.F. It seems you sent me a wrong file. It does not
> > correspond to ex1f90. There is no vecgetarraymystruct in the file you sent
> > me. Could you please send me the correct version of ex1f90 so I can test it
> > ? Many thanks in advance.
> >
> > Christophe
> >
> > CIEMAT
> > Laboratorio Nacional de Fusión por Confinamiento Magnético
> > Unidad de Materiales
> > Edificio 2 - Planta 0 - Despacho 28m
> > Avenida Complutense 40,
> > 28040 Madrid, Spain
> > Tel: +34 91496 2582
> > Fax: +34 91346 6442
> >
> > --
> > Q
> > Por favor, piense en el medio ambiente antes de imprimir este mensaje.
> > Please consider the environment before printing this email.
> >
> >
> > 2013/9/27 Barry Smith <[email protected]>
> >
> > On Sep 27, 2013, at 3:06 PM, Christophe Ortiz <[email protected]>
> > wrote:
> >
> > > The files were not attached. Could you please send them again ?
> > >
> > > Thanks.
> > > Christophe
> > >
> > > El 27/09/2013 21:28, "Barry Smith" <[email protected]> escribió:
> > >
> > > Please find attached a new version of the code. ex1.c and ex1f90.F
> > > please run it as is and send me all the output.
> > >
> > > It works correctly with gfortran and the intel 12.x ifort compiler, I
> > > believe it should now work for all fortran compilers. If it does not work
> > > as expected please send me the name and version of your fortran compiler.
> > >
> > > Barry
> > >
> > > On Sep 27, 2013, at 8:22 AM, Christophe Ortiz
> > > <[email protected]> wrote:
> > >
> > > > There was a misunderstanding. There are n entries in total in the
> > > > vector. I thought it was a vector of n points, with 3 components
> > > > associated to each point, ie 3*n entries. This is clear now.
> > > >
> > > > I've tested your example and I checked the values assigned to the
> > > > components with your piece of code:
> > > >
> > > > call VecGetArrayMyStruct(x,xarray,ierr)
> > > > do i=1,10
> > > > xarray(i)%a = i
> > > > xarray(i)%b = 100*i
> > > > xarray(i)%c = 10000*i
> > > > enddo
> > > >
> > > > call VecRestoreArrayMyStruct(x,xarray,ierr)
> > > >
> > > >
> > > > I did it with:
> > > >
> > > > call VecGetArrayMyStruct(x,xarray,ierr)
> > > > do i = 1 , 10
> > > > write(*,*) xarray(i)%a,xarray(i)%b,xarray(i)%c
> > > > end do
> > > >
> > > > call VecRestoreArrayMyStruct(x,xarray,ierr)
> > > >
> > > > I obtained the following:
> > > >
> > > > 1.00000000000000 2.00000000000000 3.00000000000000
> > > > 2.00000000000000 3.00000000000000 4.00000000000000
> > > > 3.00000000000000 4.00000000000000 5.00000000000000
> > > > 4.00000000000000 5.00000000000000 6.00000000000000
> > > > 5.00000000000000 6.00000000000000 7.00000000000000
> > > > 6.00000000000000 7.00000000000000 8.00000000000000
> > > > 7.00000000000000 8.00000000000000 9.00000000000000
> > > > 8.00000000000000 9.00000000000000 10.0000000000000
> > > > 9.00000000000000 10.0000000000000 1000.00000000000
> > > > 10.0000000000000 1000.00000000000 100000.000000000
> > > >
> > > > First column if ok (1st component), but the other two columns show
> > > > wrong numbers. For component 2 and 3, only the very last value is ok.
> > > > It seems values from component b are copied from a, and those from
> > > > component c from b, with an offset.
> > > >
> > > > What do you think ?
> > > >
> > > > Christophe
> > > >
> > > >
> > > > On Fri, Sep 27, 2013 at 2:47 PM, Barry Smith <[email protected]> wrote:
> > > >
> > > > Indices start in fortran arrays from 1 by default so don't even try
> > > > zero
> > > >
> > > > I've attached a new version that has three components,
> > > > 30 TOTAL vector entries and 10 for each component. Please run it and
> > > > send the results back.
> > > >
> > > > Barry
> > > >
> > > > On Sep 27, 2013, at 5:40 AM, Christophe Ortiz
> > > > <[email protected]> wrote:
> > > >
> > > > > Hi Barry,
> > > > >
> > > > > I made some tests with your code that allows to use struct in
> > > > > Fortran90. It compiles and does not complain with the struct with 2
> > > > > components but...it exhibits some strange behaviour. Here are the
> > > > > behaviours I observed. To make it simple, I used n=10.
> > > > >
> > > > > - After using call VecSet(x,one,ierr), I checked the values of each
> > > > > component of the vector using
> > > > >
> > > > > call VecGetArrayMyStruct(x,xarray,ierr)
> > > > > do i = 0 , n
> > > > > write(*,*) i , xarray(i)%a, xarray(i)%b
> > > > > end do
> > > > > call VecRestoreArrayMyStruct(x,xarray,ierr)
> > > > >
> > > > > I checked the values from i= 0 to n since I was not sure what was the
> > > > > convention (0- or 1-based). I obtained:
> > > > >
> > > > > 0 0.000000000000000E+000 1.00000000000000
> > > > > 1 1.00000000000000 1.00000000000000
> > > > > 2 1.00000000000000 1.00000000000000
> > > > > 3 1.00000000000000 1.00000000000000
> > > > > 4 1.00000000000000 1.00000000000000
> > > > > 5 1.00000000000000 1.00000000000000
> > > > > 6 1.00000000000000 1.00000000000000
> > > > > 7 1.00000000000000 1.00000000000000
> > > > > 8 1.00000000000000 1.00000000000000
> > > > > 9 1.00000000000000 1.00000000000000
> > > > > 10 1.00000000000000 1.996650376645798E-314
> > > > >
> > > > >
> > > > > As you can see, for component a, there is a 0 for i=0, then the
> > > > > values are correct. For component b, correct values go from i=1 to
> > > > > 9.There is garbage in second component for i=10 . I checked that if
> > > > > you define a third component c, its values go from i=-1 to n-2. It's
> > > > > like if values of second component started at the end of first
> > > > > component, and same thing for third component.
> > > > >
> > > > > Then, instead of using VecSet(), I tried to assign values directly to
> > > > > each component of the vector with
> > > > >
> > > > > call VecGetArrayMyStruct(x,xarray,ierr)
> > > > > do i = 0 , n-1
> > > > > xarray(i)%a = 2.0
> > > > > xarray(i)%b = 3.0
> > > > > end do
> > > > > call VecRestoreArrayMyStruct(x,xarray,ierr)
> > > > >
> > > > > Here I checked that indices actually start at 0 and not at 1. With a
> > > > > loop from i=1 to n I got a Memory corruption message.
> > > > >
> > > > > and I checked that the values are correctly assigned with
> > > > >
> > > > > call VecGetArrayMyStruct(x,xarray,ierr)
> > > > > do i = 0 , n-1
> > > > > write(*,*) i, xarray(i)%a, xarray(i)%b
> > > > > end do
> > > > > call VecRestoreArrayMyStruct(x,xarray,ierr)
> > > > >
> > > > >
> > > > >
> > > > > Then, I obtained:
> > > > >
> > > > > 0 2.00000000000000 2.00000000000000
> > > > > 1 2.00000000000000 2.00000000000000
> > > > > 2 2.00000000000000 2.00000000000000
> > > > > 3 2.00000000000000 2.00000000000000
> > > > > 4 2.00000000000000 2.00000000000000
> > > > > 5 2.00000000000000 2.00000000000000
> > > > > 6 2.00000000000000 2.00000000000000
> > > > > 7 2.00000000000000 2.00000000000000
> > > > > 8 2.00000000000000 2.00000000000000
> > > > > 9 2.00000000000000 3.00000000000000
> > > > >
> > > > >
> > > > > Here the problem seems more severe. Values are not correctly assigned
> > > > > to the second component. There should be 3.0 in the second column. It
> > > > > seems values of the first component are copied to the second one.
> > > > > Only for i=10 the value of xarray(i)%b is correct (3.0).
> > > > >
> > > > > Any idea where it could come from ?
> > > > > I guess it needs some fixes here and there but I think your idea is
> > > > > good and that it could work.
> > > > >
> > > > > Christophe
> > > > >
> > > > >
> > > > >
> > > > > On Thu, Sep 26, 2013 at 5:45 PM, Barry Smith <[email protected]>
> > > > > wrote:
> > > > >
> > > > > I just put it with the Fortran source code and compile it with the
> > > > > rest of the application code; here is the makefile I used
> > > > >
> > > > >
> > > > >
> > > > >
> > > > > On Sep 26, 2013, at 10:14 AM, Christophe Ortiz
> > > > > <[email protected]> wrote:
> > > > >
> > > > > > By the way, what should I do with the small .c code ? Where should
> > > > > > I put it ?
> > > > > >
> > > > > > Christophe
> > > > > >
> > > > > >
> > > > > >
> > > > > > On Thu, Sep 26, 2013 at 4:09 PM, Barry Smith <[email protected]>
> > > > > > wrote:
> > > > > >
> > > > > > Christophe,
> > > > > >
> > > > > > Despite my LinkIn Endorsement for expertise in Fortran :-) I
> > > > > > cannot pretend to be an expert in FortranXX but I have cooked up an
> > > > > > example demonstrating accessing the Vec entries as if they are in
> > > > > > an array of derived types. I've attached the example code; there
> > > > > > needs to be a small C stub that defines the functions for your
> > > > > > specific derived type name.
> > > > > > Note that it will only work I think if your N is a compile time
> > > > > > constant.
> > > > > >
> > > > > >
> > > > > >
> > > > > >
> > > > > > It worked with
> > > > > > ~/Downloads$ gfortran --version
> > > > > > GNU Fortran (GCC) 4.8.1 20130404 (prerelease)
> > > > > >
> > > > > >
> > > > > > I do not understand exactly why it works since it uses
> > > > > > F90Array1dCreate(fa,PETSC_SCALAR,1,len,ptr
> > > > > > PETSC_F90_2PTR_PARAM(ptrd)); which has a single PETSC_SCALAR as a
> > > > > > building block but … I hope it works for you. If it doesn't, let us
> > > > > > know the compiler you are using and we may be able to get it
> > > > > > working for that compiler.
> > > > > >
> > > > > > Barry
> > > > > >
> > > > > >
> > > > > >
> > > > > >
> > > > > > On Sep 26, 2013, at 4:41 AM, Christophe Ortiz
> > > > > > <[email protected]> wrote:
> > > > > >
> > > > > > > Hi all,
> > > > > > >
> > > > > > > Me again !
> > > > > > >
> > > > > > > I have read in previous posts that it is hard in Fortran to
> > > > > > > declare something similar to a typedef struct in C to manage a
> > > > > > > multicomponent problem.
> > > > > > >
> > > > > > > Is it still the case ? Has the problem been solved ?
> > > > > > >
> > > > > > > I am asking because my plan is to implement a multicomponent
> > > > > > > problem (in 1D), with many components that will be organized in
> > > > > > > arrays of two dimensions. In C I could define
> > > > > > >
> > > > > > > typedef struct{
> > > > > > > PetscScalar U;
> > > > > > > PetscScalar V;
> > > > > > > PetscScalar A[N][N];
> > > > > > > } Field;
> > > > > > >
> > > > > > > and then I could calculate the residual function with
> > > > > > >
> > > > > > > F[i].U = ...
> > > > > > > F[i].V = ...
> > > > > > > F[i].A[k][n] = ...
> > > > > > >
> > > > > > > which is quite convenient.
> > > > > > >
> > > > > > > If in Fortran it is not possible to use a struct as in C, I am
> > > > > > > afraid I'll have to deal with
> > > > > > >
> > > > > > > F(jdof,i) = ...
> > > > > > >
> > > > > > > where I will have only jdof to address U, V and A[ ][ ], which
> > > > > > > can be difficult and not very convenient I guess. Before I move
> > > > > > > all my code to C, does anyone have an alternative idea to manage
> > > > > > > this multi(many)component problem in Fortran ?
> > > > > > >
> > > > > > > Many thanks in advance for your help and suggestion !
> > > > > > >
> > > > > > > Christophe
> > > > > >
> > > > > >
> > > > > >
> > > > > >
> > > > >
> > > > >
> > > > >
> > > >
> > > >
> > > >
> > >
> >
> >
>
>
