Browse Source

first commit

mkilani 1 year ago
commit
1cf25b5a2a
10 changed files with 805 additions and 0 deletions
  1. 0 0
      README.md
  2. 15 0
      TES-LICENCE
  3. 79 0
      input.f
  4. 120 0
      lq.f
  5. 187 0
      maintes.f
  6. 48 0
      matcoef.f
  7. 40 0
      matinv.f
  8. 232 0
      quick_sort2.f
  9. 68 0
      regional.f
  10. 16 0
      structure.txt

+ 0 - 0
README.md


+ 15 - 0
TES-LICENCE

@@ -0,0 +1,15 @@
+    This file is part of TES.
+
+    TES is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    TES is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with TES.  If not, see <http://www.gnu.org/licenses/>.
+

+ 79 - 0
input.f

@@ -0,0 +1,79 @@
+      SUBROUTINE INPUTDATA(FILENAME,NOBS,TES,EMPL,EMPG,PROD,
+     &            REGIONAME,DEBUG,METHOD,DELTA,LABELS,CODES,SORT)
+      IMPLICIT NONE
+      INTEGER NOBS,DEBUG
+      CHARACTER*25 FILENAME
+      CHARACTER*20 REGIONAME
+      CHARACTER*40 LABELS(NOBS)
+      CHARACTER*5  CODES(NOBS)
+      CHARACTER*5 METHOD
+      CHARACTER*6 SORT
+      REAL*8 TES(NOBS,NOBS),EMPL(NOBS),EMPG(NOBS),PROD(NOBS)
+      REAL*8 DELTA
+      REAL*8 EMPLTOT,EMPGTOT
+      CHARACTER*40 LINE
+      CHARACTER*10 FMT
+      INTEGER N,I,J
+
+C     LECTURE DES DONNÉES
+      OPEN(UNIT=1,FILE=FILENAME)
+  1   READ(1,*) LINE
+      N=SCAN(LINE,':',.FALSE.)
+      IF(N.EQ.0) GOTO 1 
+      LINE=LINE(1:N-1)
+
+      IF(LINE.EQ.'NOBS') THEN
+         READ(1,*) NOBS
+      ELSEIF(LINE.EQ.'TES') THEN
+         READ(1,*) ((TES(I,J),J=1,NOBS),I=1,NOBS)
+      ELSEIF(LINE.EQ.'EMP') THEN
+         DO I=1,NOBS
+            READ(1,*) CODES(I),EMPL(I),EMPG(I),PROD(I)
+         ENDDO
+      ELSEIF(LINE.EQ.'EMP-PROP') THEN
+         READ(1,*) EMPLTOT,EMPGTOT
+         DO I=1,NOBS
+            READ(1,*) CODES(I),EMPL(I),EMPG(I),PROD(I)
+            EMPL(I)=EMPLTOT*EMPL(I)/100.D0
+            EMPG(I)=EMPGTOT*EMPG(I)/100.D0
+         ENDDO
+      ELSEIF(LINE.EQ.'DELTA') THEN
+         READ(1,*) DELTA
+      ELSEIF(LINE.EQ.'LABELS') THEN
+         DO I=1,NOBS
+            READ(1,'(A40)') LABELS(I)
+         ENDDO
+      ELSEIF(LINE.EQ.'METHOD') THEN
+         READ(1,*) METHOD 
+      ELSEIF(LINE.EQ.'SORTING')THEN
+         READ(1,*) SORT
+      ELSEIF(LINE.EQ.'NAME') THEN
+         READ(1,*) REGIONAME
+      ELSEIF(LINE.EQ.'DEBUG')THEN
+         READ(1,*) DEBUG
+      ELSEIF(LINE.EQ.'END') THEN 
+         GOTO 100
+      ENDIF 
+      GOTO 1
+
+ 100  CONTINUE
+      CLOSE(1) 
+      IF(DEBUG.EQ.1) THEN
+        WRITE(FMT,'("("I0,"F12.3)")') NOBS
+        PRINT*,'NOBS=',NOBS
+        PRINT*,'CODES: ',(CODES(I),I=1,NOBS)
+        PRINT*, 'EMPL'
+        WRITE(6,FMT),(EMPL(I),I=1,NOBS)
+        PRINT*,'EMPG'
+        WRITE(6,FMT),(EMPG(I),I=1,NOBS)
+
+        WRITE(FMT,'("("I0,"F9.3)")') NOBS
+        PRINT*,'TES'
+        DO I=1,NOBS
+          WRITE(*,FMT),(TES(I,J),J=1,NOBS)
+        ENDDO
+        PRINT*,'PROD'
+        WRITE(6,FMT),(PROD(I),I=1,NOBS)
+      ENDIF
+      RETURN
+      ENDSUBROUTINE

