1C Output from Public domain Ratfor, version 1.05
2      subroutine pfnb(n,p,m,a,y,q,r,b,band,m0,d,u,wn,wp, aa,yy,slo,shi,r
3     *hs,glob,ghib,nit,info)
4      integer n,p,m,kk(2),mm,m0,nit(5,m),info,sumbad
5      integer i,j,k,slo(n),shi(n),ifix,ibad
6      logical notopt
7      double precision a(p,n),y(n),r(n),q(m),d(n),u(n),b(p,m)
8      double precision wn(n,9), wp(p,(p+3)),band(n),qk(2)
9      double precision glob(p),ghib(p),aa(p,n),yy(n),rhs(p)
10      double precision zero,one,tau,beta,eps,big,fm,fn
11      parameter(zero = 0.0d0)
12      parameter(one = 1.0d0)
13      parameter(beta = 0.99995d0)
14      parameter(big = 1.0d+10)
15      parameter(eps = 1.0d-06)
16      do23000 iq = 1,m
17      notopt = .true.
18      tau = q(iq)
19      mm = m0
20      ifix = 0
21      ibad = 0
2223002 if(notopt)then
23      ibad = ibad + 1
24      fm = mm
25      fn = dble(n)
26      kk(1) = int(n * dmax1(1./fn, tau - fm/(2 * fn))) + 1
27      kk(2) = int(n * dmin1(tau + fm/(2. * fn), (fn - 1)/fn))
28      do23004 i = 1,n
29      u(i) = r(i)/band(i)
3023004 continue
3123005 continue
32      call kuantiles(kk,2,n,u)
33      qk(1) = u(kk(1))
34      qk(2) = u(kk(2))
35      call iphil(n,0,slo)
36      call iphil(n,0,shi)
37      do23006 i = 1,n
38      if(r(i) .lt. (band(i) * qk(1)))then
39      slo(i) = 1
40      else
41      if(r(i) .gt. (band(i) * qk(2)))then
42      shi(i) = 1
43      endif
44      endif
4523006 continue
4623007 continue
4723012 if(notopt)then
48      ifix = ifix + 1
49      call dphil(p,zero,glob)
50      call dphil(p,zero,ghib)
51      call dphil(n,one,d)
52      call dphil(n,one,u)
53      k = 0
54      do23014 i = 1,n
55      if(slo(i) .eq. 0 .and. shi(i) .eq. 0)then
56      k = k + 1
57      call dcopy(p,a(1,i),1,aa(1,k),1)
58      yy(k) = -y(i)
59      else
60      if(slo(i) .eq. 1)then
61      do23020 j = 1,p
62      glob(j) = glob(j) + a(j,i)
6323020 continue
6423021 continue
65      else
66      if(shi(i) .eq. 1)then
67      do23024 j = 1,p
68      ghib(j) = ghib(j) + a(j,i)
6923024 continue
7023025 continue
71      endif
72      endif
73      endif
7423014 continue
7523015 continue
76      call dcopy(p,glob,1,aa(1,k+1),1)
77      call dcopy(p,ghib,1,aa(1,k+2),1)
78      yy(k+1) = big
79      yy(k+2) = -big
80      call dgemv('N',p,k+2,one-tau,aa,p,d,1,zero,rhs,1)
81      call dscal(k+2,zero,wn,1)
82      call daxpy(k+2,one-tau,u,1,wn,1)
83      call rqfnb(k+2,p,aa,yy,rhs,d,u,beta,eps,wn,wp,nit(1,iq),info)
84      call dcopy(p,wp,1,b(1,iq),1)
85      call dcopy(n,y,1,r,1)
86      call dgemv('T',p,n,one,a,p,b(1,iq),1,one,r,1)
87      sumbad = 0
88      do23026 i = 1,n
89      if((r(i) .gt. 0) .and. slo(i) .eq. 1)then
90      slo(i) = 0
91      sumbad = sumbad + 1
92      endif
93      if((r(i) .lt. 0) .and. shi(i) .eq. 1)then
94      shi(i) = 0
95      sumbad = sumbad + 1
96      endif
9723026 continue
9823027 continue
99      if(sumbad .gt. 0)then
100      if(sumbad .gt. 0.1 * mm)then
101      mm = min(2 * mm, n)
102      goto 23013
103      endif
104      else
105      notopt = .false.
106      endif
107      goto 23012
108      endif
10923013 continue
110      nit(4,iq) = ifix
111      nit(5,iq) = ibad
112      goto 23002
113      endif
11423003 continue
11523000 continue
11623001 continue
117      return
118      end
119      subroutine iphil(n,a,v)
120      integer n,i,a,v(n)
121      do23036 i = 1,n
122      v(i) = a
12323036 continue
12423037 continue
125      return
126      end
127      subroutine dphil(n,a,v)
128      integer i,n
129      double precision a,v(n)
130      do23038 i = 1,n
131      v(i) = a
13223038 continue
13323039 continue
134      return
135      end
136