C Output from Public domain Ratfor, version 1.0 subroutine rqbr(m,nn,m5,n3,n4,a,b,t,toler,ift,x,e,s,wa,wb,nsol,nds *ol,sol,dsol,lsol,h,qn,cutoff,ci,tnmat,big,lci1) integer i,j,k,kl,kount,kr,l,lsol,m,m1,m2,m3,m4,m5,ift integer n,n1,n2,n3,nsol,ndsol,out,s(m),h(nn,nsol) integer nn,n4,idxcf logical stage,test,init,iend,lup logical lci1,lci2,skip double precision a1,aux,b1,big,d,dif,pivot,smax,t,t0,t1,tnt double precision min,max,toler,zero,half,one,two double precision b(m),sol(n3,nsol),a(m,nn),x(nn),wa(m5,n4),wb(m) double precision sum,e(m),dsol(m,ndsol) double precision qn(nn),cutoff,ci(4,nn),tnmat(4,nn),tnew,told,tn parameter( zero = 0.d0) parameter( one = 1.d0) parameter( two = 2.d0) n=nn ift = 0 wa(m+2,nn+1) = one if(m5.ne.m+5)then ift = 3 endif if(n3.ne.n+3)then ift = 4 endif if(n4.ne.n+4)then ift = 5 endif if(m.le.zero.or.n.le.zero)then ift = 6 endif if(ift.le.two)then half = one/two iend = .true. lci2 = .false. lup = .true. skip = .false. idxcf = 0 tnew = zero tn = zero m1 = m+1 n1 = n+1 n2 = n+2 m2 = m+2 m3 = m+3 m4 = m+4 do23010 j = 1,n x(j) = zero 23010 continue 23011 continue do23012 i = 1,m e(i) = zero 23012 continue 23013 continue if(t.lt.zero.or.t.gt.one)then t0 = one/(two*float(m)) t1 = one-t0 t = t0 iend = .false. lci1 = .false. endif 23016 continue do23019 i = 1,m k = 1 do23021 j = 1,nn if(k.le.nn)then if(j.eq.idxcf)then skip = .true. else skip = .false. endif if(.not.skip)then wa(i,k) = a(i,j) k = k+1 endif endif 23021 continue 23022 continue wa(i,n4) = n+i wa(i,n2) = b(i) if(idxcf .ne. 0)then wa(i,n3) = tnew*a(i,idxcf) else wa(i,n3) = zero endif wa(i,n1) = wa(i,n2)-wa(i,n3) if(wa(i,n1).lt.zero)then do23033 j = 1,n4 wa(i,j) = -wa(i,j) 23033 continue 23034 continue endif 23019 continue 23020 continue do23035 j = 1,n wa(m4,j) = j wa(m2,j) = zero wa(m3,j) = zero do23037 i = 1,m aux = sign(one,wa(m4,j))*wa(i,j) wa(m2,j) = wa(m2,j)+aux*(one-sign(one,wa(i,n4))) wa(m3,j) = wa(m3,j)+aux*sign(one,wa(i,n4)) 23037 continue 23038 continue wa(m3,j) = two*wa(m3,j) 23035 continue 23036 continue dif = zero init = .false. if(.not.lci2)then do23041 k = 1,n wa(m5,k) = zero do23043 i = 1,m wa(m5,k) = wa(m5,k)+a(i,k) 23043 continue 23044 continue wa(m5,k) = wa(m5,k)/float(m) 23041 continue 23042 continue endif lsol = 1 kount = 0 23045 continue do23048 j = 1,n wa(m1,j) = wa(m2,j)+wa(m3,j)*t 23048 continue 23049 continue if(.not.init)then stage = .true. kr = 1 kl = 1 go to 30 endif 23052 continue stage = .false. 23055 continue max = -big do23058 j = kr,n d = wa(m1,j) if(d.lt.zero)then if(d.gt.(-two))then goto 23058 endif d = -d-two endif if(d.gt.max)then max = d in = j endif 23058 continue 23059 continue if(max.le.toler)then goto 23054 endif if(wa(m1,in).le.zero)then do23070 i = 1,m4 wa(i,in) = -wa(i,in) 23070 continue 23071 continue wa(m1,in) = wa(m1,in)-two wa(m2,in) = wa(m2,in)-two endif 23072 continue k = 0 do23075 i = kl,m d = wa(i,in) if(d.gt.toler)then k = k+1 wb(k) = wa(i,n1)/d s(k) = i test = .true. endif 23075 continue 23076 continue 23079 continue if(k.le.0)then test = .false. else min = big do23084 i = 1,k if(wb(i).lt.min)then j = i min = wb(i) out = s(i) endif 23084 continue 23085 continue wb(j) = wb(k) s(j) = s(k) k = k-1 endif if(.not.test.and.stage)then goto 23081 endif if(.not.test)then goto 23047 endif pivot = wa(out,in) if(wa(m1,in)-pivot-pivot.le.toler)then go to 10 endif do23094 j = kr,n3 d = wa(out,j) wa(m1,j) = wa(m1,j)-d-d wa(m2,j) = wa(m2,j)-d-d wa(out,j) = -d 23094 continue 23095 continue wa(out,n4) = -wa(out,n4) 23080 goto 23079 23081 continue do23096 i = 1,m4 d = wa(i,kr) wa(i,kr) = wa(i,in) wa(i,in) = d 23096 continue 23097 continue kr = kr+1 go to 20 10 do23098 j = kr,n3 if(j.ne.in)then wa(out,j) = wa(out,j)/pivot endif 23098 continue 23099 continue do23102 i = 1,m3 if(i.ne.out)then d = wa(i,in) do23106 j = kr,n3 if(j.ne.in)then wa(i,j) = wa(i,j)-d*wa(out,j) endif 23106 continue 23107 continue endif 23102 continue 23103 continue do23110 i = 1,m3 if(i.ne.out)then wa(i,in) = -wa(i,in)/pivot endif 23110 continue 23111 continue wa(out,in) = one/pivot d = wa(out,n4) wa(out,n4) = wa(m4,in) wa(m4,in) = d kount = kount+1 if(.not.stage)then goto 23074 endif kl = kl+1 do23116 j = kr,n4 d = wa(out,j) wa(out,j) = wa(kount,j) wa(kount,j) = d 23116 continue 23117 continue 20 if(kount+kr.eq.n1)then goto 23057 endif 30 max = -one do23120 j = kr,n if(abs(wa(m4,j)).le.n)then d = abs(wa(m1,j)) if(d.gt.max)then max = d in = j endif endif 23120 continue 23121 continue if(wa(m1,in).lt.zero)then do23128 i = 1,m4 wa(i,in) = -wa(i,in) 23128 continue 23129 continue endif 23073 goto 23072 23074 continue 23056 goto 23055 23057 continue 23053 goto 23052 23054 continue if(kr.eq.1)then do23132 j = 1,n d = abs(wa(m1,j)) if(d.le.toler.or.two-d.le.toler)then ift = 1 wa(m2,nn+1) = zero go to 80 endif 23132 continue 23133 continue endif 80 kount = 0 sum = zero if(.not.lci2)then do23138 i = 1,kl-1 k = wa(i,n4)*sign(one,wa(i,n4)) x(k) = wa(i,n1)*sign(one,wa(i,n4)) 23138 continue 23139 continue endif do23140 i = 1,n kd = abs(wa(m4,i))-n dsol(kd,lsol) = one+wa(m1,i)/two if(wa(m4,i).lt.zero)then dsol(kd,lsol) = one-dsol(kd,lsol) endif if(.not.lci2)then sum = sum + x(i)*wa(m5,i) sol(i+3,lsol) = x(i) h(i,lsol) = kd endif 23140 continue 23141 continue do23146 i = kl,m kd = abs(wa(i,n4))-n if(wa(i,n4).lt.zero)then dsol(kd,lsol) = zero endif if(wa(i,n4).gt.zero)then dsol(kd,lsol) = one endif 23146 continue 23147 continue if(.not.lci2)then sol(1,lsol) = smax sol(2,lsol) = sum sum = zero do23154 j=kl,m d = wa(j,n1)*sign(one,wa(j,n4)) sum = sum + d*(smax + half*(sign(one,d) - one)) 23154 continue 23155 continue sol(3,lsol) = sum do23156 i=1,m dsol(i,lsol+1) = dsol(i,lsol) 23156 continue 23157 continue endif if(lci2)then a1 = zero do23160 i = 1,m a1 = a1+a(i,idxcf)*(dsol(i,lsol)+t-one) 23160 continue 23161 continue tn = a1/sqrt(qn(idxcf)*t*(one-t)) if(abs(tn).lt.cutoff)then if(lup)then smax = big else smax = -big endif do23166 i =1,kl-1 k = wa(i,n4)*sign(one,wa(i,n4)) sol(k,1) = wa(i,n2)*sign(one,wa(i,n4)) sol(k,2) = wa(i,n3)*sign(one,wa(i,n4))/tnew 23166 continue 23167 continue do23168 i = kl,m a1 = zero b1 = zero k = wa(i,n4)*sign(one,wa(i,n4))-n l = 1 do23170 j = 1,n if(j.eq.idxcf)then l = l+1 endif a1 = a1 + a(k,l)*sol(j,1) b1 = b1 + a(k,l)*sol(j,2) l = l+1 23170 continue 23171 continue tnt = (b(k)-a1)/(a(k,idxcf)-b1) if(lup)then if(tnt.gt.tnew)then if(tnt.lt.smax)then smax = tnt out = i endif endif else if(tnt.lt.tnew)then if(tnt.gt.smax)then smax = tnt out = i endif endif endif 23168 continue 23169 continue if(lup)then told = tnew tnew = smax + toler ci(3,idxcf) = told - toler tnmat(3,idxcf) = tn if(.not.(tnew .lt. big-toler))then ci(3,idxcf) = big ci(4,idxcf) = big tnmat(3,idxcf) = tn tnmat(4,idxcf) = tn lup = .false. go to 70 endif else told = tnew tnew = smax - toler ci(2,idxcf) = told + toler tnmat(2,idxcf) = tn if(.not.(tnew .gt. -big+toler))then ci(2,idxcf) = -big ci(1,idxcf) = -big tnmat(2,idxcf) = tn tnmat(1,idxcf) = tn lup = .true. go to 60 endif endif do23190 i = 1,m wa(i,n3) = wa(i,n3)/told*tnew wa(i,n1) = wa(i,n2) - wa(i,n3) 23190 continue 23191 continue do23192 j = kr,n3 d = wa(out,j) wa(m1,j) = wa(m1,j) -d -d wa(m2,j) = wa(m2,j) -d -d wa(out,j) = -d 23192 continue 23193 continue wa(out,n4) = -wa(out,n4) init = .true. else if(lup)then ci(4,idxcf) = tnew - toler tnmat(4,idxcf) = tn lup = .false. go to 70 else ci(1,idxcf) = tnew + toler tnmat(1,idxcf) = tn lup = .true. go to 60 endif endif endif if((iend).and.(.not.lci2))then go to 40 endif if(.not.lci2)then init = .true. lsol = lsol+1 do23200 i = 1,m s(i) = zero 23200 continue 23201 continue do23202 j = 1,n x(j) = zero 23202 continue 23203 continue smax = two do23204 j = 1,n b1 = wa(m3,j) a1 = (-two-wa(m2,j))/b1 b1 = -wa(m2,j)/b1 if(a1.ge.t)then if(a1.lt.smax)then smax = a1 dif = (b1-a1)/two endif endif if(b1.gt.t)then if(b1.lt.smax)then smax = b1 dif = (b1-a1)/two endif endif 23204 continue 23205 continue tnt = smax+toler*(one+abs(dif)) if(tnt.ge.t1+toler)then iend = .true. endif t = tnt if(iend)then t = t1 endif endif 23046 goto 23045 23047 continue wa(m2,nn+1) = two ift = 2 go to 50 40 if(lsol.gt.2)then sol(1,1) = zero sol(3,1) = zero sol(1,lsol) = one sol(3,lsol) = zero do23220 i = 1,m dsol(i,1) = one dsol(i,lsol) = zero dsol(i,lsol+1) = zero 23220 continue 23221 continue endif l = kl-1 do23222 i = 1,l if(wa(i,n1).lt.zero)then do23226 j = kr,n4 wa(i,j) = -wa(i,j) 23226 continue 23227 continue endif 23222 continue 23223 continue 50 sum = zero if(.not.lci2)then do23230 i = kl,m k = wa(i,n4)*sign(one,wa(i,n4)) d = wa(i,n1)*sign(one,wa(i,n4)) sum = sum+d*sign(one,d)*(half+sign(one,d)*(t-half)) k = k-n e(k) = d 23230 continue 23231 continue wa(m2,n2) = kount wa(m1,n2) = n1-kr wa(m1,n1) = sum endif if(wa(m2,nn+1).eq.two)then goto 23018 endif if(.not.lci1)then goto 23018 endif if(.not.(.not.lci2))goto 23234 lci2 = .true. n = nn-1 n1 = n+1 n2 = n+2 n3 = n+3 n4 = n+4 60 idxcf = idxcf+1 if(.not.(idxcf.gt.nn))goto 23236 goto 23018 23236 continue 70 if(.not.(lup))goto 23238 tnew = x(idxcf)+toler told = tnew ci(3,idxcf) = x(idxcf) tnmat(3,idxcf) = zero goto 23239 23238 continue tnew = x(idxcf)-toler told = tnew ci(2,idxcf) = x(idxcf) tnmat(2,idxcf) = zero 23239 continue 23234 continue 23017 goto 23016 23018 continue do23242 i=1,m dsol(i,lsol) = dsol(i,lsol+1) 23242 continue 23243 continue endif return end