+ 120 - 0
lq.f

@@ -0,0 +1,120 @@
+      SUBROUTINE LOCQ(METHOD,EMPLOC,EMPGEN,NOBS,NLD,DELTA,MATLQ,LQ,
+     +                                                  DEBUG)
+      IMPLICIT NONE
+C
+C     Author: Moez Kilani
+C     Date: March 2012
+C     Contact: moez [dot] kilani [at] univ-lille3.fr
+C
+C
+C     INPUT
+C        METHOD:  the method used to compute the coefficient
+C                 matrix. There are four methods accepted :
+C                  SLQ, CILQ, RLQ, FLQ
+C        EMPLOC:  local employment data, an vector of NOBS
+C                 elements
+C        EMPGEN:  global employment data (at nation level, in general)
+C        NOBS:    number of sectors considered.
+C     OUTPUT
+C        MATLQ:   the coefficient matrix
+C
+      INTEGER NOBS,NLD,DEBUG
+      REAL*8 EMPLOC(*),EMPGEN(*),MATLQ(NLD,*)
+      CHARACTER*5 METHOD,METHODS(5)
+      REAL*8 TOTLOC,TOTGEN,LQ(*),DELTA,LAMBDA
+      CHARACTER*10 FMT
+
+      INTEGER I,J
+
+      DATA METHODS /'SLQ','CILQ','RLQ','FLQ','NAT'/
+
+      TOTLOC=0.D0
+      TOTGEN=0.D0
+      DO I=1,NOBS
+         TOTLOC=TOTLOC+EMPLOC(I)
+         TOTGEN=TOTGEN+EMPGEN(I)
+      ENDDO
+      DO I=1,NOBS
+         LQ(I)=(EMPLOC(I)/TOTLOC)/(EMPGEN(I)/TOTGEN)
+      ENDDO
+
+      IF(METHOD.EQ.METHODS(1)) THEN
+         GOTO 1   !  method SLQ
+      ELSEIF(METHOD.EQ.METHODS(2)) THEN
+         GOTO 2   ! method CILQ
+      ELSEIF(METHOD.EQ.METHODS(3)) THEN
+         GOTO 3   ! method RLQ
+      ELSEIF(METHOD.EQ.METHODS(4)) THEN
+         GOTO 4   ! method FLQ
+      ELSEIF(METHOD.EQ.METHODS(5)) THEN
+         GOTO 5   ! method NAT
+      ENDIF
+
+C     If none of the methods is matched then print
+C     an error message here and go to the end
+      PRINT*, 'ERROR: Method ', METHOD, ' is not implemented'
+      GOTO 100
+
+   1  CONTINUE
+C     METHOD 1: SLQ
+      DO I=1,NOBS
+         DO J=1,NOBS
+            MATLQ(I,J)=MIN(1.D0,LQ(I))
+         ENDDO
+      ENDDO
+      GOTO 100
+
+
+   2  CONTINUE
+C     METHOD 2: CILQ
+      DO I=1,NOBS
+         DO J=1,NOBS
+            MATLQ(I,J)=MIN(1.D0,LQ(I)/LQ(J))
+         ENDDO
+      ENDDO
+      GOTO 100
+
+
+   3  CONTINUE
+C     METHOD 3: RLQ 
+      DO I=1,NOBS
+         DO J=1,NOBS
+            MATLQ(I,J)=MIN(1.D0,LQ(I)/(LOG(1.D0+LQ(J))/LOG(2.D0)))
+         ENDDO
+      ENDDO
+      GOTO 100
+
+   4  CONTINUE
+C     METHOD 4: FLQ
+      LAMBDA=(LOG(1.D0+TOTLOC/TOTGEN)/LOG(2.D0))**DELTA
+      DO I=1,NOBS
+         DO J=1,NOBS
+            MATLQ(I,J)=MIN(1.D0,LAMBDA*LQ(I)/LQ(J))
+         ENDDO
+      ENDDO
+      DO I=1,NOBS
+         MATLQ(I,I)=MIN(1.D0,LAMBDA*LQ(I))
+      ENDDO
+      GOTO 100 
+
+   5  CONTINUE
+C     METHOD 5: NAT
+      DO I=1,NOBS
+         DO J=1,NOBS
+            MATLQ(I,J)=1.0D0
+         ENDDO
+      ENDDO
+      GOTO 100 
+
+ 100  CONTINUE 
+
+      IF(DEBUG.EQ.1)THEN
+        WRITE(FMT,'("("I0,"F9.3)")') NOBS
+        PRINT*,'LOCATION QUOTIENT MATRIX'
+        DO I=1,NOBS
+          WRITE(*,FMT),(MATLQ(I,J),J=1,NOBS)
+        ENDDO
+      ENDIF
+
+      RETURN
+      ENDSUBROUTINE

