MKFM6-1 Source Code

MKFM6-1 source code

program MKF_ML

c 19/7/07 some errors in error handling

c 4/2/05 update

c 5 nov 2002

c KALMAN FILTER + normal ML estimation based on the pred. err. decomp.

c filter (Harvey, 1996, p. 140ff).

c hamaker(molenaar) project 7/6/01

c aio-flynn(dolan) NWO aio project

c progress (juli 2001):

c version 2 * implement kf,

c calculate loglikelihood single subject/regime

c steady state

c version 3 * add duminf to optimize based on logl function values

c version 5 * add gradients

c & change duminf to duming -> version6 - minor clean up.

c 14/6/01

c version 7 * st errs

c add c and d (levels)

c version 9 * npsol replaces duming

c version 10 * multi-case + major overhaul i/o.

c error trapping: FUN and GRD cannot generate fatal errors

c errors there are dealt with using penalty,

c i.e. function values or

c gradients are set to large value (see penalty).

c version mkfm2 * missing data (26/02/2002) #1) harvey p. 144 ok.

c * error messages

c mkfm3 * slight change in i/o (diagonal matrix spec.)

c * se= option

c * it= option

c mkfm4 * temp version voor simulation

c mkfm5 * fixed regressors nx= & Z & exo & xo

c * xo is assumed to contain no missing data

c y[t] = S a[t] + d + e[t] + Z xo[t]

