1*DECK DWNNLS
2      SUBROUTINE DWNNLS (W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE,
3     +   IWORK, WORK)
4C***BEGIN PROLOGUE  DWNNLS
5C***PURPOSE  Solve a linearly constrained least squares problem with
6C            equality constraints and nonnegativity constraints on
7C            selected variables.
8C***LIBRARY   SLATEC
9C***CATEGORY  K1A2A
10C***TYPE      DOUBLE PRECISION (WNNLS-S, DWNNLS-D)
11C***KEYWORDS  CONSTRAINED LEAST SQUARES, CURVE FITTING, DATA FITTING,
12C             EQUALITY CONSTRAINTS, INEQUALITY CONSTRAINTS,
13C             NONNEGATIVITY CONSTRAINTS, QUADRATIC PROGRAMMING
14C***AUTHOR  Hanson, R. J., (SNLA)
15C           Haskell, K. H., (SNLA)
16C***DESCRIPTION
17C
18C     Abstract
19C
20C     This subprogram solves a linearly constrained least squares
21C     problem.  Suppose there are given matrices E and A of
22C     respective dimensions ME by N and MA by N, and vectors F
23C     and B of respective lengths ME and MA.  This subroutine
24C     solves the problem
25C
26C               EX = F, (equations to be exactly satisfied)
27C
28C               AX = B, (equations to be approximately satisfied,
29C                        in the least squares sense)
30C
31C               subject to components L+1,...,N nonnegative
32C
33C     Any values ME.GE.0, MA.GE.0 and 0.LE. L .LE.N are permitted.
34C
35C     The problem is reposed as problem DWNNLS
36C
37C               (WT*E)X = (WT*F)
38C               (   A)    (   B), (least squares)
39C               subject to components L+1,...,N nonnegative.
40C
41C     The subprogram chooses the heavy weight (or penalty parameter) WT.
42C
43C     The parameters for DWNNLS are
44C
45C     INPUT.. All TYPE REAL variables are DOUBLE PRECISION
46C
47C     W(*,*),MDW,  The array W(*,*) is double subscripted with first
48C     ME,MA,N,L    dimensioning parameter equal to MDW.  For this
49C                  discussion let us call M = ME + MA.  Then MDW
50C                  must satisfy MDW.GE.M.  The condition MDW.LT.M
51C                  is an error.
52C
53C                  The array W(*,*) contains the matrices and vectors
54C
55C                       (E  F)
56C                       (A  B)
57C
58C                  in rows and columns 1,...,M and 1,...,N+1
59C                  respectively.  Columns 1,...,L correspond to
60C                  unconstrained variables X(1),...,X(L).  The
61C                  remaining variables are constrained to be
62C                  nonnegative. The condition L.LT.0 or L.GT.N is
63C                  an error.
64C
65C     PRGOPT(*)    This double precision array is the option vector.
66C                  If the user is satisfied with the nominal
67C                  subprogram features set
68C
69C                  PRGOPT(1)=1 (or PRGOPT(1)=1.0)
70C
71C                  Otherwise PRGOPT(*) is a linked list consisting of
72C                  groups of data of the following form
73C
74C                  LINK
75C                  KEY
76C                  DATA SET
77C
78C                  The parameters LINK and KEY are each one word.
79C                  The DATA SET can be comprised of several words.
80C                  The number of items depends on the value of KEY.
81C                  The value of LINK points to the first
82C                  entry of the next group of data within
83C                  PRGOPT(*).  The exception is when there are
84C                  no more options to change.  In that
85C                  case LINK=1 and the values KEY and DATA SET
86C                  are not referenced. The general layout of
87C                  PRGOPT(*) is as follows.
88C
89C               ...PRGOPT(1)=LINK1 (link to first entry of next group)
90C               .  PRGOPT(2)=KEY1 (key to the option change)
91C               .  PRGOPT(3)=DATA VALUE (data value for this change)
92C               .       .
93C               .       .
94C               .       .
95C               ...PRGOPT(LINK1)=LINK2 (link to the first entry of
96C               .                       next group)
97C               .  PRGOPT(LINK1+1)=KEY2 (key to the option change)
98C               .  PRGOPT(LINK1+2)=DATA VALUE
99C               ...     .
100C               .       .
101C               .       .
102C               ...PRGOPT(LINK)=1 (no more options to change)
103C
104C                  Values of LINK that are nonpositive are errors.
105C                  A value of LINK.GT.NLINK=100000 is also an error.
106C                  This helps prevent using invalid but positive
107C                  values of LINK that will probably extend
108C                  beyond the program limits of PRGOPT(*).
109C                  Unrecognized values of KEY are ignored.  The
110C                  order of the options is arbitrary and any number
111C                  of options can be changed with the following
112C                  restriction.  To prevent cycling in the
113C                  processing of the option array a count of the
114C                  number of options changed is maintained.
115C                  Whenever this count exceeds NOPT=1000 an error
116C                  message is printed and the subprogram returns.
117C
118C                  OPTIONS..
119C
120C                  KEY=6
121C                         Scale the nonzero columns of the
122C                  entire data matrix
123C                  (E)
124C                  (A)
125C                  to have length one. The DATA SET for
126C                  this option is a single value.  It must
127C                  be nonzero if unit length column scaling is
128C                  desired.
129C
130C                  KEY=7
131C                         Scale columns of the entire data matrix
132C                  (E)
133C                  (A)
134C                  with a user-provided diagonal matrix.
135C                  The DATA SET for this option consists
136C                  of the N diagonal scaling factors, one for
137C                  each matrix column.
138C
139C                  KEY=8
140C                         Change the rank determination tolerance from
141C                  the nominal value of SQRT(SRELPR).  This quantity
142C                  can be no smaller than SRELPR, The arithmetic-
143C                  storage precision.  The quantity used
144C                  here is internally restricted to be at
145C                  least SRELPR.  The DATA SET for this option
146C                  is the new tolerance.
147C
148C                  KEY=9
149C                         Change the blow-up parameter from the
150C                  nominal value of SQRT(SRELPR).  The reciprocal of
151C                  this parameter is used in rejecting solution
152C                  components as too large when a variable is
153C                  first brought into the active set.  Too large
154C                  means that the proposed component times the
155C                  reciprocal of the parameter is not less than
156C                  the ratio of the norms of the right-side
157C                  vector and the data matrix.
158C                  This parameter can be no smaller than SRELPR,
159C                  the arithmetic-storage precision.
160C
161C                  For example, suppose we want to provide
162C                  a diagonal matrix to scale the problem
163C                  matrix and change the tolerance used for
164C                  determining linear dependence of dropped col
165C                  vectors.  For these options the dimensions of
166C                  PRGOPT(*) must be at least N+6.  The FORTRAN
167C                  statements defining these options would
168C                  be as follows.
169C
170C                  PRGOPT(1)=N+3 (link to entry N+3 in PRGOPT(*))
171C                  PRGOPT(2)=7 (user-provided scaling key)
172C
173C                  CALL DCOPY(N,D,1,PRGOPT(3),1) (copy the N
174C                  scaling factors from a user array called D(*)
175C                  into PRGOPT(3)-PRGOPT(N+2))
176C
177C                  PRGOPT(N+3)=N+6 (link to entry N+6 of PRGOPT(*))
178C                  PRGOPT(N+4)=8 (linear dependence tolerance key)
179C                  PRGOPT(N+5)=... (new value of the tolerance)
180C
181C                  PRGOPT(N+6)=1 (no more options to change)
182C
183C
184C     IWORK(1),    The amounts of working storage actually allocated
185C     IWORK(2)     for the working arrays WORK(*) and IWORK(*),
186C                  respectively.  These quantities are compared with
187C                  the actual amounts of storage needed for DWNNLS( ).
188C                  Insufficient storage allocated for either WORK(*)
189C                  or IWORK(*) is considered an error.  This feature
190C                  was included in DWNNLS( ) because miscalculating
191C                  the storage formulas for WORK(*) and IWORK(*)
192C                  might very well lead to subtle and hard-to-find
193C                  execution errors.
194C
195C                  The length of WORK(*) must be at least
196C
197C                  LW = ME+MA+5*N
198C                  This test will not be made if IWORK(1).LE.0.
199C
200C                  The length of IWORK(*) must be at least
201C
202C                  LIW = ME+MA+N
203C                  This test will not be made if IWORK(2).LE.0.
204C
205C     OUTPUT.. All TYPE REAL variables are DOUBLE PRECISION
206C
207C     X(*)         An array dimensioned at least N, which will
208C                  contain the N components of the solution vector
209C                  on output.
210C
211C     RNORM        The residual norm of the solution.  The value of
212C                  RNORM contains the residual vector length of the
213C                  equality constraints and least squares equations.
214C
215C     MODE         The value of MODE indicates the success or failure
216C                  of the subprogram.
217C
218C                  MODE = 0  Subprogram completed successfully.
219C
220C                       = 1  Max. number of iterations (equal to
221C                            3*(N-L)) exceeded. Nearly all problems
222C                            should complete in fewer than this
223C                            number of iterations. An approximate
224C                            solution and its corresponding residual
225C                            vector length are in X(*) and RNORM.
226C
227C                       = 2  Usage error occurred.  The offending
228C                            condition is noted with the error
229C                            processing subprogram, XERMSG( ).
230C
231C     User-designated
232C     Working arrays..
233C
234C     WORK(*)      A double precision working array of length at least
235C                  M + 5*N.
236C
237C     IWORK(*)     An integer-valued working array of length at least
238C                  M+N.
239C
240C***REFERENCES  K. H. Haskell and R. J. Hanson, An algorithm for
241C                 linear least squares problems with equality and
242C                 nonnegativity constraints, Report SAND77-0552, Sandia
243C                 Laboratories, June 1978.
244C               K. H. Haskell and R. J. Hanson, Selected algorithms for
245C                 the linearly constrained least squares problem - a
246C                 users guide, Report SAND78-1290, Sandia Laboratories,
247C                 August 1979.
248C               K. H. Haskell and R. J. Hanson, An algorithm for
249C                 linear least squares problems with equality and
250C                 nonnegativity constraints, Mathematical Programming
251C                 21 (1981), pp. 98-118.
252C               R. J. Hanson and K. H. Haskell, Two algorithms for the
253C                 linearly constrained least squares problem, ACM
254C                 Transactions on Mathematical Software, September 1982.
255C               C. L. Lawson and R. J. Hanson, Solving Least Squares
256C                 Problems, Prentice-Hall, Inc., 1974.
257C***ROUTINES CALLED  DWNLSM, XERMSG
258C***REVISION HISTORY  (YYMMDD)
259C   790701  DATE WRITTEN
260C   890531  Changed all specific intrinsics to generic.  (WRB)
261C   890618  Completely restructured and revised.  (WRB & RWC)
262C   891006  REVISION DATE from Version 3.2
263C   891214  Prologue converted to Version 4.0 format.  (BAB)
264C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
265C   900510  Convert XERRWV calls to XERMSG calls, change Prologue
266C           comments to agree with WNNLS.  (RWC)
267C   920501  Reformatted the REFERENCES section.  (WRB)
268C***END PROLOGUE  DWNNLS
269      INTEGER IWORK(*), L, L1, L2, L3, L4, L5, LIW, LW, MA, MDW, ME,
270     *     MODE, N
271      DOUBLE PRECISION  PRGOPT(*), RNORM, W(MDW,*), WORK(*), X(*)
272      CHARACTER*8 XERN1
273C***FIRST EXECUTABLE STATEMENT  DWNNLS
274      MODE = 0
275      IF (MA+ME.LE.0 .OR. N.LE.0) RETURN
276C
277      IF (IWORK(1).GT.0) THEN
278         LW = ME + MA + 5*N
279         IF (IWORK(1).LT.LW) THEN
280            WRITE (XERN1, '(I8)') LW
281            CALL XERMSG ('SLATEC', 'DWNNLS', 'INSUFFICIENT STORAGE ' //
282     *         'ALLOCATED FOR WORK(*), NEED LW = ' // XERN1, 2, 1)
283            MODE = 2
284            RETURN
285         ENDIF
286      ENDIF
287C
288      IF (IWORK(2).GT.0) THEN
289         LIW = ME + MA + N
290         IF (IWORK(2).LT.LIW) THEN
291            WRITE (XERN1, '(I8)') LIW
292            CALL XERMSG ('SLATEC', 'DWNNLS', 'INSUFFICIENT STORAGE ' //
293     *         'ALLOCATED FOR IWORK(*), NEED LIW = ' // XERN1, 2, 1)
294            MODE = 2
295            RETURN
296         ENDIF
297      ENDIF
298C
299      IF (MDW.LT.ME+MA) THEN
300         CALL XERMSG ('SLATEC', 'DWNNLS',
301     *      'THE VALUE MDW.LT.ME+MA IS AN ERROR', 1, 1)
302         MODE = 2
303         RETURN
304      ENDIF
305C
306      IF (L.LT.0 .OR. L.GT.N) THEN
307         CALL XERMSG ('SLATEC', 'DWNNLS',
308     *      'L.GE.0 .AND. L.LE.N IS REQUIRED', 2, 1)
309         MODE = 2
310         RETURN
311      ENDIF
312C
313C     THE PURPOSE OF THIS SUBROUTINE IS TO BREAK UP THE ARRAYS
314C     WORK(*) AND IWORK(*) INTO SEPARATE WORK ARRAYS
315C     REQUIRED BY THE MAIN SUBROUTINE DWNLSM( ).
316C
317      L1 = N + 1
318      L2 = L1 + N
319      L3 = L2 + ME + MA
320      L4 = L3 + N
321      L5 = L4 + N
322C
323      CALL DWNLSM(W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE, IWORK,
324     *            IWORK(L1), WORK(1), WORK(L1), WORK(L2), WORK(L3),
325     *            WORK(L4), WORK(L5))
326      RETURN
327      END
328