Incorrect output from Fortran loop

Question

I have written a routine to give the column index position for the furthers right cell containing a 1, marking the right edge of a polygon in a mask array.

But when I print the index arrays (I realize the row position array is redundant), all rows are reporting a position (other than 0), which I know shouldn't be the case. I can't seem to find where the following could be incorrect.

Example mask array:

0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 0 0 0 0 0 0 
0 0 0 0 0 0 1 1 1 0 0 0 0 0 
0 0 0 0 0 1 1 1 1 1 1 0 0 0 
0 0 0 1 1 1 1 1 1 1 1 0 0 0 
0 0 0 1 1 1 1 1 1 1 1 0 0 0
0 0 0 0 0 0 1 1 1 1 0 0 0 0
0 0 0 0 0 0 0 0 0 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0

Desired output:

[0,10,10,11,11,11,9,8,0]

Is this a Fortran thing, or is my logic just off (

Routine:

  subroutine get_right_idx(mask_array, idx_x, idx_y)                        
  integer, parameter :: x = 169 ! num columns                               
  integer, parameter :: y = 124 ! num rows                                  
  integer i ! iterator for rows                                             
  integer j ! iterator for columns                                          
  integer row_x ! x position for furthest south cell in row                 
  integer row_y ! y position for furthest south cell in row                 
  integer :: idx_x(y) ! index positions for lowest lat in model grid - x    
  integer :: idx_y(y) ! index positions for lowest lat in model grid - y    
  real mask_array(y,x) ! mask array of zeros, containing polygon of ones    

  do j=1,y                                                                  
    row_x = 0                                                               
    row_y = 0                                                               
    do i=1,x                                                                
      if (mask_array(j,i).eq.1) then                                        
        row_x = i                                                           
        row_y = j                                                           
      endif                                                                 
    enddo                                                                   
    idx_x(j)=row_x                                                          
    idx_y(j)=row_y                                                          
  enddo                                                                     
  endsubroutine get_right_idx

Actual mask (zoomed out in Open Office): enter image description here

Below is the mask that I am trying to evaluate. My resulting array has a non-zero value in for all elements, where there should be zero elements at the start and end of the array, no matter which direction it is evaluated from.

Output:

125         104         104         104         104         104         104         114         114         114         114         103         103         103         108         108         103         103         103         103          97          97          97         107         107         107         107         107         107         107         107         107          97         101         101         101         101         101         111         111         111         111         111         111         101         101         100         105         105         105         105         105         105         100         115         115         104         104         104         104         104         104         104         104         104         104          98          98          98          98         108         108         108         108         108         108         108         108          98         102         102         102         102         102         112         112         112         112         112         112         101         101         101         106         106         106         101         101         101          95          95         105         105         105         105         105         105         105         105         105         105          99          99          99          99          99         109         109         109         109         109         109         109          99          99

If you want the highest index with 1, why not start from the top and work downwards with an exit command? — Kyle Kanos
– Kyle Kanos, Commented Jan 9, 2014 at 18:49
Hi Kyle, thanks, I'm not seeing how that would work. I think I did a poor job of explaining the sort of weird thing I'm trying to do, but I've given an example array and output to help. — ryanjdillon
– ryanjdillon, Commented Jan 9, 2014 at 19:03

Kyle Kanos · Accepted Answer · 2014-01-09 19:17:51Z

2

I made a quick template of your code, but with my suggestion. Since you want to find the largest index j that has a 1 in it, you can simply test this by starting your j-loop with the maximum (in the example, 14) and working downwards but exiting the inner do-loop if mask_array==1.

You could generalize this by replacing 14 and 9 with the maximum values for those dimensions.

program get_right
   implicit none
   integer, dimension(9,14) :: mask_array
   integer :: i,j,mask_out(9),j_tmp
   mask_array(1,:)=[0,0,0,0,0,0,0,0,0,0,0,0,0,0]
   mask_array(2,:)=[0,0,0,0,0,0,0,1,0,0,0,0,0,0]
   mask_array(3,:)=[0,0,0,0,0,0,1,1,1,0,0,0,0,0]
   mask_array(4,:)=[0,0,0,0,0,1,1,1,1,1,1,0,0,0]
   mask_array(5,:)=[0,0,0,1,1,1,1,1,1,1,1,0,0,0]
   mask_array(6,:)=[0,0,0,1,1,1,1,1,1,1,1,0,0,0]
   mask_array(7,:)=[0,0,0,0,0,0,1,1,1,1,0,0,0,0]
   mask_array(8,:)=[0,0,0,0,0,0,0,0,0,1,0,0,0,0]
   mask_array(9,:)=[0,0,0,0,0,0,0,0,0,0,0,0,0,0]

   mask_out=0
   do i=1,9
      j_tmp=0
      print '(14(i0,2x))',mask_array(i,:)
      do j=14,1,-1
         if(mask_array(i,j)==1) then
            j_tmp=j
            exit
         endif
      enddo
      mask_out(i)=j
   enddo

   print *,""
   print '(9(i0,2x))',mask_out

end program get_right

I got 0,8,9,11,11,11,10,10,0 as the answer (which is what you got but backwards).

answered Jan 9, 2014 at 19:17

Kyle Kanos

3,2752 gold badges25 silver badges38 bronze badges

Sign up to request clarification or add additional context in comments.

4 Comments

ryanjdillon Over a year ago

Interesting. I actually didn't try running my example array, but have been using another larger array as input. I'll take another look at this tomorrow. Thank you Kyle.

ryanjdillon Over a year ago

This seems perfectly logical, and is just about the same as mine, but I continue to get values from my mask array that are incorrect. For rows where I know no mask value (i.e. 1) exist, I'm getting something. Uff.

Kyle Kanos Over a year ago

@shootingstars: it is possible that you have your 2D array mapped differently in your head than in the code. Fortran is column major and you might be thinking your array as row major (see this wiki article).

ryanjdillon Over a year ago

It doesn't look like it Kyle. Not sure what is going on. I updated my post to give a visual of the mask I am working with and the output generated.

High Performance Mark · Accepted Answer · 2014-01-12 09:27:32Z

2

Or, if typing makes you tired you could evaluate this

maxval(mask_array*spread([(ix,ix=lbound(mask_array,dim=2),&
        ubound(mask_array,dim=2))],dim=1, ncopies=size(mask_array,dim=1)),dim=2)

As you can see this makes a temporary array (which may be undesirable if your masks are large) using spread and an implied-do loop, each element of this temporary contains its own column index number. Then it multiplies the temporary with the mask_array, performing element-wise multiplication not matrix multiplication. Finally take the maxval of each row of the result. This returns the vector in the same order as Kyle's code does.

I've edited the code to use lbound and ubound rather than 1 and size in case you want to use the code on arrays with lower bounds other than 1.

Bob's yer uncle, but don't ask whether this is faster than Kyle's code. If you are interested in execution speed test and measure.

Incidentally, since this returns just one vector and doesn't modify its arguments or have any other side effects, I'd package it as a function rather than as a subroutine.

edited Jan 12, 2014 at 9:27

answered Jan 9, 2014 at 21:55

High Performance Mark

78.6k7 gold badges109 silver badges168 bronze badges

1 Comment

ryanjdillon Over a year ago

This is a neat approach. I tend to stick to things that I can immediately grasp, so at the moment I'm still leaning on my original method. Thanks!

Collectives™ on Stack Overflow

Incorrect output from Fortran loop

2 Answers 2

4 Comments

1 Comment

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

4 Comments

1 Comment

Your Answer

Sign up or log in

Post as a guest

Related