1C Output from Public domain Ratfor, version 1.05
2      subroutine rqfnc(n1,n2,p,a1,y,a2,r,rhs,d1,d2,u,beta,eps,wn1,wn2,wp
3     *,nit,info)
4      integer n1,n2,p,info,nit(3)
5      double precision a1(p,n1),a2(p,n2),y(n1),r(n2),rhs(p),d1(n1),d2(n2
6     *),u(n1)
7      double precision wn1(n1,9),wn2(n2,6),wp(p,p+3)
8      double precision one,beta,eps
9      parameter(one = 1.0d0)
10      call lpfnc(n1,n2,p,a1,y,a2,r,rhs,d1,d2,u,beta,eps,wn1(1,1),wn2(1,1
11     *),wn1(1,2), wp(1,1),wn1(1,3),wn2(1,2),wn1(1,4),wn1(1,5),wn2(1,3),w
12     *n1(1,6), wp(1,2),wn1(1,7),wn2(1,4),wn1(1,8),wn1(1,9),wn2(1,5),wn2(
13     *1,6), wp(1,3),wp(1,4),nit,info)
14      return
15      end
16      subroutine lpfnc(n1,n2,p,a1,c1,a2,c2,b,d1,d2,u,beta,eps,x1,x2,s, y
17     *,z1,z2,w,dx1,dx2,ds,dy,dz1,dz2,dw,dr1,dr2,r2, rhs,ada,nit,info)
18      integer n1,p,i,info,nit(3),maxit
19      double precision a1(p,n1),a2(p,n2),c1(n1),c2(n2),b(p)
20      double precision zero,one,mone,big,ddot,dmax1,dmin1,dxdz1,dxdz2,ds
21     *dw
22      double precision deltap,deltad,beta,eps,mu,gap,g
23      double precision x1(n1),x2(n2),u(n1),s(n1),y(p),z1(n1),z2(n2),w(n1
24     *)
25      double precision d1(n1),d2(n2),rhs(p),ada(p,p)
26      double precision dx1(n1),dx2(n2),ds(n1),dy(p),dz1(n1),dz2(n2),dw(n
27     *1)
28      double precision dr1(n1),dr2(n2),r2(n2)
29      parameter(zero = 0.0d0)
30      parameter(one = 1.0d0)
31      parameter(mone = -1.0d0)
32      parameter(big = 1.0d+20)
33      parameter(maxit = 500)
34      nit(1)=0
35      nit(2)=0
36      nit(3)=n1
37      call dgemv('N',p,n1,one,a1,p,c1,1,zero,y,1)
38      do23000 i=1,n1
39      d1(i)=one
4023000 continue
4123001 continue
42      do23002 i=1,n2
43      d2(i)=zero
44      z2(i)=one
4523002 continue
4623003 continue
47      call stepy2(n1,n2,p,a1,d1,a2,d2,y,ada,info)
48      if(info .ne. 0)then
49      return
50      endif
51      call dcopy(n1,c1,1,s,1)
52      call dgemv('T',p,n1,mone,a1,p,y,1,one,s,1)
53      do23006 i=1,n1
54      if(dabs(s(i)) .lt. eps)then
55      z1(i)=dmax1(s(i),zero)+eps
56      w(i)=dmax1(-s(i),zero)+eps
57      else
58      z1(i)=dmax1(s(i),zero)
59      w(i)=dmax1(-s(i),zero)
60      endif
61      s(i)=u(i)-x1(i)
6223006 continue
6323007 continue
64      gap = ddot(n1,z1,1,x1,1)+ddot(n2,z2,1,x2,1)+ddot(n1,w,1,s,1)
6523010 if(gap .gt. eps .and. nit(1).lt.maxit)then
66      nit(1)=nit(1)+1
67      call dcopy(n2,c2,1,r2,1)
68      call dgemv('T',p,n2,mone,a2,p,y,1,one,r2,1)
69      call dcopy(p,b,1,dy,1)
70      call dgemv('N',p,n1,mone,a1,p,x1,1,one,dy,1)
71      call dgemv('N',p,n2,mone,a2,p,x2,1,one,dy,1)
72      do23012 i = 1,n1
73      d1(i)=one/(z1(i)/x1(i) + w(i)/s(i))
74      ds(i)=z1(i)-w(i)
75      dz1(i)=d1(i)*ds(i)
7623012 continue
7723013 continue
78      do23014 i = 1,n2
79      d2(i)=x2(i)/z2(i)
80      dz2(i)=d2(i)*r2(i)
8123014 continue
8223015 continue
83      call dgemv('N',p,n1,one,a1,p,dz1,1,one,dy,1)
84      call dgemv('N',p,n2,one,a2,p,dz2,1,one,dy,1)
85      call dcopy(p,dy,1,rhs,1)
86      call stepy2(n1,n2,p,a1,d1,a2,d2,dy,ada,info)
87      if(info .ne. 0)then
88      return
89      endif
90      call dgemv('T',p,n1,one,a1,p,dy,1,mone,ds,1)
91      deltap=big
92      deltad=big
93      do23018 i=1,n1
94      dx1(i)=d1(i)*ds(i)
95      ds(i)=-dx1(i)
96      dz1(i)=-z1(i)*(dx1(i)/x1(i) + one)
97      dw(i)=-w(i)*(ds(i)/s(i) + one)
98      if(dx1(i).lt.0)then
99      deltap=dmin1(deltap,-x1(i)/dx1(i))
100      endif
101      if(ds(i).lt.0)then
102      deltap=dmin1(deltap,-s(i)/ds(i))
103      endif
104      if(dz1(i).lt.0)then
105      deltad=dmin1(deltad,-z1(i)/dz1(i))
106      endif
107      if(dw(i).lt.0)then
108      deltad=dmin1(deltad,-w(i)/dw(i))
109      endif
11023018 continue
11123019 continue
112      call dcopy(n2,r2,1,dx2,1)
113      call dgemv('T',p,n2,one,a2,p,dy,1,mone,dx2,1)
114      do23028 i=1,n2
115      dx2(i)=d2(i)*dx2(i)
116      dz2(i)=-z2(i)*(dx2(i)/x2(i) + one)
117      if(dx2(i).lt.0)then
118      deltap=dmin1(deltap,-x2(i)/dx2(i))
119      endif
120      if(dz2(i).lt.0)then
121      deltad=dmin1(deltad,-z2(i)/dz2(i))
122      endif
12323028 continue
12423029 continue
125      deltap=dmin1(beta*deltap,one)
126      deltad=dmin1(beta*deltad,one)
127      if(min(deltap,deltad) .lt. one)then
128      nit(2)=nit(2)+1
129      mu = ddot(n1,x1,1,z1,1)+ddot(n2,x2,1,z2,1)+ddot(n1,s,1,w,1)
130      g = mu + deltap*ddot(n1,dx1,1,z1,1)+ deltad*ddot(n1,dz1,1,x1,1)+ d
131     *eltap*deltad*ddot(n1,dz1,1,dx1,1)+ deltap*ddot(n2,dx2,1,z2,1)+ del
132     *tad*ddot(n2,dz2,1,x2,1)+ deltap*deltad*ddot(n2,dz2,1,dx2,1)+ delta
133     *p*ddot(n1,ds,1,w,1)+ deltad*ddot(n1,dw,1,s,1) + deltap*deltad*ddot
134     *(n1,ds,1,dw,1)
135      mu = mu * ((g/mu)**3) /(dfloat(2*n1)+dfloat(n2))
136      do23036 i=1,n1
137      dsdw = ds(i)*dw(i)
138      dr1(i)=d1(i)*(mu*(one/s(i)-one/x1(i))+ dx1(i)*dz1(i)/x1(i)-dsdw/s(
139     *i))
14023036 continue
14123037 continue
142      do23038 i=1,n2
143      dr2(i)=d2(i)*(dx2(i)*dz2(i)/x2(i)-mu/x2(i))
14423038 continue
14523039 continue
146      call dswap(p,rhs,1,dy,1)
147      call dgemv('N',p,n1,one,a1,p,dr1,1,one,dy,1)
148      call dgemv('N',p,n2,one,a2,p,dr2,1,one,dy,1)
149      call dpotrs('U',p,1,ada,p,dy,p,info)
150      call dgemv('T',p,n1,one,a1,p,dy,1,zero,u,1)
151      deltap=big
152      deltad=big
153      do23040 i=1,n1
154      dsdw = ds(i)*dw(i)
155      dxdz1 = dx1(i)*dz1(i)
156      dx1(i) = d1(i)*(u(i)-z1(i)+w(i))-dr1(i)
157      ds(i) = -dx1(i)
158      dz1(i) = -z1(i)+(mu - z1(i)*dx1(i) - dxdz1)/x1(i)
159      dw(i) = -w(i)+(mu - w(i)*ds(i) - dsdw)/s(i)
160      if(dx1(i).lt.0)then
161      deltap=dmin1(deltap,-x1(i)/dx1(i))
162      endif
163      if(ds(i).lt.0)then
164      deltap=dmin1(deltap,-s(i)/ds(i))
165      endif
166      if(dz1(i).lt.0)then
167      deltad=dmin1(deltad,-z1(i)/dz1(i))
168      endif
169      if(dw(i).lt.0)then
170      deltad=dmin1(deltad,-w(i)/dw(i))
171      endif
17223040 continue
17323041 continue
174      call dgemv('T',p,n2,one,a2,p,dy,1,zero,u,1)
175      do23050 i=1,n2
176      dxdz2 = dx2(i)*dz2(i)
177      dx2(i) = d2(i)*(u(i)-r2(i))-dr2(i)
178      dz2(i) = -z2(i)+(mu - z2(i)*dx2(i) - dxdz2)/x2(i)
179      if(dx2(i).lt.0)then
180      deltap=dmin1(deltap,-x2(i)/dx2(i))
181      endif
182      if(dz2(i).lt.0)then
183      deltad=dmin1(deltad,-z2(i)/dz2(i))
184      endif
18523050 continue
18623051 continue
187      deltap=dmin1(beta*deltap,one)
188      deltad=dmin1(beta*deltad,one)
189      endif
190      call daxpy(n1,deltap,dx1,1,x1,1)
191      call daxpy(n2,deltap,dx2,1,x2,1)
192      call daxpy(n1,deltap,ds,1,s,1)
193      call daxpy(p,deltad,dy,1,y,1)
194      call daxpy(n1,deltad,dz1,1,z1,1)
195      call daxpy(n2,deltad,dz2,1,z2,1)
196      call daxpy(n1,deltad,dw,1,w,1)
197      gap = ddot(n1,z1,1,x1,1)+ddot(n2,z2,1,x2,1)+ddot(n1,w,1,s,1)
198      goto 23010
199      endif
20023011 continue
201      call daxpy(n1,mone,w,1,z1,1)
202      call dswap(n1,z1,1,x1,1)
203      return
204      end
205      subroutine stepy2(n1,n2,p,a1,d1,a2,d2,b,ada,info)
206      integer n1,n2,p,i,j,k,info
207      double precision a1(p,n1),a2(p,n2),b(p),d1(n1),d2(n2),ada(p,p),zer
208     *o
209      parameter(zero = 0.0d0)
210      do23056 j=1,p
211      do23058 k=1,p
212      ada(j,k)=zero
21323058 continue
21423059 continue
21523056 continue
21623057 continue
217      do23060 i=1,n1
218      call dsyr('U',p,d1(i),a1(1,i),1,ada,p)
21923060 continue
22023061 continue
221      do23062 i=1,n2
222      call dsyr('U',p,d2(i),a2(1,i),1,ada,p)
22323062 continue
22423063 continue
225      call dposv('U',p,1,ada,p,b,p,info)
226      return
227      end
228