This is a companion subprogram to DWNNLS. The documentation for DWNNLS has complete usage instructions.
In addition to the parameters discussed in the prologue to subroutine DWNNLS, the following work arrays are used in subroutine DWNLSM (they are passed through the calling sequence from DWNNLS for purposes of variable dimensioning). Their contents will in general be of no interest to the user.
IPIVOT(*)
An array of length N. Upon completion it contains the
pivoting information for the cols of W(*,*).
ITYPE(*)
An array of length M which is used to keep track
of the classification of the equations. ITYPE(I)=0
denotes equation I as an equality constraint.
ITYPE(I)=1 denotes equation I as a least squares
equation.
WD(*)
An array of length N. Upon completion it contains the
dual solution vector.
H(*)
An array of length N. Upon completion it contains the
pivot scalars of the Householder transformations performed
in the case KRANK<L.
SCALE(*)
An array of length M which is used by the subroutine
to store the diagonal matrix of weights.
These are used to apply the modified Givens
transformations.
Z(*),TEMP(*)
Working arrays of length N.
D(*)
An array of length N that contains the
column scaling for the matrix (E).
(A)
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp) | :: | w(mdw,*) | ||||
| integer | :: | mdw | ||||
| integer | :: | mme | ||||
| integer | :: | ma | ||||
| integer | :: | n | ||||
| integer | :: | l | ||||
| real(kind=wp) | :: | prgopt(*) | ||||
| real(kind=wp) | :: | x(*) | ||||
| real(kind=wp) | :: | rnorm | ||||
| integer | :: | mode | ||||
| integer | :: | ipivot(*) | ||||
| integer | :: | itype(*) | ||||
| real(kind=wp) | :: | wd(*) | ||||
| real(kind=wp) | :: | h(*) | ||||
| real(kind=wp) | :: | scale(*) | ||||
| real(kind=wp) | :: | z(*) | ||||
| real(kind=wp) | :: | temp(*) | ||||
| real(kind=wp) | :: | d(*) |
subroutine dwnlsm (w, mdw, mme, ma, n, l, prgopt, x, rnorm, mode, & ipivot, itype, wd, h, scale, z, temp, d) integer :: ipivot(*), itype(*), l, ma, mdw, mme, mode, n real(wp) :: d(*), h(*), prgopt(*), rnorm, scale(*), temp(*), & w(mdw,*), wd(*), x(*), z(*) real(wp) :: alamda, alpha, alsq, amax, blowup, bnorm, & dope(3), eanorm, fac, sm, sparam(5), t, tau, wmax, z2, & zz integer :: i, idope(3), imax, isol, itemp, iter, itmax, iwmax, j, & jcon, jp, key, krank, l1, last, link, m, me, next, niv, nlink, & nopt, nsoln, ntimes logical :: done, feasbl, hitcon, pos ! Set the nominal tolerance used in the code. tau = sqrt(drelpr) m = ma + mme me = mme mode = 2 ! To process option vector fac = 1.0e-4_wp ! Set the nominal blow up factor used in the code. blowup = tau ! The nominal column scaling used in the code is ! the identity scaling. call dcopy (n, [1.0_wp], 0, d, 1) ! Define bound for number of options to change. nopt = 1000 ! Define bound for positive value of LINK. nlink = 100000 ntimes = 0 last = 1 link = prgopt(1) if (link<=0 .or. link>nlink) then write(*,*) 'IN DWNNLS, THE OPTION VECTOR IS UNDEFINED' return endif do if (link<=1) exit ntimes = ntimes + 1 if (ntimes>nopt) then write(*,*) 'IN DWNNLS, THE LINKS IN THE OPTION VECTOR ARE CYCLING.' return endif key = prgopt(last+1) if (key==6 .and. prgopt(last+2)/=0.0_wp) then do j = 1,n t = dnrm2(m,w(1,j),1) if (t/=0.0_wp) t = 1.0_wp/t d(j) = t end do endif if (key==7) call dcopy (n, prgopt(last+2), 1, d, 1) if (key==8) tau = max(drelpr,prgopt(last+2)) if (key==9) blowup = max(drelpr,prgopt(last+2)) next = prgopt(link) if (next<=0 .or. next>nlink) then write(*,*) 'IN DWNNLS, THE OPTION VECTOR IS UNDEFINED' return endif last = link link = next end do do j = 1,n call dscal (m, d(j), w(1,j), 1) end do ! Process option vector done = .false. iter = 0 itmax = 3*(n-l) mode = 0 nsoln = l l1 = min(m,l) ! Compute scale factor to apply to equality constraint equations. do j = 1,n wd(j) = dasum(m,w(1,j),1) end do imax = idamax(n,wd,1) eanorm = wd(imax) bnorm = dasum(m,w(1,n+1),1) alamda = eanorm/(drelpr*fac) ! On machines, such as the VAXes using D floating, with a very ! limited exponent range for double precision values, the previously ! computed value of ALAMDA may cause an overflow condition. ! Therefore, this code further limits the value of ALAMDA. alamda = min(alamda,sqrt(huge(1.0_wp))) ! Define scaling diagonal matrix for modified Givens usage and ! classify equation types. alsq = alamda**2 do i = 1,m ! When equation I is heavily weighted ITYPE(I)=0, ! else ITYPE(I)=1. if (i<=me) then t = alsq itemp = 0 else t = 1.0_wp itemp = 1 endif scale(i) = t itype(i) = itemp end do ! Set the solution vector X(*) to zero and the column interchange ! matrix to the identity. call dcopy (n, [0.0_wp], 0, x, 1) do i = 1,n ipivot(i) = i end do ! Perform initial triangularization in the submatrix ! corresponding to the unconstrained variables. ! Set first L components of dual vector to zero because ! these correspond to the unconstrained variables. call dcopy (l, [0.0_wp], 0, wd, 1) ! The arrays IDOPE(*) and DOPE(*) are used to pass ! information to DWNLIT(). This was done to avoid ! a long calling sequence or the use of COMMON. idope(1) = me idope(2) = nsoln idope(3) = l1 dope(1) = alsq dope(2) = eanorm dope(3) = tau call dwnlit (w, mdw, m, n, l, ipivot, itype, h, scale, rnorm, & idope, dope, done) me = idope(1) krank = idope(2) niv = idope(3) main : do ! Perform WNNLS algorithm using the following steps. ! ! Until(DONE) ! compute search direction and feasible point ! when (HITCON) add constraints ! else perform multiplier test and drop a constraint ! fin ! Compute-Final-Solution ! ! To compute search direction and feasible point, ! solve the triangular system of currently non-active ! variables and store the solution in Z(*). ! ! To solve system ! Copy right hand side into TEMP vector to use overwriting method. if (done) exit main isol = l + 1 if (nsoln>=isol) then call dcopy (niv, w(1,n+1), 1, temp, 1) do j = nsoln,isol,-1 if (j>krank) then i = niv - nsoln + j else i = j endif if (j>krank .and. j<=l) then z(j) = 0.0_wp else z(j) = temp(i)/w(i,j) call daxpy (i-1, -z(j), w(1,j), 1, temp, 1) endif end do endif ! Increment iteration counter and check against maximum number ! of iterations. iter = iter + 1 if (iter>itmax) then mode = 1 done = .true. endif ! Check to see if any constraints have become active. ! If so, calculate an interpolation factor so that all ! active constraints are removed from the basis. alpha = 2.0_wp hitcon = .false. do j = l+1,nsoln zz = z(j) if (zz<=0.0_wp) then t = x(j)/(x(j)-zz) if (t<alpha) then alpha = t jcon = j endif hitcon = .true. endif end do ! Compute search direction and feasible point if (hitcon) then ! To add constraints, use computed ALPHA to interpolate between ! last feasible solution X(*) and current unconstrained (and ! infeasible) solution Z(*). do j = l+1,nsoln x(j) = x(j) + alpha*(z(j)-x(j)) end do feasbl = .false. do ! Remove column JCON and shift columns JCON+1 through N to the ! left. Swap column JCON into the N th position. This achieves ! upper Hessenberg form for the nonactive constraints and ! leaves an upper Hessenberg matrix to retriangularize. do i = 1,m t = w(i,jcon) call dcopy (n-jcon, w(i, jcon+1), mdw, w(i, jcon), mdw) w(i,n) = t end do ! Update permuted index vector to reflect this shift and swap. itemp = ipivot(jcon) do i = jcon,n - 1 ipivot(i) = ipivot(i+1) end do ipivot(n) = itemp ! Similarly permute X(*) vector. call dcopy (n-jcon, x(jcon+1), 1, x(jcon), 1) x(n) = 0.0_wp nsoln = nsoln - 1 niv = niv - 1 ! Retriangularize upper Hessenberg matrix after adding ! constraints. i = krank + jcon - l do j = jcon,nsoln if (itype(i)==0 .and. itype(i+1)==0) then ! Zero IP1 to I in column J if (w(i+1,j)/=0.0_wp) then call drotmg (scale(i), scale(i+1), w(i,j), w(i+1,j), & sparam) w(i+1,j) = 0.0_wp call drotm (n+1-j, w(i,j+1), mdw, w(i+1,j+1), mdw, & sparam) endif elseif (itype(i)==1 .and. itype(i+1)==1) then ! Zero IP1 to I in column J if (w(i+1,j)/=0.0_wp) then call drotmg (scale(i), scale(i+1), w(i,j), w(i+1,j), & sparam) w(i+1,j) = 0.0_wp call drotm (n+1-j, w(i,j+1), mdw, w(i+1,j+1), mdw, & sparam) endif elseif (itype(i)==1 .and. itype(i+1)==0) then call dswap (n+1, w(i,1), mdw, w(i+1,1), mdw) call dswap (1, scale(i), 1, scale(i+1), 1) itemp = itype(i+1) itype(i+1) = itype(i) itype(i) = itemp ! Swapped row was formerly a pivot element, so it will ! be large enough to perform elimination. ! Zero IP1 to I in column J. if (w(i+1,j)/=0.0_wp) then call drotmg (scale(i), scale(i+1), w(i,j), w(i+1,j), & sparam) w(i+1,j) = 0.0_wp call drotm (n+1-j, w(i,j+1), mdw, w(i+1,j+1), mdw, & sparam) endif elseif (itype(i)==0 .and. itype(i+1)==1) then if (scale(i)*w(i,j)**2/alsq>(tau*eanorm)**2) then ! Zero IP1 to I in column J if (w(i+1,j)/=0.0_wp) then call drotmg (scale(i), scale(i+1), w(i,j), & w(i+1,j), sparam) w(i+1,j) = 0.0_wp call drotm (n+1-j, w(i,j+1), mdw, w(i+1,j+1), mdw, & sparam) endif else call dswap (n+1, w(i,1), mdw, w(i+1,1), mdw) call dswap (1, scale(i), 1, scale(i+1), 1) itemp = itype(i+1) itype(i+1) = itype(i) itype(i) = itemp w(i+1,j) = 0.0_wp endif endif i = i + 1 end do ! See if the remaining coefficients in the solution set are ! feasible. They should be because of the way ALPHA was ! determined. If any are infeasible, it is due to roundoff ! error. Any that are non-positive will be set to zero and ! removed from the solution set. do jcon = l+1,nsoln if (x(jcon)<=0.0_wp) then exit else if (jcon==nsoln) feasbl = .true. end if end do if (feasbl) exit end do else ! To perform multiplier test and drop a constraint. call dcopy (nsoln, z, 1, x, 1) if (nsoln<n) call dcopy (n-nsoln, [0.0_wp], 0, x(nsoln+1), 1) ! Reclassify least squares equations as equalities as necessary. i = niv + 1 do if (i>me) exit if (itype(i)==0) then i = i + 1 else call dswap (n+1, w(i,1), mdw, w(me,1), mdw) call dswap (1, scale(i), 1, scale(me), 1) itemp = itype(i) itype(i) = itype(me) itype(me) = itemp me = me - 1 endif end do ! Form inner product vector WD(*) of dual coefficients. do j = nsoln+1,n sm = 0.0_wp do i = nsoln+1,m sm = sm + scale(i)*w(i,j)*w(i,n+1) end do wd(j) = sm end do do ! Find J such that WD(J)=WMAX is maximum. This determines ! that the incoming column J will reduce the residual vector ! and be positive. wmax = 0.0_wp iwmax = nsoln + 1 do j = nsoln+1,n if (wd(j)>wmax) then wmax = wd(j) iwmax = j endif end do if (wmax<=0.0_wp) exit main ! Set dual coefficients to zero for incoming column. wd(iwmax) = 0.0_wp ! WMAX > 0.0_wp, so okay to move column IWMAX to solution set. ! Perform transformation to retriangularize, and test for near ! linear dependence. ! ! Swap column IWMAX into NSOLN-th position to maintain upper ! Hessenberg form of adjacent columns, and add new column to ! triangular decomposition. nsoln = nsoln + 1 niv = niv + 1 if (nsoln/=iwmax) then call dswap (m, w(1,nsoln), 1, w(1,iwmax), 1) wd(iwmax) = wd(nsoln) wd(nsoln) = 0.0_wp itemp = ipivot(nsoln) ipivot(nsoln) = ipivot(iwmax) ipivot(iwmax) = itemp endif ! Reduce column NSOLN so that the matrix of nonactive constraints ! variables is triangular. do j = m,niv+1,-1 jp = j - 1 ! When operating near the ME line, test to see if the pivot ! element is near zero. If so, use the largest element above ! it as the pivot. This is to maintain the sharp interface ! between weighted and non-weighted rows in all cases. if (j==me+1) then imax = me amax = scale(me)*w(me,nsoln)**2 do jp = j - 1,niv,-1 t = scale(jp)*w(jp,nsoln)**2 if (t>amax) then imax = jp amax = t endif end do jp = imax endif if (w(j,nsoln)/=0.0_wp) then call drotmg (scale(jp), scale(j), w(jp,nsoln), w(j,nsoln), sparam) w(j,nsoln) = 0.0_wp call drotm (n+1-nsoln, w(jp,nsoln+1), mdw, w(j,nsoln+1), mdw, sparam) endif end do ! Solve for Z(NSOLN)=proposed new value for X(NSOLN). Test if ! this is nonpositive or too large. If this was true or if the ! pivot term was zero, reject the column as dependent. if (w(niv,nsoln)/=0.0_wp) then isol = niv z2 = w(isol,n+1)/w(isol,nsoln) z(nsoln) = z2 pos = z2 > 0.0_wp if (z2*eanorm>=bnorm .and. pos) then pos = .not. (blowup*z2*eanorm>=bnorm) endif elseif (niv<=me .and. w(me+1,nsoln)/=0.0_wp) then ! Try to add row ME+1 as an additional equality constraint. ! Check size of proposed new solution component. ! Reject it if it is too large. isol = me + 1 if (pos) then ! Swap rows ME+1 and NIV, and scale factors for these rows. call dswap (n+1, w(me+1,1), mdw, w(niv,1), mdw) call dswap (1, scale(me+1), 1, scale(niv), 1) itemp = itype(me+1) itype(me+1) = itype(niv) itype(niv) = itemp me = me + 1 endif else pos = .false. endif if (.not.pos) then nsoln = nsoln - 1 niv = niv - 1 endif if (pos.or.done) exit end do endif end do main ! Else perform multiplier test and drop a constraint. To compute ! final solution. Solve system, store results in X(*). ! ! Copy right hand side into TEMP vector to use overwriting method. isol = 1 if (nsoln>=isol) then call dcopy (niv, w(1,n+1), 1, temp, 1) do j = nsoln,isol,-1 if (j>krank) then i = niv - nsoln + j else i = j endif if (j>krank .and. j<=l) then z(j) = 0.0_wp else z(j) = temp(i)/w(i,j) call daxpy (i-1, -z(j), w(1,j), 1, temp, 1) endif end do endif ! Solve system. call dcopy (nsoln, z, 1, x, 1) ! Apply Householder transformations to X(*) if KRANK<L if (krank<l) then do i = 1,krank call dh12 (2, i, krank+1, l, w(i,1), mdw, h(i), x, 1, 1, 1) end do endif ! Fill in trailing zeroes for constrained variables not in solution. if (nsoln<n) call dcopy (n-nsoln, [0.0_wp], 0, x(nsoln+1), 1) ! Permute solution vector to natural order. do i = 1,n j = i do if (ipivot(j)==i) exit j = j + 1 end do ipivot(j) = ipivot(i) ipivot(i) = j call dswap (1, x(j), 1, x(i), 1) end do ! Rescale the solution using the column scaling. do j = 1,n x(j) = x(j)*d(j) end do do i = nsoln+1,m t = w(i,n+1) if (i<=me) t = t/alamda t = (scale(i)*t)*t rnorm = rnorm + t end do rnorm = sqrt(rnorm) end subroutine dwnlsm