Barrowes Consulting
Fortran consulting: code modernization, refactoring.
Matlab consulting: code development, GUIs, physics/electrical engineering
Specializing in fortran => Matlab conversion/translation/porting (see below)
The fortran goto statement allows jumps within a fortran program. Most goto's jump to the same or shallower nested if/do/while level. Jumps into a more deeply nested if/do/while were rare.
To have any hope of converting or translating fortran code to a more modern language like matlab/octave, python, etc., one has to remove every goto in the fortran code. I do this using while/break/continue combinations with a tool called remgoto.m I wrote. Remgoto analyzes the fortran code for goto's, then refactors the goto's by replacing them with do/while loops and branch statements. The resulting fortran code is logically and functionally equivalent to the original fortran code, but without the goto's. After this step, the fortran code can now be more readily converted into matlab/octave source code.
As an example, consider the following subroutine cbal.f from slatec:
subroutine CBAL (NM, N, AR, AI, LOW, IGH, SCALE) INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC REAL AR(NM,*),AI(NM,*),SCALE(*) REAL C,F,G,R,S,B2,RADIX LOGICAL NOCONV RADIX = 16 B2 = RADIX * RADIX K = 1 L = N go to 100 20 SCALE(M) = J if (J == M) go to 50 DO 30 I = 1, L F = AR(I,J) AR(I,J) = AR(I,M) AR(I,M) = F F = AI(I,J) AI(I,J) = AI(I,M) AI(I,M) = F 30 CONTINUE DO 40 I = K, N F = AR(J,I) AR(J,I) = AR(M,I) AR(M,I) = F F = AI(J,I) AI(J,I) = AI(M,I) AI(M,I) = F 40 CONTINUE 50 go to (80,130), IEXC 80 if (L == 1) go to 280 L = L - 1 100 DO 120 JJ = 1, L J = L + 1 - JJ DO 110 I = 1, L if (I == J) go to 110 if (AR(J,I) /= 0.0E0 .OR. AI(J,I) /= 0.0E0) go to 120 110 CONTINUE M = L IEXC = 1 go to 20 120 CONTINUE go to 140 130 K = K + 1 140 DO 170 J = K, L DO 150 I = K, L if (I == J) go to 150 if (AR(I,J) /= 0.0E0 .OR. AI(I,J) /= 0.0E0) go to 170 150 CONTINUE M = K IEXC = 2 go to 20 170 CONTINUE DO 180 I = K, L 180 SCALE(I) = 1.0E0 190 NOCONV = .FALSE. DO 270 I = K, L C = 0.0E0 R = 0.0E0 DO 200 J = K, L if (J == I) go to 200 C = C + ABS(AR(J,I)) + ABS(AI(J,I)) R = R + ABS(AR(I,J)) + ABS(AI(I,J)) 200 CONTINUE if (C == 0.0E0 .OR. R == 0.0E0) go to 270 G = R / RADIX F = 1.0E0 S = C + R 210 if (C >= G) go to 220 F = F * RADIX C = C * B2 go to 210 220 G = R * RADIX 230 if (C < G) go to 240 F = F / RADIX C = C / B2 go to 230 240 if ((C + R) / F >= 0.95E0 * S) go to 270 G = 1.0E0 / F SCALE(I) = SCALE(I) * F NOCONV = .TRUE. DO 250 J = K, N AR(I,J) = AR(I,J) * G AI(I,J) = AI(I,J) * G 250 CONTINUE DO 260 J = 1, L AR(J,I) = AR(J,I) * F AI(J,I) = AI(J,I) * F 260 CONTINUE 270 CONTINUE if (NOCONV) go to 190 280 LOW = K IGH = L return end
This has 19 goto statements along with being older style f77. Even with a program like Polyhedron's spag, we still get the following:
SUBROUTINE CBAL(NM, N, AR, AI, LOW, IGH, SCALE) IMPLICIT NONE INTEGER , INTENT(IN) :: NM INTEGER , INTENT(IN) :: N INTEGER , INTENT(OUT) :: LOW INTEGER , INTENT(OUT) :: IGH REAL , INTENT(INOUT) :: AR(NM,*) REAL , INTENT(INOUT) :: AI(NM,*) REAL , INTENT(INOUT) :: SCALE(*) INTEGER :: I, J, K, L, M, JJ, IEXC REAL :: C, F, G, R, S, B2, RADIX LOGICAL :: NOCONV RADIX = 16 B2 = RADIX*RADIX K = 1 L = N GO TO 100 20 SCALE(M) = J IF (J /= M) THEN DO I = 1, L F = AR(I,J) AR(I,J) = AR(I,M) AR(I,M) = F F = AI(I,J) AI(I,J) = AI(I,M) AI(I,M) = F END DO DO I = K, N F = AR(J,I) AR(J,I) = AR(M,I) AR(M,I) = F F = AI(J,I) AI(J,I) = AI(M,I) AI(M,I) = F END DO ENDIF IF (IEXC == 1) GO TO 80 IF (IEXC == 2) GO TO 130 80 IF (L == 1) GO TO 280 L = L - 1 100 L120: DO JJ = 1, L J = L + 1 - JJ DO I = 1, L IF (I == J) CYCLE IF (AR(J,I)==0.0E0 .AND. AI(J,I)==0.0E0) CYCLE CYCLE L120 END DO M = L IEXC = 1 GO TO 20 END DO L120 GO TO 140 130 K = K + 1 140 L170: DO J = K, L DO I = K, L IF (I == J) CYCLE IF (AR(I,J)==0.0E0 .AND. AI(I,J)==0.0E0) CYCLE CYCLE L170 END DO M = K IEXC = 2 GO TO 20 END DO L170 SCALE(K:L) = 1.0E0 190 NOCONV = .FALSE. DO I = K, L C = 0.0E0 R = 0.0E0 DO J = K, L IF (J == I) CYCLE C = C + ABS(AR(J,I)) + ABS(AI(J,I)) R = R + ABS(AR(I,J)) + ABS(AI(I,J)) END DO IF (C==0.0E0 .OR. R==0.0E0) CYCLE G = R/RADIX F = 1.0E0 S = C + R 210 IF (C >= G) GO TO 220 F = F*RADIX C = C*B2 GO TO 210 220 G = R*RADIX 230 IF (C < G) GO TO 240 F = F/RADIX C = C/B2 GO TO 230 240 IF ((C + R)/F >= 0.95E0*S) CYCLE G = 1.0E0/F SCALE(I) = SCALE(I)*F NOCONV = .TRUE. AR(I,K:N) = AR(I,K:N)*G AI(I,K:N) = AI(I,K:N)*G AR(:L,I) = AR(:L,I)*F AI(:L,I) = AI(:L,I)*F END DO IF (NOCONV) GO TO 190 280 LOW = K IGH = L RETURN END SUBROUTINE CBAL
This is more modern and looks better, but the routine still has 12 goto's. Using remgoto, we get the following fortran90 code with no goto's, but instead has do/exit/cycle/if statements:
SUBROUTINE CBAL(NM, N, AR, AI, LOW, IGH, SCALE) IMPLICIT NONE logical :: remg(15)=.true. INTEGER , INTENT(IN) :: NM INTEGER , INTENT(IN) :: N INTEGER , INTENT(OUT) :: LOW INTEGER , INTENT(OUT) :: IGH REAL , INTENT(INOUT) :: AR(NM,*) REAL , INTENT(INOUT) :: AI(NM,*) REAL , INTENT(INOUT) :: SCALE(*) INTEGER :: I, J, K, L, M, JJ, IEXC REAL :: C, F, G, R, S, B2, RADIX LOGICAL :: NOCONV do if (remg(11)) then if (remg(6)) then if (remg(5)) then if (remg(4)) then if (remg(3)) then if (remg(2)) then if (remg(1)) then RADIX = 16 B2 = RADIX*RADIX K = 1 L = N remg(3)=.false.;cycle endif remg(1)=.true. 20 SCALE(M) = J IF (J /= M) THEN DO I = 1, L F = AR(I,J) AR(I,J) = AR(I,M) AR(I,M) = F F = AI(I,J) AI(I,J) = AI(I,M) AI(I,M) = F enddo DO I = K, N F = AR(J,I) AR(J,I) = AR(M,I) AR(M,I) = F F = AI(J,I) AI(J,I) = AI(M,I) AI(M,I) = F enddo ENDIF IF (IEXC == 1) then remg(2)=.false.;cycle endif IF (IEXC == 2) then remg(4)=.false.;cycle endif endif remg(2)=.true. 80 IF (L == 1) then remg(11)=.false.;cycle endif L = L - 1 endif remg(3)=.true. 100 DO JJ = 1, L!L120 J = L + 1 - JJ DO I = 1, L IF (I == J) then cycle endif IF (AR(J,I)==0.0E0 .AND. AI(J,I)==0.0E0) then cycle endif remg(12)=.false.;exit enddo if (.not.(remg(12))) then remg([12])=.true. ; cycle endif M = L IEXC = 1 remg(1)=.false.;exit enddo !L120 if (.not.(remg(1))) cycle remg(5)=.false.;cycle endif remg(4)=.true. 130 K = K + 1 endif remg(5)=.true. 140 DO J = K, L!L170 DO I = K, L IF (I == J) then cycle endif IF (AR(I,J)==0.0E0 .AND. AI(I,J)==0.0E0) then cycle endif remg(13)=.false.;exit enddo if (.not.(remg(13))) then remg([13])=.true. ; cycle endif M = K IEXC = 2 remg(1)=.false.;exit enddo !L170 if (.not.(remg(1))) cycle SCALE(K:L) = 1.0E0 endif remg(6)=.true. 190 NOCONV = .FALSE. DO I = K, L do if (remg(10)) then if (remg(9)) then if (remg(8)) then if (remg(7)) then C = 0.0E0 R = 0.0E0 DO J = K, L IF (J == I) then cycle endif C = C + ABS(AR(J,I)) + ABS(AI(J,I)) R = R + ABS(AR(I,J)) + ABS(AI(I,J)) enddo IF (C==0.0E0 .OR. R==0.0E0) then remg(14)=.false.;exit endif G = R/RADIX F = 1.0E0 S = C + R endif remg(7)=.true. 210 IF (C >= G) then remg(8)=.false.;cycle endif F = F*RADIX C = C*B2 remg(7)=.false.;cycle endif remg(8)=.true. 220 G = R*RADIX endif remg(9)=.true. 230 IF (C < G) then remg(10)=.false.;cycle endif F = F/RADIX C = C/B2 remg(9)=.false.;cycle endif remg(10)=.true. 240 IF ((C + R)/F >= 0.95E0*S) then remg(15)=.false.;exit endif G = 1.0E0/F SCALE(I) = SCALE(I)*F NOCONV = .TRUE. AR(I,K:N) = AR(I,K:N)*G AI(I,K:N) = AI(I,K:N)*G AR(:L,I) = AR(:L,I)*F AI(:L,I) = AI(:L,I)*F exit enddo if (.not.(remg(14).and.remg(15))) then remg([14,15])=.true. ; cycle endif enddo IF (NOCONV) then remg(6)=.false.;cycle endif endif remg(11)=.true. 280 LOW = K IGH = L exit enddo RETURN END SUBROUTINE CBAL
This can now be converted into matlab/octave code with f2matlab:
function [nm, n, ar, ai, low, igh, scale_ml]=cbal(nm, n, ar, ai, low, igh, scale_ml); persistent b2 c f g i iexc j jj k l m noconv r radix_ml remg s ; if isempty(remg), remg([1:15])=true; end; ar_shape=size(ar);ar=reshape([ar(:).',zeros(1,ceil(numel(ar)./prod([abs(nm)])).*prod([abs(nm)])-numel(ar))],abs(nm),[]); ai_shape=size(ai);ai=reshape([ai(:).',zeros(1,ceil(numel(ai)./prod([abs(nm)])).*prod([abs(nm)])-numel(ai))],abs(nm),[]); scale_shape=size(scale_ml);scale_ml=reshape(scale_ml,1,[]); if isempty(i), i=0; end; if isempty(j), j=0; end; if isempty(k), k=0; end; if isempty(l), l=0; end; if isempty(m), m=0; end; if isempty(jj), jj=0; end; if isempty(iexc), iexc=0; end; if isempty(c), c=0; end; if isempty(f), f=0; end; if isempty(g), g=0; end; if isempty(r), r=0; end; if isempty(s), s=0; end; if isempty(b2), b2=0; end; if isempty(radix_ml), radix_ml=0; end; if isempty(noconv), noconv=false; end; while (1); if(remg(11)) if(remg(6)) if(remg(5)) if(remg(4)) if(remg(3)) if(remg(2)) if(remg(1)) radix_ml = 16; b2 = radix_ml.*radix_ml; k = 1; l = fix(n); remg(3)=false; continue; end; remg(1)=true; scale_ml(m) = j; if(j ~= m) for i = 1: l; f = ar(i,j); ar(i,j) = ar(i,m); ar(i,m) = f; f = ai(i,j); ai(i,j) = ai(i,m); ai(i,m) = f; end; i = fix(l+1); for i = k: n; f = ar(j,i); ar(j,i) = ar(m,i); ar(m,i) = f; f = ai(j,i); ai(j,i) = ai(m,i); ai(m,i) = f; end; i = fix(n+1); end; if(iexc == 1) remg(2)=false; continue; end; if(iexc == 2) remg(4)=false; continue; end; end; remg(2)=true; if(l == 1) remg(11)=false; continue; end; l = fix(l - 1); end; remg(3)=true; %L120 for jj = 1: l; j = fix(l + 1 - jj); for i = 1: l; if(i == j) continue; end; if(ar(j,i)==0.0e0 && ai(j,i)==0.0e0) continue; end; remg(12)=false; break; end; if(~(remg(12))) remg([12])=true; continue; end; m = fix(l); iexc = 1; remg(1)=false; break; %L120 end; if(~(remg(1))) continue; end; remg(5)=false; continue; end; remg(4)=true; k = fix(k + 1); end; remg(5)=true; %L170 for j = k: l; for i = k: l; if(i == j) continue; end; if(ar(i,j)==0.0e0 && ai(i,j)==0.0e0) continue; end; remg(13)=false; break; end; if(~(remg(13))) remg([13])=true; continue; end; m = fix(k); iexc = 2; remg(1)=false; break; %L170 end; if(~(remg(1))) continue; end; scale_ml([k:l]) = 1.0e0; end; remg(6)=true; noconv = false; for i = k: l; while (1); if(remg(10)) if(remg(9)) if(remg(8)) if(remg(7)) c = 0.0e0; r = 0.0e0; for j = k: l; if(j == i) continue; end; c = c + abs(ar(j,i)) + abs(ai(j,i)); r = r + abs(ar(i,j)) + abs(ai(i,j)); end; j = fix(l+1); if(c==0.0e0 || r==0.0e0) remg(14)=false; break; end; g = r./radix_ml; f = 1.0e0; s = c + r; end; remg(7)=true; if(c >= g) remg(8)=false; continue; end; f = f.*radix_ml; c = c.*b2; remg(7)=false; continue; end; remg(8)=true; g = r.*radix_ml; end; remg(9)=true; if(c < g) remg(10)=false; continue; end; f = f./radix_ml; c = c./b2; remg(9)=false; continue; end; remg(10)=true; if((c + r)./f >= 0.95e0.*s) remg(15)=false; break; end; g = 1.0e0./f; scale_ml(i) = scale_ml(i).*f; noconv = true; ar(i,[k:n]) = ar(i,[k:n]).*g; ai(i,[k:n]) = ai(i,[k:n]).*g; ar([1:l],i) = ar([1:l],i).*f; ai([1:l],i) = ai([1:l],i).*f; break; end; if(~(remg(14)&&remg(15))) remg([14,15])=true; continue; end; end; if(noconv) remg(6)=false; continue; end; end; remg(11)=true; low = fix(k); igh = fix(l); break; end; return; end %subroutine cbal
This converted matlab/octave code can now be optimized and modified as an m-file as desired.
My background is in physics and electrical engineering as per my CV, dissertation, and other publications.
References upon request.
Contact: Ben Barrowes, barrowes@alum.mit.edu
