Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

  • Congratulations strongm on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

Need some help to understand fortran code

Status
Not open for further replies.

NateRiver

Programmer
Jun 30, 2014
1
US
Right now I am translating program written in Fortran to Matlab, and I encounter 2 problems I really need help with:
1. I saw these lines of codes:
Code:
c     OPEN (UNIT=41,FILE=TRIM(TMPDIR)//LFNAME//'GN.DAT', STATUS='OLD')
c     DO 11 K=1,INT,1
c     READ (41,*) XXLIGHT(K),LIGHTSC(K)
c     11            CONTINUE
c     CLOSE (41)
I know that a capital C at the eginning make a Fortran code a comment, but I'm not quite sure if it is also true with lowercase c.

2. I saw two subroutine, which are ...
Code:
	SUBROUTINE DGBFA(ABD,LDA,N,ML,MU,IPVT,INFO)
	INTEGER LDA,N,ML,MU,IPVT(1),INFO,I,ISAMAX,IO,J,
     *		JU,JZ,JO,J1,K,KP1,L,LM,M,MM,NM1
	DOUBLE PRECISION ABD(LDA,1),T
C
	INCLUDE './misc/const.inc.f'
C
D     WRITE(6,*) ' Start:  DGBFA'
C
	M=ML+MU+1
	INFO=0
	JO=MU+2
	J1=MIN0(N,M)-1
	IF (J1.LT.JO) GO TO 30
	DO 20  JZ=JO,J1
		IO=M+1-JZ
		DO 10  I=IO,ML
			ABD(I,JZ)=ZERO
 10		CONTINUE
 20	CONTINUE
 30	CONTINUE
	JZ=J1
	JU=0
	NM1=N-1
	IF (NM1.LT.1) GO TO 130
	DO 120  K=1,NM1
		KP1=K+1
		JZ=JZ+1
		IF (JZ.GT.N) GO TO 50
		IF (ML.LT.1) GO TO 50
		DO 40  I=1,ML
			ABD(I,JZ)=ZERO
 40		CONTINUE
 50		CONTINUE
		LM=MIN0(ML,N-K)
		L=ISAMAX(LM+1,ABD(M,K),1)+M-1
		IPVT(K)=L+K-M
		IF (ABD(L,K).EQ.ZERO) GO TO 100
		IF (L.EQ.M) GO TO 60
		T=ABD(L,K)
		ABD(L,K)=ABD(M,K)
		ABD(M,K)=T
 60		CONTINUE
		T=-ONE/ABD(M,K)
		CALL DSCAL(LM,T,ABD(M+1,K),1)
		JU=MIN0(MAX0(JU,MU+IPVT(K)),N)
		MM=M
		IF (JU.LT.KP1) GO TO 90
		DO 80  J=KP1,JU
			L=L-1
			MM=MM-1
			T=ABD(L,J)
			IF (L.EQ.MM) GO TO 70
			ABD(L,J)=ABD(MM,J)
			ABD(MM,J)=T
 70			CONTINUE
			CALL DAXPY(LM,T,ABD(M+1,K),1,ABD(MM+1,J),1)
 80		CONTINUE
 90		CONTINUE
		GO TO 110
 100		CONTINUE
		INFO=K
 110		CONTINUE
 120	CONTINUE
 130	CONTINUE
	IPVT(N)=N
	IF (ABD(M,N).EQ.ZERO) INFO=N
D     WRITE(6,*) ' Stop:   DGBFA'
	RETURN
	END
and...
Code:
SUBROUTINE DAXPY(N,SA,SX,INCX,SY,INCY)
	DOUBLE PRECISION SX(1),SY(1),SA
	INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
	INCLUDE './misc/const.inc.f'
C
	IF (N.LE.0) RETURN
	IF (SA.EQ.ZERO) RETURN
	IF (INCX.EQ.1.AND.INCY.EQ.1) GO TO 20
	IX=1
	IY=1
	IF (INCX.LT.0) IX=(-N+1)*INCX+1
	IF (INCY.LT.0) IY=(-N+1)*INCY+1
	DO 10 I=1,N
		SY(IY)=SY(IY)+SA*SX(IX)
		IX=IX+INCX
		IY=IY+INCY
 10	CONTINUE
	RETURN
 20	M=MOD(N,4)
	IF (M.EQ.0) GO TO 40
	DO 30  I=1,M
		SY(I)=SY(I)+SA*SX(I)
 30	CONTINUE
	IF (N.LT.4) RETURN
 40	MP1=M+1
	DO 50  I=MP1,N,4
		SY(I)=SY(I)+SA*SX(I)
		SY(I+1)=SY(I+1)+SA*SX(I+1)
		SY(I+2)=SY(I+2)+SA*SX(I+2)
		SY(I+3)=SY(I+3)+SA*SX(I+3)
 50	CONTINUE
	RETURN
	END
C

When subroutine DGBFA calls subroutine DAXPY, I dont understand the way it pass the values of array ABD (a 2D array) to SX and SY (a 1D array). My guess is that it pass only one column or row of ABD to SX and SY, but even so, I am not sure if it uses the whole column/row or just a portion of that column/row.

I havent been able to find the answers for these problem anywhere, so I am really appreciated if someone here can help me with them.
 

do you have a way to compile the fortran program with added WRITE statements and do some inspection as to what is doing? That would be your best answer.

In any case, your question is a very good one.

Engineers programming in Fortran have had to resort to a few "tricks" here and there, some legit tricks, some abusing the condition of the language from before there was a legit way of doing something.

What you see here is a misleading declaration style of an Assumed Size Array; I say misleading because it would seem you are declaring a 1-dimensional array that is only 1 item long...but that is not necessarily the case, although if you compile this program for debugging it may certainly report out-of-bounds error.

Looking at the array entry being passed, the first thought would be that a single array entry (scalar) is being passed; but the fact that they have declared the receiving parameter to be an (assumed size) array says otherwise.

Memory storage of full matrices come in two flavors row-major and column-major. C/C++, for example, store row-major; Fortran, on the other hand, stores column-major...same as Matlab, by the way, probably because linear algebra in Matlab is done with LAPACK which in turn is written in Fortran 90.

So, it seems to me that arrays SX and SY will be traversing the matrix ABD one column at a time, from top to bottom and from left to right starting at the address given (Fortran is pass by reference unless said otherwise, rarely done).

 
Status
Not open for further replies.

Part and Inventory Search

Sponsor

Back
Top