This subprogram solves a linearly constrained least squares
problem. Suppose there are given matrices E and A of
respective dimensions ME by N and MA by N, and vectors F
and B of respective lengths ME and MA. This subroutine
solves the problem
EX = F, (equations to be exactly satisfied)AX = B, (equations to be approximately satisfied, in the least squares sense)subject to components L+1,...,N nonnegative
Any values ME>=0, MA>=0 and 0<= L <=N are permitted.
The problem is reposed as problem DWNNLS
(WT*E)X = (WT*F)
( A) ( B), (least squares)
subject to components L+1,...,N nonnegative.
The subprogram chooses the heavy weight (or penalty parameter) WT.
The parameters for DWNNLS are
INPUT.. All TYPE REAL variables are DOUBLE PRECISION
W(*,*),MDW, The array W(*,*) is double subscripted with first
ME,MA,N,L dimensioning parameter equal to MDW. For this
discussion let us call M = ME + MA. Then MDW
must satisfy MDW>=M. The condition MDW<M
is an error.
The array W(*,*) contains the matrices and vectors
(E F)
(A B)
in rows and columns 1,...,M and 1,...,N+1
respectively. Columns 1,...,L correspond to
unconstrained variables X(1),...,X(L). The
remaining variables are constrained to be
nonnegative. The condition L<0 or L>N is
an error.
PRGOPT(*) This double precision array is the option vector.
If the user is satisfied with the nominal
subprogram features set
PRGOPT(1)=1 (or PRGOPT(1)=1.0)
Otherwise PRGOPT(*) is a linked list consisting of
groups of data of the following form
LINK
KEY
DATA SET
The parameters LINK and KEY are each one word.
The DATA SET can be comprised of several words.
The number of items depends on the value of KEY.
The value of LINK points to the first
entry of the next group of data within
PRGOPT(*). The exception is when there are
no more options to change. In that
case LINK=1 and the values KEY and DATA SET
are not referenced. The general layout of
PRGOPT(*) is as follows.
...PRGOPT(1)=LINK1 (link to first entry of next group)
. PRGOPT(2)=KEY1 (key to the option change)
. PRGOPT(3)=DATA VALUE (data value for this change)
. .
. .
. .
...PRGOPT(LINK1)=LINK2 (link to the first entry of
. next group)
. PRGOPT(LINK1+1)=KEY2 (key to the option change)
. PRGOPT(LINK1+2)=DATA VALUE
... .
. .
. .
...PRGOPT(LINK)=1 (no more options to change)
Values of LINK that are nonpositive are errors.
A value of LINK>NLINK=100000 is also an error.
This helps prevent using invalid but positive
values of LINK that will probably extend
beyond the program limits of PRGOPT(*).
Unrecognized values of KEY are ignored. The
order of the options is arbitrary and any number
of options can be changed with the following
restriction. To prevent cycling in the
processing of the option array a count of the
number of options changed is maintained.
Whenever this count exceeds NOPT=1000 an error
message is printed and the subprogram returns.
OPTIONS..
KEY=6
Scale the nonzero columns of the
entire data matrix
(E)
(A)
to have length one. The DATA SET for
this option is a single value. It must
be nonzero if unit length column scaling is
desired.
KEY=7
Scale columns of the entire data matrix
(E)
(A)
with a user-provided diagonal matrix.
The DATA SET for this option consists
of the N diagonal scaling factors, one for
each matrix column.
KEY=8
Change the rank determination tolerance from
the nominal value of SQRT(SRELPR). This quantity
can be no smaller than SRELPR, The arithmetic-
storage precision. The quantity used
here is internally restricted to be at
least SRELPR. The DATA SET for this option
is the new tolerance.
KEY=9
Change the blow-up parameter from the
nominal value of SQRT(SRELPR). The reciprocal of
this parameter is used in rejecting solution
components as too large when a variable is
first brought into the active set. Too large
means that the proposed component times the
reciprocal of the parameter is not less than
the ratio of the norms of the right-side
vector and the data matrix.
This parameter can be no smaller than SRELPR,
the arithmetic-storage precision.
For example, suppose we want to provide
a diagonal matrix to scale the problem
matrix and change the tolerance used for
determining linear dependence of dropped col
vectors. For these options the dimensions of
PRGOPT(*) must be at least N+6. The FORTRAN
statements defining these options would
be as follows.
PRGOPT(1)=N+3 (link to entry N+3 in PRGOPT(*))
PRGOPT(2)=7 (user-provided scaling key)
CALL DCOPY(N,D,1,PRGOPT(3),1) (copy the N
scaling factors from a user array called D(*)
into PRGOPT(3)-PRGOPT(N+2))
PRGOPT(N+3)=N+6 (link to entry N+6 of PRGOPT(*))
PRGOPT(N+4)=8 (linear dependence tolerance key)
PRGOPT(N+5)=... (new value of the tolerance)
PRGOPT(N+6)=1 (no more options to change)
IWORK(1), The amounts of working storage actually allocated
IWORK(2) for the working arrays WORK(*) and IWORK(*),
respectively. These quantities are compared with
the actual amounts of storage needed for DWNNLS( ).
Insufficient storage allocated for either WORK(*)
or IWORK(*) is considered an error. This feature
was included in DWNNLS( ) because miscalculating
the storage formulas for WORK(*) and IWORK(*)
might very well lead to subtle and hard-to-find
execution errors.
The length of WORK(*) must be at least
LW = ME+MA+5*N
This test will not be made if IWORK(1)<=0.
The length of IWORK(*) must be at least
LIW = ME+MA+N
This test will not be made if IWORK(2)<=0.
OUTPUT.. All TYPE REAL variables are DOUBLE PRECISION
X(*) An array dimensioned at least N, which will
contain the N components of the solution vector
on output.
RNORM The residual norm of the solution. The value of
RNORM contains the residual vector length of the
equality constraints and least squares equations.
MODE The value of MODE indicates the success or failure
of the subprogram.
MODE = 0 Subprogram completed successfully.
= 1 Max. number of iterations (equal to
3*(N-L)) exceeded. Nearly all problems
should complete in fewer than this
number of iterations. An approximate
solution and its corresponding residual
vector length are in X(*) and RNORM.
= 2 Usage error occurred. The offending
condition is noted with the error
processing subprogram, XERMSG( ).
User-designated
Working arrays..
WORK(*) A double precision working array of length at least
M + 5*N.
IWORK(*) An integer-valued working array of length at least
M+N.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp) | :: | w(mdw,*) | ||||
| integer | :: | mdw | ||||
| integer | :: | me | ||||
| integer | :: | ma | ||||
| integer | :: | n | ||||
| integer | :: | l | ||||
| real(kind=wp) | :: | prgopt(*) | ||||
| real(kind=wp) | :: | x(*) | ||||
| real(kind=wp) | :: | rnorm | ||||
| integer | :: | mode | ||||
| integer | :: | iwork(*) | ||||
| real(kind=wp) | :: | work(*) |
subroutine dwnnls (w, mdw, me, ma, n, l, prgopt, x, rnorm, mode, & iwork, work) integer :: iwork(*), l, l1, l2, l3, l4, l5, liw, lw, ma, mdw, me, & mode, n real(wp) :: prgopt(*), rnorm, w(mdw,*), work(*), x(*) character(len=8) :: xern1 mode = 0 if (ma+me<=0 .or. n<=0) return if (iwork(1)>0) then lw = me + ma + 5*n if (iwork(1)<lw) then write (xern1, '(I8)') lw write(*,*) 'INSUFFICIENT STORAGE ' // & 'ALLOCATED FOR WORK(*), NEED LW = ' // xern1 mode = 2 return endif endif if (iwork(2)>0) then liw = me + ma + n if (iwork(2)<liw) then write (xern1, '(I8)') liw write(*,*) 'INSUFFICIENT STORAGE ' // & 'ALLOCATED FOR IWORK(*), NEED LIW = ' // xern1 mode = 2 return endif endif if (mdw<me+ma) then write(*,*) 'THE VALUE MDW<ME+MA IS AN ERROR' mode = 2 return endif if (l<0 .or. l>n) then write(*,*) 'L>=0 .AND. L<=N IS REQUIRED' mode = 2 return endif ! THE PURPOSE OF THIS SUBROUTINE IS TO BREAK UP THE ARRAYS ! WORK(*) AND IWORK(*) INTO SEPARATE WORK ARRAYS ! REQUIRED BY THE MAIN SUBROUTINE DWNLSM( ). l1 = n + 1 l2 = l1 + n l3 = l2 + me + ma l4 = l3 + n l5 = l4 + n call dwnlsm(w, mdw, me, ma, n, l, prgopt, x, rnorm, mode, iwork, & iwork(l1), work(1), work(l1), work(l2), work(l3), & work(l4), work(l5)) end subroutine dwnnls