+ 187 - 0
maintes.f

@@ -0,0 +1,187 @@
+      PROGRAM REGIONTES
+      IMPLICIT NONE
+      INTEGER N,NOBS,NARGS,NLQ
+      PARAMETER(N=1000)
+      REAL*8 EMPL(N),EMPG(N),PROD(N),PRODL(N)
+      REAL*8 TES(N,N),COEF(N,N),LQ(N)
+      REAL*8 MATLQ(N,N),LEONTIEF(N,N),MULTIPLIERS(N)
+      REAL*8 DELTA,BMULT,BEMPL
+      INTEGER ORDM(N),DEBUG,LMAT
+      CHARACTER*25 FILENAME
+      CHARACTER*20 REGIONAME
+      CHARACTER*40 LABELS(N)
+      CHARACTER*5 CODES(N) 
+      CHARACTER*5 METHOD
+      CHARACTER*6 SORT
+      CHARACTER*10 ARGS
+      CHARACTER*10 FMT
+      REAL*8 EMPL09(N),EMPL30(9),EMPL1(N),EMPL2(N)
+      REAL*8 PRODL09(N),PRODL30(N),PRODL1(N),PRODL2(N)
+      REAL*8 F09(N),F30(N),F1(N),F2(N) 
+      INTEGER I,J 
+********************************
+      REAL*8 F(N),Y(N),WORK(N),Z(N),RCOND,LL(N,N),L(N,N)
+      REAL*8 DET, NEWPRODL(N)
+      INTEGER IPVT(N)
+
+      NARGS=IARGC()
+      IF(NARGS.LT.1) THEN
+        PRINT*,'YOU SHOULD PROVIDE THE INPUT FILE NAME'
+        STOP
+      ENDIF 
+      CALL GETARG(1,FILENAME)
+
+      NOBS=N
+      CALL INPUTDATA(FILENAME,NOBS,TES,EMPL,EMPG,PROD,
+     &       REGIONAME,DEBUG,METHOD,DELTA,LABELS,CODES,SORT) 
+C     PRINT*,'step 1'
+C     WRITE(*,'(36F9.2)') TES(25,1:NOBS)
+
+C The command line options overwrite those 
+C given in the input file
+
+      IF(NARGS.GT.1) THEN
+      DO 10 I=1,NARGS
+         CALL GETARG(I,ARGS)
+         IF(ARGS.EQ.'-L=YES') THEN
+            LMAT=1
+         ELSEIF(ARGS.EQ.'-L=NO') THEN
+            LMAT=0
+         ELSEIF(ARGS.EQ.'-d=YES') THEN
+            DEBUG=1
+         ELSEIF(ARGS.EQ.'-d=NO') THEN
+            DEBUG=0
+         ELSEIF(ARGS.EQ.'-s=LQ') THEN
+            SORT='LQ'
+         ELSEIF(ARGS.EQ.'-s=MULT') THEN
+            SORT='METHOD'
+         ELSEIF(ARGS.EQ.'-s=NONE') THEN
+            SORT='NONE'
+         ELSEIF(ARGS.EQ.'-m=SLQ') THEN
+            METHOD='SLQ'
+         ELSEIF(ARGS.EQ.'-m=CILQ') THEN
+            METHOD='CILQ'
+         ELSEIF(ARGS.EQ.'-m=RLQ') THEN
+            METHOD='RLQ'
+         ELSEIF(ARGS.EQ.'-m=CILQ') THEN
+            METHOD='CILQ'
+         ELSEIF(ARGS.EQ.'-m=FLQ') THEN
+            METHOD='FLQ'
+         ELSEIF(ARGS.EQ.'-m=NAT') THEN
+            METHOD='NAT'
+         ENDIF
+ 10   CONTINUE
+      ENDIF
+
+      CALL MATCOEF(TES,PROD,COEF,NOBS,N,DEBUG)
+      CALL LOCQ(METHOD,EMPL,EMPG,NOBS,N,DELTA,MATLQ,LQ,DEBUG)
+C     PRINT*,'step 3'
+      CALL REGIONAL(N,NOBS,COEF,MATLQ,LEONTIEF,MULTIPLIERS,DEBUG)
+C     PRINT*,COEF(25,1:NOBS)
+C     PRINT*,'step 4'
+      NLQ=0
+      DO I=1,NOBS
+         IF(LQ(I).GT.1.D0) THEN
+                 NLQ=NLQ+1
+                 BEMPL=BEMPL+EMPL(I)
+         ENDIF
+      ENDDO
+      BMULT=SUM(EMPL(1:NOBS))/BEMPL
+      PRINT*, '======================================'
+      PRINT*, 'REGION/CITY NAME: ', REGIONAME
+      PRINT*, '--------------------------------------'
+      PRINT '(X,A31,I3)',   'NUMBER OF SECTORS WITH LQ > 1 :', NLQ
+      PRINT*, '--------------------------------------'
+      PRINT '(X,A31,F7.3)', 'BASE MULTIPLIER               :', BMULT
+      PRINT*, '--------------------------------------'
+      IF(SORT.EQ.'LQ') THEN
+         CALL SORTRX(NOBS,LQ,ORDM)
+      ELSEIF(SORT.EQ.'METHOD')THEN
+         CALL SORTRX(NOBS,MULTIPLIERS,ORDM)
+      ELSE
+         DO I=1,NOBS
+            ORDM(I)=NOBS-I+1
+         ENDDO
+      ENDIF
+      PRINT'(18X,A6)',METHOD                 
+      PRINT*,' CODES   LQ       MULTIP.   DESCR.'
+      PRINT*,'======================================'
+      DO I=1,NOBS
+         J=ORDM(NOBS-I+1)
+         PRINT '((5X)(A4)(f6.3)(3X)(F6.3)(5X)(A35))', 
+     &           CODES(J), LQ(J),MULTIPLIERS(J),LABELS(J)
+      ENDDO 
+      PRINT*,'======================================'
+      PRINT* 
+      IF(LMAT.EQ.1) THEN
+        WRITE(FMT,'("("I0,"F9.3)")') NOBS
+        PRINT*,'LEONTIEF MATRIX (I-A)^(-1)'
+        DO I=1,NOBS
+          WRITE(*,FMT),(LEONTIEF(I,J),J=1,NOBS)
+        ENDDO
+      ENDIF
+      PRINT* 
+
+************************************************
+************************************************
+* This part is added for computations to the
+* case of IDF.
+************************************************
+************************************************
+C     DO I=1,NOBS
+C        PRODL(I)=PROD(I)/EMPG(I)*EMPL(I)
+*        PRINT '(A2,3F15.3)', CODES(I),PRODL(I),PROD(I),PROD(I)/EMPG(I)
+C     ENDDO
+* LL = (I-A) and
+* LEONTIEF = (I-A)^-1
+* PRODL : the initial production
+* EMPL : the initial employment
+* NEWPRODL : the new production
+* NEWEMPL : the new employment
+* 
+C     LL=LEONTIEF
+C     CALL DGECO(LL,N,NOBS,IPVT,RCOND,Z)
+C     CALL DGEDI(LL,N,NOBS,IPVT,DET,WORK,01)
+
+* Now LL is the inverse of LEONTIEF
+      CALL DGEMV('N',NOBS,NOBS,1.0D0,LL,N,PRODL,1,0.0D0,Y,1)
+* Y contains the final consumption
+ 
+*     PRINT*, 'Emploi final à l échelle locale (M€)'
+*     DO I=1,NOBS
+*         PRINT '(F12.3)', Y(I)
+*     ENDDO
+
+* here a test of whether the computation of 
+* the inverse matrix are correct 
+*     CALL DGEMM('N','N',NOBS,NOBS,NOBS,1.D0,LL,N,LEONTIEF,N,0.D0,L,N)
+
+* the transport sector has an increase of 10% in final demand
+*     PRINT '(F12.3)',Y(20)
+C     Y(20)=1.1D0*Y(20)
+*     PRINT '(F12.3)',Y(20)
+
+* the induced change in production is
+
+C     CALL DGEMV('N',NOBS,NOBS,1.D0,LEONTIEF,N,Y,1,0.0D0,NEWPRODL,1)
+
+C     DO I=1,NOBS
+C        PRINT '(A2,",",F12.3,",",F12.3)', CODES(I),LQ(I),
+C    &   (NEWPRODL(I)-PRODL(I))/(PROD(I)/EMPG(I))
+C     ENDDO 
+
+*     PRINT*, LEONTIEF(1,1:NOBS) 
+*     DO I=1,5
+*        PRINT '(5F12.3)', LL(I,1:5)
+*     ENDDO
+*     PRINT*
+*     DO I=1,5
+*        PRINT '(5F12.3)', LEONTIEF(I,1:5)
+*     ENDDO
+*     PRINT*
+*     DO I=1,5
+*        PRINT '(5F12.3)', L(I,1:5)
+*     ENDDO
+*     Y(20)=Y(20)*1.1D0
+
+      END

