- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Dear all, I wrote a fortran code that calculates the ith-permutation of a given list {1,2,3,...,n}, without computing all the others, that are n! I needed that in order to find the ith-path of the TSP (Travelling salesman problem). When n! is big, the code gives me some error and I tested that the ith-permutation found is not the exact value. For n=10, there are not problems at all, but for n=20, the code crashes or wrong values are found. I think this is due to errors that Fortran makes operating with big numbers (sums of big numbers). I use Visual Fortran Ultimate 2013. In attached you find the subroutine I use for my goal. WeightAdjMatRete is the distance matrix between each pair of knots of the network. I hope someone can help me. Thank you in advance.
! Fattoriale RECURSIVE FUNCTION factorial(n) RESULT(n_factorial) IMPLICIT NONE REAL, INTENT(IN) :: n REAL :: n_factorial IF(n>0) THEN n_factorial=n*factorial(n-1) ELSE n_factorial=1. ENDIF ENDFUNCTION factorial ! ith-permutazione di una lista SUBROUTINE ith_permutazione(lista_iniziale,n,i,ith_permutation) IMPLICIT NONE INTEGER :: k,n REAL :: j,f REAL, INTENT(IN) :: i INTEGER, DIMENSION(1:n), INTENT(IN) :: lista_iniziale INTEGER, DIMENSION(1:n) :: lista_lavoro INTEGER, DIMENSION(1:n), INTENT(OUT) :: ith_permutation lista_lavoro=lista_iniziale j=i DO k=1,n f=factorial(REAL(n-k)) ith_permutation(k)=lista_lavoro(FLOOR(j/f)+1) lista_lavoro=PACK(lista_lavoro,MASK=lista_lavoro/=ith_permutation(k)) j=MOD(j,f) ENDDO ENDSUBROUTINE ith_permutazione ! Funzione modulo, adattata PURE FUNCTION mood(k,modulo) RESULT(ris) IMPLICIT NONE INTEGER, INTENT(IN) :: k,modulo INTEGER :: ris IF(MOD(k,modulo)/=0) THEN ris=MOD(k,modulo) ELSE ris=modulo ENDIF ENDFUNCTION mood ! Funzione quoziente, adattata PURE FUNCTION quoziente(a,p) RESULT(ris) IMPLICIT NONE INTEGER, INTENT(IN) :: a,p INTEGER :: ris IF(MOD(a,p)/=0) THEN ris=(a/p)+1 ELSE ris=a/p ENDIF ENDFUNCTION quoziente ! Vettori contenenti tutti i payoff percepiti dagli agenti allo state vector attuale e quelli ad ogni sua singola permutazione SUBROUTINE tuttipayoff(n,m,nodi,nodi_rete,sigma,bvector,MatVecSomma,VecPos,lista_iniziale,ith_permutation,lunghezze_percorso,WeightAdjMatRete,array_perceived_payoff_old,array_perceived_payoff_neg) IMPLICIT NONE INTEGER, INTENT(IN) :: n,m,nodi,nodi_rete INTEGER, DIMENSION(1:nodi), INTENT(IN) :: sigma INTEGER, DIMENSION(1:nodi), INTENT(OUT) :: bvector REAL, DIMENSION(1:m,1:n), INTENT(OUT) :: MatVecSomma REAL, DIMENSION(1:m), INTENT(OUT) :: VecPos INTEGER, DIMENSION(1:nodi_rete), INTENT(IN) :: lista_iniziale INTEGER, DIMENSION(1:nodi_rete), INTENT(OUT) :: ith_permutation REAL, DIMENSION(1:nodi_rete), INTENT(OUT) :: lunghezze_percorso REAL, DIMENSION(1:nodi_rete,1:nodi_rete), INTENT(IN) :: WeightAdjMatRete REAL, DIMENSION(1:nodi), INTENT(OUT) :: array_perceived_payoff_old,array_perceived_payoff_neg INTEGER :: i,j,k bvector=sigma FORALL(i=1:nodi,bvector(i)==-1) bvector(i)=0 ENDFORALL FORALL(i=1:m,j=1:n) MatVecSomma(i,j)=bvector(m*(j-1)+i)*(2.**REAL(n-j)) ENDFORALL FORALL(i=1:m) VecPos(i)=1.+SUM(MatVecSomma(i,:)) ENDFORALL DO k=1,nodi IF(VecPos(mood(k,m))<=factorial(REAL(nodi_rete))) THEN CALL ith_permutazione(lista_iniziale,nodi_rete,VecPos(mood(k,m))-1.,ith_permutation) FORALL(i=1:(nodi_rete-1)) lunghezze_percorso(i)=WeightAdjMatRete(ith_permutation(i),ith_permutation(i+1)) ENDFORALL lunghezze_percorso(nodi_rete)=WeightAdjMatRete(ith_permutation(nodi_rete),ith_permutation(1)) array_perceived_payoff_old(k)=(1./SUM(lunghezze_percorso)) ELSE array_perceived_payoff_old(k)=0. ENDIF IF(VecPos(mood(k,m))-SIGN(1,sigma(m*(quoziente(k,m)-1)+mood(k,m)))*2**(n-quoziente(k,m))<=factorial(REAL(nodi_rete))) THEN CALL ith_permutazione(lista_iniziale,nodi_rete,VecPos(mood(k,m))-SIGN(1,sigma(m*(quoziente(k,m)-1)+mood(k,m)))*2**(n-quoziente(k,m))-1.,ith_permutation) FORALL(i=1:(nodi_rete-1)) lunghezze_percorso(i)=WeightAdjMatRete(ith_permutation(i),ith_permutation(i+1)) ENDFORALL lunghezze_percorso(nodi_rete)=WeightAdjMatRete(ith_permutation(nodi_rete),ith_permutation(1)) array_perceived_payoff_neg(k)=(1./SUM(lunghezze_percorso)) ELSE array_perceived_payoff_neg(k)=0. ENDIF ENDDO ENDSUBROUTINE tuttipayoff
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
You can increase the precision of your REAL from the default of REAL(4) to REAL(8) or REAL(16).
If (when) you exceed the precision of the mantissa (fraction) portion of the given precision, rounding or truncation occurs. For n=20 (20!) you will require a precision of 19 digits.
REAL(4) 7.22 decimal digits (fraction)
REAL(8) 15.95 decimal digits (fraction)
REAL(16) 34.02 decimal digits (fraction)
Note, not all languages support REAL(16), and those that do, may not support the same internal format (meaning the precision may differ).
Jim Dempsey
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I solved my problem, thanks to your answer and to M Oehm advise, working by integers of KIND=18.
This is the link where the problem is exposed and fixed http://stackoverflow.com/questions/36136542/operations-with-big-real-numbers-in-fortran
I have not problems when NODI_RETE=12 or less, but have problem for bigger values. I attached my code .f90 here, written in visual studio ultimate 2013, hoping someone can help me to run the program also for NODI_RETE=20, that is my goal. For modify NODI_RETE go to row 482 and change 12 to 20.
All the files needed are attached. Thanks in advance.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Made some minor changes. Haven't verified if code is completely correct. (for example if the REAL's calculating results below need to be changed from REAL(4) to REAL(8)).
With NODI_RETE=20
Nodi_20_d_1.000_MinDist_1.000_MaxDist_10.000_Rete_1_M_7_T0_0.022_R_1 Payoff=0.321 err=3209.407 Consenso=0.162 err=1618.930 Payoff=0.314 err=0.021 Consenso=0.187 err=0.157 Payoff=0.318 err=0.013 Consenso=0.212 err=0.131 Payoff=0.310 err=0.025 Consenso=0.257 err=0.211 Payoff=0.320 err=0.030 Consenso=0.290 err=0.129 Payoff=0.311 err=0.026 Consenso=0.364 err=0.255 Payoff=0.320 err=0.028 Consenso=0.423 err=0.163 Payoff=0.332 err=0.038 Consenso=0.486 err=0.148 Payoff=0.337 err=0.015 Consenso=0.588 err=0.211 Payoff=0.322 err=0.044 Consenso=0.624 err=0.061 Payoff=0.332 err=0.031 Consenso=0.734 err=0.175 Payoff=0.341 err=0.027 Consenso=0.762 err=0.038 Payoff=0.475 err=0.393 Consenso=0.937 err=0.230 Payoff=0.485 err=0.022 Consenso=0.972 err=0.037 Payoff=0.511 err=0.053 Consenso=0.991 err=0.020 Payoff=0.512 err=0.002 Consenso=0.992 err=0.001 Percorso ottimale individuato: 17 4 15 18 13 12 10 2 14 7 16 9 6 5 19 8 20 3 1 11 con lunghezza complessiva di 66.957 contro quella minima di 34.373 Il blocco di simulazioni richieste e' terminato! Vuoi proseguire? Scrivi: S per SI N per NO
Jim Dempsey
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Thank you very much Jim. You solved with some minor changes my problem. Do you think are there any ways to use NODI_RETE>20 for 64 bit PC? This means to be able to calculate numbers bigger then 20! Best Regards
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Intel Fortran does not have an integer(16). Therefor to stay within the confines of the Fortran language you could change the code that uses integer(8) into real(16) and store the integers in there. Real(16) uses an exponent of 15 bits, sign is 1 bit, leaving 128-16=112 bits for precision. Or 48 bits (well 49 bits) more precision than integer(8). This would permit your factorial routine to handle up to 30!
To go above this, search the internet for multi-precision libraries, there are several available.
(here is one: https://www.researchgate.net/publication/220492964_A_Fortran-90_Based_Multiprecision_System)
As to if changing all INTEGER(kind=long_integer) into REAL(kind=long_kind) is sufficient, I will leave that up to you to test and verify.
Don't forget to change:
N=CEILING(LOG(REAL(factorial(NODI_RETE),kind=long_real))/LOG(REAL(2,kind=long_real)))
to
N=CEILING(LOG(factorial(NODI_RETE))/LOG(REAL(2,kind=long_real)))
Your interest is to avoid overflow of the precision necessary for your computation. I suggest you add some conditional compile statements to insert some asserts to detect overflows).
Jim Dempsey
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
The program code in #3 has the following error (at least; there may be others).
The first executable statement in BETAJSTAR uses the value of the INTENT(IN) variable index. This function is invoked in the line
DO WHILE(MAXVAL((/errore1,errore2/))>=errore_ammesso .OR. BetaJstar(A,BETA10,BETAJSTAR_TH,INDEX)<BETAJSTAR_TH)
but at this point the function argument variable index has no defined value.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Good catch.
Looks like a bug to me.
llario will have to correct for this. Possibly this should be M_INDEX or RETI_INDEX, I do not know.
Jim Dempsey
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page