How to export matrixes in Fortran?

[Ilario_D_]

Hi! I am tryng to generate a matrix(m x n), MatContribution, with the following subroutine. Its elements are 1 with a certain probability p, 0 otherwise.

    SUBROUTINE Matrice_Conoscenza(p,n,m,ConoscenzaAgente,MatContribution)
    REAL,INTENT(IN) :: p
    INTEGER,INTENT(IN) :: n,m
    REAL, DIMENSION(1:n), INTENT(INOUT) :: ConoscenzaAgente
    REAL, DIMENSION(1:m,1:n), INTENT(OUT) :: MatContribution
    INTEGER :: i,j
    DO i=1,m
        CALL RANDOM_SEED()
        CALL RANDOM_NUMBER(ConoscenzaAgente)
        DO j=1,n
             IF(ConoscenzaAgente(j)<p) THEN 
                 MatContribution(i,j)=1  
             ELSE 
                 MatContribution(i,j)=0
             ENDIF
        ENDDO
    ENDDO
    ENDSUBROUTINE Matrice_Conoscenza

 I call the subroutine in the main program by the following statements and export MAT_CONTRIBUTION.

CALL Matrice_Conoscenza(P,N,M,CONOSCENZA_AGENTE,MAT_CONTRIBUTION)
   OPEN(UNIT=0,FILE='MatContribution.DAT')
   DO i=1,M
       WRITE (0,"(100000(f2.0))") MAT_CONTRIBUTION(i,:)
   ENDDO
   CLOSE(UNIT=0)

The problem is that the MAT_CONTRIBUTION I obtain has all rows equals.

0.0.0.0.1.0.1.0.1.0.0.0.
0.0.0.0.1.0.1.0.1.0.0.0.
0.0.0.0.1.0.1.0.1.0.0.0.
0.0.0.0.1.0.1.0.1.0.0.0.
0.0.0.0.1.0.1.0.1.0.0.0.
0.0.0.0.1.0.1.0.1.0.0.0.
0.0.0.0.1.0.1.0.1.0.0.0.

Where are my errors in the code? Is it correct to export a matrix in this way? Thank you in advance.

[Arjen_Markus]

The problem is in the use of the routine random_seed. It resets the random number generator, but how it does that depends on the compiler. The Intel Fortran compiler may in fact be using the same seed at each invocation if you do not pass a seed. So the obvious solution is to leave out the call to random_seed. First of all it is not necessary, as subsequent calls to random_number generate new random numbers, secondly it may cause the kind of behaviour you are observing.

Just a few more remarks on the code:

You pass the dimensions of the array as separate arguments. This is not necessary, use assumed-shape arrays along these lines:

subroutine mysub( p, array )
    real, dimension(:) :: array

    do i = 1,size(array)
          write(*,*) i, array(i)
   enddo
end subroutine mysub

The size of the array (or size(matrix,1), size(matrix,2) for the individual dimensions) gives the number of elements in the array (or the number of rows and columns).

Using logical unit numbers like 0 is not recommended practice - there are a few magic number, usually 5 and 6 for input from the keyboard and output to the terminal, and 0 for "standard error". Use LU-numbers from 10 onwards. Or better still, use the NEWUNIT= keyword.

[Steven_L_Intel1]

arjenmarkus is correct about RANDOM_SEED. Our implementation is that if you don't call RANDOM_SEED, the same seed is used each time you run the program. If you call it with no arguments, it resets the seed based on the time-of-day, but this behavior is not specified by the standard.

Fortran 2015 adds a new RANDOM_INIT intrinsic to let you specify how RANDOM_NUMBER is initialized.

In a typical program, you call RANDOM_SEED just once at the start of the program. You DON'T call it before each call to RANDOM_NUMBER!

[andrew_4619]

following from Arjen's comments and noting also Steve's contribution. If the matrices are large some F2008 might optimise better and look a little more compact.

MODULE test
implicit none
contains
SUBROUTINE test_Matrice_Conoscenza()
    implicit none
    integer :: istat
    real    :: p = 0.4
    real    :: MatContribution(5,5)
    call Matrice_Conoscenza(p,MatContribution,istat)
END SUBROUTINE test_Matrice_Conoscenza
SUBROUTINE Matrice_Conoscenza(p,MatContribution, istat)
    implicit none
    integer, intent(out) :: istat
    REAL,INTENT(IN)      :: p
    REAL, INTENT(OUT)    :: MatContribution(:,:)
    real, allocatable    :: tmp(:,:)
    allocate(tmp, mold=MatContribution, STAT = ISTAT)
    if( istat /= 0) then
        MatContribution = 0.
        return
    endif
    CALL RANDOM_NUMBER(tmp)
    where ( tmp < p ) 
        MatContribution = 1.
    elsewhere
        MatContribution = 0.
    endwhere
end SUBROUTINE Matrice_Conoscenza
END MODULE test

[dboggs]

Re Arjen Markus #2:

