Mercurial > hg > nnls-chroma

comparison nnls.c @ 25:6d9e1ee7b35a matthiasm-plugin

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
* OK, let's revert to a C version (after discovering that the FORTRAN compiler isn't standard on the Mac)
author Chris Cannam
date Thu, 21 Oct 2010 14:38:02 +0100
parents
children da3195577172
comparison
equal deleted inserted replaced
24:568ff0daa659 25:6d9e1ee7b35a
1
2 #include "nnls.h"
3
4 typedef int integer;
5 typedef unsigned int uinteger;
6 typedef char *address;
7 typedef short int shortint;
8 typedef float real;
9 typedef double doublereal;
10
11 #define abs(x) ((x) >= 0 ? (x) : -(x))
12 #define dabs(x) (doublereal)abs(x)
13 #define min(a,b) ((a) <= (b) ? (a) : (b))
14 #define max(a,b) ((a) >= (b) ? (a) : (b))
15 #define dmin(a,b) (doublereal)min(a,b)
16 #define dmax(a,b) (doublereal)max(a,b)
17 #define bit_test(a,b) ((a) >> (b) & 1)
18 #define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
19 #define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
20
21 double d_sign(doublereal *a, doublereal *b)
22 {
23 double x;
24 x = (*a >= 0 ? *a : - *a);
25 return (*b >= 0 ? x : -x);
26 }
27
28 /* Table of constant values */
29
30 static integer c__1 = 1;
31 static integer c__0 = 0;
32 static integer c__2 = 2;
33
34 doublereal diff_(doublereal *x, doublereal *y)
35 {
36 /* System generated locals */
37 doublereal ret_val;
38
39
40 /* Function used in tests that depend on machine precision. */
41
42 /* The original version of this code was developed by */
43 /* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */
44 /* 1973 JUN 7, and published in the book */
45 /* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
46 /* Revised FEB 1995 to accompany reprinting of the book by SIAM. */
47
48 ret_val = *x - *y;
49 return ret_val;
50 } /* diff_ */
51
52 /* Subroutine */ int g1_(doublereal *a, doublereal *b, doublereal *cterm,
53 doublereal *sterm, doublereal *sig)
54 {
55 /* System generated locals */
56 doublereal d__1;
57
58 /* Builtin functions */
59 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
60
61 /* Local variables */
62 static doublereal xr, yr;
63
64
65 /* COMPUTE ORTHOGONAL ROTATION MATRIX.. */
66
67 /* The original version of this code was developed by */
68 /* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */
69 /* 1973 JUN 12, and published in the book */
70 /* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
71 /* Revised FEB 1995 to accompany reprinting of the book by SIAM. */
72
73 /* COMPUTE.. MATRIX (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2)) */
74 /* (-S,C) (-S,C)(B) ( 0 ) */
75 /* COMPUTE SIG = SQRT(A**2+B**2) */
76 /* SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT */
77 /* SIG MAY BE IN THE SAME LOCATION AS A OR B . */
78 /* ------------------------------------------------------------------ */
79 /* ------------------------------------------------------------------ */
80 if (abs(*a) > abs(*b)) {
81 xr = *b / *a;
82 /* Computing 2nd power */
83 d__1 = xr;
84 yr = sqrt(d__1 * d__1 + 1.);
85 d__1 = 1. / yr;
86 *cterm = d_sign(&d__1, a);
87 *sterm = *cterm * xr;
88 *sig = abs(*a) * yr;
89 return 0;
90 }
91 if (*b != 0.) {
92 xr = *a / *b;
93 /* Computing 2nd power */
94 d__1 = xr;
95 yr = sqrt(d__1 * d__1 + 1.);
96 d__1 = 1. / yr;
97 *sterm = d_sign(&d__1, b);
98 *cterm = *sterm * xr;
99 *sig = abs(*b) * yr;
100 return 0;
101 }
102 *sig = 0.;
103 *cterm = 0.;
104 *sterm = 1.;
105 return 0;
106 } /* g1_ */
107
108 /* SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) */
109
110 /* CONSTRUCTION AND/OR APPLICATION OF A SINGLE */
111 /* HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B */
112
113 /* The original version of this code was developed by */
114 /* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */
115 /* 1973 JUN 12, and published in the book */
116 /* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
117 /* Revised FEB 1995 to accompany reprinting of the book by SIAM. */
118 /* ------------------------------------------------------------------ */
119 /* Subroutine Arguments */
120
121 /* MODE = 1 OR 2 Selects Algorithm H1 to construct and apply a */
122 /* Householder transformation, or Algorithm H2 to apply a */
123 /* previously constructed transformation. */
124 /* LPIVOT IS THE INDEX OF THE PIVOT ELEMENT. */
125 /* L1,M IF L1 .LE. M THE TRANSFORMATION WILL BE CONSTRUCTED TO */
126 /* ZERO ELEMENTS INDEXED FROM L1 THROUGH M. IF L1 GT. M */
127 /* THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION. */
128 /* U(),IUE,UP On entry with MODE = 1, U() contains the pivot */
129 /* vector. IUE is the storage increment between elements. */
130 /* On exit when MODE = 1, U() and UP contain quantities */
131 /* defining the vector U of the Householder transformation. */
132 /* on entry with MODE = 2, U() and UP should contain */
133 /* quantities previously computed with MODE = 1. These will */
134 /* not be modified during the entry with MODE = 2. */
135 /* C() ON ENTRY with MODE = 1 or 2, C() CONTAINS A MATRIX WHICH */
136 /* WILL BE REGARDED AS A SET OF VECTORS TO WHICH THE */
137 /* HOUSEHOLDER TRANSFORMATION IS TO BE APPLIED. */
138 /* ON EXIT C() CONTAINS THE SET OF TRANSFORMED VECTORS. */
139 /* ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C(). */
140 /* ICV STORAGE INCREMENT BETWEEN VECTORS IN C(). */
141 /* NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0 */
142 /* NO OPERATIONS WILL BE DONE ON C(). */
143 /* ------------------------------------------------------------------ */
144 /* Subroutine */ int h12_(integer *mode, integer *lpivot, integer *l1,
145 integer *m, doublereal *u, integer *iue, doublereal *up, doublereal *
146 c__, integer *ice, integer *icv, integer *ncv)
147 {
148 /* System generated locals */
149 integer u_dim1, u_offset, i__1, i__2;
150 doublereal d__1, d__2;
151
152 /* Builtin functions */
153 double sqrt(doublereal);
154
155 /* Local variables */
156 static doublereal b;
157 static integer i__, j, i2, i3, i4;
158 static doublereal cl, sm;
159 static integer incr;
160 static doublereal clinv;
161
162 /* ------------------------------------------------------------------ */
163 /* double precision U(IUE,M) */
164 /* ------------------------------------------------------------------ */
165 /* Parameter adjustments */
166 u_dim1 = *iue;
167 u_offset = 1 + u_dim1;
168 u -= u_offset;
169 --c__;
170
171 /* Function Body */
172 if (0 >= *lpivot || *lpivot >= *l1 || *l1 > *m) {
173 return 0;
174 }
175 cl = (d__1 = u[*lpivot * u_dim1 + 1], abs(d__1));
176 if (*mode == 2) {
177 goto L60;
178 }
179 /* ****** CONSTRUCT THE TRANSFORMATION. ****** */
180 i__1 = *m;
181 for (j = *l1; j <= i__1; ++j) {
182 /* L10: */
183 /* Computing MAX */
184 d__2 = (d__1 = u[j * u_dim1 + 1], abs(d__1));
185 cl = max(d__2,cl);
186 }
187 if (cl <= 0.) {
188 goto L130;
189 } else {
190 goto L20;
191 }
192 L20:
193 clinv = 1. / cl;
194 /* Computing 2nd power */
195 d__1 = u[*lpivot * u_dim1 + 1] * clinv;
196 sm = d__1 * d__1;
197 i__1 = *m;
198 for (j = *l1; j <= i__1; ++j) {
199 /* L30: */
200 /* Computing 2nd power */
201 d__1 = u[j * u_dim1 + 1] * clinv;
202 sm += d__1 * d__1;
203 }
204 cl *= sqrt(sm);
205 if (u[*lpivot * u_dim1 + 1] <= 0.) {
206 goto L50;
207 } else {
208 goto L40;
209 }
210 L40:
211 cl = -cl;
212 L50:
213 *up = u[*lpivot * u_dim1 + 1] - cl;
214 u[*lpivot * u_dim1 + 1] = cl;
215 goto L70;
216 /* ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C. ****** */
217
218 L60:
219 if (cl <= 0.) {
220 goto L130;
221 } else {
222 goto L70;
223 }
224 L70:
225 if (*ncv <= 0) {
226 return 0;
227 }
228 b = *up * u[*lpivot * u_dim1 + 1];
229 /* B MUST BE NONPOSITIVE HERE. IF B = 0., RETURN. */
230
231 if (b >= 0.) {
232 goto L130;
233 } else {
234 goto L80;
235 }
236 L80:
237 b = 1. / b;
238 i2 = 1 - *icv + *ice * (*lpivot - 1);
239 incr = *ice * (*l1 - *lpivot);
240 i__1 = *ncv;
241 for (j = 1; j <= i__1; ++j) {
242 i2 += *icv;
243 i3 = i2 + incr;
244 i4 = i3;
245 sm = c__[i2] * *up;
246 i__2 = *m;
247 for (i__ = *l1; i__ <= i__2; ++i__) {
248 sm += c__[i3] * u[i__ * u_dim1 + 1];
249 /* L90: */
250 i3 += *ice;
251 }
252 if (sm != 0.) {
253 goto L100;
254 } else {
255 goto L120;
256 }
257 L100:
258 sm *= b;
259 c__[i2] += sm * *up;
260 i__2 = *m;
261 for (i__ = *l1; i__ <= i__2; ++i__) {
262 c__[i4] += sm * u[i__ * u_dim1 + 1];
263 /* L110: */
264 i4 += *ice;
265 }
266 L120:
267 ;
268 }
269 L130:
270 return 0;
271 } /* h12_ */
272
273 /* SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) */
274
275 /* Algorithm NNLS: NONNEGATIVE LEAST SQUARES */
276
277 /* The original version of this code was developed by */
278 /* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */
279 /* 1973 JUN 15, and published in the book */
280 /* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
281 /* Revised FEB 1995 to accompany reprinting of the book by SIAM. */
282
283 /* GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN */
284 /* N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM */
285
286 /* A * X = B SUBJECT TO X .GE. 0 */
287 /* ------------------------------------------------------------------ */
288 /* Subroutine Arguments */
289
290 /* A(),MDA,M,N MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE */
291 /* ARRAY, A(). ON ENTRY A() CONTAINS THE M BY N */
292 /* MATRIX, A. ON EXIT A() CONTAINS */
293 /* THE PRODUCT MATRIX, Q*A , WHERE Q IS AN */
294 /* M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY */
295 /* THIS SUBROUTINE. */
296 /* B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() CON- */
297 /* TAINS Q*B. */
298 /* X() ON ENTRY X() NEED NOT BE INITIALIZED. ON EXIT X() WILL */
299 /* CONTAIN THE SOLUTION VECTOR. */
300 /* RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE */
301 /* RESIDUAL VECTOR. */
302 /* W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN */
303 /* THE DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0. */
304 /* FOR ALL I IN SET P AND W(I) .LE. 0. FOR ALL I IN SET Z */
305 /* ZZ() AN M-ARRAY OF WORKING SPACE. */
306 /* INDEX() AN INTEGER WORKING ARRAY OF LENGTH AT LEAST N. */
307 /* ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS */
308 /* P AND Z AS FOLLOWS.. */
309
310 /* INDEX(1) THRU INDEX(NSETP) = SET P. */
311 /* INDEX(IZ1) THRU INDEX(IZ2) = SET Z. */
312 /* IZ1 = NSETP + 1 = NPP1 */
313 /* IZ2 = N */
314 /* MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING */
315 /* MEANINGS. */
316 /* 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY. */
317 /* 2 THE DIMENSIONS OF THE PROBLEM ARE BAD. */
318 /* EITHER M .LE. 0 OR N .LE. 0. */
319 /* 3 ITERATION COUNT EXCEEDED. MORE THAN 3*N ITERATIONS. */
320
321 /* ------------------------------------------------------------------ */
322 /* Subroutine */ int nnls_(doublereal *a, integer *mda, integer *m, integer *
323 n, doublereal *b, doublereal *x, doublereal *rnorm, doublereal *w,
324 doublereal *zz, integer *index, integer *mode)
325 {
326 /* System generated locals */
327 integer a_dim1, a_offset, i__1, i__2;
328 doublereal d__1, d__2;
329
330 /* Builtin functions */
331 double sqrt(doublereal);
332
333 /* Local variables */
334 static integer i__, j, l;
335 static doublereal t;
336 extern /* Subroutine */ int g1_(doublereal *, doublereal *, doublereal *,
337 doublereal *, doublereal *);
338 static doublereal cc;
339 extern /* Subroutine */ int h12_(integer *, integer *, integer *, integer
340 *, doublereal *, integer *, doublereal *, doublereal *, integer *,
341 integer *, integer *);
342 static integer ii, jj, ip;
343 static doublereal sm;
344 static integer iz, jz;
345 static doublereal up, ss;
346 static integer iz1, iz2, npp1;
347 extern doublereal diff_(doublereal *, doublereal *);
348 static integer iter;
349 static doublereal temp, wmax, alpha, asave;
350 static integer itmax, izmax, nsetp;
351 static doublereal dummy, unorm, ztest;
352 static integer rtnkey;
353
354 /* ------------------------------------------------------------------ */
355 /* integer INDEX(N) */
356 /* double precision A(MDA,N), B(M), W(N), X(N), ZZ(M) */
357 /* ------------------------------------------------------------------ */
358 /* Parameter adjustments */
359 a_dim1 = *mda;
360 a_offset = 1 + a_dim1;
361 a -= a_offset;
362 --b;
363 --x;
364 --w;
365 --zz;
366 --index;
367
368 /* Function Body */
369 *mode = 1;
370 if (*m <= 0 || *n <= 0) {
371 *mode = 2;
372 return 0;
373 }
374 iter = 0;
375 itmax = *n * 3;
376
377 /* INITIALIZE THE ARRAYS INDEX() AND X(). */
378
379 i__1 = *n;
380 for (i__ = 1; i__ <= i__1; ++i__) {
381 x[i__] = 0.;
382 /* L20: */
383 index[i__] = i__;
384 }
385
386 iz2 = *n;
387 iz1 = 1;
388 nsetp = 0;
389 npp1 = 1;
390 /* ****** MAIN LOOP BEGINS HERE ****** */
391 L30:
392 /* QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION. */
393 /* OR IF M COLS OF A HAVE BEEN TRIANGULARIZED. */
394
395 if (iz1 > iz2 || nsetp >= *m) {
396 goto L350;
397 }
398
399 /* COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W(). */
400
401 i__1 = iz2;
402 for (iz = iz1; iz <= i__1; ++iz) {
403 j = index[iz];
404 sm = 0.;
405 i__2 = *m;
406 for (l = npp1; l <= i__2; ++l) {
407 /* L40: */
408 sm += a[l + j * a_dim1] * b[l];
409 }
410 w[j] = sm;
411 /* L50: */
412 }
413 /* FIND LARGEST POSITIVE W(J). */
414 L60:
415 wmax = 0.;
416 i__1 = iz2;
417 for (iz = iz1; iz <= i__1; ++iz) {
418 j = index[iz];
419 if (w[j] > wmax) {
420 wmax = w[j];
421 izmax = iz;
422 }
423 /* L70: */
424 }
425
426 /* IF WMAX .LE. 0. GO TO TERMINATION. */
427 /* THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS. */
428
429 if (wmax <= 0.) {
430 goto L350;
431 }
432 iz = izmax;
433 j = index[iz];
434
435 /* THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P. */
436 /* BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID */
437 /* NEAR LINEAR DEPENDENCE. */
438
439 asave = a[npp1 + j * a_dim1];
440 i__1 = npp1 + 1;
441 h12_(&c__1, &npp1, &i__1, m, &a[j * a_dim1 + 1], &c__1, &up, &dummy, &
442 c__1, &c__1, &c__0);
443 unorm = 0.;
444 if (nsetp != 0) {
445 i__1 = nsetp;
446 for (l = 1; l <= i__1; ++l) {
447 /* L90: */
448 /* Computing 2nd power */
449 d__1 = a[l + j * a_dim1];
450 unorm += d__1 * d__1;
451 }
452 }
453 unorm = sqrt(unorm);
454 d__2 = unorm + (d__1 = a[npp1 + j * a_dim1], abs(d__1)) * .01;
455 if (diff_(&d__2, &unorm) > 0.) {
456
457 /* COL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE ZZ */
458 /* AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ). */
459
460 i__1 = *m;
461 for (l = 1; l <= i__1; ++l) {
462 /* L120: */
463 zz[l] = b[l];
464 }
465 i__1 = npp1 + 1;
466 h12_(&c__2, &npp1, &i__1, m, &a[j * a_dim1 + 1], &c__1, &up, &zz[1], &
467 c__1, &c__1, &c__1);
468 ztest = zz[npp1] / a[npp1 + j * a_dim1];
469
470 /* SEE IF ZTEST IS POSITIVE */
471
472 if (ztest > 0.) {
473 goto L140;
474 }
475 }
476
477 /* REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P. */
478 /* RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL */
479 /* COEFFS AGAIN. */
480
481 a[npp1 + j * a_dim1] = asave;
482 w[j] = 0.;
483 goto L60;
484
485 /* THE INDEX J=INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM */
486 /* SET Z TO SET P. UPDATE B, UPDATE INDICES, APPLY HOUSEHOLDER */
487 /* TRANSFORMATIONS TO COLS IN NEW SET Z, ZERO SUBDIAGONAL ELTS IN */
488 /* COL J, SET W(J)=0. */
489
490 L140:
491 i__1 = *m;
492 for (l = 1; l <= i__1; ++l) {
493 /* L150: */
494 b[l] = zz[l];
495 }
496
497 index[iz] = index[iz1];
498 index[iz1] = j;
499 ++iz1;
500 nsetp = npp1;
501 ++npp1;
502
503 if (iz1 <= iz2) {
504 i__1 = iz2;
505 for (jz = iz1; jz <= i__1; ++jz) {
506 jj = index[jz];
507 h12_(&c__2, &nsetp, &npp1, m, &a[j * a_dim1 + 1], &c__1, &up, &a[
508 jj * a_dim1 + 1], &c__1, mda, &c__1);
509 /* L160: */
510 }
511 }
512
513 if (nsetp != *m) {
514 i__1 = *m;
515 for (l = npp1; l <= i__1; ++l) {
516 /* L180: */
517 a[l + j * a_dim1] = 0.;
518 }
519 }
520
521 w[j] = 0.;
522 /* SOLVE THE TRIANGULAR SYSTEM. */
523 /* STORE THE SOLUTION TEMPORARILY IN ZZ(). */
524 rtnkey = 1;
525 goto L400;
526 L200:
527
528 /* ****** SECONDARY LOOP BEGINS HERE ****** */
529
530 /* ITERATION COUNTER. */
531
532 L210:
533 ++iter;
534 if (iter > itmax) {
535 *mode = 3;
536 /* write (*,'(/a)') ' NNLS quitting on iteration count.' */
537 goto L350;
538 }
539
540 /* SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE. */
541 /* IF NOT COMPUTE ALPHA. */
542
543 alpha = 2.;
544 i__1 = nsetp;
545 for (ip = 1; ip <= i__1; ++ip) {
546 l = index[ip];
547 if (zz[ip] <= 0.) {
548 t = -x[l] / (zz[ip] - x[l]);
549 if (alpha > t) {
550 alpha = t;
551 jj = ip;
552 }
553 }
554 /* L240: */
555 }
556
557 /* IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL */
558 /* STILL = 2. IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP. */
559
560 if (alpha == 2.) {
561 goto L330;
562 }
563
564 /* OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO */
565 /* INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ. */
566
567 i__1 = nsetp;
568 for (ip = 1; ip <= i__1; ++ip) {
569 l = index[ip];
570 x[l] += alpha * (zz[ip] - x[l]);
571 /* L250: */
572 }
573
574 /* MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I */
575 /* FROM SET P TO SET Z. */
576
577 i__ = index[jj];
578 L260:
579 x[i__] = 0.;
580
581 if (jj != nsetp) {
582 ++jj;
583 i__1 = nsetp;
584 for (j = jj; j <= i__1; ++j) {
585 ii = index[j];
586 index[j - 1] = ii;
587 g1_(&a[j - 1 + ii * a_dim1], &a[j + ii * a_dim1], &cc, &ss, &a[j
588 - 1 + ii * a_dim1]);
589 a[j + ii * a_dim1] = 0.;
590 i__2 = *n;
591 for (l = 1; l <= i__2; ++l) {
592 if (l != ii) {
593
594 /* Apply procedure G2 (CC,SS,A(J-1,L),A(J,L)) */
595
596 temp = a[j - 1 + l * a_dim1];
597 a[j - 1 + l * a_dim1] = cc * temp + ss * a[j + l * a_dim1]
598 ;
599 a[j + l * a_dim1] = -ss * temp + cc * a[j + l * a_dim1];
600 }
601 /* L270: */
602 }
603
604 /* Apply procedure G2 (CC,SS,B(J-1),B(J)) */
605
606 temp = b[j - 1];
607 b[j - 1] = cc * temp + ss * b[j];
608 b[j] = -ss * temp + cc * b[j];
609 /* L280: */
610 }
611 }
612
613 npp1 = nsetp;
614 --nsetp;
615 --iz1;
616 index[iz1] = i__;
617
618 /* SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE. THEY SHOULD */
619 /* BE BECAUSE OF THE WAY ALPHA WAS DETERMINED. */
620 /* IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR. ANY */
621 /* THAT ARE NONPOSITIVE WILL BE SET TO ZERO */
622 /* AND MOVED FROM SET P TO SET Z. */
623
624 i__1 = nsetp;
625 for (jj = 1; jj <= i__1; ++jj) {
626 i__ = index[jj];
627 if (x[i__] <= 0.) {
628 goto L260;
629 }
630 /* L300: */
631 }
632
633 /* COPY B( ) INTO ZZ( ). THEN SOLVE AGAIN AND LOOP BACK. */
634
635 i__1 = *m;
636 for (i__ = 1; i__ <= i__1; ++i__) {
637 /* L310: */
638 zz[i__] = b[i__];
639 }
640 rtnkey = 2;
641 goto L400;
642 L320:
643 goto L210;
644 /* ****** END OF SECONDARY LOOP ****** */
645
646 L330:
647 i__1 = nsetp;
648 for (ip = 1; ip <= i__1; ++ip) {
649 i__ = index[ip];
650 /* L340: */
651 x[i__] = zz[ip];
652 }
653 /* ALL NEW COEFFS ARE POSITIVE. LOOP BACK TO BEGINNING. */
654 goto L30;
655
656 /* ****** END OF MAIN LOOP ****** */
657
658 /* COME TO HERE FOR TERMINATION. */
659 /* COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR. */
660
661 L350:
662 sm = 0.;
663 if (npp1 <= *m) {
664 i__1 = *m;
665 for (i__ = npp1; i__ <= i__1; ++i__) {
666 /* L360: */
667 /* Computing 2nd power */
668 d__1 = b[i__];
669 sm += d__1 * d__1;
670 }
671 } else {
672 i__1 = *n;
673 for (j = 1; j <= i__1; ++j) {
674 /* L380: */
675 w[j] = 0.;
676 }
677 }
678 *rnorm = sqrt(sm);
679 return 0;
680
681 /* THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE */
682 /* TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ(). */
683
684 L400:
685 i__1 = nsetp;
686 for (l = 1; l <= i__1; ++l) {
687 ip = nsetp + 1 - l;
688 if (l != 1) {
689 i__2 = ip;
690 for (ii = 1; ii <= i__2; ++ii) {
691 zz[ii] -= a[ii + jj * a_dim1] * zz[ip + 1];
692 /* L410: */
693 }
694 }
695 jj = index[ip];
696 zz[ip] /= a[ip + jj * a_dim1];
697 /* L430: */
698 }
699 switch (rtnkey) {
700 case 1: goto L200;
701 case 2: goto L320;
702 }
703 return 0;
704 } /* nnls_ */
705