Mercurial > octave

comparison libcruft/npsol/lsadd.f @ 2329:30c606bec7a8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-07-19 01:29:05 by jwe] Initial revision
author jwe
date Fri, 19 Jul 1996 01:29:55 +0000
parents
children
comparison
equal deleted inserted replaced
2328:b44c3b2a5fce 2329:30c606bec7a8
1 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 * File LSSUBS FORTRAN
3 *
4 * LSADD LSADDS LSBNDS LSCHOL LSCORE LSCRSH LSDEL
5 * LSDFLT LSFEAS LSFILE LSGETP LSGSET LSKEY LSLOC
6 * LSMOVE LSMULS LSOPTN LSPRT LSSETX LSSOL
7 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8
9 SUBROUTINE LSADD ( UNITQ,
10 $ INFORM, IFIX, IADD, JADD,
11 $ NACTIV, NZ, NFREE, NRANK, NRES, NGQ,
12 $ N, NROWA, NQ, NROWR, NROWT,
13 $ KX, CONDMX,
14 $ A, R, T, RES, GQ, ZY, WRK1, WRK2 )
15
16 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
17 LOGICAL UNITQ
18 INTEGER KX(N)
19 DOUBLE PRECISION A(NROWA,*), R(NROWR,*), T(NROWT,*),
20 $ RES(N,*), GQ(N,*), ZY(NQ,*)
21 DOUBLE PRECISION WRK1(N), WRK2(N)
22 ************************************************************************
23 * LSADD updates the factorization, A(free) * (Z Y) = (0 T), when a
24 * constraint is added to the working set. If NRANK .gt. 0, the
25 * factorization ( R ) = PWQ is also updated, where W is the
26 * ( 0 )
27 * least squares matrix, R is upper-triangular, and P is an
28 * orthogonal matrix. The matrices W and P are not stored.
29 *
30 * There are three separate cases to consider (although each case
31 * shares code with another)...
32 *
33 * (1) A free variable becomes fixed on one of its bounds when there
34 * are already some general constraints in the working set.
35 *
36 * (2) A free variable becomes fixed on one of its bounds when there
37 * are only bound constraints in the working set.
38 *
39 * (3) A general constraint (corresponding to row IADD of A) is
40 * added to the working set.
41 *
42 * In cases (1) and (2), we assume that KX(IFIX) = JADD.
43 * In all cases, JADD is the index of the constraint being added.
44 *
45 * If there are no general constraints in the working set, the
46 * matrix Q = (Z Y) is the identity and will not be touched.
47 *
48 * If NRES .GT. 0, the row transformations are applied to the rows of
49 * the (N by NRES) matrix RES.
50 * If NGQ .GT. 0, the column transformations are applied to the
51 * columns of the (NGQ by N) matrix GQ'.
52 *
53 * Systems Optimization Laboratory, Stanford University.
54 * Original version written 31-October--1984.
55 * This version of LSADD dated 29-December-1985.
56 ************************************************************************
57 COMMON /SOL1CM/ NOUT
58 COMMON /SOL4CM/ EPSPT3, EPSPT5, EPSPT8, EPSPT9
59 COMMON /SOL5CM/ ASIZE, DTMAX, DTMIN
60
61 LOGICAL LSDBG
62 PARAMETER (LDBG = 5)
63 COMMON /LSDEBG/ ILSDBG(LDBG), LSDBG
64
65 LOGICAL BOUND , OVERFL
66 EXTERNAL DDOT , DDIV , DNRM2
67 INTRINSIC MAX , MIN
68 PARAMETER (ZERO = 0.0D+0, ONE = 1.0D+0)
69
70 * If the condition estimator of the updated factors is greater than
71 * CONDBD, a warning message is printed.
72
73 CONDBD = ONE / EPSPT9
74
75 OVERFL = .FALSE.
76 BOUND = JADD .LE. N
77 IF (BOUND) THEN
78 * ===============================================================
79 * A simple bound has entered the working set. IADD is not used.
80 * ===============================================================
81 IF (LSDBG .AND. ILSDBG(1) .GT. 0)
82 $ WRITE (NOUT, 1010) NACTIV, NZ, NFREE, IFIX, JADD, UNITQ
83 NANEW = NACTIV
84
85 IF (UNITQ) THEN
86
87 * Q is not stored, but KX defines an ordering of the columns
88 * of the identity matrix that implicitly define Q.
89 * Reorder KX so that variable IFIX is moved to position
90 * NFREE+1 and variables IFIX+1,...,NFREE+1 are moved one
91 * position to the left.
92
93 CALL DLOAD ( NFREE, (ZERO), WRK1, 1 )
94 WRK1(IFIX) = ONE
95
96 DO 100 I = IFIX, NFREE-1
97 KX(I) = KX(I+1)
98 100 CONTINUE
99 ELSE
100 * ------------------------------------------------------------
101 * Q is stored explicitly.
102 * ------------------------------------------------------------
103 * Set WRK1 = the (IFIX)-th row of Q.
104 * Move the (NFREE)-th row of Q to position IFIX.
105
106 CALL DCOPY ( NFREE, ZY(IFIX,1), NQ, WRK1, 1 )
107 IF (IFIX .LT. NFREE) THEN
108 CALL DCOPY ( NFREE, ZY(NFREE,1), NQ, ZY(IFIX,1), NQ )
109 KX(IFIX) = KX(NFREE)
110 END IF
111 END IF
112 KX(NFREE) = JADD
113 ELSE
114 * ===============================================================
115 * A general constraint has entered the working set.
116 * IFIX is not used.
117 * ===============================================================
118 IF (LSDBG .AND. ILSDBG(1) .GT. 0)
119 $ WRITE (NOUT, 1020) NACTIV, NZ, NFREE, IADD, JADD, UNITQ
120
121 NANEW = NACTIV + 1
122
123 * Transform the incoming row of A by Q'.
124
125 CALL DCOPY ( N, A(IADD,1), NROWA, WRK1, 1 )
126 CALL CMQMUL( 8, N, NZ, NFREE, NQ, UNITQ, KX, WRK1, ZY, WRK2)
127
128 * Check that the incoming row is not dependent upon those
129 * already in the working set.
130
131 DTNEW = DNRM2 ( NZ, WRK1, 1 )
132 IF (NACTIV .EQ. 0) THEN
133
134 * This is the only general constraint in the working set.
135
136 COND = DDIV ( ASIZE, DTNEW, OVERFL )
137 TDTMAX = DTNEW
138 TDTMIN = DTNEW
139 ELSE
140
141 * There are already some general constraints in the working
142 * set. Update the estimate of the condition number.
143
144 TDTMAX = MAX( DTNEW, DTMAX )
145 TDTMIN = MIN( DTNEW, DTMIN )
146 COND = DDIV ( TDTMAX, TDTMIN, OVERFL )
147 END IF
148
149 IF (COND .GT. CONDMX .OR. OVERFL) GO TO 900
150
151 IF (UNITQ) THEN
152
153 * This is the first general constraint to be added.
154 * Set Q = I.
155
156 DO 200 J = 1, NFREE
157 CALL DLOAD ( NFREE, (ZERO), ZY(1,J), 1 )
158 ZY(J,J) = ONE
159 200 CONTINUE
160 UNITQ = .FALSE.
161 END IF
162 END IF
163
164 NZERO = NZ - 1
165 IF (BOUND) NZERO = NFREE - 1
166
167 * ------------------------------------------------------------------
168 * Use a sequence of 2*2 column transformations to reduce the
169 * first NZERO elements of WRK1 to zero. This affects ZY, except
170 * when UNITQ is true. The transformations may also be applied
171 * to R, T and GQ'.
172 * ------------------------------------------------------------------
173 LROWR = N
174 NELM = 1
175 IROWT = NACTIV
176
177 DO 300 K = 1, NZERO
178
179 * Compute the transformation that reduces WRK1(K) to zero,
180 * then apply it to the relevant columns of Z and GQ'.
181
182 CALL DROT3G( WRK1(K+1), WRK1(K), CS, SN )
183 IF (.NOT. UNITQ)
184 $ CALL DROT3 ( NFREE, ZY(1,K+1), 1, ZY(1,K), 1, CS, SN )
185 IF (NGQ .GT. 0)
186 $ CALL DROT3 ( NGQ , GQ(K+1,1), N, GQ(K,1), N, CS, SN )
187
188 IF (K .GE. NZ .AND. NACTIV .GT. 0) THEN
189
190 * Apply the rotation to T.
191
192 T(IROWT,K) = ZERO
193 CALL DROT3 ( NELM, T(IROWT,K+1), 1, T(IROWT,K), 1, CS, SN )
194 NELM = NELM + 1
195 IROWT = IROWT - 1
196 END IF
197
198 IF (NRANK .GT. 0) THEN
199
200 * Apply the same transformation to the columns of R.
201 * This generates a subdiagonal element in R that must be
202 * eliminated by a row rotation.
203
204 IF (K .LT. NRANK) R(K+1,K) = ZERO
205 LCOL = MIN( K+1, NRANK )
206
207 CALL DROT3 ( LCOL, R(1,K+1), 1, R(1,K), 1, CS, SN )
208 IF (K .LT. NRANK) THEN
209 CALL DROT3G( R(K,K), R(K+1,K), CS, SN )
210 LROWR = LROWR - 1
211 CALL DROT3 ( LROWR, R(K,K+1) , NROWR,
212 $ R(K+1,K+1), NROWR, CS, SN )
213
214 IF (NRES .GT. 0)
215 $ CALL DROT3 ( NRES, RES(K,1) , N ,
216 $ RES(K+1,1), N , CS, SN )
217 END IF
218 END IF
219 300 CONTINUE
220
221 IF (BOUND) THEN
222
223 * The last row and column of ZY has been transformed to plus
224 * or minus the unit vector E(NFREE). We can reconstitute the
225 * columns of GQ and R corresponding to the new fixed variable.
226
227 IF (WRK1(NFREE) .LT. ZERO) THEN
228 NFMIN = MIN( NRANK, NFREE )
229 IF (NFMIN .GT. 0) CALL DSCAL ( NFMIN, -ONE, R(1,NFREE) , 1 )
230 IF (NGQ .GT. 0) CALL DSCAL ( NGQ , -ONE, GQ(NFREE,1), N )
231 END IF
232
233 * ---------------------------------------------------------------
234 * The diagonals of T have been altered. Recompute the
235 * largest and smallest values.
236 * ---------------------------------------------------------------
237 IF (NACTIV .GT. 0) THEN
238 CALL DCOND( NACTIV, T(NACTIV,NZ), NROWT-1, TDTMAX, TDTMIN )
239 COND = DDIV ( TDTMAX, TDTMIN, OVERFL )
240 END IF
241 ELSE
242 * ---------------------------------------------------------------
243 * General constraint. Install the new row of T.
244 * ---------------------------------------------------------------
245 CALL DCOPY ( NANEW, WRK1(NZ), 1, T(NANEW,NZ), NROWT )
246 END IF
247
248 * ==================================================================
249 * Prepare to exit. Check the magnitude of the condition estimator.
250 * ==================================================================
251 900 IF (NANEW .GT. 0) THEN
252 IF (COND .LT. CONDMX .AND. .NOT. OVERFL) THEN
253
254 * The factorization has been successfully updated.
255
256 INFORM = 0
257 DTMAX = TDTMAX
258 DTMIN = TDTMIN
259 IF (COND .GE. CONDBD) WRITE (NOUT, 2000) JADD
260 ELSE
261
262 * The proposed working set appears to be linearly dependent.
263
264 INFORM = 1
265 IF (LSDBG .AND. ILSDBG(1) .GT. 0) THEN
266 WRITE( NOUT, 3000 )
267 IF (BOUND) THEN
268 WRITE (NOUT, 3010) ASIZE, DTMAX, DTMIN
269 ELSE
270 IF (NACTIV .GT. 0) THEN
271 WRITE (NOUT, 3020) ASIZE, DTMAX, DTMIN, DTNEW
272 ELSE
273 WRITE (NOUT, 3030) ASIZE, DTNEW
274 END IF
275 END IF
276 END IF
277 END IF
278 END IF
279
280 RETURN
281
282 1010 FORMAT(/ ' //LSADD // Simple bound added.'
283 $ / ' //LSADD // NACTIV NZ NFREE IFIX JADD UNITQ'
284 $ / ' //LSADD // ', 5I6, L6 )
285 1020 FORMAT(/ ' //LSADD // General constraint added. '
286 $ / ' //LSADD // NACTIV NZ NFREE IADD JADD UNITQ'
287 $ / ' //LSADD // ', 5I6, L6 )
288 2000 FORMAT(/ ' XXX Serious ill-conditioning in the working set',
289 $ ' after adding constraint ', I5
290 $ / ' XXX Overflow may occur in subsequent iterations.'//)
291 3000 FORMAT(/ ' //LSADD // Dependent constraint rejected.' )
292 3010 FORMAT(/ ' //LSADD // ASIZE DTMAX DTMIN '
293 $ / ' //LSADD //', 1P3E10.2 )
294 3020 FORMAT(/ ' //LSADD // ASIZE DTMAX DTMIN DTNEW'
295 $ / ' //LSADD //', 1P4E10.2 )
296 3030 FORMAT(/ ' //LSADD // ASIZE DTNEW'
297 $ / ' //LSADD //', 1P2E10.2 )
298
299 * End of LSADD .
300
301 END