c mkfm6 * added parameter matrix G (premultiplies random innov. terms

c a[t]=H a[t-1] + c + G z[t] (Nov 5 2002)

c * inform (NPSOL) feedback

c * 13 okt 03 issues warning if G=0 and Q=1

c mkfm6.1* janne adolf 24/4/2010 fixed cov for latent process

c bug 3/feb05 bug in nsand - multi-group problem

c (calls IMSL svign, npsol)

c ------

implicit double precision (a-h,o-z)

implicit integer (i-n)

integer maxny,maxne,maxd,mpar,maxnt,maxns,maxnm

integer nclin,ncnln,nrowa,nrowj,nrowr,mtot,lenw,leniw

c parameter statement: here

parameter (maxny=30,maxnx=5,maxne=30,maxd=30)

c...... max n parameters

parameter (mpar=400)

c...... max timeseries length, max n of cases max n of models

parameter (maxnt=5000,maxns=10,maxnm=5)

c j_adolf

parameter (nparm=10)

c npsol ----

parameter (nclin=0,ncnln=0,nrowa=1,nrowj=1,nrowr=mpar)

parameter (mtot=mpar+nclin+ncnln)

c work array size for npsol

parameter (lenw=mpar*mpar+mpar*nclin+2*mpar*ncnln+

& 20*mpar+11*nclin+21*ncnln)

parameter (leniw=3*mpar+nclin+2*ncnln)

c npsol ---

double precision anps(1,1)

double precision xpar(mpar),grad(mpar),bl(mtot),bu(mtot)

double precision cnps(1),cjac(1,1),clamda(mtot),hess(mpar,mpar)

integer istate(mtot),iter

c --- additional array se (for st errors)

double precision sterr(mpar),scale

c --- npsol work arrays

double precision w(lenw)

integer iw(leniw)

logical dcheck,symm,diag

integer ih(4),im(4),is(4)

integer iyear,imonth,iday

character title*75,dfile(maxnm)*12,rfile(maxnm)*12

character mnames(nparm)*1,sp1*1,tag*1,tag2*2

character version*25,niter*5,plevel*1

c ------

c pattern matrices parameter matices

c maxnm - number of models, 10 matrices

common / jpat / imat(maxnm,nparm,maxd,maxd),

& idim(maxnm,nparm,2)

common / symm / symm(nparm),diag(maxnm,nparm)

common / xpat / xmat(maxnm,nparm,maxd,maxd),

& pst(maxnm,maxne,maxne),

& pend(maxnm,maxne,maxne),xmis(maxnm)

common / iscala / nt(maxns),ns(maxnm),model(maxnm),

& ibad(3),

& ipat(maxnm,nparm),nm,nns(maxnm,2),

& isfr(maxnm,nparm),isterr

common / xscala / dln2pi,penalty,xlogl(maxnm),zero,one

c time series: ts, filtered series: xest, states at t=0: xst

common / times / ts(maxns,maxnt,maxny),

& exo(maxns,maxnt,maxny),

& xest(maxns,maxnt,maxne),

& xst(maxns,maxne)

c posterior probablilities for mixtures

common / posts / fdens(maxns,maxnt)

data dln2pi / 1.83787706640935 /

data penalty / 100000.d0 /

data zero,one / 0.d0,1.d0 /

c lower and upper bounds on mixing proportions

c parameter names j_adolf

data mnames /'S','R','H','Q','P','d','c','Z','G','B'/

data version /'MKFv1 April 2010'/

data sp1 /' '/

c indicates whether S R H Q P d c Z G B symmetric

data symm

c S R H Q P d c

& /.false.,.true.,.false.,.true.,.true.,.false.,.false.,

c Z G B (j_adolf)

& .false.,.false.,.false. /

c check derivatives ? NPSOL check

c data dcheck / .true. /

c data plevel / '1' /

data dcheck / .false. /

data plevel / '0' /

c externals for npsol

external fun,grd,objfun,confun

c externals

external readi,readd,badrep,getse,writer

call getdat(iyear,imonth,iday)

call gettim(ih(1),im(1),is(1),i)

call setz(diag,xmat,imat,maxnm,nparm,maxd)

open(9,file='npsol.out')

open(8,file='pars.out')

write(*,991) maxnm,maxnt,maxns,maxny,maxnx,maxne,mpar

c start input script part

call readi(title,mnames,xpar,bu,bl,leniw,iw,

& dfile,rfile,npar,niter,ier)

call errcheck(ier)

c read data

call readd(dfile,mnames,idim,nm,nt,ns,nns,maxnm,

& ts,exo,xst,maxns,maxnt,maxny,maxnx,maxne,nparm,ier)

call errcheck(ier)

c end input

call nsand(nns,ns,nm,maxnt,maxnm)

call writer(version,title,dfile,symm,diag,

& npar,nm,ipat,mnames,idim,imat,xmat,pend,

& isfr,maxnm,maxnt,maxns,maxny,mpar,maxd,maxne,nparm,1)

call sumdat(ts,exo,xmis,ns,nt,nns,nm,idim,xst,dfile,

& maxnt,maxnm,maxns,maxne,nparm,ier)

call errcheck(ier)

if (npar.gt.0) then

c first call to filter + loglikelihood

call FUN(npar,xpar,f0)

c report penalty errors

call badrep(ibad,1)

if (niter(1:1).ne.'0') then

call setnp(niter,dcheck,plevel)

c optimize logl

call npsol(npar,nclin,ncnln,nrowa,nrowj,nrowr,

& anps,bl,bu,

& confun,objfun,

& inform,iter,istate,

& cnps,cjac,clamda,f0,grad,hess,xpar,

& iw,leniw,w,lenw)

endif

call badrep(ibad,2)

call grd(npar,xpar,grad)

call FUN(npar,xpar,f0)

call badrep(ibad,3)

if (isterr.eq.1) then

call gettim(ih(2),im(2),is(2),i)

c get st errs

call getse(hess,xpar,w,sterr,npar,mpar,scale,ier)

if (ier.ne.0) write(*,*) ' warning: error st.err ',ier

if (scale.gt.1.d0.and.ier.eq.0) write(*,990) scale

990 format(/,' eps in getse was set: eps * ',f8.1)

call gettim(ih(3),im(3),is(3),i)

endif

c hess error does not terminate program

call writep(mpar,npar,xpar,grad,sterr,hess,f0,bl,bu,isterr,

& inform)

call writer(version,title,dfile,symm,diag,

& npar,nm,ipat,mnames,idim,imat,xmat,pend,

& isfr,maxnm,maxnt,maxns,maxny,mpar,maxd,maxne,nparm,2)

call writed(rfile,fdens,xest,nm,idim,nns,nt,

& maxnt,maxnm,maxny,maxne,maxns,nparm,ier)

call errcheck(ier)

c npar=0

else

call fun(npar,xpar,f0)

call badrep(ibad,3)

call writep(mpar,npar,xpar,grad,sterr,hess,f0,bl,bu,isterr,

& inform)

call writed(rfile,fdens,xest,nm,idim,nns,nt,

& maxnt,maxnm,maxny,maxne,maxns,nparm,ier)

call errcheck(ier)

endif

call gettim(ih(4),im(4),is(4),i)

write(*,994)

&ih(1),im(1),is(1),ih(4),im(4),is(4),iday,imonth,iyear

994 format(//,'start:',1x,i2,1x,i2,1x,i2,' end:',1x,i2,1x,i2,1x,i2,

&' date:',1x,i2,1x,i2,1x,i4)

write(*,995) ih(2),im(2),is(2),ih(3),im(3),is(3)

995 format(/,'start getse:',1x,i2,1x,i2,1x,i2,

&' end getse:',1x,i2,1x,i2,1x,i2)

991 format('max nm=',i2,' nt=',i4,' ns=',i4,' ny=',i2,' nx=',i2,

& ' ne=',i2,' npar=',i3,//)

stop

end

subroutine setnp(niter,dcheck,plevel)

character niter*5,plevel*1

logical dcheck

call npoptn('der 3')

if (dcheck) then

call npoptn('verify yes')

else

call npoptn('verify no')

endif

call npoptn('print level '//plevel)

call npoptn('major iter '//niter)

call npoptn('minor iter 200 ')

return

end

subroutine setz(diag,xmat,imat,maxnm,nparm,maxd)

implicit double precision (a-h,o-z)

implicit integer (i-n)

double precision xmat(maxnm,nparm,maxd,maxd)

integer imat(maxnm,nparm,maxd,maxd)

logical diag(maxnm,nparm)

do 1 i=1,maxnm

do 1 j=1,nparm

diag(i,j)=.false.

do 1 k=1,maxd

do 1 l=1,maxd

xmat(i,j,k,l)=0.d0

imat(i,j,k,l)=0

1 continue

return

end

subroutine sumdat(ts,exo,xmis,ns,nt,nns,nm,idim,xst,dfile,

& maxnt,maxnm,maxns,maxne,nparm,ier)

c 'S','R','H','Q','P','d','c','Z',,'G' 'B'/

implicit double precision (a-h,o-z)

implicit integer (i-n)

parameter (maxny=30,maxnx=5)

double precision ts(maxns,maxnt,maxny),xmis(maxnm)

double precision exo(maxns,maxnt,maxnx)

double precision xst(maxns,maxne)

integer ns(maxnm),nns(maxnm,2),nt(maxns)

integer idim(maxnm,nparm,2)

double precision wmis(maxny),wmean(maxny),wvar(maxny)

double precision wmin(maxny),wmax(maxny)

character dfile(maxnm)*12

write(*,991)

do 1 inm=1,nm

c write(*,992)

ny=idim(inm,1,1)

ne=idim(inm,1,2)

nx=idim(inm,8,2)

write(*,992) inm,nm,ny,nx,ne,ns(inm),nns(inm,1),nns(inm,2)

do 2 ins=nns(inm,1),nns(inm,2)

do 4 i=1,ny

wmin(i)=100000000.d0

wmax(i)=-100000000.d0

wmean(i)=0.d0

wvar(i)=0.d0

4 wmis(i)=0.d0

misny=0

do 3 itime=1,nt(ins)

itmp=0

do 5 i=1,ny

if (ts(ins,itime,i).eq.xmis(inm)) then

itmp=itmp+1

wmis(i)=wmis(i)+1.

else

if (ts(ins,itime,i).gt.wmax(i)) wmax(i)=ts(ins,itime,i)

if (ts(ins,itime,i).lt.wmin(i)) wmin(i)=ts(ins,itime,i)

wmean(i)=wmean(i)+ts(ins,itime,i)

wvar(i)=wvar(i)+ts(ins,itime,i)**2

endif

5 continue

if (itmp.eq.ny) misny=misny+1

if (nx.gt.0) then

do 13 inx=1,nx

if (exo(ins,itime,inx).eq.xmis(inm)) then

ier=5230

return

endif

13 continue

endif

3 continue

do 6 i=1,ny

wmean(i)=wmean(i)/(dble(nt(ins))-wmis(i))

wvar(i)=(wvar(i)/(dble(nt(ins))-wmis(i)))-wmean(i)**2

wmis(i)=wmis(i)/dble(nt(ins))

6 continue

write(*,997) ins,nt(ins),misny,dfile(inm)

write(*,990) (xst(ins,i),i=1,ne)

nr=6

iend=0

12 ist=iend+1

iend=ist+nr

if (iend.gt.ny) iend=ny

write(*,999) (i,i=ist,iend)

write(*,996) (wmis(i),i=ist,iend)

write(*,995) (wmean(i),i=ist,iend)

write(*,994) (wvar(i),i=ist,iend)

write(*,993) (dsqrt(wvar(i)),i=ist,iend)

write(*,983) (wmin(i),i=ist,iend)

write(*,984) (wmax(i),i=ist,iend)

if (nx.gt.0) write(*,981) nx

if (iend.lt.ny) goto 12

write(*,*)

2 continue

1 continue

981 format(/,'Number of fixed regressors ',i2,/)

998 format(i2,' of ',i2,2x,i2,2x,i2,2x,i4,4x,i4,2x,i4,/)

991 format(/,'DATA SUMMARY ',/)

992 format('MODEL ',i2,' of ',i2,1x,'NY=',i2,1x,'NX=',i2,1x,

&'NE=',i2,1x,'Ncases=',i4,1x,'START=',i4,1x,'END=',i4,/)

997 format('CASE ',i4,3x,' T=',i4,' N of T missing=',i4

&,' datafile ',a12,/)

990 format(' State_0 ',8(1x,f6.2),/)

999 format(' var ',10(1x,i8))

996 format(' %miss',10(1x,f8.4))

995 format(' mean ',10(1x,f8.2))

994 format(' var ',10(1x,f8.2))

993 format(' std ',10(1x,f8.2))

983 format(' min ',10(1x,f8.2))

984 format(' max ',10(1x,f8.2))

return

end

subroutine writed(rfile,fdens,xest,nm,idim,nns,nt,

& maxnt,maxnm,maxny,maxne,maxns,nparm,ier)

implicit double precision (a-h,o-z)

implicit integer (i-n)

integer idim(maxnm,nparm,2),nns(maxnm,2),nt(maxns)

double precision xest(maxns,maxnt,maxne)

double precision fdens(maxns,maxnt)

character rfile(12)*12,lasto*12

ier=0

ifo=0

c 6300

do 1 inm=1,nm

if (rfile(inm).eq.'no') goto 1

if (inm.eq.1.or.ifo.eq.0) then

ifo=1

open(10,file=rfile(inm))

lasto=rfile(inm)

else if (ifo.eq.1.and.(rfile(inm).ne.lasto)) then

close(10,status='keep')

open(10,file=rfile(inm))

endif

do 90 ins=nns(inm,1),nns(inm,2)

do 92 itime=1,nt(ins)

write(10,822)

& inm,ins,itime,(xest(ins,itime,j),j=1,idim(inm,1,2))

92 continue

90 continue

1 continue

822 format(1x,i2,1x,i2,1x,i4,6(1x,f8.3))

821 format(5(1x,f8.4))

return

end

subroutine writep(mpar,npar,xpar,grad,sterr,hess,f0,bl,bu,

& isterr,inform)

double precision xpar(mpar),grad(mpar),sterr(mpar)

double precision f0,hess(mpar,mpar),bl(mpar),bu(mpar)

integer i

character onb*3

write(*,*)

if (npar.eq.0) then

write(*,*) ' all parameters fixed '

write(*,998) -f0,2.d0*f0

return

endif

write(*,*) 'ML parameter estimates '

iwarnl=0

iwarnu=0

do 22 i=1,npar

write(8,*) xpar(i)

onb=' '

if (xpar(i).eq.bl(i)) then

iwarnl=iwarnl+1

onb='bl!'

endif

if (xpar(i).eq.bu(i)) then

iwarnu=iwarnu+1

onb='bu!'

endif

if (isterr.eq.1)

& write(*,918) i,onb,xpar(i),grad(i),sterr(i),xpar(i)/sterr(i)

if (isterr.eq.0)

& write(*,919) i,onb,xpar(i),grad(i)

22 continue

if (iwarnl.gt.0)

& write(*,*) ' * Parameter(s) on lower bound ',iwarnl

if (iwarnu.gt.0)

& write(*,*) ' * Parameter(s) on upper bound ',iwarnu

write(*,998) -f0,2.d0*f0,inform

if (inform.ne.0) then

write(*,*) ' See NPSOL MANUAL for INFORM = ',inform

write(*,*)

endif

918 format

&(' nr ',i3,1x,a3,1x,f15.5,' g ',f12.6,' se ',f12.4,' t ',f10.2)

919 format

&(' nr ',i3,1x,a3,1x,f15.5,' g ',f12.6)

998 format(/,' Logl ',f15.3,' -2xLogL ',f15.3,

&' Inform(NPSOL) ',i2,/)

c write(*,*) ' inv. info matrix (correlations)'

c do 23 i=1,npar

c23 write(*,997) i,(hess(i,j),j=1,i-1)

c997 format(1x,' par ',i3,8(1x,f6.3))

return

end

subroutine writer(version,title,dfile,symm,diag,

& npar,nm,ipat,mnames,idim,imat,xmat,pend,

& isfr,maxnm,maxnt,maxns,maxny,mpar,maxd,maxne,nparm,indi)

implicit double precision (a-h,o-z)

implicit integer (i-n)

double precision xmat(maxnm,nparm,maxd,maxd)

integer imat(maxnm,nparm,maxd,maxd)

double precision pend(maxnm,maxne,maxne)

integer idim(maxnm,nparm,2),ipat(maxnm,nparm)

integer isfr(maxnm,nparm)

character dfile(maxnm)*12,mnames(nparm)*1,version*25

character title*75

logical symm(nparm),diag(maxnm,nparm)

c before npsol

if (indi.eq.1) write(*,999) version

if (indi.eq.1.and.npar.gt.0) then

write(*,'(a75)') title

do 2 inm=1,nm

write(*,998) inm,nm

do 6 im=1,nparm

if (ipat(inm,im).eq.0.or.isfr(inm,im).eq.0) goto 6

if (diag(inm,im)) then

write(*,991) mnames(im)

irw=idim(inm,im,2)

write(*,996)(imat(inm,im,ir,ir),ir=1,irw)

else

write(*,997) mnames(im)

do 3 icl=1,idim(inm,im,1)

irw=idim(inm,im,2)

if (symm(im)) irw=icl

3 write(*,996) (imat(inm,im,icl,ir),ir=1,irw)

endif

6 continue

2 continue

endif

if (indi.eq.2.and.npar.gt.0) then

write(*,'(a75)') title

do 12 inm=1,nm

write(*,998) inm,nm

do 16 im=1,nparm

if (ipat(inm,im).eq.0) goto 16

if (diag(inm,im)) then

write(*,992) mnames(im)

irw=idim(inm,im,2)

write(*,993)(xmat(inm,im,ir,ir),ir=1,irw)

else

write(*,994) mnames(im)

do 14 icl=1,idim(inm,im,1)

irw=idim(inm,im,2)

if (symm(im)) irw=icl

14 write(*,993) (xmat(inm,im,icl,ir),ir=1,irw)

endif

16 continue

write(*,*)

write(*,981)

981 format('P(t|t) error cov ')

do 17 i=1,idim(inm,5,1)

17 write(*,993) (pend(inm,i,j),j=1,i)

write(*,*)

12 continue

endif

999 format(/,26x,19('='),/,26x,a25,/,26x,19('='),/)

998 format(/,'Model ',i2,' of ',i2)

997 format(/,a1,1x,'fr parameters (nonzero) ')

991 format(/,a1,1x,'fr parameters (nonzero) - diagonal ')

996 format(10(1x,i3))

995 format(/,'Mix proportions free ',i3)

994 format(/,a1,1x,'parameters ')

992 format(/,a1,1x,'parameters - diagonal ')

993 format(8(1x,f8.3))

return

end

subroutine badrep(ibad,indi)

integer ibad(3),i,indi,itmp

itmp=ibad(1)+ibad(2)+ibad(3)

if (itmp.eq.0) return

write(*,*)

if (indi.eq.1) write(*,*) 'N of penalized error - on entry '

if (indi.eq.2) write(*,*) 'N of penalized error - during npsol'

if (indi.eq.3) write(*,*) 'N of penalized error - on exit '

if (ibad(1).gt.0) write(*,999) ibad(1)

if (ibad(2).gt.0) write(*,998) ibad(2)

if (ibad(3).gt.0) write(*,997) ibad(3)

999 format(' --- P matrix not posdef in KF ',i4,' times ')

998 format(' --- Density less/equal zero ',i4,' times ')

997 format(' --- Mix prob out of bouns ',i4,' times ')

write(*,*)

do 1 i=1,3

1 ibad(i)=0

return

end

subroutine errcheck(ier)

integer ier

character error*50

error=' '

if (ier.eq.0) return

call geterror(ier,error)

write(*,999) ier,error

999 format(//,' fatal...error number ',i4,/,

&' message: ',a50)

777 read(*,'(a50)',end=888) error

if (error(1:1).eq.'!'.or.error(1:10).eq.' ') goto 777

write(*,*)

write(*,*) ' next line in input file following this error '

write(*,*)

write(*,'(a50)') error

888 stop

return

end

subroutine geterror(ier,error)

character error*50

integer ier

error=' '

error=' error not documented ? '

if (ier.eq.5000)

&error='error reading "title" '

if (ier.eq.5005)

&error='eof reading "title" '

if (ier.eq.5010)

&error='error reading nm= line '

if (ier.eq.5015)

&error='eof reading nm= line '

if (ier.eq.5020)

&error='cannot find "nm=" '

if (ier.eq.5021)

&error='cannot find "se=" '

if (ier.eq.5025)

&error='error reading "nm=" '

if (ier.eq.5026)

&error='error reading "se=" '

if (ier.eq.5027)

&error='should be se=yes or se=no '

if (ier.eq.5201)

&error='nm < 0 of nm > maxnm '

if (ier.eq.5030)

&error='error reading mo line '

if (ier.eq.5035)

&error='eof reading mo line '

if (ier.eq.5040)

&error='cannot find "mo=" '

if (ier.eq.5045)

&error='error reading "mo=" '

if (ier.eq.5050)

&error='cannot find "ny=" '

if (ier.eq.5051)

&error='cannot find "nx=" '

if (ier.eq.5055)

&error='error reading "ny=" '

if (ier.eq.5205)

&error='ny le 0 or ny > maxny '

if (ier.eq.5060)

&error='cannot find "ne=" '

if (ier.eq.5065)

&error='error reading "ne=" '

if (ier.eq.5210)

&error='ne le 0 or ne > maxne '

if (ier.eq.5075)

&error='error reading data line '

if (ier.eq.5080)

&error='eof reading data line '

if (ier.eq.5085)

&error='cannot find "df=" '

if (ier.eq.5090)

&error='error reading "df=" '

if (ier.eq.5092)

&error='cannot find "rf=" '

if (ier.eq.5094)

&error='error reading "rf=" '

if (ier.eq.5095)

&error='cannot find "ns=" '

if (ier.eq.5100)

&error='error reading "ns=" '

if (ier.eq.5105)

&error='cannot find "mi=" '

if (ier.eq.5107)

&error='error reading "mi=" '

if (ier.eq.5110)

&error='error reading model specification line ("S=",etc.)'

if (ier.eq.5115)

&error='eof reading model specification line ("S=",etc.) '

if (ier.eq.5120)

&error='cannot find "S=" '

if (ier.eq.5125)

&error='error reading "S=" '

if (ier.eq.5130)

&error='cannot find "R=" '

if (ier.eq.5135)

&error='error reading "R=" '

if (ier.eq.5140)

&error='cannot find "H=" '

if (ier.eq.5145)

&error='error reading "H=" '

if (ier.eq.5150)

&error='cannot find "Q=" '

if (ier.eq.5155)

&error='error reading "Q=" '

if (ier.eq.5160)

&error='cannot find "d=" '

if (ier.eq.5165)

&error='error reading "d=" '

if (ier.eq.5170)

&error='cannot find "c=" '

if (ier.eq.5175)

&error='error reading "c=" '

if (ier.eq.5180)

&error='cannot find "P=" '

if (ier.eq.5181)

&error='cannot find "Z=" '

if (ier.eq.5182)

&error='error reading "Z=" '

c j_adolf

if (ier.eq.5280)

&error='cannot find "B=" '

if (ier.eq.5281)

&error='error reading "B=" '

if (ier.eq.5183)

&error='cannot find "G=" '

if (ier.eq.5184)

&error='error reading "G=" '

if (ier.eq.5185)

&error='error reading "P=" '

if (ier.eq.5186)

&error='error Q=1 and G=0 (or vice versa), so GQG` = zero '

if (ier.eq.5191)

&error='S ne 0 or S ne 1 '

if (ier.eq.5192)

&error='R ne 0 or R ne 1 '

if (ier.eq.5193)

&error='H ne 0 or H ne 1 '

if (ier.eq.5194)

&error='Q ne 0 or Q ne 1 '

if (ier.eq.5195)

&error='P ne 0 or P ne 1 '

if (ier.eq.5196)

&error='d ne 0 or d ne 1 '

if (ier.eq.5197)

&error='c ne 0 or c ne 1 '

if (ier.eq.5198)

&error='Z ne 0 or Z ne 1 '

if (ier.eq.5199)

&error='G ne 0 or G ne 1 '

if (ier.eq.5200)

&error='B ne 0 or B ne 1 '

if (ier.eq.6000)

&error='eof reading "S fr", etc. etc. '

if (ier.eq.6020)

&error='keyword "fr" or "fi" expected, not found '

if (ier.eq.6025)

&error='matrix name incorrect (not S,R,Q,H,P,d,c,Z,B,G) '

if (ier.eq.6009)

&error='"S" fr or "S fi" encountered while "S=0" '

if (ier.eq.6010)

&error='"R" fr or "R fi" encountered while "R=0" '