Mercurial > octave-nkf
view libcruft/npsol/lsgset.f @ 2329:30c606bec7a8
[project @ 1996-07-19 01:29:05 by jwe]
Initial revision
author | jwe |
---|---|
date | Fri, 19 Jul 1996 01:29:55 +0000 |
parents | |
children |
line wrap: on
line source
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE LSGSET( PRBTYP, LINOBJ, SINGLR, UNITGZ, UNITQ, $ N, NCLIN, NFREE, $ NROWA, NQ, NROWR, NRANK, NZ, NZ1, $ ISTATE, KX, $ BIGBND, TOLRNK, NUMINF, SUMINF, $ BL, BU, A, RES, FEATOL, $ GQ, CQ, R, X, WTINF, ZY, WRK ) IMPLICIT DOUBLE PRECISION(A-H,O-Z) CHARACTER*2 PRBTYP LOGICAL LINOBJ, SINGLR, UNITGZ, UNITQ INTEGER ISTATE(*), KX(N) DOUBLE PRECISION BL(*), BU(*), A(NROWA,*), $ RES(*), FEATOL(*) DOUBLE PRECISION GQ(N), CQ(*), R(NROWR,*), X(N), WTINF(*), $ ZY(NQ,*) DOUBLE PRECISION WRK(N) ************************************************************************ * LSGSET finds the number and weighted sum of infeasibilities for * the bounds and linear constraints. An appropriate transformed * gradient vector is returned in GQ. * * Positive values of ISTATE(j) will not be altered. These mean * the following... * * 1 2 3 * a'x = bl a'x = bu bl = bu * * Other values of ISTATE(j) will be reset as follows... * a'x lt bl a'x gt bu a'x free * - 2 - 1 0 * * If x is feasible, LSGSET computes the vector Q(free)'g(free), * where g is the gradient of the the sum of squares plus the * linear term. The matrix Q is of the form * ( Q(free) 0 ), * ( 0 I(fixed)) * where Q(free) is the orthogonal factor of A(free) and A is * the matrix of constraints in the working set. The transformed * gradients are stored in GQ. * * Systems Optimization Laboratory, Stanford University. * Original version written 31-October-1984. * Level 2 Blas added 11-June-1986. * This version of LSGSET dated 24-June-1986. ************************************************************************ EXTERNAL DDOT , IDRANK INTRINSIC ABS , MAX , MIN PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 ) BIGUPP = BIGBND BIGLOW = - BIGBND NUMINF = 0 SUMINF = ZERO CALL DLOAD ( N, ZERO, GQ, 1 ) DO 200 J = 1, N+NCLIN IF (ISTATE(J) .LE. 0) THEN FEASJ = FEATOL(J) IF (J .LE. N) THEN CTX = X(J) ELSE K = J - N CTX = DDOT ( N, A(K,1), NROWA, X, 1 ) END IF ISTATE(J) = 0 * See if the lower bound is violated. IF (BL(J) .GT. BIGLOW) THEN S = BL(J) - CTX IF (S .GT. FEASJ ) THEN ISTATE(J) = - 2 WEIGHT = - WTINF(J) GO TO 160 END IF END IF * See if the upper bound is violated. IF (BU(J) .GE. BIGUPP) GO TO 200 S = CTX - BU(J) IF (S .LE. FEASJ ) GO TO 200 ISTATE(J) = - 1 WEIGHT = WTINF(J) * Add the infeasibility. 160 NUMINF = NUMINF + 1 SUMINF = SUMINF + ABS( WEIGHT ) * S IF (J .LE. N) THEN GQ(J) = WEIGHT ELSE CALL DAXPY ( N, WEIGHT, A(K,1), NROWA, GQ, 1 ) END IF END IF 200 CONTINUE * ------------------------------------------------------------------ * Install GQ, the transformed gradient. * ------------------------------------------------------------------ SINGLR = .FALSE. UNITGZ = .TRUE. IF (NUMINF .GT. 0) THEN CALL CMQMUL( 6, N, NZ, NFREE, NQ, UNITQ, KX, GQ, ZY, WRK ) ELSE IF (NUMINF .EQ. 0 .AND. PRBTYP .EQ. 'FP') THEN CALL DLOAD ( N, ZERO, GQ, 1 ) ELSE * Ready for the Optimality Phase. * Set NZ1 so that Rz1 is nonsingular. IF (NRANK .EQ. 0) THEN IF (LINOBJ) THEN CALL DCOPY ( N, CQ, 1, GQ, 1 ) ELSE CALL DLOAD ( N, ZERO, GQ, 1 ) END IF NZ1 = 0 ELSE * Compute GQ = - R' * (transformed residual) CALL DCOPY ( NRANK, RES, 1, GQ, 1 ) CALL DSCAL ( NRANK, (-ONE), GQ, 1 ) CALL DTRMV ( 'U', 'T', 'N', NRANK, R, NROWR, GQ, 1 ) IF (NRANK .LT. N) $ CALL DGEMV( 'T', NRANK, N-NRANK, -ONE,R(1,NRANK+1),NROWR, $ RES, 1, ZERO, GQ(NRANK+1), 1 ) IF (LINOBJ) CALL DAXPY ( N, ONE, CQ, 1, GQ, 1 ) UNITGZ = .FALSE. NZ1 = IDRANK( MIN(NRANK, NZ), R, NROWR+1, TOLRNK ) END IF END IF RETURN * End of LSGSET. END