Skip to content
Browse files

WHT: whitespace clean up of some fortran code to prevent gfortran war…

…nings (convert tabs to spaces, remove some trailing whitespace)
  • Loading branch information...
1 parent 7c03d21 commit 00e546329440b77058e705cc20433de9954b14d2 Warren Weckesser committed Jan 29, 2012
Showing with 479 additions and 479 deletions.
  1. +245 −245 scipy/optimize/lbfgsb/routines.f
  2. +3 −3 scipy/stats/mvndst.f
  3. +231 −231 scipy/stats/statlib/ansari.f
View
490 scipy/optimize/lbfgsb/routines.f
@@ -339,7 +339,7 @@ subroutine mainlb(n, m, x, l, u, nbd, f, g, factr, pgtol, ws, wy,
c wy, of dimension n x m, stores Y, the matrix of y-vectors;
c sy, of dimension m x m, stores S'Y;
c ss, of dimension m x m, stores S'S;
-c yy, of dimension m x m, stores Y'Y;
+c yy, of dimension m x m, stores Y'Y;
c wt, of dimension m x m, stores the Cholesky factorization
c of (theta*S'S+LD^(-1)L'); see eq.
c (2.26) in [3].
@@ -356,7 +356,7 @@ subroutine mainlb(n, m, x, l, u, nbd, f, g, factr, pgtol, ws, wy,
c used to store the lower triangular part of
c N = [Y' ZZ'Y L_a'+R_z']
c [L_a +R_z S'AA'S ]
-c
+c
c z(n),r(n),d(n),t(n),wa(8*m) are double precision working arrays.
c z is used at different times to store the Cauchy point and
c the Newton point.
@@ -507,7 +507,7 @@ subroutine mainlb(n, m, x, l, u, nbd, f, g, factr, pgtol, ws, wy,
c Check the input arguments for errors.
- call errclb(n,m,factr,l,u,nbd,task,info,k)
+ call errclb(n,m,factr,l,u,nbd,task,info,k)
if (task(1:5) .eq. 'ERROR') then
call prn3lb(n,x,f,task,iprint,info,itfile,
+ iter,nfgv,nintol,nskip,nact,sbgnrm,
@@ -616,7 +616,7 @@ subroutine mainlb(n, m, x, l, u, nbd, f, g, factr, pgtol, ws, wy,
if (.not. cnstnd .and. col .gt. 0) then
c skip the search for GCP.
call dcopy(n,x,1,z,1)
- wrk = updatd
+ wrk = updatd
nint = 0
goto 333
endif
@@ -630,7 +630,7 @@ subroutine mainlb(n, m, x, l, u, nbd, f, g, factr, pgtol, ws, wy,
call timer(cpu1)
call cauchy(n,x,l,u,nbd,g,indx2,iwhere,t,d,z,
+ m,wy,ws,sy,wt,theta,col,head,
- + wa(1),wa(2*m+1),wa(4*m+1),wa(6*m+1),nint,
+ + wa(1),wa(2*m+1),wa(4*m+1),wa(6*m+1),nint,
+ sg,yg,iprint,sbgnrm,info,epsmch)
if (info .ne. 0) then
c singular triangular system detected; refresh the lbfgs memory.
@@ -732,7 +732,7 @@ subroutine mainlb(n, m, x, l, u, nbd, f, g, factr, pgtol, ws, wy,
c Generate the search direction d:=z-x.
do 40 i = 1, n
- d(i) = z(i) - x(i)
+ d(i) = z(i) - x(i)
40 continue
call timer(cpu1)
666 continue
@@ -773,7 +773,7 @@ subroutine mainlb(n, m, x, l, u, nbd, f, g, factr, pgtol, ws, wy,
endif
else if (task(1:5) .eq. 'FG_LN') then
c return to the driver for calculating f and g; reenter at 666.
- goto 1000
+ goto 1000
else
c calculate and print out the quantities related to the new X.
call timer(cpu2)
@@ -820,7 +820,7 @@ subroutine mainlb(n, m, x, l, u, nbd, f, g, factr, pgtol, ws, wy,
ddum = -gdold
else
dr = (gd - gdold)*stp
- call dscal(n,stp,d,1)
+ call dscal(n,stp,d,1)
ddum = -gdold*stp
endif
@@ -1003,15 +1003,15 @@ subroutine active(n, l, u, nbd, x, iwhere, iprint,
do 10 i = 1, n
if (nbd(i) .gt. 0) then
if (nbd(i) .le. 2 .and. x(i) .le. l(i)) then
- if (x(i) .lt. l(i)) then
+ if (x(i) .lt. l(i)) then
prjctd = .true.
- x(i) = l(i)
+ x(i) = l(i)
endif
nbdd = nbdd + 1
else if (nbd(i) .ge. 2 .and. x(i) .ge. u(i)) then
- if (x(i) .gt. u(i)) then
+ if (x(i) .gt. u(i)) then
prjctd = .true.
- x(i) = u(i)
+ x(i) = u(i)
endif
nbdd = nbdd + 1
endif
@@ -1024,14 +1024,14 @@ subroutine active(n, l, u, nbd, x, iwhere, iprint,
if (nbd(i) .ne. 2) boxed = .false.
if (nbd(i) .eq. 0) then
c this variable is always free
- iwhere(i) = -1
+ iwhere(i) = -1
c otherwise set x(i)=mid(x(i), u(i), l(i)).
else
- cnstnd = .true.
+ cnstnd = .true.
if (nbd(i) .eq. 2 .and. u(i) - l(i) .le. zero) then
c this variable is always fixed
- iwhere(i) = 3
+ iwhere(i) = 3
else
iwhere(i) = 0
endif
@@ -1065,9 +1065,9 @@ subroutine bmv(m, sy, wt, col, v, p, info)
c Subroutine bmv
c
c This subroutine computes the product of the 2m x 2m middle matrix
-c in the compact L-BFGS formula of B and a 2m vector v;
-c it returns the product in p.
-c
+c in the compact L-BFGS formula of B and a 2m vector v;
+c it returns the product in p.
+c
c m is an integer variable.
c On entry m is the maximum number of variable metric corrections
c used to define the limited memory matrix.
@@ -1126,7 +1126,7 @@ subroutine bmv(m, sy, wt, col, v, p, info)
c PART I: solve [ D^(1/2) O ] [ p1 ] = [ v1 ]
c [ -L*D^(-1/2) J ] [ p2 ] [ v2 ].
-c solve Jp2=v2+LD^(-1)v1.
+c solve Jp2=v2+LD^(-1)v1.
p(col + 1) = v(col + 1)
do 20 i = 2, col
i2 = col + i
@@ -1140,7 +1140,7 @@ subroutine bmv(m, sy, wt, col, v, p, info)
call dtrsl(wt,m,col,p(col+1),11,info)
if (info .ne. 0) return
-c solve D^(1/2)p1=v1.
+c solve D^(1/2)p1=v1.
do 30 i = 1, col
p(i) = v(i)/sqrt(sy(i,i))
30 continue
@@ -1375,7 +1375,7 @@ subroutine cauchy(n, x, l, u, nbd, g, iorder, iwhere, t, d, xcp,
if (sbgnrm .le. zero) then
if (iprint .ge. 0) write (6,*) 'Subgnorm = 0. GCP = X.'
call dcopy(n,x,1,xcp,1)
- return
+ return
endif
bnded = .true.
nfree = n + 1
@@ -1437,18 +1437,18 @@ subroutine cauchy(n, x, l, u, nbd, g, iorder, iwhere, t, d, xcp,
nbreak = nbreak + 1
iorder(nbreak) = i
t(nbreak) = tl/(-neggi)
- if (nbreak .eq. 1 .or. t(nbreak) .lt. bkmin) then
- bkmin = t(nbreak)
- ibkmin = nbreak
+ if (nbreak .eq. 1 .or. t(nbreak) .lt. bkmin) then
+ bkmin = t(nbreak)
+ ibkmin = nbreak
endif
else if (nbd(i) .ge. 2 .and. neggi .gt. zero) then
c x(i) + d(i) is bounded; compute t(i).
nbreak = nbreak + 1
iorder(nbreak) = i
t(nbreak) = tu/neggi
- if (nbreak .eq. 1 .or. t(nbreak) .lt. bkmin) then
- bkmin = t(nbreak)
- ibkmin = nbreak
+ if (nbreak .eq. 1 .or. t(nbreak) .lt. bkmin) then
+ bkmin = t(nbreak)
+ ibkmin = nbreak
endif
else
c x(i) + d(i) is not bounded.
@@ -1489,9 +1489,9 @@ subroutine cauchy(n, x, l, u, nbd, g, iorder, iwhere, t, d, xcp,
f2 = -theta*f1
f2_org = f2
if (col .gt. 0) then
- call bmv(m,sy,wt,col,p,v,info)
- if (info .ne. 0) return
- f2 = f2 - ddot(col2,v,1,p,1)
+ call bmv(m,sy,wt,col,p,v,info)
+ if (info .ne. 0) return
+ f2 = f2 - ddot(col2,v,1,p,1)
endif
dtm = -f1/f2
tsum = zero
@@ -1521,15 +1521,15 @@ subroutine cauchy(n, x, l, u, nbd, g, iorder, iwhere, t, d, xcp,
c Since we already have the smallest breakpoint we need not do
c heapsort yet. Often only one breakpoint is used and the
c cost of heapsort is avoided.
- tj = bkmin
- ibp = iorder(ibkmin)
+ tj = bkmin
+ ibp = iorder(ibkmin)
else
if (iter .eq. 2) then
c Replace the already used smallest breakpoint with the
c breakpoint numbered nbreak > nlast, before heapsort call.
if (ibkmin .ne. nbreak) then
t(ibkmin) = t(nbreak)
- iorder(ibkmin) = iorder(nbreak)
+ iorder(ibkmin) = iorder(nbreak)
endif
c Update heap structure of breakpoints
c (if iter=2, initialize heap).
@@ -1538,14 +1538,14 @@ subroutine cauchy(n, x, l, u, nbd, g, iorder, iwhere, t, d, xcp,
tj = t(nleft)
ibp = iorder(nleft)
endif
-
+
dt = tj - tj0
if (dt .ne. zero .and. iprint .ge. 100) then
write (6,4011) nint,f1,f2
write (6,5010) dt
write (6,6010) dtm
- endif
+ endif
c If a minimizer is within this interval, locate the GCP and return.
@@ -1560,20 +1560,20 @@ subroutine cauchy(n, x, l, u, nbd, g, iorder, iwhere, t, d, xcp,
dibp = d(ibp)
d(ibp) = zero
if (dibp .gt. zero) then
- zibp = u(ibp) - x(ibp)
- xcp(ibp) = u(ibp)
+ zibp = u(ibp) - x(ibp)
+ xcp(ibp) = u(ibp)
iwhere(ibp) = 2
else
- zibp = l(ibp) - x(ibp)
- xcp(ibp) = l(ibp)
+ zibp = l(ibp) - x(ibp)
+ xcp(ibp) = l(ibp)
iwhere(ibp) = 1
endif
if (iprint .ge. 100) write (6,*) 'Variable ',ibp,' is fixed.'
if (nleft .eq. 0 .and. nbreak .eq. n) then
c all n variables are fixed,
c return with xcp as GCP.
- dtm = dt
- goto 999
+ dtm = dt
+ goto 999
endif
c Update the derivative information.
@@ -1589,30 +1589,30 @@ subroutine cauchy(n, x, l, u, nbd, g, iorder, iwhere, t, d, xcp,
if (col .gt. 0) then
c update c = c + dt*p.
- call daxpy(col2,dt,p,1,c,1)
+ call daxpy(col2,dt,p,1,c,1)
c choose wbp,
c the row of W corresponding to the breakpoint encountered.
- pointr = head
+ pointr = head
do 70 j = 1,col
- wbp(j) = wy(ibp,pointr)
- wbp(col + j) = theta*ws(ibp,pointr)
+ wbp(j) = wy(ibp,pointr)
+ wbp(col + j) = theta*ws(ibp,pointr)
pointr = mod(pointr,m) + 1
70 continue
c compute (wbp)Mc, (wbp)Mp, and (wbp)M(wbp)'.
call bmv(m,sy,wt,col,wbp,v,info)
- if (info .ne. 0) return
- wmc = ddot(col2,c,1,v,1)
- wmp = ddot(col2,p,1,v,1)
- wmw = ddot(col2,wbp,1,v,1)
+ if (info .ne. 0) return
+ wmc = ddot(col2,c,1,v,1)
+ wmp = ddot(col2,p,1,v,1)
+ wmw = ddot(col2,wbp,1,v,1)
c update p = p - dibp*wbp.
- call daxpy(col2,-dibp,wbp,1,p,1)
+ call daxpy(col2,-dibp,wbp,1,p,1)
c complete updating f1 and f2 while col > 0.
- f1 = f1 + dibp*wmc
- f2 = f2 + 2.0d0*dibp*wmp - dibp2*wmw
+ f1 = f1 + dibp*wmc
+ f2 = f2 + 2.0d0*dibp*wmp - dibp2*wmw
endif
f2 = max(epsmch*f2_org,f2)
@@ -1621,9 +1621,9 @@ subroutine cauchy(n, x, l, u, nbd, g, iorder, iwhere, t, d, xcp,
goto 777
c to repeat the loop for unsearched intervals.
else if(bnded) then
- f1 = zero
- f2 = zero
- dtm = zero
+ f1 = zero
+ f2 = zero
+ dtm = zero
else
dtm = -f1/f2
endif
@@ -1707,27 +1707,27 @@ subroutine cmprlb(n, m, x, g, ws, wy, sy, wt, z, r, wa, index,
if (.not. cnstnd .and. col .gt. 0) then
do 26 i = 1, n
- r(i) = -g(i)
+ r(i) = -g(i)
26 continue
else
do 30 i = 1, nfree
k = index(i)
- r(i) = -theta*(z(k) - x(k)) - g(k)
+ r(i) = -theta*(z(k) - x(k)) - g(k)
30 continue
- call bmv(m,sy,wt,col,wa(2*m+1),wa(1),info)
+ call bmv(m,sy,wt,col,wa(2*m+1),wa(1),info)
if (info .ne. 0) then
info = -8
- return
+ return
endif
- pointr = head
- do 34 j = 1, col
- a1 = wa(j)
+ pointr = head
+ do 34 j = 1, col
+ a1 = wa(j)
a2 = theta*wa(col + j)
- do 32 i = 1, nfree
- k = index(i)
- r(i) = r(i) + wy(k,pointr)*a1 + ws(k,pointr)*a2
+ do 32 i = 1, nfree
+ k = index(i)
+ r(i) = r(i) + wy(k,pointr)*a1 + ws(k,pointr)*a2
32 continue
- pointr = mod(pointr,m) + 1
+ pointr = mod(pointr,m) + 1
34 continue
endif
@@ -1778,15 +1778,15 @@ subroutine errclb(n, m, factr, l, u, nbd, task, info, k)
if (nbd(i) .lt. 0 .or. nbd(i) .gt. 3) then
c return
task = 'ERROR: INVALID NBD'
- info = -6
- k = i
+ info = -6
+ k = i
endif
- if (nbd(i) .eq. 2) then
- if (l(i) .gt. u(i)) then
+ if (nbd(i) .eq. 2) then
+ if (l(i) .gt. u(i)) then
c return
task = 'ERROR: NO FEASIBLE SOLUTION'
- info = -7
- k = i
+ info = -7
+ k = i
endif
endif
10 continue
@@ -1947,90 +1947,90 @@ subroutine formk(n, nsub, ind, nenter, ileave, indx2, iupdat,
c shift old part of WN1.
do 10 jy = 1, m - 1
js = m + jy
- call dcopy(m-jy,wn1(jy+1,jy+1),1,wn1(jy,jy),1)
- call dcopy(m-jy,wn1(js+1,js+1),1,wn1(js,js),1)
- call dcopy(m-1,wn1(m+2,jy+1),1,wn1(m+1,jy),1)
+ call dcopy(m-jy,wn1(jy+1,jy+1),1,wn1(jy,jy),1)
+ call dcopy(m-jy,wn1(js+1,js+1),1,wn1(js,js),1)
+ call dcopy(m-1,wn1(m+2,jy+1),1,wn1(m+1,jy),1)
10 continue
endif
c put new rows in blocks (1,1), (2,1) and (2,2).
pbegin = 1
- pend = nsub
+ pend = nsub
dbegin = nsub + 1
- dend = n
+ dend = n
iy = col
is = m + col
ipntr = head + col - 1
- if (ipntr .gt. m) ipntr = ipntr - m
+ if (ipntr .gt. m) ipntr = ipntr - m
jpntr = head
do 20 jy = 1, col
js = m + jy
temp1 = zero
- temp2 = zero
- temp3 = zero
+ temp2 = zero
+ temp3 = zero
c compute element jy of row 'col' of Y'ZZ'Y
- do 15 k = pbegin, pend
- k1 = ind(k)
- temp1 = temp1 + wy(k1,ipntr)*wy(k1,jpntr)
+ do 15 k = pbegin, pend
+ k1 = ind(k)
+ temp1 = temp1 + wy(k1,ipntr)*wy(k1,jpntr)
15 continue
c compute elements jy of row 'col' of L_a and S'AA'S
- do 16 k = dbegin, dend
- k1 = ind(k)
- temp2 = temp2 + ws(k1,ipntr)*ws(k1,jpntr)
- temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr)
+ do 16 k = dbegin, dend
+ k1 = ind(k)
+ temp2 = temp2 + ws(k1,ipntr)*ws(k1,jpntr)
+ temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr)
16 continue
- wn1(iy,jy) = temp1
- wn1(is,js) = temp2
- wn1(is,jy) = temp3
+ wn1(iy,jy) = temp1
+ wn1(is,js) = temp2
+ wn1(is,jy) = temp3
jpntr = mod(jpntr,m) + 1
20 continue
c put new column in block (2,1).
- jy = col
+ jy = col
jpntr = head + col - 1
if (jpntr .gt. m) jpntr = jpntr - m
ipntr = head
do 30 i = 1, col
is = m + i
- temp3 = zero
+ temp3 = zero
c compute element i of column 'col' of R_z
- do 25 k = pbegin, pend
- k1 = ind(k)
- temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr)
+ do 25 k = pbegin, pend
+ k1 = ind(k)
+ temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr)
25 continue
- ipntr = mod(ipntr,m) + 1
+ ipntr = mod(ipntr,m) + 1
wn1(is,jy) = temp3
30 continue
- upcl = col - 1
+ upcl = col - 1
else
upcl = col
endif
c modify the old parts in blocks (1,1) and (2,2) due to changes
c in the set of free variables.
- ipntr = head
+ ipntr = head
do 45 iy = 1, upcl
is = m + iy
- jpntr = head
- do 40 jy = 1, iy
- js = m + jy
- temp1 = zero
- temp2 = zero
- temp3 = zero
- temp4 = zero
- do 35 k = 1, nenter
- k1 = indx2(k)
- temp1 = temp1 + wy(k1,ipntr)*wy(k1,jpntr)
- temp2 = temp2 + ws(k1,ipntr)*ws(k1,jpntr)
+ jpntr = head
+ do 40 jy = 1, iy
+ js = m + jy
+ temp1 = zero
+ temp2 = zero
+ temp3 = zero
+ temp4 = zero
+ do 35 k = 1, nenter
+ k1 = indx2(k)
+ temp1 = temp1 + wy(k1,ipntr)*wy(k1,jpntr)
+ temp2 = temp2 + ws(k1,ipntr)*ws(k1,jpntr)
35 continue
- do 36 k = ileave, n
- k1 = indx2(k)
- temp3 = temp3 + wy(k1,ipntr)*wy(k1,jpntr)
- temp4 = temp4 + ws(k1,ipntr)*ws(k1,jpntr)
+ do 36 k = ileave, n
+ k1 = indx2(k)
+ temp3 = temp3 + wy(k1,ipntr)*wy(k1,jpntr)
+ temp4 = temp4 + ws(k1,ipntr)*ws(k1,jpntr)
36 continue
- wn1(iy,jy) = wn1(iy,jy) + temp1 - temp3
- wn1(is,js) = wn1(is,js) - temp2 + temp4
- jpntr = mod(jpntr,m) + 1
+ wn1(iy,jy) = wn1(iy,jy) + temp1 - temp3
+ wn1(is,js) = wn1(is,js) - temp2 + temp4
+ jpntr = mod(jpntr,m) + 1
40 continue
ipntr = mod(ipntr,m) + 1
45 continue
@@ -2041,21 +2041,21 @@ subroutine formk(n, nsub, ind, nenter, ileave, indx2, iupdat,
jpntr = head
do 55 jy = 1, upcl
temp1 = zero
- temp3 = zero
- do 50 k = 1, nenter
- k1 = indx2(k)
- temp1 = temp1 + ws(k1,ipntr)*wy(k1,jpntr)
- 50 continue
- do 51 k = ileave, n
- k1 = indx2(k)
- temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr)
- 51 continue
+ temp3 = zero
+ do 50 k = 1, nenter
+ k1 = indx2(k)
+ temp1 = temp1 + ws(k1,ipntr)*wy(k1,jpntr)
+ 50 continue
+ do 51 k = ileave, n
+ k1 = indx2(k)
+ temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr)
+ 51 continue
if (is .le. jy + m) then
- wn1(is,jy) = wn1(is,jy) + temp1 - temp3
- else
- wn1(is,jy) = wn1(is,jy) - temp1 + temp3
- endif
- jpntr = mod(jpntr,m) + 1
+ wn1(is,jy) = wn1(is,jy) + temp1 - temp3
+ else
+ wn1(is,jy) = wn1(is,jy) - temp1 + temp3
+ endif
+ jpntr = mod(jpntr,m) + 1
55 continue
ipntr = mod(ipntr,m) + 1
60 continue
@@ -2065,21 +2065,21 @@ subroutine formk(n, nsub, ind, nenter, ileave, indx2, iupdat,
m2 = 2*m
do 70 iy = 1, col
- is = col + iy
- is1 = m + iy
- do 65 jy = 1, iy
- js = col + jy
+ is = col + iy
+ is1 = m + iy
+ do 65 jy = 1, iy
+ js = col + jy
js1 = m + jy
- wn(jy,iy) = wn1(iy,jy)/theta
- wn(js,is) = wn1(is1,js1)*theta
+ wn(jy,iy) = wn1(iy,jy)/theta
+ wn(js,is) = wn1(is1,js1)*theta
65 continue
- do 66 jy = 1, iy - 1
- wn(jy,is) = -wn1(is1,jy)
+ do 66 jy = 1, iy - 1
+ wn(jy,is) = -wn1(is1,jy)
66 continue
- do 67 jy = iy, col
- wn(jy,is) = wn1(is1,jy)
+ do 67 jy = iy, col
+ wn(jy,is) = wn1(is1,jy)
67 continue
- wn(iy,iy) = wn(iy,iy) + sy(iy,iy)
+ wn(iy,iy) = wn(iy,iy) + sy(iy,iy)
70 continue
c Form the upper triangle of WN= [ LL' L^-1(-L_a'+R_z')]
@@ -2089,8 +2089,8 @@ subroutine formk(n, nsub, ind, nenter, ileave, indx2, iupdat,
c with L' stored in the upper triangle of wn.
call dpofa(wn,m2,col,info)
if (info .ne. 0) then
- info = -1
- return
+ info = -1
+ return
endif
c then form L^-1(-L_a'+R_z') in the (1,2) block.
col2 = 2*col
@@ -2104,16 +2104,16 @@ subroutine formk(n, nsub, ind, nenter, ileave, indx2, iupdat,
do 72 is = col+1, col2
do 74 js = is, col2
- wn(is,js) = wn(is,js) + ddot(col,wn(1,is),1,wn(1,js),1)
+ wn(is,js) = wn(is,js) + ddot(col,wn(1,is),1,wn(1,js),1)
74 continue
72 continue
c Cholesky factorization of (2,2) block of wn.
call dpofa(wn(col+1,col+1),m2,col,info)
if (info .ne. 0) then
- info = -2
- return
+ info = -2
+ return
endif
return
@@ -2163,10 +2163,10 @@ subroutine formt(m, wt, sy, ss, col, theta, info)
c store T in the upper triangle of the array wt.
do 52 j = 1, col
- wt(1,j) = theta*ss(1,j)
+ wt(1,j) = theta*ss(1,j)
52 continue
do 55 i = 2, col
- do 54 j = i, col
+ do 54 j = i, col
k1 = min(i,j) - 1
ddum = zero
do 53 k = 1, k1
@@ -2241,21 +2241,21 @@ subroutine freev(n, nfree, index, nenter, ileave, indx2,
ileave = n + 1
if (iter .gt. 0 .and. cnstnd) then
c count the entering and leaving variables.
- do 20 i = 1, nfree
- k = index(i)
- if (iwhere(k) .gt. 0) then
- ileave = ileave - 1
- indx2(ileave) = k
- if (iprint .ge. 100) write (6,*)
+ do 20 i = 1, nfree
+ k = index(i)
+ if (iwhere(k) .gt. 0) then
+ ileave = ileave - 1
+ indx2(ileave) = k
+ if (iprint .ge. 100) write (6,*)
+ 'Variable ',k,' leaves the set of free variables'
endif
20 continue
- do 22 i = 1 + nfree, n
- k = index(i)
- if (iwhere(k) .le. 0) then
- nenter = nenter + 1
- indx2(nenter) = k
- if (iprint .ge. 100) write (6,*)
+ do 22 i = 1 + nfree, n
+ k = index(i)
+ if (iwhere(k) .le. 0) then
+ nenter = nenter + 1
+ indx2(nenter) = k
+ if (iprint .ge. 100) write (6,*)
+ 'Variable ',k,' enters the set of free variables'
endif
22 continue
@@ -2269,9 +2269,9 @@ subroutine freev(n, nfree, index, nenter, ileave, indx2,
nfree = 0
iact = n + 1
do 24 i = 1, n
- if (iwhere(i) .le. 0) then
- nfree = nfree + 1
- index(nfree) = i
+ if (iwhere(i) .le. 0) then
+ nfree = nfree + 1
+ index(nfree) = i
else
iact = iact - 1
index(iact) = i
@@ -2482,7 +2482,7 @@ subroutine lnsrlb(n, l, u, nbd, x, f, fold, gd, gdold, g, d, r, t,
if (iter .eq. 0 .and. .not. boxed) then
stp = min(one/dnorm, stpmx)
else
- stp = one
+ stp = one
endif
call dcopy(n,x,1,t,1)
@@ -2494,7 +2494,7 @@ subroutine lnsrlb(n, l, u, nbd, x, f, fold, gd, gdold, g, d, r, t,
556 continue
gd = ddot(n,g,1,d,1)
if (ifun .eq. 0) then
- gdold=gd
+ gdold=gd
if (gd .ge. zero) then
c the directional derivative >=0.
c Line search is impossible.
@@ -2507,15 +2507,15 @@ subroutine lnsrlb(n, l, u, nbd, x, f, fold, gd, gdold, g, d, r, t,
xstep = stp*dnorm
if (csave(1:4) .ne. 'CONV' .and. csave(1:4) .ne. 'WARN') then
- task = 'FG_LNSRCH'
- ifun = ifun + 1
+ task = 'FG_LNSRCH'
+ ifun = ifun + 1
nfgv = nfgv + 1
iback = ifun - 1
if (stp .eq. one) then
call dcopy(n,z,1,x,1)
else
do 41 i = 1, n
- x(i) = stp*d(i) + t(i)
+ x(i) = stp*d(i) + t(i)
41 continue
endif
else
@@ -2567,11 +2567,11 @@ subroutine matupd(n, m, ws, wy, sy, ss, d, r, itail,
c Set pointers for matrices WS and WY.
if (iupdat .le. m) then
- col = iupdat
- itail = mod(head+iupdat-2,m) + 1
+ col = iupdat
+ itail = mod(head+iupdat-2,m) + 1
else
- itail = mod(itail,m) + 1
- head = mod(head,m) + 1
+ itail = mod(itail,m) + 1
+ head = mod(head,m) + 1
endif
c Update matrices WS and WY.
@@ -2598,8 +2598,8 @@ subroutine matupd(n, m, ws, wy, sy, ss, d, r, itail,
c and the last column of SS:
pointr = head
do 51 j = 1, col - 1
- sy(col,j) = ddot(n,d,1,wy(1,pointr),1)
- ss(j,col) = ddot(n,ws(1,pointr),1,d,1)
+ sy(col,j) = ddot(n,d,1,wy(1,pointr),1)
+ ss(j,col) = ddot(n,ws(1,pointr),1,d,1)
pointr = mod(pointr,m) + 1
51 continue
if (stp .eq. one) then
@@ -2649,7 +2649,7 @@ subroutine prn1lb(n, m, l, u, x, iprint, itfile, epsmch)
if (iprint .ge. 1) then
write (itfile,2001) epsmch
write (itfile,*)'N = ',n,' M = ',m
- write (itfile,9001)
+ write (itfile,9001)
if (iprint .gt. 100) then
write (6,1004) 'L =',(l(i),i = 1,n)
write (6,1004) 'X0 =',(x(i),i = 1,n)
@@ -2719,20 +2719,20 @@ subroutine prn2lb(n, x, f, g, iprint, itfile, iter, nfgv, nact,
c 'word' records the status of subspace solutions.
if (iword .eq. 0) then
c the subspace minimization converged.
- word = 'con'
+ word = 'con'
else if (iword .eq. 1) then
c the subspace minimization stopped at a bound.
word = 'bnd'
else if (iword .eq. 5) then
c the truncated Newton step has been used.
- word = 'TNT'
+ word = 'TNT'
else
word = '---'
endif
if (iprint .ge. 99) then
write (6,*) 'LINE SEARCH',iback,' times; norm of step = ',xstep
write (6,2001) iter,f,sbgnrm
- if (iprint .gt. 100) then
+ if (iprint .gt. 100) then
write (6,1004) 'X =',(x(i), i = 1, n)
write (6,1004) 'G =',(g(i), i = 1, n)
endif
@@ -2811,7 +2811,7 @@ subroutine prn3lb(n, x, f, task, iprint, info, itfile,
if (info .eq. -5) write (6,9015)
if (info .eq. -6) write (6,*)' Input nbd(',k,') is invalid.'
if (info .eq. -7)
- + write (6,*)' l(',k,') > u(',k,'). No feasible solution.'
+ + write (6,*)' l(',k,') > u(',k,'). No feasible solution.'
if (info .eq. -8) write (6,9018)
if (info .eq. -9) write (6,9019)
endif
@@ -2919,15 +2919,15 @@ subroutine projgr(n, l, u, nbd, x, g, sbgnrm)
sbgnrm = zero
do 15 i = 1, n
- gi = g(i)
+ gi = g(i)
if (nbd(i) .ne. 0) then
if (gi .lt. zero) then
if (nbd(i) .ge. 2) gi = max((x(i)-u(i)),gi)
- else
+ else
if (nbd(i) .le. 2) gi = min((x(i)-l(i)),gi)
endif
endif
- sbgnrm = max(sbgnrm,abs(gi))
+ sbgnrm = max(sbgnrm,abs(gi))
15 continue
return
@@ -2951,24 +2951,24 @@ subroutine subsm(n, m, nsub, ind, l, u, nbd, x, d, ws, wy, theta,
c Subroutine subsm
c
c Given xcp, l, u, r, an index set that specifies
-c the active set at xcp, and an l-BFGS matrix B
-c (in terms of WY, WS, SY, WT, head, col, and theta),
-c this subroutine computes an approximate solution
-c of the subspace problem
+c the active set at xcp, and an l-BFGS matrix B
+c (in terms of WY, WS, SY, WT, head, col, and theta),
+c this subroutine computes an approximate solution
+c of the subspace problem
c
-c (P) min Q(x) = r'(x-xcp) + 1/2 (x-xcp)' B (x-xcp)
+c (P) min Q(x) = r'(x-xcp) + 1/2 (x-xcp)' B (x-xcp)
c
c subject to l<=x<=u
-c x_i=xcp_i for all i in A(xcp)
-c
-c along the subspace unconstrained Newton direction
-c
-c d = -(Z'BZ)^(-1) r.
+c x_i=xcp_i for all i in A(xcp)
+c
+c along the subspace unconstrained Newton direction
+c
+c d = -(Z'BZ)^(-1) r.
c
c The formula for the Newton direction, given the L-BFGS matrix
c and the Sherman-Morrison formula, is
c
-c d = (1/theta)r + (1/theta*2) Z'WK^(-1)W'Z r.
+c d = (1/theta)r + (1/theta*2) Z'WK^(-1)W'Z r.
c
c where
c K = [-D -Y'ZZ'Y/theta L_a'-R_z' ]
@@ -3106,16 +3106,16 @@ subroutine subsm(n, m, nsub, ind, l, u, nbd, x, d, ws, wy, theta,
pointr = head
do 20 i = 1, col
- temp1 = zero
- temp2 = zero
- do 10 j = 1, nsub
- k = ind(j)
- temp1 = temp1 + wy(k,pointr)*d(j)
- temp2 = temp2 + ws(k,pointr)*d(j)
+ temp1 = zero
+ temp2 = zero
+ do 10 j = 1, nsub
+ k = ind(j)
+ temp1 = temp1 + wy(k,pointr)*d(j)
+ temp2 = temp2 + ws(k,pointr)*d(j)
10 continue
- wv(i) = temp1
- wv(col + i) = theta*temp2
- pointr = mod(pointr,m) + 1
+ wv(i) = temp1
+ wv(col + i) = theta*temp2
+ pointr = mod(pointr,m) + 1
20 continue
c Compute wv:=K^(-1)wv.
@@ -3125,7 +3125,7 @@ subroutine subsm(n, m, nsub, ind, l, u, nbd, x, d, ws, wy, theta,
call dtrsl(wn,m2,col2,wv,11,info)
if (info .ne. 0) return
do 25 i = 1, col
- wv(i) = -wv(i)
+ wv(i) = -wv(i)
25 continue
call dtrsl(wn,m2,col2,wv,01,info)
if (info .ne. 0) return
@@ -3135,70 +3135,70 @@ subroutine subsm(n, m, nsub, ind, l, u, nbd, x, d, ws, wy, theta,
pointr = head
do 40 jy = 1, col
js = col + jy
- do 30 i = 1, nsub
- k = ind(i)
- d(i) = d(i) + wy(k,pointr)*wv(jy)/theta
+ do 30 i = 1, nsub
+ k = ind(i)
+ d(i) = d(i) + wy(k,pointr)*wv(jy)/theta
+ + ws(k,pointr)*wv(js)
30 continue
- pointr = mod(pointr,m) + 1
+ pointr = mod(pointr,m) + 1
40 continue
do 50 i = 1, nsub
- d(i) = d(i)/theta
+ d(i) = d(i)/theta
50 continue
c Backtrack to the feasible region.
alpha = one
- temp1 = alpha
+ temp1 = alpha
do 60 i = 1, nsub
- k = ind(i)
+ k = ind(i)
dk = d(i)
- if (nbd(k) .ne. 0) then
- if (dk .lt. zero .and. nbd(k) .le. 2) then
- temp2 = l(k) - x(k)
- if (temp2 .ge. zero) then
- temp1 = zero
- else if (dk*alpha .lt. temp2) then
- temp1 = temp2/dk
- endif
- else if (dk .gt. zero .and. nbd(k) .ge. 2) then
- temp2 = u(k) - x(k)
- if (temp2 .le. zero) then
- temp1 = zero
- else if (dk*alpha .gt. temp2) then
- temp1 = temp2/dk
- endif
+ if (nbd(k) .ne. 0) then
+ if (dk .lt. zero .and. nbd(k) .le. 2) then
+ temp2 = l(k) - x(k)
+ if (temp2 .ge. zero) then
+ temp1 = zero
+ else if (dk*alpha .lt. temp2) then
+ temp1 = temp2/dk
+ endif
+ else if (dk .gt. zero .and. nbd(k) .ge. 2) then
+ temp2 = u(k) - x(k)
+ if (temp2 .le. zero) then
+ temp1 = zero
+ else if (dk*alpha .gt. temp2) then
+ temp1 = temp2/dk
+ endif
endif
if (temp1 .lt. alpha) then
- alpha = temp1
- ibd = i
+ alpha = temp1
+ ibd = i
endif
endif
60 continue
if (alpha .lt. one) then
- dk = d(ibd)
- k = ind(ibd)
- if (dk .gt. zero) then
+ dk = d(ibd)
+ k = ind(ibd)
+ if (dk .gt. zero) then
x(k) = u(k)
d(ibd) = zero
- else if (dk .lt. zero) then
+ else if (dk .lt. zero) then
x(k) = l(k)
- d(ibd) = zero
- endif
+ d(ibd) = zero
+ endif
endif
do 70 i = 1, nsub
- k = ind(i)
- x(k) = x(k) + alpha*d(i)
+ k = ind(i)
+ x(k) = x(k) + alpha*d(i)
70 continue
if (iprint .ge. 99) then
- if (alpha .lt. one) then
+ if (alpha .lt. one) then
write (6,1002) alpha
else
write (6,*) 'SM solution inside the box'
- end if
- if (iprint .gt.100) write (6,1003) (x(i),i=1,n)
+ end if
+ if (iprint .gt.100) write (6,1003) (x(i),i=1,n)
endif
if (alpha .lt. one) then
@@ -3209,7 +3209,7 @@ subroutine subsm(n, m, nsub, ind, l, u, nbd, x, d, ws, wy, theta,
if (iprint .ge. 99) write (6,1004)
1001 format (/,'----------------SUBSM entered-----------------',/)
- 1002 format ( 'ALPHA = ',f7.5,' backtrack to the BOX')
+ 1002 format ( 'ALPHA = ',f7.5,' backtrack to the BOX')
1003 format ('Subspace solution X = ',/,(4x,1p,6(1x,d11.4)))
1004 format (/,'----------------exit SUBSM --------------------',/)
@@ -3282,13 +3282,13 @@ subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax,
c function at stp.
c On exit f is the value of the function at stp.
c
-c g is a double precision variable.
+c g is a double precision variable.
c On initial entry g is the derivative of the function at 0.
c On subsequent entries g is the derivative of the
c function at stp.
c On exit g is the derivative of the function at stp.
c
-c stp is a double precision variable.
+c stp is a double precision variable.
c On entry stp is the current estimate of a satisfactory
c step. On initial entry, a positive initial estimate
c must be provided.
@@ -3306,18 +3306,18 @@ subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax,
c curvature condition.
c On exit gtol is unchanged.
c
-c xtol is a double precision variable.
+c xtol is a double precision variable.
c On entry xtol specifies a nonnegative relative tolerance
c for an acceptable step. The subroutine exits with a
c warning if the relative difference between sty and stx
c is less than xtol.
c On exit xtol is unchanged.
c
-c stpmin is a double precision variable.
+c stpmin is a double precision variable.
c On entry stpmin is a nonnegative lower bound for the step.
c On exit stpmin is unchanged.
c
-c stpmax is a double precision variable.
+c stpmax is a double precision variable.
c On entry stpmax is a nonnegative upper bound for the step.
c On exit stpmax is unchanged.
c
@@ -3346,7 +3346,7 @@ subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax,
c
c Subprograms called
c
-c MINPACK-2 ... dcstep
+c MINPACK-2 ... dcstep
c
c MINPACK-1 Project. June 1983.
c Argonne National Laboratory.
View
6 scipy/stats/mvndst.f
@@ -1063,7 +1063,7 @@ DOUBLE PRECISION FUNCTION BVU( SH, SK, R )
BVN = -BVN/TWOPI
ENDIF
IF ( R .GT. 0 ) BVN = BVN + MVNPHI( -MAX( H, K ) )
- IF ( R .LT. 0 ) BVN = -BVN + MAX( ZERO, MVNPHI(-H)-MVNPHI(-K) )
+ IF ( R .LT. 0 ) BVN = -BVN + MAX( ZERO, MVNPHI(-H)-MVNPHI(-K) )
ENDIF
BVU = BVN
END
@@ -1089,8 +1089,8 @@ DOUBLE PRECISION FUNCTION MVNUNI()
PARAMETER ( INVMP1 = 4.656612873077392578125D-10 )
* INVMP1 = 1/(M1+1)
SAVE X10, X11, X12, X20, X21, X22
- DATA X10, X11, X12, X20, X21, X22
- & / 15485857, 17329489, 36312197, 55911127, 75906931, 96210113 /
+ DATA X10, X11, X12, X20, X21, X22
+ & / 15485857, 17329489, 36312197, 55911127, 75906931, 96210113 /
*
* Component 1
*
View
462 scipy/stats/statlib/ansari.f
@@ -21,8 +21,8 @@ subroutine wprob(test, other, astart, a1, l1, a2, a3, ifault)
nrows = 1 + (test * other)/2
sum = zero
do 10 i = 1, nrows
- sum = sum + a1(i)
- a1(i) = sum
+ sum = sum + a1(i)
+ a1(i) = sum
10 continue
do 20 i = 1, nrows
20 a1(i) = a1(i) / sum
@@ -32,280 +32,280 @@ subroutine wprob(test, other, astart, a1, l1, a2, a3, ifault)
c----------------------------------------------------------------------
- SUBROUTINE GSCALE(TEST, OTHER, ASTART, A1, L1, A2, A3, IFAULT)
+ SUBROUTINE GSCALE(TEST, OTHER, ASTART, A1, L1, A2, A3, IFAULT)
C
-C ALGORITHM AS 93 APPL. STATIST. (1976) VOL.25, NO.1
+C ALGORITHM AS 93 APPL. STATIST. (1976) VOL.25, NO.1
C
-C FROM THE SIZES OF TWO SAMPLES THE DISTRIBUTION OF THE
-C ANSARI-BRADLEY TEST FOR SCALE IS GENERATED IN ARRAY A1.
+C FROM THE SIZES OF TWO SAMPLES THE DISTRIBUTION OF THE
+C ANSARI-BRADLEY TEST FOR SCALE IS GENERATED IN ARRAY A1.
C
- REAL ASTART, A1(L1), A2(L1), A3(L1), AI, ONE, FPOINT
- INTEGER TEST, OTHER
- LOGICAL SYMM
- DATA ONE /1.0/
+ REAL ASTART, A1(L1), A2(L1), A3(L1), AI, ONE, FPOINT
+ INTEGER TEST, OTHER
+ LOGICAL SYMM
+ DATA ONE /1.0/
C
-C TYPE CONVERSION (EFFECT DEPENDS ON TYPE STATEMENT ABOVE).
+C TYPE CONVERSION (EFFECT DEPENDS ON TYPE STATEMENT ABOVE).
C
- FPOINT(I) = I
+ FPOINT(I) = I
C
-C CHECK PROBLEM SIZE AND DEFINE BASE VALUE OF THE DISTRIBUTION.
+C CHECK PROBLEM SIZE AND DEFINE BASE VALUE OF THE DISTRIBUTION.
C
- M = MIN0(TEST, OTHER)
- IFAULT = 2
- IF (M. LT. 0) RETURN
- ASTART = FPOINT((TEST + 1) / 2) * FPOINT(1 + TEST / 2)
- N = MAX0(TEST, OTHER)
+ M = MIN0(TEST, OTHER)
+ IFAULT = 2
+ IF (M. LT. 0) RETURN
+ ASTART = FPOINT((TEST + 1) / 2) * FPOINT(1 + TEST / 2)
+ N = MAX0(TEST, OTHER)
C
-C CHECK SIZE OF RESULT ARRAY.
+C CHECK SIZE OF RESULT ARRAY.
C
- IFAULT = 1
- LRES = 1 + (M * N) / 2
- IF (L1 .LT. LRES) RETURN
- SYMM = MOD(M + N, 2) .EQ. 0
+ IFAULT = 1
+ LRES = 1 + (M * N) / 2
+ IF (L1 .LT. LRES) RETURN
+ SYMM = MOD(M + N, 2) .EQ. 0
C
-C TREAT SMALL SAMPLES SEPARATELY.
+C TREAT SMALL SAMPLES SEPARATELY.
C
- MM1 = M - 1
- IF (M .GT. 2) GOTO 5
+ MM1 = M - 1
+ IF (M .GT. 2) GOTO 5
C
-C START-UP PROCEDURES ONLY NEEDED.
+C START-UP PROCEDURES ONLY NEEDED.
C
- IF (MM1) 1, 2, 3
+ IF (MM1) 1, 2, 3
C
-C ONE SAMPLE ONLY.
+C ONE SAMPLE ONLY.
C
-1 A1(1) = ONE
- GOTO 15
+1 A1(1) = ONE
+ GOTO 15
C
-C SMALLER SAMPLE SIZE = 1.
+C SMALLER SAMPLE SIZE = 1.
C
-2 CALL START1(N, A1, L1, LN1)
- GOTO 4
+2 CALL START1(N, A1, L1, LN1)
+ GOTO 4
C
-C SMALLER SAMPLE SIZE = 2.
+C SMALLER SAMPLE SIZE = 2.
C
-3 CALL START2(N, A1, L1, LN1)
+3 CALL START2(N, A1, L1, LN1)
C
-C RETURN IF A1 IS NOT IN REVERSE ORDER.
+C RETURN IF A1 IS NOT IN REVERSE ORDER.
C
-4 IF (SYMM .OR. (OTHER .GT. TEST)) GOTO 15
- GOTO 13
+4 IF (SYMM .OR. (OTHER .GT. TEST)) GOTO 15
+ GOTO 13
C
-C FULL GENERATOR NEEDED
-C SET UP INITIAL CONDITIONS (DEPENDS ON MOD(N, 2)).
+C FULL GENERATOR NEEDED
+C SET UP INITIAL CONDITIONS (DEPENDS ON MOD(N, 2)).
C
-5 NM1 = N - 1
- NM2 = N - 2
- MNOW = 3
- NC = 3
- IF (MOD(N, 2) .EQ. 1) GOTO 6
-C SET UP FOR EVEN N.
+5 NM1 = N - 1
+ NM2 = N - 2
+ MNOW = 3
+ NC = 3
+ IF (MOD(N, 2) .EQ. 1) GOTO 6
+C SET UP FOR EVEN N.
C
- N2B1 = 3
- N2B2 = 2
- CALL START2(N, A1, L1, LN1)
- CALL START2(NM2, A3, L1, LN3)
- CALL START1(NM1, A2, L1, LN2)
- GOTO 8
+ N2B1 = 3
+ N2B2 = 2
+ CALL START2(N, A1, L1, LN1)
+ CALL START2(NM2, A3, L1, LN3)
+ CALL START1(NM1, A2, L1, LN2)
+ GOTO 8
C
-C SET UP FOR ODD N.
+C SET UP FOR ODD N.
C
-6 N2B1 = 2
- N2B2 = 3
- CALL START1(N, A1, L1, LN1)
- CALL START2(NM1, A2, L1, LN2)
+6 N2B1 = 2
+ N2B2 = 3
+ CALL START1(N, A1, L1, LN1)
+ CALL START2(NM1, A2, L1, LN2)
C
-C INCREASE ORDER OF DISTRIBUTION IN A1 BY 2
-C (USING A2 AND IMPLYING A3).
+C INCREASE ORDER OF DISTRIBUTION IN A1 BY 2
+C (USING A2 AND IMPLYING A3).
C
-7 CALL FRQADD(A1, LN1, L1OUT, L1, A2, LN2, N2B1)
- LN1 = LN1 + N
- CALL IMPLY(A1, L1OUT, LN1, A3, LN3, L1, NC)
- NC = NC + 1
- IF (MNOW .EQ. M) GOTO 9
- MNOW = MNOW + 1
+7 CALL FRQADD(A1, LN1, L1OUT, L1, A2, LN2, N2B1)
+ LN1 = LN1 + N
+ CALL IMPLY(A1, L1OUT, LN1, A3, LN3, L1, NC)
+ NC = NC + 1
+ IF (MNOW .EQ. M) GOTO 9
+ MNOW = MNOW + 1
C
-C INCREASE ORDER OF DISTRIBUTION IN A2 BY 2 (USING A3).
+C INCREASE ORDER OF DISTRIBUTION IN A2 BY 2 (USING A3).
C
-8 CALL FRQADD(A2, LN2, L2OUT, L1, A3, LN3, N2B2)
- LN2 = LN2 + NM1
- CALL IMPLY(A2, L2OUT, LN2, A3, J, L1, NC)
- NC = NC + 1
- IF (MNOW .EQ. M) GOTO 9
- MNOW = MNOW + 1
- GOTO 7
+8 CALL FRQADD(A2, LN2, L2OUT, L1, A3, LN3, N2B2)
+ LN2 = LN2 + NM1
+ CALL IMPLY(A2, L2OUT, LN2, A3, J, L1, NC)
+ NC = NC + 1
+ IF (MNOW .EQ. M) GOTO 9
+ MNOW = MNOW + 1
+ GOTO 7
C
-C IF SYMMETRICAL, RESULTS IN A1 ARE COMPLETE.
+C IF SYMMETRICAL, RESULTS IN A1 ARE COMPLETE.
C
-9 IF (SYMM) GOTO 15
+9 IF (SYMM) GOTO 15
C
-C FOR A SKEW RESULT ADD A2 (OFFSET) INTO A1.
+C FOR A SKEW RESULT ADD A2 (OFFSET) INTO A1.
C
- KS = (M + 3) / 2
- J = 1
- DO 12 I = KS, LRES
- IF (I .GT. LN1) GOTO 10
- A1(I) = A1(I) + A2(J)
- GOTO 11
-10 A1(I) = A2(J)
-11 J = J + 1
-12 CONTINUE
+ KS = (M + 3) / 2
+ J = 1
+ DO 12 I = KS, LRES
+ IF (I .GT. LN1) GOTO 10
+ A1(I) = A1(I) + A2(J)
+ GOTO 11
+10 A1(I) = A2(J)
+11 J = J + 1
+12 CONTINUE
C
-C DISTRIBUTION IN A1 POSSIBLY IN REVERSE ORDER.
+C DISTRIBUTION IN A1 POSSIBLY IN REVERSE ORDER.
C
- IF (OTHER .LT. TEST) GOTO 15
+ IF (OTHER .LT. TEST) GOTO 15
C
-C REVERSE THE RESULTS IN A1.
+C REVERSE THE RESULTS IN A1.
C
-13 J = LRES
- NDO = LRES / 2
- DO 14 I = 1, NDO
- AI = A1(I)
- A1(I) =A1(J)
- A1(J) = AI
- J = J - 1
-14 CONTINUE
+13 J = LRES
+ NDO = LRES / 2
+ DO 14 I = 1, NDO
+ AI = A1(I)
+ A1(I) =A1(J)
+ A1(J) = AI
+ J = J - 1
+14 CONTINUE
C
-C FINAL RESULTS NOW IN A1.
+C FINAL RESULTS NOW IN A1.
C
-15 IFAULT = 0
- RETURN
- END
+15 IFAULT = 0
+ RETURN
+ END
- SUBROUTINE START1(N, F, L, LOUT)
-C
-C ALGORITHM AS 93.1 APPL. STATIST. (1976) VOL.25, NO.1
-C
-C GENERATES A 1,N ANSARI-BRADLEY DISTRIBUTION IN F.
+ SUBROUTINE START1(N, F, L, LOUT)
+C
+C ALGORITHM AS 93.1 APPL. STATIST. (1976) VOL.25, NO.1
+C
+C GENERATES A 1,N ANSARI-BRADLEY DISTRIBUTION IN F.
C
- REAL F(L), ONE, TWO
- DATA ONE, TWO /1.0, 2.0/
- LOUT = 1 + N / 2
- DO 1 I = 1, LOUT
-1 F(I) = TWO
- IF (MOD(N, 2) .EQ. 0) F(LOUT) = ONE
- RETURN
- END
+ REAL F(L), ONE, TWO
+ DATA ONE, TWO /1.0, 2.0/
+ LOUT = 1 + N / 2
+ DO 1 I = 1, LOUT
+1 F(I) = TWO
+ IF (MOD(N, 2) .EQ. 0) F(LOUT) = ONE
+ RETURN
+ END
C
- SUBROUTINE START2(N, F, L, LOUT)
+ SUBROUTINE START2(N, F, L, LOUT)
C
-C ALGORITHM AS 93.2 APPL. STATIST. (1976) VOL.25, NO.1
+C ALGORITHM AS 93.2 APPL. STATIST. (1976) VOL.25, NO.1
C
-C GENERATES A 2,N ANSARI-BRADLEY DISTRIBUTION IN F.
+C GENERATES A 2,N ANSARI-BRADLEY DISTRIBUTION IN F.
C
- REAL F(L), ONE, TWO, THREE, FOUR
- DATA ONE, TWO, THREE, FOUR /1.0, 2.0, 3.0, 4.0/
+ REAL F(L), ONE, TWO, THREE, FOUR
+ DATA ONE, TWO, THREE, FOUR /1.0, 2.0, 3.0, 4.0/
C
-C DERIVE F FOR 2, NU, WHERE NU IS HIGHEST EVEN INTEGER
-C LESS THAN OR EQUAL TO N.
-C DEFINE NU AND ARRAY LIMITS.
-C
- NU = N - MOD(N, 2)
- J = NU + 1
- LOUT = J
- LT1 = LOUT + 1
- NDO = LT1 / 2
- A = ONE
- B = THREE
-C
-C GENERATE THE SYMMETRICAL 2,NU DISTRIBUTION.
-C
- DO 1 I = 1, NDO
- F(I) = A
- F(J) = A
- J = J - 1
- A = A + B
- B = FOUR - B
-1 CONTINUE
- IF (NU .EQ. N) RETURN
-C
-C ADD AN OFFSET 1,N DISTRIBUTION INTO F TO GIVE 2,N RESULT.
-C
- NU = NDO + 1
- DO 2 I = NU, LOUT
-2 F(I) = F(I) + TWO
- F(LT1) = TWO
- LOUT = LT1
- RETURN
- END
-C
- SUBROUTINE FRQADD(F1, L1IN, L1OUT, L1, F2, L2, NSTART)
-C
-C ALGORITHM AS 93.3 APPL. STATIST. (1976) VOL.25, NO.1
-C
-C ARRAY F1 HAS TWICE THE CONTENTS OF ARRAY F2 ADDED INTO IT
-C STARTING WITH ELEMENTS NSTART AND 1 IN F1 AND F2 RESPECTIVELY.
-C
- REAL F1(L1), F2(L2), MUL2
- DATA MUL2 /2.0/
- I2 = 1
- DO 1 I1 = NSTART, L1IN
- F1(I1) = F1(I1) + MUL2 * F2(I2)
- I2 = I2 + 1
-1 CONTINUE
- NXT = L1IN + 1
- L1OUT = L2 + NSTART - 1
- DO 2 I1 = NXT, L1OUT
- F1(I1) = MUL2 * F2(I2)
- I2 = I2 + 1
-2 CONTINUE
- NSTART = NSTART + 1
- RETURN
- END
-C
- SUBROUTINE IMPLY(F1, L1IN, L1OUT, F2, L2, L2MAX, NOFF)
-C
-C ALGORITHM AS 93.4 APPL. STATIST. (1976) VOL.25, NO.1
-C
-C GIVEN L1IN ELEMENTS OF AN ARRAY F1, A SYMMETRICAL
-C ARRAY F2 IS DERIVED AND ADDED ONTO F1, LEAVING THE
-C FIRST NOFF ELEMENTS OF F1 UNCHANGED AND GIVING A
-C SYMMETRICAL RESULT OF L1OUT ELEMENTS IN F1.
-C
- REAL F1(L1OUT), F2(L2MAX), SUM, DIFF
-C
-C SET-UP SUBSCRIPTS AND LOOP COUNTER.
-C
- I2 = 1 - NOFF
- J1 = L1OUT
- J2 = L1OUT - NOFF
- L2 = J2
- J2MIN = (J2 + 1) / 2
- NDO = (L1OUT + 1) / 2
-C
-C DERIVE AND IMPLY NEW VALUES FROM OUTSIDE INWARDS.
-C
- DO 6 I1 = 1, NDO
-C
-C GET NEW F1 VALUE FROM SUM OF L/H ELEMENTS OF
-C F1 + F2 (IF F2 IS IN RANGE).
-C
- IF (I2 .GT. 0) GOTO 1
- SUM = F1(I1)
- GOTO 2
-1 SUM = F1(I1) + F2(I2)
-C
-C REVISE LEFT ELEMENT OF F1.
+C DERIVE F FOR 2, NU, WHERE NU IS HIGHEST EVEN INTEGER
+C LESS THAN OR EQUAL TO N.
+C DEFINE NU AND ARRAY LIMITS.
+C
+ NU = N - MOD(N, 2)
+ J = NU + 1
+ LOUT = J
+ LT1 = LOUT + 1
+ NDO = LT1 / 2
+ A = ONE
+ B = THREE
+C
+C GENERATE THE SYMMETRICAL 2,NU DISTRIBUTION.
+C
+ DO 1 I = 1, NDO
+ F(I) = A
+ F(J) = A
+ J = J - 1
+ A = A + B
+ B = FOUR - B
+1 CONTINUE
+ IF (NU .EQ. N) RETURN
+C
+C ADD AN OFFSET 1,N DISTRIBUTION INTO F TO GIVE 2,N RESULT.
+C
+ NU = NDO + 1
+ DO 2 I = NU, LOUT
+2 F(I) = F(I) + TWO
+ F(LT1) = TWO
+ LOUT = LT1
+ RETURN
+ END
+C
+ SUBROUTINE FRQADD(F1, L1IN, L1OUT, L1, F2, L2, NSTART)
+C
+C ALGORITHM AS 93.3 APPL. STATIST. (1976) VOL.25, NO.1
+C
+C ARRAY F1 HAS TWICE THE CONTENTS OF ARRAY F2 ADDED INTO IT
+C STARTING WITH ELEMENTS NSTART AND 1 IN F1 AND F2 RESPECTIVELY.
+C
+ REAL F1(L1), F2(L2), MUL2
+ DATA MUL2 /2.0/
+ I2 = 1
+ DO 1 I1 = NSTART, L1IN
+ F1(I1) = F1(I1) + MUL2 * F2(I2)
+ I2 = I2 + 1
+1 CONTINUE
+ NXT = L1IN + 1
+ L1OUT = L2 + NSTART - 1
+ DO 2 I1 = NXT, L1OUT
+ F1(I1) = MUL2 * F2(I2)
+ I2 = I2 + 1
+2 CONTINUE
+ NSTART = NSTART + 1
+ RETURN
+ END
+C
+ SUBROUTINE IMPLY(F1, L1IN, L1OUT, F2, L2, L2MAX, NOFF)
+C
+C ALGORITHM AS 93.4 APPL. STATIST. (1976) VOL.25, NO.1
+C
+C GIVEN L1IN ELEMENTS OF AN ARRAY F1, A SYMMETRICAL
+C ARRAY F2 IS DERIVED AND ADDED ONTO F1, LEAVING THE
+C FIRST NOFF ELEMENTS OF F1 UNCHANGED AND GIVING A
+C SYMMETRICAL RESULT OF L1OUT ELEMENTS IN F1.
+C
+ REAL F1(L1OUT), F2(L2MAX), SUM, DIFF
+C
+C SET-UP SUBSCRIPTS AND LOOP COUNTER.
+C
+ I2 = 1 - NOFF
+ J1 = L1OUT
+ J2 = L1OUT - NOFF
+ L2 = J2
+ J2MIN = (J2 + 1) / 2
+ NDO = (L1OUT + 1) / 2
+C
+C DERIVE AND IMPLY NEW VALUES FROM OUTSIDE INWARDS.
+C
+ DO 6 I1 = 1, NDO
+C
+C GET NEW F1 VALUE FROM SUM OF L/H ELEMENTS OF
+C F1 + F2 (IF F2 IS IN RANGE).
+C
+ IF (I2 .GT. 0) GOTO 1
+ SUM = F1(I1)
+ GOTO 2
+1 SUM = F1(I1) + F2(I2)
+C
+C REVISE LEFT ELEMENT OF F1.
C
- F1(I1) = SUM
-C
-C IF F2 NOT COMPLETE IMPLY AND ASSIGN F2 VALUES
-C AND REVISE SUBSCRIPTS.
-C
-2 I2 = I2 + 1
- IF (J2 .LT. J2MIN) GOTO 5
- IF (J1 .LE. L1IN) GOTO 3
- DIFF = SUM
- GOTO 4
-3 DIFF = SUM - F1(J1)
-4 F2(I1) = DIFF
- F2(J2) = DIFF
- J2 = J2 - 1
-C
-C ASSIGN R/H ELEMENT OF F1 AND REVISE SUBSCRIPT.
-C
-5 F1(J1) = SUM
- J1 = J1 - 1
-6 CONTINUE
- RETURN
- END
+ F1(I1) = SUM
+C
+C IF F2 NOT COMPLETE IMPLY AND ASSIGN F2 VALUES
+C AND REVISE SUBSCRIPTS.
+C
+2 I2 = I2 + 1
+ IF (J2 .LT. J2MIN) GOTO 5
+ IF (J1 .LE. L1IN) GOTO 3
+ DIFF = SUM
+ GOTO 4
+3 DIFF = SUM - F1(J1)
+4 F2(I1) = DIFF
+ F2(J2) = DIFF
+ J2 = J2 - 1
+C
+C ASSIGN R/H ELEMENT OF F1 AND REVISE SUBSCRIPT.
+C
+5 F1(J1) = SUM
+ J1 = J1 - 1
+6 CONTINUE
+ RETURN
+ END

0 comments on commit 00e5463

Please sign in to comment.
Something went wrong with that request. Please try again.