
forums.silverfrost.com Welcome to the Silverfrost forums

View previous topic :: View next topic 
Author 
Message 
AntonioGJerez
Joined: 01 Jun 2015 Posts: 9 Location: Spain

Posted: Mon Jun 01, 2015 1:51 pm Post subject: error in array construction 


Dear Friends
This code produces a runtime Access Violation Error when trying to construct an array which includes elements from an arrayvalued function
Moreover, If I define REAL kk(5), then kk=(/ 0.0, diag4(matrix) /) gives a compilation error "Nonconformant array shapes infirst Rank of an array expression (5 and 2)"
It seems that diag4(matrix) is interpreted as a real instead of as a 4 elements array.
Any sugestion?
Thank you in advance
Antonio
program bug
interface diag4
! Extracts the diagonal from a 4x4 matrix (assumed to be square)
function diag4(A)
real,dimension(4,4),intent(in)::A
real,dimension(4)::diag4
end function diag4
end interface
real matrix(4,4)
matrix=reshape((/1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,16./),(/4,4/))
print*,diag4(matrix) ! this runs ok
print*,(/ 0.0, diag4(matrix) /) ! this leads to Access Violation Error
end program
function diag4(A)
implicit none
real,dimension(4,4),intent(in)::A
real,dimension(4)::diag4
integer index
do index=1,size(diag4)
diag4(index)=A(index,index)
enddo
return
end function


Back to top 


JohnCampbell
Joined: 16 Feb 2006 Posts: 2239 Location: Sydney

Posted: Tue Jun 02, 2015 8:34 am Post subject: 


I find the following works:
real out(5)
...
out = 0
out(2:5) = diag4(matrix)
print*, out
I am not sure if out = (/ 0., diag4(matrix) /) should work, but I suspect not, although
out = (/ 0., (/ 1., 6., 11., 16./) /) or
out = (/ (/ 0., 1./), (/ 6., 11., 16./) /)
does appear to work.
John 

Back to top 


AntonioGJerez
Joined: 01 Jun 2015 Posts: 9 Location: Spain

Posted: Tue Jun 02, 2015 8:00 pm Post subject: 


Hi John
thanks, nevertheless, I would like to include the array (/ 0., diag4(matrix) /) into a more complicated expression,
e.g. ... * (/ 0., diag4(matrix) /)+...
it should be better to avoid storing (/ 0., diag4(matrix) /) in a variable in several steps.
I have checked that the above array constructor does work a scalar function defined as
function diag1(A)
real diag1
...
nevertheless it does not work indeed with a one element array:
function diag1(A)
real,dimension(1)::diag1
...
regards 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 1429

Posted: Wed Jun 03, 2015 11:31 pm Post subject: 


Your post touches multiple issues, and I think that we can tackle them better if we take them one at a time.
1. Matlab and Fortran are different. Matlab is an interpreted language and specialized for one thing: manipulating matrices. It has rather weak notions of types and is not suited for general computing. Fortran is a compiled language and newer versions of Fortran (e.g. F95) have stricter rules regarding types and mixing types. A discussion of the Matlab way visavis the Fortran way is probably not appropriate in the Silverfrost forums.
2. When the compiler rejects a piece of code, we may have to consult the relevant Fortran standard to interpret the error messages given by the compiler. If the reading of the standard indicates that the code is conforming, file a compiler bug report with the shortest example code that demonstrates the bug, and cite the relevant sections of the standard. Your code demonstrates such a bug, and you may proceed to file a bug report. When you do, refrain from describing alternative codes that work, but do state if you have tried other compilers and they compile the code correctly even though FTN95 does not.
3. If there is a certain idiom that best expresses a calculation, but the language standard in effect (F95 for FTN95) does not allow it, describe the idiom and ask for alternative ways of expressing the idiom in standard Fortran.
The following program falls into the third category that I just described. It uses the "implied DO loop" construct to assemble a 1D array using one or more constants and the diagonal elements of a 2D array. Note that it avoids using a function with a return value that is an array. Such functions are problematic because they require an explicit interface wherever they are used and, if no such interface is provided, cause erratic run time behaviour.
Code:  program nobug
implicit none
integer, parameter :: M = 4
real, parameter :: Pi = 3.1415927
real matrix(M,M)
real :: diag(M+1)
integer :: i
matrix=reshape((/ 1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,16. /), &
(/ M,M /))
diag=(/ Pi, (matrix(i,i), i=1,M) /)
print *,diag
print *,(/ Pi, (matrix(i,i), i=1,M) /)
end program



Back to top 


AntonioGJerez
Joined: 01 Jun 2015 Posts: 9 Location: Spain

Posted: Thu Jun 04, 2015 12:48 am Post subject: Thanks. Nevertheless... 


Ok, thanks for your code (even though the proposal of supporting an arbitrary function giving a vector would be much more powerful).
Nevertheless, you say "Such functions are problematic because they require an explicit interface wherever they are used and, if no such interface is provided, cause erratic run time behaviour"
Note that in the original example of this post, the interface is provided, the function works but the array constructor desnīt work
Antonio 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 1429

Posted: Thu Jun 04, 2015 4:08 am Post subject: 


It is possible to use an implied DO loop to write a function that gives the diagonal of an arbitrary matrix. Here is the code. Note that I had to provide an alternative syntax in Line6 to work around the FTN95 bug with array constructors.
This code is open to criticism for allocating an array for the function result every time that the function is used. The rules of Fortran 2003 may help to overcome this type of "memory leak" with the "allocate on assignment" feature.
Code:  program getdiag
implicit none
real matrix(2,2),d3(3)
!
matrix=reshape((/ 1.,2.,3.,4. /),(/2,2/))
d3(1) = 0; d3(2:3)=diagf(matrix)
! d3 = (/ 0.0, diagf(matrix) /) !should work, but FTN95 bug is in the way
print*,d3
contains
function diagf(A)
real, dimension(:,:), intent(in) :: A
real, dimension(:),pointer :: diagf
integer :: i,m
m=min(size(A,1),size(A,2))
!
allocate(diagf(m))
diagf = (/ (A(i,i), i=1,m) /)
return
end function
end program



Back to top 




You cannot post new topics in this forum You cannot reply to topics in this forum You cannot edit your posts in this forum You cannot delete your posts in this forum You cannot vote in polls in this forum

Powered by phpBB © 2001, 2005 phpBB Group
