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