Mercurial > hg > nnls-chroma

comparison nnls.f @ 22:444c344681f3 matthiasm-plugin

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
* Rather than worry about provenance of C versions, why not just use the FORTRAN? (We can back this out if it doesn't go well with the build scripts)
author Chris Cannam
date Thu, 21 Oct 2010 11:58:28 +0100
parents
children 6d9e1ee7b35a
comparison
equal deleted inserted replaced
21:1e2d0db15920 22:444c344681f3
1
2 double precision FUNCTION DIFF(X,Y)
3 c
4 c Function used in tests that depend on machine precision.
5 c
6 c The original version of this code was developed by
7 c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
8 c 1973 JUN 7, and published in the book
9 c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
10 c Revised FEB 1995 to accompany reprinting of the book by SIAM.
11 C
12 double precision X, Y
13 DIFF=X-Y
14 RETURN
15 END
16
17
18 SUBROUTINE G1 (A,B,CTERM,STERM,SIG)
19 c
20 C COMPUTE ORTHOGONAL ROTATION MATRIX..
21 c
22 c The original version of this code was developed by
23 c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
24 c 1973 JUN 12, and published in the book
25 c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
26 c Revised FEB 1995 to accompany reprinting of the book by SIAM.
27 C
28 C COMPUTE.. MATRIX (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2))
29 C (-S,C) (-S,C)(B) ( 0 )
30 C COMPUTE SIG = SQRT(A**2+B**2)
31 C SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT
32 C SIG MAY BE IN THE SAME LOCATION AS A OR B .
33 C ------------------------------------------------------------------
34 double precision A, B, CTERM, ONE, SIG, STERM, XR, YR, ZERO
35 parameter(ONE = 1.0d0, ZERO = 0.0d0)
36 C ------------------------------------------------------------------
37 if (abs(A) .gt. abs(B)) then
38 XR=B/A
39 YR=sqrt(ONE+XR**2)
40 CTERM=sign(ONE/YR,A)
41 STERM=CTERM*XR
42 SIG=abs(A)*YR
43 RETURN
44 endif
45
46 if (B .ne. ZERO) then
47 XR=A/B
48 YR=sqrt(ONE+XR**2)
49 STERM=sign(ONE/YR,B)
50 CTERM=STERM*XR
51 SIG=abs(B)*YR
52 RETURN
53 endif
54
55 SIG=ZERO
56 CTERM=ZERO
57 STERM=ONE
58 RETURN
59 END
60
61
62 C SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)
63 C
64 C CONSTRUCTION AND/OR APPLICATION OF A SINGLE
65 C HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B
66 C
67 c The original version of this code was developed by
68 c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
69 c 1973 JUN 12, and published in the book
70 c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
71 c Revised FEB 1995 to accompany reprinting of the book by SIAM.
72 C ------------------------------------------------------------------
73 c Subroutine Arguments
74 c
75 C MODE = 1 OR 2 Selects Algorithm H1 to construct and apply a
76 c Householder transformation, or Algorithm H2 to apply a
77 c previously constructed transformation.
78 C LPIVOT IS THE INDEX OF THE PIVOT ELEMENT.
79 C L1,M IF L1 .LE. M THE TRANSFORMATION WILL BE CONSTRUCTED TO
80 C ZERO ELEMENTS INDEXED FROM L1 THROUGH M. IF L1 GT. M
81 C THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION.
82 C U(),IUE,UP On entry with MODE = 1, U() contains the pivot
83 c vector. IUE is the storage increment between elements.
84 c On exit when MODE = 1, U() and UP contain quantities
85 c defining the vector U of the Householder transformation.
86 c on entry with MODE = 2, U() and UP should contain
87 c quantities previously computed with MODE = 1. These will
88 c not be modified during the entry with MODE = 2.
89 C C() ON ENTRY with MODE = 1 or 2, C() CONTAINS A MATRIX WHICH
90 c WILL BE REGARDED AS A SET OF VECTORS TO WHICH THE
91 c HOUSEHOLDER TRANSFORMATION IS TO BE APPLIED.
92 c ON EXIT C() CONTAINS THE SET OF TRANSFORMED VECTORS.
93 C ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C().
94 C ICV STORAGE INCREMENT BETWEEN VECTORS IN C().
95 C NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0
96 C NO OPERATIONS WILL BE DONE ON C().
97 C ------------------------------------------------------------------
98 SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)
99 C ------------------------------------------------------------------
100 integer I, I2, I3, I4, ICE, ICV, INCR, IUE, J
101 integer L1, LPIVOT, M, MODE, NCV
102 double precision B, C(*), CL, CLINV, ONE, SM
103 c double precision U(IUE,M)
104 double precision U(IUE,*)
105 double precision UP
106 parameter(ONE = 1.0d0)
107 C ------------------------------------------------------------------
108 IF (0.GE.LPIVOT.OR.LPIVOT.GE.L1.OR.L1.GT.M) RETURN
109 CL=abs(U(1,LPIVOT))
110 IF (MODE.EQ.2) GO TO 60
111 C ****** CONSTRUCT THE TRANSFORMATION. ******
112 DO 10 J=L1,M
113 10 CL=MAX(abs(U(1,J)),CL)
114 IF (CL) 130,130,20
115 20 CLINV=ONE/CL
116 SM=(U(1,LPIVOT)*CLINV)**2
117 DO 30 J=L1,M
118 30 SM=SM+(U(1,J)*CLINV)**2
119 CL=CL*SQRT(SM)
120 IF (U(1,LPIVOT)) 50,50,40
121 40 CL=-CL
122 50 UP=U(1,LPIVOT)-CL
123 U(1,LPIVOT)=CL
124 GO TO 70
125 C ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C. ******
126 C
127 60 IF (CL) 130,130,70
128 70 IF (NCV.LE.0) RETURN
129 B= UP*U(1,LPIVOT)
130 C B MUST BE NONPOSITIVE HERE. IF B = 0., RETURN.
131 C
132 IF (B) 80,130,130
133 80 B=ONE/B
134 I2=1-ICV+ICE*(LPIVOT-1)
135 INCR=ICE*(L1-LPIVOT)
136 DO 120 J=1,NCV
137 I2=I2+ICV
138 I3=I2+INCR
139 I4=I3
140 SM=C(I2)*UP
141 DO 90 I=L1,M
142 SM=SM+C(I3)*U(1,I)
143 90 I3=I3+ICE
144 IF (SM) 100,120,100
145 100 SM=SM*B
146 C(I2)=C(I2)+SM*UP
147 DO 110 I=L1,M
148 C(I4)=C(I4)+SM*U(1,I)
149 110 I4=I4+ICE
150 120 CONTINUE
151 130 RETURN
152 END
153
154
155
156 C SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE)
157 C
158 C Algorithm NNLS: NONNEGATIVE LEAST SQUARES
159 C
160 c The original version of this code was developed by
161 c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
162 c 1973 JUN 15, and published in the book
163 c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
164 c Revised FEB 1995 to accompany reprinting of the book by SIAM.
165 c
166 C GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN
167 C N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM
168 C
169 C A * X = B SUBJECT TO X .GE. 0
170 C ------------------------------------------------------------------
171 c Subroutine Arguments
172 c
173 C A(),MDA,M,N MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE
174 C ARRAY, A(). ON ENTRY A() CONTAINS THE M BY N
175 C MATRIX, A. ON EXIT A() CONTAINS
176 C THE PRODUCT MATRIX, Q*A , WHERE Q IS AN
177 C M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY
178 C THIS SUBROUTINE.
179 C B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() CON-
180 C TAINS Q*B.
181 C X() ON ENTRY X() NEED NOT BE INITIALIZED. ON EXIT X() WILL
182 C CONTAIN THE SOLUTION VECTOR.
183 C RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE
184 C RESIDUAL VECTOR.
185 C W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN
186 C THE DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0.
187 C FOR ALL I IN SET P AND W(I) .LE. 0. FOR ALL I IN SET Z
188 C ZZ() AN M-ARRAY OF WORKING SPACE.
189 C INDEX() AN INTEGER WORKING ARRAY OF LENGTH AT LEAST N.
190 C ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS
191 C P AND Z AS FOLLOWS..
192 C
193 C INDEX(1) THRU INDEX(NSETP) = SET P.
194 C INDEX(IZ1) THRU INDEX(IZ2) = SET Z.
195 C IZ1 = NSETP + 1 = NPP1
196 C IZ2 = N
197 C MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING
198 C MEANINGS.
199 C 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY.
200 C 2 THE DIMENSIONS OF THE PROBLEM ARE BAD.
201 C EITHER M .LE. 0 OR N .LE. 0.
202 C 3 ITERATION COUNT EXCEEDED. MORE THAN 3*N ITERATIONS.
203 C
204 C ------------------------------------------------------------------
205 SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE)
206 C ------------------------------------------------------------------
207 integer I, II, IP, ITER, ITMAX, IZ, IZ1, IZ2, IZMAX, J, JJ, JZ, L
208 integer M, MDA, MODE,N, NPP1, NSETP, RTNKEY
209 c integer INDEX(N)
210 c double precision A(MDA,N), B(M), W(N), X(N), ZZ(M)
211 integer INDEX(*)
212 double precision A(MDA,*), B(*), W(*), X(*), ZZ(*)
213 double precision ALPHA, ASAVE, CC, DIFF, DUMMY, FACTOR, RNORM
214 double precision SM, SS, T, TEMP, TWO, UNORM, UP, WMAX
215 double precision ZERO, ZTEST
216 parameter(FACTOR = 0.01d0)
217 parameter(TWO = 2.0d0, ZERO = 0.0d0)
218 C ------------------------------------------------------------------
219 MODE=1
220 IF (M .le. 0 .or. N .le. 0) then
221 MODE=2
222 RETURN
223 endif
224 ITER=0
225 ITMAX=3*N
226 C
227 C INITIALIZE THE ARRAYS INDEX() AND X().
228 C
229 DO 20 I=1,N
230 X(I)=ZERO
231 20 INDEX(I)=I
232 C
233 IZ2=N
234 IZ1=1
235 NSETP=0
236 NPP1=1
237 C ****** MAIN LOOP BEGINS HERE ******
238 30 CONTINUE
239 C QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION.
240 C OR IF M COLS OF A HAVE BEEN TRIANGULARIZED.
241 C
242 IF (IZ1 .GT.IZ2.OR.NSETP.GE.M) GO TO 350
243 C
244 C COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W().
245 C
246 DO 50 IZ=IZ1,IZ2
247 J=INDEX(IZ)
248 SM=ZERO
249 DO 40 L=NPP1,M
250 40 SM=SM+A(L,J)*B(L)
251 W(J)=SM
252 50 continue
253 C FIND LARGEST POSITIVE W(J).
254 60 continue
255 WMAX=ZERO
256 DO 70 IZ=IZ1,IZ2
257 J=INDEX(IZ)
258 IF (W(J) .gt. WMAX) then
259 WMAX=W(J)
260 IZMAX=IZ
261 endif
262 70 CONTINUE
263 C
264 C IF WMAX .LE. 0. GO TO TERMINATION.
265 C THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS.
266 C
267 IF (WMAX .le. ZERO) go to 350
268 IZ=IZMAX
269 J=INDEX(IZ)
270 C
271 C THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P.
272 C BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID
273 C NEAR LINEAR DEPENDENCE.
274 C
275 ASAVE=A(NPP1,J)
276 CALL H12 (1,NPP1,NPP1+1,M,A(1,J),1,UP,DUMMY,1,1,0)
277 UNORM=ZERO
278 IF (NSETP .ne. 0) then
279 DO 90 L=1,NSETP
280 90 UNORM=UNORM+A(L,J)**2
281 endif
282 UNORM=sqrt(UNORM)
283 IF (DIFF(UNORM+ABS(A(NPP1,J))*FACTOR,UNORM) .gt. ZERO) then
284 C
285 C COL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE ZZ
286 C AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ).
287 C
288 DO 120 L=1,M
289 120 ZZ(L)=B(L)
290 CALL H12 (2,NPP1,NPP1+1,M,A(1,J),1,UP,ZZ,1,1,1)
291 ZTEST=ZZ(NPP1)/A(NPP1,J)
292 C
293 C SEE IF ZTEST IS POSITIVE
294 C
295 IF (ZTEST .gt. ZERO) go to 140
296 endif
297 C
298 C REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P.
299 C RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL
300 C COEFFS AGAIN.
301 C
302 A(NPP1,J)=ASAVE
303 W(J)=ZERO
304 GO TO 60
305 C
306 C THE INDEX J=INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM
307 C SET Z TO SET P. UPDATE B, UPDATE INDICES, APPLY HOUSEHOLDER
308 C TRANSFORMATIONS TO COLS IN NEW SET Z, ZERO SUBDIAGONAL ELTS IN
309 C COL J, SET W(J)=0.
310 C
311 140 continue
312 DO 150 L=1,M
313 150 B(L)=ZZ(L)
314 C
315 INDEX(IZ)=INDEX(IZ1)
316 INDEX(IZ1)=J
317 IZ1=IZ1+1
318 NSETP=NPP1
319 NPP1=NPP1+1
320 C
321 IF (IZ1 .le. IZ2) then
322 DO 160 JZ=IZ1,IZ2
323 JJ=INDEX(JZ)
324 CALL H12 (2,NSETP,NPP1,M,A(1,J),1,UP,A(1,JJ),1,MDA,1)
325 160 continue
326 endif
327 C
328 IF (NSETP .ne. M) then
329 DO 180 L=NPP1,M
330 180 A(L,J)=ZERO
331 endif
332 C
333 W(J)=ZERO
334 C SOLVE THE TRIANGULAR SYSTEM.
335 C STORE THE SOLUTION TEMPORARILY IN ZZ().
336 RTNKEY = 1
337 GO TO 400
338 200 CONTINUE
339 C
340 C ****** SECONDARY LOOP BEGINS HERE ******
341 C
342 C ITERATION COUNTER.
343 C
344 210 continue
345 ITER=ITER+1
346 IF (ITER .gt. ITMAX) then
347 MODE=3
348 write (*,'(/a)') ' NNLS quitting on iteration count.'
349 GO TO 350
350 endif
351 C
352 C SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE.
353 C IF NOT COMPUTE ALPHA.
354 C
355 ALPHA=TWO
356 DO 240 IP=1,NSETP
357 L=INDEX(IP)
358 IF (ZZ(IP) .le. ZERO) then
359 T=-X(L)/(ZZ(IP)-X(L))
360 IF (ALPHA .gt. T) then
361 ALPHA=T
362 JJ=IP
363 endif
364 endif
365 240 CONTINUE
366 C
367 C IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL
368 C STILL = 2. IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP.
369 C
370 IF (ALPHA.EQ.TWO) GO TO 330
371 C
372 C OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO
373 C INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ.
374 C
375 DO 250 IP=1,NSETP
376 L=INDEX(IP)
377 X(L)=X(L)+ALPHA*(ZZ(IP)-X(L))
378 250 continue
379 C
380 C MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I
381 C FROM SET P TO SET Z.
382 C
383 I=INDEX(JJ)
384 260 continue
385 X(I)=ZERO
386 C
387 IF (JJ .ne. NSETP) then
388 JJ=JJ+1
389 DO 280 J=JJ,NSETP
390 II=INDEX(J)
391 INDEX(J-1)=II
392 CALL G1 (A(J-1,II),A(J,II),CC,SS,A(J-1,II))
393 A(J,II)=ZERO
394 DO 270 L=1,N
395 IF (L.NE.II) then
396 c
397 c Apply procedure G2 (CC,SS,A(J-1,L),A(J,L))
398 c
399 TEMP = A(J-1,L)
400 A(J-1,L) = CC*TEMP + SS*A(J,L)
401 A(J,L) =-SS*TEMP + CC*A(J,L)
402 endif
403 270 CONTINUE
404 c
405 c Apply procedure G2 (CC,SS,B(J-1),B(J))
406 c
407 TEMP = B(J-1)
408 B(J-1) = CC*TEMP + SS*B(J)
409 B(J) =-SS*TEMP + CC*B(J)
410 280 continue
411 endif
412 c
413 NPP1=NSETP
414 NSETP=NSETP-1
415 IZ1=IZ1-1
416 INDEX(IZ1)=I
417 C
418 C SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE. THEY SHOULD
419 C BE BECAUSE OF THE WAY ALPHA WAS DETERMINED.
420 C IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR. ANY
421 C THAT ARE NONPOSITIVE WILL BE SET TO ZERO
422 C AND MOVED FROM SET P TO SET Z.
423 C
424 DO 300 JJ=1,NSETP
425 I=INDEX(JJ)
426 IF (X(I) .le. ZERO) go to 260
427 300 CONTINUE
428 C
429 C COPY B( ) INTO ZZ( ). THEN SOLVE AGAIN AND LOOP BACK.
430 C
431 DO 310 I=1,M
432 310 ZZ(I)=B(I)
433 RTNKEY = 2
434 GO TO 400
435 320 CONTINUE
436 GO TO 210
437 C ****** END OF SECONDARY LOOP ******
438 C
439 330 continue
440 DO 340 IP=1,NSETP
441 I=INDEX(IP)
442 340 X(I)=ZZ(IP)
443 C ALL NEW COEFFS ARE POSITIVE. LOOP BACK TO BEGINNING.
444 GO TO 30
445 C
446 C ****** END OF MAIN LOOP ******
447 C
448 C COME TO HERE FOR TERMINATION.
449 C COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR.
450 C
451 350 continue
452 SM=ZERO
453 IF (NPP1 .le. M) then
454 DO 360 I=NPP1,M
455 360 SM=SM+B(I)**2
456 else
457 DO 380 J=1,N
458 380 W(J)=ZERO
459 endif
460 RNORM=sqrt(SM)
461 RETURN
462 C
463 C THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE
464 C TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ().
465 C
466 400 continue
467 DO 430 L=1,NSETP
468 IP=NSETP+1-L
469 IF (L .ne. 1) then
470 DO 410 II=1,IP
471 ZZ(II)=ZZ(II)-A(II,JJ)*ZZ(IP+1)
472 410 continue
473 endif
474 JJ=INDEX(IP)
475 ZZ(IP)=ZZ(IP)/A(IP,JJ)
476 430 continue
477 go to (200, 320), RTNKEY
478 END
479