Line data Source code
1 : !!****m* ABINIT/m_slsqp
2 : !! NAME
3 : !! m_slsqp
4 : !!
5 : !! FUNCTION
6 : !! This module contains the SLSQP optimization routine.
7 : !! This was taken from scipy.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (c) 2001-2002 Enthought, Inc. 2003-2024, SciPy Developers.
11 : !! All rights reserved.
12 : !!
13 : !! Redistribution and use in source and binary forms, with or without
14 : !! modification, are permitted provided that the following conditions
15 : !! are met:
16 : !!
17 : !! 1. Redistributions of source code must retain the above copyright
18 : !! notice, this list of conditions and the following disclaimer.
19 : !!
20 : !! 2. Redistributions in binary form must reproduce the above
21 : !! copyright notice, this list of conditions and the following
22 : !! disclaimer in the documentation and/or other materials provided
23 : !! with the distribution.
24 : !!
25 : !! 3. Neither the name of the copyright holder nor the names of its
26 : !! contributors may be used to endorse or promote products derived
27 : !! from this software without specific prior written permission.
28 : !!
29 : !! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
30 : !! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
31 : !! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
32 : !! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
33 : !! OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
34 : !! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
35 : !! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
36 : !! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
37 : !! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
38 : !! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
39 : !! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
40 : !!
41 : !! This file is distributed under the terms of the
42 : !! GNU General Public License, see ~abinit/COPYING
43 : !! or http://www.gnu.org/copyleft/gpl.txt .
44 : !!
45 : !! SOURCE
46 :
47 : #if defined HAVE_CONFIG_H
48 : #include "config.h"
49 : #endif
50 :
51 : #include "abi_common.h"
52 :
53 : MODULE m_slsqp
54 :
55 : implicit none
56 :
57 : private
58 :
59 : public :: slsqp
60 :
61 : CONTAINS !========================================================================================
62 : !!***
63 :
64 : !C
65 : !C ALGORITHM 733, COLLECTED ALGORITHMS FROM ACM.
66 : !C TRANSACTIONS ON MATHEMATICAL SOFTWARE,
67 : !C VOL. 20, NO. 3, SEPTEMBER, 1994, PP. 262-281.
68 : !C https://doi.org/10.1145/192115.192124
69 : !C
70 : !C
71 : !C https://web.archive.org/web/20170106155705/http://permalink.gmane.org/gmane.comp.python.scientific.devel/6725
72 : !C ------
73 : !C From: Deborah Cotton <cotton@hq.acm.org>
74 : !C Date: Fri, 14 Sep 2007 12:35:55 -0500
75 : !C Subject: RE: Algorithm License requested
76 : !C To: Alan Isaac
77 : !C
78 : !C Prof. Issac,
79 : !C
80 : !C In that case, then because the author consents to [the ACM] releasing
81 : !C the code currently archived at http://www.netlib.org/toms/733 under the
82 : !C BSD license, the ACM hereby releases this code under the BSD license.
83 : !C
84 : !C Regards,
85 : !C
86 : !C Deborah Cotton, Copyright & Permissions
87 : !C ACM Publications
88 : !C 2 Penn Plaza, Suite 701**
89 : !C New York, NY 10121-0701
90 : !C permissions@acm.org
91 : !C 212.869.7440 ext. 652
92 : !C Fax. 212.869.0481
93 : !C ------
94 : !C
95 :
96 : !************************************************************************
97 : !* optimizer *
98 : !************************************************************************
99 :
100 0 : SUBROUTINE slsqp (m, meq, la, n, x, xl, xu, f, c, g, a, &
101 0 : acc, iter, mode, w, l_w, jw, l_jw, &
102 : alpha, f0, gs, h1, h2, h3, h4, t, t0, tol, &
103 : iexact, incons, ireset, itermx, line, &
104 : n1, n2, n3)
105 :
106 : !C SLSQP S EQUENTIAL L EAST SQ UARES P ROGRAMMING
107 : !C TO SOLVE GENERAL NONLINEAR OPTIMIZATION PROBLEMS
108 :
109 : !C***********************************************************************
110 : !C* *
111 : !C* *
112 : !C* A NONLINEAR PROGRAMMING METHOD WITH *
113 : !C* QUADRATIC PROGRAMMING SUBPROBLEMS *
114 : !C* *
115 : !C* *
116 : !C* THIS SUBROUTINE SOLVES THE GENERAL NONLINEAR PROGRAMMING PROBLEM *
117 : !C* *
118 : !C* MINIMIZE F(X) *
119 : !C* *
120 : !C* SUBJECT TO C (X) .EQ. 0 , J = 1,...,MEQ *
121 : !C* J *
122 : !C* *
123 : !C* C (X) .GE. 0 , J = MEQ+1,...,M *
124 : !C* J *
125 : !C* *
126 : !C* XL .LE. X .LE. XU , I = 1,...,N. *
127 : !C* I I I *
128 : !C* *
129 : !C* THE ALGORITHM IMPLEMENTS THE METHOD OF HAN AND POWELL *
130 : !C* WITH BFGS-UPDATE OF THE B-MATRIX AND L1-TEST FUNCTION *
131 : !C* WITHIN THE STEPLENGTH ALGORITHM. *
132 : !C* *
133 : !C* PARAMETER DESCRIPTION: *
134 : !C* ( * MEANS THIS PARAMETER WILL BE CHANGED DURING CALCULATION ) *
135 : !C* *
136 : !C* M IS THE TOTAL NUMBER OF CONSTRAINTS, M .GE. 0 *
137 : !C* MEQ IS THE NUMBER OF EQUALITY CONSTRAINTS, MEQ .GE. 0 *
138 : !C* LA SEE A, LA .GE. MAX(M,1) *
139 : !C* N IS THE NUMBER OF VARIBLES, N .GE. 1 *
140 : !C* * X() X() STORES THE CURRENT ITERATE OF THE N VECTOR X *
141 : !C* ON ENTRY X() MUST BE INITIALIZED. ON EXIT X() *
142 : !C* STORES THE SOLUTION VECTOR X IF MODE = 0. *
143 : !C* XL() XL() STORES AN N VECTOR OF LOWER BOUNDS XL TO X. *
144 : !C* ELEMENTS MAY BE NAN TO INDICATE NO LOWER BOUND. *
145 : !C* XU() XU() STORES AN N VECTOR OF UPPER BOUNDS XU TO X. *
146 : !C* ELEMENTS MAY BE NAN TO INDICATE NO UPPER BOUND. *
147 : !C* F IS THE VALUE OF THE OBJECTIVE FUNCTION. *
148 : !C* C() C() STORES THE M VECTOR C OF CONSTRAINTS, *
149 : !C* EQUALITY CONSTRAINTS (IF ANY) FIRST. *
150 : !C* DIMENSION OF C MUST BE GREATER OR EQUAL LA, *
151 : !C* which must be GREATER OR EQUAL MAX(1,M). *
152 : !C* G() G() STORES THE N VECTOR G OF PARTIALS OF THE *
153 : !C* OBJECTIVE FUNCTION; DIMENSION OF G MUST BE *
154 : !C* GREATER OR EQUAL N+1. *
155 : !C* A(),LA,M,N THE LA BY N + 1 ARRAY A() STORES *
156 : !C* THE M BY N MATRIX A OF CONSTRAINT NORMALS. *
157 : !C* A() HAS FIRST DIMENSIONING PARAMETER LA, *
158 : !C* WHICH MUST BE GREATER OR EQUAL MAX(1,M). *
159 : !C* F,C,G,A MUST ALL BE SET BY THE USER BEFORE EACH CALL. *
160 : !C* * ACC ABS(ACC) CONTROLS THE FINAL ACCURACY. *
161 : !C* IF ACC .LT. ZERO AN EXACT LINESEARCH IS PERFORMED,*
162 : !C* OTHERWISE AN ARMIJO-TYPE LINESEARCH IS USED. *
163 : !C* * ITER PRESCRIBES THE MAXIMUM NUMBER OF ITERATIONS. *
164 : !C* ON EXIT ITER INDICATES THE NUMBER OF ITERATIONS. *
165 : !C* * MODE MODE CONTROLS CALCULATION: *
166 : !C* REVERSE COMMUNICATION IS USED IN THE SENSE THAT *
167 : !C* THE PROGRAM IS INITIALIZED BY MODE = 0; THEN IT IS*
168 : !C* TO BE CALLED REPEATEDLY BY THE USER UNTIL A RETURN*
169 : !C* WITH MODE .NE. IABS(1) TAKES PLACE. *
170 : !C* IF MODE = -1 GRADIENTS HAVE TO BE CALCULATED, *
171 : !C* WHILE WITH MODE = 1 FUNCTIONS HAVE TO BE CALCULATED
172 : !C* MODE MUST NOT BE CHANGED BETWEEN SUBSEQUENT CALLS *
173 : !C* OF SQP. *
174 : !C* EVALUATION MODES: *
175 : !C* MODE = -1: GRADIENT EVALUATION, (G&A) *
176 : !C* 0: ON ENTRY: INITIALIZATION, (F,G,C&A) *
177 : !C* ON EXIT : REQUIRED ACCURACY FOR SOLUTION OBTAINED *
178 : !C* 1: FUNCTION EVALUATION, (F&C) *
179 : !C* *
180 : !C* FAILURE MODES: *
181 : !C* 2: NUMBER OF EQUALITY CONSTRAINTS LARGER THAN N *
182 : !C* 3: MORE THAN 3*N ITERATIONS IN LSQ SUBPROBLEM *
183 : !C* 4: INEQUALITY CONSTRAINTS INCOMPATIBLE *
184 : !C* 5: SINGULAR MATRIX E IN LSQ SUBPROBLEM *
185 : !C* 6: SINGULAR MATRIX C IN LSQ SUBPROBLEM *
186 : !C* 7: RANK-DEFICIENT EQUALITY CONSTRAINT SUBPROBLEM HFTI*
187 : !C* 8: POSITIVE DIRECTIONAL DERIVATIVE FOR LINESEARCH *
188 : !C* 9: MORE THAN ITER ITERATIONS IN SQP *
189 : !C* >=10: WORKING SPACE W OR JW TOO SMALL, *
190 : !C* W SHOULD BE ENLARGED TO L_W=MODE/1000 *
191 : !C* JW SHOULD BE ENLARGED TO L_JW=MODE-1000*L_W *
192 : !C* * W(), L_W W() IS A ONE DIMENSIONAL WORKING SPACE, *
193 : !C* THE LENGTH L_W OF WHICH SHOULD BE AT LEAST *
194 : !C* (3*N1+M)*(N1+1) for LSQ *
195 : !C* +(N1-MEQ+1)*(MINEQ+2) + 2*MINEQ for LSI *
196 : !C* +(N1+MINEQ)*(N1-MEQ) + 2*MEQ + N1 for LSEI *
197 : !C* + N1*N/2 + 2*M + 3*N + 3*N1 + 1 for SLSQPB *
198 : !C* with MINEQ = M - MEQ + 2*N1 & N1 = N+1 *
199 : !C* NOTICE: FOR PROPER DIMENSIONING OF W IT IS RECOMMENDED TO *
200 : !C* COPY THE FOLLOWING STATEMENTS INTO THE HEAD OF *
201 : !C* THE CALLING PROGRAM (AND REMOVE THE COMMENT C) *
202 : !c#######################################################################
203 : !C INTEGER LEN_W, LEN_JW, M, N, N1, MEQ, MINEQ
204 : !C PARAMETER (M=... , MEQ=... , N=... )
205 : !C PARAMETER (N1= N+1, MINEQ= M-MEQ+N1+N1)
206 : !C PARAMETER (LEN_W=
207 : !c $ (3*N1+M)*(N1+1)
208 : !c $ +(N1-MEQ+1)*(MINEQ+2) + 2*MINEQ
209 : !c $ +(N1+MINEQ)*(N1-MEQ) + 2*MEQ + N1
210 : !c $ +(N+1)*N/2 + 2*M + 3*N + 3*N1 + 1,
211 : !c $ LEN_JW=MINEQ)
212 : !C DOUBLE PRECISION W(LEN_W)
213 : !C INTEGER JW(LEN_JW)
214 : !c#######################################################################
215 : !C* THE FIRST M+N+N*N1/2 ELEMENTS OF W MUST NOT BE *
216 : !C* CHANGED BETWEEN SUBSEQUENT CALLS OF SLSQP. *
217 : !C* ON RETURN W(1) ... W(M) CONTAIN THE MULTIPLIERS *
218 : !C* ASSOCIATED WITH THE GENERAL CONSTRAINTS, WHILE *
219 : !C* W(M+1) ... W(M+N(N+1)/2) STORE THE CHOLESKY FACTOR*
220 : !C* L*D*L(T) OF THE APPROXIMATE HESSIAN OF THE *
221 : !C* LAGRANGIAN COLUMNWISE DENSE AS LOWER TRIANGULAR *
222 : !C* UNIT MATRIX L WITH D IN ITS 'DIAGONAL' and *
223 : !C* W(M+N(N+1)/2+N+2 ... W(M+N(N+1)/2+N+2+M+2N) *
224 : !C* CONTAIN THE MULTIPLIERS ASSOCIATED WITH ALL *
225 : !C* ALL CONSTRAINTS OF THE QUADRATIC PROGRAM FINDING *
226 : !C* THE SEARCH DIRECTION TO THE SOLUTION X* *
227 : !C* * JW(), L_JW JW() IS A ONE DIMENSIONAL INTEGER WORKING SPACE *
228 : !C* THE LENGTH L_JW OF WHICH SHOULD BE AT LEAST *
229 : !C* MINEQ *
230 : !C* with MINEQ = M - MEQ + 2*N1 & N1 = N+1 *
231 : !C* *
232 : !C* THE USER HAS TO PROVIDE THE FOLLOWING SUBROUTINES: *
233 : !C* LDL(N,A,Z,SIG,W) : UPDATE OF THE LDL'-FACTORIZATION. *
234 : !C* LINMIN(A,B,F,TOL) : LINESEARCH ALGORITHM IF EXACT = 1 *
235 : !C* LSQ(M,MEQ,LA,N,NC,C,D,A,B,XL,XU,X,LAMBDA,W,....) : *
236 : !C* *
237 : !C* SOLUTION OF THE QUADRATIC PROGRAM *
238 : !C* QPSOL IS RECOMMENDED: *
239 : !C* PE GILL, W MURRAY, MA SAUNDERS, MH WRIGHT: *
240 : !C* USER'S GUIDE FOR SOL/QPSOL: *
241 : !C* A FORTRAN PACKAGE FOR QUADRATIC PROGRAMMING, *
242 : !C* TECHNICAL REPORT SOL 83-7, JULY 1983 *
243 : !C* DEPARTMENT OF OPERATIONS RESEARCH, STANFORD UNIVERSITY *
244 : !C* STANFORD, CA 94305 *
245 : !C* QPSOL IS THE MOST ROBUST AND EFFICIENT QP-SOLVER *
246 : !C* AS IT ALLOWS WARM STARTS WITH PROPER WORKING SETS *
247 : !C* *
248 : !C* IF IT IS NOT AVAILABLE USE LSEI, A CONSTRAINT LINEAR LEAST *
249 : !C* SQUARES SOLVER IMPLEMENTED USING THE SOFTWARE HFTI, LDP, NNLS *
250 : !C* FROM C.L. LAWSON, R.J.HANSON: SOLVING LEAST SQUARES PROBLEMS, *
251 : !C* PRENTICE HALL, ENGLEWOOD CLIFFS, 1974. *
252 : !C* LSEI COMES WITH THIS PACKAGE, together with all necessary SR's. *
253 : !C* *
254 : !C* TOGETHER WITH A COUPLE OF SUBROUTINES FROM BLAS LEVEL 1 *
255 : !C* *
256 : !C* SQP IS HEAD SUBROUTINE FOR BODY SUBROUTINE SQPBDY *
257 : !C* IN WHICH THE ALGORITHM HAS BEEN IMPLEMENTED. *
258 : !C* *
259 : !C* IMPLEMENTED BY: DIETER KRAFT, DFVLR OBERPFAFFENHOFEN *
260 : !C* as described in Dieter Kraft: A Software Package for *
261 : !C* Sequential Quadratic Programming *
262 : !C* DFVLR-FB 88-28, 1988 *
263 : !C* which should be referenced if the user publishes results of SLSQP *
264 : !C* *
265 : !C* DATE: APRIL - OCTOBER, 1981. *
266 : !C* STATUS: DECEMBER, 31-ST, 1984. *
267 : !C* STATUS: MARCH , 21-ST, 1987, REVISED TO FORTRAN 77 *
268 : !C* STATUS: MARCH , 20-th, 1989, REVISED TO MS-FORTRAN *
269 : !C* STATUS: APRIL , 14-th, 1989, HESSE in-line coded *
270 : !C* STATUS: FEBRUARY, 28-th, 1991, FORTRAN/2 Version 1.04 *
271 : !C* accepts Statement Functions *
272 : !C* STATUS: MARCH , 1-st, 1991, tested with SALFORD *
273 : !C* FTN77/386 COMPILER VERS 2.40*
274 : !C* in protected mode *
275 : !C* *
276 : !C***********************************************************************
277 : !C* *
278 : !C* Copyright 1991: Dieter Kraft, FHM *
279 : !C* *
280 : !C***********************************************************************
281 :
282 : INTEGER il, im, ir, is, iter, iu, iv, iw, ix, l_w, l_jw, &
283 : jw(l_jw), la, m, meq, mineq, mode, n
284 :
285 : DOUBLE PRECISION acc, a(la,n+1), c(la), f, g(n+1), &
286 : x(n), xl(n), xu(n), w(l_w)
287 :
288 : INTEGER iexact, incons, ireset, itermx, line, n1, n2, n3
289 :
290 : DOUBLE PRECISION alpha, f0, gs, h1, h2, h3, h4, t, t0, tol
291 :
292 : !c dim(W) = N1*(N1+1) + MEQ*(N1+1) + MINEQ*(N1+1) for LSQ
293 : !c +(N1-MEQ+1)*(MINEQ+2) + 2*MINEQ for LSI
294 : !c +(N1+MINEQ)*(N1-MEQ) + 2*MEQ + N1 for LSEI
295 : !c + N1*N/2 + 2*M + 3*N +3*N1 + 1 for SLSQPB
296 : !c with MINEQ = M - MEQ + 2*N1 & N1 = N+1
297 :
298 : !C CHECK LENGTH OF WORKING ARRAYS
299 :
300 0 : n1 = n+1
301 0 : mineq = m-meq+n1+n1
302 : il = (3*n1+m)*(n1+1) + &
303 : (n1-meq+1)*(mineq+2) + 2*mineq + &
304 : (n1+mineq)*(n1-meq) + 2*meq + &
305 0 : n1*n/2 + 2*m + 3*n + 4*n1 + 1
306 0 : im = MAX(mineq, n1-meq)
307 0 : IF (l_w .LT. il .OR. l_jw .LT. im) THEN
308 0 : mode = 1000*MAX(10,il)
309 0 : mode = mode+MAX(10,im)
310 0 : RETURN
311 : ENDIF
312 :
313 : !C PREPARE DATA FOR CALLING SQPBDY - INITIAL ADDRESSES IN W
314 :
315 0 : im = 1
316 0 : il = im + MAX(1,m)
317 0 : il = im + la
318 0 : ix = il + n1*n/2 + 1
319 0 : ir = ix + n
320 0 : is = ir + n + n + MAX(1,m)
321 0 : is = ir + n + n + la
322 0 : iu = is + n1
323 0 : iv = iu + n1
324 0 : iw = iv + n1
325 :
326 : CALL slsqpb (m, meq, la, n, x, xl, xu, f, c, g, a, acc, iter, &
327 : mode, w(ir), w(il), w(ix), w(im), w(is), w(iu), w(iv), w(iw), jw, &
328 : alpha, f0, gs, h1, h2, h3, h4, t, t0, tol, &
329 : iexact, incons, ireset, itermx, line, &
330 0 : n1, n2, n3)
331 :
332 : END SUBROUTINE slsqp
333 :
334 0 : SUBROUTINE slsqpb (m, meq, la, n, x, xl, xu, f, c, g, a, acc, &
335 0 : iter, mode, r, l, x0, mu, s, u, v, w, iw, &
336 : alpha, f0, gs, h1, h2, h3, h4, t, t0, tol, &
337 : iexact, incons, ireset, itermx, line, &
338 : n1, n2, n3)
339 :
340 : !C NONLINEAR PROGRAMMING BY SOLVING SEQUENTIALLY QUADRATIC PROGRAMS
341 :
342 : !C - L1 - LINE SEARCH, POSITIVE DEFINITE BFGS UPDATE -
343 :
344 : !C BODY SUBROUTINE FOR SLSQP
345 :
346 : INTEGER iw(*), i, iexact, incons, ireset, iter, itermx, &
347 : k, j, la, line, m, meq, mode, n, n1, n2, n3
348 : LOGICAL badlin
349 :
350 : DOUBLE PRECISION a(la,n+1), c(la), g(n+1), l((n+1)*(n+2)/2), &
351 : mu(la), r(m+n+n+2), s(n+1), u(n+1), v(n+1), w(*), &
352 : x(n), xl(n), xu(n), x0(n), &
353 : acc, alfmin, alpha, f, f0, gs, h1, h2, h3, h4, &
354 : hun, one, t, t0, ten, tol, two, ZERO
355 :
356 : !c dim(W) = N1*(N1+1) + MEQ*(N1+1) + MINEQ*(N1+1) for LSQ
357 : !c +(N1-MEQ+1)*(MINEQ+2) + 2*MINEQ
358 : !c +(N1+MINEQ)*(N1-MEQ) + 2*MEQ + N1 for LSEI
359 : !c with MINEQ = M - MEQ + 2*N1 & N1 = N+1
360 :
361 : DATA ZERO /0.0d0/, one /1.0d0/, alfmin /1.0d-1/, &
362 : hun /1.0d+2/, ten /1.0d+1/, two /2.0d0/
363 :
364 : !C The badlin flag keeps track whether the SQP problem on the current
365 : !C iteration was inconsistent or not.
366 0 : badlin = .false.
367 :
368 0 : IF (mode.LT.0) GO TO 260
369 0 : IF (mode.EQ.0) GO TO 100
370 : IF (mode.GT.0) GO TO 220
371 :
372 0 : 100 itermx = iter
373 0 : IF (acc.GE.ZERO) THEN
374 0 : iexact = 0
375 : ELSE
376 0 : iexact = 1
377 : ENDIF
378 0 : acc = ABS(acc)
379 0 : tol = ten*acc
380 0 : iter = 0
381 0 : ireset = 0
382 0 : n1 = n + 1
383 0 : n2 = n1*n/2
384 0 : n3 = n2 + 1
385 0 : s(1) = ZERO
386 0 : mu(1) = ZERO
387 0 : CALL dcopy_(n, s(1), 0, s, 1)
388 0 : CALL dcopy_(m, mu(1), 0, mu, 1)
389 :
390 : !C RESET BFGS MATRIX
391 :
392 0 : 110 ireset = ireset + 1
393 0 : IF (ireset.GT.5) GO TO 255
394 0 : l(1) = ZERO
395 0 : CALL dcopy_(n2, l(1), 0, l, 1)
396 0 : j = 1
397 0 : DO 120 i=1,n
398 0 : l(j) = one
399 0 : j = j + n1 - i
400 0 : 120 CONTINUE
401 :
402 : !C MAIN ITERATION : SEARCH DIRECTION, STEPLENGTH, LDL'-UPDATE
403 :
404 0 : 130 iter = iter + 1
405 0 : mode = 9
406 0 : IF (iter.GT.itermx) GO TO 330
407 :
408 : !C SEARCH DIRECTION AS SOLUTION OF QP - SUBPROBLEM
409 :
410 0 : CALL dcopy_(n, xl, 1, u, 1)
411 0 : CALL dcopy_(n, xu, 1, v, 1)
412 0 : CALL daxpy_sl(n, -one, x, 1, u, 1)
413 0 : CALL daxpy_sl(n, -one, x, 1, v, 1)
414 0 : h4 = one
415 0 : CALL lsq (m, meq, n , n3, la, l, g, a, c, u, v, s, r, w, iw, mode)
416 :
417 : !C AUGMENTED PROBLEM FOR INCONSISTENT LINEARIZATION
418 : !C
419 : !C If it turns out that the original SQP problem is inconsistent,
420 : !C disallow termination with convergence on this iteration,
421 : !C even if the augmented problem was solved.
422 :
423 0 : badlin = .false.
424 0 : IF (mode.EQ.6) THEN
425 0 : IF (n.EQ.meq) THEN
426 0 : mode = 4
427 : ENDIF
428 : ENDIF
429 0 : IF (mode.EQ.4) THEN
430 0 : badlin = .true.
431 0 : DO 140 j=1,m
432 0 : IF (j.LE.meq) THEN
433 0 : a(j,n1) = -c(j)
434 : ELSE
435 0 : a(j,n1) = MAX(-c(j),ZERO)
436 : ENDIF
437 0 : 140 CONTINUE
438 0 : s(1) = ZERO
439 0 : CALL dcopy_(n, s(1), 0, s, 1)
440 0 : h3 = ZERO
441 0 : g(n1) = ZERO
442 0 : l(n3) = hun
443 0 : s(n1) = one
444 0 : u(n1) = ZERO
445 0 : v(n1) = one
446 0 : incons = 0
447 : 150 CALL lsq (m, meq, n1, n3, la, l, g, a, c, u, v, s, r, &
448 0 : w, iw, mode)
449 0 : h4 = one - s(n1)
450 0 : IF (mode.EQ.4) THEN
451 0 : l(n3) = ten*l(n3)
452 0 : incons = incons + 1
453 0 : IF (incons.GT.5) GO TO 330
454 : GOTO 150
455 0 : ELSE IF (mode.NE.1) THEN
456 : GOTO 330
457 : ENDIF
458 0 : ELSE IF (mode.NE.1) THEN
459 : GOTO 330
460 : ENDIF
461 :
462 : !C UPDATE MULTIPLIERS FOR L1-TEST
463 :
464 0 : DO 160 i=1,n
465 0 : v(i) = g(i) - ddot_sl(m,a(1,i),1,r,1)
466 0 : 160 CONTINUE
467 0 : f0 = f
468 0 : CALL dcopy_(n, x, 1, x0, 1)
469 0 : gs = ddot_sl(n, g, 1, s, 1)
470 0 : h1 = ABS(gs)
471 0 : h2 = ZERO
472 0 : DO 170 j=1,m
473 0 : IF (j.LE.meq) THEN
474 0 : h3 = c(j)
475 : ELSE
476 0 : h3 = ZERO
477 : ENDIF
478 0 : h2 = h2 + MAX(-c(j),h3)
479 0 : h3 = ABS(r(j))
480 0 : mu(j) = MAX(h3,(mu(j)+h3)/two)
481 0 : h1 = h1 + h3*ABS(c(j))
482 0 : 170 CONTINUE
483 :
484 : !C CHECK CONVERGENCE
485 :
486 0 : mode = 0
487 : IF (h1.LT.acc .AND. h2.LT.acc .AND. .NOT. badlin &
488 0 : .AND. f .EQ. f) GO TO 330
489 0 : h1 = ZERO
490 0 : DO 180 j=1,m
491 0 : IF (j.LE.meq) THEN
492 0 : h3 = c(j)
493 : ELSE
494 0 : h3 = ZERO
495 : ENDIF
496 0 : h1 = h1 + mu(j)*MAX(-c(j),h3)
497 0 : 180 CONTINUE
498 0 : t0 = f + h1
499 0 : h3 = gs - h1*h4
500 0 : mode = 8
501 0 : IF (h3.GE.ZERO) GO TO 110
502 :
503 : !C LINE SEARCH WITH AN L1-TESTFUNCTION
504 :
505 0 : line = 0
506 0 : alpha = one
507 0 : IF (iexact.EQ.1) GOTO 210
508 :
509 : !C INEXACT LINESEARCH
510 :
511 0 : 190 line = line + 1
512 0 : h3 = alpha*h3
513 0 : CALL dscal_sl(n, alpha, s, 1)
514 0 : CALL dcopy_(n, x0, 1, x, 1)
515 0 : CALL daxpy_sl(n, one, s, 1, x, 1)
516 0 : mode = 1
517 0 : GO TO 330
518 0 : 200 IF (h1.LE.h3/ten .OR. line.GT.10) GO TO 240
519 0 : alpha = MAX(h3/(two*(h3-h1)),alfmin)
520 0 : GO TO 190
521 :
522 : !C EXACT LINESEARCH
523 :
524 0 : 210 IF (line.NE.3) THEN
525 0 : alpha = linmin(line,alfmin,one,t,tol)
526 0 : CALL dcopy_(n, x0, 1, x, 1)
527 0 : CALL daxpy_sl(n, alpha, s, 1, x, 1)
528 0 : mode = 1
529 0 : GOTO 330
530 : ENDIF
531 0 : CALL dscal_sl(n, alpha, s, 1)
532 0 : GOTO 240
533 :
534 : !C CALL FUNCTIONS AT CURRENT X
535 :
536 0 : 220 t = f
537 0 : DO 230 j=1,m
538 0 : IF (j.LE.meq) THEN
539 0 : h1 = c(j)
540 : ELSE
541 0 : h1 = ZERO
542 : ENDIF
543 0 : t = t + mu(j)*MAX(-c(j),h1)
544 0 : 230 CONTINUE
545 0 : h1 = t - t0
546 0 : GOTO (200, 210) iexact+1
547 :
548 : !C CHECK CONVERGENCE
549 :
550 0 : 240 h3 = ZERO
551 0 : DO 250 j=1,m
552 0 : IF (j.LE.meq) THEN
553 0 : h1 = c(j)
554 : ELSE
555 0 : h1 = ZERO
556 : ENDIF
557 0 : h3 = h3 + MAX(-c(j),h1)
558 0 : 250 CONTINUE
559 : IF ((ABS(f-f0).LT.acc .OR. dnrm2_(n,s,1).LT.acc) .AND. h3.LT.acc &
560 0 : .AND. .NOT. badlin .AND. f .EQ. f) &
561 : THEN
562 0 : mode = 0
563 : ELSE
564 0 : mode = -1
565 : ENDIF
566 0 : GO TO 330
567 :
568 : !C CHECK relaxed CONVERGENCE in case of positive directional derivative
569 :
570 : 255 CONTINUE
571 0 : h3 = ZERO
572 0 : DO 256 j=1,m
573 0 : IF (j.LE.meq) THEN
574 0 : h1 = c(j)
575 : ELSE
576 0 : h1 = ZERO
577 : ENDIF
578 0 : h3 = h3 + MAX(-c(j),h1)
579 0 : 256 CONTINUE
580 : IF ((ABS(f-f0).LT.tol .OR. dnrm2_(n,s,1).LT.tol) .AND. h3.LT.tol &
581 0 : .AND. .NOT. badlin .AND. f .EQ. f) &
582 : THEN
583 0 : mode = 0
584 : ELSE
585 0 : mode = 8
586 : ENDIF
587 0 : GO TO 330
588 :
589 : !C CALL JACOBIAN AT CURRENT X
590 :
591 : !C UPDATE CHOLESKY-FACTORS OF HESSIAN MATRIX BY MODIFIED BFGS FORMULA
592 :
593 0 : 260 DO 270 i=1,n
594 0 : u(i) = g(i) - ddot_sl(m,a(1,i),1,r,1) - v(i)
595 0 : 270 CONTINUE
596 :
597 : !C L'*S
598 :
599 : k = 0
600 0 : DO 290 i=1,n
601 0 : h1 = ZERO
602 0 : k = k + 1
603 0 : DO 280 j=i+1,n
604 0 : k = k + 1
605 0 : h1 = h1 + l(k)*s(j)
606 0 : 280 CONTINUE
607 0 : v(i) = s(i) + h1
608 0 : 290 CONTINUE
609 :
610 : !C D*L'*S
611 :
612 : k = 1
613 0 : DO 300 i=1,n
614 0 : v(i) = l(k)*v(i)
615 0 : k = k + n1 - i
616 0 : 300 CONTINUE
617 :
618 : !C L*D*L'*S
619 :
620 0 : DO 320 i=n,1,-1
621 0 : h1 = ZERO
622 0 : k = i
623 0 : DO 310 j=1,i - 1
624 0 : h1 = h1 + l(k)*v(j)
625 0 : k = k + n - j
626 0 : 310 CONTINUE
627 0 : v(i) = v(i) + h1
628 0 : 320 CONTINUE
629 :
630 0 : h1 = ddot_sl(n,s,1,u,1)
631 0 : h2 = ddot_sl(n,s,1,v,1)
632 0 : h3 = 0.2d0*h2
633 0 : IF (h1.LT.h3) THEN
634 0 : h4 = (h2-h3)/(h2-h1)
635 0 : h1 = h3
636 0 : CALL dscal_sl(n, h4, u, 1)
637 0 : CALL daxpy_sl(n, one-h4, v, 1, u, 1)
638 : ENDIF
639 0 : IF (h1.EQ.0 .or. h2.EQ.0) THEN
640 : !C Singular update: reset hessian.
641 : GO TO 110
642 : end if
643 0 : CALL ldl(n, l, u, +one/h1, v)
644 0 : CALL ldl(n, l, v, -one/h2, u)
645 :
646 : !C END OF MAIN ITERATION
647 :
648 0 : GO TO 130
649 :
650 : !C END OF SLSQPB
651 :
652 : 330 continue
653 0 : END SUBROUTINE slsqpb
654 :
655 :
656 0 : SUBROUTINE lsq(m,meq,n,nl,la,l,g,a,b,xl,xu,x,y,w,jw,mode)
657 :
658 : !C MINIMIZE with respect to X
659 :
660 : !C ||E*X - F||
661 : !C 1/2 T
662 : !C WITH UPPER TRIANGULAR MATRIX E = +D *L ,
663 :
664 : !C -1/2 -1
665 : !C AND VECTOR F = -D *L *G,
666 :
667 : !C WHERE THE UNIT LOWER TRIDIANGULAR MATRIX L IS STORED COLUMNWISE
668 : !C DENSE IN THE N*(N+1)/2 ARRAY L WITH VECTOR D STORED IN ITS
669 : !C 'DIAGONAL' THUS SUBSTITUTING THE ONE-ELEMENTS OF L
670 :
671 : !C SUBJECT TO
672 :
673 : !C A(J)*X - B(J) = 0 , J=1,...,MEQ,
674 : !C A(J)*X - B(J) >=0, J=MEQ+1,...,M,
675 : !C XL(I) <= X(I) <= XU(I), I=1,...,N,
676 : !C ON ENTRY, THE USER HAS TO PROVIDE THE ARRAYS L, G, A, B, XL, XU.
677 : !C WITH DIMENSIONS: L(N*(N+1)/2), G(N), A(LA,N), B(M), XL(N), XU(N)
678 : !C THE WORKING ARRAY W MUST HAVE AT LEAST THE FOLLOWING DIMENSION:
679 : !c DIM(W) = (3*N+M)*(N+1) for LSQ
680 : !c +(N-MEQ+1)*(MINEQ+2) + 2*MINEQ for LSI
681 : !c +(N+MINEQ)*(N-MEQ) + 2*MEQ + N for LSEI
682 : !c with MINEQ = M - MEQ + 2*N
683 : !C ON RETURN, NO ARRAY WILL BE CHANGED BY THE SUBROUTINE.
684 : !C X STORES THE N-DIMENSIONAL SOLUTION VECTOR
685 : !C Y STORES THE VECTOR OF LAGRANGE MULTIPLIERS OF DIMENSION
686 : !C M+N+N (CONSTRAINTS+LOWER+UPPER BOUNDS)
687 : !C MODE IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING MEANINGS:
688 : !C MODE=1: SUCCESSFUL COMPUTATION
689 : !C 2: ERROR RETURN BECAUSE OF WRONG DIMENSIONS (N<1)
690 : !C 3: ITERATION COUNT EXCEEDED BY NNLS
691 : !C 4: INEQUALITY CONSTRAINTS INCOMPATIBLE
692 : !C 5: MATRIX E IS NOT OF FULL RANK
693 : !C 6: MATRIX C IS NOT OF FULL RANK
694 : !C 7: RANK DEFECT IN HFTI
695 :
696 : !c coded Dieter Kraft, april 1987
697 : !c revised march 1989
698 :
699 : DOUBLE PRECISION l,g,a,b,w,xl,xu,x,y, &
700 : diag,ZERO,one,xnorm
701 :
702 : INTEGER jw(*),i,ic,id,ie,IF,ig,ih,il,ip,iw,&
703 : i1,i2,i3,i4,la,m,meq,mineq,mode,m1,n,nl,n1,n2,n3,&
704 : nancnt,j
705 :
706 : DIMENSION a(la,n), b(la), g(n), l(nl),&
707 : w(*), x(n), xl(n), xu(n), y(m+n+n)
708 :
709 : DATA ZERO/0.0d0/, one/1.0d0/
710 :
711 0 : n1 = n + 1
712 0 : mineq = m - meq
713 0 : m1 = mineq + n + n
714 :
715 : !c determine whether to solve problem
716 : !c with inconsistent linerarization (n2=1)
717 : !c or not (n2=0)
718 :
719 0 : n2 = n1*n/2 + 1
720 0 : IF (n2.EQ.nl) THEN
721 : n2 = 0
722 : ELSE
723 0 : n2 = 1
724 : ENDIF
725 0 : n3 = n-n2
726 :
727 : !C RECOVER MATRIX E AND VECTOR F FROM L AND G
728 :
729 0 : i2 = 1
730 0 : i3 = 1
731 0 : i4 = 1
732 0 : ie = 1
733 0 : IF = n*n+1
734 0 : DO 10 i=1,n3
735 0 : i1 = n1-i
736 0 : diag = SQRT (l(i2))
737 0 : w(i3) = ZERO
738 0 : CALL dcopy_ (i1 , w(i3), 0, w(i3), 1)
739 0 : CALL dcopy_ (i1-n2, l(i2), 1, w(i3), n)
740 0 : CALL dscal_sl (i1-n2, diag, w(i3), n)
741 0 : w(i3) = diag
742 0 : w(IF-1+i) = (g(i) - ddot_sl (i-1, w(i4), 1, w(IF), 1))/diag
743 0 : i2 = i2 + i1 - n2
744 0 : i3 = i3 + n1
745 0 : i4 = i4 + n
746 0 : 10 CONTINUE
747 0 : IF (n2.EQ.1) THEN
748 0 : w(i3) = l(nl)
749 0 : w(i4) = ZERO
750 0 : CALL dcopy_ (n3, w(i4), 0, w(i4), 1)
751 0 : w(IF-1+n) = ZERO
752 : ENDIF
753 0 : CALL dscal_sl (n, - one, w(IF), 1)
754 :
755 0 : ic = IF + n
756 0 : id = ic + meq*n
757 :
758 0 : IF (meq .GT. 0) THEN
759 :
760 : !C RECOVER MATRIX C FROM UPPER PART OF A
761 :
762 0 : DO 20 i=1,meq
763 0 : CALL dcopy_ (n, a(i,1), la, w(ic-1+i), meq)
764 0 : 20 CONTINUE
765 :
766 : !C RECOVER VECTOR D FROM UPPER PART OF B
767 :
768 0 : CALL dcopy_ (meq, b(1), 1, w(id), 1)
769 0 : CALL dscal_sl (meq, - one, w(id), 1)
770 :
771 : ENDIF
772 :
773 0 : ig = id + meq
774 :
775 : !C RECOVER MATRIX G FROM LOWER PART OF A
776 : !C The matrix G(mineq+2*n,m1) is stored at w(ig)
777 : !C Not all rows will be filled if some of the upper/lower
778 : !C bounds are unbounded.
779 :
780 0 : IF (mineq .GT. 0) THEN
781 :
782 0 : DO 30 i=1,mineq
783 0 : CALL dcopy_ (n, a(meq+i,1), la, w(ig-1+i), m1)
784 0 : 30 CONTINUE
785 :
786 : ENDIF
787 :
788 0 : ih = ig + m1*n
789 0 : iw = ih + mineq + 2*n
790 :
791 0 : IF (mineq .GT. 0) THEN
792 :
793 : !C RECOVER H FROM LOWER PART OF B
794 : !C The vector H(mineq+2*n) is stored at w(ih)
795 :
796 0 : CALL dcopy_ (mineq, b(meq+1), 1, w(ih), 1)
797 0 : CALL dscal_sl (mineq, - one, w(ih), 1)
798 :
799 : ENDIF
800 :
801 : !C AUGMENT MATRIX G BY +I AND -I, AND,
802 : !C AUGMENT VECTOR H BY XL AND XU
803 : !C NaN value indicates no bound
804 :
805 0 : ip = ig + mineq
806 0 : il = ih + mineq
807 0 : nancnt = 0
808 :
809 0 : DO 40 i=1,n
810 0 : if (xl(i).eq.xl(i)) then
811 0 : w(il) = xl(i)
812 0 : do 41 j=1,n
813 0 : w(ip + m1*(j-1)) = 0
814 0 : 41 continue
815 0 : w(ip + m1*(i-1)) = 1
816 0 : ip = ip + 1
817 0 : il = il + 1
818 : else
819 0 : nancnt = nancnt + 1
820 : end if
821 0 : 40 CONTINUE
822 :
823 0 : DO 50 i=1,n
824 0 : if (xu(i).eq.xu(i)) then
825 0 : w(il) = -xu(i)
826 0 : do 51 j=1,n
827 0 : w(ip + m1*(j-1)) = 0
828 0 : 51 continue
829 0 : w(ip + m1*(i-1)) = -1
830 0 : ip = ip + 1
831 0 : il = il + 1
832 : else
833 0 : nancnt = nancnt + 1
834 : end if
835 0 : 50 CONTINUE
836 :
837 : CALL lsei (w(ic), w(id), w(ie), w(IF), w(ig), w(ih), MAX(1,meq),&
838 0 : meq, n, n, m1, m1-nancnt, n, x, xnorm, w(iw), jw, mode)
839 :
840 0 : IF (mode .EQ. 1) THEN
841 :
842 : !c restore Lagrange multipliers (only for user-defined variables)
843 :
844 0 : CALL dcopy_ (m, w(iw), 1, y(1), 1)
845 :
846 : !c set rest of the multipliers to nan (they are not used)
847 :
848 0 : IF (n3 .GT. 0) THEN
849 0 : y(m+1) = 0
850 0 : y(m+1) = 0 / y(m+1)
851 0 : do 60 i=m+2,m+n3+n3
852 0 : y(i) = y(m+1)
853 0 : 60 continue
854 : ENDIF
855 :
856 : ENDIF
857 0 : call bound(n, x, xl, xu)
858 :
859 : !C END OF SUBROUTINE LSQ
860 :
861 0 : END SUBROUTINE lsq
862 :
863 :
864 0 : SUBROUTINE lsei(c,d,e,f,g,h,lc,mc,LE,me,lg,mg,n,x,xnrm,w,jw,mode)
865 :
866 : !C FOR MODE=1, THE SUBROUTINE RETURNS THE SOLUTION X OF
867 : !C EQUALITY & INEQUALITY CONSTRAINED LEAST SQUARES PROBLEM LSEI :
868 :
869 : !C MIN ||E*X - F||
870 : !C X
871 :
872 : !C S.T. C*X = D,
873 : !C G*X >= H.
874 :
875 : !C USING QR DECOMPOSITION & ORTHOGONAL BASIS OF NULLSPACE OF C
876 : !C CHAPTER 23.6 OF LAWSON & HANSON: SOLVING LEAST SQUARES PROBLEMS.
877 :
878 : !C THE FOLLOWING DIMENSIONS OF THE ARRAYS DEFINING THE PROBLEM
879 : !C ARE NECESSARY
880 : !C DIM(E) : FORMAL (LE,N), ACTUAL (ME,N)
881 : !C DIM(F) : FORMAL (LE ), ACTUAL (ME )
882 : !C DIM(C) : FORMAL (LC,N), ACTUAL (MC,N)
883 : !C DIM(D) : FORMAL (LC ), ACTUAL (MC )
884 : !C DIM(G) : FORMAL (LG,N), ACTUAL (MG,N)
885 : !C DIM(H) : FORMAL (LG ), ACTUAL (MG )
886 : !C DIM(X) : FORMAL (N ), ACTUAL (N )
887 : !C DIM(W) : 2*MC+ME+(ME+MG)*(N-MC) for LSEI
888 : !C +(N-MC+1)*(MG+2)+2*MG for LSI
889 : !C DIM(JW): MAX(MG,L)
890 : !C ON ENTRY, THE USER HAS TO PROVIDE THE ARRAYS C, D, E, F, G, AND H.
891 : !C ON RETURN, ALL ARRAYS WILL BE CHANGED BY THE SUBROUTINE.
892 : !C X STORES THE SOLUTION VECTOR
893 : !C XNORM STORES THE RESIDUUM OF THE SOLUTION IN EUCLIDIAN NORM
894 : !C W STORES THE VECTOR OF LAGRANGE MULTIPLIERS IN ITS FIRST
895 : !C MC+MG ELEMENTS
896 : !C MODE IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING MEANINGS:
897 : !C MODE=1: SUCCESSFUL COMPUTATION
898 : !C 2: ERROR RETURN BECAUSE OF WRONG DIMENSIONS (N<1)
899 : !C 3: ITERATION COUNT EXCEEDED BY NNLS
900 : !C 4: INEQUALITY CONSTRAINTS INCOMPATIBLE
901 : !C 5: MATRIX E IS NOT OF FULL RANK
902 : !C 6: MATRIX C IS NOT OF FULL RANK
903 : !C 7: RANK DEFECT IN HFTI
904 :
905 : !C 18.5.1981, DIETER KRAFT, DFVLR OBERPFAFFENHOFEN
906 : !C 20.3.1987, DIETER KRAFT, DFVLR OBERPFAFFENHOFEN
907 :
908 : INTEGER jw(*),i,ie,IF,ig,iw,j,k,krank,l,lc,LE,lg,&
909 : mc,mc1,me,mg,mode,n
910 : DOUBLE PRECISION c(lc,n),e(LE,n),g(lg,n),d(lc),f(LE),h(lg),x(n),&
911 : w(*),t,xnrm,rnorm(1),epmach,ZERO
912 : DATA epmach/2.22d-16/,ZERO/0.0d+00/
913 :
914 0 : mode=2
915 0 : IF(mc.GT.n) GOTO 75
916 0 : l=n-mc
917 0 : mc1=mc+1
918 0 : iw=(l+1)*(mg+2)+2*mg+mc
919 0 : ie=iw+mc+1
920 0 : IF=ie+me*l
921 0 : ig=IF+me
922 :
923 : !C TRIANGULARIZE C AND APPLY FACTORS TO E AND G
924 :
925 0 : DO 10 i=1,mc
926 0 : j=MIN(i+1,lc)
927 0 : CALL h12(1,i,i+1,n,c(i,1),lc,w(iw+i),c(j,1),lc,1,mc-i)
928 0 : CALL h12(2,i,i+1,n,c(i,1),lc,w(iw+i),e ,LE,1,me)
929 0 : CALL h12(2,i,i+1,n,c(i,1),lc,w(iw+i),g ,lg,1,mg)
930 0 : 10 CONTINUE
931 :
932 : !C SOLVE C*X=D AND MODIFY F
933 :
934 0 : mode=6
935 0 : DO 15 i=1,mc
936 0 : IF(ABS(c(i,i)).LT.epmach) GOTO 75
937 0 : x(i)=(d(i)-ddot_sl(i-1,c(i,1),lc,x,1))/c(i,i)
938 0 : 15 CONTINUE
939 0 : mode=1
940 0 : w(mc1) = ZERO
941 0 : CALL dcopy_ (mg-mc,w(mc1),0,w(mc1),1)
942 :
943 0 : IF(mc.EQ.n) GOTO 50
944 :
945 0 : DO 20 i=1,me
946 0 : w(IF-1+i)=f(i)-ddot_sl(mc,e(i,1),LE,x,1)
947 0 : 20 CONTINUE
948 : !C STORE TRANSFORMED E & G
949 :
950 0 : DO 25 i=1,me
951 0 : CALL dcopy_(l,e(i,mc1),LE,w(ie-1+i),me)
952 0 : 25 CONTINUE
953 0 : DO 30 i=1,mg
954 0 : CALL dcopy_(l,g(i,mc1),lg,w(ig-1+i),mg)
955 0 : 30 CONTINUE
956 :
957 0 : IF(mg.GT.0) GOTO 40
958 :
959 : !C SOLVE LS WITHOUT INEQUALITY CONSTRAINTS
960 :
961 0 : mode=7
962 0 : k=MAX(LE,n)
963 0 : t=SQRT(epmach)
964 0 : CALL hfti (w(ie),me,me,l,w(IF),k,1,t,krank,rnorm,w,w(l+1),jw)
965 : !C HFTI IS MORE GENERIC, BUT WE ONLY CALL IT WITH NB=1, SO RETRIEVE THE
966 : !C SINGLE VALUE WE NEED FROM RNORM HERE
967 0 : xnrm = rnorm(1)
968 0 : CALL dcopy_(l,w(IF),1,x(mc1),1)
969 0 : IF(krank.NE.l) GOTO 75
970 0 : mode=1
971 0 : GOTO 50
972 : !C MODIFY H AND SOLVE INEQUALITY CONSTRAINED LS PROBLEM
973 :
974 0 : 40 DO 45 i=1,mg
975 0 : h(i)=h(i)-ddot_sl(mc,g(i,1),lg,x,1)
976 0 : 45 CONTINUE
977 : CALL lsi &
978 0 : (w(ie),w(IF),w(ig),h,me,me,mg,mg,l,x(mc1),xnrm,w(mc1),jw,mode)
979 0 : IF(mc.EQ.0) GOTO 75
980 0 : t=dnrm2_(mc,x,1)
981 0 : xnrm=SQRT(xnrm*xnrm+t*t)
982 0 : IF(mode.NE.1) GOTO 75
983 :
984 : !C SOLUTION OF ORIGINAL PROBLEM AND LAGRANGE MULTIPLIERS
985 :
986 0 : 50 DO 55 i=1,me
987 0 : f(i)=ddot_sl(n,e(i,1),LE,x,1)-f(i)
988 0 : 55 CONTINUE
989 0 : DO 60 i=1,mc
990 0 : d(i)=ddot_sl(me,e(1,i),1,f,1)-ddot_sl(mg,g(1,i),1,w(mc1),1)
991 0 : 60 CONTINUE
992 0 : DO 65 i=mc,1,-1
993 0 : CALL h12(2,i,i+1,n,c(i,1),lc,w(iw+i),x,1,1,1)
994 0 : 65 CONTINUE
995 0 : DO 70 i=mc,1,-1
996 0 : j=MIN(i+1,lc)
997 0 : w(i)=(d(i)-ddot_sl(mc-i,c(j,i),1,w(j),1))/c(i,i)
998 0 : 70 CONTINUE
999 :
1000 : !C END OF SUBROUTINE LSEI
1001 :
1002 : 75 CONTINUE
1003 :
1004 0 : END SUBROUTINE lsei
1005 :
1006 0 : SUBROUTINE lsi(e,f,g,h,LE,me,lg,mg,n,x,xnorm,w,jw,mode)
1007 :
1008 : !C FOR MODE=1, THE SUBROUTINE RETURNS THE SOLUTION X OF
1009 : !C INEQUALITY CONSTRAINED LINEAR LEAST SQUARES PROBLEM:
1010 :
1011 : !C MIN ||E*X-F||
1012 : !C X
1013 :
1014 : !C S.T. G*X >= H
1015 :
1016 : !C THE ALGORITHM IS BASED ON QR DECOMPOSITION AS DESCRIBED IN
1017 : !C CHAPTER 23.5 OF LAWSON & HANSON: SOLVING LEAST SQUARES PROBLEMS
1018 :
1019 : !C THE FOLLOWING DIMENSIONS OF THE ARRAYS DEFINING THE PROBLEM
1020 : !C ARE NECESSARY
1021 : !C DIM(E) : FORMAL (LE,N), ACTUAL (ME,N)
1022 : !C DIM(F) : FORMAL (LE ), ACTUAL (ME )
1023 : !C DIM(G) : FORMAL (LG,N), ACTUAL (MG,N)
1024 : !C DIM(H) : FORMAL (LG ), ACTUAL (MG )
1025 : !C DIM(X) : N
1026 : !C DIM(W) : (N+1)*(MG+2) + 2*MG
1027 : !C DIM(JW): LG
1028 : !C ON ENTRY, THE USER HAS TO PROVIDE THE ARRAYS E, F, G, AND H.
1029 : !C ON RETURN, ALL ARRAYS WILL BE CHANGED BY THE SUBROUTINE.
1030 : !C X STORES THE SOLUTION VECTOR
1031 : !C XNORM STORES THE RESIDUUM OF THE SOLUTION IN EUCLIDIAN NORM
1032 : !C W STORES THE VECTOR OF LAGRANGE MULTIPLIERS IN ITS FIRST
1033 : !C MG ELEMENTS
1034 : !C MODE IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING MEANINGS:
1035 : !C MODE=1: SUCCESSFUL COMPUTATION
1036 : !C 2: ERROR RETURN BECAUSE OF WRONG DIMENSIONS (N<1)
1037 : !C 3: ITERATION COUNT EXCEEDED BY NNLS
1038 : !C 4: INEQUALITY CONSTRAINTS INCOMPATIBLE
1039 : !C 5: MATRIX E IS NOT OF FULL RANK
1040 :
1041 : !C 03.01.1980, DIETER KRAFT: CODED
1042 : !C 20.03.1987, DIETER KRAFT: REVISED TO FORTRAN 77
1043 :
1044 : INTEGER i,j,LE,lg,me,mg,mode,n,jw(lg)
1045 : DOUBLE PRECISION e(LE,n),f(LE),g(lg,n),h(lg),x(n),w(*),&
1046 : xnorm,epmach,t,one
1047 : DATA epmach/2.22d-16/,one/1.0d+00/
1048 :
1049 : !C QR-FACTORS OF E AND APPLICATION TO F
1050 :
1051 0 : DO 10 i=1,n
1052 0 : j=MIN(i+1,n)
1053 0 : CALL h12(1,i,i+1,me,e(1,i),1,t,e(1,j),1,LE,n-i)
1054 0 : CALL h12(2,i,i+1,me,e(1,i),1,t,f ,1,1 ,1 )
1055 0 : 10 CONTINUE
1056 : !C TRANSFORM G AND H TO GET LEAST DISTANCE PROBLEM
1057 :
1058 0 : mode=5
1059 0 : DO 30 i=1,mg
1060 0 : DO 20 j=1,n
1061 0 : IF (.NOT.(ABS(e(j,j)).GE.epmach)) GOTO 50
1062 0 : g(i,j)=(g(i,j)-ddot_sl(j-1,g(i,1),lg,e(1,j),1))/e(j,j)
1063 0 : 20 CONTINUE
1064 0 : h(i)=h(i)-ddot_sl(n,g(i,1),lg,f,1)
1065 0 : 30 CONTINUE
1066 :
1067 : !C SOLVE LEAST DISTANCE PROBLEM
1068 :
1069 0 : CALL ldp(g,lg,mg,n,h,x,xnorm,w,jw,mode)
1070 0 : IF (mode.NE.1) GOTO 50
1071 :
1072 : !C SOLUTION OF ORIGINAL PROBLEM
1073 :
1074 0 : CALL daxpy_sl(n,one,f,1,x,1)
1075 0 : DO 40 i=n,1,-1
1076 0 : j=MIN(i+1,n)
1077 0 : x(i)=(x(i)-ddot_sl(n-i,e(i,j),LE,x(j),1))/e(i,i)
1078 0 : 40 CONTINUE
1079 0 : j=MIN(n+1,me)
1080 0 : t=dnrm2_(me-n,f(j),1)
1081 0 : xnorm=SQRT(xnorm*xnorm+t*t)
1082 :
1083 : !C END OF SUBROUTINE LSI
1084 :
1085 : 50 CONTINUE
1086 :
1087 0 : END SUBROUTINE lsi
1088 :
1089 0 : SUBROUTINE ldp(g,mg,m,n,h,x,xnorm,w,INDEX,mode)
1090 :
1091 : !C T
1092 : !C MINIMIZE 1/2 X X SUBJECT TO G * X >= H.
1093 :
1094 : !C C.L. LAWSON, R.J. HANSON: 'SOLVING LEAST SQUARES PROBLEMS'
1095 : !C PRENTICE HALL, ENGLEWOOD CLIFFS, NEW JERSEY, 1974.
1096 :
1097 : !C PARAMETER DESCRIPTION:
1098 :
1099 : !C G(),MG,M,N ON ENTRY G() STORES THE M BY N MATRIX OF
1100 : !C LINEAR INEQUALITY CONSTRAINTS. G() HAS FIRST
1101 : !C DIMENSIONING PARAMETER MG
1102 : !C H() ON ENTRY H() STORES THE M VECTOR H REPRESENTING
1103 : !C THE RIGHT SIDE OF THE INEQUALITY SYSTEM
1104 :
1105 : !C REMARK: G(),H() WILL NOT BE CHANGED DURING CALCULATIONS BY LDP
1106 :
1107 : !C X() ON ENTRY X() NEED NOT BE INITIALIZED.
1108 : !C ON EXIT X() STORES THE SOLUTION VECTOR X IF MODE=1.
1109 : !C XNORM ON EXIT XNORM STORES THE EUCLIDIAN NORM OF THE
1110 : !C SOLUTION VECTOR IF COMPUTATION IS SUCCESSFUL
1111 : !C W() W IS A ONE DIMENSIONAL WORKING SPACE, THE LENGTH
1112 : !C OF WHICH SHOULD BE AT LEAST (M+2)*(N+1) + 2*M
1113 : !C ON EXIT W() STORES THE LAGRANGE MULTIPLIERS
1114 : !C ASSOCIATED WITH THE CONSTRAINTS
1115 : !C AT THE SOLUTION OF PROBLEM LDP
1116 : !C INDEX() INDEX() IS A ONE DIMENSIONAL INTEGER WORKING SPACE
1117 : !C OF LENGTH AT LEAST M
1118 : !C MODE MODE IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING
1119 : !C MEANINGS:
1120 : !C MODE=1: SUCCESSFUL COMPUTATION
1121 : !C 2: ERROR RETURN BECAUSE OF WRONG DIMENSIONS (N.LE.0)
1122 : !C 3: ITERATION COUNT EXCEEDED BY NNLS
1123 : !C 4: INEQUALITY CONSTRAINTS INCOMPATIBLE
1124 :
1125 : DOUBLE PRECISION g,h,x,xnorm,w,u,v,&
1126 : ZERO,one,fac,rnorm,diff
1127 : INTEGER INDEX,i,IF,iw,iwdual,iy,iz,j,m,mg,mode,n,n1
1128 : DIMENSION g(mg,n),h(m),x(n),w(*),INDEX(m)
1129 : diff(u,v)= u-v
1130 : DATA ZERO,one/0.0d0,1.0d0/
1131 :
1132 0 : mode=2
1133 0 : IF(n.LE.0) GOTO 50
1134 :
1135 : !C STATE DUAL PROBLEM
1136 :
1137 0 : mode=1
1138 0 : x(1)=ZERO
1139 0 : CALL dcopy_(n,x(1),0,x,1)
1140 0 : xnorm=ZERO
1141 0 : IF(m.EQ.0) GOTO 50
1142 : iw=0
1143 0 : DO 20 j=1,m
1144 0 : DO 10 i=1,n
1145 0 : iw=iw+1
1146 0 : w(iw)=g(j,i)
1147 0 : 10 CONTINUE
1148 0 : iw=iw+1
1149 0 : w(iw)=h(j)
1150 0 : 20 CONTINUE
1151 : IF=iw+1
1152 0 : DO 30 i=1,n
1153 0 : iw=iw+1
1154 0 : w(iw)=ZERO
1155 0 : 30 CONTINUE
1156 0 : w(iw+1)=one
1157 0 : n1=n+1
1158 0 : iz=iw+2
1159 0 : iy=iz+n1
1160 0 : iwdual=iy+m
1161 :
1162 : !C SOLVE DUAL PROBLEM
1163 :
1164 0 : CALL nnls (w,n1,n1,m,w(IF),w(iy),rnorm,w(iwdual),w(iz),INDEX,mode)
1165 :
1166 0 : IF(mode.NE.1) GOTO 50
1167 0 : mode=4
1168 0 : IF(rnorm.LE.ZERO) GOTO 50
1169 :
1170 : !C COMPUTE SOLUTION OF PRIMAL PROBLEM
1171 :
1172 0 : fac=one-ddot_sl(m,h,1,w(iy),1)
1173 0 : IF(.NOT.(diff(one+fac,one).GT.ZERO)) GOTO 50
1174 0 : mode=1
1175 0 : fac=one/fac
1176 0 : DO 40 j=1,n
1177 0 : x(j)=fac*ddot_sl(m,g(1,j),1,w(iy),1)
1178 0 : 40 CONTINUE
1179 0 : xnorm=dnrm2_(n,x,1)
1180 :
1181 : !C COMPUTE LAGRANGE MULTIPLIERS FOR PRIMAL PROBLEM
1182 :
1183 0 : w(1)=ZERO
1184 0 : CALL dcopy_(m,w(1),0,w,1)
1185 0 : CALL daxpy_sl(m,fac,w(iy),1,w,1)
1186 :
1187 : !C END OF SUBROUTINE LDP
1188 :
1189 : 50 CONTINUE
1190 :
1191 0 : END SUBROUTINE ldp
1192 :
1193 :
1194 0 : SUBROUTINE nnls (a, mda, m, n, b, x, rnorm, w, z, INDEX, mode)
1195 :
1196 : !C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY:
1197 : !C 'SOLVING LEAST SQUARES PROBLEMS'. PRENTICE-HALL.1974
1198 :
1199 : !C ********** NONNEGATIVE LEAST SQUARES **********
1200 :
1201 : !C GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN
1202 : !C N-VECTOR, X, WHICH SOLVES THE LEAST SQUARES PROBLEM
1203 :
1204 : !C A*X = B SUBJECT TO X >= 0
1205 :
1206 : !C A(),MDA,M,N
1207 : !C MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE ARRAY,A().
1208 : !C ON ENTRY A() CONTAINS THE M BY N MATRIX,A.
1209 : !C ON EXIT A() CONTAINS THE PRODUCT Q*A,
1210 : !C WHERE Q IS AN M BY M ORTHOGONAL MATRIX GENERATED
1211 : !C IMPLICITLY BY THIS SUBROUTINE.
1212 : !C EITHER M>=N OR M<N IS PERMISSIBLE.
1213 : !C THERE IS NO RESTRICTION ON THE RANK OF A.
1214 : !C B() ON ENTRY B() CONTAINS THE M-VECTOR, B.
1215 : !C ON EXIT B() CONTAINS Q*B.
1216 : !C X() ON ENTRY X() NEED NOT BE INITIALIZED.
1217 : !C ON EXIT X() WILL CONTAIN THE SOLUTION VECTOR.
1218 : !C RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE
1219 : !C RESIDUAL VECTOR.
1220 : !C W() AN N-ARRAY OF WORKING SPACE.
1221 : !C ON EXIT W() WILL CONTAIN THE DUAL SOLUTION VECTOR.
1222 : !C W WILL SATISFY W(I)=0 FOR ALL I IN SET P
1223 : !C AND W(I)<=0 FOR ALL I IN SET Z
1224 : !C Z() AN M-ARRAY OF WORKING SPACE.
1225 : !C INDEX()AN INTEGER WORKING ARRAY OF LENGTH AT LEAST N.
1226 : !C ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS
1227 : !C P AND Z AS FOLLOWS:
1228 : !C INDEX(1) THRU INDEX(NSETP) = SET P.
1229 : !C INDEX(IZ1) THRU INDEX (IZ2) = SET Z.
1230 : !C IZ1=NSETP + 1 = NPP1, IZ2=N.
1231 : !C MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING MEANING:
1232 : !C 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY.
1233 : !C 2 THE DIMENSIONS OF THE PROBLEM ARE WRONG,
1234 : !C EITHER M <= 0 OR N <= 0.
1235 : !C 3 ITERATION COUNT EXCEEDED, MORE THAN 3*N ITERATIONS.
1236 :
1237 : INTEGER i,ii,ip,iter,itmax,iz,izmax,iz1,iz2,j,jj,jz,&
1238 : k,l,m,mda,mode,n,npp1,nsetp,INDEX(n)
1239 :
1240 : DOUBLE PRECISION a(mda,n),b(m),x(n),w(n),z(m),asave,diff,&
1241 : factor,ZERO,one,wmax,alpha,&
1242 : c,s,t,u,v,up,rnorm,unorm
1243 :
1244 : diff(u,v)= u-v
1245 :
1246 : DATA ZERO,one,factor/0.0d0,1.0d0,1.0d-2/
1247 :
1248 : !c revised Dieter Kraft, March 1983
1249 :
1250 0 : mode=2
1251 0 : IF(m.LE.0.OR.n.LE.0) GOTO 290
1252 0 : mode=1
1253 0 : iter=0
1254 0 : itmax=3*n
1255 :
1256 : !C STEP ONE (INITIALIZE)
1257 :
1258 0 : DO 100 i=1,n
1259 0 : INDEX(i)=i
1260 0 : 100 CONTINUE
1261 0 : iz1=1
1262 0 : iz2=n
1263 0 : nsetp=0
1264 0 : npp1=1
1265 0 : x(1)=ZERO
1266 0 : CALL dcopy_(n,x(1),0,x,1)
1267 :
1268 : !C STEP TWO (COMPUTE DUAL VARIABLES)
1269 : !C .....ENTRY LOOP A
1270 :
1271 0 : 110 IF(iz1.GT.iz2.OR.nsetp.GE.m) GOTO 280
1272 0 : DO 120 iz=iz1,iz2
1273 0 : j=INDEX(iz)
1274 0 : w(j)=ddot_sl(m-nsetp,a(npp1,j),1,b(npp1),1)
1275 0 : 120 CONTINUE
1276 :
1277 : !C STEP THREE (TEST DUAL VARIABLES)
1278 :
1279 0 : 130 wmax=ZERO
1280 0 : DO 140 iz=iz1,iz2
1281 0 : j=INDEX(iz)
1282 0 : IF(w(j).LE.wmax) GOTO 140
1283 0 : wmax=w(j)
1284 0 : izmax=iz
1285 0 : 140 CONTINUE
1286 :
1287 : !C .....EXIT LOOP A
1288 :
1289 0 : IF(wmax.LE.ZERO) GOTO 280
1290 0 : iz=izmax
1291 0 : j=INDEX(iz)
1292 :
1293 : !C STEP FOUR (TEST INDEX J FOR LINEAR DEPENDENCY)
1294 :
1295 0 : asave=a(npp1,j)
1296 0 : CALL h12(1,npp1,npp1+1,m,a(1,j),1,up,z,1,1,0)
1297 0 : unorm=dnrm2_(nsetp,a(1,j),1)
1298 0 : t=factor*ABS(a(npp1,j))
1299 0 : IF(diff(unorm+t,unorm).LE.ZERO) GOTO 150
1300 0 : CALL dcopy_(m,b,1,z,1)
1301 0 : CALL h12(2,npp1,npp1+1,m,a(1,j),1,up,z,1,1,1)
1302 0 : IF(z(npp1)/a(npp1,j).GT.ZERO) GOTO 160
1303 0 : 150 a(npp1,j)=asave
1304 0 : w(j)=ZERO
1305 0 : GOTO 130
1306 : !C STEP FIVE (ADD COLUMN)
1307 :
1308 0 : 160 CALL dcopy_(m,z,1,b,1)
1309 0 : INDEX(iz)=INDEX(iz1)
1310 0 : INDEX(iz1)=j
1311 0 : iz1=iz1+1
1312 0 : nsetp=npp1
1313 0 : npp1=npp1+1
1314 0 : DO 170 jz=iz1,iz2
1315 0 : jj=INDEX(jz)
1316 0 : CALL h12(2,nsetp,npp1,m,a(1,j),1,up,a(1,jj),1,mda,1)
1317 0 : 170 CONTINUE
1318 0 : k=MIN(npp1,mda)
1319 0 : w(j)=ZERO
1320 0 : CALL dcopy_(m-nsetp,w(j),0,a(k,j),1)
1321 :
1322 : !C STEP SIX (SOLVE LEAST SQUARES SUB-PROBLEM)
1323 : !C .....ENTRY LOOP B
1324 :
1325 0 : 180 DO 200 ip=nsetp,1,-1
1326 0 : IF(ip.EQ.nsetp) GOTO 190
1327 0 : CALL daxpy_sl(ip,-z(ip+1),a(1,jj),1,z,1)
1328 0 : 190 jj=INDEX(ip)
1329 0 : z(ip)=z(ip)/a(ip,jj)
1330 0 : 200 CONTINUE
1331 0 : iter=iter+1
1332 0 : IF(iter.LE.itmax) GOTO 220
1333 0 : 210 mode=3
1334 0 : GOTO 280
1335 : !C STEP SEVEN TO TEN (STEP LENGTH ALGORITHM)
1336 :
1337 0 : 220 alpha=one
1338 0 : jj=0
1339 0 : DO 230 ip=1,nsetp
1340 0 : IF(z(ip).GT.ZERO) GOTO 230
1341 0 : l=INDEX(ip)
1342 0 : t=-x(l)/(z(ip)-x(l))
1343 0 : IF(alpha.LT.t) GOTO 230
1344 0 : alpha=t
1345 0 : jj=ip
1346 0 : 230 CONTINUE
1347 0 : DO 240 ip=1,nsetp
1348 0 : l=INDEX(ip)
1349 0 : x(l)=(one-alpha)*x(l) + alpha*z(ip)
1350 0 : 240 CONTINUE
1351 :
1352 : !C .....EXIT LOOP B
1353 :
1354 0 : IF(jj.EQ.0) GOTO 110
1355 :
1356 : !C STEP ELEVEN (DELETE COLUMN)
1357 :
1358 0 : i=INDEX(jj)
1359 0 : 250 x(i)=ZERO
1360 0 : jj=jj+1
1361 0 : DO 260 j=jj,nsetp
1362 0 : ii=INDEX(j)
1363 0 : INDEX(j-1)=ii
1364 0 : CALL dsrotg(a(j-1,ii),a(j,ii),c,s)
1365 0 : t=a(j-1,ii)
1366 0 : CALL dsrot(n,a(j-1,1),mda,a(j,1),mda,c,s)
1367 0 : a(j-1,ii)=t
1368 0 : a(j,ii)=ZERO
1369 0 : CALL dsrot(1,b(j-1),1,b(j),1,c,s)
1370 0 : 260 CONTINUE
1371 0 : npp1=nsetp
1372 0 : nsetp=nsetp-1
1373 0 : iz1=iz1-1
1374 0 : INDEX(iz1)=i
1375 0 : IF(nsetp.LE.0) GOTO 210
1376 0 : DO 270 jj=1,nsetp
1377 0 : i=INDEX(jj)
1378 0 : IF(x(i).LE.ZERO) GOTO 250
1379 0 : 270 CONTINUE
1380 0 : CALL dcopy_(m,b,1,z,1)
1381 0 : GOTO 180
1382 : !C STEP TWELVE (SOLUTION)
1383 :
1384 0 : 280 k=MIN(npp1,m)
1385 0 : rnorm=dnrm2_(m-nsetp,b(k),1)
1386 0 : IF(npp1.GT.m) THEN
1387 0 : w(1)=ZERO
1388 0 : CALL dcopy_(n,w(1),0,w,1)
1389 : ENDIF
1390 :
1391 : !C END OF SUBROUTINE NNLS
1392 :
1393 : 290 CONTINUE
1394 :
1395 0 : END SUBROUTINE nnls
1396 :
1397 0 : SUBROUTINE hfti(a,mda,m,n,b,mdb,nb,tau,krank,rnorm,h,g,ip)
1398 :
1399 : !C RANK-DEFICIENT LEAST SQUARES ALGORITHM AS DESCRIBED IN:
1400 : !C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12
1401 : !C TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974
1402 :
1403 : !C A(*,*),MDA,M,N THE ARRAY A INITIALLY CONTAINS THE M x N MATRIX A
1404 : !C OF THE LEAST SQUARES PROBLEM AX = B.
1405 : !C THE FIRST DIMENSIONING PARAMETER MDA MUST SATISFY
1406 : !C MDA >= M. EITHER M >= N OR M < N IS PERMITTED.
1407 : !C THERE IS NO RESTRICTION ON THE RANK OF A.
1408 : !C THE MATRIX A WILL BE MODIFIED BY THE SUBROUTINE.
1409 : !C B(*,*),MDB,NB IF NB = 0 THE SUBROUTINE WILL MAKE NO REFERENCE
1410 : !C TO THE ARRAY B. IF NB > 0 THE ARRAY B() MUST
1411 : !C INITIALLY CONTAIN THE M x NB MATRIX B OF THE
1412 : !C THE LEAST SQUARES PROBLEM AX = B AND ON RETURN
1413 : !C THE ARRAY B() WILL CONTAIN THE N x NB SOLUTION X.
1414 : !C IF NB>1 THE ARRAY B() MUST BE DOUBLE SUBSCRIPTED
1415 : !C WITH FIRST DIMENSIONING PARAMETER MDB>=MAX(M,N),
1416 : !C IF NB=1 THE ARRAY B() MAY BE EITHER SINGLE OR
1417 : !C DOUBLE SUBSCRIPTED.
1418 : !C TAU ABSOLUTE TOLERANCE PARAMETER FOR PSEUDORANK
1419 : !C DETERMINATION, PROVIDED BY THE USER.
1420 : !C KRANK PSEUDORANK OF A, SET BY THE SUBROUTINE.
1421 : !C RNORM ON EXIT, RNORM(J) WILL CONTAIN THE EUCLIDIAN
1422 : !C NORM OF THE RESIDUAL VECTOR FOR THE PROBLEM
1423 : !C DEFINED BY THE J-TH COLUMN VECTOR OF THE ARRAY B.
1424 : !C H(), G() ARRAYS OF WORKING SPACE OF LENGTH >= N.
1425 : !C IP() INTEGER ARRAY OF WORKING SPACE OF LENGTH >= N
1426 : !C RECORDING PERMUTATION INDICES OF COLUMN VECTORS
1427 :
1428 : INTEGER i,j,jb,k,kp1,krank,l,ldiag,lmax,m,&
1429 : mda,mdb,n,nb,ip(n)
1430 : DOUBLE PRECISION a(mda,n),b(mdb,nb),h(n),g(n),rnorm(nb),factor,&
1431 : tau,ZERO,hmax,diff,tmp,u,v
1432 : diff(u,v)= u-v
1433 : DATA ZERO/0.0d0/, factor/1.0d-3/
1434 :
1435 0 : k=0
1436 0 : ldiag=MIN(m,n)
1437 0 : IF(ldiag.LE.0) GOTO 270
1438 :
1439 : !C COMPUTE LMAX
1440 :
1441 0 : DO 80 j=1,ldiag
1442 0 : IF(j.EQ.1) GOTO 20
1443 : lmax=j
1444 0 : DO 10 l=j,n
1445 0 : h(l)=h(l)-a(j-1,l)**2
1446 0 : IF(h(l).GT.h(lmax)) lmax=l
1447 0 : 10 CONTINUE
1448 0 : IF(diff(hmax+factor*h(lmax),hmax).GT.ZERO)&
1449 : GOTO 50
1450 0 : 20 lmax=j
1451 0 : DO 40 l=j,n
1452 0 : h(l)=ZERO
1453 0 : DO 30 i=j,m
1454 0 : h(l)=h(l)+a(i,l)**2
1455 0 : 30 CONTINUE
1456 0 : IF(h(l).GT.h(lmax)) lmax=l
1457 0 : 40 CONTINUE
1458 0 : hmax=h(lmax)
1459 :
1460 : !C COLUMN INTERCHANGES IF NEEDED
1461 :
1462 0 : 50 ip(j)=lmax
1463 0 : IF(ip(j).EQ.j) GOTO 70
1464 0 : DO 60 i=1,m
1465 0 : tmp=a(i,j)
1466 0 : a(i,j)=a(i,lmax)
1467 0 : a(i,lmax)=tmp
1468 0 : 60 CONTINUE
1469 0 : h(lmax)=h(j)
1470 :
1471 : !C J-TH TRANSFORMATION AND APPLICATION TO A AND B
1472 :
1473 0 : 70 i=MIN(j+1,n)
1474 0 : CALL h12(1,j,j+1,m,a(1,j),1,h(j),a(1,i),1,mda,n-j)
1475 0 : CALL h12(2,j,j+1,m,a(1,j),1,h(j),b,1,mdb,nb)
1476 0 : 80 CONTINUE
1477 :
1478 : !C DETERMINE PSEUDORANK
1479 :
1480 0 : DO 90 j=1,ldiag
1481 0 : IF(ABS(a(j,j)).LE.tau) GOTO 100
1482 0 : 90 CONTINUE
1483 : k=ldiag
1484 0 : GOTO 110
1485 0 : 100 k=j-1
1486 0 : 110 kp1=k+1
1487 :
1488 : !C NORM OF RESIDUALS
1489 :
1490 0 : DO 130 jb=1,nb
1491 0 : rnorm(jb)=dnrm2_(m-k,b(kp1,jb),1)
1492 0 : 130 CONTINUE
1493 0 : IF(k.GT.0) GOTO 160
1494 0 : DO 150 jb=1,nb
1495 0 : DO 149 i=1,n
1496 0 : b(i,jb)=ZERO
1497 0 : 149 CONTINUE
1498 0 : 150 CONTINUE
1499 0 : GOTO 270
1500 0 : 160 IF(k.EQ.n) GOTO 180
1501 :
1502 : !C HOUSEHOLDER DECOMPOSITION OF FIRST K ROWS
1503 :
1504 0 : DO 170 i=k,1,-1
1505 0 : CALL h12(1,i,kp1,n,a(i,1),mda,g(i),a,mda,1,i-1)
1506 0 : 170 CONTINUE
1507 0 : 180 DO 250 jb=1,nb
1508 :
1509 : !C SOLVE K*K TRIANGULAR SYSTEM
1510 :
1511 0 : DO 210 i=k,1,-1
1512 0 : j=MIN(i+1,n)
1513 0 : b(i,jb)=(b(i,jb)-ddot_sl(k-i,a(i,j),mda,b(j,jb),1))/a(i,i)
1514 0 : 210 CONTINUE
1515 :
1516 : !C COMPLETE SOLUTION VECTOR
1517 :
1518 0 : IF(k.EQ.n) GOTO 240
1519 0 : DO 220 j=kp1,n
1520 0 : b(j,jb)=ZERO
1521 0 : 220 CONTINUE
1522 0 : DO 230 i=1,k
1523 0 : CALL h12(2,i,kp1,n,a(i,1),mda,g(i),b(1,jb),1,mdb,1)
1524 0 : 230 CONTINUE
1525 :
1526 : !C REORDER SOLUTION ACCORDING TO PREVIOUS COLUMN INTERCHANGES
1527 :
1528 0 : 240 DO 249 j=ldiag,1,-1
1529 0 : IF(ip(j).EQ.j) GOTO 249
1530 0 : l=ip(j)
1531 0 : tmp=b(l,jb)
1532 0 : b(l,jb)=b(j,jb)
1533 0 : b(j,jb)=tmp
1534 0 : 249 CONTINUE
1535 0 : 250 CONTINUE
1536 0 : 270 krank=k
1537 0 : END SUBROUTINE hfti
1538 :
1539 0 : SUBROUTINE h12 (mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv)
1540 :
1541 : !C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12
1542 : !C TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974
1543 :
1544 : !C CONSTRUCTION AND/OR APPLICATION OF A SINGLE
1545 : !C HOUSEHOLDER TRANSFORMATION Q = I + U*(U**T)/B
1546 :
1547 : !C MODE = 1 OR 2 TO SELECT ALGORITHM H1 OR H2 .
1548 : !C LPIVOT IS THE INDEX OF THE PIVOT ELEMENT.
1549 : !C L1,M IF L1 <= M THE TRANSFORMATION WILL BE CONSTRUCTED TO
1550 : !C ZERO ELEMENTS INDEXED FROM L1 THROUGH M.
1551 : !C IF L1 > M THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION.
1552 : !C U(),IUE,UP
1553 : !C ON ENTRY TO H1 U() STORES THE PIVOT VECTOR.
1554 : !C IUE IS THE STORAGE INCREMENT BETWEEN ELEMENTS.
1555 : !C ON EXIT FROM H1 U() AND UP STORE QUANTITIES DEFINING
1556 : !C THE VECTOR U OF THE HOUSEHOLDER TRANSFORMATION.
1557 : !C ON ENTRY TO H2 U() AND UP
1558 : !C SHOULD STORE QUANTITIES PREVIOUSLY COMPUTED BY H1.
1559 : !C THESE WILL NOT BE MODIFIED BY H2.
1560 : !C C() ON ENTRY TO H1 OR H2 C() STORES A MATRIX WHICH WILL BE
1561 : !C REGARDED AS A SET OF VECTORS TO WHICH THE HOUSEHOLDER
1562 : !C TRANSFORMATION IS TO BE APPLIED.
1563 : !C ON EXIT C() STORES THE SET OF TRANSFORMED VECTORS.
1564 : !C ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C().
1565 : !C ICV STORAGE INCREMENT BETWEEN VECTORS IN C().
1566 : !C NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED.
1567 : !C IF NCV <= 0 NO OPERATIONS WILL BE DONE ON C().
1568 :
1569 : INTEGER incr, ice, icv, iue, lpivot, l1, mode, ncv
1570 : INTEGER i, i2, i3, i4, j, m
1571 : DOUBLE PRECISION u,up,c,cl,clinv,b,sm,one,ZERO
1572 : DIMENSION u(iue,*), c(*)
1573 : DATA one/1.0d+00/, ZERO/0.0d+00/
1574 :
1575 0 : IF (0.GE.lpivot.OR.lpivot.GE.l1.OR.l1.GT.m) GOTO 80
1576 0 : cl=ABS(u(1,lpivot))
1577 0 : IF (mode.EQ.2) GOTO 30
1578 :
1579 : !C ****** CONSTRUCT THE TRANSFORMATION ******
1580 :
1581 0 : DO 10 j=l1,m
1582 0 : sm=ABS(u(1,j))
1583 0 : cl=MAX(sm,cl)
1584 0 : 10 CONTINUE
1585 0 : IF (cl.LE.ZERO) GOTO 80
1586 0 : clinv=one/cl
1587 0 : sm=(u(1,lpivot)*clinv)**2
1588 0 : DO 20 j=l1,m
1589 0 : sm=sm+(u(1,j)*clinv)**2
1590 0 : 20 CONTINUE
1591 0 : cl=cl*SQRT(sm)
1592 0 : IF (u(1,lpivot).GT.ZERO) cl=-cl
1593 0 : up=u(1,lpivot)-cl
1594 0 : u(1,lpivot)=cl
1595 0 : GOTO 40
1596 : !C ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C ******
1597 :
1598 0 : 30 IF (cl.LE.ZERO) GOTO 80
1599 0 : 40 IF (ncv.LE.0) GOTO 80
1600 0 : b=up*u(1,lpivot)
1601 0 : IF (b.GE.ZERO) GOTO 80
1602 0 : b=one/b
1603 0 : i2=1-icv+ice*(lpivot-1)
1604 0 : incr=ice*(l1-lpivot)
1605 0 : DO 70 j=1,ncv
1606 0 : i2=i2+icv
1607 0 : i3=i2+incr
1608 0 : i4=i3
1609 0 : sm=c(i2)*up
1610 0 : DO 50 i=l1,m
1611 0 : sm=sm+c(i3)*u(1,i)
1612 0 : i3=i3+ice
1613 0 : 50 CONTINUE
1614 0 : IF (sm.EQ.ZERO) GOTO 70
1615 0 : sm=sm*b
1616 0 : c(i2)=c(i2)+sm*up
1617 0 : DO 60 i=l1,m
1618 0 : c(i4)=c(i4)+sm*u(1,i)
1619 0 : i4=i4+ice
1620 0 : 60 CONTINUE
1621 0 : 70 CONTINUE
1622 : 80 CONTINUE
1623 :
1624 0 : END SUBROUTINE h12
1625 :
1626 0 : SUBROUTINE ldl (n,a,z,sigma,w)
1627 : !C LDL LDL' - RANK-ONE - UPDATE
1628 :
1629 : !C PURPOSE:
1630 : !C UPDATES THE LDL' FACTORS OF MATRIX A BY RANK-ONE MATRIX
1631 : !C SIGMA*Z*Z'
1632 :
1633 : !C INPUT ARGUMENTS: (* MEANS PARAMETERS ARE CHANGED DURING EXECUTION)
1634 : !C N : ORDER OF THE COEFFICIENT MATRIX A
1635 : !C * A : POSITIVE DEFINITE MATRIX OF DIMENSION N;
1636 : !C ONLY THE LOWER TRIANGLE IS USED AND IS STORED COLUMN BY
1637 : !C COLUMN AS ONE DIMENSIONAL ARRAY OF DIMENSION N*(N+1)/2.
1638 : !C * Z : VECTOR OF DIMENSION N OF UPDATING ELEMENTS
1639 : !C SIGMA : SCALAR FACTOR BY WHICH THE MODIFYING DYADE Z*Z' IS
1640 : !C MULTIPLIED
1641 :
1642 : !C OUTPUT ARGUMENTS:
1643 : !C A : UPDATED LDL' FACTORS
1644 :
1645 : !C WORKING ARRAY:
1646 : !C W : VECTOR OP DIMENSION N (USED ONLY IF SIGMA .LT. ZERO)
1647 :
1648 : !C METHOD:
1649 : !C THAT OF FLETCHER AND POWELL AS DESCRIBED IN :
1650 : !C FLETCHER,R.,(1974) ON THE MODIFICATION OF LDL' FACTORIZATION.
1651 : !C POWELL,M.J.D. MATH.COMPUTATION 28, 1067-1078.
1652 :
1653 : !C IMPLEMENTED BY:
1654 : !C KRAFT,D., DFVLR - INSTITUT FUER DYNAMIK DER FLUGSYSTEME
1655 : !C D-8031 OBERPFAFFENHOFEN
1656 :
1657 : !C STATUS: 15. JANUARY 1980
1658 :
1659 : !C SUBROUTINES REQUIRED: NONE
1660 :
1661 : INTEGER i, ij, j, n
1662 : DOUBLE PRECISION a(*), t, v, w(*), z(*), u, tp, one, beta, four,&
1663 : ZERO, alpha, delta, gamma, sigma, epmach
1664 : DATA ZERO, one, four, epmach /0.0d0, 1.0d0, 4.0d0, 2.22d-16/
1665 :
1666 0 : IF(sigma.EQ.ZERO) GOTO 280
1667 0 : ij=1
1668 0 : t=one/sigma
1669 0 : IF(sigma.GT.ZERO) GOTO 220
1670 : !C PREPARE NEGATIVE UPDATE
1671 0 : DO 150 i=1,n
1672 0 : w(i)=z(i)
1673 0 : 150 CONTINUE
1674 0 : DO 170 i=1,n
1675 0 : v=w(i)
1676 0 : t=t+v*v/a(ij)
1677 0 : DO 160 j=i+1,n
1678 0 : ij=ij+1
1679 0 : w(j)=w(j)-v*a(ij)
1680 0 : 160 CONTINUE
1681 0 : ij=ij+1
1682 0 : 170 CONTINUE
1683 0 : IF(t.GE.ZERO) t=epmach/sigma
1684 0 : DO 210 i=1,n
1685 0 : j=n+1-i
1686 0 : ij=ij-i
1687 0 : u=w(j)
1688 0 : w(j)=t
1689 0 : t=t-u*u/a(ij)
1690 0 : 210 CONTINUE
1691 : 220 CONTINUE
1692 : !C HERE UPDATING BEGINS
1693 0 : DO 270 i=1,n
1694 0 : v=z(i)
1695 0 : delta=v/a(ij)
1696 0 : IF(sigma.LT.ZERO) tp=w(i)
1697 0 : IF(sigma.GT.ZERO) tp=t+delta*v
1698 0 : alpha=tp/t
1699 0 : a(ij)=alpha*a(ij)
1700 0 : IF(i.EQ.n) GOTO 280
1701 0 : beta=delta/tp
1702 0 : IF(alpha.GT.four) GOTO 240
1703 0 : DO 230 j=i+1,n
1704 0 : ij=ij+1
1705 0 : z(j)=z(j)-v*a(ij)
1706 0 : a(ij)=a(ij)+beta*z(j)
1707 0 : 230 CONTINUE
1708 0 : GOTO 260
1709 0 : 240 gamma=t/tp
1710 0 : DO 250 j=i+1,n
1711 0 : ij=ij+1
1712 0 : u=a(ij)
1713 0 : a(ij)=gamma*u+beta*z(j)
1714 0 : z(j)=z(j)-v*u
1715 0 : 250 CONTINUE
1716 0 : 260 ij=ij+1
1717 0 : t=tp
1718 0 : 270 CONTINUE
1719 0 : 280 RETURN
1720 : !C END OF LDL
1721 : END SUBROUTINE ldl
1722 :
1723 0 : DOUBLE PRECISION FUNCTION linmin (mode, ax, bx, f, tol)
1724 : !C LINMIN LINESEARCH WITHOUT DERIVATIVES
1725 :
1726 : !C PURPOSE:
1727 :
1728 : !C TO FIND THE ARGUMENT LINMIN WHERE THE FUNCTION F TAKES IT'S MINIMUM
1729 : !C ON THE INTERVAL AX, BX.
1730 : !C COMBINATION OF GOLDEN SECTION AND SUCCESSIVE QUADRATIC INTERPOLATION.
1731 :
1732 : !C INPUT ARGUMENTS: (* MEANS PARAMETERS ARE CHANGED DURING EXECUTION)
1733 :
1734 : !C *MODE SEE OUTPUT ARGUMENTS
1735 : !C AX LEFT ENDPOINT OF INITIAL INTERVAL
1736 : !C BX RIGHT ENDPOINT OF INITIAL INTERVAL
1737 : !C F FUNCTION VALUE AT LINMIN WHICH IS TO BE BROUGHT IN BY
1738 : !C REVERSE COMMUNICATION CONTROLLED BY MODE
1739 : !C TOL DESIRED LENGTH OF INTERVAL OF UNCERTAINTY OF FINAL RESULT
1740 :
1741 : !C OUTPUT ARGUMENTS:
1742 :
1743 : !C LINMIN ABSCISSA APPROXIMATING THE POINT WHERE F ATTAINS A MINIMUM
1744 : !C MODE CONTROLS REVERSE COMMUNICATION
1745 : !C MUST BE SET TO 0 INITIALLY, RETURNS WITH INTERMEDIATE
1746 : !C VALUES 1 AND 2 WHICH MUST NOT BE CHANGED BY THE USER,
1747 : !C ENDS WITH CONVERGENCE WITH VALUE 3.
1748 :
1749 : !C WORKING ARRAY:
1750 :
1751 : !C NONE
1752 :
1753 : !C METHOD:
1754 :
1755 : !C THIS FUNCTION SUBPROGRAM IS A SLIGHTLY MODIFIED VERSION OF THE
1756 : !C ALGOL 60 PROCEDURE LOCALMIN GIVEN IN
1757 : !C R.P. BRENT: ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES,
1758 : !C PRENTICE-HALL (1973).
1759 :
1760 : !C IMPLEMENTED BY:
1761 :
1762 : !C KRAFT, D., DFVLR - INSTITUT FUER DYNAMIK DER FLUGSYSTEME
1763 : !C D-8031 OBERPFAFFENHOFEN
1764 :
1765 : !C STATUS: 31. AUGUST 1984
1766 :
1767 : !C SUBROUTINES REQUIRED: NONE
1768 :
1769 : INTEGER mode
1770 : DOUBLE PRECISION f, tol, a, b, c, d, e, p, q, r, u, v, w, x, m,&
1771 : & fu, fv, fw, fx, eps, tol1, tol2, ZERO, ax, bx
1772 : DATA c /0.381966011d0/, eps /1.5d-8/, ZERO /0.0d0/
1773 :
1774 : !C EPS = SQUARE - ROOT OF MACHINE PRECISION
1775 : !C C = GOLDEN SECTION RATIO = (3-SQRT(5))/2
1776 :
1777 0 : GOTO (10, 55), mode
1778 :
1779 : !C INITIALIZATION
1780 :
1781 0 : a = ax
1782 0 : b = bx
1783 0 : e = ZERO
1784 0 : v = a + c*(b - a)
1785 0 : w = v
1786 0 : x = w
1787 0 : linmin = x
1788 0 : mode = 1
1789 0 : GOTO 100
1790 :
1791 : !C MAIN LOOP STARTS HERE
1792 :
1793 0 : 10 fx = f
1794 0 : fv = fx
1795 0 : fw = fv
1796 0 : 20 m = 0.5d0*(a + b)
1797 0 : tol1 = eps*ABS(x) + tol
1798 0 : tol2 = tol1 + tol1
1799 :
1800 : !C TEST CONVERGENCE
1801 :
1802 0 : IF (ABS(x - m) .LE. tol2 - 0.5d0*(b - a)) GOTO 90
1803 0 : r = ZERO
1804 0 : q = r
1805 0 : p = q
1806 0 : IF (ABS(e) .LE. tol1) GOTO 30
1807 :
1808 : !C FIT PARABOLA
1809 :
1810 0 : r = (x - w)*(fx - fv)
1811 0 : q = (x - v)*(fx - fw)
1812 0 : p = (x - v)*q - (x - w)*r
1813 0 : q = q - r
1814 0 : q = q + q
1815 0 : IF (q .GT. ZERO) p = -p
1816 0 : IF (q .LT. ZERO) q = -q
1817 : r = e
1818 0 : e = d
1819 :
1820 : !C IS PARABOLA ACCEPTABLE
1821 :
1822 : 30 IF (ABS(p) .GE. 0.5d0*ABS(q*r) .OR.&
1823 0 : & p .LE. q*(a - x) .OR. p .GE. q*(b-x)) GOTO 40
1824 :
1825 : !C PARABOLIC INTERPOLATION STEP
1826 :
1827 0 : d = p/q
1828 :
1829 : !C F MUST NOT BE EVALUATED TOO CLOSE TO A OR B
1830 :
1831 0 : IF (u - a .LT. tol2) d = SIGN(tol1, m - x)
1832 0 : IF (b - u .LT. tol2) d = SIGN(tol1, m - x)
1833 0 : GOTO 50
1834 :
1835 : !C GOLDEN SECTION STEP
1836 :
1837 0 : 40 IF (x .GE. m) e = a - x
1838 0 : IF (x .LT. m) e = b - x
1839 0 : d = c*e
1840 :
1841 : !C F MUST NOT BE EVALUATED TOO CLOSE TO X
1842 :
1843 0 : 50 IF (ABS(d) .LT. tol1) d = SIGN(tol1, d)
1844 0 : u = x + d
1845 0 : linmin = u
1846 0 : mode = 2
1847 0 : GOTO 100
1848 0 : 55 fu = f
1849 :
1850 : !C UPDATE A, B, V, W, AND X
1851 :
1852 0 : IF (fu .GT. fx) GOTO 60
1853 : IF (u .GE. x) a = x
1854 : IF (u .LT. x) b = x
1855 : v = w
1856 : fv = fw
1857 : w = x
1858 : fw = fx
1859 : x = u
1860 : fx = fu
1861 0 : GOTO 85
1862 : 60 IF (u .LT. x) a = u
1863 : IF (u .GE. x) b = u
1864 0 : IF (fu .LE. fw .OR. w .EQ. x) GOTO 70
1865 0 : IF (fu .LE. fv .OR. v .EQ. x .OR. v .EQ. w) GOTO 80
1866 0 : GOTO 85
1867 : 70 v = w
1868 : fv = fw
1869 : w = u
1870 : fw = fu
1871 0 : GOTO 85
1872 0 : 80 v = u
1873 0 : fv = fu
1874 0 : 85 GOTO 20
1875 :
1876 : !C END OF MAIN LOOP
1877 :
1878 0 : 90 linmin = x
1879 0 : mode = 3
1880 : 100 RETURN
1881 :
1882 : !C END OF LINMIN
1883 :
1884 : END FUNCTION linmin
1885 :
1886 : !C## Following a selection from BLAS Level 1
1887 :
1888 0 : SUBROUTINE daxpy_sl(n,da,dx,incx,dy,incy)
1889 :
1890 : !C CONSTANT TIMES A VECTOR PLUS A VECTOR.
1891 : !C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
1892 : !C JACK DONGARRA, LINPACK, 3/11/78.
1893 :
1894 : DOUBLE PRECISION dx(*),dy(*),da
1895 : INTEGER i,incx,incy,ix,iy,m,mp1,n
1896 :
1897 0 : IF(n.LE.0)RETURN
1898 0 : IF(da.EQ.0.0d0)RETURN
1899 0 : IF(incx.EQ.1.AND.incy.EQ.1)GO TO 20
1900 :
1901 : !C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
1902 : !C NOT EQUAL TO 1
1903 :
1904 0 : ix = 1
1905 0 : iy = 1
1906 0 : IF(incx.LT.0)ix = (-n+1)*incx + 1
1907 0 : IF(incy.LT.0)iy = (-n+1)*incy + 1
1908 0 : DO 10 i = 1,n
1909 0 : dy(iy) = dy(iy) + da*dx(ix)
1910 0 : ix = ix + incx
1911 0 : iy = iy + incy
1912 0 : 10 CONTINUE
1913 0 : RETURN
1914 :
1915 : !C CODE FOR BOTH INCREMENTS EQUAL TO 1
1916 :
1917 : !C CLEAN-UP LOOP
1918 :
1919 0 : 20 m = MOD(n,4)
1920 0 : IF( m .EQ. 0 ) GO TO 40
1921 0 : DO 30 i = 1,m
1922 0 : dy(i) = dy(i) + da*dx(i)
1923 0 : 30 CONTINUE
1924 0 : IF( n .LT. 4 ) RETURN
1925 0 : 40 mp1 = m + 1
1926 0 : DO 50 i = mp1,n,4
1927 0 : dy(i) = dy(i) + da*dx(i)
1928 0 : dy(i + 1) = dy(i + 1) + da*dx(i + 1)
1929 0 : dy(i + 2) = dy(i + 2) + da*dx(i + 2)
1930 0 : dy(i + 3) = dy(i + 3) + da*dx(i + 3)
1931 0 : 50 CONTINUE
1932 : RETURN
1933 : END SUBROUTINE daxpy_sl
1934 :
1935 0 : SUBROUTINE dcopy_(n,dx,incx,dy,incy)
1936 :
1937 : !C COPIES A VECTOR, X, TO A VECTOR, Y.
1938 : !C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
1939 : !C JACK DONGARRA, LINPACK, 3/11/78.
1940 :
1941 : DOUBLE PRECISION dx(*),dy(*)
1942 : INTEGER i,incx,incy,ix,iy,m,mp1,n
1943 :
1944 0 : IF(n.LE.0)RETURN
1945 0 : IF(incx.EQ.1.AND.incy.EQ.1)GO TO 20
1946 :
1947 : !C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
1948 : !C NOT EQUAL TO 1
1949 :
1950 0 : ix = 1
1951 0 : iy = 1
1952 0 : IF(incx.LT.0)ix = (-n+1)*incx + 1
1953 0 : IF(incy.LT.0)iy = (-n+1)*incy + 1
1954 0 : DO 10 i = 1,n
1955 0 : dy(iy) = dx(ix)
1956 0 : ix = ix + incx
1957 0 : iy = iy + incy
1958 0 : 10 CONTINUE
1959 0 : RETURN
1960 :
1961 : !C CODE FOR BOTH INCREMENTS EQUAL TO 1
1962 :
1963 : !C CLEAN-UP LOOP
1964 :
1965 0 : 20 m = MOD(n,7)
1966 0 : IF( m .EQ. 0 ) GO TO 40
1967 0 : DO 30 i = 1,m
1968 0 : dy(i) = dx(i)
1969 0 : 30 CONTINUE
1970 0 : IF( n .LT. 7 ) RETURN
1971 0 : 40 mp1 = m + 1
1972 0 : DO 50 i = mp1,n,7
1973 0 : dy(i) = dx(i)
1974 0 : dy(i + 1) = dx(i + 1)
1975 0 : dy(i + 2) = dx(i + 2)
1976 0 : dy(i + 3) = dx(i + 3)
1977 0 : dy(i + 4) = dx(i + 4)
1978 0 : dy(i + 5) = dx(i + 5)
1979 0 : dy(i + 6) = dx(i + 6)
1980 0 : 50 CONTINUE
1981 : RETURN
1982 : END SUBROUTINE dcopy_
1983 :
1984 0 : DOUBLE PRECISION FUNCTION ddot_sl(n,dx,incx,dy,incy)
1985 :
1986 : !C FORMS THE DOT PRODUCT OF TWO VECTORS.
1987 : !C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
1988 : !C JACK DONGARRA, LINPACK, 3/11/78.
1989 :
1990 : DOUBLE PRECISION dx(*),dy(*),dtemp
1991 : INTEGER i,incx,incy,ix,iy,m,mp1,n
1992 :
1993 0 : ddot_sl = 0.0d0
1994 0 : dtemp = 0.0d0
1995 0 : IF(n.LE.0)RETURN
1996 0 : IF(incx.EQ.1.AND.incy.EQ.1)GO TO 20
1997 :
1998 : !C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
1999 : !C NOT EQUAL TO 1
2000 :
2001 0 : ix = 1
2002 0 : iy = 1
2003 0 : IF(incx.LT.0)ix = (-n+1)*incx + 1
2004 0 : IF(incy.LT.0)iy = (-n+1)*incy + 1
2005 0 : DO 10 i = 1,n
2006 0 : dtemp = dtemp + dx(ix)*dy(iy)
2007 0 : ix = ix + incx
2008 0 : iy = iy + incy
2009 0 : 10 CONTINUE
2010 0 : ddot_sl = dtemp
2011 0 : RETURN
2012 :
2013 : !C CODE FOR BOTH INCREMENTS EQUAL TO 1
2014 :
2015 : !C CLEAN-UP LOOP
2016 :
2017 0 : 20 m = MOD(n,5)
2018 0 : IF( m .EQ. 0 ) GO TO 40
2019 0 : DO 30 i = 1,m
2020 0 : dtemp = dtemp + dx(i)*dy(i)
2021 0 : 30 CONTINUE
2022 0 : IF( n .LT. 5 ) GO TO 60
2023 0 : 40 mp1 = m + 1
2024 0 : DO 50 i = mp1,n,5
2025 : dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + &
2026 0 : dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
2027 0 : 50 CONTINUE
2028 : 60 ddot_sl = dtemp
2029 : RETURN
2030 : END FUNCTION ddot_sl
2031 :
2032 : DOUBLE PRECISION FUNCTION dnrm1(n,x,i,j)
2033 : INTEGER n, i, j, k
2034 : DOUBLE PRECISION snormx, sum, x(n), ZERO, one, scale, temp
2035 : DATA ZERO/0.0d0/, one/1.0d0/
2036 :
2037 : !C DNRM1 - COMPUTES THE I-NORM OF A VECTOR
2038 : !C BETWEEN THE ITH AND THE JTH ELEMENTS
2039 :
2040 : !C INPUT -
2041 : !C N LENGTH OF VECTOR
2042 : !C X VECTOR OF LENGTH N
2043 : !C I INITIAL ELEMENT OF VECTOR TO BE USED
2044 : !C J FINAL ELEMENT TO USE
2045 :
2046 : !C OUTPUT -
2047 : !C DNRM1 NORM
2048 :
2049 : snormx=ZERO
2050 : DO 10 k=i,j
2051 : snormx=MAX(snormx,ABS(x(k)))
2052 : 10 CONTINUE
2053 : dnrm1 = snormx
2054 : IF (snormx.EQ.ZERO) RETURN
2055 : scale = snormx
2056 : IF (snormx.GE.one) scale=SQRT(snormx)
2057 : sum=ZERO
2058 : DO 20 k=i,j
2059 : temp=ZERO
2060 : IF (ABS(x(k))+scale .NE. scale) temp = x(k)/snormx
2061 : IF (one+temp.NE.one) sum = sum+temp*temp
2062 : 20 CONTINUE
2063 : sum=SQRT(sum)
2064 : dnrm1=snormx*sum
2065 : RETURN
2066 : END FUNCTION dnrm1
2067 :
2068 0 : DOUBLE PRECISION FUNCTION dnrm2_ ( n, dx, incx)
2069 : INTEGER n, i, j, nn, next, incx
2070 : DOUBLE PRECISION dx(*), cutlo, cuthi, hitest, sum, xmax, ZERO, one
2071 : DATA ZERO, one /0.0d0, 1.0d0/
2072 :
2073 : !C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE
2074 : !C INCREMENT INCX .
2075 : !C IF N .LE. 0 RETURN WITH RESULT = 0.
2076 : !C IF N .GE. 1 THEN INCX MUST BE .GE. 1
2077 :
2078 : !C C.L.LAWSON, 1978 JAN 08
2079 :
2080 : !C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE
2081 : !C HOPEFULLY APPLICABLE TO ALL MACHINES.
2082 : !C CUTLO = MAXIMUM OF SQRT(U/EPS) OVER ALL KNOWN MACHINES.
2083 : !C CUTHI = MINIMUM OF SQRT(V) OVER ALL KNOWN MACHINES.
2084 : !C WHERE
2085 : !C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
2086 : !C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT)
2087 : !C V = LARGEST NO. (OVERFLOW LIMIT)
2088 :
2089 : !C BRIEF OUTLINE OF ALGORITHM..
2090 :
2091 : !C PHASE 1 SCANS ZERO COMPONENTS.
2092 : !C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
2093 : !C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
2094 : !C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
2095 : !C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
2096 :
2097 : !C VALUES FOR CUTLO AND CUTHI..
2098 : !C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
2099 : !C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
2100 : !C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE
2101 : !C UNIVAC AND DEC AT 2**(-103)
2102 : !C THUS CUTLO = 2**(-51) = 4.44089E-16
2103 : !C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
2104 : !C THUS CUTHI = 2**(63.5) = 1.30438E19
2105 : !C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
2106 : !C THUS CUTLO = 2**(-33.5) = 8.23181D-11
2107 : !C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19
2108 : !C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 /
2109 : !C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 /
2110 : DATA cutlo, cuthi / 8.232d-11, 1.304d19 /
2111 :
2112 0 : IF(n .GT. 0) GO TO 10
2113 0 : dnrm2_ = ZERO
2114 0 : GO TO 300
2115 :
2116 0 : 10 next = 30
2117 0 : sum = ZERO
2118 0 : nn = n * incx
2119 : !C BEGIN MAIN LOOP
2120 0 : i = 1
2121 0 : 20 IF( next .EQ. 30) GO TO 30
2122 0 : IF( next .EQ. 50) GO TO 50
2123 0 : IF( next .EQ. 70) GO TO 70
2124 : IF( next .EQ. 110) GO TO 110
2125 0 : 30 IF( ABS(dx(i)) .GT. cutlo) GO TO 85
2126 : next = 50
2127 0 : xmax = ZERO
2128 :
2129 : !C PHASE 1. SUM IS ZERO
2130 :
2131 0 : 50 IF( dx(i) .EQ. ZERO) GO TO 200
2132 0 : IF( ABS(dx(i)) .GT. cutlo) GO TO 85
2133 :
2134 : !C PREPARE FOR PHASE 2.
2135 :
2136 : next = 70
2137 0 : GO TO 105
2138 :
2139 : !C PREPARE FOR PHASE 4.
2140 :
2141 0 : 100 i = j
2142 0 : next = 110
2143 0 : sum = (sum / dx(i)) / dx(i)
2144 0 : 105 xmax = ABS(dx(i))
2145 0 : GO TO 115
2146 :
2147 : !C PHASE 2. SUM IS SMALL.
2148 : !C SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
2149 :
2150 0 : 70 IF( ABS(dx(i)) .GT. cutlo ) GO TO 75
2151 :
2152 : !C COMMON CODE FOR PHASES 2 AND 4.
2153 : !C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW.
2154 :
2155 0 : 110 IF( ABS(dx(i)) .LE. xmax ) GO TO 115
2156 0 : sum = one + sum * (xmax / dx(i))**2
2157 0 : xmax = ABS(dx(i))
2158 0 : GO TO 200
2159 :
2160 0 : 115 sum = sum + (dx(i)/xmax)**2
2161 0 : GO TO 200
2162 :
2163 : !C PREPARE FOR PHASE 3.
2164 :
2165 0 : 75 sum = (sum * xmax) * xmax
2166 :
2167 : !C FOR REAL OR D.P. SET HITEST = CUTHI/N
2168 : !C FOR COMPLEX SET HITEST = CUTHI/(2*N)
2169 :
2170 0 : 85 hitest = cuthi/float( n )
2171 :
2172 : !C PHASE 3. SUM IS MID-RANGE. NO SCALING.
2173 :
2174 0 : DO 95 j =i,nn,incx
2175 0 : IF(ABS(dx(j)) .GE. hitest) GO TO 100
2176 0 : sum = sum + dx(j)**2
2177 0 : 95 CONTINUE
2178 0 : dnrm2_ = SQRT( sum )
2179 0 : GO TO 300
2180 :
2181 : 200 CONTINUE
2182 0 : i = i + incx
2183 0 : IF ( i .LE. nn ) GO TO 20
2184 :
2185 : !C END OF MAIN LOOP.
2186 :
2187 : !C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
2188 :
2189 0 : dnrm2_ = xmax * SQRT(sum)
2190 : 300 CONTINUE
2191 : RETURN
2192 : END FUNCTION dnrm2_
2193 :
2194 0 : SUBROUTINE dsrot (n,dx,incx,dy,incy,c,s)
2195 :
2196 : !C APPLIES A PLANE ROTATION.
2197 : !C JACK DONGARRA, LINPACK, 3/11/78.
2198 :
2199 : DOUBLE PRECISION dx(*),dy(*),dtemp,c,s
2200 : INTEGER i,incx,incy,ix,iy,n
2201 :
2202 0 : IF(n.LE.0)RETURN
2203 0 : IF(incx.EQ.1.AND.incy.EQ.1)GO TO 20
2204 :
2205 : !C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
2206 : !C TO 1
2207 :
2208 0 : ix = 1
2209 0 : iy = 1
2210 0 : IF(incx.LT.0)ix = (-n+1)*incx + 1
2211 0 : IF(incy.LT.0)iy = (-n+1)*incy + 1
2212 0 : DO 10 i = 1,n
2213 0 : dtemp = c*dx(ix) + s*dy(iy)
2214 0 : dy(iy) = c*dy(iy) - s*dx(ix)
2215 0 : dx(ix) = dtemp
2216 0 : ix = ix + incx
2217 0 : iy = iy + incy
2218 0 : 10 CONTINUE
2219 0 : RETURN
2220 :
2221 : !C CODE FOR BOTH INCREMENTS EQUAL TO 1
2222 :
2223 0 : 20 DO 30 i = 1,n
2224 0 : dtemp = c*dx(i) + s*dy(i)
2225 0 : dy(i) = c*dy(i) - s*dx(i)
2226 0 : dx(i) = dtemp
2227 0 : 30 CONTINUE
2228 : RETURN
2229 : END SUBROUTINE dsrot
2230 :
2231 0 : SUBROUTINE dsrotg(da,db,c,s)
2232 :
2233 : !C CONSTRUCT GIVENS PLANE ROTATION.
2234 : !C JACK DONGARRA, LINPACK, 3/11/78.
2235 : !C MODIFIED 9/27/86.
2236 :
2237 : DOUBLE PRECISION da,db,c,s,roe,scale,r,z,one,ZERO
2238 : DATA one, ZERO /1.0d+00, 0.0d+00/
2239 :
2240 0 : roe = db
2241 0 : IF( ABS(da) .GT. ABS(db) ) roe = da
2242 0 : scale = ABS(da) + ABS(db)
2243 0 : IF( scale .NE. ZERO ) GO TO 10
2244 0 : c = one
2245 0 : s = ZERO
2246 0 : r = ZERO
2247 0 : GO TO 20
2248 0 : 10 r = scale*SQRT((da/scale)**2 + (db/scale)**2)
2249 0 : r = SIGN(one,roe)*r
2250 0 : c = da/r
2251 0 : s = db/r
2252 0 : 20 z = s
2253 0 : IF( ABS(c) .GT. ZERO .AND. ABS(c) .LE. s ) z = one/c
2254 0 : da = r
2255 0 : db = z
2256 0 : RETURN
2257 : END SUBROUTINE dsrotg
2258 :
2259 0 : SUBROUTINE dscal_sl(n,da,dx,incx)
2260 :
2261 : !C SCALES A VECTOR BY A CONSTANT.
2262 : !C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
2263 : !C JACK DONGARRA, LINPACK, 3/11/78.
2264 :
2265 : DOUBLE PRECISION da,dx(*)
2266 : INTEGER i,incx,m,mp1,n,nincx
2267 :
2268 0 : IF(n.LE.0)RETURN
2269 0 : IF(incx.EQ.1)GO TO 20
2270 :
2271 :
2272 : !C CODE FOR INCREMENT NOT EQUAL TO 1
2273 :
2274 0 : nincx = n*incx
2275 0 : DO 10 i = 1,nincx,incx
2276 0 : dx(i) = da*dx(i)
2277 0 : 10 CONTINUE
2278 0 : RETURN
2279 :
2280 : !C CODE FOR INCREMENT EQUAL TO 1
2281 :
2282 : !C CLEAN-UP LOOP
2283 :
2284 0 : 20 m = MOD(n,5)
2285 0 : IF( m .EQ. 0 ) GO TO 40
2286 0 : DO 30 i = 1,m
2287 0 : dx(i) = da*dx(i)
2288 0 : 30 CONTINUE
2289 0 : IF( n .LT. 5 ) RETURN
2290 0 : 40 mp1 = m + 1
2291 0 : DO 50 i = mp1,n,5
2292 0 : dx(i) = da*dx(i)
2293 0 : dx(i + 1) = da*dx(i + 1)
2294 0 : dx(i + 2) = da*dx(i + 2)
2295 0 : dx(i + 3) = da*dx(i + 3)
2296 0 : dx(i + 4) = da*dx(i + 4)
2297 0 : 50 CONTINUE
2298 : RETURN
2299 : END SUBROUTINE dscal_sl
2300 :
2301 0 : subroutine bound(n, x, xl, xu)
2302 : integer n, i
2303 : double precision x(n), xl(n), xu(n)
2304 0 : do i = 1, n
2305 : !C Note that xl(i) and xu(i) may be NaN to indicate no bound
2306 0 : if(xl(i).eq.xl(i).and.x(i) < xl(i))then
2307 0 : x(i) = xl(i)
2308 0 : else if(xu(i).eq.xu(i).and.x(i) > xu(i))then
2309 0 : x(i) = xu(i)
2310 : end if
2311 : end do
2312 0 : end subroutine bound
2313 :
2314 : !----------------------------------------------------------------------
2315 :
2316 : END MODULE m_slsqp
2317 : !!***
2318 :
|
