c program lp(input=130,output=130,tape1,tape2,tape3) c c c For general details see: "An Upper Bound for the c Entropy and its Applications to Maximal Entropy Problem" by c Y. Alhassid, N. Agmon and R.D. Levine, Chem. Phys. Lett. c 53,22 (1978). c For details on the algorithm see: "An Algorithm for c Finding the Distribution of Maximal Entropy", by N. Agmon, c Y. Alhassid and R.D. Levine, J. Comp. Phys. 30, 250 (1979). c c function - given an equation Ax=b ( a(m*n),b(m),x(n) ) find: c (1) Are A rows dependent? c (2) Is there a feasible solution? c c Input data is as follows: (All are in free format) c 1) m,n,iflag,delta c m : No of constraints c n : No of states c iflag=-1: end of data c iflag=-2: only averages are read in this run c delta: accuracy parameter c 2) minimum value possible in A, maximum value possible,theta c if both 0 defaults of (-1e36,1e36) given by the program. c theta - dependecy of constraints minimal measure c 3) rows of A c 4) mode=0,1,2 or 3: c mode 0: read only averages c mode 1: read only experimental distribution c mode 2: read averages and prior distribution c mode 3: read experimental and prior distributions c if mode=3 add: exp. distribution (pexp(i),i=1,n) c prior distribution (prio(i),i=1,n) c 5) iguess - initial guess of Lagrange parameters: c iguess=0:don't read lagrange multipliers (use same values c as in previous run) c iguess=1: set initial values to zero c iguess=2: read initial values c c output - a,b, and error code if any on output file. c print out the feasible solution if any,and dependency. c print out any error occured through the run. c c calls - zx3lp from imsl library with minor changes. c c comments - a is variable dimension,up to 20 rows. c memory is given externally according to A's dimension. c in order to read tape1 mentioned above use c read(7)((a(i,j),i=1,i1),j=1,j1),(b(i),i=1,i1) c ier=70: a rows are linearly dependent. c =133: no feasible solution. c =0: no error. c implicit real*8 (a-h,o-z) common/e/g(210),f,limn,delta dimension a(16,200) common/run/nrun c c read matrix dimension nrun=0 open(unit=7,form='unformatted',status='scratch') open(unit=8,form='unformatted',status='scratch') open(unit=9,form='unformatted',status='scratch') 10 read *,me,n,iflag,delta if(iflag.eq.-1) stop 1 if(delta.le.1.e-12) delta=1.e-06 nrun=nrun+1 print *,' - - - - - - - - - - - - - - - - - - - - - - -' print *,' -- run no.',nrun,' --' if(me.gt.20) then print*,'No. of constraints greater than 20 caused abortion of j 1ob' stop 2 endif if(n.gt.200) then print*,'No. of states greater than 200 caused abortion of job' stop 2 endif if(n.eq.0) then print*,'No. of states n=0 caused abortion of job' stop 2 endif rewind 7 print*,'There are',me,' constraints and',n,' states in this run' print*,'Accuracy test:',delta call lp1(a,me+1,me,me+3,n,iflag) goto 10 end subroutine lp1(a,m,me,ia,n,iflag) implicit real*8 (a-h,o-z) c c function - read and print input matrix.call routine to find c feasible solution c variable list c m - me+1 (an additional row for 1...1) c m - me+2 (1 for inequality,1 for row of 1....1) c ia - m+2(2 rows for working storage) c n - number of unkowns(input card) c a - (ia*n) matrix of the equation to be solved(ax=b) c b - (ia) right hand size. c m1 = number of inequalities(0) c m2 - number of equalities. c dimension a(ia,n),b(23),x(20),c(21) common/mod/mode,iguess,q(23,21) common/run/nrun common/fail/nfail c c ... iflag eq -2 means that only the b vector is changed c so read the original a from tape1 if(iflag.ne.-2)goto 5 rewind 8 read(8) ((a(i,j),i=2,m),j=1,n) do 6 j=1,n 6 a(1,j)=1. b(1)=1. goto 41 c c ... read limit card for matrix variable. c if both min and max are 0- set them to default value 5 read *,bmin,bmax,theta if(bmin.ne.0..or.bmax.ne.0.) go to 15 c set default bmin=-1.e38 bmax=1.e38 theta=0.001 15 print *,'Limits for this run ',bmin,bmax c c ... read and print a,b c set row 1 for one inequality and row 2 to all ones do 25 j=1,n 25 a(1,j)=1. b(1)=1. print *,'A matrix:' do 40 i=2,m read *,(a(i,j),j=1,n) k=min0(n,5) print 34,i-1,(a(i,j),j=1,k) 34 format('row',i4,1x,5e12.4) if(k.lt.n)print 35,(a(i,j),j=k+1,n) 35 format(8x,5e12.4) c check limits do 37 j=1,n if(a(i,j).lt.bmin.or.a(i,j).gt.bmax)print 36,i,j,a(i,j) 36 format('*** a(',i3,',',i3,')=',e12.3,' so over/underflows limit') 37 continue 40 continue 41 read *,mode if(mode.eq.1.or.mode.eq.3) then read *,(a(m+1,j),j=1,n) print *,'Experimental distribution is' print 42,(a(m+1,j),j=1,n) 42 format(/(7f10.6)) do 47 j=2,m b(j)=0. do 47 i=1,n b(j)=b(j)+a(m+1,i)*a(j,i) 47 continue else c read vector b read *,(b(i),i=2,m) endif if(mode.eq.2.or.mode.eq.3) then read *,(a(m+2,j),j=1,n) print *,'Prior distribution is' print 42,(a(m+2,j),j=1,n) else do 59 i=1,n 59 a(m+2,i)=1. endif if(mode.eq.1.or.mode.eq.3) then c entropy of exp dist=pexp*ln(pexp/p(o)) szb=0. do 224 j=1,n 224 szb=szb + a(m+1,j)*dlog(a(m+1,j)/a(m+2,j)) print *,'Entropy of experimental distribution=',szb endif print *,'Given averages are:' print 35,(b(i),i=2,m) rewind 8 write(8) ((a(i,j),i=2,m),j=1,n),(b(i),i=2,m) read*,iguess if(iguess.ne.2) go to 60 read *,(x(i),i=1,me) print *,'guess: ',(x(i),i=1,me) c c dump a,b to tape1 c 60 i1=m-1 c if(nrun.eq.1.and.iguess.ne.1) iguess=1 if(nfail.eq.1.and.iguess.ne.1) iguess=1 if(iflag.eq.-2) go to 79 call gs(a,q,ia,m,n,ier,kk,theta,psi) if(ier.eq.70) go to 120 y=sqrt(float(n)) do 121 i=1,m do 121 j=1,n 121 a(i,j)=a(i,j)*y do 122 i=1,m do 122 j=1,m 122 q(i,j)=q(i,j)*y print *,'The orthogonolized constraints are:' do 124 i=2,m k=min0(n,5) print 34,i-1,(a(i,j),j=1,k) if(k.lt.n) print 35,(a(i,j),j=k+1,n) 124 continue print *,'The Orthogonalization Matrix is:' do 126 i=1,m k=min0(m,5) print 34,i-1,(q(i,j),j=1,k) if(k.lt.m) print 35,(q(i,j),j=k+1,m) 126 continue rewind 9 write(9) ((q(i,j),i=1,m),j=1,m) go to 80 79 if(ier.eq.70) go to 120 rewind 7 read(7) ((a(i,j),i=2,m),j=1,n) 80 do 90 i=1,m c(i)=0. do 90 j=1,m 90 c(i)=c(i)+q(i,j)*b(j) do 91 i=1,m 91 b(i)=c(i) print *,'The averages of the orthogonolized constraints:' print 35,(b(i),i=2,m) rewind 7 write(7) ((a(i,j),i=2,m),j=1,n),(b(i),i=2,m) write(7) (a(m+2,j),j=1,n) if(ier.eq.0) go to 81 if(ier.eq.1) go to 82 120 print*,'Constraints are linearly dependent; constraint No.',kk-1 print*,'is a linear combination of the previous ones' nfail=1 return 82 print 111,theta,kk-1,psi 111 format(/,'*** Warning: Constraints are nearly linearly dependent; c'/'angle less than ',e10.2,' first occured in line ',i2,' Minimal cangle = ',e10.2,' degrees.') go to 112 81 print *,'Minimal angle = ',psi,' degrees.' 112 if(iguess.eq.2) write(8) (x(i),i=1,me),(a(m+2,j),j=1,n) c c ... find solution and check dependency c m1=0, m2=m-m1 call lp2(a,m,n,0,m,ia,b) return end subroutine lp2(a,m,n,m1,m2,ia,b) implicit real*8 (a-h,o-z) c c function - (1) check dependency of equation. c (2) find feasible solution c c input - calling sequence parms. c output - ier=0,70,133 etc. error explanation (see IMSL c manual under ZX3LP) c an error message is printed and ier is written on tape1. c if feassible solution exists-print it. dimension a(ia,n),b(ia) dimension c(200),psol(200),dsol(26),rw(700),iw(55) common/fail/nfail do 2 j=1,n 2 c(j)=1. call zx3lp(a,ia,b,c,n,m1,m2,s,psol,dsol,rw,iw,ier) c c handle error if any. nfail=0 if(ier.eq.0) go to 30 nfail=1 if(ier.ne.70) go to 10 print 5 5 format('*** Warning: You may have a redundant constraint; you may 1be on the boundry',/,'of the simplex') nfail=0 go to 30 10 if(ier.ne.133) go to 20 print *,'***** No feasible solution exists.' return 20 print *,'***** Error No. ',ier return c feassible solution. 30 continue c check solution. rewind 7 read(7) ((a(i,j),i=2,m),j=1,n),(b(i),i=2,m) do 40 i=2,m x=0. do 35 j=1,n 35 x=x+a(i,j)*psol(j) if(abs(x-b(i)).gt.1.e-7)print 37,i 37 format('***** Solution does not match on row ',i5) 40 continue c end of check print 42,(psol(j),j=1,n) 42 format('Feasible solution is'/,(8x,6g12.5)) jok=m*n-n+1 if(ier.eq.0.or.ier.eq.70)call entpy(a,b,m-1,n,a(jok,1),a(jok+n,1)) return end subroutine entpy(a,b,n,imn,pr,p0) implicit real*8 (a-h,o-z) dimension a(n,imn),b(n),pr(imn),p0(imn),u(20),x(23),dx(20),e(21) dimension ind(200) common/e/ g(210),f,limn,delta common/j/nut common/mod/mode,iguess,q(23,21) c nut=1 error in exp at sub prob ends run. limn=imn n1=n+1 sa=1. print *,'Enter phase 3 - solution of maximal entropy problem' if(iguess.ne.2) go to 80 rewind 8 read(8) ((a(i,j),i=1,n),j=1,limn),(b(i),i=1,n) read(8) (x(i),i=1,n),(p0(i),i=1,limn) print 61 print 3,(x(i),i=1,n) do 90 i=1,n do 90 j=1,limn 90 a(i,j)=a(i,j)-b(i) r=fn(x,a,p0,n,pr) do 100 i=1,n 100 r=r-x(i)*b(i) do 81 i=1,n jj=n+1-i 81 x(jj+1)=x(jj) x(1)=r rewind 9 read(9) ((q(i,j),j=1,n1),i=1,n1) call leqt1f(q,1,n1,23,x,0,e,ie ) rewind 9 read(9) ((q(i,j),i=1,n1),j=1,n1) do 82 i=1,n 82 x(i)=x(i+1) 80 rewind 7 read(7) ((a(i,j),i=1,n),j=1,limn),(b(i),i=1,n) read(7) (p0(i),i=1,limn) if(iguess.eq.2) go to 17 if(iguess.eq.1) go to 30 print 61 print 3,(e(i),i=1,n) 17 do 56 j=1,n do 55 i=1,limn 55 a(j,i)=a(j,i)-b(j) 56 continue print*,'Subtracting the averages from every row of A yielded' print*,'the new A matrix' do 59 i=1,n k=min0(limn,5) print 34,i,(a(i,j),j=1,k) if(k.lt.limn) print 35,(a(i,j),j=k+1,limn) 59 continue nut=0 call rafi(x,dx,sa,z,a,iter,n,pr,p0,u,ier) if(nut.eq.1) return if(ier.ne.0) go to 25 16 dz=f do 18 i=1,n 18 dz=dz-x(i)*b(i) do 83 i=1,n jj=n+1-i 83 x(jj+1)=x(jj) x(1)=dz do 101 i=1,n1 e(i)=0. do 101 j=1,n1 101 e(i)=e(i)+x(j)*q(j ,i) dz=e(1) do 84 i=1,n x(i)=x(i+1) 84 e(i)=e(i+1) z=exp(dz) print 110,(x(i),i=1,n) print 7,(e(i),i=1,n) print *,'The distribution of maximal entropy is:' print 8,(pr(i),i=1,limn) open(19,file='me1',form='unformatted') write(19) (pr(i),i=1,limn) job=0 do 15 i=1,limn if(pr(i).gt.delta) goto15 job=job+1 ind(job)=i 15 continue print *,'Warning: The following calculated probablities are less t 1han',delta,' :' print 2,(ind(i),i=1,job) 2 format('p(i),i=',(20i3)) rewind 8 read(8) ((a(i,j),i=1,n),j=1,limn),(b(i),i=1,n) do 85 j=1,n do 85 i=1,limn 85 a(j,i)=a(j,i)-b(j) call grdint (e,zz,a,n,pr,p0,u) print 9,(u(i),i=1,n) if(mode.eq.3)f=-f print 11,z,dz,f nn=n*(n+1)/2 print 12 l=1 k=1 do 21 m=1,n print 3,(g(j),j=l,k) l=k+1 21 k=k+m +1 if(job.eq.1) print 10 return 25 if(ier.eq.130) print 13 if(iter.eq.0) return print 14 go to 16 30 do 31 i=1,n 31 x(i)=0. print 61 print 3 ,(x(i),i=1,n) go to 17 3 format(/,6g13.4) 110 format(/,'The Lagrange parameters of the orthogonolized constraint 1s are '/,6g13.4) 7 format(/'The Lagrange parameters are',/,(5e15.6)) 8 format(8f10.6) 9 format(/'Differences between given and calculated averages' 1,/,8f10.7) 10 format(/,'We recommend repeating this run after deleting ' 1 ,'states with too small probabilities.') 11 format(/,'Z=',e12.4,' ln(Z)=',f10.5,' S=',f10.5) 12 format(/'The correlation matrix') 13 format(/'**** Iterative improvement failed to converge.' c,/,'Hessian is too ill-conditioned. (IMSL sub LEQ2TP)') 14 format(/,'Results up to this stage are') 34 format('row',i4,1x,5g12.4) 35 format(8x,5g12.4) 61 format(/,'The initial guess for the Lagrange multipliers is') end subroutine grdint(x,z,a,n,pr,p0,u) implicit real*8 (a-h,o-z) dimension p0(1),a(n,1),pr(1),u(n),x(n) common/j/nut common/e/ g(210),f,limn call prob(x,z,a,p0,n,pr) if(nut.eq.1) return nn=n*(n+1)/2 do 7 m=1,n 7 u(m)=0. do 2 j=1,nn 2 g(j)=0. do 3 i=1,limn j=0 do 3 m=1,n v=pr(i)*a(m,i) u(m)=u(m)+v do 3 k=1,m j=j+1 3 g(j)=g(j)+v*a(k,i) j=0 do 1 m=1,n do 1 k=1,m j=j+1 1 g(j)=g(j)-u(m)*u(k) return end subroutine minima(x,dx,p,n,u,a,p0,pr) implicit real*8 (a-h,o-z) dimension p0(1),a(n,1),pr(1),x(n),dx(n),u(n) common/j/nut common/e/g(210),f,limn 10 xr=1. fr=fu(n,x,u,dx,xr ,a,p0, pr) if(nut.eq.1) return if(fr.ge.f) go to 50 20 xm=xr fm=fr xr=xr+xr fr=fu(n,x,u,dx,xr ,a,p0, pr) if(nut.eq.1) return if(fr.ge.fm) go to 40 if(abs(xr).lt.10.) go to 20 p=xr f=fr return 40 p=xm-(0.5*xm*(fr-f))/(fr+f-fm-fm) f=fu(n,x,u,dx,xr,a,p0,pr) if(nut.eq.1) return if(f.le.fm) return p=xm f=fm return 50 if(abs(f).lt.1.e-10) go to 80 if(fr-f.lt.1000.*abs(f)) go to 80 print *,'U is devided by Thousend' do 60 i=1,n 60 u(i)=.01*u(i) go to 10 80 xm=xr/2. fm=fu(n,x,u,dx,xm ,a,p0, pr) if(nut.eq.1) return if(fm.lt.f) go to 40 if(xm.lt.0.001) go to 90 xr=xm fr=fm go to 80 90 xr=-xm fr=fu(n,x,u,dx,xr ,a,p0, pr) if(nut.eq.1) return if(fr.lt.f) go to 20 p=-(0.5*xm*(fm-fr))/(fm+fr-f-f ) f=fu(n,x,u,dx,p,a,p0,pr) return end function fu(n,x,u,dx,t,a,p0,pr) implicit real*8 (a-h,o-z) dimension p0(1),a(n,1),pr(1),x(n),u(n),dx(n) do 10 i=1,n 10 dx(i)=x(i)+t*u(i) fu=fn(dx,a,p0,n,pr) return end function fn(x,a,p0, n,pr) implicit real*8 (a-h,o-z) dimension p0(1),a(n,1),pr(1),x(1) common/e/ g(210),f,limn common/j/nut call prob(x,z,a,p0,n,pr) if(nut.eq.1) return fn=dlog(z) return end subroutine prob(x,z,a,p0, n,pr) implicit real*8 (a-h,o-z) dimension a(n,1),p0(1),pr(1),x(1) common/j/nut common/e/ g(210),f,limn z=0. do 2 i=1,limn pr(i)=0. do 1 k=1,n 1 pr(i)=pr(i)+x(k)*a(k,i) if(abs(pr(i)).ge.500.) go to 5 pr(i)=exp(-pr(i))*p0(i) 2 z=z+pr(i) do 3 i=1,limn 3 pr(i)=pr(i)/z return 5 print 6 nut=1 return 6 format (/,'***** Abs value of argument greater than 500. Error c detected by exp at sub prob',/,'End of run') end subroutine gs(a,b,ma,m,n,ier,kk,psilim,psi) implicit real*8 (a-h,o-z) c c calling: c a-ma*n matrix.on exit a contains m orthonormal vectors c b-ma*m conversion matrix (e=ba where e is the orthonormal vectors) c ma-actual number of rows in a(from the dimension) c m-number of rows in a(in). c n-number of coloumns.(in). c ier-if error occurded it contains err number(out). c 70-dependency. c k-row number on which err first occured(out). c psilim-minimal angle in degrees(in) c psi-minimal angle(out) c c this program computes from a given matrix,the orthnormal c matrix,and the transformation matrix. c the method in which it works: c c let a(k) be the k-s row in a c let e(k) be the orthonormal row composed from a(k) through c the grahm-smidt process. c let s(k) be 1/norm(a(k)) c let s(k,j) be -(a(k),a(j))(inner product). c then we compose a set of matrixes the multiplication of c which will transform a into an orthonormal matrix e c e.g. m,n=3.the set of matarixes is: c c 1 0 0 1 0 0 s1 0 0 c c 0 1 0 s2*s21 s2 0 0 1 0 c c s3*s31 s3*s32 s3 0 0 1 0 0 1 c c when we multiply(right one first) we get the transformation c matrix b c in this program-in order to get efficiency we build: c c s1 0 0 c c s2*s21 s2 0 c c s3*s31*s1 s3*s32 s3 c c and then multiply row and coloumns,and get the trans matrix. dimension a(ma,n),b(23,m) data delta/1.e-20/ ier=0 psi=500. kk=0 do 5 i=1,m do 5 j=1,m 5 b(i,j)=0 c row 1. do 10 j=1,n 10 b(1,1)=b(1,1)+a(1,j)*a(1,j) c norm(a(1)) if(b(1,1).ge.delta) go to 12 ier=70 kk=1 return 12 x=sqrt(b(1,1)) b(1,1)=1./x do 15 j=1,n 15 a(1,j)=a(1,j)/x do 50 k=2,m do 20 l=1,k-1 do 20 j=1,n 20 b(k,l)=b(k,l)-a(l,j)*a(k,j) b(k,k)=1. ck=0. do 32 j=1,n x=0. do 30 i=1,k-1 30 x=x+b(k,i)*a(i,j) a(k,j)=a(k,j)+x ck=ck+x*x 32 continue c normalize a(k) x=0. do 40 j=1,n 40 x=x+a(k,j)*a(k,j) if(x.gt.delta) go to 410 ier=70 kk=k return 410 ck=atan(x/ck)*180./3.14 if(ck.lt.psi)psi=ck if(ck.ge.psilim)goto 41 if(kk.ne.0)goto 41 kk=k ier=1 41 x=sqrt(x) do 42 j=1,n 42 a(k,j)=a(k,j)/x c set b(k,j) to s(k,j)*s(k) do 44 j=1,k 44 b(k,j)=b(k,j)/x 50 continue c now we multiply rows and coloumns in b,so that it will c contain the trans matrix do 60 k=2,m do 60 j=1,k-1 x=0. do 57 i=j,k-1 57 x=x+b(k,i)*b(i,j) 60 b(k,j)=x return end subroutine rafi (x,dx,sa,z,a,iter,n,pr,p0,u,ier) implicit real*8 (a-h,o-z) dimension p0(1),a(n,1),pr(1),u(n),x(n),dx(n),wa(250) common/e/g(210),f,limn,delta common/j/nut iter=0 f=fn(x,a,p0,n,pr) if(nut.eq.1) return 44 call grdint (x,z,a,n,pr,p0,u) if(nut.eq.1) return if(iter.gt.20) go to 77 s=0. do 70 i=1,n 70 s=s+u(i)*u(i) s=s/sa if(s.lt.delta) go to 78 call leqt2p (g,1,n,20,u,id,d1,d2,wa,ier) if(ier.ne.0) go to 80 61 call minima (x,dx,p,n,u,a,p0, pr) if(nut.eq.1) return iter=iter+1 do 50 i=1,n dx(i)=p*u(i) 50 x(i)=x(i)+dx(i) go to 44 77 print *,'Over 20 iterations used' return 78 print*,iter,' iterations were needed to assure that grad2 is lt.' c,delta return 80 print*,'Error detected by IMSL sub LEQT2P at iteration No.',iter if(ier.eq.130) return print *,'***** Hessian algorithmically not positive definite' print *,'(IMSL sub LEQT2P).' if (iter.gt.3) go to 90 print *,'We try the direction of the gradient' call grdint(x,z,a,n,pr,p0,u) go to 61 90 print *,'Give up.' return end 3 20 0 1.e-7 3*0. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 0 15. 25. 35. 1 3 20 0 1.e-7 3*0. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 4 9 16 25 26 29 34 41 50 51 54 59 66 75 76 79 84 91 100 1 4 9 16 25 36 49 64 81 100 101 104 109 116 125 136 149 164 181 200 0 15. 25. 25. 1 3 20 -2 1.e-7 0 10 50 100 0 3 20 -2 1.e-7 0 1.31 2.07 2.08 1 0 0 -1 0 - - - - - - - - - - - - - - - - - - - - - - - -- run no. 1 -- There are 3 constraints and 20 states in this run Accuracy test: 1.0000000000000d-07 Limits for this run -1.0000000000000d+38 1.0000000000000d+38 A matrix: row 1 0.1000e+01 0.2000e+01 0.3000e+01 0.4000e+01 0.5000e+01 0.6000e+01 0.7000e+01 0.8000e+01 0.9000e+01 0.1000e+02 0.1200e+02 0.1300e+02 0.1400e+02 0.1500e+02 0.1600e+02 0.1800e+02 0.1900e+02 0.2000e+02 row 2 0.1100e+02 0.1200e+02 0.1300e+02 0.1400e+02 0.1500e+02 0.1600e+02 0.1700e+02 0.1800e+02 0.1900e+02 0.2000e+02 0.2200e+02 0.2300e+02 0.2400e+02 0.2500e+02 0.2600e+02 0.2800e+02 0.2900e+02 0.3000e+02 row 3 0.2100e+02 0.2200e+02 0.2300e+02 0.2400e+02 0.2500e+02 0.2600e+02 0.2700e+02 0.2800e+02 0.2900e+02 0.3000e+02 0.3200e+02 0.3300e+02 0.3400e+02 0.3500e+02 0.3600e+02 0.3800e+02 0.3900e+02 0.4000e+02 Given averages are: 0.1500e+02 0.2500e+02 0.3500e+02 Constraints are linearly dependent; constraint No. 2 is a linear combination of the previous ones - - - - - - - - - - - - - - - - - - - - - - - -- run no. 2 -- There are 3 constraints and 20 states in this run Accuracy test: 1.0000000000000d-07 Limits for this run -1.0000000000000d+38 1.0000000000000d+38 A matrix: row 1 0.1000e+01 0.2000e+01 0.3000e+01 0.4000e+01 0.5000e+01 0.6000e+01 0.7000e+01 0.8000e+01 0.9000e+01 0.1000e+02 0.1200e+02 0.1300e+02 0.1400e+02 0.1500e+02 0.1600e+02 0.1800e+02 0.1900e+02 0.2000e+02 row 2 0.1000e+01 0.4000e+01 0.9000e+01 0.1600e+02 0.2500e+02 0.2600e+02 0.2900e+02 0.3400e+02 0.4100e+02 0.5000e+02 0.5400e+02 0.5900e+02 0.6600e+02 0.7500e+02 0.7600e+02 0.8400e+02 0.9100e+02 0.1000e+03 row 3 0.1000e+01 0.4000e+01 0.9000e+01 0.1600e+02 0.2500e+02 0.3600e+02 0.4900e+02 0.6400e+02 0.8100e+02 0.1000e+03 0.1040e+03 0.1090e+03 0.1160e+03 0.1250e+03 0.1360e+03 0.1640e+03 0.1810e+03 0.2000e+03 Given averages are: 0.1500e+02 0.2500e+02 0.2500e+02 The orthogonolized constraints are: row 1 -0.1648e+01 -0.1474e+01 -0.1301e+01 -0.1127e+01 -0.9538e+00 -0.7804e+00 -0.6070e+00 -0.4336e+00 -0.2601e+00 -0.8671e-01 0.2601e+00 0.4336e+00 0.6070e+00 0.7804e+00 0.9538e+00 0.1301e+01 0.1474e+01 0.1648e+01 row 2 0.2642e+00 -0.6882e+00 -0.7160e+00 0.1807e+00 0.2002e+01 0.1251e+00 -0.8272e+00 -0.8550e+00 0.4171e-01 0.1863e+01 -0.9662e+00 -0.9940e+00 -0.9732e-01 0.1724e+01 -0.1529e+00 -0.1133e+01 -0.2363e+00 0.1585e+01 row 3 0.1258e+01 0.5362e+00 -0.1637e+00 -0.8420e+00 -0.1499e+01 -0.8939e+00 -0.2677e+00 0.3803e+00 0.1050e+01 0.1741e+01 0.2763e+00 -0.4236e+00 -0.1102e+01 -0.1758e+01 -0.1154e+01 0.1204e+00 0.7900e+00 0.1481e+01 The Orthogonalization Matrix is: row 0 0.1000e+01 0. e+00 0. e+00 0. e+00 row 1 -0.1821e+01 0.1734e+00 0. e+00 0. e+00 row 2 0.2141e+01 -0.2339e+01 0.4623e+00 0. e+00 row 3 0.2001e+01 -0.7541e+00 -0.1240e+00 0.1348e+00 The averages of the orthogonolized constraints: 0.7804e+00 -0.2139e+02 -0.9039e+01 Minimal angle = 8.3737396517982d-02 degrees. *** TERMINAL ERROR (IER = 133) FROM IMSL ROUTINE ZX3LP ***** No feasible solution exists. - - - - - - - - - - - - - - - - - - - - - - - -- run no. 3 -- There are 3 constraints and 20 states in this run Accuracy test: 1.0000000000000d-07 Given averages are: 0.1000e+02 0.5000e+02 0.1000e+03 The averages of the orthogonolized constraints: -0.8671e-01 0.1863e+01 0.1741e+01 Minimal angle = 8.3737396517982d-02 degrees. Feasible solution is 0. 0.33647e-08 0. 0. 0.34258e-08 0. 0. 0. 1.0000 0. 0. 0. 0. 0. 0. 0.64235e-08 0. 0. Enter phase 3 - solution of maximal entropy problem The initial guess for the Lagrange multipliers is 0. 0. 0. Subtracting the averages from every row of A yielded the new A matrix row 1 -1.561 -1.387 -1.214 -1.041 -0.8671 -0.6937 -0.5203 -0.3468 -0.1734 0.1735e-17 0.3468 0.5203 0.6937 0.8671 1.041 1.387 1.561 1.734 row 2 -1.599 -2.551 -2.579 -1.682 0.1390 -1.738 -2.690 -2.718 -1.821 0.2220e-15 -2.829 -2.857 -1.960 -0.1390 -2.016 -2.996 -2.099 -0.2781 row 3 -0.4834 -1.205 -1.905 -2.583 -3.240 -2.635 -2.009 -1.361 -0.6912 0.1665e-15 -1.465 -2.165 -2.843 -3.500 -2.895 -1.621 -0.9511 -0.2599 1 iterations were needed to assure that grad2 is lt. 1.0000000000000d-07 The Lagrange parameters of the orthogonolized constraints are 1.387 -29.81 -27.86 The Lagrange parameters are 0.909702e+02 -0.103257e+02 -0.375482e+01 The distribution of maximal entropy is: 0.000000 0. 0. 0. 0.000000 0. 0. 0. 1.000000 0.000000 0. 0. 0. 0. 0. 0. 0. 0.000000 Warning: The following calculated probablities are less than 1.0000000000000d-07 : p(i),i= 1 2 3 4 5 6 7 8 9 11 12 13 14 15 16 17 18 19 20 Differences between given and calculated averages 0.0000002 0.0000008 0.0000016 Z= 0.1628e-07 ln(Z)= -17.93345 S= 0.00000 The correlation matrix 0.1628e-05 0.8139e-05 0.4070e-04 0.1628e-04 0.8139e-04 0.1628e-03 - - - - - - - - - - - - - - - - - - - - - - - -- run no. 4 -- There are 3 constraints and 20 states in this run Accuracy test: 1.0000000000000d-07 Given averages are: 0.1310e+01 0.2070e+01 0.2080e+01 The averages of the orthogonolized constraints: -0.1594e+01 0.3364e-01 0.1037e+01 Minimal angle = 8.3737396517982d-02 degrees. Feasible solution is 0.72550 0.26300 0. 0. 0.11300e-01 0. 0. 0. 0.20001e-03 0. 0. 0. 0. 0. 0. 0. 0. 0. Enter phase 3 - solution of maximal entropy problem The initial guess for the Lagrange multipliers is 0. 0. 0. Subtracting the averages from every row of A yielded the new A matrix row 1 -0.5376e-01 0.1197 0.2931 0.4665 0.6399 0.8133 0.9868 1.160 1.334 1.507 1.854 2.027 2.201 2.374 2.548 2.894 3.068 3.241 row 2 0.2305 -0.7218 -0.7496 0.1471 1.968 0.9148e-01 -0.8609 -0.8887 0.8064e-02 1.829 -0.9999 -1.028 -0.1310 1.690 -0.1866 -1.167 -0.2700 1.551 row 3 0.2208 -0.5007 -1.201 -1.879 -2.536 -1.931 -1.305 -0.6567 0.1293e-01 0.7042 -0.7606 -1.461 -2.139 -2.795 -2.191 -0.9165 -0.2469 0.4443 6 iterations were needed to assure that grad2 is lt. 1.0000000000000d-07 The Lagrange parameters of the orthogonolized constraints are 7.890 0.2791 -0.3045 The Lagrange parameters are 0.945042e+00 0.166769e+00 -0.410456e-01 The distribution of maximal entropy is: 0.750344 0.200000 0.041457 0.006683 0.000838 0.000433 0.000174 0.000054 0.000002 0.000001 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 Warning: The following calculated probablities are less than 1.0000000000000d-07 : p(i),i= 13 14 15 16 17 18 19 20 Differences between given and calculated averages 0.0000423 0.0002008 0.0003489 Z= 0.4568e+00 ln(Z)= -0.78354 S= 0.71430 The correlation matrix 0.3642 1.414 5.790 1.471 6.069 6.540