C Output from Public domain Ratfor, version 1.0 subroutine powell(n,p,p2,a,b,c,x,tau,h,f,u,s,g,d,xh,maxit,nflag) integer n,p,p2 double precision x(p),a(n,p),b(n),c(n) double precision f(n),u(p,p),s(n),g(p2),d(p),xh(p) double precision zero, one, mone, step, tau,pow integer h(p),hin,hout,k,it,inset,maxit,nflag parameter(zero = 0.0d0, one = 1.d0, mone = -1.d0) it = 0 23000 continue it = it + 1 if(it .gt. 1)then call pivot(n,p,h,hin,hout,a,u,d,xh,nflag) endif if(nflag .gt. 0)then nflag = nflag + 2 return endif do23007 i = 1,p xh(i) = b(h(i)) 23007 continue 23008 continue call dgemv('N',p,p,one,u,p,xh,1,zero,x,1) call dgemv('N',n,p,one,a,n,x,1,zero,f,1) do23009 i = 1,n if(inset(p,i,h) .gt. 0 .or. f(i) .gt. c(i))then s(i) = zero else if(b(i) .lt. f(i))then s(i) = one - tau else s(i) = - tau endif endif 23009 continue 23010 continue call dgemv('T',n,p,one,a,n,s,1,zero,xh,1) call dgemv('T',p,p,one,u,p,xh,1,zero,g,1) do23015 i = 1,p if(f(h(i)) .lt. c(h(i)))then if(b(h(i)) .lt. c(h(i)))then g(i + p) = - g(i) + one - tau else g(i + p) = - g(i) - tau endif else g(i + p) = - g(i) + tau endif g(i) = g(i) + one - tau 23015 continue 23016 continue k = idmin(p2,g,1) if(g(k) .ge. 0 .or. it .gt. maxit)then goto 23002 endif call dscal(p,zero,d,1) if(k .le. p)then call daxpy(p,one,u(1,k),1,d,1) else k = k - p call daxpy(p,mone,u(1,k),1,d,1) endif call dgemv('N',n,p,one,a,n,d,1,zero,s,1) do23025 i = 1,n call dcopy(p,x,1,xh,1) step = (b(i) - f(i))/s(i) call daxpy(p,step,d,1,xh,1) s(i) = pow(n,p,xh,a,b,c,tau) 23025 continue 23026 continue hin = idmin(n,s,1) if(inset(p,hin,h) .gt. 0)then nflag = 2 goto 23002 endif hout = h(k) 23001 goto 23000 23002 continue if(it .gt. maxit)then nflag = 1 endif return end subroutine pivot(n,p,h,hin,hout,a,b,u,v,eflag) integer n,p,h(p),hin,hout,inset,k,eflag double precision a(n,p),b(p,p),u(p),v(p) double precision zero,one parameter(zero = 0.d0, one = 1.d0) eflag = 0 k = inset(p,hout,h) if(k .eq. 0)then eflag = 1 return endif if(inset(p,hin,h) .gt. 0)then eflag = 2 return endif if(hin .lt. 1 .or. hin .gt. n)then eflag = 3 return endif call dcopy(p,a(hin,1),n,v,1) call dgemv('T',p,p,one,b,p,v,1,zero,u,1) call dcopy(p,b(1,k),1,v,1) do23037 j = 1,p do23039 i = 1,p if(j .eq. k)then b(i,j) = b(i,j)/u(k) else b(i,j) = b(i,j) - (u(j)/u(k)) * v(i) endif 23039 continue 23040 continue 23037 continue 23038 continue h(k) = hin return end integer function inset(p,k,h) integer p,k,h(p) do23043 inset = 1,p if(h(inset) .eq. k)then return endif 23043 continue 23044 continue inset = 0 return end double precision function pow(n,p,x,a,b,c,tau) integer n,p double precision x(p),a(n,p),b(n),c(n) double precision tau,u,zero,rho,fit,ddot parameter(zero= 0.d0) pow = zero do23047 i = 1,n fit = ddot(p,a(i,1),n,x,1) u = b(i) - min(fit,c(i)) pow = pow + rho(u, tau) 23047 continue 23048 continue return end double precision function rho(u,tau) double precision u,tau,one parameter(one = 1.d0) if(u .lt. 0)then rho = u * (tau - one) else rho = u * tau endif return end