c  [this is src/GLMnet.f90 from glmnet-1.1-3 by Jerome Friedman, Trevor Hastie, Rob Tibshirani]
c  [from: http://cran.r-project.org/web/packages/glmnet/index.html]
c  [License: GPL-2]
c
c                          GLMnet (5/17/08)
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,ne,nx,nlam,flmin,ulam,thr,isd,
c            lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
c
c   x(no,ni) = predictor data matrix flat file
c
c
c sparse predictor matrix:
c
c call spelnet(ka,parm,no,ni,x,ix,jx,y,w,jd,vp,ne,nx,nlam,flmin,ulam,thr,
c             isd,lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
c
c   x, ix, jx = predictor data matrix in compressed sparse column format
c
c
c other inputs:
c
c   ka = algorithm flag
c      ka=1 => covariance updating algorithm
c      ka=2 => naive algorithm
c   parm = family 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
c   w(no)= observation weights
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   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 standardized coefficient
c      change from the previous iteration is less than thr
c      (suggested value, thr=1.0e-4)
c   isd = 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
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 = total passes over the data summed over 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
c
c Note: x, y and w are overwritten by programs
c
c
c least-squares utility routines:
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          Symmetric binomial/multinomial logistic elastic net
c
c
c dense predictor matrix:
c
c call lognet (parm,no,ni,nc,x,y,jd,vp,ne,nx,nlam,flmin,ulam,thr,isd,
c              maxit,kopt,lmu,a0,ca,ia,nin,dev,alm,nlp,jerr)
c
c   x(no,ni) = predictor data matrix flat file
c
c
c sparse predictor matrix:
c
c call splognet (parm,no,ni,nc,x,ix,jx,y,jd,vp,ne,nx,nlam,flmin,
c             ulam,thr,isd,maxit,kopt,lmu,a0,ca,ia,nin,dev,alm,nlp,jerr)
c
c   x, ix, jx = predictor data matrix in compressed sparse column format
c
c
c other inputs:
c
c   parm, no, ni, jd, vp, ne, nx, nlam, flmin, ulam, thr, isd, same as 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(overwritten)
c   maxit = maximum number of iterations allowed for any lamda value
c           (suggested value, maxit = 100)
c   kopt = optimization flag
c      kopt = 0 => Newton-Raphson
c      kpot = 1 => modified Newton-Raphson (recommended)
c
c
c output:
c
c   lmu, ia, nin, alm, nlp, same as 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   dev(lmu) = fraction of explained devience for 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 < 0 => non fatal error - partial output returned
c         jerr = -k => convergence for kth lamda value not reached
c            after maxit (see above) iterations. Solutions for
c            larger lamdas returned
c         jerr = -10000-k => number of non zero coefficients along path
c            exceeds nx (see above) at kth lamda value. Solutions for
c            larger lamdas returned
c
c
c
c logistic utilitity routines:
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              
      subroutine elnet  (ka,parm,no,ni,x,y,w,jd,vp,ne,nx,nlam,flmin,ulam    352 
     *,thr,isd,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
      real x(no,ni),y(no),w(no),vp(ni),ca(nx,nlam)                          353
      real ulam(nlam),a0(nlam),rsq(nlam),alm(nlam)                          354
      integer jd(*),ia(nx),nin(nlam)                                        355
      real, dimension (:), allocatable :: vq;                                   
      if(maxval(vp) .gt. 0.0)goto 10021                                     358
      jerr=10000                                                            358
      return                                                                358
10021 continue                                                              359
      allocate(vq(1:ni),stat=jerr)                                          359
      if(jerr.ne.0) return                                                  360
      vq=max(0.0,vp)                                                        360
      vq=vq*ni/sum(vq)                                                      361
      if(ka .ne. 1)goto 10041                                               362
      call elnetu  (parm,no,ni,x,y,w,jd,vq,ne,nx,nlam,flmin,ulam,thr,isd    365 
     *,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
      goto 10051                                                            366
10041 continue                                                              367
      call elnetn (parm,no,ni,x,y,w,jd,vq,ne,nx,nlam,flmin,ulam,thr,isd,    370 
     *  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
10051 continue                                                              371
10031 continue                                                              371
      deallocate(vq)                                                        372
      return                                                                373
      end                                                                   374
      subroutine elnetu  (parm,no,ni,x,y,w,jd,vp,ne,nx,nlam,flmin,ulam,t    377 
     *hr,isd,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
      real x(no,ni),y(no),w(no),vp(ni),ulam(nlam)                           378
      real ca(nx,nlam),a0(nlam),rsq(nlam),alm(nlam)                         379
      integer jd(*),ia(nx),nin(nlam)                                        380
      real, dimension (:), allocatable :: xm,xs,g,xv,vlam                       
      integer, dimension (:), allocatable :: ju                                 
      allocate(g(1:ni),stat=jerr)                                           385
      allocate(xm(1:ni),stat=ierr)                                          385
      jerr=jerr+ierr                                                        386
      allocate(xs(1:ni),stat=ierr)                                          386
      jerr=jerr+ierr                                                        387
      allocate(ju(1:ni),stat=ierr)                                          387
      jerr=jerr+ierr                                                        388
      allocate(xv(1:ni),stat=ierr)                                          388
      jerr=jerr+ierr                                                        389
      allocate(vlam(1:nlam),stat=ierr)                                      389
      jerr=jerr+ierr                                                        390
      if(jerr.ne.0) return                                                  391
      call chkvars(no,ni,x,ju)                                              392
      if(jd(1).gt.0) ju(jd(2:(jd(1)+1)))=0                                  393
      if(maxval(ju) .gt. 0)goto 10071                                       393
      jerr=7777                                                             393
      return                                                                393
10071 continue                                                              394
      call standard(no,ni,x,y,w,isd,ju,g,xm,xs,ym,ys,xv,jerr)               395
      if(jerr.ne.0) return                                                  396
      if(flmin.ge.1.0) vlam=ulam/ys                                         397
      call elnet1(parm,ni,ju,vp,g,no,ne,nx,x,nlam,flmin,vlam,thr,xv,  lm    399 
     *u,ca,ia,nin,rsq,alm,nlp,jerr)
      if(jerr.ne.0) return                                                  400
10080 do 10081 k=1,lmu                                                      400
      alm(k)=ys*alm(k)                                                      400
      nk=nin(k)                                                             401
10090 do 10091 l=1,nk                                                       401
      ca(l,k)=ys*ca(l,k)/xs(ia(l))                                          401
10091 continue                                                              402
10092 continue                                                              402
      a0(k)=ym-dot_product(ca(1:nk,k),xm(ia(1:nk)))                         403
10081 continue                                                              404
10082 continue                                                              404
      deallocate(xm,xs,g,ju,xv,vlam)                                        405
      return                                                                406
      end                                                                   407
      subroutine standard (no,ni,x,y,w,isd,ju,g,xm,xs,ym,ys,xv,jerr)        408
      real x(no,ni),y(no),w(no),g(ni),xm(ni),xs(ni),xv(ni)                  408
      integer ju(ni)                                                        409
      real, dimension (:), allocatable :: v                                     
      allocate(v(1:no),stat=jerr)                                           412
      if(jerr.ne.0) return                                                  413
      w=w/sum(w)                                                            413
      v=sqrt(w)                                                             414
10100 do 10101 j=1,ni                                                       414
      if(ju(j).eq.0)goto 10101                                              415
      xm(j)=dot_product(w,x(:,j))                                           415
      x(:,j)=v*(x(:,j)-xm(j))                                               416
      xv(j)=dot_product(x(:,j),x(:,j))                                      417
10101 continue                                                              418
10102 continue                                                              418
      if(isd .ne. 0)goto 10121                                              418
      xs=1.0                                                                418
      goto 10131                                                            419
10121 continue                                                              420
10140 do 10141 j=1,ni                                                       420
      if(ju(j).eq.0)goto 10141                                              420
      xs(j)=sqrt(xv(j))                                                     420
      x(:,j)=x(:,j)/xs(j)                                                   420
10141 continue                                                              421
10142 continue                                                              421
      xv=1.0                                                                422
10131 continue                                                              423
10111 continue                                                              423
      ym=dot_product(w,y)                                                   423
      y=v*(y-ym)                                                            423
      ys=sqrt(dot_product(y,y))                                             423
      y=y/ys                                                                423
      g=0.0                                                                 424
10150 do 10151 j=1,ni                                                       424
      if(ju(j).ne.0) g(j)=dot_product(y,x(:,j))                             424
10151 continue                                                              425
10152 continue                                                              425
      deallocate(v)                                                         426
      return                                                                427
      end                                                                   428
      subroutine elnet1 (beta,ni,ju,vp,g,no,ne,nx,x,nlam,flmin,ulam,thr,    430 
     *xv,  lmu,ao,ia,kin,rsqo,almo,nlp,jerr)
      parameter(sml=1.0e-5, eps=1.0e-6, big=9.9e30, mnlam=5, rsqmax=0.99    431 
     *9)
      real vp(ni),g(ni),x(no,ni),ulam(nlam),ao(nx,nlam),rsqo(nlam),almo(    432 
     *nlam),xv(ni)
      integer ju(ni),ia(nx),kin(nlam)                                       433
      real, dimension (:), allocatable :: a,da                                  
      integer, dimension (:), allocatable :: mm                                 
      real, dimension (:,:), allocatable :: c                                   
      allocate(c(1:ni,1:nx),stat=jerr)                                          
      allocate(a(1:ni),stat=ierr)                                           439
      jerr=jerr+ierr                                                        440
      allocate(mm(1:ni),stat=ierr)                                          440
      jerr=jerr+ierr                                                        441
      allocate(da(1:ni),stat=ierr)                                          441
      jerr=jerr+ierr                                                        442
      if(jerr.ne.0) return                                                  443
      bta=max(beta,1.0e-3)                                                  443
      omb=1.0-bta                                                           444
      if(flmin .ge. 1.0)goto 10171                                          444
      eqs=max(eps,flmin)                                                    444
      alf=eqs**(1.0/(nlam-1))                                               444
10171 continue                                                              445
      rsq=0.0                                                               445
      a=0.0                                                                 445
      mm=0                                                                  445
      nlp=0                                                                 445
      nin=nlp                                                               445
      iz=0                                                                  445
      mnl=min(mnlam,nlam)                                                   446
10180 do 10181 m=1,nlam                                                     447
      if(flmin .lt. 1.0)goto 10201                                          447
      alm=ulam(m)                                                           447
      goto 10191                                                            448
10201 if(m .le. 2)goto 10211                                                448
      alm=alm*alf                                                           448
      goto 10191                                                            449
10211 if(m .ne. 1)goto 10221                                                449
      alm=big                                                               449
      goto 10231                                                            450
10221 continue                                                              450
      alm=0.0                                                               451
10240 do 10241 j=1,ni                                                       451
      if(ju(j).eq.0)goto 10241                                              451
      if(vp(j).le.0.0)goto 10241                                            452
      alm=max(alm,abs(g(j))/vp(j))                                          453
10241 continue                                                              454
10242 continue                                                              454
      alm=alf*alm/bta                                                       455
10231 continue                                                              456
10191 continue                                                              456
      dem=alm*omb                                                           456
      ab=alm*bta                                                            456
      rsq0=rsq                                                              456
      jz=1                                                                  457
10250 continue                                                              457
10251 continue                                                              457
      if(iz*jz.ne.0) go to 10260                                            457
      nlp=nlp+1                                                             457
      dlx=0.0                                                               458
10270 do 10271 k=1,ni                                                       458
      if(ju(k).eq.0)goto 10271                                              459
      ak=a(k)                                                               459
      u=g(k)+ak*xv(k)                                                       459
      v=abs(u)-vp(k)*ab                                                     459
      a(k)=0.0                                                              460
      if(v.gt.0.0) a(k)=sign(v,u)/(xv(k)+vp(k)*dem)                         461
      if(a(k).eq.ak)goto 10271                                              462
      if(mm(k) .ne. 0)goto 10291                                            462
      nin=nin+1                                                             462
      if(nin.gt.nx)goto 10272                                               463
10300 do 10301 j=1,ni                                                       463
      if(ju(j).eq.0)goto 10301                                              464
      if(mm(j) .eq. 0)goto 10321                                            464
      c(j,nin)=c(k,mm(j))                                                   464
      goto 10301                                                            464
10321 continue                                                              465
      if(j .ne. k)goto 10341                                                465
      c(j,nin)=xv(j)                                                        465
      goto 10301                                                            465
10341 continue                                                              466
      c(j,nin)=dot_product(x(:,j),x(:,k))                                   467
10301 continue                                                              468
10302 continue                                                              468
      mm(k)=nin                                                             468
      ia(nin)=k                                                             469
10291 continue                                                              470
      del=a(k)-ak                                                           470
      rsq=rsq+del*(2.0*g(k)-del*xv(k))                                      471
      dlx=max(abs(del)/sqrt(xv(k)),dlx)                                     472
10350 do 10351 j=1,ni                                                       472
      if(ju(j).ne.0) g(j)=g(j)-c(j,mm(k))*del                               472
10351 continue                                                              473
10352 continue                                                              473
10271 continue                                                              474
10272 continue                                                              474
      if(dlx.lt.thr)goto 10252                                              474
      if(nin.gt.nx)goto 10252                                               475
10260 continue                                                              475
      iz=1                                                                  475
      da(1:nin)=a(ia(1:nin))                                                476
10360 continue                                                              476
10361 continue                                                              476
      nlp=nlp+1                                                             476
      dlx=0.0                                                               477
10370 do 10371 l=1,nin                                                      477
      k=ia(l)                                                               477
      ak=a(k)                                                               477
      u=g(k)+ak*xv(k)                                                       477
      v=abs(u)-vp(k)*ab                                                     478
      a(k)=0.0                                                              479
      if(v.gt.0.0) a(k)=sign(v,u)/(xv(k)+vp(k)*dem)                         480
      if(a(k).eq.ak)goto 10371                                              481
      del=a(k)-ak                                                           481
      rsq=rsq+del*(2.0*g(k)-del*xv(k))                                      482
      dlx=max(abs(del)/sqrt(xv(k)),dlx)                                     483
10380 do 10381 j=1,nin                                                      483
      g(ia(j))=g(ia(j))-c(ia(j),mm(k))*del                                  483
10381 continue                                                              484
10382 continue                                                              484
10371 continue                                                              485
10372 continue                                                              485
      if(dlx.lt.thr)goto 10362                                              485
      goto 10361                                                            486
10362 continue                                                              486
      da(1:nin)=a(ia(1:nin))-da(1:nin)                                      487
10390 do 10391 j=1,ni                                                       487
      if(mm(j).ne.0)goto 10391                                              488
      if(ju(j).ne.0) g(j)=g(j)-dot_product(da(1:nin),c(j,1:nin))            489
10391 continue                                                              490
10392 continue                                                              490
      jz=0                                                                  491
      goto 10251                                                            492
10252 continue                                                              492
      if(nin.gt.nx)goto 10182                                               493
      if(nin.gt.0) ao(1:nin,m)=a(ia(1:nin))                                 493
      kin(m)=nin                                                            494
      rsqo(m)=rsq                                                           494
      almo(m)=alm                                                           494
      lmu=m                                                                 495
      if(m.lt.mnl)goto 10181                                                495
      if(flmin.ge.1.0)goto 10181                                            496
      me=0                                                                  496
10400 do 10401 j=1,nin                                                      496
      if(ao(j,m).ne.0.0) me=me+1                                            496
10401 continue                                                              496
10402 continue                                                              496
      if(me.gt.ne)goto 10182                                                497
      if(rsq-rsq0.lt.sml*rsq)goto 10182                                     497
      if(rsq.gt.rsqmax)goto 10182                                           498
10181 continue                                                              499
10182 continue                                                              499
      deallocate(a,mm,c,da)                                                 500
      return                                                                501
      end                                                                   502
      subroutine elnetn (parm,no,ni,x,y,w,jd,vp,ne,nx,nlam,flmin,ulam,th    504 
     *r,isd,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
      real vp(ni),x(no,ni),y(no),w(no),ulam(nlam)                           505
      real ca(nx,nlam),a0(nlam),rsq(nlam),alm(nlam)                         506
      integer jd(*),ia(nx),nin(nlam)                                        507
      real, dimension (:), allocatable :: xm,xs,xv,vlam                         
      integer, dimension (:), allocatable :: ju                                 
      allocate(xm(1:ni),stat=jerr)                                          512
      allocate(xs(1:ni),stat=ierr)                                          512
      jerr=jerr+ierr                                                        513
      allocate(ju(1:ni),stat=ierr)                                          513
      jerr=jerr+ierr                                                        514
      allocate(xv(1:ni),stat=ierr)                                          514
      jerr=jerr+ierr                                                        515
      allocate(vlam(1:nlam),stat=ierr)                                      515
      jerr=jerr+ierr                                                        516
      if(jerr.ne.0) return                                                  517
      call chkvars(no,ni,x,ju)                                              518
      if(jd(1).gt.0) ju(jd(2:(jd(1)+1)))=0                                  519
      if(maxval(ju) .gt. 0)goto 10421                                       519
      jerr=7777                                                             519
      return                                                                519
10421 continue                                                              520
      call standard1(no,ni,x,y,w,isd,ju,xm,xs,ym,ys,xv,jerr)                521
      if(jerr.ne.0) return                                                  522
      if(flmin.ge.1.0) vlam=ulam/ys                                         523
      call elnet2(parm,ni,ju,vp,y,no,ne,nx,x,nlam,flmin,vlam,thr,xv,  lm    525 
     *u,ca,ia,nin,rsq,alm,nlp,jerr)
      if(jerr.ne.0) return                                                  526
10430 do 10431 k=1,lmu                                                      526
      alm(k)=ys*alm(k)                                                      526
      nk=nin(k)                                                             527
10440 do 10441 l=1,nk                                                       527
      ca(l,k)=ys*ca(l,k)/xs(ia(l))                                          527
10441 continue                                                              528
10442 continue                                                              528
      a0(k)=ym-dot_product(ca(1:nk,k),xm(ia(1:nk)))                         529
10431 continue                                                              530
10432 continue                                                              530
      deallocate(xm,xs,ju,xv,vlam)                                          531
      return                                                                532
      end                                                                   533
      subroutine standard1 (no,ni,x,y,w,isd,ju,xm,xs,ym,ys,xv,jerr)         534
      real x(no,ni),y(no),w(no),xm(ni),xs(ni),xv(ni)                        534
      integer ju(ni)                                                        535
      real, dimension (:), allocatable :: v                                     
      allocate(v(1:no),stat=jerr)                                           538
      if(jerr.ne.0) return                                                  539
      w=w/sum(w)                                                            539
      v=sqrt(w)                                                             540
10450 do 10451 j=1,ni                                                       540
      if(ju(j).eq.0)goto 10451                                              541
      xm(j)=dot_product(w,x(:,j))                                           541
      x(:,j)=v*(x(:,j)-xm(j))                                               542
      xv(j)=dot_product(x(:,j),x(:,j))                                      543
10451 continue                                                              544
10452 continue                                                              544
      if(isd .ne. 0)goto 10471                                              544
      xs=1.0                                                                544
      goto 10481                                                            545
10471 continue                                                              545
10490 do 10491 j=1,ni                                                       545
      if(ju(j).eq.0)goto 10491                                              545
      xs(j)=sqrt(xv(j))                                                     545
      x(:,j)=x(:,j)/xs(j)                                                   545
10491 continue                                                              546
10492 continue                                                              546
      xv=1.0                                                                547
10481 continue                                                              548
10461 continue                                                              548
      ym=dot_product(w,y)                                                   548
      y=v*(y-ym)                                                            548
      ys=sqrt(dot_product(y,y))                                             548
      y=y/ys                                                                549
      deallocate(v)                                                         550
      return                                                                551
      end                                                                   552
      subroutine elnet2(beta,ni,ju,vp,y,no,ne,nx,x,nlam,flmin,ulam,thr,x    554 
     *v,  lmu,ao,ia,kin,rsqo,almo,nlp,jerr)
      parameter(sml=1.0e-5, eps=1.0e-6, big=9.9e30, mnlam=5, rsqmax=0.99    555 
     *9)
      real vp(ni),y(no),x(no,ni),ulam(nlam),ao(nx,nlam),rsqo(nlam),almo(    556 
     *nlam),xv(ni)
      integer ju(ni),ia(nx),kin(nlam)                                       557
      real, dimension (:), allocatable :: a                                     
      integer, dimension (:), allocatable :: mm                                 
      allocate(a(1:ni),stat=jerr)                                           562
      allocate(mm(1:ni),stat=ierr)                                          562
      jerr=jerr+ierr                                                        563
      if(jerr.ne.0) return                                                  564
      bta=max(beta,1.0e-3)                                                  564
      omb=1.0-bta                                                           565
      if(flmin .ge. 1.0)goto 10511                                          565
      eqs=max(eps,flmin)                                                    565
      alf=eqs**(1.0/(nlam-1))                                               565
10511 continue                                                              566
      rsq=0.0                                                               566
      a=0.0                                                                 566
      mm=0                                                                  566
      nlp=0                                                                 566
      nin=nlp                                                               566
      iz=0                                                                  566
      mnl=min(mnlam,nlam)                                                   567
10520 do 10521 m=1,nlam                                                     568
      if(flmin .lt. 1.0)goto 10541                                          568
      alm=ulam(m)                                                           568
      goto 10531                                                            569
10541 if(m .le. 2)goto 10551                                                569
      alm=alm*alf                                                           569
      goto 10531                                                            570
10551 if(m .ne. 1)goto 10561                                                570
      alm=big                                                               570
      goto 10571                                                            571
10561 continue                                                              571
      alm=0.0                                                               572
10580 do 10581 j=1,ni                                                       572
      if(ju(j).eq.0)goto 10581                                              572
      if(vp(j).le.0.0)goto 10581                                            573
      alm=max(alm,abs(dot_product(y,x(:,j)))/vp(j))                         574
10581 continue                                                              575
10582 continue                                                              575
      alm=alf*alm/bta                                                       576
10571 continue                                                              577
10531 continue                                                              577
      dem=alm*omb                                                           577
      ab=alm*bta                                                            577
      rsq0=rsq                                                              577
      jz=1                                                                  578
10590 continue                                                              578
10591 continue                                                              578
      if(iz*jz.ne.0) go to 10260                                            578
      nlp=nlp+1                                                             578
      dlx=0.0                                                               579
10600 do 10601 k=1,ni                                                       579
      if(ju(k).eq.0)goto 10601                                              579
      gk=dot_product(y,x(:,k))                                              580
      ak=a(k)                                                               580
      u=gk+ak*xv(k)                                                         580
      v=abs(u)-vp(k)*ab                                                     580
      a(k)=0.0                                                              581
      if(v.gt.0.0) a(k)=sign(v,u)/(xv(k)+vp(k)*dem)                         582
      if(a(k).eq.ak)goto 10601                                              583
      if(mm(k) .ne. 0)goto 10621                                            583
      nin=nin+1                                                             583
      if(nin.gt.nx)goto 10602                                               584
      mm(k)=nin                                                             584
      ia(nin)=k                                                             585
10621 continue                                                              586
      del=a(k)-ak                                                           586
      rsq=rsq+del*(2.0*gk-del*xv(k))                                        587
      y=y-del*x(:,k)                                                        587
      dlx=max(abs(del)/sqrt(xv(k)),dlx)                                     588
10601 continue                                                              589
10602 continue                                                              589
      if(dlx.lt.thr)goto 10592                                              589
      if(nin.gt.nx)goto 10592                                               590
10260 continue                                                              590
      iz=1                                                                  591
10630 continue                                                              591
10631 continue                                                              591
      nlp=nlp+1                                                             591
      dlx=0.0                                                               592
10640 do 10641 l=1,nin                                                      592
      k=ia(l)                                                               592
      gk=dot_product(y,x(:,k))                                              593
      ak=a(k)                                                               593
      u=gk+ak*xv(k)                                                         593
      v=abs(u)-vp(k)*ab                                                     593
      a(k)=0.0                                                              594
      if(v.gt.0.0) a(k)=sign(v,u)/(xv(k)+vp(k)*dem)                         595
      if(a(k).eq.ak)goto 10641                                              596
      del=a(k)-ak                                                           596
      rsq=rsq+del*(2.0*gk-del*xv(k))                                        597
      y=y-del*x(:,k)                                                        597
      dlx=max(abs(del)/sqrt(xv(k)),dlx)                                     598
10641 continue                                                              599
10642 continue                                                              599
      if(dlx.lt.thr)goto 10632                                              599
      goto 10631                                                            600
10632 continue                                                              600
      jz=0                                                                  601
      goto 10591                                                            602
10592 continue                                                              602
      if(nin.gt.nx)goto 10522                                               603
      if(nin.gt.0) ao(1:nin,m)=a(ia(1:nin))                                 603
      kin(m)=nin                                                            604
      rsqo(m)=rsq                                                           604
      almo(m)=alm                                                           604
      lmu=m                                                                 605
      if(m.lt.mnl)goto 10521                                                605
      if(flmin.ge.1.0)goto 10521                                            606
      me=0                                                                  606
10650 do 10651 j=1,nin                                                      606
      if(ao(j,m).ne.0.0) me=me+1                                            606
10651 continue                                                              606
10652 continue                                                              606
      if(me.gt.ne)goto 10522                                                607
      if(rsq-rsq0.lt.sml*rsq)goto 10522                                     607
      if(rsq.gt.rsqmax)goto 10522                                           608
10521 continue                                                              609
10522 continue                                                              609
      deallocate(a,mm)                                                      610
      return                                                                611
      end                                                                   612
      subroutine chkvars(no,ni,x,ju)                                        613
      real x(no,ni)                                                         613
      integer ju(ni)                                                        614
10660 do 10661 j=1,ni                                                       614
      ju(j)=0                                                               614
      t=x(1,j)                                                              615
10670 do 10671 i=2,no                                                       615
      if(x(i,j).eq.t)goto 10671                                             615
      ju(j)=1                                                               615
      goto 10672                                                            615
10671 continue                                                              616
10672 continue                                                              616
10661 continue                                                              617
10662 continue                                                              617
      return                                                                618
      end                                                                   619
      subroutine uncomp(ni,ca,ia,nin,a)                                     620
      real ca(*),a(ni)                                                      620
      integer ia(*)                                                         621
      a=0.0                                                                 621
      if(nin.gt.0) a(ia(1:nin))=ca(1:nin)                                   622
      return                                                                623
      end                                                                   624
      subroutine modval(a0,ca,ia,nin,n,x,f)                                 625
      real ca(nin),x(n,*),f(n)                                              625
      integer ia(nin)                                                       626
      f=a0                                                                  626
      if(nin.le.0) return                                                   627
10680 do 10681 i=1,n                                                        627
      f(i)=f(i)+dot_product(ca(1:nin),x(i,ia(1:nin)))                       627
10681 continue                                                              628
10682 continue                                                              628
      return                                                                629
      end                                                                   630
      subroutine spelnet  (ka,parm,no,ni,x,ix,jx,y,w,jd,vp,ne,nx,nlam,fl    633 
     *min,ulam,thr,isd,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
      real x(*),y(no),w(no),vp(ni),ulam(nlam)                               634
      real ca(nx,nlam),a0(nlam),rsq(nlam),alm(nlam)                         635
      integer ix(*),jx(*),jd(*),ia(nx),nin(nlam)                            636
      real, dimension (:), allocatable :: vq;                                   
      if(maxval(vp) .gt. 0.0)goto 10701                                     639
      jerr=10000                                                            639
      return                                                                639
10701 continue                                                              640
      allocate(vq(1:ni),stat=jerr)                                          640
      if(jerr.ne.0) return                                                  641
      vq=max(0.0,vp)                                                        641
      vq=vq*ni/sum(vq)                                                      642
      if(ka .ne. 1)goto 10721                                               643
      call spelnetu  (parm,no,ni,x,ix,jx,y,w,jd,vq,ne,nx,nlam,flmin,ulam    646 
     *,thr,isd,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
      goto 10731                                                            647
10721 continue                                                              648
      call spelnetn (parm,no,ni,x,ix,jx,y,w,jd,vq,ne,nx,nlam,flmin,ulam,    651 
     *thr,isd,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
10731 continue                                                              652
10711 continue                                                              652
      deallocate(vq)                                                        653
      return                                                                654
      end                                                                   655
      subroutine spelnetu  (parm,no,ni,x,ix,jx,y,w,jd,vp,ne,nx,nlam,flmi    658 
     *n,ulam,thr,isd,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
      real x(*),y(no),w(no),vp(ni),ulam(nlam)                               659
      real ca(nx,nlam),a0(nlam),rsq(nlam),alm(nlam)                         660
      integer ix(*),jx(*),jd(*),ia(nx),nin(nlam)                            661
      real, dimension (:), allocatable :: xm,xs,g,xv,vlam                       
      integer, dimension (:), allocatable :: ju                                 
      allocate(g(1:ni),stat=jerr)                                           666
      allocate(xm(1:ni),stat=ierr)                                          666
      jerr=jerr+ierr                                                        667
      allocate(xs(1:ni),stat=ierr)                                          667
      jerr=jerr+ierr                                                        668
      allocate(ju(1:ni),stat=ierr)                                          668
      jerr=jerr+ierr                                                        669
      allocate(xv(1:ni),stat=ierr)                                          669
      jerr=jerr+ierr                                                        670
      allocate(vlam(1:nlam),stat=ierr)                                      670
      jerr=jerr+ierr                                                        671
      if(jerr.ne.0) return                                                  672
      call spchkvars(no,ni,x,ix,ju)                                         673
      if(jd(1).gt.0) ju(jd(2:(jd(1)+1)))=0                                  674
      if(maxval(ju) .gt. 0)goto 10751                                       674
      jerr=7777                                                             674
      return                                                                674
10751 continue                                                              675
      call spstandard(no,ni,x,ix,jx,y,w,ju,isd,g,xm,xs,ym,ys,xv,jerr)       676
      if(jerr.ne.0) return                                                  677
      if(flmin.ge.1.0) vlam=ulam/ys                                         678
      call spelnet1(parm,ni,g,no,w,ne,nx,x,ix,jx,ju,vp,nlam,flmin,vlam,t    680 
     *hr,  xm,xs,xv,lmu,ca,ia,nin,rsq,alm,nlp,jerr)
      if(jerr.ne.0) return                                                  681
10760 do 10761 k=1,lmu                                                      681
      alm(k)=ys*alm(k)                                                      681
      nk=nin(k)                                                             682
10770 do 10771 l=1,nk                                                       682
      ca(l,k)=ys*ca(l,k)/xs(ia(l))                                          682
10771 continue                                                              683
10772 continue                                                              683
      a0(k)=ym-dot_product(ca(1:nk,k),xm(ia(1:nk)))                         684
10761 continue                                                              685
10762 continue                                                              685
      deallocate(xm,xs,g,ju,xv,vlam)                                        686
      return                                                                687
      end                                                                   688
      subroutine spstandard (no,ni,x,ix,jx,y,w,ju,isd,g,xm,xs,ym,ys,xv,j    689 
     *err)
      real x(*),y(no),w(no),g(ni),xm(ni),xs(ni),xv(ni)                      689
      integer ix(*),jx(*),ju(ni)                                            690
      w=w/sum(w)                                                            691
10780 do 10781 j=1,ni                                                       691
      if(ju(j).eq.0)goto 10781                                              692
      jb=ix(j)                                                              692
      je=ix(j+1)-1                                                          692
      xm(j)=dot_product(w(jx(jb:je)),x(jb:je))                              693
      xv(j)=dot_product(w(jx(jb:je)),x(jb:je)**2)-xm(j)**2                  694
10781 continue                                                              695
10782 continue                                                              695
      if(isd .ne. 0)goto 10801                                              695
      xs=1.0                                                                695
      goto 10811                                                            696
10801 continue                                                              696
10820 do 10821 j=1,ni                                                       696
      if(ju(j).ne.0) xs(j)=sqrt(xv(j))                                      696
10821 continue                                                              696
10822 continue                                                              696
      xv=1.0                                                                696
10811 continue                                                              697
10791 continue                                                              697
      ym=dot_product(w,y)                                                   697
      y=y-ym                                                                697
      ys=sqrt(dot_product(w,y**2))                                          697
      y=y/ys                                                                697
      g=0.0                                                                 698
10830 do 10831 j=1,ni                                                       698
      if(ju(j).eq.0)goto 10831                                              698
      jb=ix(j)                                                              698
      je=ix(j+1)-1                                                          699
      g(j)=dot_product(w(jx(jb:je))*y(jx(jb:je)),x(jb:je))/xs(j)            700
10831 continue                                                              701
10832 continue                                                              701
      return                                                                702
      end                                                                   703
      subroutine spelnet1(beta,ni,g,no,w,ne,nx,x,ix,jx,ju,vp,nlam,flmin,    705 
     *ulam,  thr,xm,xs,xv,lmu,ao,ia,kin,rsqo,almo,nlp,jerr)
      parameter(sml=1.0e-5, eps=1.0e-6, big=9.9e30, mnlam=5, rsqmax=0.99    706 
     *9)
      real g(ni),vp(ni),x(*),ulam(nlam),w(no)                               707
      real ao(nx,nlam),rsqo(nlam),almo(nlam),xm(ni),xs(ni),xv(ni)           708
      integer ix(*),jx(*),ju(ni),ia(nx),kin(nlam)                           709
      real, dimension (:), allocatable :: a,da                                  
      integer, dimension (:), allocatable :: mm                                 
      real, dimension (:,:), allocatable :: c                                   
      allocate(c(1:ni,1:nx),stat=jerr)                                          
      allocate(a(1:ni),stat=ierr)                                           715
      jerr=jerr+ierr                                                        716
      allocate(mm(1:ni),stat=ierr)                                          716
      jerr=jerr+ierr                                                        717
      allocate(da(1:ni),stat=ierr)                                          717
      jerr=jerr+ierr                                                        718
      if(jerr.ne.0) return                                                  719
      bta=max(beta,1.0e-3)                                                  719
      omb=1.0-bta                                                           720
      if(flmin .ge. 1.0)goto 10851                                          720
      eqs=max(eps,flmin)                                                    720
      alf=eqs**(1.0/(nlam-1))                                               720
10851 continue                                                              721
      rsq=0.0                                                               721
      a=0.0                                                                 721
      mm=0                                                                  721
      nlp=0                                                                 721
      nin=nlp                                                               721
      iz=0                                                                  721
      mnl=min(mnlam,nlam)                                                   722
10860 do 10861 m=1,nlam                                                     723
      if(flmin .lt. 1.0)goto 10881                                          723
      alm=ulam(m)                                                           723
      goto 10871                                                            724
10881 if(m .le. 2)goto 10891                                                724
      alm=alm*alf                                                           724
      goto 10871                                                            725
10891 if(m .ne. 1)goto 10901                                                725
      alm=big                                                               725
      goto 10911                                                            726
10901 continue                                                              726
      alm=0.0                                                               727
10920 do 10921 j=1,ni                                                       727
      if(ju(j).eq.0)goto 10921                                              727
      if(vp(j).le.0.0)goto 10921                                            728
      alm=max(alm,abs(g(j))/vp(j))                                          729
10921 continue                                                              730
10922 continue                                                              730
      alm=alf*alm/bta                                                       731
10911 continue                                                              732
10871 continue                                                              732
      dem=alm*omb                                                           732
      ab=alm*bta                                                            732
      rsq0=rsq                                                              732
      jz=1                                                                  733
10930 continue                                                              733
10931 continue                                                              733
      if(iz*jz.ne.0) go to 10260                                            733
      nlp=nlp+1                                                             733
      dlx=0.0                                                               734
10940 do 10941 k=1,ni                                                       734
      if(ju(k).eq.0)goto 10941                                              735
      ak=a(k)                                                               735
      u=g(k)+ak*xv(k)                                                       735
      v=abs(u)-vp(k)*ab                                                     735
      a(k)=0.0                                                              736
      if(v.gt.0.0) a(k)=sign(v,u)/(xv(k)+vp(k)*dem)                         737
      if(a(k).eq.ak)goto 10941                                              738
      if(mm(k) .ne. 0)goto 10961                                            738
      nin=nin+1                                                             738
      if(nin.gt.nx)goto 10942                                               739
10970 do 10971 j=1,ni                                                       739
      if(ju(j).eq.0)goto 10971                                              740
      if(mm(j) .eq. 0)goto 10991                                            740
      c(j,nin)=c(k,mm(j))                                                   740
      goto 10971                                                            740
10991 continue                                                              741
      if(j .ne. k)goto 11011                                                741
      c(j,nin)=xv(j)                                                        741
      goto 10971                                                            741
11011 continue                                                              742
      c(j,nin)=  (row_prod(j,k,ix,jx,x,w)-xm(j)*xm(k))/(xs(j)*xs(k))        744
10971 continue                                                              745
10972 continue                                                              745
      mm(k)=nin                                                             745
      ia(nin)=k                                                             746
10961 continue                                                              747
      del=a(k)-ak                                                           747
      rsq=rsq+del*(2.0*g(k)-del*xv(k))                                      748
      dlx=max(abs(del)/sqrt(xv(k)),dlx)                                     749
11020 do 11021 j=1,ni                                                       749
      if(ju(j).ne.0) g(j)=g(j)-c(j,mm(k))*del                               749
11021 continue                                                              750
11022 continue                                                              750
10941 continue                                                              751
10942 continue                                                              751
      if(dlx.lt.thr)goto 10932                                              751
      if(nin.gt.nx)goto 10932                                               752
10260 continue                                                              752
      iz=1                                                                  752
      da(1:nin)=a(ia(1:nin))                                                753
11030 continue                                                              753
11031 continue                                                              753
      nlp=nlp+1                                                             753
      dlx=0.0                                                               754
11040 do 11041 l=1,nin                                                      754
      k=ia(l)                                                               755
      ak=a(k)                                                               755
      u=g(k)+ak*xv(k)                                                       755
      v=abs(u)-vp(k)*ab                                                     755
      a(k)=0.0                                                              756
      if(v.gt.0.0) a(k)=sign(v,u)/(xv(k)+vp(k)*dem)                         757
      if(a(k).eq.ak)goto 11041                                              758
      del=a(k)-ak                                                           758
      rsq=rsq+del*(2.0*g(k)-del*xv(k))                                      759
      dlx=max(abs(del)/sqrt(xv(k)),dlx)                                     760
11050 do 11051 j=1,nin                                                      760
      g(ia(j))=g(ia(j))-c(ia(j),mm(k))*del                                  760
11051 continue                                                              761
11052 continue                                                              761
11041 continue                                                              762
11042 continue                                                              762
      if(dlx.lt.thr)goto 11032                                              762
      goto 11031                                                            763
11032 continue                                                              763
      da(1:nin)=a(ia(1:nin))-da(1:nin)                                      764
11060 do 11061 j=1,ni                                                       764
      if(mm(j).ne.0)goto 11061                                              765
      if(ju(j).ne.0) g(j)=g(j)-dot_product(da(1:nin),c(j,1:nin))            766
11061 continue                                                              767
11062 continue                                                              767
      jz=0                                                                  768
      goto 10931                                                            769
10932 continue                                                              769
      if(nin.gt.nx)goto 10862                                               770
      if(nin.gt.0) ao(1:nin,m)=a(ia(1:nin))                                 770
      kin(m)=nin                                                            771
      rsqo(m)=rsq                                                           771
      almo(m)=alm                                                           771
      lmu=m                                                                 772
      if(m.lt.mnl)goto 10861                                                772
      if(flmin.ge.1.0)goto 10861                                            773
      me=0                                                                  773
11070 do 11071 j=1,nin                                                      773
      if(ao(j,m).ne.0.0) me=me+1                                            773
11071 continue                                                              773
11072 continue                                                              773
      if(me.gt.ne)goto 10862                                                774
      if(rsq-rsq0.lt.sml*rsq)goto 10862                                     774
      if(rsq.gt.rsqmax)goto 10862                                           775
10861 continue                                                              776
10862 continue                                                              776
      deallocate(a,mm,c,da)                                                 777
      return                                                                778
      end                                                                   779
      subroutine spelnetn(parm,no,ni,x,ix,jx,y,w,jd,vp,ne,nx,nlam,flmin,    781 
     *ulam,  thr,isd,lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
      real x(*),vp(ni),y(no),w(no),ulam(nlam)                               782
      real ca(nx,nlam),a0(nlam),rsq(nlam),alm(nlam)                         783
      integer ix(*),jx(*),jd(*),ia(nx),nin(nlam)                            784
      real, dimension (:), allocatable :: xm,xs,xv,vlam                         
      integer, dimension (:), allocatable :: ju                                 
      allocate(xm(1:ni),stat=jerr)                                          789
      allocate(xs(1:ni),stat=ierr)                                          789
      jerr=jerr+ierr                                                        790
      allocate(ju(1:ni),stat=ierr)                                          790
      jerr=jerr+ierr                                                        791
      allocate(xv(1:ni),stat=ierr)                                          791
      jerr=jerr+ierr                                                        792
      allocate(vlam(1:nlam),stat=ierr)                                      792
      jerr=jerr+ierr                                                        793
      if(jerr.ne.0) return                                                  794
      call spchkvars(no,ni,x,ix,ju)                                         795
      if(jd(1).gt.0) ju(jd(2:(jd(1)+1)))=0                                  796
      if(maxval(ju) .gt. 0)goto 11091                                       796
      jerr=7777                                                             796
      return                                                                796
11091 continue                                                              797
      call spstandard1(no,ni,x,ix,jx,y,w,ju,isd,xm,xs,ym,ys,xv,jerr)        798
      if(jerr.ne.0) return                                                  799
      if(flmin.ge.1.0) vlam=ulam/ys                                         800
      call spelnet2(parm,ni,y,w,no,ne,nx,x,ix,jx,ju,vp,nlam,flmin,vlam,t    802 
     *hr,xm,xs,xv,  lmu,ca,ia,nin,rsq,alm,nlp,jerr)
      if(jerr.ne.0) return                                                  803
11100 do 11101 k=1,lmu                                                      803
      alm(k)=ys*alm(k)                                                      803
      nk=nin(k)                                                             804
11110 do 11111 l=1,nk                                                       804
      ca(l,k)=ys*ca(l,k)/xs(ia(l))                                          804
11111 continue                                                              805
11112 continue                                                              805
      a0(k)=ym-dot_product(ca(1:nk,k),xm(ia(1:nk)))                         806
11101 continue                                                              807
11102 continue                                                              807
      deallocate(xm,xs,ju,xv,vlam)                                          808
      return                                                                809
      end                                                                   810
      subroutine spstandard1 (no,ni,x,ix,jx,y,w,ju,isd,xm,xs,ym,ys,xv,je    811 
     *rr)
      real x(*),y(no),w(no),xm(ni),xs(ni),xv(ni)                            811
      integer ix(*),jx(*),ju(ni)                                            812
      w=w/sum(w)                                                            813
11120 do 11121 j=1,ni                                                       813
      if(ju(j).eq.0)goto 11121                                              814
      jb=ix(j)                                                              814
      je=ix(j+1)-1                                                          814
      xm(j)=dot_product(w(jx(jb:je)),x(jb:je))                              815
      xv(j)=dot_product(w(jx(jb:je)),x(jb:je)**2)-xm(j)**2                  816
11121 continue                                                              817
11122 continue                                                              817
      if(isd .ne. 0)goto 11141                                              817
      xs=1.0                                                                817
      goto 11151                                                            818
11141 continue                                                              818
11160 do 11161 j=1,ni                                                       818
      if(ju(j).ne.0) xs(j)=sqrt(xv(j))                                      818
11161 continue                                                              818
11162 continue                                                              818
      xv=1.0                                                                818
11151 continue                                                              819
11131 continue                                                              819
      ym=dot_product(w,y)                                                   819
      y=y-ym                                                                819
      ys=sqrt(dot_product(w,y**2))                                          819
      y=y/ys                                                                820
      return                                                                821
      end                                                                   822
      subroutine spelnet2(beta,ni,y,w,no,ne,nx,x,ix,jx,ju,vp,nlam,flmin,    824 
     *ulam,  thr,xm,xs,xv,lmu,ao,ia,kin,rsqo,almo,nlp,jerr)
      parameter(sml=1.0e-5, eps=1.0e-6, big=9.9e30, mnlam=5, rsqmax=0.99    825 
     *9)
      real y(no),w(no),x(*),vp(ni),ulam(nlam)                               826
      real ao(nx,nlam),rsqo(nlam),almo(nlam),xm(ni),xs(ni),xv(ni)           827
      integer ix(*),jx(*),ju(ni),ia(nx),kin(nlam)                           828
      real, dimension (:), allocatable :: a                                     
      integer, dimension (:), allocatable :: mm                                 
      allocate(a(1:ni),stat=jerr)                                           833
      allocate(mm(1:ni),stat=ierr)                                          833
      jerr=jerr+ierr                                                        834
      if(jerr.ne.0) return                                                  835
      bta=max(beta,1.0e-3)                                                  835
      omb=1.0-bta                                                           836
      if(flmin .ge. 1.0)goto 11181                                          836
      eqs=max(eps,flmin)                                                    836
      alf=eqs**(1.0/(nlam-1))                                               836
11181 continue                                                              837
      rsq=0.0                                                               837
      a=0.0                                                                 837
      mm=0                                                                  837
      o=0.0                                                                 837
      nlp=0                                                                 837
      nin=nlp                                                               837
      iz=0                                                                  837
      mnl=min(mnlam,nlam)                                                   838
11190 do 11191 m=1,nlam                                                     839
      if(flmin .lt. 1.0)goto 11211                                          839
      alm=ulam(m)                                                           839
      goto 11201                                                            840
11211 if(m .le. 2)goto 11221                                                840
      alm=alm*alf                                                           840
      goto 11201                                                            841
11221 if(m .ne. 1)goto 11231                                                841
      alm=big                                                               841
      goto 11241                                                            842
11231 continue                                                              842
      alm=0.0                                                               843
11250 do 11251 j=1,ni                                                       843
      if(ju(j).eq.0)goto 11251                                              843
      if(vp(j).le.0.0)goto 11251                                            844
      jb=ix(j)                                                              844
      je=ix(j+1)-1                                                          845
      alm=max(alm,abs(dot_product(y(jx(jb:je))+o,w(jx(jb:je))*x(jb:je))     847 
     * /(vp(j)*xs(j))))
11251 continue                                                              848
11252 continue                                                              848
      alm=alf*alm/bta                                                       849
11241 continue                                                              850
11201 continue                                                              850
      dem=alm*omb                                                           850
      ab=alm*bta                                                            850
      rsq0=rsq                                                              850
      jz=1                                                                  851
11260 continue                                                              851
11261 continue                                                              851
      if(iz*jz.ne.0) go to 10260                                            851
      nlp=nlp+1                                                             851
      dlx=0.0                                                               852
11270 do 11271 k=1,ni                                                       852
      if(ju(k).eq.0)goto 11271                                              852
      jb=ix(k)                                                              852
      je=ix(k+1)-1                                                          853
      gk=dot_product(y(jx(jb:je))+o,w(jx(jb:je))*x(jb:je))/xs(k)            854
      ak=a(k)                                                               854
      u=gk+ak*xv(k)                                                         854
      v=abs(u)-vp(k)*ab                                                     854
      a(k)=0.0                                                              855
      if(v.gt.0.0) a(k)=sign(v,u)/(xv(k)+vp(k)*dem)                         856
      if(a(k).eq.ak)goto 11271                                              857
      if(mm(k) .ne. 0)goto 11291                                            857
      nin=nin+1                                                             857
      if(nin.gt.nx)goto 11272                                               858
      mm(k)=nin                                                             858
      ia(nin)=k                                                             859
11291 continue                                                              860
      del=a(k)-ak                                                           860
      rsq=rsq+del*(2.0*gk-del*xv(k))                                        861
      y(jx(jb:je))=y(jx(jb:je))-del*x(jb:je)/xs(k)                          862
      o=o+del*xm(k)/xs(k)                                                   862
      dlx=max(abs(del)/sqrt(xv(k)),dlx)                                     863
11271 continue                                                              864
11272 continue                                                              864
      if(dlx.lt.thr)goto 11262                                              864
      if(nin.gt.nx)goto 11262                                               865
10260 continue                                                              865
      iz=1                                                                  866
11300 continue                                                              866
11301 continue                                                              866
      nlp=nlp+1                                                             866
      dlx=0.0                                                               867
11310 do 11311 l=1,nin                                                      867
      k=ia(l)                                                               867
      jb=ix(k)                                                              867
      je=ix(k+1)-1                                                          868
      gk=dot_product(y(jx(jb:je))+o,w(jx(jb:je))*x(jb:je))/xs(k)            869
      ak=a(k)                                                               869
      u=gk+ak*xv(k)                                                         869
      v=abs(u)-vp(k)*ab                                                     869
      a(k)=0.0                                                              870
      if(v.gt.0.0) a(k)=sign(v,u)/(xv(k)+vp(k)*dem)                         871
      if(a(k).eq.ak)goto 11311                                              872
      del=a(k)-ak                                                           872
      rsq=rsq+del*(2.0*gk-del*xv(k))                                        873
      y(jx(jb:je))=y(jx(jb:je))-del*x(jb:je)/xs(k)                          874
      o=o+del*xm(k)/xs(k)                                                   874
      dlx=max(abs(del)/sqrt(xv(k)),dlx)                                     875
11311 continue                                                              876
11312 continue                                                              876
      if(dlx.lt.thr)goto 11302                                              876
      goto 11301                                                            877
11302 continue                                                              877
      jz=0                                                                  878
      goto 11261                                                            879
11262 continue                                                              879
      if(nin.gt.nx)goto 11192                                               880
      if(nin.gt.0) ao(1:nin,m)=a(ia(1:nin))                                 880
      kin(m)=nin                                                            881
      rsqo(m)=rsq                                                           881
      almo(m)=alm                                                           881
      lmu=m                                                                 882
      if(m.lt.mnl)goto 11191                                                882
      if(flmin.ge.1.0)goto 11191                                            883
      me=0                                                                  883
11320 do 11321 j=1,nin                                                      883
      if(ao(j,m).ne.0.0) me=me+1                                            883
11321 continue                                                              883
11322 continue                                                              883
      if(me.gt.ne)goto 11192                                                884
      if(rsq-rsq0.lt.sml*rsq)goto 11192                                     884
      if(rsq.gt.rsqmax)goto 11192                                           885
11191 continue                                                              886
11192 continue                                                              886
      deallocate(a,mm)                                                      887
      return                                                                888
      end                                                                   889
      subroutine spchkvars(no,ni,x,ix,ju)                                   890
      real x(*)                                                             890
      integer ix(*),ju(ni)                                                  891
11330 do 11331 j=1,ni                                                       891
      ju(j)=0                                                               891
      jb=ix(j)                                                              891
      nj=ix(j+1)-jb                                                         891
      if(nj.eq.0)goto 11331                                                 892
      je=ix(j+1)-1                                                          893
      if(nj .ge. no)goto 11351                                              893
11360 do 11361 i=jb,je                                                      893
      if(x(i).eq.0.0)goto 11361                                             893
      ju(j)=1                                                               893
      goto 11362                                                            893
11361 continue                                                              893
11362 continue                                                              893
      goto 11371                                                            894
11351 continue                                                              894
      t=x(jb)                                                               894
11380 do 11381 i=jb+1,je                                                    894
      if(x(i).eq.t)goto 11381                                               894
      ju(j)=1                                                               894
      goto 11382                                                            894
11381 continue                                                              894
11382 continue                                                              894
11371 continue                                                              895
11341 continue                                                              895
11331 continue                                                              896
11332 continue                                                              896
      return                                                                897
      end                                                                   898
      subroutine cmodval(a0,ca,ia,nin,x,ix,jx,n,f)                          899
      real ca(*),x(*),f(n)                                                  899
      integer ia(*),ix(*),jx(*)                                             900
      f=a0                                                                  901
11390 do 11391 j=1,nin                                                      901
      k=ia(j)                                                               901
      kb=ix(k)                                                              901
      ke=ix(k+1)-1                                                          902
      f(jx(kb:ke))=f(jx(kb:ke))+ca(j)*x(kb:ke)                              903
11391 continue                                                              904
11392 continue                                                              904
      return                                                                905
      end                                                                   906
      function row_prod(i,j,ia,ja,ra,w)                                     907
      integer ia(*),ja(*)                                                   907
      real ra(*),w(*)                                                       908
      row_prod=dot(ra(ia(i)),ra(ia(j)),ja(ia(i)),ja(ia(j)),  ia(i+1)-ia(    910 
     *i),ia(j+1)-ia(j),w)
      return                                                                911
      end                                                                   912
      function dot(x,y,mx,my,nx,ny,w)                                       913
      real x(*),y(*),w(*)                                                   913
      integer mx(*),my(*)                                                   914
      i=1                                                                   914
      j=i                                                                   914
      s=0.0                                                                 915
11400 continue                                                              915
11401 continue                                                              915
11410 continue                                                              916
11411 if(mx(i).ge.my(j))goto 11412                                          916
      i=i+1                                                                 916
      if(i.gt.nx) go to 11420                                               916
      goto 11411                                                            917
11412 continue                                                              917
      if(mx(i).eq.my(j)) go to 11430                                        918
11440 continue                                                              918
11441 if(my(j).ge.mx(i))goto 11442                                          918
      j=j+1                                                                 918
      if(j.gt.ny) go to 11420                                               918
      goto 11441                                                            919
11442 continue                                                              919
      if(mx(i).eq.my(j)) go to 11430                                        919
      goto 11401                                                            920
11430 continue                                                              920
      s=s+w(mx(i))*x(i)*y(j)                                                921
      i=i+1                                                                 921
      if(i.gt.nx)goto 11402                                                 921
      j=j+1                                                                 921
      if(j.gt.ny)goto 11402                                                 922
      goto 11401                                                            923
11402 continue                                                              923
11420 continue                                                              923
      dot=s                                                                 924
      return                                                                925
      end                                                                   926
      subroutine lognet (parm,no,ni,nc,x,y,jd,vp,ne,nx,nlam,flmin,ulam,t    928 
     *hr,  isd,maxit,kopt,lmu,a0,ca,ia,nin,dev,alm,nlp,jerr)
      real x(no,ni),y(no,max(2,nc)),vp(ni),ulam(nlam)                       929
      real ca(nx,nc,nlam),a0(nc,nlam),dev(nlam),alm(nlam)                   930
      integer jd(*),ia(nx),nin(nlam)                                        931
      real, dimension (:), allocatable :: xm,xs,ww,vq                           
      integer, dimension (:), allocatable :: ju                                 
      if(maxval(vp) .gt. 0.0)goto 11461                                     935
      jerr=10000                                                            935
      return                                                                935
11461 continue                                                              936
      allocate(ww(1:no),stat=jerr)                                          937
      allocate(ju(1:ni),stat=ierr)                                          937
      jerr=jerr+ierr                                                        938
      allocate(vq(1:ni),stat=ierr)                                          938
      jerr=jerr+ierr                                                        939
      allocate(xm(1:ni),stat=ierr)                                          939
      jerr=jerr+ierr                                                        940
      if(isd .le. 0)goto 11481                                              940
      allocate(xs(1:ni),stat=ierr)                                          940
      jerr=jerr+ierr                                                        940
11481 continue                                                              941
      if(jerr.ne.0) return                                                  942
      call chkvars(no,ni,x,ju)                                              943
      if(jd(1).gt.0) ju(jd(2:(jd(1)+1)))=0                                  944
      if(maxval(ju) .gt. 0)goto 11501                                       944
      jerr=7777                                                             944
      return                                                                944
11501 continue                                                              945
      vq=max(0.0,vp)                                                        945
      vq=vq*ni/sum(vq)                                                      946
11510 do 11511 i=1,no                                                       946
      ww(i)=sum(y(i,:))                                                     946
      y(i,:)=y(i,:)/ww(i)                                                   946
11511 continue                                                              946
11512 continue                                                              946
      ww=ww/sum(ww)                                                         947
      call lstandard1(no,ni,x,ww,ju,isd,xm,xs)                              948
      if(nc .ne. 1)goto 11531                                               949
      call lognet2n(parm,no,ni,x,y(:,1),ww,ju,vq,ne,nx,nlam,flmin,ulam,t    951 
     *hr,  isd,maxit,kopt,lmu,a0,ca,ia,nin,dev,alm,nlp,jerr)
      goto 11541                                                            952
11531 continue                                                              953
      call lognetn(parm,no,ni,nc,x,y,ww,ju,vq,ne,nx,nlam,flmin,ulam,thr,    955 
     *  isd,maxit,kopt,lmu,a0,ca,ia,nin,dev,alm,nlp,jerr)
11541 continue                                                              956
11521 continue                                                              956
      if(jerr.gt.0) return                                                  957
11550 do 11551 k=1,lmu                                                      957
      nk=nin(k)                                                             958
11560 do 11561 ic=1,nc                                                      958
      if(isd .le. 0)goto 11581                                              958
11590 do 11591 l=1,nk                                                       958
      ca(l,ic,k)=ca(l,ic,k)/xs(ia(l))                                       958
11591 continue                                                              958
11592 continue                                                              958
11581 continue                                                              959
      a0(ic,k)=a0(ic,k)-dot_product(ca(1:nk,ic,k),xm(ia(1:nk)))             960
11561 continue                                                              961
11562 continue                                                              961
11551 continue                                                              962
11552 continue                                                              962
      deallocate(ww,ju,vq,xm)                                               962
      if(isd.gt.0) deallocate(xs)                                           963
      return                                                                964
      end                                                                   965
      subroutine lstandard1 (no,ni,x,w,ju,isd,xm,xs)                        966
      real x(no,ni),w(no),xm(ni),xs(ni)                                     966
      integer ju(ni)                                                        967
11600 do 11601 j=1,ni                                                       967
      if(ju(j).eq.0)goto 11601                                              968
      xm(j)=dot_product(w,x(:,j))                                           968
      x(1:no,j)=x(1:no,j)-xm(j)                                             969
      if(isd .le. 0)goto 11621                                              970
      xs(j)=sqrt(dot_product(w*x(:,j),x(:,j)))                              971
      x(1:no,j)=x(1:no,j)/xs(j)                                             972
11621 continue                                                              973
11601 continue                                                              974
11602 continue                                                              974
      return                                                                975
      end                                                                   976
      subroutine lognet2n(parm,no,ni,x,y,w,ju,vp,ne,nx,nlam,flmin,ulam,s    978 
     *hr,  isd,maxit,kopt,lmu,a0,a,m,kin,dev,alm,nlp,jerr)
      parameter(sml=1.0e-5, eps=1.0e-6, pmin=1.0e-5,  big=9.9e30, mnlam=    980 
     *5, devmax=0.999)
      real x(no,ni),y(no),w(no),vp(ni),ulam(nlam)                           981
      real a(nx,nlam),a0(nlam),dev(nlam),alm(nlam)                          982
      integer ju(ni),m(nx),kin(nlam)                                        983
      real, dimension (:), allocatable :: b,bs,v,r,xv,q                         
      integer, dimension (:), allocatable :: mm                                 
      allocate(b(0:ni),stat=jerr)                                           988
      allocate(xv(1:ni),stat=ierr)                                          988
      jerr=jerr+ierr                                                        989
      allocate(bs(0:ni),stat=ierr)                                          989
      jerr=jerr+ierr                                                        990
      allocate(mm(1:ni),stat=ierr)                                          990
      jerr=jerr+ierr                                                        991
      allocate(r(1:no),stat=ierr)                                           991
      jerr=jerr+ierr                                                        992
      allocate(v(1:no),stat=ierr)                                           992
      jerr=jerr+ierr                                                        993
      allocate(q(1:no),stat=ierr)                                           993
      jerr=jerr+ierr                                                        994
      if(jerr.ne.0) return                                                  995
      fmax=log(1.0/pmin-1.0)                                                995
      fmin=-fmax                                                            995
      vmin=(1.0+pmin)*pmin*(1.0-pmin)                                       996
      bta=max(parm,1.0e-3)                                                  996
      omb=1.0-bta                                                           997
      q0=dot_product(w,y)                                                   997
      if(q0 .gt. pmin)goto 11641                                            997
      jerr=8001                                                             997
      return                                                                997
11641 continue                                                              998
      if(q0 .lt. 1.0-pmin)goto 11661                                        998
      jerr=9001                                                             998
      return                                                                998
11661 continue                                                              999
      vi=q0*(1.0-q0)                                                        999
      b(1:ni)=0.0                                                           999
      b(0)=log(q0/(1.0-q0))                                                 999
      v=vi*w                                                               1000
      r=w*(y-q0)                                                           1000
      dev1=-(b(0)*q0+log(1.0-q0))                                          1000
      q=q0                                                                 1001
      if(isd .le. 0)goto 11681                                             1001
      xv=0.25                                                              1001
      goto 11691                                                           1002
11681 continue                                                             1002
11700 do 11701 j=1,ni                                                      1002
      if(ju(j).ne.0) xv(j)=0.25*dot_product(w,x(:,j)**2)                   1002
11701 continue                                                             1002
11702 continue                                                             1002
11691 continue                                                             1003
11671 continue                                                             1003
      xmz=vi                                                               1003
      bs=0.0                                                               1004
      if(flmin .ge. 1.0)goto 11721                                         1004
      eqs=max(eps,flmin)                                                   1004
      alf=eqs**(1.0/(nlam-1))                                              1004
11721 continue                                                             1005
      m=0                                                                  1005
      mm=0                                                                 1005
      nlp=0                                                                1005
      nin=nlp                                                              1005
      mnl=min(mnlam,nlam)                                                  1006
11730 do 11731 ilm=1,nlam                                                  1007
      if(flmin .lt. 1.0)goto 11751                                         1007
      al=ulam(ilm)                                                         1007
      goto 11741                                                           1008
11751 if(ilm .le. 2)goto 11761                                             1008
      al=al*alf                                                            1008
      goto 11741                                                           1009
11761 if(ilm .ne. 1)goto 11771                                             1009
      al=big                                                               1009
      goto 11781                                                           1010
11771 continue                                                             1010
      al=0.0                                                               1011
11790 do 11791 j=1,ni                                                      1011
      if(ju(j).eq.0)goto 11791                                             1011
      if(vp(j).le.0.0)goto 11791                                           1012
      al=max(al,abs(dot_product(r,x(:,j)))/vp(j))                          1013
11791 continue                                                             1014
11792 continue                                                             1014
      al=alf*al/bta                                                        1015
11781 continue                                                             1016
11741 continue                                                             1016
      al2=al*omb                                                           1016
      al1=al*bta                                                           1016
      nit=0                                                                1017
11800 continue                                                             1017
11801 continue                                                             1017
      bs(0)=b(0)                                                           1017
      if(nin.gt.0) bs(m(1:nin))=b(m(1:nin))                                1018
11810 continue                                                             1018
11811 continue                                                             1018
      nlp=nlp+1                                                            1018
      dlx=0.0                                                              1019
11820 do 11821 k=1,ni                                                      1019
      if(ju(k).eq.0)goto 11821                                             1020
      bk=b(k)                                                              1020
      gk=dot_product(r,x(:,k))                                             1021
      u=gk+xv(k)*b(k)                                                      1021
      au=abs(u)-vp(k)*al1                                                  1022
      if(au .gt. 0.0)goto 11841                                            1022
      b(k)=0.0                                                             1022
      goto 11851                                                           1023
11841 continue                                                             1023
      b(k)=sign(au,u)/(xv(k)+vp(k)*al2)                                    1023
11851 continue                                                             1024
11831 continue                                                             1024
      d=b(k)-bk                                                            1024
      if(abs(d).le.0.0)goto 11821                                          1024
      dlx=max(dlx,abs(d))                                                  1025
      r=r-d*v*x(:,k)                                                       1026
      if(mm(k) .ne. 0)goto 11871                                           1026
      nin=nin+1                                                            1026
      if(nin.gt.nx)goto 11822                                              1027
      mm(k)=nin                                                            1027
      m(nin)=k                                                             1028
11871 continue                                                             1029
11821 continue                                                             1030
11822 continue                                                             1030
      if(nin.gt.nx)goto 11812                                              1031
      d=sum(r)/xmz                                                         1032
      if(d .eq. 0.0)goto 11891                                             1032
      b(0)=b(0)+d                                                          1032
      dlx=max(dlx,abs(d))                                                  1032
      r=r-d*v                                                              1032
11891 continue                                                             1033
      if(dlx.lt.shr)goto 11812                                             1034
11900 continue                                                             1034
11901 continue                                                             1034
      nlp=nlp+1                                                            1034
      dlx=0.0                                                              1035
11910 do 11911 l=1,nin                                                     1035
      k=m(l)                                                               1035
      bk=b(k)                                                              1035
      gk=dot_product(r,x(:,k))                                             1036
      u=gk+xv(k)*b(k)                                                      1036
      au=abs(u)-vp(k)*al1                                                  1037
      if(au .gt. 0.0)goto 11931                                            1037
      b(k)=0.0                                                             1037
      goto 11941                                                           1038
11931 continue                                                             1038
      b(k)=sign(au,u)/(xv(k)+vp(k)*al2)                                    1038
11941 continue                                                             1039
11921 continue                                                             1039
      d=b(k)-bk                                                            1039
      if(abs(d).le.0.0)goto 11911                                          1039
      dlx=max(dlx,abs(d))                                                  1040
      r=r-d*v*x(:,k)                                                       1041
11911 continue                                                             1042
11912 continue                                                             1042
      d=sum(r)/xmz                                                         1043
      if(d .eq. 0.0)goto 11961                                             1043
      b(0)=b(0)+d                                                          1043
      dlx=max(dlx,abs(d))                                                  1043
      r=r-d*v                                                              1043
11961 continue                                                             1045
      if(dlx.lt.shr)goto 11902                                             1045
      goto 11901                                                           1046
11902 continue                                                             1046
      goto 11811                                                           1047
11812 continue                                                             1047
      if(nin.gt.nx)goto 11802                                              1048
      if(abs(b(0)-bs(0)) .ge. shr)goto 11981                               1048
      ix=0                                                                 1049
11990 do 11991 j=1,nin                                                     1049
      if(abs(b(m(j))-bs(m(j))).lt.shr)goto 11991                           1049
      ix=1                                                                 1049
      goto 11992                                                           1049
11991 continue                                                             1050
11992 continue                                                             1050
      if(ix.eq.0)goto 11802                                                1051
11981 continue                                                             1052
12000 do 12001 i=1,no                                                      1052
      fi=b(0)                                                              1053
      if(nin.gt.0) fi=fi+dot_product(b(m(1:nin)),x(i,m(1:nin)))            1054
      if(fi .ge. fmin)goto 12021                                           1054
      q(i)=0.0                                                             1054
      goto 12011                                                           1054
12021 if(fi .le. fmax)goto 12031                                           1054
      q(i)=1.0                                                             1054
      goto 12041                                                           1055
12031 continue                                                             1055
      q(i)=1.0/(1.0+exp(-fi))                                              1055
12041 continue                                                             1056
12011 continue                                                             1056
12001 continue                                                             1057
12002 continue                                                             1057
      v=w*q*(1.0-q)                                                        1057
      xmz=sum(v)                                                           1057
      if(xmz.le.vmin)goto 11802                                            1057
      r=w*(y-q)                                                            1058
      if(kopt .ne. 0)goto 12061                                            1059
12070 do 12071 j=1,nin                                                     1059
      xv(m(j))=dot_product(v,x(:,m(j))**2)                                 1059
12071 continue                                                             1060
12072 continue                                                             1060
12061 continue                                                             1061
      nit=nit+1                                                            1061
      if(nit .le. maxit)goto 12091                                         1061
      jerr=-ilm                                                            1061
      return                                                               1061
12091 continue                                                             1062
      goto 11801                                                           1063
11802 continue                                                             1063
      if(nin .le. nx)goto 12111                                            1063
      jerr=-10000-ilm                                                      1063
      goto 11732                                                           1063
12111 continue                                                             1064
      if(nin.gt.0) a(1:nin,ilm)=b(m(1:nin))                                1064
      kin(ilm)=nin                                                         1065
      a0(ilm)=b(0)                                                         1065
      alm(ilm)=al                                                          1065
      lmu=ilm                                                              1066
      devi=dev2(no,w,y,q,pmin)                                             1067
      dev(ilm)=(dev1-devi)/dev1                                            1067
      if(xmz.le.vmin)goto 11732                                            1068
      if(ilm.lt.mnl)goto 11731                                             1068
      if(flmin.ge.1.0)goto 11731                                           1069
      me=0                                                                 1069
12120 do 12121 j=1,nin                                                     1069
      if(a(j,ilm).ne.0.0) me=me+1                                          1069
12121 continue                                                             1069
12122 continue                                                             1069
      if(me.gt.ne)goto 11732                                               1070
      if(dev(ilm).gt.devmax)goto 11732                                     1070
      if(dev(ilm)-dev(ilm-1).lt.sml)goto 11732                             1071
11731 continue                                                             1072
11732 continue                                                             1072
      deallocate(b,bs,v,r,xv,q,mm)                                         1073
      return                                                               1074
      end                                                                  1075
      function dev2(n,w,y,p,pmin)                                          1076
      real w(n),y(n),p(n)                                                  1077
      pmax=1.0-pmin                                                        1077
      s=0.0                                                                1078
12130 do 12131 i=1,n                                                       1078
      pi=min(max(pmin,p(i)),pmax)                                          1079
      s=s-w(i)*(y(i)*log(pi)+(1.0-y(i))*log(1.0-pi))                       1080
12131 continue                                                             1081
12132 continue                                                             1081
      dev2=s                                                               1082
      return                                                               1083
      end                                                                  1084
      subroutine lognetn(parm,no,ni,nc,x,y,w,ju,vp,ne,nx,nlam,flmin,ulam   1086 
     *,shr,  isd,maxit,kopt,lmu,a0,a,m,kin,dev,alm,nlp,jerr)
      parameter(sml=1.0e-5, eps=1.0e-6, pmin=1.0e-5,  big=9.9e30, mnlam=   1088 
     *5, devmax=0.999, exmx=250.0, exmn=-exmx)
      real x(no,ni),y(no,nc),w(no),vp(ni),ulam(nlam)                       1089
      real a(nx,nc,nlam),a0(nc,nlam),dev(nlam),alm(nlam)                   1090
      integer ju(ni),m(nx),kin(nlam)                                       1091
      double precision, dimension (:,:), allocatable :: q                       
      double precision, dimension (:), allocatable :: sxp                       
      real, dimension (:), allocatable :: di,v,r                                
      real, dimension (:,:), allocatable :: b,bs,xv                             
      integer, dimension (:), allocatable :: mm,is                              
      allocate(b(0:ni,1:nc),stat=jerr)                                          
      allocate(xv(1:ni,1:nc),stat=ierr); jerr=jerr+ierr                         
      allocate(bs(0:ni,1:nc),stat=ierr); jerr=jerr+ierr                         
      allocate(q(1:no,1:nc),stat=ierr); jerr=jerr+ierr                          
      allocate(r(1:no),stat=ierr)                                          1102
      jerr=jerr+ierr                                                       1103
      allocate(v(1:no),stat=ierr)                                          1103
      jerr=jerr+ierr                                                       1104
      allocate(mm(1:ni),stat=ierr)                                         1104
      jerr=jerr+ierr                                                       1105
      allocate(is(1:max(nc,ni)),stat=ierr)                                 1105
      jerr=jerr+ierr                                                       1106
      allocate(sxp(1:no),stat=ierr)                                        1106
      jerr=jerr+ierr                                                       1107
      allocate(di(1:no),stat=ierr)                                         1107
      jerr=jerr+ierr                                                       1108
      if(jerr.ne.0) return                                                 1109
      pmax=1.0-pmin                                                        1109
      emin=pmin/pmax                                                       1109
      emax=1.0/emin                                                        1110
      pfm=(1.0+pmin)*pmin                                                  1110
      pfx=(1.0-pmin)*pmax                                                  1110
      vmin=pfm*pmax                                                        1111
      bta=max(parm,1.0e-3)                                                 1111
      omb=1.0-bta                                                          1111
      dev1=0.0                                                             1112
12140 do 12141 ic=1,nc                                                     1112
      q0=dot_product(w,y(:,ic))                                            1113
      if(q0 .gt. pmin)goto 12161                                           1113
      jerr =8000+ic                                                        1113
      return                                                               1113
12161 continue                                                             1114
      if(q0 .lt. 1.0-pmin)goto 12181                                       1114
      jerr =9000+ic                                                        1114
      return                                                               1114
12181 continue                                                             1115
      vi=q0*(1.0-q0)                                                       1115
      v=vi*w                                                               1115
      b(1:ni,ic)=0.0                                                       1116
      b(0,ic)=log(q0)                                                      1116
      dev1=dev1-q0*b(0,ic)                                                 1117
12141 continue                                                             1118
12142 continue                                                             1118
      if(isd .le. 0)goto 12201                                             1118
      xv=0.25                                                              1118
      goto 12211                                                           1119
12201 continue                                                             1119
12220 do 12221 j=1,ni                                                      1119
      if(ju(j).ne.0) xv(j,:)=0.25*dot_product(w,x(:,j)**2)                 1119
12221 continue                                                             1119
12222 continue                                                             1119
12211 continue                                                             1120
12191 continue                                                             1120
      b(0,:)=b(0,:)-sum(b(0,:))/nc                                         1120
      sxp=0.0                                                              1121
12230 do 12231 ic=1,nc                                                     1121
      q(:,ic)=exp(b(0,ic))                                                 1121
      sxp=sxp+q(:,ic)                                                      1121
12231 continue                                                             1122
12232 continue                                                             1122
      if(flmin .ge. 1.0)goto 12251                                         1122
      eqs=max(eps,flmin)                                                   1122
      alf=eqs**(1.0/(nlam-1))                                              1122
12251 continue                                                             1123
      m=0                                                                  1123
      mm=0                                                                 1123
      nin=0                                                                1123
      nlp=0                                                                1123
      mnl=min(mnlam,nlam)                                                  1123
      bs=0.0                                                               1124
12260 do 12261 ilm=1,nlam                                                  1125
      if(flmin .lt. 1.0)goto 12281                                         1125
      al=ulam(ilm)                                                         1125
      goto 12271                                                           1126
12281 if(ilm .le. 2)goto 12291                                             1126
      al=al*alf                                                            1126
      goto 12271                                                           1127
12291 if(ilm .ne. 1)goto 12301                                             1127
      al=big                                                               1127
      goto 12311                                                           1128
12301 continue                                                             1128
      al=0.0                                                               1129
12320 do 12321 ic=1,nc                                                     1129
      r=w*(y(:,ic)-q(:,ic)/sxp)                                            1130
12330 do 12331 j=1,ni                                                      1130
      if(ju(j).eq.0)goto 12331                                             1130
      if(vp(j).le.0.0)goto 12331                                           1131
      al=max(al,abs(dot_product(r,x(:,j)))/vp(j))                          1132
12331 continue                                                             1133
12332 continue                                                             1133
12321 continue                                                             1134
12322 continue                                                             1134
      al=alf*al/bta                                                        1135
12311 continue                                                             1136
12271 continue                                                             1136
      al2=al*omb                                                           1136
      al1=al*bta                                                           1136
      nit=0                                                                1137
12340 continue                                                             1137
12341 continue                                                             1137
      ix=0                                                                 1137
      jx=ix                                                                1137
      ig=0                                                                 1138
12350 do 12351 ic=1,nc                                                     1138
      bs(0,ic)=b(0,ic)                                                     1139
      if(nin.gt.0) bs(m(1:nin),ic)=b(m(1:nin),ic)                          1140
      xmz=0.0                                                              1141
12360 do 12361 i=1,no                                                      1141
      pic=q(i,ic)/sxp(i)                                                   1142
      if(pic .ge. pfm)goto 12381                                           1142
      pic=0.0                                                              1142
      v(i)=0.0                                                             1142
      goto 12371                                                           1143
12381 if(pic .le. pfx)goto 12391                                           1143
      pic=1.0                                                              1143
      v(i)=0.0                                                             1143
      goto 12401                                                           1144
12391 continue                                                             1144
      v(i)=w(i)*pic*(1.0-pic)                                              1144
      xmz=xmz+v(i)                                                         1144
12401 continue                                                             1145
12371 continue                                                             1145
      r(i)=w(i)*(y(i,ic)-pic)                                              1146
12361 continue                                                             1147
12362 continue                                                             1147
      if(xmz.le.vmin)goto 12351                                            1147
      ig=1                                                                 1148
      if(kopt .ne. 0)goto 12421                                            1149
12430 do 12431 j=1,nin                                                     1149
      xv(m(j),ic)=dot_product(v,x(:,m(j))**2)                              1149
12431 continue                                                             1150
12432 continue                                                             1150
12421 continue                                                             1151
12440 continue                                                             1151
12441 continue                                                             1151
      nlp=nlp+1                                                            1151
      dlx=0.0                                                              1152
12450 do 12451 k=1,ni                                                      1152
      if(ju(k).eq.0)goto 12451                                             1153
      bk=b(k,ic)                                                           1153
      gk=dot_product(r,x(:,k))                                             1154
      u=gk+xv(k,ic)*b(k,ic)                                                1154
      au=abs(u)-vp(k)*al1                                                  1155
      if(au .gt. 0.0)goto 12471                                            1155
      b(k,ic)=0.0                                                          1155
      goto 12481                                                           1156
12471 continue                                                             1156
      b(k,ic)=sign(au,u)/(xv(k,ic)+vp(k)*al2)                              1156
12481 continue                                                             1157
12461 continue                                                             1157
      d=b(k,ic)-bk                                                         1157
      if(abs(d).le.0.0)goto 12451                                          1157
      dlx=max(dlx,abs(d))                                                  1158
      r=r-d*v*x(:,k)                                                       1159
      if(mm(k) .ne. 0)goto 12501                                           1159
      nin=nin+1                                                            1160
      if(nin .le. nx)goto 12521                                            1160
      jx=1                                                                 1160
      goto 12452                                                           1160
12521 continue                                                             1161
      mm(k)=nin                                                            1161
      m(nin)=k                                                             1162
12501 continue                                                             1163
12451 continue                                                             1164
12452 continue                                                             1164
      if(jx.gt.0)goto 12442                                                1165
      d=sum(r)/xmz                                                         1166
      if(d .eq. 0.0)goto 12541                                             1166
      b(0,ic)=b(0,ic)+d                                                    1166
      dlx=max(dlx,abs(d))                                                  1166
      r=r-d*v                                                              1166
12541 continue                                                             1167
      if(dlx.lt.shr)goto 12442                                             1168
12550 continue                                                             1168
12551 continue                                                             1168
      nlp=nlp+1                                                            1168
      dlx=0.0                                                              1169
12560 do 12561 l=1,nin                                                     1169
      k=m(l)                                                               1169
      bk=b(k,ic)                                                           1170
      gk=dot_product(r,x(:,k))                                             1171
      u=gk+xv(k,ic)*b(k,ic)                                                1171
      au=abs(u)-vp(k)*al1                                                  1172
      if(au .gt. 0.0)goto 12581                                            1172
      b(k,ic)=0.0                                                          1172
      goto 12591                                                           1173
12581 continue                                                             1173
      b(k,ic)=sign(au,u)/(xv(k,ic)+vp(k)*al2)                              1173
12591 continue                                                             1174
12571 continue                                                             1174
      d=b(k,ic)-bk                                                         1174
      if(abs(d).le.0.0)goto 12561                                          1175
      dlx=max(dlx,abs(d))                                                  1175
      r=r-d*v*x(:,k)                                                       1176
12561 continue                                                             1177
12562 continue                                                             1177
      d=sum(r)/xmz                                                         1178
      if(d .eq. 0.0)goto 12611                                             1178
      b(0,ic)=b(0,ic)+d                                                    1179
      dlx=max(dlx,abs(d))                                                  1179
      r=r-d*v                                                              1180
12611 continue                                                             1181
      if(dlx.lt.shr)goto 12552                                             1181
      goto 12551                                                           1182
12552 continue                                                             1182
      goto 12441                                                           1183
12442 continue                                                             1183
      if(jx.gt.0)goto 12352                                                1184
      if(abs(b(0,ic)-bs(0,ic)).gt.shr) ix=1                                1185
      if(ix .ne. 0)goto 12631                                              1186
12640 do 12641 j=1,nin                                                     1187
      if(abs(b(m(j),ic)-bs(m(j),ic)) .le. shr)goto 12661                   1187
      ix=1                                                                 1187
      goto 12642                                                           1187
12661 continue                                                             1188
12641 continue                                                             1189
12642 continue                                                             1189
12631 continue                                                             1190
12670 do 12671 i=1,no                                                      1190
      fi=b(0,ic)                                                           1192
      if(nin.gt.0) fi=fi+dot_product(b(m(1:nin),ic),x(i,m(1:nin)))         1193
      fi=min(max(exmn,fi),exmx)                                            1193
      sxp(i)=sxp(i)-q(i,ic)                                                1194
      q(i,ic)=min(max(emin*sxp(i),exp(dble(fi))),emax*sxp(i))              1195
      sxp(i)=sxp(i)+q(i,ic)                                                1196
12671 continue                                                             1197
12672 continue                                                             1197
12351 continue                                                             1198
12352 continue                                                             1198
      if(jx.gt.0)goto 12342                                                1198
      if(ix.eq.0)goto 12342                                                1198
      if(ig.eq.0)goto 12342                                                1199
      s=-sum(b(0,:))/nc                                                    1199
      b(0,:)=b(0,:)+s                                                      1199
      di=s                                                                 1200
12680 do 12681 j=1,nin                                                     1200
      l=m(j)                                                               1201
      if(vp(l) .gt. 0.0)goto 12701                                         1201
      s=sum(b(l,:))/nc                                                     1201
      goto 12711                                                           1202
12701 continue                                                             1202
      s=elc(parm,nc,b(l,:),is)                                             1202
12711 continue                                                             1203
12691 continue                                                             1203
      b(l,:)=b(l,:)-s                                                      1203
      di=di-s*x(:,l)                                                       1204
12681 continue                                                             1205
12682 continue                                                             1205
      di=exp(di)                                                           1205
      sxp=sxp*di                                                           1205
12720 do 12721 ic=1,nc                                                     1205
      q(:,ic)=q(:,ic)*di                                                   1205
12721 continue                                                             1206
12722 continue                                                             1206
      nit=nit+1                                                            1206
      if(nit .le. maxit)goto 12741                                         1206
      jerr=-ilm                                                            1206
      return                                                               1206
12741 continue                                                             1207
      goto 12341                                                           1208
12342 continue                                                             1208
      if(jx .le. 0)goto 12761                                              1208
      jerr=-10000-ilm                                                      1208
      goto 12262                                                           1208
12761 continue                                                             1208
      devi=0.0                                                             1209
12770 do 12771 ic=1,nc                                                     1210
      if(nin.gt.0) a(1:nin,ic,ilm)=b(m(1:nin),ic)                          1210
      a0(ic,ilm)=b(0,ic)                                                   1211
12780 do 12781 i=1,no                                                      1211
      if(y(i,ic).le.0.0)goto 12781                                         1212
      devi=devi-w(i)*y(i,ic)*log(q(i,ic)/sxp(i))                           1213
12781 continue                                                             1214
12782 continue                                                             1214
12771 continue                                                             1215
12772 continue                                                             1215
      kin(ilm)=nin                                                         1215
      alm(ilm)=al                                                          1215
      lmu=ilm                                                              1216
      dev(ilm)=(dev1-devi)/dev1                                            1216
      if(ig.eq.0)goto 12262                                                1217
      if(ilm.lt.mnl)goto 12261                                             1217
      if(flmin.ge.1.0)goto 12261                                           1218
      if(nintot(ni,nx,nc,a(1,1,ilm),m,nin,is).gt.ne)goto 12262             1219
      if(dev(ilm).gt.devmax)goto 12262                                     1219
      if(dev(ilm)-dev(ilm-1).lt.sml)goto 12262                             1220
12261 continue                                                             1221
12262 continue                                                             1221
      deallocate(sxp,b,bs,v,r,xv,q,mm,is)                                  1222
      return                                                               1223
      end                                                                  1224
      function elc(parm,n,a,m)                                             1225
      real a(n)                                                            1225
      integer m(n)                                                         1226
      fn=n                                                                 1226
      am=sum(a)/fn                                                         1227
      if((parm .ne. 0.0) .and. (n .ne. 2))goto 12801                       1227
      elc=am                                                               1227
      return                                                               1227
12801 continue                                                             1228
12810 do 12811 i=1,n                                                       1228
      m(i)=i                                                               1228
12811 continue                                                             1228
12812 continue                                                             1228
      call psort7(a,m,1,n)                                                 1229
      if(a(m(1)) .ne. a(m(n)))goto 12831                                   1229
      elc=a(1)                                                             1229
      return                                                               1229
12831 continue                                                             1230
      if(mod(n,2) .ne. 1)goto 12851                                        1230
      ad=a(m(n/2+1))                                                       1230
      goto 12861                                                           1231
12851 continue                                                             1231
      ad=0.5*(a(m(n/2+1))+a(m(n/2)))                                       1231
12861 continue                                                             1232
12841 continue                                                             1232
      if(parm .ne. 1.0)goto 12881                                          1232
      elc=ad                                                               1232
      return                                                               1232
12881 continue                                                             1233
      b1=min(am,ad)                                                        1233
      b2=max(am,ad)                                                        1233
      k2=1                                                                 1234
12890 continue                                                             1234
12891 if(a(m(k2)).gt.b1)goto 12892                                         1234
      k2=k2+1                                                              1234
      goto 12891                                                           1234
12892 continue                                                             1234
      k1=k2-1                                                              1235
12900 continue                                                             1235
12901 if(a(m(k2)).ge.b2)goto 12902                                         1235
      k2=k2+1                                                              1235
      goto 12901                                                           1236
12902 continue                                                             1236
      r=parm/((1.0-parm)*fn)                                               1236
      is=0                                                                 1236
      sm=n-2*(k1-1)                                                        1237
12910 do 12911 k=k1,k2-1                                                   1237
      sm=sm-2.0                                                            1237
      s=r*sm+am                                                            1238
      if(s .le. a(m(k)) .or. s .gt. a(m(k+1)))goto 12931                   1238
      is=k                                                                 1238
      goto 12912                                                           1238
12931 continue                                                             1239
12911 continue                                                             1240
12912 continue                                                             1240
      if(is .eq. 0)goto 12951                                              1240
      elc=s                                                                1240
      return                                                               1240
12951 continue                                                             1240
      r2=2.0*r                                                             1240
      s1=a(m(k1))                                                          1240
      am2=2.0*am                                                           1241
      cri=r2*sum(abs(a-s1))+s1*(s1-am2)                                    1241
      elc=s1                                                               1242
12960 do 12961 k=k1+1,k2                                                   1242
      s=a(m(k))                                                            1242
      if(s.eq.s1)goto 12961                                                1243
      c=r2*sum(abs(a-s))+s*(s-am2)                                         1244
      if(c .ge. cri)goto 12981                                             1244
      cri=c                                                                1244
      elc=s                                                                1244
12981 continue                                                             1244
      s1=s                                                                 1245
12961 continue                                                             1246
12962 continue                                                             1246
      return                                                               1247
      end                                                                  1248
      function nintot(ni,nx,nc,a,m,nin,is)                                 1249
      real a(nx,nc)                                                        1249
      integer m(nx),is(ni)                                                 1250
      is=0                                                                 1250
      nintot=0                                                             1251
12990 do 12991 ic=1,nc                                                     1251
13000 do 13001 j=1,nin                                                     1251
      k=m(j)                                                               1251
      if(is(k).ne.0)goto 13001                                             1252
      if(a(j,ic).eq.0.0)goto 13001                                         1252
      is(k)=k                                                              1252
      nintot=nintot+1                                                      1253
13001 continue                                                             1253
13002 continue                                                             1253
12991 continue                                                             1254
12992 continue                                                             1254
      return                                                               1255
      end                                                                  1256
      subroutine luncomp(ni,nx,nc,ca,ia,nin,a)                             1257
      real ca(nx,nc),a(ni,nc)                                              1257
      integer ia(nx)                                                       1258
      a=0.0                                                                1259
13010 do 13011 ic=1,nc                                                     1259
      if(nin.gt.0) a(ia(1:nin),ic)=ca(1:nin,ic)                            1259
13011 continue                                                             1260
13012 continue                                                             1260
      return                                                               1261
      end                                                                  1262
      subroutine lmodval(nt,x,nc,nx,a0,ca,ia,nin,ans)                      1263
      real a0(nc),ca(nx,nc),x(nt,*),ans(nc,nt)                             1263
      integer ia(nx)                                                       1264
13020 do 13021 i=1,nt                                                      1264
13030 do 13031 ic=1,nc                                                     1264
      ans(ic,i)=a0(ic)                                                     1266
      if(nin.gt.0) ans(ic,i)=ans(ic,i)+dot_product(ca(1:nin,ic),x(i,ia(1   1267 
     *:nin)))
13031 continue                                                             1267
13032 continue                                                             1267
13021 continue                                                             1268
13022 continue                                                             1268
      return                                                               1269
      end                                                                  1270
      subroutine splognet (parm,no,ni,nc,x,ix,jx,y,jd,vp,ne,nx,nlam,flmi   1272 
     *n,  ulam,thr,isd,maxit,kopt,lmu,a0,ca,ia,nin,dev,alm,nlp,jerr)
      real x(*),y(no,max(2,nc)),vp(ni),ulam(nlam)                          1273
      real ca(nx,nc,nlam),a0(nc,nlam),dev(nlam),alm(nlam)                  1274
      integer ix(*),jx(*),jd(*),ia(nx),nin(nlam)                           1275
      real, dimension (:), allocatable :: xm,xs,ww,vq                           
      integer, dimension (:), allocatable :: ju                                 
      if(maxval(vp) .gt. 0.0)goto 13051                                    1279
      jerr=10000                                                           1279
      return                                                               1279
13051 continue                                                             1280
      allocate(ww(1:no),stat=jerr)                                         1281
      allocate(ju(1:ni),stat=ierr)                                         1281
      jerr=jerr+ierr                                                       1282
      allocate(vq(1:ni),stat=ierr)                                         1282
      jerr=jerr+ierr                                                       1283
      allocate(xm(1:ni),stat=ierr)                                         1283
      jerr=jerr+ierr                                                       1284
      allocate(xs(1:ni),stat=ierr)                                         1284
      jerr=jerr+ierr                                                       1285
      if(jerr.ne.0) return                                                 1286
      call spchkvars(no,ni,x,ix,ju)                                        1287
      if(jd(1).gt.0) ju(jd(2:(jd(1)+1)))=0                                 1288
      if(maxval(ju) .gt. 0)goto 13071                                      1288
      jerr=7777                                                            1288
      return                                                               1288
13071 continue                                                             1289
      vq=max(0.0,vp)                                                       1289
      vq=vq*ni/sum(vq)                                                     1290
13080 do 13081 i=1,no                                                      1290
      ww(i)=sum(y(i,:))                                                    1290
      y(i,:)=y(i,:)/ww(i)                                                  1290
13081 continue                                                             1290
13082 continue                                                             1290
      ww=ww/sum(ww)                                                        1291
      call splstandard2(no,ni,x,ix,jx,ww,ju,isd,xm,xs)                     1292
      if(nc .ne. 1)goto 13101                                              1293
      call sprlognet2n(parm,no,ni,x,ix,jx,y(:,1),ww,ju,vq,ne,nx,nlam,flm   1295 
     *in,  ulam,thr,isd,maxit,kopt,xm,xs,lmu,a0,ca,ia,nin,dev,alm,nlp,je
     *rr)
      goto 13111                                                           1296
13101 continue                                                             1297
      call sprlognetn(parm,no,ni,nc,x,ix,jx,y,ww,ju,vq,ne,nx,nlam,flmin,   1299 
     *ulam,thr,  isd,maxit,kopt,xm,xs,lmu,a0,ca,ia,nin,dev,alm,nlp,jerr)
13111 continue                                                             1300
13091 continue                                                             1300
      if(jerr.gt.0) return                                                 1301
13120 do 13121 k=1,lmu                                                     1301
      nk=nin(k)                                                            1302
13130 do 13131 ic=1,nc                                                     1302
      if(isd .le. 0)goto 13151                                             1302
13160 do 13161 l=1,nk                                                      1302
      ca(l,ic,k)=ca(l,ic,k)/xs(ia(l))                                      1302
13161 continue                                                             1302
13162 continue                                                             1302
13151 continue                                                             1303
      a0(ic,k)=a0(ic,k)-dot_product(ca(1:nk,ic,k),xm(ia(1:nk)))            1304
13131 continue                                                             1305
13132 continue                                                             1305
13121 continue                                                             1306
13122 continue                                                             1306
      deallocate(ww,ju,vq,xm,xs)                                           1307
      return                                                               1308
      end                                                                  1309
      subroutine splstandard2(no,ni,x,ix,jx,w,ju,isd,xm,xs)                1310
      real x(*),w(no),xm(ni),xs(ni)                                        1310
      integer ix(*),jx(*),ju(ni)                                           1311
13170 do 13171 j=1,ni                                                      1311
      if(ju(j).eq.0)goto 13171                                             1311
      jb=ix(j)                                                             1311
      je=ix(j+1)-1                                                         1312
      xm(j)=dot_product(w(jx(jb:je)),x(jb:je))                             1313
      if(isd.gt.0) xs(j)=sqrt(dot_product(w(jx(jb:je)),x(jb:je)**2)-xm(j   1314 
     *)**2)
13171 continue                                                             1315
13172 continue                                                             1315
      if(isd.eq.0) xs=1.0                                                  1316
      return                                                               1317
      end                                                                  1318
      subroutine sprlognet2n (parm,no,ni,x,ix,jx,y,w,ju,vp,ne,nx,nlam,     1320 
     *flmin,ulam,shr,isd,maxit,kopt,xb,xs,lmu,a0,a,m,kin,dev,alm,nlp,jer
     *r)
      parameter(sml=1.0e-5, eps=1.0e-6, pmin=1.0e-5,  big=9.9e30, mnlam=   1322 
     *5, devmax=0.999)
      real x(*),y(no),w(no),vp(ni),ulam(nlam)                              1323
      real a(nx,nlam),a0(nlam),dev(nlam),alm(nlam)                         1324
      real xb(ni),xs(ni)                                                   1324
      integer ix(*),jx(*),ju(ni),m(nx),kin(nlam)                           1325
      real, dimension (:), allocatable :: xm,b,bs,v,r,sc,xv,q                   
      integer, dimension (:), allocatable :: mm                                 
      allocate(b(0:ni),stat=jerr)                                          1330
      allocate(xm(0:ni),stat=ierr)                                         1330
      jerr=jerr+ierr                                                       1331
      allocate(xv(1:ni),stat=ierr)                                         1331
      jerr=jerr+ierr                                                       1332
      allocate(bs(0:ni),stat=ierr)                                         1332
      jerr=jerr+ierr                                                       1333
      allocate(mm(1:ni),stat=ierr)                                         1333
      jerr=jerr+ierr                                                       1334
      allocate(q(1:no),stat=ierr)                                          1334
      jerr=jerr+ierr                                                       1335
      allocate(r(1:no),stat=ierr)                                          1335
      jerr=jerr+ierr                                                       1336
      allocate(v(1:no),stat=ierr)                                          1336
      jerr=jerr+ierr                                                       1337
      allocate(sc(1:no),stat=ierr)                                         1337
      jerr=jerr+ierr                                                       1338
      if(jerr.ne.0) return                                                 1339
      fmax=log(1.0/pmin-1.0)                                               1339
      fmin=-fmax                                                           1339
      vmin=(1.0+pmin)*pmin*(1.0-pmin)                                      1340
      bta=max(parm,1.0e-3)                                                 1340
      omb=1.0-bta                                                          1341
      q0=dot_product(w,y)                                                  1341
      if(q0 .gt. pmin)goto 13191                                           1341
      jerr=8001                                                            1341
      return                                                               1341
13191 continue                                                             1342
      if(q0 .lt. 1.0-pmin)goto 13211                                       1342
      jerr=9001                                                            1342
      return                                                               1342
13211 continue                                                             1343
      vi=q0*(1.0-q0)                                                       1343
      b(1:ni)=0.0                                                          1343
      b(0)=log(q0/(1.0-q0))                                                1343
      v=vi*w                                                               1344
      r=w*(y-q0)                                                           1344
      dev1=-(b(0)*q0+log(1.0-q0))                                          1344
      q=q0                                                                 1345
      if(isd .le. 0)goto 13231                                             1345
      xv=0.25                                                              1345
      goto 13241                                                           1346
13231 continue                                                             1347
13250 do 13251 j=1,ni                                                      1347
      if(ju(j).eq.0)goto 13251                                             1347
      jb=ix(j)                                                             1347
      je=ix(j+1)-1                                                         1348
      xv(j)=0.25*(dot_product(w(jx(jb:je)),x(jb:je)**2)-xb(j)**2)          1349
13251 continue                                                             1350
13252 continue                                                             1350
13241 continue                                                             1351
13221 continue                                                             1351
      xm(0)=vi                                                             1351
      nlp=0                                                                1351
      nin=nlp                                                              1352
      if(flmin .ge. 1.0)goto 13271                                         1352
      eqs=max(eps,flmin)                                                   1352
      alf=eqs**(1.0/(nlam-1))                                              1352
13271 continue                                                             1353
      m=0                                                                  1353
      mm=0                                                                 1353
      nin=0                                                                1353
      o=0.0                                                                1353
      svr=o                                                                1353
      mnl=min(mnlam,nlam)                                                  1353
      bs=0.0                                                               1354
13280 do 13281 ilm=1,nlam                                                  1355
      if(flmin .lt. 1.0)goto 13301                                         1355
      al=ulam(ilm)                                                         1355
      goto 13291                                                           1356
13301 if(ilm .le. 2)goto 13311                                             1356
      al=al*alf                                                            1356
      goto 13291                                                           1357
13311 if(ilm .ne. 1)goto 13321                                             1357
      al=big                                                               1357
      goto 13331                                                           1358
13321 continue                                                             1358
      al=0.0                                                               1359
13340 do 13341 j=1,ni                                                      1359
      if(ju(j).eq.0)goto 13341                                             1359
      if(vp(j).le.0.0)goto 13341                                           1360
      jb=ix(j)                                                             1360
      je=ix(j+1)-1                                                         1360
      jn=ix(j+1)-ix(j)                                                     1361
      sc(1:jn)=r(jx(jb:je))+v(jx(jb:je))*o                                 1362
      gj=dot_product(sc(1:jn),x(jb:je))                                    1363
      gj=(gj-svr*xb(j))/xs(j)                                              1364
      al=max(al,abs(gj)/vp(j))                                             1365
13341 continue                                                             1366
13342 continue                                                             1366
      al=alf*al/bta                                                        1367
13331 continue                                                             1368
13291 continue                                                             1368
      al2=al*omb                                                           1368
      al1=al*bta                                                           1368
      nit=0                                                                1369
13350 continue                                                             1369
13351 continue                                                             1369
      bs(0)=b(0)                                                           1369
      if(nin.gt.0) bs(m(1:nin))=b(m(1:nin))                                1370
13360 continue                                                             1370
13361 continue                                                             1370
      nlp=nlp+1                                                            1370
      dlx=0.0                                                              1371
13370 do 13371 k=1,ni                                                      1371
      if(ju(k).eq.0)goto 13371                                             1372
      jb=ix(k)                                                             1372
      je=ix(k+1)-1                                                         1372
      jn=ix(k+1)-ix(k)                                                     1372
      bk=b(k)                                                              1373
      sc(1:jn)=r(jx(jb:je))+v(jx(jb:je))*o                                 1374
      gk=dot_product(sc(1:jn),x(jb:je))                                    1375
      gk=(gk-svr*xb(k))/xs(k)                                              1376
      u=gk+xv(k)*b(k)                                                      1376
      au=abs(u)-vp(k)*al1                                                  1377
      if(au .gt. 0.0)goto 13391                                            1377
      b(k)=0.0                                                             1377
      goto 13401                                                           1378
13391 continue                                                             1378
      b(k)=sign(au,u)/(xv(k)+vp(k)*al2)                                    1378
13401 continue                                                             1379
13381 continue                                                             1379
      d=b(k)-bk                                                            1379
      if(abs(d).le.0.0)goto 13371                                          1379
      dlx=max(dlx,abs(d))                                                  1380
      if(mm(k) .ne. 0)goto 13421                                           1380
      nin=nin+1                                                            1380
      if(nin.gt.nx)goto 13372                                              1381
      mm(k)=nin                                                            1381
      m(nin)=k                                                             1381
      sc(1:jn)=v(jx(jb:je))                                                1382
      xm(k)=dot_product(sc(1:jn),x(jb:je))                                 1383
13421 continue                                                             1384
      r(jx(jb:je))=r(jx(jb:je))-d*v(jx(jb:je))*x(jb:je)/xs(k)              1385
      o=o+d*(xb(k)/xs(k))                                                  1386
      svr=svr-d*(xm(k)-xb(k)*xm(0))/xs(k)                                  1387
13371 continue                                                             1388
13372 continue                                                             1388
      if(nin.gt.nx)goto 13362                                              1389
      d=svr/xm(0)                                                          1390
      if(d .eq. 0.0)goto 13441                                             1390
      b(0)=b(0)+d                                                          1390
      dlx=max(dlx,abs(d))                                                  1390
      r=r-d*v                                                              1390
13441 continue                                                             1391
      svr=svr-d*xm(0)                                                      1391
      if(dlx.lt.shr)goto 13362                                             1392
13450 continue                                                             1392
13451 continue                                                             1392
      nlp=nlp+1                                                            1392
      dlx=0.0                                                              1393
13460 do 13461 l=1,nin                                                     1393
      k=m(l)                                                               1393
      jb=ix(k)                                                             1393
      je=ix(k+1)-1                                                         1394
      jn=ix(k+1)-ix(k)                                                     1394
      bk=b(k)                                                              1395
      sc(1:jn)=r(jx(jb:je))+v(jx(jb:je))*o                                 1396
      gk=dot_product(sc(1:jn),x(jb:je))                                    1397
      gk=(gk-svr*xb(k))/xs(k)                                              1398
      u=gk+xv(k)*b(k)                                                      1398
      au=abs(u)-vp(k)*al1                                                  1399
      if(au .gt. 0.0)goto 13481                                            1399
      b(k)=0.0                                                             1399
      goto 13491                                                           1400
13481 continue                                                             1400
      b(k)=sign(au,u)/(xv(k)+vp(k)*al2)                                    1400
13491 continue                                                             1401
13471 continue                                                             1401
      d=b(k)-bk                                                            1401
      if(abs(d).le.0.0)goto 13461                                          1401
      dlx=max(dlx,abs(d))                                                  1402
      r(jx(jb:je))=r(jx(jb:je))-d*v(jx(jb:je))*x(jb:je)/xs(k)              1403
      o=o+d*(xb(k)/xs(k))                                                  1404
      svr=svr-d*(xm(k)-xb(k)*xm(0))/xs(k)                                  1405
13461 continue                                                             1406
13462 continue                                                             1406
      d=svr/xm(0)                                                          1407
      if(d .eq. 0.0)goto 13511                                             1407
      b(0)=b(0)+d                                                          1407
      dlx=max(dlx,abs(d))                                                  1407
      r=r-d*v                                                              1407
13511 continue                                                             1408
      svr=svr-d*xm(0)                                                      1409
      if(dlx.lt.shr)goto 13452                                             1409
      goto 13451                                                           1410
13452 continue                                                             1410
      goto 13361                                                           1411
13362 continue                                                             1411
      if(nin.gt.nx)goto 13352                                              1412
      if(abs(b(0)-bs(0)) .ge. shr)goto 13531                               1412
      kx=0                                                                 1413
13540 do 13541 j=1,nin                                                     1413
      if(abs(b(m(j))-bs(m(j))).lt.shr)goto 13541                           1413
      kx=1                                                                 1413
      goto 13542                                                           1413
13541 continue                                                             1414
13542 continue                                                             1414
      if(kx.eq.0)goto 13352                                                1415
13531 continue                                                             1416
      sc=b(0)                                                              1416
      b0=0.0                                                               1417
13550 do 13551 j=1,nin                                                     1417
      l=m(j)                                                               1417
      jb=ix(l)                                                             1417
      je=ix(l+1)-1                                                         1418
      sc(jx(jb:je))=sc(jx(jb:je))+b(l)*x(jb:je)/xs(l)                      1419
      b0=b0-b(l)*xb(l)/xs(l)                                               1420
13551 continue                                                             1421
13552 continue                                                             1421
      sc=sc+b0                                                             1422
13560 do 13561 i=1,no                                                      1422
      fi=sc(i)                                                             1423
      if(fi .ge. fmin)goto 13581                                           1423
      q(i)=0.0                                                             1423
      goto 13571                                                           1423
13581 if(fi .le. fmax)goto 13591                                           1423
      q(i)=1.0                                                             1423
      goto 13601                                                           1424
13591 continue                                                             1424
      q(i)=1.0/(1.0+exp(-fi))                                              1424
13601 continue                                                             1425
13571 continue                                                             1425
13561 continue                                                             1426
13562 continue                                                             1426
      v=w*q*(1.0-q)                                                        1426
      r=w*(y-q)                                                            1427
      xm(0)=sum(v)                                                         1427
      if(xm(0).lt.vmin)goto 13352                                          1427
      svr=sum(r)                                                           1427
      o=0.0                                                                1428
13610 do 13611 l=1,nin                                                     1428
      j=m(l)                                                               1429
      jb=ix(j)                                                             1429
      je=ix(j+1)-1                                                         1429
      jn=ix(j+1)-ix(j)                                                     1430
      sc(1:jn)=v(jx(jb:je))                                                1431
      xm(j)=dot_product(sc(1:jn),x(jb:je))                                 1432
      if(kopt .ne. 0)goto 13631                                            1433
      xv(j)=dot_product(sc(1:jn),x(jb:je)**2)                              1434
      xv(j)=(xv(j)-2.0*xb(j)*xm(j)+xm(0)*xb(j)**2)/xs(j)**2                1436
13631 continue                                                             1437
13611 continue                                                             1438
13612 continue                                                             1438
      nit=nit+1                                                            1438
      if(nit .le. maxit)goto 13651                                         1438
      jerr=-ilm                                                            1438
      return                                                               1438
13651 continue                                                             1439
      goto 13351                                                           1440
13352 continue                                                             1440
      if(nin .le. nx)goto 13671                                            1440
      jerr=-10000-ilm                                                      1440
      goto 13282                                                           1440
13671 continue                                                             1441
      if(nin.gt.0) a(1:nin,ilm)=b(m(1:nin))                                1441
      kin(ilm)=nin                                                         1442
      a0(ilm)=b(0)                                                         1442
      alm(ilm)=al                                                          1442
      lmu=ilm                                                              1443
      devi=dev2(no,w,y,q,pmin)                                             1444
      dev(ilm)=(dev1-devi)/dev1                                            1445
      if(ilm.lt.mnl)goto 13281                                             1445
      if(flmin.ge.1.0)goto 13281                                           1446
      me=0                                                                 1446
13680 do 13681 j=1,nin                                                     1446
      if(a(j,ilm).ne.0.0) me=me+1                                          1446
13681 continue                                                             1446
13682 continue                                                             1446
      if(me.gt.ne)goto 13282                                               1447
      if(dev(ilm).gt.devmax)goto 13282                                     1447
      if(dev(ilm)-dev(ilm-1).lt.sml)goto 13282                             1448
      if(xm(0).lt.vmin)goto 13282                                          1449
13281 continue                                                             1450
13282 continue                                                             1450
      deallocate(xm,b,bs,v,r,sc,xv,q,mm)                                   1451
      return                                                               1452
      end                                                                  1453
      subroutine sprlognetn(parm,no,ni,nc,x,ix,jx,y,w,ju,vp,ne,nx,nlam,f   1455 
     *lmin,ulam,  shr,isd,maxit,kopt,xb,xs,lmu,a0,a,m,kin,dev,alm,nlp,je
     *rr)
      parameter(sml=1.0e-5, eps=1.0e-6, pmin=1.0e-5,  big=9.9e30, mnlam=   1457 
     *5, devmax=0.999, exmx=250.0, exmn=-exmx)
      real x(*),y(no,nc),w(no),vp(ni),ulam(nlam),xb(ni),xs(ni)             1458
      real a(nx,nc,nlam),a0(nc,nlam),dev(nlam),alm(nlam)                   1459
      integer ix(*),jx(*),ju(ni),m(nx),kin(nlam)                           1460
      double precision, dimension (:,:), allocatable :: q                       
      double precision, dimension (:), allocatable :: sxp                       
      real, dimension (:), allocatable :: sc,xm,v,r                             
      real, dimension (:,:), allocatable :: b,bs,xv                             
      integer, dimension (:), allocatable :: mm,is                              
      allocate(b(0:ni,1:nc),stat=jerr)                                          
      allocate(xv(1:ni,1:nc),stat=ierr); jerr=jerr+ierr                         
      allocate(bs(0:ni,1:nc),stat=ierr); jerr=jerr+ierr                         
      allocate(q(1:no,1:nc),stat=ierr); jerr=jerr+ierr                          
      allocate(xm(0:ni),stat=ierr)                                         1471
      jerr=jerr+ierr                                                       1472
      allocate(r(1:no),stat=ierr)                                          1472
      jerr=jerr+ierr                                                       1473
      allocate(v(1:no),stat=ierr)                                          1473
      jerr=jerr+ierr                                                       1474
      allocate(mm(1:ni),stat=ierr)                                         1474
      jerr=jerr+ierr                                                       1475
      allocate(is(1:max(nc,ni)),stat=ierr)                                 1475
      jerr=jerr+ierr                                                       1476
      allocate(sxp(1:no),stat=ierr)                                        1476
      jerr=jerr+ierr                                                       1477
      allocate(sc(1:no),stat=ierr)                                         1477
      jerr=jerr+ierr                                                       1478
      if(jerr.ne.0) return                                                 1479
      pmax=1.0-pmin                                                        1479
      emin=pmin/pmax                                                       1479
      emax=1.0/emin                                                        1480
      pfm=(1.0+pmin)*pmin                                                  1480
      pfx=(1.0-pmin)*pmax                                                  1480
      vmin=pfm*pmax                                                        1481
      bta=max(parm,1.0e-3)                                                 1481
      omb=1.0-bta                                                          1481
      dev1=0.0                                                             1482
13690 do 13691 ic=1,nc                                                     1482
      q0=dot_product(w,y(:,ic))                                            1483
      if(q0 .gt. pmin)goto 13711                                           1483
      jerr =8000+ic                                                        1483
      return                                                               1483
13711 continue                                                             1484
      if(q0 .lt. 1.0-pmin)goto 13731                                       1484
      jerr =9000+ic                                                        1484
      return                                                               1484
13731 continue                                                             1485
      vi=q0*(1.0-q0)                                                       1485
      v=vi*w                                                               1485
      b(1:ni,ic)=0.0                                                       1486
      b(0,ic)=log(q0)                                                      1486
      xm(0)=vi                                                             1486
      dev1=dev1-q0*b(0,ic)                                                 1487
13691 continue                                                             1488
13692 continue                                                             1488
      if(isd .le. 0)goto 13751                                             1488
      xv=0.25                                                              1488
      goto 13761                                                           1489
13751 continue                                                             1490
13770 do 13771 j=1,ni                                                      1490
      if(ju(j).eq.0)goto 13771                                             1490
      jb=ix(j)                                                             1490
      je=ix(j+1)-1                                                         1491
      xv(j,:)=0.25*(dot_product(w(jx(jb:je)),x(jb:je)**2)-xb(j)**2)        1492
13771 continue                                                             1493
13772 continue                                                             1493
13761 continue                                                             1494
13741 continue                                                             1494
      b(0,:)=b(0,:)-sum(b(0,:))/nc                                         1494
      sxp=0.0                                                              1495
13780 do 13781 ic=1,nc                                                     1495
      q(:,ic)=exp(b(0,ic))                                                 1495
      sxp=sxp+q(:,ic)                                                      1495
13781 continue                                                             1496
13782 continue                                                             1496
      if(flmin .ge. 1.0)goto 13801                                         1496
      eqs=max(eps,flmin)                                                   1496
      alf=eqs**(1.0/(nlam-1))                                              1496
13801 continue                                                             1497
      m=0                                                                  1497
      mm=0                                                                 1497
      nin=0                                                                1497
      nlp=0                                                                1497
      mnl=min(mnlam,nlam)                                                  1497
      bs=0.0                                                               1497
      svr=0.0                                                              1497
      o=0.0                                                                1498
13810 do 13811 ilm=1,nlam                                                  1499
      if(flmin .lt. 1.0)goto 13831                                         1499
      al=ulam(ilm)                                                         1499
      goto 13821                                                           1500
13831 if(ilm .le. 2)goto 13841                                             1500
      al=al*alf                                                            1500
      goto 13821                                                           1501
13841 if(ilm .ne. 1)goto 13851                                             1501
      al=big                                                               1501
      goto 13861                                                           1502
13851 continue                                                             1502
      al=0.0                                                               1503
13870 do 13871 ic=1,nc                                                     1503
      v=q(:,ic)/sxp                                                        1503
      r=w*(y(:,ic)-v)                                                      1503
      v=w*v*(1.0-v)                                                        1504
13880 do 13881 j=1,ni                                                      1504
      if(ju(j).eq.0)goto 13881                                             1504
      if(vp(j).le.0.0)goto 13881                                           1505
      jb=ix(j)                                                             1505
      je=ix(j+1)-1                                                         1505
      jn=ix(j+1)-ix(j)                                                     1506
      sc(1:jn)=r(jx(jb:je))+o*v(jx(jb:je))                                 1507
      gj=dot_product(sc(1:jn),x(jb:je))                                    1508
      gj=(gj-svr*xb(j))/xs(j)                                              1509
      al=max(al,abs(gj)/vp(j))                                             1510
13881 continue                                                             1511
13882 continue                                                             1511
13871 continue                                                             1512
13872 continue                                                             1512
      al=alf*al/bta                                                        1513
13861 continue                                                             1514
13821 continue                                                             1514
      al2=al*omb                                                           1514
      al1=al*bta                                                           1514
      nit=0                                                                1515
13890 continue                                                             1515
13891 continue                                                             1515
      ixx=0                                                                1515
      jxx=ixx                                                              1515
      ig=0                                                                 1516
13900 do 13901 ic=1,nc                                                     1516
      bs(0,ic)=b(0,ic)                                                     1517
      if(nin.gt.0) bs(m(1:nin),ic)=b(m(1:nin),ic)                          1518
      xm(0)=0.0                                                            1518
      svr=0.0                                                              1518
      o=0.0                                                                1519
13910 do 13911 i=1,no                                                      1519
      pic=q(i,ic)/sxp(i)                                                   1520
      if(pic .ge. pfm)goto 13931                                           1520
      pic=0.0                                                              1520
      v(i)=0.0                                                             1520
      goto 13921                                                           1521
13931 if(pic .le. pfx)goto 13941                                           1521
      pic=1.0                                                              1521
      v(i)=0.0                                                             1521
      goto 13951                                                           1522
13941 continue                                                             1522
      v(i)=w(i)*pic*(1.0-pic)                                              1522
      xm(0)=xm(0)+v(i)                                                     1522
13951 continue                                                             1523
13921 continue                                                             1523
      r(i)=w(i)*(y(i,ic)-pic)                                              1523
      svr=svr+r(i)                                                         1524
13911 continue                                                             1525
13912 continue                                                             1525
      if(xm(0).le.vmin)goto 13901                                          1525
      ig=1                                                                 1526
13960 do 13961 l=1,nin                                                     1526
      j=m(l)                                                               1527
      jb=ix(j)                                                             1527
      je=ix(j+1)-1                                                         1528
      xm(j)=dot_product(v(jx(jb:je)),x(jb:je))                             1529
      if(kopt .ne. 0)goto 13981                                            1530
      xv(j,ic)=dot_product(v(jx(jb:je)),x(jb:je)**2)                       1531
      xv(j,ic)=(xv(j,ic)-2.0*xb(j)*xm(j)+xm(0)*xb(j)**2)/xs(j)**2          1532
13981 continue                                                             1533
13961 continue                                                             1534
13962 continue                                                             1534
13990 continue                                                             1534
13991 continue                                                             1534
      nlp=nlp+1                                                            1534
      dlx=0.0                                                              1535
14000 do 14001 k=1,ni                                                      1535
      if(ju(k).eq.0)goto 14001                                             1536
      jb=ix(k)                                                             1536
      je=ix(k+1)-1                                                         1536
      jn=ix(k+1)-ix(k)                                                     1536
      bk=b(k,ic)                                                           1537
      sc(1:jn)=r(jx(jb:je))+o*v(jx(jb:je))                                 1538
      gk=dot_product(sc(1:jn),x(jb:je))                                    1539
      gk=(gk-svr*xb(k))/xs(k)                                              1540
      u=gk+xv(k,ic)*b(k,ic)                                                1540
      au=abs(u)-vp(k)*al1                                                  1541
      if(au .gt. 0.0)goto 14021                                            1541
      b(k,ic)=0.0                                                          1541
      goto 14031                                                           1542
14021 continue                                                             1542
      b(k,ic)=sign(au,u)/(xv(k,ic)+vp(k)*al2)                              1542
14031 continue                                                             1543
14011 continue                                                             1543
      d=b(k,ic)-bk                                                         1543
      if(abs(d).le.0.0)goto 14001                                          1543
      dlx=max(dlx,abs(d))                                                  1544
      if(mm(k) .ne. 0)goto 14051                                           1544
      nin=nin+1                                                            1545
      if(nin .le. nx)goto 14071                                            1545
      jxx=1                                                                1545
      goto 14002                                                           1545
14071 continue                                                             1546
      mm(k)=nin                                                            1546
      m(nin)=k                                                             1547
      xm(k)=dot_product(v(jx(jb:je)),x(jb:je))                             1548
14051 continue                                                             1549
      r(jx(jb:je))=r(jx(jb:je))-d*v(jx(jb:je))*x(jb:je)/xs(k)              1550
      o=o+d*(xb(k)/xs(k))                                                  1551
      svr=svr-d*(xm(k)-xb(k)*xm(0))/xs(k)                                  1552
14001 continue                                                             1553
14002 continue                                                             1553
      if(jxx.gt.0)goto 13992                                               1554
      d=svr/xm(0)                                                          1555
      if(d .eq. 0.0)goto 14091                                             1555
      b(0,ic)=b(0,ic)+d                                                    1555
      dlx=max(dlx,abs(d))                                                  1556
      r=r-d*v                                                              1556
      svr=svr-d*xm(0)                                                      1557
14091 continue                                                             1558
      if(dlx.lt.shr)goto 13992                                             1559
14100 continue                                                             1559
14101 continue                                                             1559
      nlp=nlp+1                                                            1559
      dlx=0.0                                                              1560
14110 do 14111 l=1,nin                                                     1560
      k=m(l)                                                               1560
      jb=ix(k)                                                             1560
      je=ix(k+1)-1                                                         1561
      jn=ix(k+1)-ix(k)                                                     1561
      bk=b(k,ic)                                                           1562
      sc(1:jn)=r(jx(jb:je))+o*v(jx(jb:je))                                 1563
      gk=dot_product(sc(1:jn),x(jb:je))                                    1564
      gk=(gk-svr*xb(k))/xs(k)                                              1565
      u=gk+xv(k,ic)*b(k,ic)                                                1565
      au=abs(u)-vp(k)*al1                                                  1566
      if(au .gt. 0.0)goto 14131                                            1566
      b(k,ic)=0.0                                                          1566
      goto 14141                                                           1567
14131 continue                                                             1567
      b(k,ic)=sign(au,u)/(xv(k,ic)+vp(k)*al2)                              1567
14141 continue                                                             1568
14121 continue                                                             1568
      d=b(k,ic)-bk                                                         1568
      if(abs(d).le.0.0)goto 14111                                          1569
      dlx=max(dlx,abs(d))                                                  1570
      r(jx(jb:je))=r(jx(jb:je))-d*v(jx(jb:je))*x(jb:je)/xs(k)              1571
      o=o+d*(xb(k)/xs(k))                                                  1572
      svr=svr-d*(xm(k)-xb(k)*xm(0))/xs(k)                                  1573
14111 continue                                                             1574
14112 continue                                                             1574
      d=svr/xm(0)                                                          1575
      if(d .eq. 0.0)goto 14161                                             1575
      b(0,ic)=b(0,ic)+d                                                    1575
      dlx=max(dlx,abs(d))                                                  1576
      r=r-d*v                                                              1576
      svr=svr-d*xm(0)                                                      1577
14161 continue                                                             1578
      if(dlx.lt.shr)goto 14102                                             1578
      goto 14101                                                           1579
14102 continue                                                             1579
      goto 13991                                                           1580
13992 continue                                                             1580
      if(jxx.gt.0)goto 13902                                               1581
      if(abs(b(0,ic)-bs(0,ic)).gt.shr) ixx=1                               1582
      if(ixx .ne. 0)goto 14181                                             1583
14190 do 14191 j=1,nin                                                     1584
      if(abs(b(m(j),ic)-bs(m(j),ic)) .le. shr)goto 14211                   1584
      ixx=1                                                                1584
      goto 14192                                                           1584
14211 continue                                                             1585
14191 continue                                                             1586
14192 continue                                                             1586
14181 continue                                                             1587
      sc=b(0,ic)                                                           1587
      b0=0.0                                                               1588
14220 do 14221 j=1,nin                                                     1588
      l=m(j)                                                               1588
      jb=ix(l)                                                             1588
      je=ix(l+1)-1                                                         1589
      sc(jx(jb:je))=sc(jx(jb:je))+b(l,ic)*x(jb:je)/xs(l)                   1590
      b0=b0-b(l,ic)*xb(l)/xs(l)                                            1591
14221 continue                                                             1592
14222 continue                                                             1592
      sc=min(max(exmn,sc+b0),exmx)                                         1593
      sxp=sxp-q(:,ic)                                                      1594
      q(:,ic)=min(max(emin*sxp,exp(dble(sc))),emax*sxp)                    1595
      sxp=sxp+q(:,ic)                                                      1596
13901 continue                                                             1597
13902 continue                                                             1597
      if(jxx.gt.0)goto 13892                                               1597
      if(ixx.eq.0)goto 13892                                               1597
      if(ig.eq.0)goto 13892                                                1598
      s=-sum(b(0,:))/nc                                                    1598
      b(0,:)=b(0,:)+s                                                      1598
      sc=s                                                                 1598
      b0=0.0                                                               1599
14230 do 14231 j=1,nin                                                     1599
      l=m(j)                                                               1600
      if(vp(l) .gt. 0.0)goto 14251                                         1600
      s=sum(b(l,:))/nc                                                     1600
      goto 14261                                                           1601
14251 continue                                                             1601
      s=elc(parm,nc,b(l,:),is)                                             1601
14261 continue                                                             1602
14241 continue                                                             1602
      b(l,:)=b(l,:)-s                                                      1603
      jb=ix(l)                                                             1603
      je=ix(l+1)-1                                                         1604
      sc(jx(jb:je))=sc(jx(jb:je))-s*x(jb:je)/xs(l)                         1605
      b0=b0+s*xb(l)/xs(l)                                                  1606
14231 continue                                                             1607
14232 continue                                                             1607
      sc=sc+b0                                                             1607
      sc=exp(sc)                                                           1607
      sxp=sxp*sc                                                           1607
14270 do 14271 ic=1,nc                                                     1607
      q(:,ic)=q(:,ic)*sc                                                   1607
14271 continue                                                             1608
14272 continue                                                             1608
      nit=nit+1                                                            1608
      if(nit .le. maxit)goto 14291                                         1608
      jerr=-ilm                                                            1608
      return                                                               1608
14291 continue                                                             1609
      goto 13891                                                           1610
13892 continue                                                             1610
      if(jxx .le. 0)goto 14311                                             1610
      jerr=-10000-ilm                                                      1610
      goto 13812                                                           1610
14311 continue                                                             1610
      devi=0.0                                                             1611
14320 do 14321 ic=1,nc                                                     1612
      if(nin.gt.0) a(1:nin,ic,ilm)=b(m(1:nin),ic)                          1612
      a0(ic,ilm)=b(0,ic)                                                   1613
14330 do 14331 i=1,no                                                      1613
      if(y(i,ic).le.0.0)goto 14331                                         1614
      devi=devi-w(i)*y(i,ic)*log(q(i,ic)/sxp(i))                           1615
14331 continue                                                             1616
14332 continue                                                             1616
14321 continue                                                             1617
14322 continue                                                             1617
      kin(ilm)=nin                                                         1617
      alm(ilm)=al                                                          1617
      lmu=ilm                                                              1618
      dev(ilm)=(dev1-devi)/dev1                                            1618
      if(ig.eq.0)goto 13812                                                1619
      if(ilm.lt.mnl)goto 13811                                             1619
      if(flmin.ge.1.0)goto 13811                                           1620
      if(nintot(ni,nx,nc,a(1,1,ilm),m,nin,is).gt.ne)goto 13812             1621
      if(dev(ilm).gt.devmax)goto 13812                                     1621
      if(dev(ilm)-dev(ilm-1).lt.sml)goto 13812                             1622
13811 continue                                                             1623
13812 continue                                                             1623
      deallocate(sxp,b,bs,v,r,xv,q,mm,is,xm,sc)                            1624
      return                                                               1625
      end                                                                  1626
      subroutine lcmodval(nc,nx,a0,ca,ia,nin,x,ix,jx,n,f)                  1627
      real a0(nc),ca(nx,nc),x(*),f(nc,n)                                   1627
      integer ia(*),ix(*),jx(*)                                            1628
14340 do 14341 ic=1,nc                                                     1628
      f(ic,:)=a0(ic)                                                       1628
14341 continue                                                             1629
14342 continue                                                             1629
14350 do 14351 j=1,nin                                                     1629
      k=ia(j)                                                              1629
      kb=ix(k)                                                             1629
      ke=ix(k+1)-1                                                         1630
14360 do 14361 ic=1,nc                                                     1630
      f(ic,jx(kb:ke))=f(ic,jx(kb:ke))+ca(j,ic)*x(kb:ke)                    1630
14361 continue                                                             1631
14362 continue                                                             1631
14351 continue                                                             1632
14352 continue                                                             1632
      return                                                               1633
      end                                                                  1634
      subroutine psort7 (v,a,ii,jj)                                             
c                                                                               
c     puts into a the permutation vector which sorts v into                     
c     increasing order. the array v is not modified.                            
c     only elements from ii to jj are considered.                               
c     arrays iu(k) and il(k) permit sorting up to 2**(k+1)-1 elements           
c                                                                               
c     this is a modification of cacm algorithm #347 by r. c. singleton,         
c     which is a modified hoare quicksort.                                      
c                                                                               
      dimension a(jj),v(jj),iu(20),il(20)                                       
      integer t,tt                                                              
      integer a                                                                 
      real v                                                                    
      m=1                                                                       
      i=ii                                                                      
      j=jj                                                                      
 10   if (i.ge.j) go to 80                                                      
 20   k=i                                                                       
      ij=(j+i)/2                                                                
      t=a(ij)                                                                   
      vt=v(t)                                                                   
      if (v(a(i)).le.vt) go to 30                                               
      a(ij)=a(i)                                                                
      a(i)=t                                                                    
      t=a(ij)                                                                   
      vt=v(t)                                                                   
 30   l=j                                                                       
      if (v(a(j)).ge.vt) go to 50                                               
      a(ij)=a(j)                                                                
      a(j)=t                                                                    
      t=a(ij)                                                                   
      vt=v(t)                                                                   
      if (v(a(i)).le.vt) go to 50                                               
      a(ij)=a(i)                                                                
      a(i)=t                                                                    
      t=a(ij)                                                                   
      vt=v(t)                                                                   
      go to 50                                                                  
 40   a(l)=a(k)                                                                 
      a(k)=tt                                                                   
 50   l=l-1                                                                     
      if (v(a(l)).gt.vt) go to 50                                               
      tt=a(l)                                                                   
      vtt=v(tt)                                                                 
 60   k=k+1                                                                     
      if (v(a(k)).lt.vt) go to 60                                               
      if (k.le.l) go to 40                                                      
      if (l-i.le.j-k) go to 70                                                  
      il(m)=i                                                                   
      iu(m)=l                                                                   
      i=k                                                                       
      m=m+1                                                                     
      go to 90                                                                  
 70   il(m)=k                                                                   
      iu(m)=j                                                                   
      j=l                                                                       
      m=m+1                                                                     
      go to 90                                                                  
 80   m=m-1                                                                     
      if (m.eq.0) return                                                        
      i=il(m)                                                                   
      j=iu(m)                                                                   
 90   if (j-i.gt.10) go to 20                                                   
      if (i.eq.ii) go to 10                                                     
      i=i-1                                                                     
 100  i=i+1                                                                     
      if (i.eq.j) go to 80                                                      
      t=a(i+1)                                                                  
      vt=v(t)                                                                   
      if (v(a(i)).le.vt) go to 100                                              
      k=i                                                                       
 110  a(k+1)=a(k)                                                               
      k=k-1                                                                     
      if (vt.lt.v(a(k))) go to 110                                              
      a(k+1)=t                                                                  
      go to 100                                                                 
      end                                                                       