Program g13befe ! G13BEF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: g13bef, nag_wp, x04abf, x04caf ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: iset = 1, nin = 5, nout = 6 ! .. Local Scalars .. Real (Kind=nag_wp) :: d, s Integer :: dp, i, ifail, imwa, inc, isttf, itc, & iwa, kef, kfc, kpriv, kzef, kzsp, & ldcm, ldxxy, mx, nadv, ncd, nce, & ncf, ncg, ndf, ndv, nis, nit, npara, & nser, nsttf, nxxy, qp, qx, smx ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: cm(:,:), para(:), res(:), sd(:), & sttf(:), wa(:), xxy(:,:) Real (Kind=nag_wp) :: zsp(4) Integer :: mr(7) Integer, Allocatable :: mt(:,:), mwa(:) ! .. Intrinsic Procedures .. Intrinsic :: max, sum ! .. Executable Statements .. Write (nout,*) 'G13BEF Example Program Results' Write (nout,*) ! Skip heading in data file Read (nin,*) ! Read in problem size Read (nin,*) kzef, kfc, nxxy, nser, kef, nit, kzsp, kpriv If (kzsp/=0) Then Read (nin,*) zsp End If ! Number of input series nis = nser - 1 ! Set the advisory channel to NOUT for monitoring information If (kpriv/=0) Then nadv = nout Call x04abf(iset,nadv) End If Allocate (mt(4,nser)) ! Read in orders Read (nin,*) mr(1:7) ! Read in transfer function Do i = 1, nis Read (nin,*) mt(1:4,i) End Do ! Calculate NPARA and various other quantities required ! for calculate array sizes npara = 0 ncg = 0 qx = 0 smx = 0 ncf = nser inc = 1 Do i = 1, nis npara = npara + mt(2,i) + mt(3,i) If (mt(4,i)>1) Then ncg = ncg + sum(mt(1:3,i)) If (mt(4,i)==3) Then mx = max(mt(1,i)+mt(2,i),mt(3,i)) qx = max(qx,mx) smx = smx + mx End If Else If (mt(4,i)==1 .And. kef==3) Then If (mt(3,i)>0) Then ncf = ncf + 1 End If inc = inc + 1 End If End Do npara = npara + mr(1) + mr(3) + mr(4) + mr(6) + nser ! Calculate size of arrays isttf = mr(4)*mr(7) + mr(2) + mr(5)*mr(7) + mr(3) + & max(mr(1),mr(6)*mr(7)) + ncg ldxxy = nxxy ldcm = npara qp = mr(3) + mr(6)*mr(7) dp = mr(2) + mr(5)*mr(7) If (mr(3)>0 .And. kef>1) Then inc = inc + 1 End If If (kfc>0 .And. kef==3) Then inc = inc + 1 End If qx = qp ncd = npara + kfc + smx If (mr(1)>0) Then ncf = ncf + inc End If If (mr(3)>0) Then ncf = ncf + inc End If If (mr(4)>0) Then ncf = ncf + inc End If If (mr(6)>0) Then ncf = ncf + inc End If If (qx>0) Then ncf = ncf + 1 End If If (kfc>0) Then ncf = ncf + 1 End If ncd = ncd + qx nce = nxxy + dp + 6*qx iwa = 2*ncd**2 + nce*(ncf+4) iwa = 2*iwa imwa = 16*nser + 7*ncd + 3*npara + 3*kfc + 27 Allocate (xxy(ldxxy,nser),para(npara),sd(npara),cm(ldcm,npara), & res(nxxy),sttf(isttf),wa(iwa),mwa(imwa)) ! Read in rest of data Read (nin,*) para(1:npara) Read (nin,*)(xxy(i,1:nser),i=1,nxxy) ifail = -1 Call g13bef(mr,nser,mt,para,npara,kfc,nxxy,xxy,ldxxy,kef,nit,kzsp,zsp, & itc,sd,cm,ldcm,s,d,ndf,kzef,res,sttf,isttf,nsttf,wa,iwa,mwa,imwa, & kpriv,ifail) If (ifail/=0) Then If (ifail/=8 .And. ifail/=9) Then Go To 100 End If End If ! Display results Write (nout,99999) 'The number of iterations carried out is', itc Write (nout,*) Write (nout,*) & 'Final values of the parameters and their standard deviations' Write (nout,*) Write (nout,*) ' I PARA(I) SD' Write (nout,*) Write (nout,99998)(i,para(i),sd(i),i=1,npara) Write (nout,*) Flush (nout) ifail = 0 Call x04caf('General',' ',npara,npara,cm,ldcm, & 'The correlation matrix is',ifail) Write (nout,*) Write (nout,*) 'The residuals and the z and n values are' Write (nout,*) Write (nout,*) ' I RES(I) z(t) n(t)' Write (nout,*) ndv = nxxy - mr(2) - mr(5)*mr(7) Do i = 1, nxxy If (i<=ndv) Then Write (nout,99997) i, res(i), xxy(i,1:nser) Else Write (nout,99996) i, xxy(i,1:nser) End If End Do If (mr(2)/=0 .Or. mr(5)/=0) Then Write (nout,*) Write (nout,*) & '** Note that the residuals relate to differenced values **' End If Write (nout,*) Write (nout,99995) 'The state set consists of', nsttf, ' values' Write (nout,*) Write (nout,99994) sttf(1:nsttf) Write (nout,*) Write (nout,99999) 'The number of degrees of freedom is', ndf 100 Continue 99999 Format (1X,A,I4) 99998 Format (1X,I4,2F20.6) 99997 Format (1X,I4,3F15.3) 99996 Format (1X,I4,F30.3,F15.3) 99995 Format (1X,A,I4,A) 99994 Format (1X,6F10.4) End Program g13befe