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
|
C SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE)
C
C Algorithm NNLS: NONNEGATIVE LEAST SQUARES
C
c The original version of this code was developed by
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
c 1973 JUN 15, and published in the book
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
c
C GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN
C N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM
C
C A * X = B SUBJECT TO X .GE. 0
C ------------------------------------------------------------------
c Subroutine Arguments
c
C A(),MDA,M,N MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE
C ARRAY, A(). ON ENTRY A() CONTAINS THE M BY N
C MATRIX, A. ON EXIT A() CONTAINS
C THE PRODUCT MATRIX, Q*A , WHERE Q IS AN
C M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY
C THIS SUBROUTINE.
C B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() CON-
C TAINS Q*B.
C X() ON ENTRY X() NEED NOT BE INITIALIZED. ON EXIT X() WILL
C CONTAIN THE SOLUTION VECTOR.
C RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE
C RESIDUAL VECTOR.
C W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN
C THE DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0.
C FOR ALL I IN SET P AND W(I) .LE. 0. FOR ALL I IN SET Z
C ZZ() AN M-ARRAY OF WORKING SPACE.
C INDEX() AN INTEGER WORKING ARRAY OF LENGTH AT LEAST N.
C ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS
C P AND Z AS FOLLOWS..
C
C INDEX(1) THRU INDEX(NSETP) = SET P.
C INDEX(IZ1) THRU INDEX(IZ2) = SET Z.
C IZ1 = NSETP + 1 = NPP1
C IZ2 = N
C MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING
C MEANINGS.
C 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY.
C 2 THE DIMENSIONS OF THE PROBLEM ARE BAD.
C EITHER M .LE. 0 OR N .LE. 0.
C 3 ITERATION COUNT EXCEEDED. MORE THAN 3*N ITERATIONS.
C
C ------------------------------------------------------------------
SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE)
C ------------------------------------------------------------------
integer I, II, IP, ITER, ITMAX, IZ, IZ1, IZ2, IZMAX, J, JJ, JZ, L
integer M, MDA, MODE,N, NPP1, NSETP, RTNKEY
c integer INDEX(N)
c double precision A(MDA,N), B(M), W(N), X(N), ZZ(M)
integer INDEX(*)
double precision A(MDA,*), B(*), W(*), X(*), ZZ(*)
double precision ALPHA, ASAVE, CC, DIFF, DUMMY, FACTOR, RNORM
double precision SM, SS, T, TEMP, TWO, UNORM, UP, WMAX
double precision ZERO, ZTEST
parameter(FACTOR = 0.01d0)
parameter(TWO = 2.0d0, ZERO = 0.0d0)
C ------------------------------------------------------------------
MODE=1
IF (M .le. 0 .or. N .le. 0) then
MODE=2
RETURN
endif
ITER=0
ITMAX=3*N
C
C INITIALIZE THE ARRAYS INDEX() AND X().
C
DO 20 I=1,N
X(I)=ZERO
20 INDEX(I)=I
C
IZ2=N
IZ1=1
NSETP=0
NPP1=1
C ****** MAIN LOOP BEGINS HERE ******
30 CONTINUE
C QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION.
C OR IF M COLS OF A HAVE BEEN TRIANGULARIZED.
C
IF (IZ1 .GT.IZ2.OR.NSETP.GE.M) GO TO 350
C
C COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W().
C
DO 50 IZ=IZ1,IZ2
J=INDEX(IZ)
SM=ZERO
DO 40 L=NPP1,M
40 SM=SM+A(L,J)*B(L)
W(J)=SM
50 continue
C FIND LARGEST POSITIVE W(J).
60 continue
WMAX=ZERO
DO 70 IZ=IZ1,IZ2
J=INDEX(IZ)
IF (W(J) .gt. WMAX) then
WMAX=W(J)
IZMAX=IZ
endif
70 CONTINUE
C
C IF WMAX .LE. 0. GO TO TERMINATION.
C THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS.
C
IF (WMAX .le. ZERO) go to 350
IZ=IZMAX
J=INDEX(IZ)
C
C THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P.
C BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID
C NEAR LINEAR DEPENDENCE.
C
ASAVE=A(NPP1,J)
CALL H12 (1,NPP1,NPP1+1,M,A(1,J),1,UP,DUMMY,1,1,0)
UNORM=ZERO
IF (NSETP .ne. 0) then
DO 90 L=1,NSETP
90 UNORM=UNORM+A(L,J)**2
endif
UNORM=sqrt(UNORM)
IF (DIFF(UNORM+ABS(A(NPP1,J))*FACTOR,UNORM) .gt. ZERO) then
C
C COL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE ZZ
C AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ).
C
DO 120 L=1,M
120 ZZ(L)=B(L)
CALL H12 (2,NPP1,NPP1+1,M,A(1,J),1,UP,ZZ,1,1,1)
ZTEST=ZZ(NPP1)/A(NPP1,J)
C
C SEE IF ZTEST IS POSITIVE
C
IF (ZTEST .gt. ZERO) go to 140
endif
C
C REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P.
C RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL
C COEFFS AGAIN.
C
A(NPP1,J)=ASAVE
W(J)=ZERO
GO TO 60
C
C THE INDEX J=INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM
C SET Z TO SET P. UPDATE B, UPDATE INDICES, APPLY HOUSEHOLDER
C TRANSFORMATIONS TO COLS IN NEW SET Z, ZERO SUBDIAGONAL ELTS IN
C COL J, SET W(J)=0.
C
140 continue
DO 150 L=1,M
150 B(L)=ZZ(L)
C
INDEX(IZ)=INDEX(IZ1)
INDEX(IZ1)=J
IZ1=IZ1+1
NSETP=NPP1
NPP1=NPP1+1
C
IF (IZ1 .le. IZ2) then
DO 160 JZ=IZ1,IZ2
JJ=INDEX(JZ)
CALL H12 (2,NSETP,NPP1,M,A(1,J),1,UP,A(1,JJ),1,MDA,1)
160 continue
endif
C
IF (NSETP .ne. M) then
DO 180 L=NPP1,M
180 A(L,J)=ZERO
endif
C
W(J)=ZERO
C SOLVE THE TRIANGULAR SYSTEM.
C STORE THE SOLUTION TEMPORARILY IN ZZ().
RTNKEY = 1
GO TO 400
200 CONTINUE
C
C ****** SECONDARY LOOP BEGINS HERE ******
C
C ITERATION COUNTER.
C
210 continue
ITER=ITER+1
IF (ITER .gt. ITMAX) then
MODE=3
write (*,'(/a)') ' NNLS quitting on iteration count.'
GO TO 350
endif
C
C SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE.
C IF NOT COMPUTE ALPHA.
C
ALPHA=TWO
DO 240 IP=1,NSETP
L=INDEX(IP)
IF (ZZ(IP) .le. ZERO) then
T=-X(L)/(ZZ(IP)-X(L))
IF (ALPHA .gt. T) then
ALPHA=T
JJ=IP
endif
endif
240 CONTINUE
C
C IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL
C STILL = 2. IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP.
C
IF (ALPHA.EQ.TWO) GO TO 330
C
C OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO
C INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ.
C
DO 250 IP=1,NSETP
L=INDEX(IP)
X(L)=X(L)+ALPHA*(ZZ(IP)-X(L))
250 continue
C
C MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I
C FROM SET P TO SET Z.
C
I=INDEX(JJ)
260 continue
X(I)=ZERO
C
IF (JJ .ne. NSETP) then
JJ=JJ+1
DO 280 J=JJ,NSETP
II=INDEX(J)
INDEX(J-1)=II
CALL G1 (A(J-1,II),A(J,II),CC,SS,A(J-1,II))
A(J,II)=ZERO
DO 270 L=1,N
IF (L.NE.II) then
c
c Apply procedure G2 (CC,SS,A(J-1,L),A(J,L))
c
TEMP = A(J-1,L)
A(J-1,L) = CC*TEMP + SS*A(J,L)
A(J,L) =-SS*TEMP + CC*A(J,L)
endif
270 CONTINUE
c
c Apply procedure G2 (CC,SS,B(J-1),B(J))
c
TEMP = B(J-1)
B(J-1) = CC*TEMP + SS*B(J)
B(J) =-SS*TEMP + CC*B(J)
280 continue
endif
c
NPP1=NSETP
NSETP=NSETP-1
IZ1=IZ1-1
INDEX(IZ1)=I
C
C SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE. THEY SHOULD
C BE BECAUSE OF THE WAY ALPHA WAS DETERMINED.
C IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR. ANY
C THAT ARE NONPOSITIVE WILL BE SET TO ZERO
C AND MOVED FROM SET P TO SET Z.
C
DO 300 JJ=1,NSETP
I=INDEX(JJ)
IF (X(I) .le. ZERO) go to 260
300 CONTINUE
C
C COPY B( ) INTO ZZ( ). THEN SOLVE AGAIN AND LOOP BACK.
C
DO 310 I=1,M
310 ZZ(I)=B(I)
RTNKEY = 2
GO TO 400
320 CONTINUE
GO TO 210
C ****** END OF SECONDARY LOOP ******
C
330 continue
DO 340 IP=1,NSETP
I=INDEX(IP)
340 X(I)=ZZ(IP)
C ALL NEW COEFFS ARE POSITIVE. LOOP BACK TO BEGINNING.
GO TO 30
C
C ****** END OF MAIN LOOP ******
C
C COME TO HERE FOR TERMINATION.
C COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR.
C
350 continue
SM=ZERO
IF (NPP1 .le. M) then
DO 360 I=NPP1,M
360 SM=SM+B(I)**2
else
DO 380 J=1,N
380 W(J)=ZERO
endif
RNORM=sqrt(SM)
RETURN
C
C THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE
C TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ().
C
400 continue
DO 430 L=1,NSETP
IP=NSETP+1-L
IF (L .ne. 1) then
DO 410 II=1,IP
ZZ(II)=ZZ(II)-A(II,JJ)*ZZ(IP+1)
410 continue
endif
JJ=INDEX(IP)
ZZ(IP)=ZZ(IP)/A(IP,JJ)
430 continue
go to (200, 320), RTNKEY
END
double precision FUNCTION DIFF(X,Y)
c
c Function used in tests that depend on machine precision.
c
c The original version of this code was developed by
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
c 1973 JUN 7, and published in the book
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
C
double precision X, Y
DIFF=X-Y
RETURN
END
C SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)
C
C CONSTRUCTION AND/OR APPLICATION OF A SINGLE
C HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B
C
c The original version of this code was developed by
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
c 1973 JUN 12, and published in the book
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
C ------------------------------------------------------------------
c Subroutine Arguments
c
C MODE = 1 OR 2 Selects Algorithm H1 to construct and apply a
c Householder transformation, or Algorithm H2 to apply a
c previously constructed transformation.
C LPIVOT IS THE INDEX OF THE PIVOT ELEMENT.
C L1,M IF L1 .LE. M THE TRANSFORMATION WILL BE CONSTRUCTED TO
C ZERO ELEMENTS INDEXED FROM L1 THROUGH M. IF L1 GT. M
C THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION.
C U(),IUE,UP On entry with MODE = 1, U() contains the pivot
c vector. IUE is the storage increment between elements.
c On exit when MODE = 1, U() and UP contain quantities
c defining the vector U of the Householder transformation.
c on entry with MODE = 2, U() and UP should contain
c quantities previously computed with MODE = 1. These will
c not be modified during the entry with MODE = 2.
C C() ON ENTRY with MODE = 1 or 2, C() CONTAINS A MATRIX WHICH
c WILL BE REGARDED AS A SET OF VECTORS TO WHICH THE
c HOUSEHOLDER TRANSFORMATION IS TO BE APPLIED.
c ON EXIT C() CONTAINS THE SET OF TRANSFORMED VECTORS.
C ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C().
C ICV STORAGE INCREMENT BETWEEN VECTORS IN C().
C NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0
C NO OPERATIONS WILL BE DONE ON C().
C ------------------------------------------------------------------
SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)
C ------------------------------------------------------------------
integer I, I2, I3, I4, ICE, ICV, INCR, IUE, J
integer L1, LPIVOT, M, MODE, NCV
double precision B, C(*), CL, CLINV, ONE, SM
c double precision U(IUE,M)
double precision U(IUE,*)
double precision UP
parameter(ONE = 1.0d0)
C ------------------------------------------------------------------
IF (0.GE.LPIVOT.OR.LPIVOT.GE.L1.OR.L1.GT.M) RETURN
CL=abs(U(1,LPIVOT))
IF (MODE.EQ.2) GO TO 60
C ****** CONSTRUCT THE TRANSFORMATION. ******
DO 10 J=L1,M
10 CL=MAX(abs(U(1,J)),CL)
IF (CL) 130,130,20
20 CLINV=ONE/CL
SM=(U(1,LPIVOT)*CLINV)**2
DO 30 J=L1,M
30 SM=SM+(U(1,J)*CLINV)**2
CL=CL*SQRT(SM)
IF (U(1,LPIVOT)) 50,50,40
40 CL=-CL
50 UP=U(1,LPIVOT)-CL
U(1,LPIVOT)=CL
GO TO 70
C ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C. ******
C
60 IF (CL) 130,130,70
70 IF (NCV.LE.0) RETURN
B= UP*U(1,LPIVOT)
C B MUST BE NONPOSITIVE HERE. IF B = 0., RETURN.
C
IF (B) 80,130,130
80 B=ONE/B
I2=1-ICV+ICE*(LPIVOT-1)
INCR=ICE*(L1-LPIVOT)
DO 120 J=1,NCV
I2=I2+ICV
I3=I2+INCR
I4=I3
SM=C(I2)*UP
DO 90 I=L1,M
SM=SM+C(I3)*U(1,I)
90 I3=I3+ICE
IF (SM) 100,120,100
100 SM=SM*B
C(I2)=C(I2)+SM*UP
DO 110 I=L1,M
C(I4)=C(I4)+SM*U(1,I)
110 I4=I4+ICE
120 CONTINUE
130 RETURN
END
SUBROUTINE G1 (A,B,CTERM,STERM,SIG)
c
C COMPUTE ORTHOGONAL ROTATION MATRIX..
c
c The original version of this code was developed by
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
c 1973 JUN 12, and published in the book
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
C
C COMPUTE.. MATRIX (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2))
C (-S,C) (-S,C)(B) ( 0 )
C COMPUTE SIG = SQRT(A**2+B**2)
C SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT
C SIG MAY BE IN THE SAME LOCATION AS A OR B .
C ------------------------------------------------------------------
double precision A, B, CTERM, ONE, SIG, STERM, XR, YR, ZERO
parameter(ONE = 1.0d0, ZERO = 0.0d0)
C ------------------------------------------------------------------
if (abs(A) .gt. abs(B)) then
XR=B/A
YR=sqrt(ONE+XR**2)
CTERM=sign(ONE/YR,A)
STERM=CTERM*XR
SIG=abs(A)*YR
RETURN
endif
if (B .ne. ZERO) then
XR=A/B
YR=sqrt(ONE+XR**2)
STERM=sign(ONE/YR,B)
CTERM=STERM*XR
SIG=abs(B)*YR
RETURN
endif
SIG=ZERO
CTERM=ZERO
STERM=ONE
RETURN
END
|
