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