Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
c
c newGLMnet (5/12/14)
c
c
c Elastic net with squared-error loss
c
c dense predictor matrix:
c
c call elnet(ka,parm,no,ni,x,y,w,jd,vp,cl,ne,nx,nlam,flmin,ulam,thr,isd,
c intr,maxit,lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
c
c x(no,ni) = predictor data matrix flat file (overwritten)
c
c
c sparse predictor matrix:
c
c call spelnet(ka,parm,no,ni,x,ix,jx,y,w,jd,vp,cl,ne,nx,nlam,flmin,ulam,thr,
c isd,intr,maxit,lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
c
c x, ix, jx = predictor data matrix in compressed sparse row format
c
c
c other inputs:
c
c ka = algorithm flag
c ka=1 => covariance updating algorithm
c ka=2 => naive algorithm
c parm = penalty member index (0 <= parm <= 1)
c = 0.0 => ridge
c = 1.0 => lasso
c no = number of observations
c ni = number of predictor variables
c y(no) = response vector (overwritten)
c w(no)= observation weights (overwritten)
c jd(jd(1)+1) = predictor variable deletion flag
c jd(1) = 0 => use all variables
c jd(1) != 0 => do not use variables jd(2)...jd(jd(1)+1)
c vp(ni) = relative penalties for each predictor variable
c vp(j) = 0 => jth variable unpenalized
c cl(2,ni) = interval constraints on coefficient values (overwritten)
c cl(1,j) = lower bound for jth coefficient value (<= 0.0)
c cl(2,j) = upper bound for jth coefficient value (>= 0.0)
c ne = maximum number of variables allowed to enter largest model
c (stopping criterion)
c nx = maximum number of variables allowed to enter all models
c along path (memory allocation, nx > ne).
c nlam = (maximum) number of lamda values
c flmin = user control of lamda values (>=0)
c flmin < 1.0 => minimum lamda = flmin*(largest lamda value)
c flmin >= 1.0 => use supplied lamda values (see below)
c ulam(nlam) = user supplied lamda values (ignored if flmin < 1.0)
c thr = convergence threshold for each lamda solution.
c iterations stop when the maximum reduction in the criterion value
c as a result of each parameter update over a single pass
c is less than thr times the null criterion value.
c (suggested value, thr=1.0e-5)
c isd = predictor variable standarization flag:
c isd = 0 => regression on original predictor variables
c isd = 1 => regression on standardized predictor variables
c Note: output solutions always reference original
c variables locations and scales.
c intr = intercept flag
c intr = 0/1 => don't/do include intercept in model
c maxit = maximum allowed number of passes over the data for all lambda
c values (suggested values, maxit = 100000)
c
c output:
c
c lmu = actual number of lamda values (solutions)
c a0(lmu) = intercept values for each solution
c ca(nx,lmu) = compressed coefficient values for each solution
c ia(nx) = pointers to compressed coefficients
c nin(lmu) = number of compressed coefficients for each solution
c rsq(lmu) = R**2 values for each solution
c alm(lmu) = lamda values corresponding to each solution
c nlp = actual number of passes over the data for all lamda values
c jerr = error flag:
c jerr = 0 => no error
c jerr > 0 => fatal error - no output returned
c jerr < 7777 => memory allocation error
c jerr = 7777 => all used predictors have zero variance
c jerr = 10000 => maxval(vp) <= 0.0
C jerr < 0 => non fatal error - partial output:
c Solutions for larger lamdas (1:(k-1)) returned.
c jerr = -k => convergence for kth lamda value not reached
c after maxit (see above) iterations.
c jerr = -10000-k => number of non zero coefficients along path
c exceeds nx (see above) at kth lamda value.
c
c
c
c least-squares utility routines:
c
c
c uncompress coefficient vectors for all solutions:
c
c call solns(ni,nx,lmu,ca,ia,nin,b)
c
c input:
c
c ni,nx = input to elnet
c lmu,ca,ia,nin = output from elnet
c
c output:
c
c b(ni,lmu) = all elnet returned solutions in uncompressed format
c
c
c uncompress coefficient vector for particular solution:
c
c call uncomp(ni,ca,ia,nin,a)
c
c input:
c
c ni = total number of predictor variables
c ca(nx) = compressed coefficient values for the solution
c ia(nx) = pointers to compressed coefficients
c nin = number of compressed coefficients for the solution
c
c output:
c
c a(ni) = uncompressed coefficient vector
c referencing original variables
c
c
c evaluate linear model from compressed coefficients and
c uncompressed predictor matrix:
c
c call modval(a0,ca,ia,nin,n,x,f);
c
c input:
c
c a0 = intercept
c ca(nx) = compressed coefficient values for a solution
c ia(nx) = pointers to compressed coefficients
c nin = number of compressed coefficients for solution
c n = number of predictor vectors (observations)
c x(n,ni) = full (uncompressed) predictor matrix
c
c output:
c
c f(n) = model predictions
c
c
c evaluate linear model from compressed coefficients and
c compressed predictor matrix:
c
c call cmodval(a0,ca,ia,nin,x,ix,jx,n,f);
c
c input:
c
c a0 = intercept
c ca(nx) = compressed coefficient values for a solution
c ia(nx) = pointers to compressed coefficients
c nin = number of compressed coefficients for solution
c x, ix, jx = predictor matrix in compressed sparse row format
c n = number of predictor vectors (observations)
c
c output:
c
c f(n) = model predictions
c
c
c
c
c Multiple response
c elastic net with squared-error loss
c
c dense predictor matrix:
c
c call multelnet(parm,no,ni,nr,x,y,w,jd,vp,cl,ne,nx,nlam,flmin,ulam,thr,isd,
c jsd,intr,maxit,lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
c
c x(no,ni) = predictor data matrix flat file (overwritten)
c
c
c sparse predictor matrix:
c
c call multspelnet(parm,no,ni,nr,x,ix,jx,y,w,jd,vp,cl,ne,nx,nlam,flmin,ulam,thr,
c isd,jsd,intr,maxit,lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
c
c x, ix, jx = predictor data matrix in compressed sparse row format
c
c other inputs:
c
c nr = number of response variables
c y(no,nr) = response data matrix (overwritten)
c jsd = response variable standardization flag
c jsd = 0 => regression using original response variables
c jsd = 1 => regression using standardized response variables
c Note: output solutions always reference original
c variables locations and scales.
c all other inputs same as elnet/spelnet above
c
c output:
c
c a0(nr,lmu) = intercept values for each solution
c ca(nx,nr,lmu) = compressed coefficient values for each solution
c all other outputs same as elnet/spelnet above
c (jerr = 90000 => bounds adjustment non convergence)
c
c
c
c multiple response least-squares utility routines:
c
c
c uncompress coefficient matrix for all solutions:
c
c call multsolns(ni,nx,nr,lmu,ca,ia,nin,b)
c
c input:
c
c ni,nx,nr = input to multelnet
c lmu,ca,ia,nin = output from multelnet
c
c output:
c
c b(ni,nr,lmu) = all multelnet returned solutions in uncompressed format
c
c
c uncompress coefficient matrix for particular solution:
c
c call multuncomp(ni,nr,nx,ca,ia,nin,a)
c
c input:
c
c ni,nr,nx = input to multelnet
c ca(nx,nr) = compressed coefficient values for the solution
c ia(nx) = pointers to compressed coefficients
c nin = number of compressed coefficients for the solution
c
c output:
c
c a(ni,nr) = uncompressed coefficient matrix
c referencing original variables
c
c
c evaluate linear model from compressed coefficients and
c uncompressed predictor matrix:
c
c call multmodval(nx,nr,a0,ca,ia,nin,n,x,f);
c
c input:
c
c nx,nr = input to multelnet
c a0(nr) = intercepts
c ca(nx,nr) = compressed coefficient values for a solution
c ia(nx) = pointers to compressed coefficients
c nin = number of compressed coefficients for solution
c n = number of predictor vectors (observations)
c x(n,ni) = full (uncompressed) predictor matrix
c
c output:
c
c f(nr,n) = model predictions
c
c
c evaluate linear model from compressed coefficients and
c compressed predictor matrix:
c
c call multcmodval(nx,nr,a0,ca,ia,nin,x,ix,jx,n,f);
c
c input:
c
c nx,nr = input to multelnet
c a0(nr) = intercepts
c ca(nx,nr) = compressed coefficient values for a solution
c ia(nx) = pointers to compressed coefficients
c nin = number of compressed coefficients for solution
c x, ix, jx = predictor matrix in compressed sparse row format
c n = number of predictor vectors (observations)
c
c output:
c
c f(nr,n) = model predictions
c
c
c
c
c Symmetric binomial/multinomial logistic elastic net
c
c
c dense predictor matrix:
c
c call lognet (parm,no,ni,nc,x,y,o,jd,vp,cl,ne,nx,nlam,flmin,ulam,thr,isd,
c intr,maxit,kopt,lmu,a0,ca,ia,nin,dev0,fdev,alm,nlp,jerr)
c
c x(no,ni) = predictor data matrix flat file (overwritten)
c
c
c sparse predictor matrix:
c
c call splognet (parm,no,ni,nc,x,ix,jx,y,o,jd,vp,cl,ne,nx,nlam,flmin,
c ulam,thr,isd,intr,maxit,kopt,lmu,a0,ca,ia,nin,dev0,fdev,alm,nlp,jerr)
c
c x, ix, jx = predictor data matrix in compressed sparse row format
c
c
c other inputs:
c
c parm,no,ni,jd,vp,cl,ne,nx,nlam,flmin,ulam,thr,isd,intr,maxit
c = same as elnet above.
c
c nc = number of classes (distinct outcome values)
c nc=1 => binomial two-class logistic regression
c (all output references class 1)
c y(no,max(2,nc)) = number of each class at each design point
c entries may have fractional values or all be zero (overwritten)
c o(no,nc) = observation off-sets for each class
c kopt = optimization flag
c kopt = 0 => Newton-Raphson (recommended)
c kpot = 1 => modified Newton-Raphson (sometimes faster)
c kpot = 2 => nonzero coefficients same for each class (nc > 1)
c
c
c output:
c
c lmu,ia,nin,alm,nlp = same as elent above
c
c a0(nc,lmu) = intercept values for each class at each solution
c ca(nx,nc,lmu) = compressed coefficient values for each class at
c each solution
c dev0 = null deviance (intercept only model)
c fdev(lmu) = fraction of devience explained by each solution
c jerr = error flag
c jerr = 0 => no error
c jerr > 0 => fatal error - no output returned
c jerr < 7777 => memory allocation error
c jerr = 7777 => all used predictors have zero variance
c jerr = 8000 + k => null probability < 1.0e-5 for class k
c jerr = 9000 + k => null probability for class k
c > 1.0 - 1.0e-5
c jerr = 10000 => maxval(vp) <= 0.0
c jerr = 90000 => bounds adjustment non convergence
C jerr < 0 => non fatal error - partial output:
c Solutions for larger lamdas (1:(k-1)) returned.
c jerr = -k => convergence for kth lamda value not reached
c after maxit (see above) iterations.
c jerr = -10000-k => number of non zero coefficients along path
c exceeds nx (see above) at kth lamda value.
c jerr = -20000-k => max(p*(1-p)) < 1.0e-6 at kth lamda value.
c o(no,nc) = training data values for last (lmu_th) solution linear
c combination.
c
c
c
c logistic/multinomial utilitity routines:
c
c
c uncompress coefficient vectors for all solutions:
c
c call lsolns(ni,nx,nc,lmu,ca,ia,nin,b)
c
c input:
c
c ni,nx,nc = input to lognet
c lmu,ca,ia,nin = output from lognet
c
c output:
c
c b(ni,nc,lmu) = all lognet returned solutions in uncompressed format
c
c
c uncompress coefficient vector for particular solution:
c
c call luncomp(ni,nx,nc,ca,ia,nin,a)
c
c input:
c
c ni, nx, nc = same as above
c ca(nx,nc) = compressed coefficient values (for each class)
c ia(nx) = pointers to compressed coefficients
c nin = number of compressed coefficients
c
c output:
c
c a(ni,nc) = uncompressed coefficient vectors
c referencing original variables
c
c
c evaluate linear model from compressed coefficients and
c uncompressed predictor vectors:
c
c call lmodval(nt,x,nc,nx,a0,ca,ia,nin,ans);
c
c input:
c
c nt = number of observations
c x(nt,ni) = full (uncompressed) predictor vectors
c nc, nx = same as above
c a0(nc) = intercepts
c ca(nx,nc) = compressed coefficient values (for each class)
c ia(nx) = pointers to compressed coefficients
c nin = number of compressed coefficients
c
c output:
c
c ans(nc,nt) = model predictions
c
c
c evaluate linear model from compressed coefficients and
c compressed predictor matrix:
c
c call lcmodval(nc,nx,a0,ca,ia,nin,x,ix,jx,n,f);
c
c input:
c
c nc, nx = same as above
c a0(nc) = intercept
c ca(nx,nc) = compressed coefficient values for a solution
c ia(nx) = pointers to compressed coefficients
c nin = number of compressed coefficients for solution
c x, ix, jx = predictor matrix in compressed sparse row format
c n = number of predictor vectors (observations)
c
c output:
c
c f(nc,n) = model predictions
c
c
c
c
c Poisson elastic net
c
c
c dense predictor matrix:
c
c call fishnet (parm,no,ni,x,y,o,w,jd,vp,cl,ne,nx,nlam,flmin,ulam,thr,
c isd,intr,maxit,lmu,a0,ca,ia,nin,dev0,fdev,alm,nlp,jerr)
c
c x(no,ni) = predictor data matrix flat file (overwritten)
c
c sparse predictor matrix:
c
c call spfishnet (parm,no,ni,x,ix,jx,y,o,w,jd,vp,cl,ne,nx,nlam,flmin,ulam,thr,
c isd,intr,maxit,lmu,a0,ca,ia,nin,dev0,fdev,alm,nlp,jerr)
c
c x, ix, jx = predictor data matrix in compressed sparse row format
c
c other inputs:
c
c y(no) = observation response counts
c o(no) = observation off-sets
c parm,no,ni,w,jd,vp,cl,ne,nx,nlam,flmin,ulam,thr,isd,intr,maxit
c = same as elnet above
c
c output:
c
c lmu,a0,ca,ia,nin,alm = same as elnet above
c dev0,fdev = same as lognet above
c nlp = total number of passes over predictor variables
c jerr = error flag
c jerr = 0 => no error
c jerr > 0 => fatal error - no output returned
c jerr < 7777 => memory allocation error
c jerr = 7777 => all used predictors have zero variance
c jerr = 8888 => negative response count y values
c jerr = 9999 => no positive observations weights
c jerr = 10000 => maxval(vp) <= 0.0
C jerr < 0 => non fatal error - partial output:
c Solutions for larger lamdas (1:(k-1)) returned.
c jerr = -k => convergence for kth lamda value not reached
c after maxit (see above) iterations.
c jerr = -10000-k => number of non zero coefficients along path
c exceeds nx (see above) at kth lamda value.
c o(no) = training data values for last (lmu_th) solution linear
c combination.
c
c
c Poisson utility routines:
c
c
c same as elnet above:
c
c call solns(ni,nx,lmu,ca,ia,nin,b)
c call uncomp(ni,ca,ia,nin,a)
c call modval(a0,ca,ia,nin,n,x,f);
c call cmodval(a0,ca,ia,nin,x,ix,jx,n,f);
c
c compute deviance for given uncompressed data and set of uncompressed
c solutions
c
c call deviance(no,ni,x,y,o,w,nsol,a0,a,flog,jerr)
c
c input:
c
c no = number of observations
c ni = number of predictor variables
c x(no,ni) = predictor data matrix flat file
c y(no) = observation response counts
c o(no) = observation off-sets
c w(no)= observation weights
c nsol = number of solutions
c a0(nsol) = intercept for each solution
c a(ni,nsol) = solution coefficient vectors (uncompressed)
c
c output:
c
c flog(nsol) = respective deviance values minus null deviance
c jerr = error flag - see above
c
c
c compute deviance for given compressed data and set of uncompressed solutions
c
c call spdeviance(no,ni,x,ix,jx,y,o,w,nsol,a0,a,flog,jerr)
c
c input:
c
c no = number of observations
c ni = number of predictor variables
c x, ix, jx = predictor data matrix in compressed sparse row format
c y(no) = observation response counts
c o(no) = observation off-sets
c w(no)= observation weights
c nsol = number of solutions
c a0(nsol) = intercept for each solution
c a(ni,nsol) = solution coefficient vectors (uncompressed)
c
c output
c
c flog(nsol) = respective deviance values minus null deviance
c jerr = error flag - see above
c
c
c compute deviance for given compressed data and compressed solutions
c
c call cspdeviance(no,x,ix,jx,y,o,w,nx,lmu,a0,ca,ia,nin,flog,jerr)
c
c input:
c
c no = number of observations
c x, ix, jx = predictor data matrix in compressed sparse row format
c y(no) = observation response counts
c o(no) = observation off-sets
c w(no)= observation weights
c nx = input to spfishnet
c lmu,a0(lmu),ca(nx,lmu),ia(nx),nin(lmu) = output from spfishnet
c
c output
c
c flog(lmu) = respective deviance values minus null deviance
c jerr = error flag - see above
c
c
c
c Elastic net with Cox proportional hazards model
c
c
c dense predictor matrix:
c
c call coxnet (parm,no,ni,x,y,d,o,w,jd,vp,cl,ne,nx,nlam,flmin,ulam,thr,
c maxit,isd,lmu,ca,ia,nin,dev0,fdev,alm,nlp,jerr)
c
c input:
c
c x(no,ni) = predictor data matrix flat file (overwritten)
c y(no) = observation times
c d(no) = died/censored indicator
c d(i)=0.0 => y(i) = censoring time
c d(i)=1.0 => y(i) = death time
c o(no) = observation off-sets
c parm,no,ni,w,jd,vp,cl,ne,nx,nlam,flmin,ulam,thr,maxit
c = same as fishnet above
c
c output:
c
c lmu,ca,ia,nin,dev0,fdev,alm,nlp = same as fishnet above
c jerr = error flag
c jerr = 0 => no error - output returned
c jerr > 0 => fatal error - no output returned
c jerr < 7777 => memory allocation error
c jerr = 7777 => all used predictors have zero variance
c jerr = 8888 => all observations censored (d(i)=0.0)
c jerr = 9999 => no positive observations weights
c jerr = 10000 => maxval(vp) <= 0.0
c jerr = 20000, 30000 => initialization numerical error
C jerr < 0 => non fatal error - partial output:
c Solutions for larger lamdas (1:(k-1)) returned.
c jerr = -k => convergence for kth lamda value not reached
c after maxit (see above) iterations.
c jerr = -10000-k => number of non zero coefficients along path
c exceeds nx (see above) at kth lamda value.
c jerr = -30000-k => numerical error at kth lambda value
c o(no) = training data values for last (lmu_th) solution linear
c combination.
c
c
c
c coxnet utility routines:
c
c
c same as elnet above:
c
c call solns(ni,nx,lmu,ca,ia,nin,b)
c call uncomp(ni,ca,ia,nin,a)
c
c
c evaluate linear model from compressed coefficients and
c uncompressed predictor matrix:
c
c call cxmodval(ca,ia,nin,n,x,f);
c
c input:
c
c ca(nx) = compressed coefficient values for a solution
c ia(nx) = pointers to compressed coefficients
c nin = number of compressed coefficients for solution
c n = number of predictor vectors (observations)
c x(n,ni) = full (uncompressed) predictor matrix
c
c output:
c
c f(n) = model predictions
c
c
c compute log-likelihood for given data set and vectors of coefficients
c
c call loglike(no,ni,x,y,d,o,w,nvec,a,flog,jerr)
c
c input:
c
c no = number of observations
c ni = number of predictor variables
c x(no,ni) = predictor data matrix flat file
c y(no) = observation times
c d(no) = died/censored indicator
c d(i)=0.0 => y(i) = censoring time
c d(i)=1.0 => y(i) = death time
c o(no) = observation off-sets
c w(no)= observation weights
c nvec = number of coefficient vectors
c a(ni,nvec) = coefficient vectors (uncompressed)
c
c output
c
c flog(nvec) = respective log-likelihood values
c jerr = error flag - see coxnet above
c
c
c
c
c Changing internal parameter values
c
c
c call chg_fract_dev(fdev)
c fdev = minimum fractional change in deviance for stopping path
c default = 1.0e-5
c
c call chg_dev_max(devmax)
c devmax = maximum fraction of explained deviance for stopping path
c default = 0.999
c
c call chg_min_flmin(eps)
c eps = minimum value of flmin (see above). default= 1.0e-6
c
c call chg_big(big)
c big = large floating point number. default = 9.9e35
c
c call chg_min_lambdas(mnlam)
c mnlam = minimum number of path points (lambda values) allowed
c default = 5
c
c call chg_min_null_prob(pmin)
c pmin = minimum null probability for any class. default = 1.0e-9
c
c call chg _max_exp(exmx)
c exmx = maximum allowed exponent. default = 250.0
c
c call chg_bnorm(prec,mxit)
c prec = convergence threshold for multi response bounds adjustment
c solution. default = 1.0e-10.
c mxit = maximum iterations for multiresponse bounds adjustment solution
c default = 100.
c
c
c Obtain current internal parameter values
c
c call get_int_parms(fdev,eps,big,mnlam,devmax,pmin,exmx)
c call get_bnorm(prec,mxit);
c
c
c
subroutine get_int_parms(sml,eps,big,mnlam,rsqmax,pmin,exmx) 772
data sml0,eps0,big0,mnlam0,rsqmax0,pmin0,exmx0 /1.0e-5,1.0e-6,9.9 774
*e35,5,0.999,1.0e-9,250.0/
sml=sml0 774
eps=eps0 774
big=big0 774
mnlam=mnlam0 774
rsqmax=rsqmax0 775
pmin=pmin0 775
exmx=exmx0 776
return 777
entry chg_fract_dev(arg) 777
sml0=arg 777
return 778
entry chg_dev_max(arg) 778
rsqmax0=arg 778
return 779
entry chg_min_flmin(arg) 779
eps0=arg 779
return 780
entry chg_big(arg) 780
big0=arg 780
return 781
entry chg_min_lambdas(irg) 781
mnlam0=irg 781
return 782
entry chg_min_null_prob(arg) 782
pmin0=arg 782
return 783
entry chg_max_exp(arg) 783
exmx0=arg 783
return 784
end 785
subroutine elnet (ka,parm,no,ni,x,y,w,jd,vp,cl,ne,nx,nlam,flmin,u 788
*lam,thr,isd,intr,maxit, lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
real x(no,ni),y(no),w(no),vp(ni),ca(nx,nlam),cl(2,ni) 789
real ulam(nlam),a0(nlam),rsq(nlam),alm(nlam) 790
integer jd(*),ia(nx),nin(nlam) 791
real, dimension (:), allocatable :: vq;
if(maxval(vp) .gt. 0.0)goto 10021 794
jerr=10000 794
return 794
10021 continue 795
allocate(vq(1:ni),stat=jerr) 795
if(jerr.ne.0) return 796
vq=max(0.0,vp) 796
vq=vq*ni/sum(vq) 797
if(ka .ne. 1)goto 10041 798
call elnetu (parm,no,ni,x,y,w,jd,vq,cl,ne,nx,nlam,flmin,ulam,thr, 801
*isd,intr,maxit, lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
goto 10051 802
10041 continue 803
call elnetn (parm,no,ni,x,y,w,jd,vq,cl,ne,nx,nlam,flmin,ulam,thr,i 806
*sd,intr,maxit, lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
10051 continue 807
10031 continue 807
deallocate(vq) 808
return 809
end 810
subroutine elnetu (parm,no,ni,x,y,w,jd,vp,cl,ne,nx,nlam,flmin,ula 813
*m,thr,isd,intr,maxit, lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
real x(no,ni),y(no),w(no),vp(ni),ulam(nlam),cl(2,ni) 814
real ca(nx,nlam),a0(nlam),rsq(nlam),alm(nlam) 815
integer jd(*),ia(nx),nin(nlam) 816
real, dimension (:), allocatable :: xm,xs,g,xv,vlam
integer, dimension (:), allocatable :: ju
allocate(g(1:ni),stat=jerr) 821
allocate(xm(1:ni),stat=ierr) 821
jerr=jerr+ierr 822
allocate(xs(1:ni),stat=ierr) 822
jerr=jerr+ierr 823
allocate(ju(1:ni),stat=ierr) 823
jerr=jerr+ierr 824
allocate(xv(1:ni),stat=ierr) 824
jerr=jerr+ierr 825
allocate(vlam(1:nlam),stat=ierr) 825
jerr=jerr+ierr 826
if(jerr.ne.0) return 827
call chkvars(no,ni,x,ju) 828
if(jd(1).gt.0) ju(jd(2:(jd(1)+1)))=0 829
if(maxval(ju) .gt. 0)goto 10071 829
jerr=7777 829
return 829
10071 continue 830
call standard(no,ni,x,y,w,isd,intr,ju,g,xm,xs,ym,ys,xv,jerr) 831
if(jerr.ne.0) return 832
cl=cl/ys 832
if(isd .le. 0)goto 10091 832
10100 do 10101 j=1,ni 832
cl(:,j)=cl(:,j)*xs(j) 832
10101 continue 832
10102 continue 832
10091 continue 833
if(flmin.ge.1.0) vlam=ulam/ys 834
call elnet1(parm,ni,ju,vp,cl,g,no,ne,nx,x,nlam,flmin,vlam,thr,maxi 836
*t,xv, lmu,ca,ia,nin,rsq,alm,nlp,jerr)
if(jerr.gt.0) return 837
10110 do 10111 k=1,lmu 837
alm(k)=ys*alm(k) 837
nk=nin(k) 838
10120 do 10121 l=1,nk 838
ca(l,k)=ys*ca(l,k)/xs(ia(l)) 838
10121 continue 838
10122 continue 838
a0(k)=0.0 839
if(intr.ne.0) a0(k)=ym-dot_product(ca(1:nk,k),xm(ia(1:nk))) 840
10111 continue 841
10112 continue 841
deallocate(xm,xs,g,ju,xv,vlam) 842
return 843
end 844
subroutine standard (no,ni,x,y,w,isd,intr,ju,g,xm,xs,ym,ys,xv,jerr 845
*)
real x(no,ni),y(no),w(no),g(ni),xm(ni),xs(ni),xv(ni) 845
integer ju(ni) 846
real, dimension (:), allocatable :: v
allocate(v(1:no),stat=jerr) 849
if(jerr.ne.0) return 850
w=w/sum(w) 850
v=sqrt(w) 851
if(intr .ne. 0)goto 10141 851
ym=0.0 851
y=v*y 852
ys=sqrt(dot_product(y,y)-dot_product(v,y)**2) 852
y=y/ys 853
10150 do 10151 j=1,ni 853
if(ju(j).eq.0)goto 10151 853
xm(j)=0.0 853
x(:,j)=v*x(:,j) 854
xv(j)=dot_product(x(:,j),x(:,j)) 855
if(isd .eq. 0)goto 10171 855
xbq=dot_product(v,x(:,j))**2 855
vc=xv(j)-xbq 856
xs(j)=sqrt(vc) 856
x(:,j)=x(:,j)/xs(j) 856
xv(j)=1.0+xbq/vc 857
goto 10181 858
10171 continue 858
xs(j)=1.0 858
10181 continue 859
10161 continue 859
10151 continue 860
10152 continue 860
goto 10191 861
10141 continue 862
10200 do 10201 j=1,ni 862
if(ju(j).eq.0)goto 10201 863
xm(j)=dot_product(w,x(:,j)) 863
x(:,j)=v*(x(:,j)-xm(j)) 864
xv(j)=dot_product(x(:,j),x(:,j)) 864
if(isd.gt.0) xs(j)=sqrt(xv(j)) 865
10201 continue 866
10202 continue 866
if(isd .ne. 0)goto 10221 866
xs=1.0 866
goto 10231 867
10221 continue 868
10240 do 10241 j=1,ni 868
if(ju(j).eq.0)goto 10241 868
x(:,j)=x(:,j)/xs(j) 868
10241 continue 869
10242 continue 869
xv=1.0 870
10231 continue 871
10211 continue 871
ym=dot_product(w,y) 871
y=v*(y-ym) 871
ys=sqrt(dot_product(y,y)) 871
y=y/ys 872
10191 continue 873
10131 continue 873
g=0.0 873
10250 do 10251 j=1,ni 873
if(ju(j).ne.0) g(j)=dot_product(y,x(:,j)) 873
10251 continue 874
10252 continue 874
deallocate(v) 875
return 876
end 877
subroutine elnet1 (beta,ni,ju,vp,cl,g,no,ne,nx,x,nlam,flmin,ulam,t 879
*hr,maxit,xv, lmu,ao,ia,kin,rsqo,almo,nlp,jerr)
real vp(ni),g(ni),x(no,ni),ulam(nlam),ao(nx,nlam),rsqo(nlam),almo( 880
*nlam),xv(ni)
real cl(2,ni) 881
integer ju(ni),ia(nx),kin(nlam) 882
real, dimension (:), allocatable :: a,da
integer, dimension (:), allocatable :: mm
real, dimension (:,:), allocatable :: c
allocate(c(1:ni,1:nx),stat=jerr)
call get_int_parms(sml,eps,big,mnlam,rsqmax,pmin,exmx) 889
allocate(a(1:ni),stat=ierr) 889
jerr=jerr+ierr 890
allocate(mm(1:ni),stat=ierr) 890
jerr=jerr+ierr 891
allocate(da(1:ni),stat=ierr) 891
jerr=jerr+ierr 892
if(jerr.ne.0) return 893
bta=beta 893
omb=1.0-bta 894
if(flmin .ge. 1.0)goto 10271 894
eqs=max(eps,flmin) 894
alf=eqs**(1.0/(nlam-1)) 894
10271 continue 895
rsq=0.0 895
a=0.0 895
mm=0 895
nlp=0 895
nin=nlp 895
iz=0 895
mnl=min(mnlam,nlam) 896
10280 do 10281 m=1,nlam 897
if(flmin .lt. 1.0)goto 10301 897
alm=ulam(m) 897
goto 10291 898
10301 if(m .le. 2)goto 10311 898
alm=alm*alf 898
goto 10291 899
10311 if(m .ne. 1)goto 10321 899
alm=big 899
goto 10331 900
10321 continue 900
alm=0.0 901
10340 do 10341 j=1,ni 901
if(ju(j).eq.0)goto 10341 901
if(vp(j).le.0.0)goto 10341 902
alm=max(alm,abs(g(j))/vp(j)) 903
10341 continue 904
10342 continue 904
alm=alf*alm/max(bta,1.0e-3) 905
10331 continue 906
10291 continue 906
dem=alm*omb 906
ab=alm*bta 906
rsq0=rsq 906
jz=1 907
10350 continue 907
10351 continue 907
if(iz*jz.ne.0) go to 10360 907
nlp=nlp+1 907
dlx=0.0 908
10370 do 10371 k=1,ni 908
if(ju(k).eq.0)goto 10371 909
ak=a(k) 909
u=g(k)+ak*xv(k) 909
v=abs(u)-vp(k)*ab 909
a(k)=0.0 911
if(v.gt.0.0) a(k)=max(cl(1,k),min(cl(2,k),sign(v,u)/(xv(k)+vp(k)*d 912
*em)))
if(a(k).eq.ak)goto 10371 913
if(mm(k) .ne. 0)goto 10391 913
nin=nin+1 913
if(nin.gt.nx)goto 10372 914
10400 do 10401 j=1,ni 914
if(ju(j).eq.0)goto 10401 915
if(mm(j) .eq. 0)goto 10421 915
c(j,nin)=c(k,mm(j)) 915
goto 10401 915
10421 continue 916
if(j .ne. k)goto 10441 916
c(j,nin)=xv(j) 916
goto 10401 916
10441 continue 917
c(j,nin)=dot_product(x(:,j),x(:,k)) 918
10401 continue 919
10402 continue 919
mm(k)=nin 919
ia(nin)=k 920
10391 continue 921
del=a(k)-ak 921
rsq=rsq+del*(2.0*g(k)-del*xv(k)) 922
dlx=max(xv(k)*del**2,dlx) 923
10450 do 10451 j=1,ni 923
if(ju(j).ne.0) g(j)=g(j)-c(j,mm(k))*del 923
10451 continue 924
10452 continue 924
10371 continue 925
10372 continue 925
if(dlx.lt.thr)goto 10352 925
if(nin.gt.nx)goto 10352 926
if(nlp .le. maxit)goto 10471 926
jerr=-m 926
return 926
10471 continue 927
10360 continue 927
iz=1 927
da(1:nin)=a(ia(1:nin)) 928
10480 continue 928
10481 continue 928
nlp=nlp+1 928
dlx=0.0 929
10490 do 10491 l=1,nin 929
k=ia(l) 929
ak=a(k) 929
u=g(k)+ak*xv(k) 929
v=abs(u)-vp(k)*ab 930
a(k)=0.0 932
if(v.gt.0.0) a(k)=max(cl(1,k),min(cl(2,k),sign(v,u)/(xv(k)+vp(k)*d 933
*em)))
if(a(k).eq.ak)goto 10491 934
del=a(k)-ak 934
rsq=rsq+del*(2.0*g(k)-del*xv(k)) 935
dlx=max(xv(k)*del**2,dlx) 936
10500 do 10501 j=1,nin 936
g(ia(j))=g(ia(j))-c(ia(j),mm(k))*del 936
10501 continue 937
10502 continue 937
10491 continue 938
10492 continue 938
if(dlx.lt.thr)goto 10482 938
if(nlp .le. maxit)goto 10521 938
jerr=-m 938
return 938
10521 continue 939
goto 10481 940
10482 continue 940
da(1:nin)=a(ia(1:nin))-da(1:nin) 941
10530 do 10531 j=1,ni 941
if(mm(j).ne.0)goto 10531 942
if(ju(j).ne.0) g(j)=g(j)-dot_product(da(1:nin),c(j,1:nin)) 943
10531 continue 944
10532 continue 944
jz=0 945
goto 10351 946
10352 continue 946
if(nin .le. nx)goto 10551 946
jerr=-10000-m 946
goto 10282 946
10551 continue 947
if(nin.gt.0) ao(1:nin,m)=a(ia(1:nin)) 947
kin(m)=nin 948
rsqo(m)=rsq 948
almo(m)=alm 948
lmu=m 949
if(m.lt.mnl)goto 10281 949
if(flmin.ge.1.0)goto 10281 950
me=0 950
10560 do 10561 j=1,nin 950
if(ao(j,m).ne.0.0) me=me+1 950
10561 continue 950
10562 continue 950
if(me.gt.ne)goto 10282 951
if(rsq-rsq0.lt.sml*rsq)goto 10282 951
if(rsq.gt.rsqmax)goto 10282 952
10281 continue 953
10282 continue 953
deallocate(a,mm,c,da) 954
return 955
end 956
subroutine elnetn (parm,no,ni,x,y,w,jd,vp,cl,ne,nx,nlam,flmin,ulam 958
*,thr,isd, intr,maxit,lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
real vp(ni),x(no,ni),y(no),w(no),ulam(nlam),cl(2,ni) 959
real ca(nx,nlam),a0(nlam),rsq(nlam),alm(nlam) 960
integer jd(*),ia(nx),nin(nlam) 961
real, dimension (:), allocatable :: xm,xs,xv,vlam
integer, dimension (:), allocatable :: ju
allocate(xm(1:ni),stat=jerr) 966
allocate(xs(1:ni),stat=ierr) 966
jerr=jerr+ierr 967
allocate(ju(1:ni),stat=ierr) 967
jerr=jerr+ierr 968
allocate(xv(1:ni),stat=ierr) 968
jerr=jerr+ierr 969
allocate(vlam(1:nlam),stat=ierr) 969
jerr=jerr+ierr 970
if(jerr.ne.0) return 971
call chkvars(no,ni,x,ju) 972
if(jd(1).gt.0) ju(jd(2:(jd(1)+1)))=0 973
if(maxval(ju) .gt. 0)goto 10581 973
jerr=7777 973
return 973
10581 continue 974
call standard1(no,ni,x,y,w,isd,intr,ju,xm,xs,ym,ys,xv,jerr) 975
if(jerr.ne.0) return 976
cl=cl/ys 976
if(isd .le. 0)goto 10601 976
10610 do 10611 j=1,ni 976
cl(:,j)=cl(:,j)*xs(j) 976
10611 continue 976
10612 continue 976
10601 continue 977
if(flmin.ge.1.0) vlam=ulam/ys 978
call elnet2(parm,ni,ju,vp,cl,y,no,ne,nx,x,nlam,flmin,vlam,thr,maxi 980
*t,xv, lmu,ca,ia,nin,rsq,alm,nlp,jerr)
if(jerr.gt.0) return 981
10620 do 10621 k=1,lmu 981
alm(k)=ys*alm(k) 981
nk=nin(k) 982
10630 do 10631 l=1,nk 982
ca(l,k)=ys*ca(l,k)/xs(ia(l)) 982
10631 continue 982
10632 continue 982
a0(k)=0.0 983
if(intr.ne.0) a0(k)=ym-dot_product(ca(1:nk,k),xm(ia(1:nk))) 984
10621 continue 985
10622 continue 985
deallocate(xm,xs,ju,xv,vlam) 986
return 987
end 988
subroutine standard1 (no,ni,x,y,w,isd,intr,ju,xm,xs,ym,ys,xv,jerr) 989
real x(no,ni),y(no),w(no),xm(ni),xs(ni),xv(ni) 989
integer ju(ni) 990
real, dimension (:), allocatable :: v
allocate(v(1:no),stat=jerr) 993
if(jerr.ne.0) return 994
w=w/sum(w) 994
v=sqrt(w) 995
if(intr .ne. 0)goto 10651 995
ym=0.0 995
y=v*y 996
ys=sqrt(dot_product(y,y)-dot_product(v,y)**2) 996
y=y/ys 997
10660 do 10661 j=1,ni 997
if(ju(j).eq.0)goto 10661 997
xm(j)=0.0 997
x(:,j)=v*x(:,j) 998
xv(j)=dot_product(x(:,j),x(:,j)) 999
if(isd .eq. 0)goto 10681 999
xbq=dot_product(v,x(:,j))**2 999
vc=xv(j)-xbq 1000
xs(j)=sqrt(vc) 1000
x(:,j)=x(:,j)/xs(j) 1000
xv(j)=1.0+xbq/vc 1001
goto 10691 1002
10681 continue 1002
xs(j)=1.0 1002
10691 continue 1003
10671 continue 1003
10661 continue 1004
10662 continue 1004
go to 10700 1005
10651 continue 1006
10710 do 10711 j=1,ni 1006
if(ju(j).eq.0)goto 10711 1007
xm(j)=dot_product(w,x(:,j)) 1007
x(:,j)=v*(x(:,j)-xm(j)) 1008
xv(j)=dot_product(x(:,j),x(:,j)) 1008
if(isd.gt.0) xs(j)=sqrt(xv(j)) 1009
10711 continue 1010
10712 continue 1010
if(isd .ne. 0)goto 10731 1010
xs=1.0 1010
goto 10741 1011
10731 continue 1011
10750 do 10751 j=1,ni 1011
if(ju(j).eq.0)goto 10751 1011
x(:,j)=x(:,j)/xs(j) 1011
10751 continue 1012
10752 continue 1012
xv=1.0 1013
10741 continue 1014
10721 continue 1014
ym=dot_product(w,y) 1014
y=v*(y-ym) 1014
ys=sqrt(dot_product(y,y)) 1014
y=y/ys 1015
10700 continue 1015
deallocate(v) 1016
return 1017
end 1018
subroutine elnet2(beta,ni,ju,vp,cl,y,no,ne,nx,x,nlam,flmin,ulam,th 1020
*r,maxit,xv, lmu,ao,ia,kin,rsqo,almo,nlp,jerr)
real vp(ni),y(no),x(no,ni),ulam(nlam),ao(nx,nlam),rsqo(nlam),almo( 1021
*nlam),xv(ni)
real cl(2,ni) 1022
integer ju(ni),ia(nx),kin(nlam) 1023
real, dimension (:), allocatable :: a,g
integer, dimension (:), allocatable :: mm,ix
call get_int_parms(sml,eps,big,mnlam,rsqmax,pmin,exmx) 1028
allocate(a(1:ni),stat=jerr) 1029
allocate(mm(1:ni),stat=ierr) 1029
jerr=jerr+ierr 1030
allocate(g(1:ni),stat=ierr) 1030
jerr=jerr+ierr 1031
allocate(ix(1:ni),stat=ierr) 1031
jerr=jerr+ierr 1032
if(jerr.ne.0) return 1033
bta=beta 1033
omb=1.0-bta 1033
ix=0 1034
if(flmin .ge. 1.0)goto 10771 1034
eqs=max(eps,flmin) 1034
alf=eqs**(1.0/(nlam-1)) 1034
10771 continue 1035
rsq=0.0 1035
a=0.0 1035
mm=0 1035
nlp=0 1035
nin=nlp 1035
iz=0 1035
mnl=min(mnlam,nlam) 1035
alm=0.0 1036
10780 do 10781 j=1,ni 1036
if(ju(j).eq.0)goto 10781 1036
g(j)=abs(dot_product(y,x(:,j))) 1036
10781 continue 1037
10782 continue 1037
10790 do 10791 m=1,nlam 1037
alm0=alm 1038
if(flmin .lt. 1.0)goto 10811 1038
alm=ulam(m) 1038
goto 10801 1039
10811 if(m .le. 2)goto 10821 1039
alm=alm*alf 1039
goto 10801 1040
10821 if(m .ne. 1)goto 10831 1040
alm=big 1040
goto 10841 1041
10831 continue 1041
alm0=0.0 1042
10850 do 10851 j=1,ni 1042
if(ju(j).eq.0)goto 10851 1042
if(vp(j).gt.0.0) alm0=max(alm0,g(j)/vp(j)) 1042
10851 continue 1043
10852 continue 1043
alm0=alm0/max(bta,1.0e-3) 1043
alm=alf*alm0 1044
10841 continue 1045
10801 continue 1045
dem=alm*omb 1045
ab=alm*bta 1045
rsq0=rsq 1045
jz=1 1046
tlam=bta*(2.0*alm-alm0) 1047
10860 do 10861 k=1,ni 1047
if(ix(k).eq.1)goto 10861 1047
if(ju(k).eq.0)goto 10861 1048
if(g(k).gt.tlam*vp(k)) ix(k)=1 1049
10861 continue 1050
10862 continue 1050
10870 continue 1050
10871 continue 1050
if(iz*jz.ne.0) go to 10360 1051
10880 continue 1051
nlp=nlp+1 1051
dlx=0.0 1052
10890 do 10891 k=1,ni 1052
if(ix(k).eq.0)goto 10891 1052
gk=dot_product(y,x(:,k)) 1053
ak=a(k) 1053
u=gk+ak*xv(k) 1053
v=abs(u)-vp(k)*ab 1053
a(k)=0.0 1055
if(v.gt.0.0) a(k)=max(cl(1,k),min(cl(2,k),sign(v,u)/(xv(k)+vp(k)*d 1056
*em)))
if(a(k).eq.ak)goto 10891 1057
if(mm(k) .ne. 0)goto 10911 1057
nin=nin+1 1057
if(nin.gt.nx)goto 10892 1058
mm(k)=nin 1058
ia(nin)=k 1059
10911 continue 1060
del=a(k)-ak 1060
rsq=rsq+del*(2.0*gk-del*xv(k)) 1061
y=y-del*x(:,k) 1061
dlx=max(xv(k)*del**2,dlx) 1062
10891 continue 1063
10892 continue 1063
if(nin.gt.nx)goto 10872 1064
if(dlx .ge. thr)goto 10931 1064
ixx=0 1065
10940 do 10941 k=1,ni 1065
if(ix(k).eq.1)goto 10941 1065
if(ju(k).eq.0)goto 10941 1066
g(k)=abs(dot_product(y,x(:,k))) 1067
if(g(k) .le. ab*vp(k))goto 10961 1067
ix(k)=1 1067
ixx=1 1067
10961 continue 1068
10941 continue 1069
10942 continue 1069
if(ixx.eq.1) go to 10880 1070
goto 10872 1071
10931 continue 1072
if(nlp .le. maxit)goto 10981 1072
jerr=-m 1072
return 1072
10981 continue 1073
10360 continue 1073
iz=1 1074
10990 continue 1074
10991 continue 1074
nlp=nlp+1 1074
dlx=0.0 1075
11000 do 11001 l=1,nin 1075
k=ia(l) 1075
gk=dot_product(y,x(:,k)) 1076
ak=a(k) 1076
u=gk+ak*xv(k) 1076
v=abs(u)-vp(k)*ab 1076
a(k)=0.0 1078
if(v.gt.0.0) a(k)=max(cl(1,k),min(cl(2,k),sign(v,u)/(xv(k)+vp(k)*d 1079
*em)))
if(a(k).eq.ak)goto 11001 1080
del=a(k)-ak 1080
rsq=rsq+del*(2.0*gk-del*xv(k)) 1081
y=y-del*x(:,k) 1081
dlx=max(xv(k)*del**2,dlx) 1082
11001 continue 1083
11002 continue 1083
if(dlx.lt.thr)goto 10992 1083
if(nlp .le. maxit)goto 11021 1083
jerr=-m 1083
return 1083
11021 continue 1084
goto 10991 1085
10992 continue 1085
jz=0 1086
goto 10871 1087
10872 continue 1087
if(nin .le. nx)goto 11041 1087
jerr=-10000-m 1087
goto 10792 1087
11041 continue 1088
if(nin.gt.0) ao(1:nin,m)=a(ia(1:nin)) 1088
kin(m)=nin 1089
rsqo(m)=rsq 1089
almo(m)=alm 1089
lmu=m 1090
if(m.lt.mnl)goto 10791 1090
if(flmin.ge.1.0)goto 10791 1091
me=0 1091
11050 do 11051 j=1,nin 1091
if(ao(j,m).ne.0.0) me=me+1 1091
11051 continue 1091
11052 continue 1091
if(me.gt.ne)goto 10792 1092
if(rsq-rsq0.lt.sml*rsq)goto 10792 1092
if(rsq.gt.rsqmax)goto 10792 1093
10791 continue 1094
10792 continue 1094
deallocate(a,mm,g,ix) 1095
return 1096
end 1097
subroutine chkvars(no,ni,x,ju) 1098
real x(no,ni) 1098
integer ju(ni) 1099
11060 do 11061 j=1,ni 1099
ju(j)=0 1099
t=x(1,j) 1100
11070 do 11071 i=2,no 1100
if(x(i,j).eq.t)goto 11071 1100
ju(j)=1 1100
goto 11072 1100
11071 continue 1101
11072 continue 1101
11061 continue 1102
11062 continue 1102
return 1103
end 1104
subroutine uncomp(ni,ca,ia,nin,a) 1105
real ca(*),a(ni) 1105
integer ia(*) 1106
a=0.0 1106
if(nin.gt.0) a(ia(1:nin))=ca(1:nin) 1107
return 1108
end 1109
subroutine modval(a0,ca,ia,nin,n,x,f) 1110
real ca(nin),x(n,*),f(n) 1110
integer ia(nin) 1111
f=a0 1111
if(nin.le.0) return 1112
11080 do 11081 i=1,n 1112
f(i)=f(i)+dot_product(ca(1:nin),x(i,ia(1:nin))) 1112
11081 continue 1113
11082 continue 1113
return 1114
end 1115
subroutine spelnet (ka,parm,no,ni,x,ix,jx,y,w,jd,vp,cl,ne,nx,nlam 1118
*,flmin,ulam,thr,isd,intr, maxit,lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr
*)
real x(*),y(no),w(no),vp(ni),ulam(nlam),cl(2,ni) 1119
real ca(nx,nlam),a0(nlam),rsq(nlam),alm(nlam) 1120
integer ix(*),jx(*),jd(*),ia(nx),nin(nlam) 1121
real, dimension (:), allocatable :: vq;
if(maxval(vp) .gt. 0.0)goto 11101 1124
jerr=10000 1124
return 1124
11101 continue 1125
allocate(vq(1:ni),stat=jerr) 1125
if(jerr.ne.0) return 1126
vq=max(0.0,vp) 1126
vq=vq*ni/sum(vq) 1127
if(ka .ne. 1)goto 11121 1128
call spelnetu (parm,no,ni,x,ix,jx,y,w,jd,vq,cl,ne,nx,nlam,flmin,u 1131
*lam,thr,isd, intr,maxit,lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
goto 11131 1132
11121 continue 1133
call spelnetn (parm,no,ni,x,ix,jx,y,w,jd,vq,cl,ne,nx,nlam,flmin,ul 1136
*am,thr,isd,intr, maxit,lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
11131 continue 1137
11111 continue 1137
deallocate(vq) 1138
return 1139
end 1140
subroutine spelnetu (parm,no,ni,x,ix,jx,y,w,jd,vp,cl,ne,nx,nlam,f 1143
*lmin,ulam,thr,isd,intr, maxit,lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
real x(*),y(no),w(no),vp(ni),ulam(nlam),cl(2,ni) 1144
real ca(nx,nlam),a0(nlam),rsq(nlam),alm(nlam) 1145
integer ix(*),jx(*),jd(*),ia(nx),nin(nlam) 1146
real, dimension (:), allocatable :: xm,xs,g,xv,vlam
integer, dimension (:), allocatable :: ju
allocate(g(1:ni),stat=jerr) 1151
allocate(xm(1:ni),stat=ierr) 1151
jerr=jerr+ierr 1152
allocate(xs(1:ni),stat=ierr) 1152
jerr=jerr+ierr 1153
allocate(ju(1:ni),stat=ierr) 1153
jerr=jerr+ierr 1154
allocate(xv(1:ni),stat=ierr) 1154
jerr=jerr+ierr 1155
allocate(vlam(1:nlam),stat=ierr) 1155
jerr=jerr+ierr 1156
if(jerr.ne.0) return 1157
call spchkvars(no,ni,x,ix,ju) 1158
if(jd(1).gt.0) ju(jd(2:(jd(1)+1)))=0 1159
if(maxval(ju) .gt. 0)goto 11151 1159
jerr=7777 1159
return 1159
11151 continue 1160
call spstandard(no,ni,x,ix,jx,y,w,ju,isd,intr,g,xm,xs,ym,ys,xv,jer 1161
*r)
if(jerr.ne.0) return 1162
cl=cl/ys 1162
if(isd .le. 0)goto 11171 1162
11180 do 11181 j=1,ni 1162
cl(:,j)=cl(:,j)*xs(j) 1162
11181 continue 1162
11182 continue 1162
11171 continue 1163
if(flmin.ge.1.0) vlam=ulam/ys 1164
call spelnet1(parm,ni,g,no,w,ne,nx,x,ix,jx,ju,vp,cl,nlam,flmin,vla 1166
*m,thr,maxit, xm,xs,xv,lmu,ca,ia,nin,rsq,alm,nlp,jerr)
if(jerr.gt.0) return 1167
11190 do 11191 k=1,lmu 1167
alm(k)=ys*alm(k) 1167
nk=nin(k) 1168
11200 do 11201 l=1,nk 1168
ca(l,k)=ys*ca(l,k)/xs(ia(l)) 1168
11201 continue 1168
11202 continue 1168
a0(k)=0.0 1169
if(intr.ne.0) a0(k)=ym-dot_product(ca(1:nk,k),xm(ia(1:nk))) 1170
11191 continue 1171
11192 continue 1171
deallocate(xm,xs,g,ju,xv,vlam) 1172
return 1173
end 1174
subroutine spstandard (no,ni,x,ix,jx,y,w,ju,isd,intr,g,xm,xs,ym,ys 1175
*,xv,jerr)
real x(*),y(no),w(no),g(ni),xm(ni),xs(ni),xv(ni) 1175
integer ix(*),jx(*),ju(ni) 1176
w=w/sum(w) 1177
if(intr .ne. 0)goto 11221 1177
ym=0.0 1178
ys=sqrt(dot_product(w,y**2)-dot_product(w,y)**2) 1178
y=y/ys 1179
11230 do 11231 j=1,ni 1179
if(ju(j).eq.0)goto 11231 1179
xm(j)=0.0 1179
jb=ix(j) 1179
je=ix(j+1)-1 1180
xv(j)=dot_product(w(jx(jb:je)),x(jb:je)**2) 1181
if(isd .eq. 0)goto 11251 1181
xbq=dot_product(w(jx(jb:je)),x(jb:je))**2 1181
vc=xv(j)-xbq 1182
xs(j)=sqrt(vc) 1182
xv(j)=1.0+xbq/vc 1183
goto 11261 1184
11251 continue 1184
xs(j)=1.0 1184
11261 continue 1185
11241 continue 1185
11231 continue 1186
11232 continue 1186
goto 11271 1187
11221 continue 1188
11280 do 11281 j=1,ni 1188
if(ju(j).eq.0)goto 11281 1189
jb=ix(j) 1189
je=ix(j+1)-1 1189
xm(j)=dot_product(w(jx(jb:je)),x(jb:je)) 1190
xv(j)=dot_product(w(jx(jb:je)),x(jb:je)**2)-xm(j)**2 1191
if(isd.gt.0) xs(j)=sqrt(xv(j)) 1192
11281 continue 1193
11282 continue 1193
if(isd .ne. 0)goto 11301 1193
xs=1.0 1193
goto 11311 1193
11301 continue 1193
xv=1.0 1193
11311 continue 1194
11291 continue 1194
ym=dot_product(w,y) 1194
y=y-ym 1194
ys=sqrt(dot_product(w,y**2)) 1194
y=y/ys 1195
11271 continue 1196
11211 continue 1196
g=0.0 1197
11320 do 11321 j=1,ni 1197
if(ju(j).eq.0)goto 11321 1197
jb=ix(j) 1197
je=ix(j+1)-1 1198
g(j)=dot_product(w(jx(jb:je))*y(jx(jb:je)),x(jb:je))/xs(j) 1199
11321 continue 1200
11322 continue 1200
return 1201
end 1202
subroutine spelnet1(beta,ni,g,no,w,ne,nx,x,ix,jx,ju,vp,cl,nlam,flm 1204
*in,ulam, thr,maxit,xm,xs,xv,lmu,ao,ia,kin,rsqo,almo,nlp,jerr)
real g(ni),vp(ni),x(*),ulam(nlam),w(no) 1205
real ao(nx,nlam),rsqo(nlam),almo(nlam),xm(ni),xs(ni),xv(ni),cl(2,n 1206
*i)
integer ix(*),jx(*),ju(ni),ia(nx),kin(nlam) 1207
real, dimension (:), allocatable :: a,da
integer, dimension (:), allocatable :: mm
real, dimension (:,:), allocatable :: c
allocate(c(1:ni,1:nx),stat=jerr)
call get_int_parms(sml,eps,big,mnlam,rsqmax,pmin,exmx) 1214
allocate(a(1:ni),stat=ierr) 1214
jerr=jerr+ierr 1215
allocate(mm(1:ni),stat=ierr) 1215
jerr=jerr+ierr 1216
allocate(da(1:ni),stat=ierr) 1216
jerr=jerr+ierr 1217
if(jerr.ne.0) return 1218
bta=beta 1218
omb=1.0-bta 1219
if(flmin .ge. 1.0)goto 11341 1219
eqs=max(eps,flmin) 1219
alf=eqs**(1.0/(nlam-1)) 1219
11341 continue 1220
rsq=0.0 1220
a=0.0 1220
mm=0 1220
nlp=0 1220
nin=nlp 1220
iz=0 1220
mnl=min(mnlam,nlam) 1221
11350 do 11351 m=1,nlam 1222
if(flmin .lt. 1.0)goto 11371 1222
alm=ulam(m) 1222
goto 11361 1223
11371 if(m .le. 2)goto 11381 1223
alm=alm*alf 1223
goto 11361 1224
11381 if(m .ne. 1)goto 11391 1224
alm=big 1224
goto 11401 1225
11391 continue 1225
alm=0.0 1226
11410 do 11411 j=1,ni 1226
if(ju(j).eq.0)goto 11411 1226
if(vp(j).le.0.0)goto 11411 1227
alm=max(alm,abs(g(j))/vp(j)) 1228
11411 continue 1229
11412 continue 1229
alm=alf*alm/max(bta,1.0e-3) 1230
11401 continue 1231
11361 continue 1231
dem=alm*omb 1231
ab=alm*bta 1231
rsq0=rsq 1231
jz=1 1232
11420 continue 1232
11421 continue 1232
if(iz*jz.ne.0) go to 10360 1232
nlp=nlp+1 1232
dlx=0.0 1233
11430 do 11431 k=1,ni 1233
if(ju(k).eq.0)goto 11431 1234
ak=a(k) 1234
u=g(k)+ak*xv(k) 1234
v=abs(u)-vp(k)*ab 1234
a(k)=0.0 1236
if(v.gt.0.0) a(k)=max(cl(1,k),min(cl(2,k),sign(v,u)/(xv(k)+vp(k)*d 1237
*em)))
if(a(k).eq.ak)goto 11431 1238
if(mm(k) .ne. 0)goto 11451 1238
nin=nin+1 1238
if(nin.gt.nx)goto 11432 1239
11460 do 11461 j=1,ni 1239
if(ju(j).eq.0)goto 11461 1240
if(mm(j) .eq. 0)goto 11481 1240
c(j,nin)=c(k,mm(j)) 1240
goto 11461 1240
11481 continue 1241
if(j .ne. k)goto 11501 1241
c(j,nin)=xv(j) 1241
goto 11461 1241
11501 continue 1242
c(j,nin)= (row_prod(j,k,ix,jx,x,w)-xm(j)*xm(k))/(xs(j)*xs(k)) 1244
11461 continue 1245
11462 continue 1245
mm(k)=nin 1245
ia(nin)=k 1246
11451 continue 1247
del=a(k)-ak 1247
rsq=rsq+del*(2.0*g(k)-del*xv(k)) 1248
dlx=max(xv(k)*del**2,dlx) 1249
11510 do 11511 j=1,ni 1249
if(ju(j).ne.0) g(j)=g(j)-c(j,mm(k))*del 1249
11511 continue 1250
11512 continue 1250
11431 continue 1251
11432 continue 1251
if(dlx.lt.thr)goto 11422 1251
if(nin.gt.nx)goto 11422 1252
if(nlp .le. maxit)goto 11531 1252
jerr=-m 1252
return 1252
11531 continue 1253
10360 continue 1253
iz=1 1253
da(1:nin)=a(ia(1:nin)) 1254
11540 continue 1254
11541 continue 1254
nlp=nlp+1 1254
dlx=0.0 1255
11550 do 11551 l=1,nin 1255
k=ia(l) 1256
ak=a(k) 1256
u=g(k)+ak*xv(k) 1256
v=abs(u)-vp(k)*ab 1256
a(k)=0.0 1258
if(v.gt.0.0) a(k)=max(cl(1,k),min(cl(2,k),sign(v,u)/(xv(k)+vp(k)*d 1259
*em)))
if(a(k).eq.ak)goto 11551 1260
del=a(k)-ak 1260
rsq=rsq+del*(2.0*g(k)-del*xv(k)) 1261
dlx=max(xv(k)*del**2,dlx) 1262
11560 do 11561 j=1,nin 1262
g(ia(j))=g(ia(j))-c(ia(j),mm(k))*del 1262
11561 continue 1263
11562 continue 1263
11551 continue 1264
11552 continue 1264
if(dlx.lt.thr)goto 11542 1264
if(nlp .le. maxit)goto 11581 1264
jerr=-m 1264
return 1264
11581 continue 1265
goto 11541 1266
11542 continue 1266
da(1:nin)=a(ia(1:nin))-da(1:nin) 1267
11590 do 11591 j=1,ni 1267
if(mm(j).ne.0)goto 11591 1268
if(ju(j).ne.0) g(j)=g(j)-dot_product(da(1:nin),c(j,1:nin)) 1269
11591 continue 1270
11592 continue 1270
jz=0 1271
goto 11421 1272
11422 continue 1272
if(nin .le. nx)goto 11611 1272
jerr=-10000-m 1272
goto 11352 1272
11611 continue 1273
if(nin.gt.0) ao(1:nin,m)=a(ia(1:nin)) 1273
kin(m)=nin 1274
rsqo(m)=rsq 1274
almo(m)=alm 1274
lmu=m 1275
if(m.lt.mnl)goto 11351 1275
if(flmin.ge.1.0)goto 11351 1276
me=0 1276
11620 do 11621 j=1,nin 1276
if(ao(j,m).ne.0.0) me=me+1 1276
11621 continue 1276
11622 continue 1276
if(me.gt.ne)goto 11352 1277
if(rsq-rsq0.lt.sml*rsq)goto 11352 1277
if(rsq.gt.rsqmax)goto 11352 1278
11351 continue 1279
11352 continue 1279
deallocate(a,mm,c,da) 1280
return 1281
end 1282
subroutine spelnetn(parm,no,ni,x,ix,jx,y,w,jd,vp,cl,ne,nx,nlam,flm 1284
*in,ulam, thr,isd,intr,maxit,lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
real x(*),vp(ni),y(no),w(no),ulam(nlam),cl(2,ni) 1285
real ca(nx,nlam),a0(nlam),rsq(nlam),alm(nlam) 1286
integer ix(*),jx(*),jd(*),ia(nx),nin(nlam) 1287
real, dimension (:), allocatable :: xm,xs,xv,vlam
integer, dimension (:), allocatable :: ju
allocate(xm(1:ni),stat=jerr) 1292
allocate(xs(1:ni),stat=ierr) 1292
jerr=jerr+ierr 1293
allocate(ju(1:ni),stat=ierr) 1293
jerr=jerr+ierr 1294
allocate(xv(1:ni),stat=ierr) 1294
jerr=jerr+ierr 1295
allocate(vlam(1:nlam),stat=ierr) 1295
jerr=jerr+ierr 1296
if(jerr.ne.0) return 1297
call spchkvars(no,ni,x,ix,ju) 1298
if(jd(1).gt.0) ju(jd(2:(jd(1)+1)))=0 1299
if(maxval(ju) .gt. 0)goto 11641 1299
jerr=7777 1299
return 1299
11641 continue 1300
call spstandard1(no,ni,x,ix,jx,y,w,ju,isd,intr,xm,xs,ym,ys,xv,jerr 1301
*)
if(jerr.ne.0) return 1302
cl=cl/ys 1302
if(isd .le. 0)goto 11661 1302
11670 do 11671 j=1,ni 1302
cl(:,j)=cl(:,j)*xs(j) 1302
11671 continue 1302
11672 continue 1302
11661 continue 1303
if(flmin.ge.1.0) vlam=ulam/ys 1304
call spelnet2(parm,ni,y,w,no,ne,nx,x,ix,jx,ju,vp,cl,nlam,flmin,vla 1306
*m,thr,maxit, xm,xs,xv,lmu,ca,ia,nin,rsq,alm,nlp,jerr)
if(jerr.gt.0) return 1307
11680 do 11681 k=1,lmu 1307
alm(k)=ys*alm(k) 1307
nk=nin(k) 1308
11690 do 11691 l=1,nk 1308
ca(l,k)=ys*ca(l,k)/xs(ia(l)) 1308
11691 continue 1308
11692 continue 1308
a0(k)=0.0 1309
if(intr.ne.0) a0(k)=ym-dot_product(ca(1:nk,k),xm(ia(1:nk))) 1310
11681 continue 1311
11682 continue 1311
deallocate(xm,xs,ju,xv,vlam) 1312
return 1313
end 1314
subroutine spstandard1 (no,ni,x,ix,jx,y,w,ju,isd,intr,xm,xs,ym,ys, 1315
*xv,jerr)
real x(*),y(no),w(no),xm(ni),xs(ni),xv(ni) 1315
integer ix(*),jx(*),ju(ni) 1316
w=w/sum(w) 1317
if(intr .ne. 0)goto 11711 1317
ym=0.0 1318
ys=sqrt(dot_product(w,y**2)-dot_product(w,y)**2) 1318
y=y/ys 1319
11720 do 11721 j=1,ni 1319
if(ju(j).eq.0)goto 11721 1319
xm(j)=0.0 1319
jb=ix(j) 1319
je=ix(j+1)-1 1320
xv(j)=dot_product(w(jx(jb:je)),x(jb:je)**2) 1321
if(isd .eq. 0)goto 11741 1321
xbq=dot_product(w(jx(jb:je)),x(jb:je))**2 1321
vc=xv(j)-xbq 1322
xs(j)=sqrt(vc) 1322
xv(j)=1.0+xbq/vc 1323
goto 11751 1324
11741 continue 1324
xs(j)=1.0 1324
11751 continue 1325
11731 continue 1325
11721 continue 1326
11722 continue 1326
return 1327
11711 continue 1328
11760 do 11761 j=1,ni 1328
if(ju(j).eq.0)goto 11761 1329
jb=ix(j) 1329
je=ix(j+1)-1 1329
xm(j)=dot_product(w(jx(jb:je)),x(jb:je)) 1330
xv(j)=dot_product(w(jx(jb:je)),x(jb:je)**2)-xm(j)**2 1331
if(isd.gt.0) xs(j)=sqrt(xv(j)) 1332
11761 continue 1333
11762 continue 1333
if(isd .ne. 0)goto 11781 1333
xs=1.0 1333
goto 11791 1333
11781 continue 1333
xv=1.0 1333
11791 continue 1334
11771 continue 1334
ym=dot_product(w,y) 1334
y=y-ym 1334
ys=sqrt(dot_product(w,y**2)) 1334
y=y/ys 1335
return 1336
end 1337
subroutine spelnet2(beta,ni,y,w,no,ne,nx,x,ix,jx,ju,vp,cl,nlam,flm 1339
*in,ulam, thr,maxit,xm,xs,xv,lmu,ao,ia,kin,rsqo,almo,nlp,jerr)
real y(no),w(no),x(*),vp(ni),ulam(nlam),cl(2,ni) 1340
real ao(nx,nlam),rsqo(nlam),almo(nlam),xm(ni),xs(ni),xv(ni) 1341
integer ix(*),jx(*),ju(ni),ia(nx),kin(nlam) 1342
real, dimension (:), allocatable :: a,g
integer, dimension (:), allocatable :: mm,iy
call get_int_parms(sml,eps,big,mnlam,rsqmax,pmin,exmx) 1347
allocate(a(1:ni),stat=jerr) 1348
allocate(mm(1:ni),stat=ierr) 1348
jerr=jerr+ierr 1349
allocate(g(1:ni),stat=ierr) 1349
jerr=jerr+ierr 1350
allocate(iy(1:ni),stat=ierr) 1350
jerr=jerr+ierr 1351
if(jerr.ne.0) return 1352
bta=beta 1352
omb=1.0-bta 1352
alm=0.0 1352
iy=0 1353
if(flmin .ge. 1.0)goto 11811 1353
eqs=max(eps,flmin) 1353
alf=eqs**(1.0/(nlam-1)) 1353
11811 continue 1354
rsq=0.0 1354
a=0.0 1354
mm=0 1354
o=0.0 1354
nlp=0 1354
nin=nlp 1354
iz=0 1354
mnl=min(mnlam,nlam) 1355
11820 do 11821 j=1,ni 1355
if(ju(j).eq.0)goto 11821 1356
jb=ix(j) 1356
je=ix(j+1)-1 1357
g(j)=abs(dot_product(y(jx(jb:je))+o,w(jx(jb:je))*x(jb:je))/xs(j)) 1358
11821 continue 1359
11822 continue 1359
11830 do 11831 m=1,nlam 1359
alm0=alm 1360
if(flmin .lt. 1.0)goto 11851 1360
alm=ulam(m) 1360
goto 11841 1361
11851 if(m .le. 2)goto 11861 1361
alm=alm*alf 1361
goto 11841 1362
11861 if(m .ne. 1)goto 11871 1362
alm=big 1362
goto 11881 1363
11871 continue 1363
alm0=0.0 1364
11890 do 11891 j=1,ni 1364
if(ju(j).eq.0)goto 11891 1364
if(vp(j).gt.0.0) alm0=max(alm0,g(j)/vp(j)) 1364
11891 continue 1365
11892 continue 1365
alm0=alm0/max(bta,1.0e-3) 1365
alm=alf*alm0 1366
11881 continue 1367
11841 continue 1367
dem=alm*omb 1367
ab=alm*bta 1367
rsq0=rsq 1367
jz=1 1368
tlam=bta*(2.0*alm-alm0) 1369
11900 do 11901 k=1,ni 1369
if(iy(k).eq.1)goto 11901 1369
if(ju(k).eq.0)goto 11901 1370
if(g(k).gt.tlam*vp(k)) iy(k)=1 1371
11901 continue 1372
11902 continue 1372
11910 continue 1372
11911 continue 1372
if(iz*jz.ne.0) go to 10360 1373
10880 continue 1373
nlp=nlp+1 1373
dlx=0.0 1374
11920 do 11921 k=1,ni 1374
if(iy(k).eq.0)goto 11921 1374
jb=ix(k) 1374
je=ix(k+1)-1 1375
gk=dot_product(y(jx(jb:je))+o,w(jx(jb:je))*x(jb:je))/xs(k) 1376
ak=a(k) 1376
u=gk+ak*xv(k) 1376
v=abs(u)-vp(k)*ab 1376
a(k)=0.0 1378
if(v.gt.0.0) a(k)=max(cl(1,k),min(cl(2,k),sign(v,u)/(xv(k)+vp(k)*d 1379
*em)))
if(a(k).eq.ak)goto 11921 1380
if(mm(k) .ne. 0)goto 11941 1380
nin=nin+1 1380
if(nin.gt.nx)goto 11922 1381
mm(k)=nin 1381
ia(nin)=k 1382
11941 continue 1383
del=a(k)-ak 1383
rsq=rsq+del*(2.0*gk-del*xv(k)) 1384
y(jx(jb:je))=y(jx(jb:je))-del*x(jb:je)/xs(k) 1385
o=o+del*xm(k)/xs(k) 1385
dlx=max(xv(k)*del**2,dlx) 1386
11921 continue 1387
11922 continue 1387
if(nin.gt.nx)goto 11912 1388
if(dlx .ge. thr)goto 11961 1388
ixx=0 1389
11970 do 11971 j=1,ni 1389
if(iy(j).eq.1)goto 11971 1389
if(ju(j).eq.0)goto 11971 1390
jb=ix(j) 1390
je=ix(j+1)-1 1391
g(j)=abs(dot_product(y(jx(jb:je))+o,w(jx(jb:je))*x(jb:je))/xs(j)) 1392
if(g(j) .le. ab*vp(j))goto 11991 1392
iy(j)=1 1392
ixx=1 1392
11991 continue 1393
11971 continue 1394
11972 continue 1394
if(ixx.eq.1) go to 10880 1395
goto 11912 1396
11961 continue 1397
if(nlp .le. maxit)goto 12011 1397
jerr=-m 1397
return 1397
12011 continue 1398
10360 continue 1398
iz=1 1399
12020 continue 1399
12021 continue 1399
nlp=nlp+1 1399
dlx=0.0 1400
12030 do 12031 l=1,nin 1400
k=ia(l) 1400
jb=ix(k) 1400
je=ix(k+1)-1 1401
gk=dot_product(y(jx(jb:je))+o,w(jx(jb:je))*x(jb:je))/xs(k) 1402
ak=a(k) 1402
u=gk+ak*xv(k) 1402
v=abs(u)-vp(k)*ab 1402
a(k)=0.0 1404
if(v.gt.0.0) a(k)=max(cl(1,k),min(cl(2,k),sign(v,u)/(xv(k)+vp(k)*d 1405
*em)))
if(a(k).eq.ak)goto 12031 1406
del=a(k)-ak 1406
rsq=rsq+del*(2.0*gk-del*xv(k)) 1407
y(jx(jb:je))=y(jx(jb:je))-del*x(jb:je)/xs(k) 1408
o=o+del*xm(k)/xs(k) 1408
dlx=max(xv(k)*del**2,dlx) 1409
12031 continue 1410
12032 continue 1410
if(dlx.lt.thr)goto 12022 1410
if(nlp .le. maxit)goto 12051 1410
jerr=-m 1410
return 1410
12051 continue 1411
goto 12021 1412
12022 continue 1412
jz=0 1413
goto 11911 1414
11912 continue 1414
if(nin .le. nx)goto 12071 1414
jerr=-10000-m 1414
goto 11832 1414
12071 continue 1415
if(nin.gt.0) ao(1:nin,m)=a(ia(1:nin)) 1415
kin(m)=nin 1416
rsqo(m)=rsq 1416
almo(m)=alm 1416
lmu=m 1417
if(m.lt.mnl)goto 11831 1417
if(flmin.ge.1.0)goto 11831 1418
me=0 1418
12080 do 12081 j=1,nin 1418
if(ao(j,m).ne.0.0) me=me+1 1418
12081 continue 1418
12082 continue 1418
if(me.gt.ne)goto 11832 1419
if(rsq-rsq0.lt.sml*rsq)goto 11832 1419
if(rsq.gt.rsqmax)goto 11832 1420
11831 continue 1421
11832 continue 1421
deallocate(a,mm,g,iy) 1422
return 1423
end 1424
subroutine spchkvars(no,ni,x,ix,ju) 1425
real x(*) 1425
integer ix(*),ju(ni) 1426
12090 do 12091 j=1,ni 1426
ju(j)=0 1426
jb=ix(j) 1426
nj=ix(j+1)-jb 1426
if(nj.eq.0)goto 12091 1427
je=ix(j+1)-1 1428
if(nj .ge. no)goto 12111 1428
12120 do 12121 i=jb,je 1428
if(x(i).eq.0.0)goto 12121 1428
ju(j)=1 1428
goto 12122 1428
12121 continue 1428
12122 continue 1428
goto 12131 1429
12111 continue 1429
t=x(jb) 1429
12140 do 12141 i=jb+1,je 1429
if(x(i).eq.t)goto 12141 1429
ju(j)=1 1429
goto 12142 1429
12141 continue 1429
12142 continue 1429
12131 continue 1430
12101 continue 1430
12091 continue 1431
12092 continue 1431
return 1432
end 1433
subroutine cmodval(a0,ca,ia,nin,x,ix,jx,n,f) 1434
real ca(*),x(*),f(n) 1434
integer ia(*),ix(*),jx(*) 1435
f=a0 1436
12150 do 12151 j=1,nin 1436
k=ia(j) 1436
kb=ix(k) 1436
ke=ix(k+1)-1 1437
f(jx(kb:ke))=f(jx(kb:ke))+ca(j)*x(kb:ke) 1438
12151 continue 1439
12152 continue 1439
return 1440
end 1441
function row_prod(i,j,ia,ja,ra,w) 1442
integer ia(*),ja(*) 1442
real ra(*),w(*) 1443
row_prod=dot(ra(ia(i)),ra(ia(j)),ja(ia(i)),ja(ia(j)), ia(i+1)-ia( 1445
*i),ia(j+1)-ia(j),w)
return 1446
end 1447
function dot(x,y,mx,my,nx,ny,w) 1448
real x(*),y(*),w(*) 1448
integer mx(*),my(*) 1449
i=1 1449
j=i 1449
s=0.0 1450
12160 continue 1450
12161 continue 1450
12170 continue 1451
12171 if(mx(i).ge.my(j))goto 12172 1451
i=i+1 1451
if(i.gt.nx) go to 12180 1451
goto 12171 1452
12172 continue 1452
if(mx(i).eq.my(j)) go to 12190 1453
12200 continue 1453
12201 if(my(j).ge.mx(i))goto 12202 1453
j=j+1 1453
if(j.gt.ny) go to 12180 1453
goto 12201 1454
12202 continue 1454
if(mx(i).eq.my(j)) go to 12190 1454
goto 12161 1455
12190 continue 1455
s=s+w(mx(i))*x(i)*y(j) 1456
i=i+1 1456
if(i.gt.nx)goto 12162 1456
j=j+1 1456
if(j.gt.ny)goto 12162 1457
goto 12161 1458
12162 continue 1458
12180 continue 1458
dot=s 1459
return 1460
end 1461
subroutine lognet (parm,no,ni,nc,x,y,g,jd,vp,cl,ne,nx,nlam,flmin,u 1463
*lam,thr, isd,intr,maxit,kopt,lmu,a0,ca,ia,nin,dev0,dev,alm,nlp,je
*rr)
real x(no,ni),y(no,max(2,nc)),g(no,nc),vp(ni),ulam(nlam) 1464
real ca(nx,nc,nlam),a0(nc,nlam),dev(nlam),alm(nlam),cl(2,ni) 1465
integer jd(*),ia(nx),nin(nlam) 1466
real, dimension (:), allocatable :: xm,xs,ww,vq,xv
integer, dimension (:), allocatable :: ju
if(maxval(vp) .gt. 0.0)goto 12221 1470
jerr=10000 1470
return 1470
12221 continue 1471
allocate(ww(1:no),stat=jerr) 1472
allocate(ju(1:ni),stat=ierr) 1472
jerr=jerr+ierr 1473
allocate(vq(1:ni),stat=ierr) 1473
jerr=jerr+ierr 1474
allocate(xm(1:ni),stat=ierr) 1474
jerr=jerr+ierr 1475
if(kopt .ne. 2)goto 12241 1475
allocate(xv(1:ni),stat=ierr) 1475
jerr=jerr+ierr 1475
12241 continue 1476
if(isd .le. 0)goto 12261 1476
allocate(xs(1:ni),stat=ierr) 1476
jerr=jerr+ierr 1476
12261 continue 1477
if(jerr.ne.0) return 1478
call chkvars(no,ni,x,ju) 1479
if(jd(1).gt.0) ju(jd(2:(jd(1)+1)))=0 1480
if(maxval(ju) .gt. 0)goto 12281 1480
jerr=7777 1480
return 1480
12281 continue 1481
vq=max(0.0,vp) 1481
vq=vq*ni/sum(vq) 1482
12290 do 12291 i=1,no 1482
ww(i)=sum(y(i,:)) 1482
if(ww(i).gt.0.0) y(i,:)=y(i,:)/ww(i) 1482
12291 continue 1483
12292 continue 1483
sw=sum(ww) 1483
ww=ww/sw 1484
if(nc .ne. 1)goto 12311 1484
call lstandard1(no,ni,x,ww,ju,isd,intr,xm,xs) 1485
if(isd .le. 0)goto 12331 1485
12340 do 12341 j=1,ni 1485
cl(:,j)=cl(:,j)*xs(j) 1485
12341 continue 1485
12342 continue 1485
12331 continue 1486
call lognet2n(parm,no,ni,x,y(:,1),g(:,1),ww,ju,vq,cl,ne,nx,nlam,fl 1488
*min,ulam, thr,isd,intr,maxit,kopt,lmu,a0,ca,ia,nin,dev0,dev,alm,n
*lp,jerr)
goto 12301 1489
12311 if(kopt .ne. 2)goto 12351 1489
call multlstandard1(no,ni,x,ww,ju,isd,intr,xm,xs,xv) 1490
if(isd .le. 0)goto 12371 1490
12380 do 12381 j=1,ni 1490
cl(:,j)=cl(:,j)*xs(j) 1490
12381 continue 1490
12382 continue 1490
12371 continue 1491
call multlognetn(parm,no,ni,nc,x,y,g,ww,ju,vq,cl,ne,nx,nlam,flmin, 1493
*ulam,thr, intr,maxit,xv,lmu,a0,ca,ia,nin,dev0,dev,alm,nlp,jerr)
goto 12391 1494
12351 continue 1494
call lstandard1(no,ni,x,ww,ju,isd,intr,xm,xs) 1495
if(isd .le. 0)goto 12411 1495
12420 do 12421 j=1,ni 1495
cl(:,j)=cl(:,j)*xs(j) 1495
12421 continue 1495
12422 continue 1495
12411 continue 1496
call lognetn(parm,no,ni,nc,x,y,g,ww,ju,vq,cl,ne,nx,nlam,flmin,ulam 1498
*,thr, isd,intr,maxit,kopt,lmu,a0,ca,ia,nin,dev0,dev,alm,nlp,jerr)
12391 continue 1499
12301 continue 1499
if(jerr.gt.0) return 1499
dev0=2.0*sw*dev0 1500
12430 do 12431 k=1,lmu 1500
nk=nin(k) 1501
12440 do 12441 ic=1,nc 1501
if(isd .le. 0)goto 12461 1501
12470 do 12471 l=1,nk 1501
ca(l,ic,k)=ca(l,ic,k)/xs(ia(l)) 1501
12471 continue 1501
12472 continue 1501
12461 continue 1502
if(intr .ne. 0)goto 12491 1502
a0(ic,k)=0.0 1502
goto 12501 1503
12491 continue 1503
a0(ic,k)=a0(ic,k)-dot_product(ca(1:nk,ic,k),xm(ia(1:nk))) 1503
12501 continue 1504
12481 continue 1504
12441 continue 1505
12442 continue 1505
12431 continue 1506
12432 continue 1506
deallocate(ww,ju,vq,xm) 1506
if(isd.gt.0) deallocate(xs) 1507
if(kopt.eq.2) deallocate(xv) 1508
return 1509
end 1510
subroutine lstandard1 (no,ni,x,w,ju,isd,intr,xm,xs) 1511
real x(no,ni),w(no),xm(ni),xs(ni) 1511
integer ju(ni) 1512
if(intr .ne. 0)goto 12521 1513
12530 do 12531 j=1,ni 1513
if(ju(j).eq.0)goto 12531 1513
xm(j)=0.0 1514
if(isd .eq. 0)goto 12551 1514
vc=dot_product(w,x(:,j)**2)-dot_product(w,x(:,j))**2 1515
xs(j)=sqrt(vc) 1515
x(:,j)=x(:,j)/xs(j) 1516
12551 continue 1517
12531 continue 1518
12532 continue 1518
return 1519
12521 continue 1520
12560 do 12561 j=1,ni 1520
if(ju(j).eq.0)goto 12561 1521
xm(j)=dot_product(w,x(:,j)) 1521
x(:,j)=x(:,j)-xm(j) 1522
if(isd .le. 0)goto 12581 1522
xs(j)=sqrt(dot_product(w,x(:,j)**2)) 1522
x(:,j)=x(:,j)/xs(j) 1522
12581 continue 1523
12561 continue 1524
12562 continue 1524
return 1525
end 1526
subroutine multlstandard1 (no,ni,x,w,ju,isd,intr,xm,xs,xv) 1527
real x(no,ni),w(no),xm(ni),xs(ni),xv(ni) 1527
integer ju(ni) 1528
if(intr .ne. 0)goto 12601 1529
12610 do 12611 j=1,ni 1529
if(ju(j).eq.0)goto 12611 1529
xm(j)=0.0 1530
xv(j)=dot_product(w,x(:,j)**2) 1531
if(isd .eq. 0)goto 12631 1531
xbq=dot_product(w,x(:,j))**2 1531
vc=xv(j)-xbq 1532
xs(j)=sqrt(vc) 1532
x(:,j)=x(:,j)/xs(j) 1532
xv(j)=1.0+xbq/vc 1533
12631 continue 1534
12611 continue 1535
12612 continue 1535
return 1536
12601 continue 1537
12640 do 12641 j=1,ni 1537
if(ju(j).eq.0)goto 12641 1538
xm(j)=dot_product(w,x(:,j)) 1538
x(:,j)=x(:,j)-xm(j) 1539
xv(j)=dot_product(w,x(:,j)**2) 1540
if(isd .le. 0)goto 12661 1540
xs(j)=sqrt(xv(j)) 1540
x(:,j)=x(:,j)/xs(j) 1540
xv(j)=1.0 1540
12661 continue 1541
12641 continue 1542
12642 continue 1542
return 1543
end 1544
subroutine lognet2n(parm,no,ni,x,y,g,w,ju,vp,cl,ne,nx,nlam,flmin,u 1546
*lam,shri, isd,intr,maxit,kopt,lmu,a0,a,m,kin,dev0,dev,alm,nlp,jer
*r)
real x(no,ni),y(no),g(no),w(no),vp(ni),ulam(nlam),cl(2,ni) 1547
real a(nx,nlam),a0(nlam),dev(nlam),alm(nlam) 1548
integer ju(ni),m(nx),kin(nlam) 1549
real, dimension (:), allocatable :: b,bs,v,r,xv,q,ga
integer, dimension (:), allocatable :: mm,ixx
call get_int_parms(sml,eps,big,mnlam,devmax,pmin,exmx) 1554
allocate(b(0:ni),stat=jerr) 1555
allocate(xv(1:ni),stat=ierr) 1555
jerr=jerr+ierr 1556
allocate(ga(1:ni),stat=ierr) 1556
jerr=jerr+ierr 1557
allocate(bs(0:ni),stat=ierr) 1557
jerr=jerr+ierr 1558
allocate(mm(1:ni),stat=ierr) 1558
jerr=jerr+ierr 1559
allocate(ixx(1:ni),stat=ierr) 1559
jerr=jerr+ierr 1560
allocate(r(1:no),stat=ierr) 1560
jerr=jerr+ierr 1561
allocate(v(1:no),stat=ierr) 1561
jerr=jerr+ierr 1562
allocate(q(1:no),stat=ierr) 1562
jerr=jerr+ierr 1563
if(jerr.ne.0) return 1564
fmax=log(1.0/pmin-1.0) 1564
fmin=-fmax 1564
vmin=(1.0+pmin)*pmin*(1.0-pmin) 1565
bta=parm 1565
omb=1.0-bta 1566
q0=dot_product(w,y) 1566
if(q0 .gt. pmin)goto 12681 1566
jerr=8001 1566
return 1566
12681 continue 1567
if(q0 .lt. 1.0-pmin)goto 12701 1567
jerr=9001 1567
return 1567
12701 continue 1568
if(intr.eq.0.0) q0=0.5 1569
ixx=0 1569
al=0.0 1569
bz=0.0 1569
if(intr.ne.0) bz=log(q0/(1.0-q0)) 1570
if(nonzero(no,g) .ne. 0)goto 12721 1570
vi=q0*(1.0-q0) 1570
b(0)=bz 1570
v=vi*w 1571
r=w*(y-q0) 1571
q=q0 1571
xmz=vi 1571
dev1=-(bz*q0+log(1.0-q0)) 1572
goto 12731 1573
12721 continue 1573
b(0)=0.0 1574
if(intr .eq. 0)goto 12751 1574
b(0)=azero(no,y,g,w,jerr) 1574
if(jerr.ne.0) return 1574
12751 continue 1575
q=1.0/(1.0+exp(-b(0)-g)) 1575
v=w*q*(1.0-q) 1575
r=w*(y-q) 1575
xmz=sum(v) 1576
dev1=-(b(0)*q0+dot_product(w,y*g+log(1.0-q))) 1577
12731 continue 1578
12711 continue 1578
if(kopt .le. 0)goto 12771 1579
if(isd .le. 0 .or. intr .eq. 0)goto 12791 1579
xv=0.25 1579
goto 12801 1580
12791 continue 1580
12810 do 12811 j=1,ni 1580
if(ju(j).ne.0) xv(j)=0.25*dot_product(w,x(:,j)**2) 1580
12811 continue 1580
12812 continue 1580
12801 continue 1581
12781 continue 1581
12771 continue 1582
dev0=dev1 1583
12820 do 12821 i=1,no 1583
if(y(i).gt.0.0) dev0=dev0+w(i)*y(i)*log(y(i)) 1584
if(y(i).lt.1.0) dev0=dev0+w(i)*(1.0-y(i))*log(1.0-y(i)) 1585
12821 continue 1586
12822 continue 1586
if(flmin .ge. 1.0)goto 12841 1586
eqs=max(eps,flmin) 1586
alf=eqs**(1.0/(nlam-1)) 1586
12841 continue 1587
m=0 1587
mm=0 1587
nlp=0 1587
nin=nlp 1587
mnl=min(mnlam,nlam) 1587
bs=0.0 1587
b(1:ni)=0.0 1588
shr=shri*dev0 1589
12850 do 12851 j=1,ni 1589
if(ju(j).eq.0)goto 12851 1589
ga(j)=abs(dot_product(r,x(:,j))) 1589
12851 continue 1590
12852 continue 1590
12860 do 12861 ilm=1,nlam 1590
al0=al 1591
if(flmin .lt. 1.0)goto 12881 1591
al=ulam(ilm) 1591
goto 12871 1592
12881 if(ilm .le. 2)goto 12891 1592
al=al*alf 1592
goto 12871 1593
12891 if(ilm .ne. 1)goto 12901 1593
al=big 1593
goto 12911 1594
12901 continue 1594
al0=0.0 1595
12920 do 12921 j=1,ni 1595
if(ju(j).eq.0)goto 12921 1595
if(vp(j).gt.0.0) al0=max(al0,ga(j)/vp(j)) 1595
12921 continue 1596
12922 continue 1596
al0=al0/max(bta,1.0e-3) 1596
al=alf*al0 1597
12911 continue 1598
12871 continue 1598
al2=al*omb 1598
al1=al*bta 1598
tlam=bta*(2.0*al-al0) 1599
12930 do 12931 k=1,ni 1599
if(ixx(k).eq.1)goto 12931 1599
if(ju(k).eq.0)goto 12931 1600
if(ga(k).gt.tlam*vp(k)) ixx(k)=1 1601
12931 continue 1602
12932 continue 1602
10880 continue 1603
12940 continue 1603
12941 continue 1603
bs(0)=b(0) 1603
if(nin.gt.0) bs(m(1:nin))=b(m(1:nin)) 1604
if(kopt .ne. 0)goto 12961 1605
12970 do 12971 j=1,ni 1605
if(ixx(j).gt.0) xv(j)=dot_product(v,x(:,j)**2) 1605
12971 continue 1606
12972 continue 1606
12961 continue 1607
12980 continue 1607
12981 continue 1607
nlp=nlp+1 1607
dlx=0.0 1608
12990 do 12991 k=1,ni 1608
if(ixx(k).eq.0)goto 12991 1609
bk=b(k) 1609
gk=dot_product(r,x(:,k)) 1610
u=gk+xv(k)*b(k) 1610
au=abs(u)-vp(k)*al1 1611
if(au .gt. 0.0)goto 13011 1611
b(k)=0.0 1611
goto 13021 1612
13011 continue 1613
b(k)=max(cl(1,k),min(cl(2,k),sign(au,u)/(xv(k)+vp(k)*al2))) 1614
13021 continue 1615
13001 continue 1615
d=b(k)-bk 1615
if(abs(d).le.0.0)goto 12991 1615
dlx=max(dlx,xv(k)*d**2) 1616
r=r-d*v*x(:,k) 1617
if(mm(k) .ne. 0)goto 13041 1617
nin=nin+1 1617
if(nin.gt.nx)goto 12992 1618
mm(k)=nin 1618
m(nin)=k 1619
13041 continue 1620
12991 continue 1621
12992 continue 1621