On Oct 1, 2013, at 2:18 AM, Christophe Ortiz <[email protected]> wrote:
> Hi Barry,
>
> Thanks for the piece of code. I will try to implement it. My problem is 1D,
> so should be ok. If I have some pb I'll let you know.
>
> If this work, do you plan to include it in next version of PETSc ?
I put the example with VecGetArray() into petsc-dev and if you send me one
with DMDAVecGetArray() I'll add that as well. Since the derived type has to
appear explicitly in the code I can't put something in the library only in
examples that people can copy and change.
Barry
> Christophe
>
> On Mon, Sep 30, 2013 at 2:04 PM, Barry Smith <[email protected]> wrote:
>
> 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
> > > > > > >
> > > > > > >
> > > > > > >
> > > > > > >
> > > > > >
> > > > > >
> > > > > >
> > > > >
> > > > >
> > > > >
> > > >
> > >
> > >
> >
> >
>
>