+ 48 - 0
matcoef.f

@@ -0,0 +1,48 @@
+      SUBROUTINE MATCOEF(TES,PROD,COEF,NOBS,NLD,DEBUG)
+* ***************************************************************
+*                                                               *
+* Summary:                                                      *
+* --------                                                      *
+* Transforms an initial input-output table                      *
+* to a table of technical coefficients                          *
+*                                                               *
+* Variables:                                                    *
+* ----------                                                    *
+* TES (input)   : (NOBSxNOBS)-matrix of the input-output table  *
+* PROD (input)  : (1xNOBS)-vector of production values          *
+* COEF (output) : the table of technical coefficients           *
+* NOBS (input)  : the number of obesrvations in the problem     *
+* NLD (input)   : the leading dimension of the tables           *
+* DEBUG (input) : a binary variable that indicates whether      *
+*                 intermediate computations are shown or not    *
+*                                                               *
+* Author : Moez Kilani                                          *
+* Date : April 2012; revised June 2012                          *
+*                                                               *
+*****************************************************************
+      IMPLICIT NONE
+      INTEGER NOBS,NLD,DEBUG
+      REAL*8 TES(NLD,*),PROD(*),COEF(NLD,*) 
+C     LOCAL VARIABLES
+      INTEGER I,J
+      REAL*8 TOT
+      CHARACTER*10 FMT
+
+      DO I=1,NOBS
+         DO J=1,NOBS
+            COEF(I,J)=TES(I,J)/PROD(J)
+         ENDDO
+      ENDDO
+
+
+      IF(DEBUG.EQ.1)THEN
+        WRITE(FMT,'("("I0,"F9.3)")') NOBS
+        PRINT*,'TECHNICAL COEFICIENTS'
+        DO I=1,NOBS
+          WRITE(*,FMT),(COEF(I,J),J=1,NOBS)
+        ENDDO
+      ENDIF
+
+      RETURN
+      ENDSUBROUTINE
+

