Attachment: ex1f.F
Description: Binary data

Attachment: ex1.c
Description: Binary data

Attachment: makefile
Description: Binary data

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
> > > >
> > > >
> > > >
> > > >
> > >
> > >
> > >
> >
> >
> >
> 

Reply via email to