diff options
Diffstat (limited to 'Utilib/src/ALNNLS.f')
| -rw-r--r-- | Utilib/src/ALNNLS.f | 301 |
1 files changed, 301 insertions, 0 deletions
diff --git a/Utilib/src/ALNNLS.f b/Utilib/src/ALNNLS.f new file mode 100644 index 0000000..aa25a16 --- /dev/null +++ b/Utilib/src/ALNNLS.f @@ -0,0 +1,301 @@ +*DECK ALNNLS + SUBROUTINE ALNNLS(A,M,N,MP,NP,B,X,RNORM,MODE) +* +*----------------------------------------------------------------------- +* +*Purpose: +* implementation of the nonnegative least square algorithm. +* +*Author(s): A. Miller +* +*Parameters: input +* A input rectangular matrix. +* M,N first/second mathematical dimension of matrix A. +* MP,NP first/second physical dimension of matrix A. +* B source M-vector. +* +*Parameters: output +* A product matrix, Q*A, where Q IS AN M x M orthogonal matrix +* generated implicitly by this subroutine. +* B product Q*B. +* X solution N-vector. +* RNORM Euclidean norm of the residual vector. +* MODE success-failure flag (1: the solution has been computed +* successfully; 2: the dimensions of the problem are +* inconsistent, either m.LE.0 or n.LE.0; 3: iteration count +* exceeded with more than 10*n iterations). +* +*Reference: +* Chen, Donghui; Plemmons, Robert J. (2009). Nonnegativity constraints +* in numerical analysis. Symposium on the Birth of Numerical Analysis. +* +*----------------------------------------------------------------------- +* + IMPLICIT NONE +*---- +* Subroutine arguments +*---- + INTEGER M, N, MP, NP, MODE + DOUBLE PRECISION A(MP,NP), B(M), X(N), RNORM +*---- +* Local variables +*---- + INTEGER I, II, IP, ITER, ITMAX, IZ, IZ1, IZ2, IZMAX, J, JJ, JZ, + > L, MDA, NPP1, NSETP + DOUBLE PRECISION DUMMY(1), ALPHA, ASAVE, CC, FACTOR, SM, SS, T, + > TEMP, UNORM, UP, WMAX, ZTEST, XR, YR + INTEGER, ALLOCATABLE, DIMENSION(:) :: INDX + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: W, ZZ + PARAMETER(FACTOR = 0.01D0) +* + MODE = 1 + IF(M.LE.0 .OR. N.LE.0) THEN + MODE = 2 + RETURN + ENDIF + ITER = 0 + ITMAX = 10*N + ALLOCATE(INDX(N),W(N), ZZ(M)) +*---- +* Initialize the arrays INDX() and X(). +*---- + DO I = 1,N + X(I) = 0.0D0 + INDX(I) = I + ENDDO + IZ2 = N + IZ1 = 1 + NSETP = 0 + NPP1 = 1 +*---- +* ****** MAIN LOOP BEGINS HERE ****** +* Quit if all coefficients are already in the solution or if M cols of +* A have been triangularized. +*---- + 30 IF(IZ1.GT.IZ2 .OR. NSETP.GE.M) GO TO 350 +*---- +* Compute components of the dual (negative gradient) vector W. +*---- + DO IZ = IZ1,IZ2 + J = INDX(IZ) + W(J) = DOT_PRODUCT(A(NPP1:M,J), B(NPP1:M)) + ENDDO +*---- +* Find largest positive W(J). +*---- + 60 IZMAX = 0 + WMAX = 0.0D0 + DO IZ = IZ1,IZ2 + J = INDX(IZ) + IF(W(J).GT.WMAX) THEN + WMAX = W(J) + IZMAX = IZ + ENDIF + ENDDO +*---- +* If WMAX.le.0. go to termination. This indicates satisfaction of the +* Kuhn-Tucker conditions. +*---- + IF(WMAX.LE.0.0D0) GO TO 350 + IZ = IZMAX + J = INDX(IZ) +*---- +* The sign of W(J) is ok for J to be moved to set P. Begin the +* transformation and check new diagonal element to avoid near linear +* dependence. +*---- + ASAVE = A(NPP1,J) + CALL ALH12(1, NPP1, NPP1+1, M, A(:,J), UP, DUMMY(1), 1, 1, 0) + UNORM = 0.0D0 + IF(NSETP.NE.0) UNORM = SUM( A(1:NSETP,J)**2 ) + UNORM = SQRT(UNORM) + IF(UNORM + ABS(A(NPP1,J))*FACTOR - UNORM .GT. 0.0D0) THEN +*---- +* Col J is sufficiently independent. Copy B into ZZ, update ZZ +* and solve for ZTEST ( = proposed new value for X(J) ). +*---- + ZZ(1:M) = B(1:M) + CALL ALH12(2, NPP1, NPP1+1, M, A(:,J), UP, ZZ(1), 1, 1, 1) + ZTEST = ZZ(NPP1)/A(NPP1,J) +*---- +* See if ZTEST is positive. +*---- + IF(ZTEST.GT.0.0D0) GO TO 140 + ENDIF +*---- +* Reject J as a candidate to be moved from set Z to set P. Restore +* A(NPP1,J), set W(J) = 0., and loop back to test dual coeffs again. +*---- + A(NPP1,J) = ASAVE + W(J) = 0.0D0 + GO TO 60 +*---- +* The index J = INDX(IZ) has been selected to be moved from +* set Z to set P. Update B, update indices, apply Householder +* transformations to cols in new set Z, sero subdiagonal elts in +* col J, set W(J) = 0. +*---- + 140 B(1:M) = ZZ(1:M) + INDX(IZ) = INDX(IZ1) + INDX(IZ1) = J + IZ1 = IZ1+1 + NSETP = NPP1 + NPP1 = NPP1+1 + MDA = SIZE(A,1) + JJ = 0 + IF(IZ1.LE.IZ2) THEN + DO JZ = IZ1,IZ2 + JJ = INDX(JZ) + CALL ALH12(2, NSETP, NPP1, M, A(:,J), UP, A(1,JJ), 1, MDA, 1) + ENDDO + ENDIF + IF(NSETP.NE.M) A(NPP1:M,J) = 0.0D0 + W(J) = 0.0D0 +*---- +* Solve the triangular system. Store the solution temporarily in ZZ(). +*---- + DO L = 1, NSETP + IP = NSETP+1-L + IF(L .NE. 1) ZZ(1:IP) = ZZ(1:IP) - A(1:IP,JJ)*ZZ(IP+1) + JJ = INDX(IP) + ZZ(IP) = ZZ(IP) / A(IP,JJ) + ENDDO +*---- +* ****** SECONDARY LOOP BEGINS HERE ****** +*---- + 210 ITER = ITER+1 + IF(ITER.GT.ITMAX) THEN + MODE = 3 + WRITE (*,'(/A)') ' NNLS QUITTING ON ITERATION COUNT.' + GO TO 350 + ENDIF +*---- +* See if all new constrained coeffs are feasible; if not compute ALPHA. +*---- + ALPHA = 2.0D0 + DO IP = 1,NSETP + L = INDX(IP) + IF(ZZ(IP).LE.0.0D0) THEN + T = -X(L)/(ZZ(IP)-X(L)) + IF(ALPHA.GT.T) THEN + ALPHA = T + JJ = IP + ENDIF + ENDIF + ENDDO +*---- +* If all new constrained coeffs are feasible then ALPHA will still be +* equal to 2. If so exit from secondary loop to main loop. +*---- + IF(ALPHA == 2.0D0) GO TO 330 +*---- +* Otherwise use ALPHA which will be between 0. and 1. to interpolate +* between the old X and the new ZZ. +*---- + DO IP = 1,NSETP + L = INDX(IP) + X(L) = X(L) + ALPHA*(ZZ(IP)-X(L)) + ENDDO +*---- +* Modify A and B and the index arrays to move coefficient I from set +* P to set Z. +*---- + I = INDX(JJ) + 260 X(I) = 0.0D0 +*---- +* Compute.. matrix (C, S) so that (C, S)(A) = (SQRT(A**2+B**2)) +* (-S,C) (-S,C)(B) ( 0 ) +* Compute SIG = SQRT(A**2+B**2) +* SIG is computed last to allow for the possibility that SIG +* may be in the same location as A OR B . +*---- + IF(JJ.NE.NSETP) THEN + JJ = JJ+1 + DO J = JJ,NSETP + II = INDX(J) + INDX(J-1) = II + IF(ABS(A(J-1,II)).GT.ABS(A(J,II))) THEN + XR = A(J,II) / A(J-1,II) + YR = SQRT(1.0D0 + XR**2) + CC = SIGN(1.0D0/YR, A(J-1,II)) + SS = CC * XR + A(J-1,II) = ABS(A(J-1,II)) * YR + ELSE IF(A(J,II).NE.0.D0) THEN + XR = A(J-1,II) / A(J,II) + YR = SQRT(1.0D0 + XR**2) + SS = SIGN(1.0D0/YR, A(J,II)) + CC = SS * XR + A(J-1,II) = ABS(A(J,II)) * YR + ELSE + CC = 0.0D0 + SS = 1.0D0 + ENDIF + A(J,II) = 0.0D0 + DO L = 1,N + IF(L.NE.II) THEN +*---- +* Apply procedure G2 (CC,SS,A(J-1,L),A(J,L)) +*---- + TEMP = A(J-1,L) + A(J-1,L) = CC*TEMP + SS*A(J,L) + A(J,L) = -SS*TEMP + CC*A(J,L) + ENDIF + ENDDO +*---- +* Apply procedure G2 (CC,SS,B(J-1),B(J)) +*---- + TEMP = B(J-1) + B(J-1) = CC*TEMP + SS*B(J) + B(J) = -SS*TEMP + CC*B(J) + ENDDO + ENDIF + NPP1 = NSETP + NSETP = NSETP-1 + IZ1 = IZ1-1 + INDX(IZ1) = I +*---- +* See if the remaining coeffs in set P are feasible. They should be +* because of the way ALPHA was determined. if any are infeasible it is +* due to round-off error. any that are nonpositive will be set to 0.0d0 +* and moved from set P to set Z. +*---- + DO JJ = 1,NSETP + I = INDX(JJ) + IF(X(I).LE.0.0D0) GO TO 260 + ENDDO +*---- +* Copy B( ) into ZZ( ). Then solve again and loop back. +*---- + ZZ(1:M) = B(1:M) + DO L = 1, NSETP + IP = NSETP+1-L + IF(L .NE. 1) ZZ(1:IP) = ZZ(1:IP) - A(1:IP,JJ)*ZZ(IP+1) + JJ = INDX(IP) + ZZ(IP) = ZZ(IP) / A(IP,JJ) + ENDDO + GO TO 210 +*---- +* ****** End of secondary loop ****** +*---- + 330 DO IP = 1,NSETP + I = INDX(IP) + X(I) = ZZ(IP) + ENDDO +*---- +* All new coeffs are positive. Loop back to beginning. +*---- + GO TO 30 +*---- +* Come to here for termination. Compute the norm of the final residual +* vector. +*---- + 350 SM = 0.0D0 + IF(NPP1.LE.M) THEN + SM = SUM( B(NPP1:M)**2 ) + ELSE + W(1:N) = 0.0D0 + ENDIF + RNORM = SQRT(SM) + DEALLOCATE(ZZ,W,INDX) + RETURN + END |