+ 40 - 0
matinv.f

@@ -0,0 +1,40 @@
+      SUBROUTINE MATINV(A,NLD,NB)
+C     ======================================================
+C     COMPUTES THE INVERSE OF A MATRIX A
+C
+C     PARAMETERS:
+C     -----------
+C       A (INPUT/OUTPUT): ON INPUT THE MATRIX TO INVERSE
+C                         AND ON OUTPUT THE RESULT
+C       NB (INPUT): THE DIMENSION OF MATRXI A IS NB x NB
+C       NLD: LEADING DIMENSION OF A
+C
+C     NOTE : USES LAPACK SUBROUTINES "DGETRF" AND "DGETRI"      
+C     ======================================================
+      INTEGER IPIV,INFO,LDA
+      REAL*8 WORK
+      INTEGER I,J,NB
+      REAL*8 A(NLD,*)
+      PARAMETER(IWORK=1000)
+      DIMENSION WORK(IWORK)
+      DIMENSION IPIV(NB)
+      LDA=NB
+C     ****************
+C     LU decomposition
+C     DGETRF : a lapack subroutine
+C     ****************
+      CALL DGETRF(NB,NB,A,NLD,IPIV,INFO)
+      IF (INFO.NE.0) THEN
+          PRINT *, "Echec dans la décomposition LU"
+      ENDIF
+C     *****************************************
+C     Invert from LU decomposition
+C     DGETRI: a lapack subroutine
+C     *****************************************
+      CALL DGETRI(NB,A,NLD,IPIV,WORK,IWORK,INFO)
+      IF (INFO.NE.0) THEN
+          PRINT *, "Echec dans l'inversion de la matrice"
+      ENDIF
+      RETURN
+      ENDSUBROUTINE
+

