Prev: class
Next: Intel fortran (em64t) very low-accuracy results on x86_64
From: Luka Djigas on 26 Dec 2009 14:29
Upon translating some matlab code, I came across this :

M1x(end+1,:) = (-h/c) * (A3*V2+A4*V1); %N

(don't mind the right side, it was just copy pasted)

The left side is what interests me. As I understand it lengthens the
array by one element, and gives it the value of the result of right
side.

Now I've never done anything similar in fortran. What would be (can it
be done?) the easiest way to do it (since I have half a page of those
to translate) ?


with regards,
and wishes for a better New Year !

Luka
From: Richard Maine on 26 Dec 2009 15:17
Luka Djigas <ldigas@___gmail___.com> wrote:

> The left side is what interests me. As I understand it lengthens the
> array by one element, and gives it the value of the result of right
> side.
>
> Now I've never done anything similar in fortran. What would be (can it
> be done?) the easiest way to do it (since I have half a page of those
> to translate) ?

In f2003, that could be written, with array allocatable, as

array = [array, new_element]

I'm basing this just on your description; I don't know enough Matlab to
tell you how accurately this mimics the actual Matlab.

The RHS is an array constuctor, with the elements of array, followed by
new_element. That part you can do in f90 (though you have to use the
ugly (/... /) instead of the [...] form). F2003 allows the
auto-reallocation of array with the assignment.

Without f2003, you pretty have to have the allocation as a separate
step, which also means that you have to have 2 arrays because the
allocation will destroy the original. It is most easily done with
pointers instead of allocatables. Doing it with allocatables ends up
requiring twice as many data copies and allocations, in addition to more
lines of code. F2003 adds the move_alloc intrinsic, which facilitates
doing it wil allocatables, but if you have f2003, the auto-reallocation
form seems simpler.

A pure f90 form with pointers looks something like

allocate(tmp_array(size(array+1))
tmp_array = (/ array, new_element /)
deallocate(array)
array -> tmp_array

Be aware that if you are building a large array by concatenating a lot
of elements like this, it is likely to be *FAR* more efficient to
reallocate in bigger chunks than each element.

--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain
From: dpb on 26 Dec 2009 15:41
Richard Maine wrote:
> Luka Djigas <ldigas@___gmail___.com> wrote:
>
>> The left side is what interests me. As I understand it lengthens the
>> array by one element, and gives it the value of the result of right
>> side.
>>
>> Now I've never done anything similar in fortran. What would be (can it
>> be done?) the easiest way to do it (since I have half a page of those
>> to translate) ?
>
> In f2003, that could be written, with array allocatable, as
>
> array = [array, new_element]
>
> I'm basing this just on your description; I don't know enough Matlab to
> tell you how accurately this mimics the actual Matlab.
....

In Matlab, the LHS ":" array subscript is the entire length of the
corresponding dimension.

So, the append of the loop is adding/appending a row vector not simply a
single element in all likelihood. That would depend on the length of
the RHS; it would have to be conformant to the size of the array in that
dimension on the LHS to avoid a size mismatch runtime error.

--
From: dpb on 26 Dec 2009 15:42
Richard Maine wrote:
....

> Be aware that if you are building a large array by concatenating a lot
> of elements like this, it is likely to be *FAR* more efficient to
> reallocate in bigger chunks than each element.

The same is true in Matlab of course; one can speed up ML quite a lot in
many instances if it is possible to "preallocate".

--


From: James Van Buskirk on 26 Dec 2009 16:52
"dpb" <none(a)non.net> wrote in message
news:hh5seb$du8$1(a)news.eternal-september.org...

> In Matlab, the LHS ":" array subscript is the entire length of the
> corresponding dimension.

> So, the append of the loop is adding/appending a row vector not simply a
> single element in all likelihood. That would depend on the length of the
> RHS; it would have to be conformant to the size of the array in that
> dimension on the LHS to avoid a size mismatch runtime error.

I have an example which should work whether or not the RHS is a
vector of the right length of a scalar:

C:\gfortran\clf\concat_array>type concat_array.f90
program concat_array
implicit none
real, allocatable :: M1x(:,:)
integer N
real h,c
real A3,V2,A4,V1
integer i
character(20) fmt
real, allocatable :: temp(:,:)

h = 6.626
c = 2.998
A3 = 3
V2 = 0.7
A4 = -4
V1 = 0.5
N = 4
write(*,'(a)') 'First attempt didn''t work because gfortran'// &
' doesn''t have allocate on assignment'
allocate(M1x(3,N))
M1x = reshape([(i,i=1,size(M1x))],shape(M1x))
write(fmt,'(a,i0,a)') '(',size(M1x,2),'f5.1)'
write(*,'(a)') 'Before concatenation M1x ='
write(*,fmt) transpose(M1x)
! M1x = reshape(transpose(M1x),shape(M1x)+[1,0], &
! pad=[(-h/c)*(A3*V2+A4*V1)],order=[2,1])
write(fmt,'(a,i0,a)') '(',size(M1x,2),'f5.1)'
write(*,'(a)') 'After concatenation M1x ='
! write(*,fmt) transpose(M1x)
write(*,fmt) transpose(reshape(transpose(M1x),shape(M1x)+[1,0], &
pad=[(-h/c)*(A3*V2+A4*V1)],order=[2,1]))
write(*,'(a)') 'Second attempt uses MOVE_ALLOC'
deallocate(M1x)
allocate(M1x(3,N))
M1x = reshape([(i,i=1,size(M1x))],shape(M1x))
write(fmt,'(a,i0,a)') '(',size(M1x,2),'f5.1)'
write(*,'(a)') 'Before concatenation M1x ='
write(*,fmt) transpose(M1x)
call MOVE_ALLOC(M1x,temp)
allocate(M1x(size(temp,1)+1,size(temp,2)))
M1x(:size(temp,1),:) = temp
M1x(size(M1x,1),:) = (-h/c)*(A3*V2+A4*V1)
write(fmt,'(a,i0,a)') '(',size(M1x,2),'f5.1)'
write(*,'(a)') 'After concatenation M1x ='
write(*,fmt) transpose(M1x)
end program concat_array

C:\gfortran\clf\concat_array>gfortran concat_array.f90 -oconcat_array

C:\gfortran\clf\concat_array>concat_array
First attempt didn't work because gfortran doesn't have allocate on
assignment
Before concatenation M1x =
1.0 4.0 7.0 10.0
2.0 5.0 8.0 11.0
3.0 6.0 9.0 12.0
After concatenation M1x =
1.0 4.0 7.0 10.0
2.0 5.0 8.0 11.0
3.0 6.0 9.0 12.0
-0.2 -0.2 -0.2 -0.2
Second attempt uses MOVE_ALLOC
Before concatenation M1x =
1.0 4.0 7.0 10.0
2.0 5.0 8.0 11.0
3.0 6.0 9.0 12.0
After concatenation M1x =
1.0 4.0 7.0 10.0
2.0 5.0 8.0 11.0
3.0 6.0 9.0 12.0
-0.2 -0.2 -0.2 -0.2

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end


 |  Next  |  Last
Pages: 1 2
Prev: class
Next: Intel fortran (em64t) very low-accuracy results on x86_64