I can't completely agree with your advice to use an assumed shape array instead of passing dimensions. Mainly, because you don't mention that an explicit interface is required to use assumed-shape arrays, and Ilario is liable to think that doing this is more trouble than passing the dimensions. I happen to agree.

One major reason to use explicit interface is to utilize such things as optional arguments, keyword arguments, generic functions and other advanced argument things; these are not involved in this case. Another major reason relates to a variety of security, portability, and compiler checking abilities. I get this. But, for simplicity and clarity, whenever the dimensions are conveniently available in the calling procedure, I still prefer to pass them to the subroutine.

Call me old fashioned, and I'm sure there are readers here who disagree, and I would love to hear some discussion on this matter.

[mecej4]

In general, whether to pass array sizes as extra arguments or depend on the properties of assumed shape array arguments: It Depends!

In simple cases, where only one or two array arguments are passed, passing them explicitly may be fine. Example: finding the weighted scalar product of two vectors: Real Function wDot(x,y,w,n). However, when two or more rectangular (or triangular) matrices are being passed, and their shapes have to conform, the F77 way can quickly get confusing and error-prone.

Consider the original code of this thread. The matrix is of size m X n, yet the argument list contains n first and then m. This discrepancy is obvious in this instance, but may not be so obvious in more complex calls. For most algorithms, such array size arguments should be INTENT(IN), but the old-fashioned way of writing code exposes them to inadvertent change. Secondly, human readers tend to assign affinity between array arguments and array size arguments based on relative placement in the argument list. Not respecting this expectation when writing code creates problems.

A large number of routines in BLAS and Lapack involve this issue. I think that BLAS95 and Lapack95 provide simplified interfaces to rather complex subroutines, which can be really beneficial to casual users. Consider this example from a recent MKL Forum thread:

F77: 
   CALL DHSEQR('E','N',N,ILO,IHI,H,LDH,WR,WI,Z,LDZ,WORK,LWORK,INFO)
   CALL DHSEIN('R','Q','N',SELECT,N,H,LDH,WR,WI,VL,LDVL,VR,LDVR,MM,M,WORK,IFAILL,IFAILR,INFO)
F95:
   hsav = h
   call hseqr(h,wr,wi,info=info)
   call hsein(hsav,wr,wi,select,vr=vr,eigsrc='n',info=info)

The user's problem was traced to forgetting to save the argument H, which is overwritten in DHSEQR. Consider the large number of integer arguments in the F77 version. With that version, it takes much more reading to find out what to pass in all those arguments, even when many of the arguments may never be used and are merely place-holders.

To sum up, "it depends", but for me the new way works better most often. In fact, and this may be an indication of my habits, I rewrote the OP's F77 code in that thread (https://software.intel.com/en-us/forums/intel-math-kernel-library/topic/606702) using F95 calls after I had given up on finding the bug in the F77 version. 

Andrew_4619's code in #4 is clean and very readable, and conveys the programmer's intent very clearly with the WHERE block.

[dboggs]

Thank you for these comments mecej4. This is valuable to me.

It should also be noted that the code in #4, where SUBROUTINE Matrice_Conoscenza uses an assumed-shape array argument , works only because it is contained in a module (thus deriving an explicit interface). Casual readers may miss this. It is a significant style change for older F77 programmers to put all their subroutines inside a module. It takes some work and dedication to appreciate the benefits.

[mecej4]

During the process of upgrading older codes to use assumed-shape array arguments, the /warn:interfaces compiler option can be very helpful. For example, with the following version (which is in error!) of the code of #1 and #4, using that option tells the user about the missing interface.

PROGRAM test_Matrice_Conoscenza
    implicit none
    integer :: i,istat
    real    :: p = 0.4
    real    :: MatContribution(5,5)

    call Matrice_Conoscenza(p,MatContribution,istat)
    open(unit=10,file='mat.dat',status='replace')
    write(10,'(5E12.4)')(matContribution(i,:),i=1,5)
    close(10)
    
END PROGRAM test_Matrice_Conoscenza

SUBROUTINE Matrice_Conoscenza(p,MatContribution, istat)
    implicit none
    integer, intent(out) :: istat
    REAL,INTENT(IN)      :: p
    REAL, INTENT(OUT)    :: MatContribution(:,:)
    real, allocatable    :: tmp(:,:)
    allocate(tmp, mold=MatContribution, STAT = ISTAT)
    if( istat /= 0) then
        MatContribution = 0.
        return
    endif
    CALL RANDOM_NUMBER(tmp)
    where ( tmp < p ) 
        MatContribution = 1.
    elsewhere
        MatContribution = 0.
    endwhere
end SUBROUTINE Matrice_Conoscenza

Without the interface checking, the program ran and output just zeros into the file. Recompiling with the warn option gave

$ ifort -warn interfaces matr.f90 
matr.f90(7): error #7978: Required interface for passing assumed shape array is missing from original source.   [MATCONTRIBUTION]
    call Matrice_Conoscenza(p,MatContribution,istat)
------------------------------^
compilation aborted for matr.f90 (code 1)