+ 232 - 0
quick_sort2.f

@@ -0,0 +1,232 @@
+C From Leonard J. Moss of SLAC:
+
+C Here's a hybrid QuickSort I wrote a number of years ago.  It's
+C based on suggestions in Knuth, Volume 3, and performs much better
+C than a pure QuickSort on short or partially ordered input arrays.  
+
+      SUBROUTINE SORTRX(N,DATA,INDEX)
+C===================================================================
+C
+C     SORTRX -- SORT, Real input, indeX output
+C
+C
+C     Input:  N     INTEGER
+C             DATA  REAL
+C
+C     Output: INDEX INTEGER (DIMENSION N)
+C
+C This routine performs an in-memory sort of the first N elements of
+C array DATA, returning into array INDEX the indices of elements of
+C DATA arranged in ascending order.  Thus,
+C
+C    DATA(INDEX(1)) will be the smallest number in array DATA;
+C    DATA(INDEX(N)) will be the largest number in DATA.
+C
+C The original data is not physically rearranged.  The original order
+C of equal input values is not necessarily preserved.
+C
+C===================================================================
+C
+C SORTRX uses a hybrid QuickSort algorithm, based on several
+C suggestions in Knuth, Volume 3, Section 5.2.2.  In particular, the
+C "pivot key" [my term] for dividing each subsequence is chosen to be
+C the median of the first, last, and middle values of the subsequence;
+C and the QuickSort is cut off when a subsequence has 9 or fewer
+C elements, and a straight insertion sort of the entire array is done
+C at the end.  The result is comparable to a pure insertion sort for
+C very short arrays, and very fast for very large arrays (of order 12
+C micro-sec/element on the 3081K for arrays of 10K elements).  It is
+C also not subject to the poor performance of the pure QuickSort on
+C partially ordered data.
+C
+C Created:  15 Jul 1986  Len Moss
+C
+C===================================================================
+ 
+      INTEGER   N,INDEX(N)
+      REAL*8    DATA(N)
+ 
+      INTEGER   LSTK(31),RSTK(31),ISTK
+      INTEGER   L,R,I,J,P,INDEXP,INDEXT
+      REAL*8    DATAP
+ 
+C     QuickSort Cutoff
+C
+C     Quit QuickSort-ing when a subsequence contains M or fewer
+C     elements and finish off at end with straight insertion sort.
+C     According to Knuth, V.3, the optimum value of M is around 9.
+ 
+      INTEGER   M
+      PARAMETER (M=9)
+ 
+C===================================================================
+C
+C     Make initial guess for INDEX
+ 
+      DO 50 I=1,N
+         INDEX(I)=I
+   50    CONTINUE
+ 
+C     If array is short, skip QuickSort and go directly to
+C     the straight insertion sort.
+ 
+      IF (N.LE.M) GOTO 900
+ 
+C===================================================================
+C
+C     QuickSort
+C
+C     The "Qn:"s correspond roughly to steps in Algorithm Q,
+C     Knuth, V.3, PP.116-117, modified to select the median
+C     of the first, last, and middle elements as the "pivot
+C     key" (in Knuth's notation, "K").  Also modified to leave
+C     data in place and produce an INDEX array.  To simplify
+C     comments, let DATA[I]=DATA(INDEX(I)).
+ 
+C Q1: Initialize
+      ISTK=0
+      L=1
+      R=N
+ 
+  200 CONTINUE
+ 
+C Q2: Sort the subsequence DATA[L]..DATA[R].
+C
+C     At this point, DATA[l] <= DATA[m] <= DATA[r] for all l < L,
+C     r > R, and L <= m <= R.  (First time through, there is no
+C     DATA for l < L or r > R.)
+ 
+      I=L
+      J=R
+ 
+C Q2.5: Select pivot key
+C
+C     Let the pivot, P, be the midpoint of this subsequence,
+C     P=(L+R)/2; then rearrange INDEX(L), INDEX(P), and INDEX(R)
+C     so the corresponding DATA values are in increasing order.
+C     The pivot key, DATAP, is then DATA[P].
+ 
+      P=(L+R)/2
+      INDEXP=INDEX(P)
+      DATAP=DATA(INDEXP)
+ 
+      IF (DATA(INDEX(L)) .GT. DATAP) THEN
+         INDEX(P)=INDEX(L)
+         INDEX(L)=INDEXP
+         INDEXP=INDEX(P)
+         DATAP=DATA(INDEXP)
+      ENDIF
+ 
+      IF (DATAP .GT. DATA(INDEX(R))) THEN
+         IF (DATA(INDEX(L)) .GT. DATA(INDEX(R))) THEN
+            INDEX(P)=INDEX(L)
+            INDEX(L)=INDEX(R)
+         ELSE
+            INDEX(P)=INDEX(R)
+         ENDIF
+         INDEX(R)=INDEXP
+         INDEXP=INDEX(P)
+         DATAP=DATA(INDEXP)
+      ENDIF
+ 
+C     Now we swap values between the right and left sides and/or
+C     move DATAP until all smaller values are on the left and all
+C     larger values are on the right.  Neither the left or right
+C     side will be internally ordered yet; however, DATAP will be
+C     in its final position.
+ 
+  300 CONTINUE
+ 
+C Q3: Search for datum on left >= DATAP
+C
+C     At this point, DATA[L] <= DATAP.  We can therefore start scanning
+C     up from L, looking for a value >= DATAP (this scan is guaranteed
+C     to terminate since we initially placed DATAP near the middle of
+C     the subsequence).
+ 
+         I=I+1
+         IF (DATA(INDEX(I)).LT.DATAP) GOTO 300
+ 
+  400 CONTINUE
+ 
+C Q4: Search for datum on right <= DATAP
+C
+C     At this point, DATA[R] >= DATAP.  We can therefore start scanning
+C     down from R, looking for a value <= DATAP (this scan is guaranteed
+C     to terminate since we initially placed DATAP near the middle of
+C     the subsequence).
+ 
+         J=J-1
+         IF (DATA(INDEX(J)).GT.DATAP) GOTO 400
+ 
+C Q5: Have the two scans collided?
+ 
+      IF (I.LT.J) THEN
+ 
+C Q6: No, interchange DATA[I] <--> DATA[J] and continue
+ 
+         INDEXT=INDEX(I)
+         INDEX(I)=INDEX(J)
+         INDEX(J)=INDEXT
+         GOTO 300
+      ELSE
+ 
+C Q7: Yes, select next subsequence to sort
+C
+C     At this point, I >= J and DATA[l] <= DATA[I] == DATAP <= DATA[r],
+C     for all L <= l < I and J < r <= R.  If both subsequences are
+C     more than M elements long, push the longer one on the stack and
+C     go back to QuickSort the shorter; if only one is more than M
+C     elements long, go back and QuickSort it; otherwise, pop a
+C     subsequence off the stack and QuickSort it.
+ 
+         IF (R-J .GE. I-L .AND. I-L .GT. M) THEN
+            ISTK=ISTK+1
+            LSTK(ISTK)=J+1
+            RSTK(ISTK)=R
+            R=I-1
+         ELSE IF (I-L .GT. R-J .AND. R-J .GT. M) THEN
+            ISTK=ISTK+1
+            LSTK(ISTK)=L
+            RSTK(ISTK)=I-1
+            L=J+1
+         ELSE IF (R-J .GT. M) THEN
+            L=J+1
+         ELSE IF (I-L .GT. M) THEN
+            R=I-1
+         ELSE
+C Q8: Pop the stack, or terminate QuickSort if empty
+            IF (ISTK.LT.1) GOTO 900
+            L=LSTK(ISTK)
+            R=RSTK(ISTK)
+            ISTK=ISTK-1
+         ENDIF
+         GOTO 200
+      ENDIF
+ 
+  900 CONTINUE
+ 
+C===================================================================
+C
+C Q9: Straight Insertion sort
+ 
+      DO 950 I=2,N
+         IF (DATA(INDEX(I-1)) .GT. DATA(INDEX(I))) THEN
+            INDEXP=INDEX(I)
+            DATAP=DATA(INDEXP)
+            P=I-1
+  920       CONTINUE
+               INDEX(P+1) = INDEX(P)
+               P=P-1
+               IF (P.GT.0) THEN
+                  IF (DATA(INDEX(P)).GT.DATAP) GOTO 920
+               ENDIF
+            INDEX(P+1) = INDEXP
+         ENDIF
+  950    CONTINUE
+ 
+C===================================================================
+C
+C     All done
+ 
+      END

