Branch data 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 : :
|
