-
Notifications
You must be signed in to change notification settings - Fork 36
/
GLMnet.f
2892 lines (2892 loc) · 215 KB
/
GLMnet.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
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