+ 68 - 0
regional.f

@@ -0,0 +1,68 @@
+      SUBROUTINE REGIONAL(NLD,NOBS,COEF,MATLQ,LEONTIEF,MULTIPLIERS,
+     +                                                     DEBUG)
+C     ------------------- 
+C     TAKES A MATRIX "A" AND PRODUCES A LEONTIEF MATRIX
+C     AND COMPUTES MULTIPLIERS
+      IMPLICIT NONE
+      INTEGER NOBS,NLD,DEBUG
+      REAL*8 COEF(NLD,*),MATLQ(NLD,*),LEONTIEF(NLD,*)
+      REAL*8 MULTIPLIERS(NOBS)
+      CHARACTER*10 FMT
+
+C     LOCAL DATA
+      INTEGER I,J
+      REAL*8 MATID(NOBS,NOBS), MATTMP(NOBS,NOBS)
+
+      DO I=1,NOBS
+        DO J=1,NOBS
+           MATID(I,J)=0.D0
+        ENDDO
+      ENDDO
+      DO I=1,NOBS
+         MATID(I,I)=1.D0
+      ENDDO 
+
+      DO 10 I=1,NOBS
+      DO 10 J=1,NOBS
+         COEF(I,J)=COEF(I,J)*MATLQ(I,J)
+         LEONTIEF(I,J)=MATID(I,J)-COEF(I,J)
+
+  10  CONTINUE
+
+      IF(DEBUG.EQ.1)THEN
+        WRITE(FMT,'("("I0,"F9.3)")') NOBS
+        PRINT*,'MATRIX A'
+        DO I=1,NOBS
+          WRITE(*,FMT),(COEF(I,J),J=1,NOBS)
+        ENDDO
+
+        PRINT*,'MATRIX I-A'
+        DO I=1,NOBS
+          WRITE(*,FMT),(LEONTIEF(I,J),J=1,NOBS)
+        ENDDO
+      ENDIF
+
+      CALL MATINV(LEONTIEF,NLD,NOBS)
+
+      IF(DEBUG.EQ.1)THEN
+        WRITE(FMT,'("("I0,"F9.3)")') NOBS
+        PRINT*,'LEONTIEF MATRIX (I-A)^(-1)'
+        DO I=1,NOBS
+          WRITE(*,FMT),(LEONTIEF(I,J),J=1,NOBS)
+        ENDDO
+      ENDIF
+
+      MULTIPLIERS=0.D0
+      DO J=1,NOBS
+         DO I=1,NOBS
+            MULTIPLIERS(J)=MULTIPLIERS(J)+LEONTIEF(I,J)
+         ENDDO
+      ENDDO
+
+
+      RETURN
+      ENDSUBROUTINE
+
+
+
+

+ 16 - 0
structure.txt

@@ -0,0 +1,16 @@
+
+
+rtes
+ |______main
+ |______lq
+ |______regional
+ |______matinv
+ |______input
+
+
+to add :
+- input manipulation
+- multiple regions in single file
+- separate input files (main script, national tes tables, labour distribution)
+- ras technique
+- docmentation