c----------------- HP version, August 1995 --------- program nnreg implicit double precision(a-h,o-z) parameter(kmax=10,jdmax=16,nxmax=3000) parameter(jdmaxo=16,nxmaxo=3000) parameter(npmax=1+kmax*(jdmax+2)) dimension ydat(nxmax),xmat(nxmax,jdmax),xxmat(nxmax,jdmax), * oydat(nxmaxo),oxmat(nxmaxo,jdmaxo),oxxmat(nxmaxo,jdmaxo), * theta(npmax),ctheta(npmax),pre2(nxmax),opre2(nxmaxo), * tsav(200,npmax),rsav(200), * yhat(nxmax),oyhat(nxmaxo) dimension h0(npmax,npmax), xsd(jdmax), xmean(jdmax) dimension oxsd(jdmaxo), oxmean(jdmaxo) character*35 fname,outname,ofname,ooutname data glow, ghigh, scale / -1.26, 1.26, 0.5 /, * maxstep1, maxstep2 / 50, 250 /, * ibrent, fltol / 0, 0.0 / common /size/ k,jd,nx,m,jforce common /xdata/ ydat,xmat common /oxdata/ oydat,oxmat external objfun cgcv=1.0 outsamp=0 outend=0 c----------- Data statement initializes many parameters c----------- GCV cost=2, ibrent and fltol=0 so never used c----------- Set initial Jacobian to identity matrix open (1,file='test1.text') call setidmat(h0,npmax) c----------- Read in parameters and minimization routine options open(11,file='nnrego.par',status='old') read(11,10) fname read(11,10) ofname read(11,10) outname read(11,10) ooutname 10 format(a35) read(11,*) nx,nxc,nxo,nxco if( (nx.gt.nxmax).or.(nxc.gt.jdmax)) then write(*,*) 'dataset too large!' stop endif if( (nxo.gt.nxmaxo).or.(nxco.gt.jdmaxo)) then write(*,*) 'dataset too large!' stop endif read(11,*) ngrid, ntries, npol, iprout, igreed read(11,*) ngrido c----------- Initialize random number generator read(11,*) iseed iseed=abs(mod(iseed,125000)) rn=rng(iseed) c----------- if ngrid is negative set switch to read in start values if( ngrid.lt.0) then istart=1 npol=1 ntries=1 ngrid=1 else istart=0 endif c if( (npol.gt.200).or.(npol.gt.ngrid)) then write(*,*)'too many values to polish' stop endif if( (igreed.gt.0) .and.(iprout.gt.0)) then write(*,*) 'Greedy algorithm can only output the best fit' write(*,*) ' Not all the polsihed results' stop endif read(11,*) ftol1, ftol2, itmax1, itmax2 iprint=itmax1+itmax2 read(11,*) k1,k2 close(11) if( ( 1+ (nxc+2)*k2).gt.npmax) then write(*,*) 'too many parameters in the model for array size' stop endif if( ( k1.gt.kmax).or.( k2.gt.kmax)) then write(*,*) ' too many hidden units' stop endif m= nx jd=nxc mo= nxo jdo=nxco c----------- read data series into xdat; set nx=total series length call getxy(fname,xmat,ydat,nxmax,nx,nxc) c----------- read data series into oxdat; set nxo=total series length call getxy(ofname,oxmat,oydat,nxmaxo,nxo,nxco) c1025 c----------- Set up the output file open(12,file=outname,status='unknown',access='append') c----- open(13,file=ydat,status='unknown',access='append') write(12,*) '# first row of data X, Y:' write(12,*) '#',(xmat(1,k), k=1,nxc), ydat(1) write(12,*) '# first row of data oX, oY:' write(12,*) '#',(oxmat(1,k), k=1,nxco), oydat(1) write(12,*) '#',' last row of data oX, oY:' write(12,*) '#',(oxmat(nxo,k), k=1,nxco), oydat( nxo) write(12,*) '#',' last row of data X, Y:' write(12,*) '#',(xmat(nx,k), k=1,nxc), ydat( nx) write(12,*) '# number of observations:', nx write(12,*) '# number of columns of X', nxc if( istart.eq.0) then write(12,*) '# number of points for grid search', ngrid write(12,*)'#', glow, ghigh, ngrid else write(12,*)'# starting values read in ' write(12,*)'# a grid search has been omitted !' endif if( igreed.gt.0) then write(12,*) '# greedy algorithm used to fit net' write(12,*) '# in steps of ', igreed endif write(12,*) '# d k RMS BIC GCV par iter' close(12) 5000 continue c------------ Standardize values to predict: SD=1, mean=0 do 500 ic=1,nxc call sdev(xmat(1,ic),1,nx,tempsd, tempm) xmean(ic)=tempm xsd(ic)=tempsd do 501 jj=1,nx xmat(jj,ic)= (xmat(jj,ic)-tempm)/tempsd 501 continue 500 continue c1025a c if(outsamp.eq.1) then do 510 ic=1,nxco call sdev(oxmat(1,ic),1,nxo,otempsd, otempm) oxmean(ic)=otempm oxsd(ic)=otempsd do 511 jj=1,nxo oxmat(jj,ic)= (oxmat(jj,ic)-otempm)/otempsd 511 continue 510 continue c endif c call sdev(ydat,1,nx,ysd, ymean) do 520 jj=1,nx ydat(jj) = (ydat(jj)-ymean)/ysd 520 continue call sdev(oydat,1,nxo,oysd, oymean) do 530 jj=1,nxo oydat(jj) = (oydat(jj)-oymean)/oysd 530 continue csept20 open(12,file=outname,status='unknown',access='append') write(12,*) 'standardized ydat sept　20　1', (ydat(j) ,j=1,nx) close(12) c----------- beginning of output for S function c----------- open file of start values if(( istart.eq.1).and.(outsamp.ne.1)) then open(2,file='nnreg.str',status='unknown') endif 5999 continue c----------- START LOOP on # units do 1004 ik=k1,k2 k=ik c----- reset number of hidden units to one for the greedy algorithm if( (igreed.gt.0).and.( ik.gt.k1)) then k=igreed endif c----------- compute number of parameters npar=1+k*(jd+2+2*jforce) if (2*npar.gt.m) goto 1004 if (npar.gt.npmax) goto 1004 c----------- if istart=1 then read in start values cprob c 5999 continue if(( istart.eq.1).and.(outsamp.ne.1)) then c if( istart.eq.1) then read(11,*) ( tsav(1,jpar), jpar=1,npar) endif do 15 i=1,200 rsav(i)=10000. 15 continue c----------- START LOOP on repeated fits c----------- if start=0 then do random search for starting values if( istart.eq.0) then c----------- to match lenns glow= -1.26 ghigh =1.26 ngrid = 256 step = (ghigh- glow)/(ngrid-1) do 99 irep=1,ngrid c----------- initialize parameter choices eps=scale*10**( step*irep + glow ) call tinit(theta,npar,eps,ntries) c----------- do a high-tolerance fit call bfgsfm(objfun,npar,npmax,theta,h0,ftol1,fltol, * itmax1,maxstep1,fnew,iter,iprint,inform,ibrent) xm=float(m) bic=fnew*sqrt( 1.0+npar*log(xm)/(xm-npar) ) c fnew=rms error; xm = sample size ( dim(X)[1] ) ; c This is actually the sqrt of what is usually called BIC. c old bic bic=1.419+dlog(fnew)+0.5*float(npar)*dlog(xm)/xm imax=idamax(npol,rsav,1) rmax=rsav(imax) if (fnew.lt.rmax) then rsav(imax)=fnew do 85 jpar=1,npar tsav(imax,jpar)=theta(jpar) 85 continue 9000 format( 5e16.8) endif 99 continue endif c----------- END LOOP on repeated fits if(outsamp.ne.1) then do 999 i=1,npol do 997 jpar=1,npar theta(jpar)=tsav(i,jpar) 997 continue c----- entering subroutine theta are starting values from grinding call bfgsfm(objfun,npar,npmax,theta,h0,ftol2,fltol, * itmax2,maxstep2,fnew,iter,iprint,inform,ibrent) c----------- convert the paramteres to canonical form to avoid c----------- identification problems call canpar(theta,k, nxc, ctheta) c----------- save new parameter estimates. do 998 jpar=1,npar tsav(i,jpar)= ctheta(jpar) 998 continue xm=float(m) bic=fnew*sqrt( 1.0+npar*log(xm)/(xm-npar) ) c fnew=rms error; xm = sample size ( dim(X)[1] ) ; c This is actually the sqrt of what is usually called BIC. gcv=(fnew/(1.-(cgcv*float(npar)/xm)))**2 if( (fnew.lt.ftest).or.(i.eq.1)) then ibest=i ftest=fnew endif c----------- Write out results if( iprout.gt.0) then write(*,8900) nxc, k, fnew write(*,9000) (xmean(jj),jj=1,nxc) write(*,9000) (xsd(jj), jj=1,nxc) write(*,9000) ymean, ysd write(*,9000) (ctheta(jj), jj=1,npar) endif open(12,file=outname,status='old',access='append') write(12,98) jd,k,fnew,bic,gcv,npar,iter, * inform close(12) 97 format(1x,2i3,2x,f7.3,2(2x,f9.6),2x,i5,i3) 98 format(1x,2i3,2x,f9.6,2(2x,f7.4),i3,i7,i3,2(2x,f7.3)) 999 continue endif c loop if(outsamp.ne.1) then do 1002 jpar=1,npar ctheta(jpar)=tsav(ibest,jpar) 1002 continue endif c----- write out the fit with smallest rms out of polished results if( iprout.eq.0) then csept20 8900 format( i8, i8, e15.8) csept20 open(12,file=outname,status='unknown',access='append') write(12,*) 'nxc, k, ftest, sept　20　1' write(12,8900) nxc, k, ftest if(outsamp.eq.1) then write(12,*) 'nxco, k, ftest, sept　20　1' write(12,8900) nxco, k, ftest endif do 9001 ii=1,nx pre1=0 do 8999 kk=1, k ggg= 0 do 8998 jj=1,nxc c xmat(jj,ic)= dexp(xmat(jj,ic)-tempm)/tempsd xxmat(kk,jj)= ctheta(2*k+1+jj+nxc*(kk-1))*xmat(ii,jj) ggg= xxmat(kk,jj) + ggg 8998 continue uum = ggg+ ctheta(k+kk+1) hide = dexp(uum)/(1+dexp(uum)) bhide= hide*ctheta(kk+1) pre1= pre1+bhide 8999 continue pre2(ii)= pre1+ctheta(1) 9001 continue write(12,*) 'prediction ' write(12,*) (pre2(ii), ii=1,nx) write(12,*) 'ctheta(jj), jj=1,npar sept　20　1' write(12,9000) (ctheta(jj), jj=1,npar) c xmean(ic)=tempm c xsd(ic)=tempsd if(outsamp.eq.1) then write(12,*) 'oxmat(jj,ic), jj=1,nxo oct 25　' write(12,*) ( (oxmat(jj,ic), jj=1,nxo) ,ic=1,nxco) write(12,*) 'oymean, oysd oct　25　1' write(12,9000) oymean, oysd else write(12,*) 'xmat(jj,ic), jj=1,nx sept　20　1' write(12,*) ( (xmat(jj,ic), jj=1,nx) ,ic=1,nxc) write(12,*) 'ymean sept　20　1' write(12,9000) ymean, ysd endif c write(*,9000) (xmean(jj),jj=1,nxc) c write(*,9000) (xsd(jj), jj=1,nxc) endif close(12) c----- calculate residual vector and substitute for y if this c----- is the greedy algorithm csept28 if( igreed.ge.0) then call netev( npar,ctheta,k,xmat, nxmax, nx,nxc,yhat) call netev( npar,ctheta,k,oxmat, nxmaxo, nxo,nxco,oyhat) if(outsamp.eq.1) then do 9100 jj=1,nxo write(1,*) 'k',k,'oydat' ,oydat(jj),'oyhat' ,oyhat(jj) 9100 continue else do 9110 jj=1,nx write(1,*) 'k',k,'ydat' ,ydat(jj),'yhat' ,yhat(jj) 9110 continue endif c if( igreed.gt.0) then do 9115 jj=1,nx c----- find residuals on raw scale ydat(jj)= (ydat(jj)- yhat(jj))*ysd 9115 continue c----- reset the scale for next fit call sdev(ydat,1,nx,ysd, ymean) do 9200 jj=1,nx ydat(jj) = (ydat(jj)-ymean)/ysd 9200 continue endif endif c goto 5999 1004 continue outend=outend+1 istart=1 outsamp=1 if(outend.gt.1) goto 1014 goto 5999 1014 continue c----------- END LOOP on #units if( istart.eq.1) close(2) open(12,file=outname,status='old',access='append') write(12,*) '@end ' close(12) close(1) c close(13) c stop end c ------------------------------------------------------------------------- c Single-layer net, "quadratic" squasher c Inputs: X (n-d) x d data matrix (d lags, n data points) c gama K x d matrix (K=# of hidden units) c beta K+1 x 1 column vector c mu K x 1 row vector c Output: gk= beta * G(X*gama' .+ mu') c ------------------------------------------------------------------------- subroutine net(xout,xin,beta,amu,gama,nx,jd,k) implicit real*8(a-h,o-z) parameter(kmax=8,jdmax=16,nxmax=3000, * hf=0.5,one=1.d0,two=2.d0) dimension xout(nxmax),xin(nxmax,jdmax),uin(nxmax,kmax) dimension gama(kmax,jdmax),beta(kmax+1),amu(kmax) do 10 it=1,nx xout(it)=beta(1) 10 continue do 50 jun=1,k u0=amu(jun) do 20 it=1,nx uin(it,jun)=u0 20 continue do 40 lag=1,jd gjunlag=gama(jun,lag) do 30 it=1,nx uin(it,jun)=uin(it,jun)+gjunlag*xin(it,lag) 30 continue 40 continue 50 continue do 100 jun=1,k do 90 it=1,nx uinij=uin(it,jun) auin=abs(uinij) xout(it)=xout(it)+beta(1+jun)*uinij * *(one+hf*auin)/(two+auin+hf*auin**2) 90 continue 100 continue return end C============================================================= subroutine objfun(mode,npar,theta,f,g,nstate) c c npar = # of variables (input integer) c theta = evaluation point (input vector of length nvar) c f = function value (output scalar) c g = gradient value (output vector of length nvar) c mode == 0 then compute f, leave g unchanged c == 1 then compute g, leave f unchanged c == 2 then compute both f and g c nstate = not used; included for compatibility with NPSOL implicit real*8(a-h,o-z) parameter(kmax=8,jdmax=16,nxmax=3000) parameter(zero=0.d0,one=1.d0,hf=0.5d0,two=2.d0) common /size/ k,jd,nx,m,jforce common /xdata/ xt,xlag dimension gama(kmax,jdmax) dimension xt(nxmax),xout(nxmax),xlag(nxmax,jdmax) dimension theta(npar),g(npar) m0=k+1 j0=2*k+1 jdf=jd+2*jforce c ------ extract gamma matrix from theta vector do 30 iunit=1,k loc0=j0+(iunit-1)*jdf do 20 lag=1,jdf gama(iunit,lag)=theta(loc0+lag) 20 continue 30 continue c ---------- compute predicted values (xout) call net(xout,xlag,theta(1),theta(m0+1),gama,nx,jdf,k) c----------- compute objective function f (mean square prediction error) sse=zero if (mode.eq.0) then do 50 i=1,m sse=sse+(xout(i)-xt(i))**2 50 continue f=sqrt(sse/dfloat(m)) else c ---------- compute objective function & its gradient vector do 55 j=1,npar g(j)=zero 55 continue do 90 i=1,m fi=xout(i)-xt(i) sse=sse+fi**2 g(1)=g(1)+fi do 80 j=1,k uin=theta(m0+j) do 60 lag=1,jdf uin=uin+gama(j,lag)*xlag(i,lag) 60 continue auin=abs(uin) dij=(two+auin +hf*(auin**2)) g(j+1)=g(j+1)+fi*uin*(one + hf*auin)/dij gpij=two*(one+auin)/(dij**2) g(m0+j)=g(m0+j)+fi*theta(1+j)*gpij loc0=j0+(j-1)*jdf do 70 lag=1,jdf g(loc0+lag)=g(loc0+lag) & +fi*theta(1+j)*xlag(i,lag)*gpij 70 continue 80 continue 90 continue fac=one/sqrt(sse*dfloat(m)) do 120 j=1,npar g(j)=fac*g(j) 120 continue if (mode.eq.2) f=sqrt(sse/dfloat(m)) endif return end subroutine bfgsfm(fcn,nvar,np,p,h0,ftol,fltol,maxiter, * maxstep,fnew,iter,iprint,inform,ibrent) implicit real*8(a-h,o-z) character*1 trans parameter(nvmax=100,eps=1.d-16,one=1.0,zero=0.0,trans='N') dimension p(np),pnew(nvmax),h0(np,np),g(nvmax), * gnew(nvmax),h(nvmax,nvmax),xi(nvmax),work(nvmax), * dx(nvmax),shy(nvmax),shys(nvmax,nvmax), * pc(nvmax),gc(nvmax) external fcn sftol=sqrt(ftol) cftol=ftol**0.35 ymin=-0.001 inform=0 jset=1 call fcn(2,nvar,p,fp,g,0) do 10 i=1,nvar do 5 j=1,nvar h(i,j)=h0(i,j) 5 continue 10 continue c---------------- set xi = -h*g call dgemv(trans,nvar,nvar,-1.d0,h,nvmax,g,1,zero,xi,1) nfail=0 jhct=0 do 99 iter=1,maxiter 12 call bhhhstp(fcn,nvar,p,fp,xi,g,maxstep,s,fnew,iters,info) c---------------- check for failed step; reset hessian & search direction if (fnew.gt.fp) then if (jset.eq.1) then inform=1 fnew=fp goto 999 endif nfail=nfail+1 jset=1 do 15 i=1,nvar do 14 j=1,nvar h(i,j)=h0(i,j) 14 continue 15 continue call dgemv(trans,nvar,nvar,-1.d0,h,nvmax,g,1,zero,xi,1) goto 12 else c---------------- set pnew=p+s*xi and compute new gradient call dcopy(nvar,p,1,pnew,1) call daxpy(nvar,s,xi,1,pnew,1) call fcn(1,nvar,pnew,fjunk,gnew,0) if (ibrent.eq.1) then c-------------- try to bracket minimum ax=0. bx=s fb=fnew step=s do 16 iterb=1,5 cx=bx+step call dcopy(nvar,p,1,pc,1) call daxpy(nvar,cx,xi,1,pc,1) call fcn(0,nvar,pc,fc,gc,0) if (fc.gt.fb) then goto 17 else ax=bx fb=fc bx=cx step=2.*step endif 16 continue 17 continue c------------- If bracketing succeeded, use Brent's method to find minimum if (fc.gt.fb) then fnew=ffmin(fcn,nvar,p,xi,ax,bx,cx,fb,fltol,xmin) call dcopy(nvar,p,1,pnew,1) call daxpy(nvar,xmin,xi,1,pnew,1) jbr=jbr+1 c--------------If bracketing failed, use lowest value found else s=cx fnew=fc call dcopy(nvar,pc,1,pnew,1) endif endif c------------------ print out results every iprint iteration itermod=mod(iter,iprint) if (itermod.eq.0) then imax=idamax(nvar,gnew,1) gnorm=ddot(nvar,gnew,1,gnew,1) write(*,21)iter,fnew,(fp-fnew)/abs(fp),gnorm, * nfail,jhct,jbr 21 format(1x,i8,2x,f12.8,2x,e13.6,2x,f12.8,2x,3i4) nfail=0 jhct=0 jbr=0 endif c----------------- Check if convergence criteria satisfied imax=idamax(nvar,gnew,1) gmax=abs(gnew(imax)) if(gmax.lt.eps) then goto 999 else if ( (fp-fnew).lt.ftol*(1.+fnew)) then if(gmax.lt.cftol*(1.+abs(fnew))) then imax=idamax(nvar,pnew,1) pmax=abs(pnew(imax)) call dcopy(nvar,pnew,1,work,1) call daxpy(nvar,-1.d0,p,1,work,1) imax=idamax(nvar,work,1) dpmax=abs(work(imax)) if (dpmax.lt.sftol*(1.+pmax)) then goto 999 endif endif endif endif c--------------------- If not, prepare for the next iteration c set dx=pnew-p call dcopy(nvar,pnew,1,dx,1) call daxpy(nvar,-1.d0,p,1,dx,1) c set shy=dx-h*gnew+h*g call dcopy(nvar,dx,1,shy,1) call dgemv(trans,nvar,nvar,-1.d0,h,nvmax,gnew,1,1.d0,shy,1) call dgemv(trans,nvar,nvar,1.d0,h,nvmax,g,1,1.d0,shy,1) c set shys=(shy*dx')+(shy*dx')' do 30 i=1,nvar do 25 j=1,nvar shys(i,j)=shy(i)*dx(j) + shy(j)*dx(i) 25 continue 30 continue c set yshy=gnew'*shy-g'*shy y1=ddot(nvar,gnew,1,shy,1) y2=ddot(nvar,g,1,shy,1) yshy=y1-y2 c set ys=gnew'*dx-g'*dx y1=ddot(nvar,gnew,1,dx,1) y2=ddot(nvar,g,1,dx,1) ys=y1-y2 c update h if ys>0, otherwise reset h if (ys.gt.0.) then jset=0 call dger(nvar,nvar,-yshy/ys,dx,1,dx,1,shys,nvmax) do 40 i=1,nvar do 35 j=1,nvar h(i,j)=h(i,j)+shys(i,j)/ys 35 continue 40 continue else jset=1 do 50 i=1,nvar do 45 j=1,nvar h(i,j)=h0(i,j) 45 continue 50 continue jhct=jhct+1 endif c update f,g,p, and search direction call dcopy(nvar,gnew,1,g,1) call dcopy(nvar,pnew,1,p,1) fp=fnew call dgemv(trans,nvar,nvar,-1.d0,h,nvmax,g,1,0.d0,xi,1) c------------------- If next step "points up", reset Hessian y1=ddot(nvar,xi,1,g,1) y2=ddot(nvar,xi,1,xi,1) if (y1.ge.ymin*y2) then jset=1 do 60 i=1,nvar do 55 j=1,nvar h(i,j)=h0(i,j) 55 continue 60 continue call dgemv(trans,nvar,nvar,-1.d0,h,nvmax,g,1,0.d0,xi,1) jhct=jhct+1 endif endif 99 continue 999 return end subroutine bhhhstp(fcn,nvar,x0,f0,d,g,itermax,s,fs, * iters,info) implicit real*8(a-h,o-z) parameter(delta=0.25,pwr=0.618,one=1.0,nvmax=100) dimension x0(nvar),d(nvar),g(nvar),work(nvmax) c--------------------- Initializations info=1 iter=0 iup=0 idown=0 factor=2. dg=ddot(nvar,d,1,g,1) downfact=dg*delta upfact=dg*(1.-delta) alam=1.0 if (itermax.eq.0) itermax=25 do 10 iter=1,itermax c----------------- Set g= alam*d + x0, vof=fcn(g) call dcopy(nvar,x0,1,work,1) call daxpy(nvar,alam,d,1,work,1) call fcn(0,nvar,work,vof,g,1) c------------------ vof above cone, hence reduce step size if ((vof-f0).gt.(downfact*alam)) then idown=1 if (iup.eq.1) then factor=factor**pwr iup=0 endif alam=alam/factor c------------------ vof below cone, hence increase step size else if ((vof-f0).lt.(upfact*alam)) then iup=1 if (idown.eq.1) then factor=factor**pwr idown=0 endif alam=alam*factor c------------------ vof within cone, hence stop iterations else info=0 goto 20 endif 10 continue 20 continue c------------------- compute returned values iters=iter s=alam fs=vof return end C function ffmin(ax,bx,cx,f,tol) function ffmin(fcn,nvar,x0,dx,AX,BX,CX,fb,TOL,XMIN) implicit real*8(a-h,o-z) parameter(nvmax=100,itmax=100) dimension x0(nvar),dx(nvar),work(nvmax),g(nvmax) c c an approximation x to the point where f attains a minimum on c the interval (ax,bx) is determined. c c c input.. c c ax left endpoint of initial interval c bx right endpoint of initial interval c f function subprogram which evaluates f(x) for any x c in the interval (ax,bx) c tol desired length of the interval of uncertainty of the final c result ( .ge. 0.0d0) c c c output.. c c fmin abcissa approximating the point where f attains a minimum c c c the method used is a combination of golden section search and c successive parabolic interpolation. convergence is never much slower c than that for a fibonacci search. if f has a continuous second c derivative which is positive at the minimum (which is not at ax or c bx), then convergence is superlinear, and usually of the order of c about 1.324.... c the function f is never evaluated at two points closer together c than eps*abs(fmin) + (tol/3), where eps is approximately the square c root of the relative machine precision. if f is a unimodal c function and the computed values of f are always unimodal when c separated by at least eps*abs(x) + (tol/3), then fmin approximates c the abcissa of the global minimum of f on the interval ax,bx with c an error less than 3*eps*abs(fmin) + tol. if f is not unimodal, c then fmin may approximate a local, but perhaps non-global, minimum to c the same accuracy. c this function subprogram is a slightly modified version of the c algol 60 procedure localmin given in richard brent, algorithms for c minimization without derivatives, prentice - hall, inc. (1973). c C C Modified by S Ellner to use relative error tolerance, tol*(1.+dabs(x)), C instead of absolute error tolerance; and to assume a machine C epsilon of 10d-16. July 1992 C C This version is adapted for calling by bfgsi, and uses the C same conventions concerning the objective function. c c c c is the squared inverse of the golden ratio c c = 0.5d0*(3. - dsqrt(5.0d0)) c c eps is approximately the square root of the relative machine c precision. c eps=1.d-8 c c initialization c a = min(ax,cx) b = max(ax,cx) v = bx w = v x = v e = 0.0d0 fx = fb fv = fx fw = fx c c main loop starts here c do 86 iter=1,itmax xm = 0.5d0*(a + b) tol1 = eps + tol*(1.+dabs(x)) tol2 = 2.0d0*tol1 c c check stopping criterion c if (dabs(x - xm) .le. (tol2 - 0.5d0*(b - a))) go to 90 c c is golden-section necessary c if (dabs(e) .le. tol1) go to 40 c c fit parabola c r = (x - w)*(fx - fv) q = (x - v)*(fx - fw) p = (x - v)*q - (x - w)*r q = 2.0d00*(q - r) if (q .gt. 0.0d0) p = -p q = dabs(q) r = e e = d c c is parabola acceptable? c 30 if (dabs(p) .ge. dabs(0.5d0*q*r)) go to 40 if (p .le. q*(a - x)) go to 40 if (p .ge. q*(b - x)) go to 40 c c a parabolic interpolation step c d = p/q u = x + d c c f must not be evaluated too close to ax or bx c if ((u - a) .lt. tol2) d = dsign(tol1, xm - x) if ((b - u) .lt. tol2) d = dsign(tol1, xm - x) go to 50 c c a golden-section step c 40 if (x .ge. xm) then e = a - x else e = b - x endif d = c*e c c f must not be evaluated too close to x c 50 if (dabs(d) .ge. tol1) then u = x + d else u = x + dsign(tol1, d) endif call dcopy(nvar,x0,1,work,1) call daxpy(nvar,u,dx,1,work,1) call fcn(0,nvar,work,fu,g,1) C fu = f(u) c c update a, b, v, w, and x c if (fu .gt. fx) go to 60 if (u .ge. x) then a = x else b=x endif v = w fv = fw w = x fw = fx x = u fx = fu go to 86 60 if (u .lt. x) then a = u else b = u endif if (fu .le. fw) go to 70 if (w .eq. x) go to 70 if (fu .le. fv) go to 80 if (v .eq. x) go to 80 if (v .eq. w) go to 80 go to 86 70 v = w fv = fw w = u fw = fu go to 86 80 v = u fv = fu go to 86 c c end of main loop c 86 CONTINUE 90 xmin = x ffmin=fx return end double precision function dasum(n,dx,incx) c c takes the sum of the absolute values. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dtemp integer i,incx,m,mp1,n,nincx c dasum = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dtemp = dtemp + dabs(dx(i)) 10 continue dasum = dtemp return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,6) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dabs(dx(i)) 30 continue if( n .lt. 6 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,6 dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2)) * + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5)) 50 continue 60 dasum = dtemp return end SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, LDA, M, N CHARACTER*1 TRANS * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * DGEMV performs one of the matrix-vector operations * * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, * * where alpha and beta are scalars, x and y are vectors and A is an * m by n matrix. * * Parameters * ========== * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' y := alpha*A*x + beta*y. * * TRANS = 'T' or 't' y := alpha*A'*x + beta*y. * * TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * X - DOUBLE PRECISION array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. * Before entry, the incremented array X must contain the * vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * BETA - DOUBLE PRECISION. * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - DOUBLE PRECISION array of DIMENSION at least * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. * Before entry with BETA non-zero, the incremented array Y * must contain the vector y. On exit, Y is overwritten by the * updated vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 1 ELSE IF( M.LT.0 )THEN INFO = 2 ELSE IF( N.LT.0 )THEN INFO = 3 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 ELSE IF( INCY.EQ.0 )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * Set LENX and LENY, the lengths of the vectors x and y, and set * up the start points in X and Y. * IF( LSAME( TRANS, 'N' ) )THEN LENX = N LENY = M ELSE LENX = M LENY = N END IF IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( LENX - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( LENY - 1 )*INCY END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * * First form y := beta*y. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )THEN IF( BETA.EQ.ZERO )THEN DO 10, I = 1, LENY Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, LENY Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE IY = KY IF( BETA.EQ.ZERO )THEN DO 30, I = 1, LENY Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE ELSE DO 40, I = 1, LENY Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE END IF END IF END IF IF( ALPHA.EQ.ZERO ) $ RETURN IF( LSAME( TRANS, 'N' ) )THEN * * Form y := alpha*A*x + y. * JX = KX IF( INCY.EQ.1 )THEN DO 60, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) DO 50, I = 1, M Y( I ) = Y( I ) + TEMP*A( I, J ) 50 CONTINUE END IF JX = JX + INCX 60 CONTINUE ELSE DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IY = KY DO 70, I = 1, M Y( IY ) = Y( IY ) + TEMP*A( I, J ) IY = IY + INCY 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF ELSE * * Form y := alpha*A'*x + y. * JY = KY IF( INCX.EQ.1 )THEN DO 100, J = 1, N TEMP = ZERO DO 90, I = 1, M TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 100 CONTINUE ELSE DO 120, J = 1, N TEMP = ZERO IX = KX DO 110, I = 1, M TEMP = TEMP + A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 120 CONTINUE END IF END IF * RETURN * * End of DGEMV . * END SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA INTEGER INCX, INCY, LDA, M, N * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * DGER performs the rank 1 operation * * A := alpha*x*y' + A, * * where alpha is a scalar, x is an m element vector, y is an n element * vector and A is an m by n matrix. * * Parameters * ========== * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - DOUBLE PRECISION array of dimension at least * ( 1 + ( m - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the m * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * Y - DOUBLE PRECISION array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. * Unchanged on exit. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. On exit, A is * overwritten by the updated matrix. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JY, KX * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( M.LT.0 )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( INCY.EQ.0 )THEN INFO = 7 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGER ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) $ RETURN * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( INCY.GT.0 )THEN JY = 1 ELSE JY = 1 - ( N - 1 )*INCY END IF IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) DO 10, I = 1, M A( I, J ) = A( I, J ) + X( I )*TEMP 10 CONTINUE END IF JY = JY + INCY 20 CONTINUE ELSE IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( M - 1 )*INCX END IF DO 40, J = 1, N IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) IX = KX DO 30, I = 1, M A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE END IF JY = JY + INCY 40 CONTINUE END IF * RETURN * * End of DGER . * END subroutine daxpy(n,da,dx,incx,dy,incy) c c constant times a vector plus a vector. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1),da integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (da .eq. 0.0d0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dy(i) + da*dx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 dy(i) = dy(i) + da*dx(i) dy(i + 1) = dy(i + 1) + da*dx(i + 1) dy(i + 2) = dy(i + 2) + da*dx(i + 2) dy(i + 3) = dy(i + 3) + da*dx(i + 3) 50 continue return end subroutine dcopy(n,dx,incx,dy,incy) c c copies a vector, x, to a vector, y. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1) integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,7) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dx(i) 30 continue if( n .lt. 7 ) return 40 mp1 = m + 1 do 50 i = mp1,n,7 dy(i) = dx(i) dy(i + 1) = dx(i + 1) dy(i + 2) = dx(i + 2) dy(i + 3) = dx(i + 3) dy(i + 4) = dx(i + 4) dy(i + 5) = dx(i + 5) dy(i + 6) = dx(i + 6) 50 continue return end double precision function ddot(n,dx,incx,dy,incy) c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1),dtemp integer i,incx,incy,ix,iy,m,mp1,n c ddot = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy 10 continue ddot = dtemp return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dx(i)*dy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 50 continue 60 ddot = dtemp return end double precision function dnrm2 ( n, dx, incx) integer next double precision dx(1), cutlo, cuthi, hitest, sum, xmax,zero,one data zero, one /0.0d0, 1.0d0/ c c euclidean norm of the n-vector stored in dx() with storage c increment incx . c if n .le. 0 return with result = 0. c if n .ge. 1 then incx must be .ge. 1 c c c.l.lawson, 1978 jan 08 c c four phase method using two built-in constants that are c hopefully applicable to all machines. c cutlo = maximum of dsqrt(u/eps) over all known machines. c cuthi = minimum of dsqrt(v) over all known machines. c where c eps = smallest no. such that eps + 1. .gt. 1. c u = smallest positive no. (underflow limit) c v = largest no. (overflow limit) c c brief outline of algorithm.. c c phase 1 scans zero components. c move to phase 2 when a component is nonzero and .le. cutlo c move to phase 3 when a component is .gt. cutlo c move to phase 4 when a component is .ge. cuthi/m c where m = n for x() real and m = 2*n for complex. c c values for cutlo and cuthi.. c from the environmental parameters listed in the imsl converter c document the limiting values are as follows.. c cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are c univac and dec at 2**(-103) c thus cutlo = 2**(-51) = 4.44089e-16 c cuthi, s.p. v = 2**127 for univac, honeywell, and dec. c thus cuthi = 2**(63.5) = 1.30438e19 c cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. c thus cutlo = 2**(-33.5) = 8.23181d-11 c cuthi, d.p. same as s.p. cuthi = 1.30438d19 c data cutlo, cuthi / 8.232d-11, 1.304d19 / c data cutlo, cuthi / 4.441e-16, 1.304e19 / data cutlo, cuthi / 8.232d-11, 1.304d19 / c if(n .gt. 0) go to 10 dnrm2 = zero go to 300 c 10 assign 30 to next sum = zero nn = n * incx c begin main loop i = 1 20 go to next,(30, 50, 70, 110) 30 if( dabs(dx(i)) .gt. cutlo) go to 85 assign 50 to next xmax = zero c c phase 1. sum is zero c 50 if( dx(i) .eq. zero) go to 200 if( dabs(dx(i)) .gt. cutlo) go to 85 c c prepare for phase 2. assign 70 to next go to 105 c c prepare for phase 4. c 100 i = j assign 110 to next sum = (sum / dx(i)) / dx(i) 105 xmax = dabs(dx(i)) go to 115 c c phase 2. sum is small. c scale to avoid destructive underflow. c 70 if( dabs(dx(i)) .gt. cutlo ) go to 75 c c common code for phases 2 and 4. c in phase 4 sum is large. scale to avoid overflow. c 110 if( dabs(dx(i)) .le. xmax ) go to 115 sum = one + sum * (xmax / dx(i))**2 xmax = dabs(dx(i)) go to 200 c 115 sum = sum + (dx(i)/xmax)**2 go to 200 c c c prepare for phase 3. c 75 sum = (sum * xmax) * xmax c c c for real or d.p. set hitest = cuthi/n c for complex set hitest = cuthi/(2*n) c 85 hitest = cuthi/float( n ) c c phase 3. sum is mid-range. no scaling. c do 95 j =i,nn,incx if(dabs(dx(j)) .ge. hitest) go to 100 95 sum = sum + dx(j)**2 dnrm2 = dsqrt( sum ) go to 300 c 200 continue i = i + incx if ( i .le. nn ) go to 20 c c end of main loop. c c compute square root and adjust for scaling. c dnrm2 = xmax * dsqrt(sum) 300 continue return end subroutine drot (n,dx,incx,dy,incy,c,s) c c applies a plane rotation. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1),dtemp,c,s integer i,incx,incy,ix,iy,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = c*dx(ix) + s*dy(iy) dy(iy) = c*dy(iy) - s*dx(ix) dx(ix) = dtemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c 20 do 30 i = 1,n dtemp = c*dx(i) + s*dy(i) dy(i) = c*dy(i) - s*dx(i) dx(i) = dtemp 30 continue return end subroutine drotg(da,db,c,s) c c construct givens plane rotation. c jack dongarra, linpack, 3/11/78. c modified 9/27/86. c double precision da,db,c,s,roe,scale,r,z c roe = db if( dabs(da) .gt. dabs(db) ) roe = da scale = dabs(da) + dabs(db) if( scale .ne. 0.0d0 ) go to 10 c = 1.0d0 s = 0.0d0 r = 0.0d0 go to 20 10 r = scale*dsqrt((da/scale)**2 + (db/scale)**2) r = dsign(1.0d0,roe)*r c = da/r s = db/r 20 z = s if( dabs(c) .gt. 0.0d0 .and. dabs(c) .le. s ) z = 1.0d0/c da = r db = z return end subroutine dscal(n,da,dx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c double precision da,dx(1) integer i,incx,m,mp1,n,nincx c if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dx(i) = da*dx(i) 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end integer function idamax(n,dx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dmax integer i,incx,ix,n c idamax = 0 if( n .lt. 1 ) return idamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 dmax = dabs(dx(1)) ix = ix + incx do 10 i = 2,n if(dabs(dx(ix)).le.dmax) go to 5 idamax = i dmax = dabs(dx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 dmax = dabs(dx(1)) do 30 i = 2,n if(dabs(dx(i)).le.dmax) go to 30 idamax = i dmax = dabs(dx(i)) 30 continue return end SUBROUTINE XERBLA ( SRNAME, INFO ) * .. Scalar Arguments .. INTEGER INFO CHARACTER*6 SRNAME * .. * * Purpose * ======= * * XERBLA is an error handler for the Level 2 BLAS routines. * * It is called by the Level 2 BLAS routines if an input parameter is * invalid. * * Installers should consider modifying the STOP statement in order to * call system-specific exception-handling facilities. * * Parameters * ========== * * SRNAME - CHARACTER*6. * On entry, SRNAME specifies the name of the routine which * called XERBLA. * * INFO - INTEGER. * On entry, INFO specifies the position of the invalid * parameter in the parameter-list of the calling routine. * * * Auxiliary routine for Level 2 Blas. * * Written on 20-July-1986. * * .. Executable Statements .. * WRITE (*,99999) SRNAME, INFO * STOP * 99999 FORMAT ( ' ** On entry to ', A6, ' parameter number ', I2, $ ' had an illegal value' ) * * End of XERBLA. * END LOGICAL FUNCTION LSAME ( CA, CB ) * .. Scalar Arguments .. CHARACTER*1 CA, CB * .. * * Purpose * ======= * * LSAME tests if CA is the same letter as CB regardless of case. * CB is assumed to be an upper case letter. LSAME returns .TRUE. if * CA is either the same as CB or the equivalent lower case letter. * * N.B. This version of the routine is only correct for ASCII code. * Installers must modify the routine for other character-codes. * * For EBCDIC systems the constant IOFF must be changed to -64. * For CDC systems using 6-12 bit representations, the system- * specific code in comments must be activated. * * Parameters * ========== * * CA - CHARACTER*1 * CB - CHARACTER*1 * On entry, CA and CB specify characters to be compared. * Unchanged on exit. * * * Auxiliary routine for Level 2 Blas. * * -- Written on 20-July-1986 * Richard Hanson, Sandia National Labs. * Jeremy Du Croz, Nag Central Office. * * .. Parameters .. INTEGER IOFF PARAMETER ( IOFF=32 ) * .. Intrinsic Functions .. INTRINSIC ICHAR * .. Executable Statements .. * * Test if the characters are equal * LSAME = CA .EQ. CB * * Now test for equivalence * IF ( .NOT.LSAME ) THEN LSAME = ICHAR(CA) - IOFF .EQ. ICHAR(CB) END IF * RETURN * * The following comments contain code for CDC systems using 6-12 bit * representations. * * .. Parameters .. * INTEGER ICIRFX * PARAMETER ( ICIRFX=62 ) * .. Scalar Arguments .. * CHARACTER*1 CB * .. Array Arguments .. * CHARACTER*1 CA(*) * .. Local Scalars .. * INTEGER IVAL * .. Intrinsic Functions .. * INTRINSIC ICHAR, CHAR * .. Executable Statements .. * * See if the first character in string CA equals string CB. * * LSAME = CA(1) .EQ. CB .AND. CA(1) .NE. CHAR(ICIRFX) * * IF (LSAME) RETURN * * The characters are not identical. Now check them for equivalence. * Look for the 'escape' character, circumflex, followed by the * letter. * * IVAL = ICHAR(CA(2)) * IF (IVAL.GE.ICHAR('A') .AND. IVAL.LE.ICHAR('Z')) THEN * LSAME = CA(1) .EQ. CHAR(ICIRFX) .AND. CA(2) .EQ. CB * END IF * * RETURN * * End of LSAME. * END subroutine dswap (n,dx,incx,dy,incy) c c interchanges two vectors. c uses unrolled loops for increments equal one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1),dtemp integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dx(ix) dx(ix) = dy(iy) dy(iy) = dtemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,3) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp 30 continue if( n .lt. 3 ) return 40 mp1 = m + 1 do 50 i = mp1,n,3 dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp dtemp = dx(i + 1) dx(i + 1) = dy(i + 1) dy(i + 1) = dtemp dtemp = dx(i + 2) dx(i + 2) = dy(i + 2) dy(i + 2) = dtemp 50 continue return end SUBROUTINE DQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,JOB) INTEGER LDX,N,P,JOB INTEGER JPVT(1) DOUBLE PRECISION X(LDX,1),QRAUX(1),WORK(1) C C DQRDC USES HOUSEHOLDER TRANSFORMATIONS TO COMPUTE THE QR C FACTORIZATION OF AN N BY P MATRIX X. COLUMN PIVOTING C BASED ON THE 2-NORMS OF THE REDUCED COLUMNS MAY BE C PERFORMED AT THE USERS OPTION. C C ON ENTRY C C X DOUBLE PRECISION(LDX,P), WHERE LDX .GE. N. C X CONTAINS THE MATRIX WHOSE DECOMPOSITION IS TO BE C COMPUTED. C C LDX INTEGER. C LDX IS THE LEADING DIMENSION OF THE ARRAY X. C C N INTEGER. C N IS THE NUMBER OF ROWS OF THE MATRIX X. C C P INTEGER. C P IS THE NUMBER OF COLUMNS OF THE MATRIX X. C C JPVT INTEGER(P). C JPVT CONTAINS INTEGERS THAT CONTROL THE SELECTION C OF THE PIVOT COLUMNS. THE K-TH COLUMN X(K) OF X C IS PLACED IN ONE OF THREE CLASSES ACCORDING TO THE C VALUE OF JPVT(K). C C IF JPVT(K) .GT. 0, THEN X(K) IS AN INITIAL C COLUMN. C C IF JPVT(K) .EQ. 0, THEN X(K) IS A FREE COLUMN. C C IF JPVT(K) .LT. 0, THEN X(K) IS A FINAL COLUMN. C C BEFORE THE DECOMPOSITION IS COMPUTED, INITIAL COLUMNS C ARE MOVED TO THE BEGINNING OF THE ARRAY X AND FINAL C COLUMNS TO THE END. BOTH INITIAL AND FINAL COLUMNS C ARE FROZEN IN PLACE DURING THE COMPUTATION AND ONLY C FREE COLUMNS ARE MOVED. AT THE K-TH STAGE OF THE C REDUCTION, IF X(K) IS OCCUPIED BY A FREE COLUMN C IT IS INTERCHANGED WITH THE FREE COLUMN OF LARGEST C REDUCED NORM. JPVT IS NOT REFERENCED IF C JOB .EQ. 0. C C WORK DOUBLE PRECISION(P). C WORK IS A WORK ARRAY. WORK IS NOT REFERENCED IF C JOB .EQ. 0. C C JOB INTEGER. C JOB IS AN INTEGER THAT INITIATES COLUMN PIVOTING. C IF JOB .EQ. 0, NO PIVOTING IS DONE. C IF JOB .NE. 0, PIVOTING IS DONE. C C ON RETURN C C X X CONTAINS IN ITS UPPER TRIANGLE THE UPPER C TRIANGULAR MATRIX R OF THE QR FACTORIZATION. C BELOW ITS DIAGONAL X CONTAINS INFORMATION FROM C WHICH THE ORTHOGONAL PART OF THE DECOMPOSITION C CAN BE RECOVERED. NOTE THAT IF PIVOTING HAS C BEEN REQUESTED, THE DECOMPOSITION IS NOT THAT C OF THE ORIGINAL MATRIX X BUT THAT OF X C WITH ITS COLUMNS PERMUTED AS DESCRIBED BY JPVT. C C QRAUX DOUBLE PRECISION(P). C QRAUX CONTAINS FURTHER INFORMATION REQUIRED TO RECOVER C THE ORTHOGONAL PART OF THE DECOMPOSITION. C C JPVT JPVT(K) CONTAINS THE INDEX OF THE COLUMN OF THE C ORIGINAL MATRIX THAT HAS BEEN INTERCHANGED INTO C THE K-TH COLUMN, IF PIVOTING WAS REQUESTED. C C LINPACK. THIS VERSION DATED 08/14/78 . C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C C DQRDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS. C C BLAS DAXPY,DDOT,DSCAL,DSWAP,DNRM2 C FORTRAN DABS,DMAX1,MIN0,DSQRT C C INTERNAL VARIABLES C INTEGER J,JP,L,LP1,LUP,MAXJ,PL,PU DOUBLE PRECISION MAXNRM,DNRM2,TT DOUBLE PRECISION DDOT,NRMXL,T LOGICAL NEGJ,SWAPJ C C PL = 1 PU = 0 IF (JOB .EQ. 0) GO TO 60 C C PIVOTING HAS BEEN REQUESTED. REARRANGE THE COLUMNS C ACCORDING TO JPVT. C DO 20 J = 1, P SWAPJ = JPVT(J) .GT. 0 NEGJ = JPVT(J) .LT. 0 JPVT(J) = J IF (NEGJ) JPVT(J) = -J IF (.NOT.SWAPJ) GO TO 10 IF (J .NE. PL) CALL DSWAP(N,X(1,PL),1,X(1,J),1) JPVT(J) = JPVT(PL) JPVT(PL) = J PL = PL + 1 10 CONTINUE 20 CONTINUE PU = P DO 50 JJ = 1, P J = P - JJ + 1 IF (JPVT(J) .GE. 0) GO TO 40 JPVT(J) = -JPVT(J) IF (J .EQ. PU) GO TO 30 CALL DSWAP(N,X(1,PU),1,X(1,J),1) JP = JPVT(PU) JPVT(PU) = JPVT(J) JPVT(J) = JP 30 CONTINUE PU = PU - 1 40 CONTINUE 50 CONTINUE 60 CONTINUE C C COMPUTE THE NORMS OF THE FREE COLUMNS. C IF (PU .LT. PL) GO TO 80 DO 70 J = PL, PU QRAUX(J) = DNRM2(N,X(1,J),1) WORK(J) = QRAUX(J) 70 CONTINUE 80 CONTINUE C C PERFORM THE HOUSEHOLDER REDUCTION OF X. C LUP = MIN0(N,P) DO 200 L = 1, LUP IF (L .LT. PL .OR. L .GE. PU) GO TO 120 C C LOCATE THE COLUMN OF LARGEST NORM AND BRING IT C INTO THE PIVOT POSITION. C MAXNRM = 0.0D0 MAXJ = L DO 100 J = L, PU IF (QRAUX(J) .LE. MAXNRM) GO TO 90 MAXNRM = QRAUX(J) MAXJ = J 90 CONTINUE 100 CONTINUE IF (MAXJ .EQ. L) GO TO 110 CALL DSWAP(N,X(1,L),1,X(1,MAXJ),1) QRAUX(MAXJ) = QRAUX(L) WORK(MAXJ) = WORK(L) JP = JPVT(MAXJ) JPVT(MAXJ) = JPVT(L) JPVT(L) = JP 110 CONTINUE 120 CONTINUE QRAUX(L) = 0.0D0 IF (L .EQ. N) GO TO 190 C C COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L. C NRMXL = DNRM2(N-L+1,X(L,L),1) IF (NRMXL .EQ. 0.0D0) GO TO 180 IF (X(L,L) .NE. 0.0D0) NRMXL = DSIGN(NRMXL,X(L,L)) CALL DSCAL(N-L+1,1.0D0/NRMXL,X(L,L),1) X(L,L) = 1.0D0 + X(L,L) C C APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS, C UPDATING THE NORMS. C LP1 = L + 1 IF (P .LT. LP1) GO TO 170 DO 160 J = LP1, P T = -DDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L) CALL DAXPY(N-L+1,T,X(L,L),1,X(L,J),1) IF (J .LT. PL .OR. J .GT. PU) GO TO 150 IF (QRAUX(J) .EQ. 0.0D0) GO TO 150 TT = 1.0D0 - (DABS(X(L,J))/QRAUX(J))**2 TT = DMAX1(TT,0.0D0) T = TT TT = 1.0D0 + 0.05D0*TT*(QRAUX(J)/WORK(J))**2 IF (TT .EQ. 1.0D0) GO TO 130 QRAUX(J) = QRAUX(J)*DSQRT(T) GO TO 140 130 CONTINUE QRAUX(J) = DNRM2(N-L,X(L+1,J),1) WORK(J) = QRAUX(J) 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE C C SAVE THE TRANSFORMATION. C QRAUX(L) = X(L,L) X(L,L) = -NRMXL 180 CONTINUE 190 CONTINUE 200 CONTINUE RETURN END SUBROUTINE DSVDC(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,JOB,INFO) INTEGER LDX,N,P,LDU,LDV,JOB,INFO DOUBLE PRECISION X(LDX,1),S(1),E(1),U(LDU,1),V(LDV,1),WORK(1) C C C DSVDC IS A SUBROUTINE TO REDUCE A DOUBLE PRECISION NXP MATRIX X C BY ORTHOGONAL TRANSFORMATIONS U AND V TO DIAGONAL FORM. THE C DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X. THE C COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS, C AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS. C C ON ENTRY C C X DOUBLE PRECISION(LDX,P), WHERE LDX.GE.N. C X CONTAINS THE MATRIX WHOSE SINGULAR VALUE C DECOMPOSITION IS TO BE COMPUTED. X IS C DESTROYED BY DSVDC. C C LDX INTEGER. C LDX IS THE LEADING DIMENSION OF THE ARRAY X. C C N INTEGER. C N IS THE NUMBER OF COLUMNS OF THE MATRIX X. C C P INTEGER. C P IS THE NUMBER OF ROWS OF THE MATRIX X. C C LDU INTEGER. C LDU IS THE LEADING DIMENSION OF THE ARRAY U. C (SEE BELOW). C C LDV INTEGER. C LDV IS THE LEADING DIMENSION OF THE ARRAY V. C (SEE BELOW). C C WORK DOUBLE PRECISION(N). C WORK IS A SCRATCH ARRAY. C C JOB INTEGER. C JOB CONTROLS THE COMPUTATION OF THE SINGULAR C VECTORS. IT HAS THE DECIMAL EXPANSION AB C WITH THE FOLLOWING MEANING C C A.EQ.0 DO NOT COMPUTE THE LEFT SINGULAR C VECTORS. C A.EQ.1 RETURN THE N LEFT SINGULAR VECTORS C IN U. C A.GE.2 RETURN THE FIRST MIN(N,P) SINGULAR C VECTORS IN U. C B.EQ.0 DO NOT COMPUTE THE RIGHT SINGULAR C VECTORS. C B.EQ.1 RETURN THE RIGHT SINGULAR VECTORS C IN V. C C ON RETURN C C S DOUBLE PRECISION(MM), WHERE MM=MIN(N+1,P). C THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE C SINGULAR VALUES OF X ARRANGED IN DESCENDING C ORDER OF MAGNITUDE. C C E DOUBLE PRECISION(P). C E ORDINARILY CONTAINS ZEROS. HOWEVER SEE THE C DISCUSSION OF INFO FOR EXCEPTIONS. C C U DOUBLE PRECISION(LDU,K), WHERE LDU.GE.N. IF C JOBA.EQ.1 THEN K.EQ.N, IF JOBA.GE.2 C THEN K.EQ.MIN(N,P). C U CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS. C U IS NOT REFERENCED IF JOBA.EQ.0. IF N.LE.P C OR IF JOBA.EQ.2, THEN U MAY BE IDENTIFIED WITH X C IN THE SUBROUTINE CALL. C C V DOUBLE PRECISION(LDV,P), WHERE LDV.GE.P. C V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS. C V IS NOT REFERENCED IF JOB.EQ.0. IF P.LE.N, C THEN V MAY BE IDENTIFIED WITH X IN THE C SUBROUTINE CALL. C C INFO INTEGER. C THE SINGULAR VALUES (AND THEIR CORRESPONDING C SINGULAR VECTORS) S(INFO+1),S(INFO+2),...,S(M) C ARE CORRECT (HERE M=MIN(N,P)). THUS IF C INFO.EQ.0, ALL THE SINGULAR VALUES AND THEIR C VECTORS ARE CORRECT. IN ANY EVENT, THE MATRIX C B = TRANS(U)*X*V IS THE BIDIAGONAL MATRIX C WITH THE ELEMENTS OF S ON ITS DIAGONAL AND THE C ELEMENTS OF E ON ITS SUPER-DIAGONAL (TRANS(U) C IS THE TRANSPOSE OF U). THUS THE SINGULAR C VALUES OF X AND B ARE THE SAME. C C LINPACK. THIS VERSION DATED 03/19/79 . C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C C DSVDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS. C C EXTERNAL DROT C BLAS DAXPY,DDOT,DSCAL,DSWAP,DNRM2,DROTG C FORTRAN DABS,DMAX1,MAX0,MIN0,MOD,DSQRT C C INTERNAL VARIABLES C INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT, * MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1 DOUBLE PRECISION DDOT,T DOUBLE PRECISION B,C,CS,EL,EMM1,F,G,DNRM2,SCALE,SHIFT,SL,SM,SN, * SMM1,T1,TEST,ZTEST LOGICAL WANTU,WANTV C C C SET THE MAXIMUM NUMBER OF ITERATIONS. C MAXIT = 30 C C DETERMINE WHAT IS TO BE COMPUTED. C WANTU = .FALSE. WANTV = .FALSE. JOBU = MOD(JOB,100)/10 NCU = N IF (JOBU .GT. 1) NCU = MIN0(N,P) IF (JOBU .NE. 0) WANTU = .TRUE. IF (MOD(JOB,10) .NE. 0) WANTV = .TRUE. C C REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS C IN S AND THE SUPER-DIAGONAL ELEMENTS IN E. C INFO = 0 NCT = MIN0(N-1,P) NRT = MAX0(0,MIN0(P-2,N)) LU = MAX0(NCT,NRT) IF (LU .LT. 1) GO TO 170 DO 160 L = 1, LU LP1 = L + 1 IF (L .GT. NCT) GO TO 20 C C COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND C PLACE THE L-TH DIAGONAL IN S(L). C S(L) = DNRM2(N-L+1,X(L,L),1) IF (S(L) .EQ. 0.0D0) GO TO 10 IF (X(L,L) .NE. 0.0D0) S(L) = DSIGN(S(L),X(L,L)) CALL DSCAL(N-L+1,1.0D0/S(L),X(L,L),1) X(L,L) = 1.0D0 + X(L,L) 10 CONTINUE S(L) = -S(L) 20 CONTINUE IF (P .LT. LP1) GO TO 50 DO 40 J = LP1, P IF (L .GT. NCT) GO TO 30 IF (S(L) .EQ. 0.0D0) GO TO 30 C C APPLY THE TRANSFORMATION. C T = -DDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L) CALL DAXPY(N-L+1,T,X(L,L),1,X(L,J),1) 30 CONTINUE C C PLACE THE L-TH ROW OF X INTO E FOR THE C SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION. C E(J) = X(L,J) 40 CONTINUE 50 CONTINUE IF (.NOT.WANTU .OR. L .GT. NCT) GO TO 70 C C PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK C MULTIPLICATION. C DO 60 I = L, N U(I,L) = X(I,L) 60 CONTINUE 70 CONTINUE IF (L .GT. NRT) GO TO 150 C C COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE C L-TH SUPER-DIAGONAL IN E(L). C E(L) = DNRM2(P-L,E(LP1),1) IF (E(L) .EQ. 0.0D0) GO TO 80 IF (E(LP1) .NE. 0.0D0) E(L) = DSIGN(E(L),E(LP1)) CALL DSCAL(P-L,1.0D0/E(L),E(LP1),1) E(LP1) = 1.0D0 + E(LP1) 80 CONTINUE E(L) = -E(L) IF (LP1 .GT. N .OR. E(L) .EQ. 0.0D0) GO TO 120 C C APPLY THE TRANSFORMATION. C DO 90 I = LP1, N WORK(I) = 0.0D0 90 CONTINUE DO 100 J = LP1, P CALL DAXPY(N-L,E(J),X(LP1,J),1,WORK(LP1),1) 100 CONTINUE DO 110 J = LP1, P CALL DAXPY(N-L,-E(J)/E(LP1),WORK(LP1),1,X(LP1,J),1) 110 CONTINUE 120 CONTINUE IF (.NOT.WANTV) GO TO 140 C C PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT C BACK MULTIPLICATION. C DO 130 I = LP1, P V(I,L) = E(I) 130 CONTINUE 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE C C SET UP THE FINAL BIDIAGONAL MATRIX OR ORDER M. C M = MIN0(P,N+1) NCTP1 = NCT + 1 NRTP1 = NRT + 1 IF (NCT .LT. P) S(NCTP1) = X(NCTP1,NCTP1) IF (N .LT. M) S(M) = 0.0D0 IF (NRTP1 .LT. M) E(NRTP1) = X(NRTP1,M) E(M) = 0.0D0 C C IF REQUIRED, GENERATE U. C IF (.NOT.WANTU) GO TO 300 IF (NCU .LT. NCTP1) GO TO 200 DO 190 J = NCTP1, NCU DO 180 I = 1, N U(I,J) = 0.0D0 180 CONTINUE U(J,J) = 1.0D0 190 CONTINUE 200 CONTINUE IF (NCT .LT. 1) GO TO 290 DO 280 LL = 1, NCT L = NCT - LL + 1 IF (S(L) .EQ. 0.0D0) GO TO 250 LP1 = L + 1 IF (NCU .LT. LP1) GO TO 220 DO 210 J = LP1, NCU T = -DDOT(N-L+1,U(L,L),1,U(L,J),1)/U(L,L) CALL DAXPY(N-L+1,T,U(L,L),1,U(L,J),1) 210 CONTINUE 220 CONTINUE CALL DSCAL(N-L+1,-1.0D0,U(L,L),1) U(L,L) = 1.0D0 + U(L,L) LM1 = L - 1 IF (LM1 .LT. 1) GO TO 240 DO 230 I = 1, LM1 U(I,L) = 0.0D0 230 CONTINUE 240 CONTINUE GO TO 270 250 CONTINUE DO 260 I = 1, N U(I,L) = 0.0D0 260 CONTINUE U(L,L) = 1.0D0 270 CONTINUE 280 CONTINUE 290 CONTINUE 300 CONTINUE C C IF IT IS REQUIRED, GENERATE V. C IF (.NOT.WANTV) GO TO 350 DO 340 LL = 1, P L = P - LL + 1 LP1 = L + 1 IF (L .GT. NRT) GO TO 320 IF (E(L) .EQ. 0.0D0) GO TO 320 DO 310 J = LP1, P T = -DDOT(P-L,V(LP1,L),1,V(LP1,J),1)/V(LP1,L) CALL DAXPY(P-L,T,V(LP1,L),1,V(LP1,J),1) 310 CONTINUE 320 CONTINUE DO 330 I = 1, P V(I,L) = 0.0D0 330 CONTINUE V(L,L) = 1.0D0 340 CONTINUE 350 CONTINUE C C MAIN ITERATION LOOP FOR THE SINGULAR VALUES. C MM = M ITER = 0 360 CONTINUE C C QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND. C C ...EXIT IF (M .EQ. 0) GO TO 620 C C IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET C FLAG AND RETURN. C IF (ITER .LT. MAXIT) GO TO 370 INFO = M C ......EXIT GO TO 620 370 CONTINUE C C THIS SECTION OF THE PROGRAM INSPECTS FOR C NEGLIGIBLE ELEMENTS IN THE S AND E ARRAYS. ON C COMPLETION THE VARIABLES KASE AND L ARE SET AS FOLLOWS. C C KASE = 1 IF S(M) AND E(L-1) ARE NEGLIGIBLE AND L.LT.M C KASE = 2 IF S(L) IS NEGLIGIBLE AND L.LT.M C KASE = 3 IF E(L-1) IS NEGLIGIBLE, L.LT.M, AND C S(L), ..., S(M) ARE NOT NEGLIGIBLE (QR STEP). C KASE = 4 IF E(M-1) IS NEGLIGIBLE (CONVERGENCE). C DO 390 LL = 1, M L = M - LL C ...EXIT IF (L .EQ. 0) GO TO 400 TEST = DABS(S(L)) + DABS(S(L+1)) ZTEST = TEST + DABS(E(L)) IF (ZTEST .NE. TEST) GO TO 380 E(L) = 0.0D0 C ......EXIT GO TO 400 380 CONTINUE 390 CONTINUE 400 CONTINUE IF (L .NE. M - 1) GO TO 410 KASE = 4 GO TO 480 410 CONTINUE LP1 = L + 1 MP1 = M + 1 DO 430 LLS = LP1, MP1 LS = M - LLS + LP1 C ...EXIT IF (LS .EQ. L) GO TO 440 TEST = 0.0D0 IF (LS .NE. M) TEST = TEST + DABS(E(LS)) IF (LS .NE. L + 1) TEST = TEST + DABS(E(LS-1)) ZTEST = TEST + DABS(S(LS)) IF (ZTEST .NE. TEST) GO TO 420 S(LS) = 0.0D0 C ......EXIT GO TO 440 420 CONTINUE 430 CONTINUE 440 CONTINUE IF (LS .NE. L) GO TO 450 KASE = 3 GO TO 470 450 CONTINUE IF (LS .NE. M) GO TO 460 KASE = 1 GO TO 470 460 CONTINUE KASE = 2 L = LS 470 CONTINUE 480 CONTINUE L = L + 1 C C PERFORM THE TASK INDICATED BY KASE. C GO TO (490,520,540,570), KASE C C DEFLATE NEGLIGIBLE S(M). C 490 CONTINUE MM1 = M - 1 F = E(M-1) E(M-1) = 0.0D0 DO 510 KK = L, MM1 K = MM1 - KK + L T1 = S(K) CALL DROTG(T1,F,CS,SN) S(K) = T1 IF (K .EQ. L) GO TO 500 F = -SN*E(K-1) E(K-1) = CS*E(K-1) 500 CONTINUE IF (WANTV) CALL DROT(P,V(1,K),1,V(1,M),1,CS,SN) 510 CONTINUE GO TO 610 C C SPLIT AT NEGLIGIBLE S(L). C 520 CONTINUE F = E(L-1) E(L-1) = 0.0D0 DO 530 K = L, M T1 = S(K) CALL DROTG(T1,F,CS,SN) S(K) = T1 F = -SN*E(K) E(K) = CS*E(K) IF (WANTU) CALL DROT(N,U(1,K),1,U(1,L-1),1,CS,SN) 530 CONTINUE GO TO 610 C C PERFORM ONE QR STEP. C 540 CONTINUE C C CALCULATE THE SHIFT. C SCALE = DMAX1(DABS(S(M)),DABS(S(M-1)),DABS(E(M-1)), * DABS(S(L)),DABS(E(L))) SM = S(M)/SCALE SMM1 = S(M-1)/SCALE EMM1 = E(M-1)/SCALE SL = S(L)/SCALE EL = E(L)/SCALE B = ((SMM1 + SM)*(SMM1 - SM) + EMM1**2)/2.0D0 C = (SM*EMM1)**2 SHIFT = 0.0D0 IF (B .EQ. 0.0D0 .AND. C .EQ. 0.0D0) GO TO 550 SHIFT = DSQRT(B**2+C) IF (B .LT. 0.0D0) SHIFT = -SHIFT SHIFT = C/(B + SHIFT) 550 CONTINUE C C I changed this line from - shift to + shift C This bug was communicated by Ake Bjork C Mike Meyer jan 9 1984 C F = (SL + SM)*(SL - SM) - SHIFT C F = (SL + SM)*(SL - SM) + SHIFT G = SL*EL C C CHASE ZEROS. C MM1 = M - 1 DO 560 K = L, MM1 CALL DROTG(F,G,CS,SN) IF (K .NE. L) E(K-1) = F F = CS*S(K) + SN*E(K) E(K) = CS*E(K) - SN*S(K) G = SN*S(K+1) S(K+1) = CS*S(K+1) IF (WANTV) CALL DROT(P,V(1,K),1,V(1,K+1),1,CS,SN) CALL DROTG(F,G,CS,SN) S(K) = F F = CS*E(K) + SN*S(K+1) S(K+1) = -SN*E(K) + CS*S(K+1) G = SN*E(K+1) E(K+1) = CS*E(K+1) IF (WANTU .AND. K .LT. N) * CALL DROT(N,U(1,K),1,U(1,K+1),1,CS,SN) 560 CONTINUE E(M-1) = F ITER = ITER + 1 GO TO 610 C C CONVERGENCE. C 570 CONTINUE C C MAKE THE SINGULAR VALUE POSITIVE. C IF (S(L) .GE. 0.0D0) GO TO 580 S(L) = -S(L) IF (WANTV) CALL DSCAL(P,-1.0D0,V(1,L),1) 580 CONTINUE C C ORDER THE SINGULAR VALUE. C 590 IF (L .EQ. MM) GO TO 600 C ...EXIT IF (S(L) .GE. S(L+1)) GO TO 600 T = S(L) S(L) = S(L+1) S(L+1) = T IF (WANTV .AND. L .LT. P) * CALL DSWAP(P,V(1,L),1,V(1,L+1),1) IF (WANTU .AND. L .LT. N) * CALL DSWAP(N,U(1,L),1,U(1,L+1),1) L = L + 1 GO TO 590 600 CONTINUE ITER = 0 M = M - 1 610 CONTINUE GO TO 360 620 CONTINUE RETURN END C============================================================= subroutine netev(npar,theta, k,xmat,nxmax,nx,nxc, yhat) c c npar = # of variables (input integer) c theta = evaluation point (input vector of length npar) c f = function value (output scalar) c g = gradient value (output vector of length npar) c mode == 0 then compute f, leave g unchanged c == 1 then compute g, leave f unchanged c == 2 then compute both f and g c nstate is not used (optional flag for initialization) implicit double precision(a-h,o-z) parameter(kmax=8,jdmax=16) dimension gama(kmax,jdmax) dimension yhat(nxmax),xmat(nxmax,jdmax) dimension theta(npar) jd= nxc m0=k+1 j0=2*k+1 do 30 iunit=1,k loc0=j0+(iunit-1)*jd do 20 lag=1,jd gama(iunit,lag)=theta(loc0+lag) 20 continue 30 continue call net(yhat,xmat,theta(1),theta(m0+1),gama,nx,jd,k) return end subroutine getxy(fname,xmat,y,nxmax,nx,nxc) implicit double precision(a-h,o-z) character*35 fname dimension xmat(nxmax,10), y(nxmax) integer nx,nxc open(unit=11,file=fname,status='old') rewind(11) m=0 do 50 i=1,nxmax read(11,*,end=60) y(i), (xmat(i,j), j=1,nxc) m=m+1 50 continue 60 continue close(11) if (m.lt.nxmax) then nx=m else nx=nxmax endif return end subroutine canpar( theta, nu,nxc,ctheta) implicit double precision (a-h, o-z) integer ind(500),loc1,loc2,loc0 real*8 theta(1) real*8 ctheta(1) , sgn(500) np= nu*(nxc+2) +1 if( np.gt.500) then write(*,*) ' too m any parameters for canpar' stop endif ctheta(1) = theta(1) do j=1,nu ind(j)=j ctheta(j+ 1)= abs(theta(j+1)) if( theta(j+1).le.0) then sgn(j)=-1 else sgn(j)=1 endif enddo loc1= 1+nu call msort( ctheta( 2),ind,nu) do j=1, nu jp= ind(j) ctheta(j +loc1)= theta( jp+ loc1)*sgn(jp) enddo loc0= 2*nu+1 do j=1, nu jp= ind(j) loc1= (jp-1)* nxc + loc0 loc2= (j-1)*nxc + loc0 do jj=1,nxc ctheta(loc2+jj)= theta( loc1+ jj)*sgn(jp) enddo enddo return end subroutine msort( k,ki,n) C HEAPSORT ALGORITHM FOR SORTING ON VECTOR OF KEYS K OF LENGTH N C J F MONAHAN TRANSCRIBED FROM KNUTH, VOL 2, PP 146-7. C integer array ki is permuted along with K REAL*8 K(1),KK integer ki(1),kki INTEGER R IF(N.LE.1) RETURN L=N/2+1 R=N 2 IF(L.GT.1) GO TO 1 KK=K(R) kki= ki(R) K(R)=K(1) ki(R)=ki(1) R=R-1 IF(R.EQ.1) GO TO 9 GO TO 3 1 L=L-1 KK=K(L) kki=ki(L) 3 J=L 4 I=J J=2*J IF(J-R) 5,6,8 5 IF(K(J).LT.K(J+1)) J=J+1 6 IF(KK.GT.K(J)) GO TO 8 7 K(I)=K(J) ki(I)=ki(J) GO TO 4 8 K(I)=KK ki(I)=kki GO TO 2 9 K(1)=KK ki(1)=kki RETURN END subroutine tinit(theta,npar,eps,nreps) implicit double precision(a-h,o-z) parameter(kmax=8,jdmax=16) parameter(npmax=1+kmax*(jdmax+2)) dimension theta(npmax),ti(npmax),g(npmax) twoeps=2.*eps fmin=10000. do 30 i=1,nreps do 10 j=1,npar ti(j)=twoeps*(rng(1)-0.5) 10 continue call objfun(0,npar,ti,fi,g,1) if(fi.lt.fmin) then call dcopy(npar,ti,1,theta,1) fmin=fi endif 30 continue return end subroutine setidmat(h,n) implicit real*8(a-h,o-z) dimension h(n,n) do 5 i=1,n do 3 j=1,n h(i,j)=0. 3 continue h(i,i)=1. 5 continue return end FUNCTION RNG(IXX) implicit double precision (a-h,o-z) C UNIFORM PSEUDORANDOM NUMBER GENERATOR C FORTRAN VERSION OF LEWIS, GOODMAN, MILLER C SCHRAGE, ACM TOMS V.5 (1979) P132 C FIRST CALL SETS SEED TO IXX, LATER IXX IGNORED INTEGER A,P,IX,B15,B16,XHI,XALO,LEFTLO,FHI,K DATA A/16807/,B15/32768/,B16/65536/,P/2147483647/ DATA IX/0/ IF(IX.EQ.0) IX=IXX XHI=IX/B16 XALO=(IX-XHI*B16)*A LEFTLO=XALO/B16 FHI=XHI*A+LEFTLO K=FHI/B15 IX=(((XALO-LEFTLO*B16)-P)+(FHI-K*B15)*B16)+K IF(IX.LT.0) IX=IX+P RNG=FLOAT(IX)*4.656612875E-10 RETURN END C------ CONFUN & MAKEOPT are for use by NPSOL, otherwise not used SUBROUTINE CONFUN(MODE,NCNLN,N,NROWJ,NEEDC,X,C,CJAC,NSTATE) IMPLICIT double precision (A-H,O-Z) INTEGER*4 MODE,NCNLN,N,NROWJ INTEGER*4 NEEDC(*) REAL*8 X(N),C(*),CJAC(NROWJ,*) SAVE RETURN END subroutine makeopt(bigbnd,itmax1,itmax2,ftol1,ftol2,msglvl, * mverify) implicit double precision(a-h,o-z) OPEN(7,FILE='options1.dat',status='unknown') rewind(7) WRITE(7,7000) write(7,7010) WRITE(7,7001) BIGBND WRITE(7,7002) ITMAX1 WRITE(7,7004) FTOL1 WRITE(7,7005) MSGLVL WRITE(7,7006) MVERIFY WRITE(7,7900) close(7) OPEN(7,FILE='options2.dat',status='unknown') rewind(7) WRITE(7,7000) write(7,7010) WRITE(7,7001) BIGBND WRITE(7,7002) ITMAX2 WRITE(7,7004) FTOL2 WRITE(7,7005) MSGLVL WRITE(7,7006) MVERIFY WRITE(7,7900) close(7) ITMAX3=100 OPEN(7,FILE='options3.dat',status='unknown') rewind(7) WRITE(7,7000) write(7,7010) WRITE(7,7001) BIGBND WRITE(7,7002) ITMAX3 WRITE(7,7004) FTOL2 WRITE(7,7005) MSGLVL WRITE(7,7006) MVERIFY WRITE(7,7900) close(7) 7000 FORMAT('BEGIN',74X,' ') 7010 format(' NOLIST ',48x,' ') 7001 FORMAT(' INFINITE BOUND ',D15.3,35X,' ') 7002 FORMAT(' MAJOR ITERATION LIMIT ',I15 ,35X,' ') 7004 FORMAT(' OPTIMALITY TOLERANCE ',D15.3,35X,' ') 7005 FORMAT(' MAJOR PRINT LEVEL ',I15 ,35X,' ') 7006 FORMAT(' VERIFY LEVEL ',I15 ,35X,' ') 7008 FORMAT(' LINESEARCH TOLERANCE ',D15.3,35X,' ') 7900 FORMAT('END ',74X,' ') return end subroutine sdev(x,n1,n,sd,ave) c c Purpose: calculates standard deviation of values n1,n2,....n c in array x of length n c implicit double precision(a-h,o-z) dimension x(n) xn=float(n-n1+1) ave=0. do 10 i=n1,n ave=ave+x(i) 10 continue ave=ave/xn var=0. do 20 i=n1,n var=var+(x(i)-ave)**2 20 continue var=var/(xn-1.) sd=sqrt(var) return end integer*4 function igcdiv(i1,i2) implicit integer*4(i-n) igcd=0 imin=min(i1,i2) do 10 i=1,imin ir1=mod(i1,i) ir2=mod(i2,i) if (ir1.eq.0.and.ir2.eq.0) igcd=i 10 continue igcdiv=igcd return end subroutine recontp(xt,xdat,xlag,nxmax,jdmax,nx,ldelay, * jd,jt1,m,jtp) implicit double precision(a-h,o-z) dimension xt(nxmax),xdat(nxmax),xlag(nxmax,jdmax) m=0 do 10 i=1,nx-jt1 m=m+1 ipred=i+jt1 xt(i)=xdat(ipred) if (jd.gt.0) then do 5 j=1,jd iback=jtp+(j-1)*ldelay xlag(i,j)=xdat(ipred-iback) 5 continue endif 10 continue return end