Newbie Question...
Andrea
I'm not an expert with fortran, but I would like to ask 2
questions:
1) Is there a fortran function/subroutine that returns the indices of
a vector (matrix) X that are non-zero? (for whom that know matlab,
like the "find" function).
2) If this is not the case, could you please tell me what is wrong in
this routine? When debugging, it still keeps telling me that "out" is
an "undefined pointer/array" while it's not...
!!!!!!!!!!!!!!!!!!!! CODE BEGIN !!!!!!!!!!!!!!!!!!!!!
program finder
integer m,n
integer, pointer, dimension(:,:) :: out
real in_1(20,1)
in_1(1:20,1) = 0.e0 ! Vector of zeros
in_1(4:9,1) = 1.e0 ! Put some non-zero elements on it
m = size(in_1,dim = 1)
n = size(in_1,dim = 2)
call findf_r(out,in_1,m,n) ! Here it gives me an error
write(*,*) size(out)
end program finder
subroutine findf_r(out,in_1,m,n)
real, dimension(m,n) :: in_1
integer, dimension(size(in_1)) :: temp
integer, pointer, dimension(:,:) :: out
integer in_1_m, in_1_n
integer i, j, num
in_1_m = size(in_1,dim = 1)
in_1_n = size(in_1,dim = 2)
num = 0
if ((in_1_m == 1).and.(in_1_n == 1)) then
! IN THIS CASE IS A SCALAR, NOT VECTOR NOR MATRIX
if (in_1(1,1) /= 0) then
num = 1
allocate(out(num,1))
out(1,1) = 1
else
allocate(out(num,1))
endif
elseif (in_1_m == 1) then
! IN THIS CASE IS A ROW VECTOR
do i = 1, in_1_n
if (in_1(1,i) /= 0) then
num = num + 1
temp(num) = i
endif
enddo
allocate(out(1,num))
out(1,1:num) = temp
elseif (in_1_n == 1) then
! IN THIS CASE IS A COLUMN VECTOR
do i = 1, in_1_m
if (in_1(i,1) /= 0) then
num = num + 1
temp(num) = i
endif
enddo
allocate(out(num,1))
out(1:num,1) = temp
else
! IN THIS CASE IS A MATRIX
do j = 1, in_1_n
do i = 1, in_1_m
if (in_1(i,j) /= 0) then
num = num + 1
temp(num) = i + (j - 1)*in_1_m
endif
enddo
enddo
allocate(out(num,1))
out(1:num,1) = temp
endif
return
end subroutine findf_r
!!!!!!!!!!!!!!!!!!!! CODE END !!!!!!!!!!!!!!!!!!!!!
Thanks for any suggestion.
Andrea.
Jugoslav Dujic
| Hello Fortran Gurus,
|
| I'm not an expert with fortran, but I would like to ask 2
| questions:
|
| 1) Is there a fortran function/subroutine that returns the indices of
| a vector (matrix) X that are non-zero? (for whom that know matlab,
| like the "find" function).
|
| 2) If this is not the case, could you please tell me what is wrong in
| this routine? When debugging, it still keeps telling me that "out" is
| an "undefined pointer/array" while it's not...
1) is not the case. I didn't look deeper into the contents of findf_r,
but the basic mistake is that any routine which has a POINTER argument
must have an explicit interface, which is not the case with your code.
As it is now, PROGRAM doesn't "know" that findf_r expects a POINTER as
1st argument, so it assumes that it expects a normal array, so it
"sends" the target of out, which is undefined at that moment (probably
NULLified). Your code should crash with a Access violation/Segmentation
fault, doesn't it?
To achieve an explicit interface, you can:
1) (best) place findf_r in a MODULE and USE that module within main
program, or
2) make findf_r a contained routine of the program (move end program
to the very end and place a CONTAINS statement instead). Note that
findf_r will "inherit" all the local variables from the program, or
3) put an INTERFACE block to findf_r in the PROGRAM (complicated &
error prone).
Oh, and yes, use IMPLICIT NONE :-).
--
Jugoslav
___________
www.geocities.com/jdujic
Dave Seaman
> Hello Fortran Gurus,
> I'm not an expert with fortran, but I would like to ask 2
> questions:
> 1) Is there a fortran function/subroutine that returns the indices of
> a vector (matrix) X that are non-zero? (for whom that know matlab,
> like the "find" function).
In what form do you need the indices? Will a logical array help? There are
lots of ways to use such arrays in vector operations.
program bittest
implicit none
integer, parameter :: n = 5
real :: a(n) = (/ 1.0, 0.0, 2.0, 0.0, 0.0 /)
logical :: flag(n)
flag = a /= 0.0 ! the "/=" operator is sometimes spelled ".ne."
print *, flag
end program bittest
This program prints
T F T F F
--
Dave Seaman
Judge Yohn's mistakes revealed in Mumia Abu-Jamal ruling.
<http://www.commoncouragepress.com/index.cfm?action=book&bookid=228>
David Ham
> On 10 Sep 2003 04:22:55 -0700, Andrea wrote:
> > Hello Fortran Gurus,
>
> > I'm not an expert with fortran, but I would like to ask 2
> > questions:
>
> > 1) Is there a fortran function/subroutine that returns the indices of
> > a vector (matrix) X that are non-zero? (for whom that know matlab,
> > like the "find" function).
>
> In what form do you need the indices? Will a logical array help? There are
> lots of ways to use such arrays in vector operations.
>
> program bittest
>
> implicit none
>
> integer, parameter :: n = 5
> real :: a(n) = (/ 1.0, 0.0, 2.0, 0.0, 0.0 /)
> logical :: flag(n)
>
> flag = a /= 0.0 ! the "/=" operator is sometimes spelled ".ne."
>
> print *, flag
>
> end program bittest
>
> This program prints
>
> T F T F F
>
This idea can be pretty easily turned into a function to do what the
OP wants:
function find(vector)
!--------------------------------------------------
! Return a pointer to an array containing the
! indices of the non-zero elements of vector.
!
! The usual warnings about pointer-valued functions
! and memory leaks apply.
!
! Remember that since vector is an assumed shape
! array, it is always indexed starting from 1.
!--------------------------------------------------
implicit none
real, dimension(:), intent(in) :: vector
integer, dimension(:), pointer :: find
integer :: i, j
logical, dimension(size(vector)) :: flag
flag = vector/=0
allocate(find(count(flag)))
j=1
do i=1, size(vector)
if (flag(i)) then
find(j)=i
j=j+1
end if
end do
end function find
If you don't fancy mucking around with pointers, replace the
declaration of find with:
integer, dimension(size(vector)) :: find
and the allocate statement with:
find=0
Then the function returns a vector the same size as vector with the
indices of non-zero entries first and trailing zeros. Note that zero
is an unamiguously out of band value.
David
Paul van Delst
>
> function find(vector)
> implicit none
>
> real, dimension(:), intent(in) :: vector
> integer, dimension(:), pointer :: find
>
> integer :: i, j
> logical, dimension(size(vector)) :: flag
>
> flag = vector/=0
>
> allocate(find(count(flag)))
find = PACK( (/(i,i=1,size(vector))/), flag )
> end function find
Would this work with a pointer argument? I use PACK a lot for finding indices based upon
some mask condition.
cheers,
paulv
--
Paul van Delst
CIMSS @ NOAA/NCEP/EMC
Ph: (301)763-8000 x7748
Fax:(301)763-8545
Michel OLAGNON
What about
pack ( (/ (i, i=1,size(array)) /), (array /= 0) )
Richard Maine
> find = PACK( (/(i,i=1,size(vector))/), flag )
...
> Would this work with a pointer argument? I use PACK a lot for
> finding indices based upon some mask condition.
Yes. Why wouldn't it?
I presume you are referring to "find" being a pointer. Note that
it is *NOT* an argument to PACK (though it is a dummy argument in
the function from which this was extracted). Indeed, "find" here
basically has nothing to do with PACK. The "find" is just the
left hand side of an assignment statement. The details of how
the right-hand-side is constructed are quite irelevant; all that
matters is that whatever expression is on the right-hand-size
evaluates to something appropriate to assign to the left-hand-side.
In this case, the right-hand-side evaluates to a rank 1 array of
type default integer. The previous line allocated "find" to
an appropriate size. Thus the left-hand-side is a definable variable
and is the same type, kind, and shape as the expression. I don't see
anything to question.
Remember that a fundamental point about Fortran pointers is that they
are automatically dereferenced. In most places, the name of the
pointer refers to its target; no explicit dereference operator is
required. Special syntax is required to deal with the pointer itself
as opposed to the target. You can use "find" pretty much anywhere
that you could use a nonpointer array of the same type, kind, and
shape. You can write expressions like find+find, and you can use
"find" on the left-hand-side of assignment statements. This is an
ordinary assignment statement, *NOT* a pointer assignment.
--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain