Mantid
Loading...
Searching...
No Matches
SLSQPMinimizer.cpp
Go to the documentation of this file.
1#include "MantidKernel/Math/Optimization/SLSQPMinimizer.h"
2#include "MantidKernel/Exception.h"
3
4#include <algorithm>
5#include <cassert>
6#include <cmath>
7#include <sstream>
8
9namespace Mantid::Kernel::Math {
10namespace {
12// Forward-declaration of minimizer
13int slsqp_(int *m, int *meq, int *la, int *n, double *x, double *xl, double *xu, double *f, double *c__, double *g,
14 double *a, double *acc, int *iter, int *mode, double *w, const int *l_w__, int *jw, const int *l_jw__);
16} // namespace
17
23std::vector<double> SLSQPMinimizer::minimize(const std::vector<double> &x0) const {
24 assert(numParameters() == x0.size());
25
26 // slsqp parameters
27 auto m = static_cast<int>(numEqualityConstraints() + numInequalityConstraints());
28 auto meq = static_cast<int>(numEqualityConstraints());
29 int la = std::max(1, m);
30 auto n = static_cast<int>(numParameters());
31
32 std::vector<double> x = x0;
33 std::vector<double> xl(n, -1e12), xu(n, 1e12); // bounds
34 // Function value holder
35 double fx(0.0);
36 // Constraint value holder
37 std::vector<double> constraintValues(m, 0.0);
38
39 // slsqp requires gradient array to be one longer than number of parameters
40 std::vector<double> gradient(m_nparams + 1, 0.0);
41
42 double acc = 1e-06; // acceptance
43 int majiter = 100;
44 int mode = 0; // Starts in mode 0
45
46 // workspace spaces for slsqp
47 int n1 = n + 1;
48 int mineq = m - meq + n1 + n1;
49 int len_w = (3 * n1 + m) * (n1 + 1) + (n1 - meq + 1) * (mineq + 2) + 2 * mineq + (n1 + mineq) * (n1 - meq) + 2 * meq +
50 n1 + (n + 1) * n / 2 + 2 * m + 3 * n + 3 * n1 + 1;
51 int len_jw = mineq;
52 std::vector<double> vec_w(len_w, 0.0);
53 double *w = vec_w.data();
54 std::vector<int> vec_jw(len_jw, 0);
55 int *jw = vec_jw.data();
56
57 while (true) {
58 if (mode == 0 || mode == 1) // objective and constraint evaluation required
59 {
60 // Compute objective function
61 fx = fvalue(x);
62 evaluateConstraints(constraintValues, x);
63 }
64 if (mode == 0 || mode == -1) // gradient evaluation required
65 {
66 // Compute the derivatives of the objective function
67 fprime(gradient, x);
68 }
69 // Call SLSQP
70 slsqp_(&m, &meq, &la, &n, x.data(), xl.data(), xu.data(), &fx, constraintValues.data(), gradient.data(),
71 m_constraintNorms.data(), &acc, &majiter, &mode, w, &len_w, jw, &len_jw);
72
73 // If exit mode is not -1 or 1, slsqp has completed
74 if (mode != 1 && mode != -1)
75 break;
76 }
77 return x;
78}
79
85void SLSQPMinimizer::fprime(std::vector<double> &grad, const std::vector<double> &x) const {
86 assert(grad.size() > numParameters()); // > for slsqp
87
88 const double epsilon(1e-08);
89
90 double f0 = fvalue(x);
91 std::vector<double> tmp = x;
92 for (size_t i = 0; i < numParameters(); ++i) {
93 double &xi = tmp[i];
94 const double curx = xi; // copy
95 xi += epsilon;
96 grad[i] = (fvalue(tmp) - f0) / epsilon;
97 xi = curx;
98 }
99}
100
106void SLSQPMinimizer::evaluateConstraints(std::vector<double> &constrValues, const std::vector<double> &x) const {
107 assert(numEqualityConstraints() + numInequalityConstraints() == constrValues.size());
108
109 // Need to compute position in flat array
110 const size_t numConstrs(constrValues.size());
111 for (size_t i = 0; i < numConstrs; ++i) {
112 double value(0.0);
113 for (size_t j = 0; j < numParameters(); ++j) {
114 value += m_constraintNorms[j * numConstrs + i] * x[j];
115 }
116 constrValues[i] = value;
117 }
118}
119
129void SLSQPMinimizer::initializeConstraints(const DblMatrix &equality, const DblMatrix &inequality) {
130 const size_t totalNumConstr = numEqualityConstraints() + numInequalityConstraints();
131 if (totalNumConstr == 0)
132 return;
133
134 // Sanity checks on matrix sizes
135 for (size_t i = 0; i < 2; ++i) {
136 size_t ncols(0);
137 std::string matrix;
138 if (i == 0) {
139 ncols = equality.numCols();
140 matrix = "equality";
141 } else {
142 ncols = inequality.numCols();
143 matrix = "inequality";
144 }
145
146 if (ncols > 0 && ncols != numParameters()) {
147 std::ostringstream os;
148 os << "SLSQPMinimizer::initializeConstraints - Invalid " << matrix
149 << " constraint matrix. Number of columns must match number of "
150 "parameters. ncols="
151 << ncols << ", nparams=" << numParameters();
152 throw std::invalid_argument(os.str());
153 }
154 }
155
156 // --- Constraint Vector Layout As Required By SLSQP -----
157 // The equality constraints are specified first, followed by the inequality
158 // constraints. The constraints
159 // should be written column-wise such that to get move between successive
160 // values of the same constraint you
161 // jump (numEqualityConstraints()+numInequalityConstraints()) in the array.
162 // The array is then padded
163 // by (numEqualityConstraints()+numInEqualityConstraints()) zeroes
164 const size_t constrVecSize = totalNumConstr * numParameters() + totalNumConstr;
165 m_constraintNorms.resize(constrVecSize, 0.0);
166
167 size_t constrCounter(0);
168 while (constrCounter < totalNumConstr) {
169 const DblMatrix *constrMatrix(nullptr);
170 if (constrCounter < numEqualityConstraints())
171 constrMatrix = &equality;
172 else
173 constrMatrix = &inequality;
174
175 for (size_t i = 0; i < constrMatrix->numRows(); ++i, ++constrCounter) {
176 const auto matrixRow = (*constrMatrix)[i];
177 for (size_t j = 0; j < constrMatrix->numCols(); ++j) {
178 m_constraintNorms[j * totalNumConstr + constrCounter] = matrixRow[j];
179 }
180 }
181 }
182 assert(totalNumConstr == constrCounter);
183}
184
185namespace {
187
188//---------------------------------------------------------------------------------------------------//
189//
190// The following code was translated from fortran using f2c and then
191// hand-cleaned so that it will
192// compile without the f2c headers.
193//
194// The main entry point is the slsqp_ function. All others are helpers.
195//
196// By hand changes:
197// - doublereal --> double
198// - integer --> int
199// - real --> float
200//
201// Line 2639: Commented out next_fmt declaration and corresponding fmt_
202// statements that cause compiler
203// warnings and are unused
204//---------------------------------------------------------------------------------------------------//
205
206/* Table of constant values */
207static int c__0 = 0;
208static int c__1 = 1;
209static int c__2 = 2;
210
211double d_sign(const double *a, const double *b) {
212 double x;
213 x = (*a >= 0 ? *a : -*a);
214 return (*b >= 0 ? x : -x);
215}
216
217// --------------------- Forward declarations of helpers ---------------------
218int slsqpb_(int * /*m*/, int * /*meq*/, int * /*la*/, int * /*n*/, double * /*x*/, double * /*xl*/, double * /*xu*/,
219 const double * /*f*/, double * /*c__*/, double * /*g*/, double * /*a*/, double * /*acc*/, int * /*iter*/,
220 int * /*mode*/, double * /*r__*/, double * /*l*/, double * /*x0*/, double * /*mu*/, double * /*s*/,
221 double * /*u*/, double * /*v*/, double * /*w*/, int * /*iw*/);
222int dcopy___(const int *n, double *dx, const int *incx, double *dy, const int *incy);
223int daxpy_sl__(const int *n, const double *da, double *dx, const int *incx, double *dy, const int *incy);
224int lsq_(int *m, int *meq, int *n, const int *nl, int *la, double *l, double *g, double *a, double *b, double *xl,
225 double *xu, double *x, double *y, double *w, int *jw, int *mode);
226double ddot_sl__(const int *n, double *dx, const int *incx, double *dy, const int *incy);
227int dscal_sl__(const int *n, const double *da, double *dx, const int *incx);
228double linmin_(int *mode, const double *ax, const double *bx, const double *f, const double *tol);
229double dnrm2___(const int *n, double *dx, const int *incx);
230int ldl_(const int *n, double *a, double *z__, const double *sigma, double *w);
231int lsei_(double *c__, double *d__, double *e, double *f, double *g, double *h__, int *lc, int *mc, int *le, int *me,
232 int *lg, int *mg, int *n, double *x, double *xnrm, double *w, int *jw, int *mode);
233int h12_(const int *mode, const int *lpivot, const int *l1, const int *m, double *u, const int *iue, double *up,
234 double *c__, const int *ice, const int *icv, const int *ncv);
235int hfti_(double *a, int *mda, int *m, int *n, double *b, int *mdb, int *nb, const double *tau, int *krank,
236 double *rnorm, double *h__, double *g, int *ip);
237int lsi_(double *e, double *f, double *g, double *h__, int *le, int *me, int *lg, int *mg, int *n, double *x,
238 double *xnorm, double *w, int *jw, int *mode);
239int ldp_(double *g, const int *mg, int *m, int *n, double *h__, double *x, double *xnorm, double *w, int *index,
240 int *mode);
241int nnls_(double *a, int *mda, int *m, int *n, double *b, double *x, double *rnorm, double *w, double *z__, int *index,
242 int *mode);
243int dsrotg_(double *da, double *db, double *c__, double *s);
244int dsrot_(const int *n, double *dx, const int *incx, double *dy, const int *incy, const double *c__, const double *s);
245//----------------------------------------------------------------------------
246
247/* ALGORITHM 733, COLLECTED ALGORITHMS FROM ACM. */
248/* TRANSACTIONS ON MATHEMATICAL SOFTWARE, */
249/* VOL. 20, NO. 3, SEPTEMBER, 1994, PP. 262-281. */
250/* http://doi.acm.org/10.1145/192115.192124 */
251
252/* http://permalink.gmane.org/gmane.comp.python.scientific.devel/6725 */
253/* ------ */
254/* From: Deborah Cotton <cotton@hq.acm.org> */
255/* Date: Fri, 14 Sep 2007 12:35:55 -0500 */
256/* Subject: RE: Algorithm License requested */
257/* To: Alan Isaac */
258
259/* Prof. Issac, */
260
261/* In that case, then because the author consents to [the ACM] releasing */
262/* the code currently archived at http://www.netlib.org/toms/733 under the
263 */
264/* BSD license, the ACM hereby releases this code under the BSD license. */
265
266/* Regards, */
267
268/* Deborah Cotton, Copyright & Permissions */
269/* ACM Publications */
270/* 2 Penn Plaza, Suite 701** */
271/* New York, NY 10121-0701 */
272/* permissions@acm.org */
273/* 212.869.7440 ext. 652 */
274/* Fax. 212.869.0481 */
275/* ------ */
276
277/* *********************************************************************** */
278/* optimizer * */
279/* *********************************************************************** */
280int slsqp_(int *m, int *meq, int *la, int *n, double *x, double *xl, double *xu, double *f, double *c__, double *g,
281 double *a, double *acc, int *iter, int *mode, double *w, const int *l_w__, int *jw, const int *l_jw__) {
282 /* System generated locals */
283 int a_dim1, a_offset, i__1, i__2;
284
285 /* Local variables */
286 static int n1, il, im, ir, is, iu, iv, iw, ix, mineq;
287
288 /* SLSQP S EQUENTIAL L EAST SQ UARES P ROGRAMMING */
289 /* TO SOLVE GENERAL NONLINEAR OPTIMIZATION PROBLEMS */
290 /* *********************************************************************** */
291 /* * * */
292 /* * * */
293 /* * A NONLINEAR PROGRAMMING METHOD WITH * */
294 /* * QUADRATIC PROGRAMMING SUBPROBLEMS * */
295 /* * * */
296 /* * * */
297 /* * THIS SUBROUTINE SOLVES THE GENERAL NONLINEAR PROGRAMMING PROBLEM * */
298 /* * * */
299 /* * MINIMIZE F(X) * */
300 /* * * */
301 /* * SUBJECT TO C (X) .EQ. 0 , J = 1,...,MEQ * */
302 /* * J * */
303 /* * * */
304 /* * C (X) .GE. 0 , J = MEQ+1,...,M * */
305 /* * J * */
306 /* * * */
307 /* * XL .LE. X .LE. XU , I = 1,...,N. * */
308 /* * I I I * */
309 /* * * */
310 /* * THE ALGORITHM IMPLEMENTS THE METHOD OF HAN AND POWELL * */
311 /* * WITH BFGS-UPDATE OF THE B-MATRIX AND L1-TEST FUNCTION * */
312 /* * WITHIN THE STEPLENGTH ALGORITHM. * */
313 /* * * */
314 /* * PARAMETER DESCRIPTION: * */
315 /* * ( * MEANS THIS PARAMETER WILL BE CHANGED DURING CALCULATION ) * */
316 /* * * */
317 /* * M IS THE TOTAL NUMBER OF CONSTRAINTS, M .GE. 0 * */
318 /* * MEQ IS THE NUMBER OF EQUALITY CONSTRAINTS, MEQ .GE. 0 * */
319 /* * LA SEE A, LA .GE. MAX(M,1) * */
320 /* * N IS THE NUMBER OF VARIBLES, N .GE. 1 * */
321 /* * * X() X() STORES THE CURRENT ITERATE OF THE N VECTOR X * */
322 /* * ON ENTRY X() MUST BE INITIALIZED. ON EXIT X() * */
323 /* * STORES THE SOLUTION VECTOR X IF MODE = 0. * */
324 /* * XL() XL() STORES AN N VECTOR OF LOWER BOUNDS XL TO X. * */
325 /* * XU() XU() STORES AN N VECTOR OF UPPER BOUNDS XU TO X. * */
326 /* * F IS THE VALUE OF THE OBJECTIVE FUNCTION. * */
327 /* * C() C() STORES THE M VECTOR C OF CONSTRAINTS, * */
328 /* * EQUALITY CONSTRAINTS (IF ANY) FIRST. * */
329 /* * DIMENSION OF C MUST BE GREATER OR EQUAL LA, * */
330 /* * which must be GREATER OR EQUAL MAX(1,M). * */
331 /* * G() G() STORES THE N VECTOR G OF PARTIALS OF THE * */
332 /* * OBJECTIVE FUNCTION; DIMENSION OF G MUST BE * */
333 /* * GREATER OR EQUAL N+1. * */
334 /* * A(),LA,M,N THE LA BY N + 1 ARRAY A() STORES * */
335 /* * THE M BY N MATRIX A OF CONSTRAINT NORMALS. * */
336 /* * A() HAS FIRST DIMENSIONING PARAMETER LA, * */
337 /* * WHICH MUST BE GREATER OR EQUAL MAX(1,M). * */
338 /* * F,C,G,A MUST ALL BE SET BY THE USER BEFORE EACH CALL. * */
339 /* * * ACC ABS(ACC) CONTROLS THE FINAL ACCURACY. * */
340 /* * IF ACC .LT. ZERO AN EXACT LINESEARCH IS PERFORMED,* */
341 /* * OTHERWISE AN ARMIJO-TYPE LINESEARCH IS USED. * */
342 /* * * ITER PRESCRIBES THE MAXIMUM NUMBER OF ITERATIONS. * */
343 /* * ON EXIT ITER INDICATES THE NUMBER OF ITERATIONS. * */
344 /* * * MODE MODE CONTROLS CALCULATION: * */
345 /* * REVERSE COMMUNICATION IS USED IN THE SENSE THAT * */
346 /* * THE PROGRAM IS INITIALIZED BY MODE = 0; THEN IT IS* */
347 /* * TO BE CALLED REPEATEDLY BY THE USER UNTIL A RETURN* */
348 /* * WITH MODE .NE. IABS(1) TAKES PLACE. * */
349 /* * IF MODE = -1 GRADIENTS HAVE TO BE CALCULATED, * */
350 /* * WHILE WITH MODE = 1 FUNCTIONS HAVE TO BE CALCULATED */
351 /* * MODE MUST NOT BE CHANGED BETWEEN SUBSEQUENT CALLS * */
352 /* * OF SQP. * */
353 /* * EVALUATION MODES: * */
354 /* * MODE = -1: GRADIENT EVALUATION, (G&A) * */
355 /* * 0: ON ENTRY: INITIALIZATION, (F,G,C&A) * */
356 /* * ON EXIT : REQUIRED ACCURACY FOR SOLUTION OBTAINED * */
357 /* * 1: FUNCTION EVALUATION, (F&C) * */
358 /* * * */
359 /* * FAILURE MODES: * */
360 /* * 2: NUMBER OF EQUALITY CONTRAINTS LARGER THAN N * */
361 /* * 3: MORE THAN 3*N ITERATIONS IN LSQ SUBPROBLEM * */
362 /* * 4: INEQUALITY CONSTRAINTS INCOMPATIBLE * */
363 /* * 5: SINGULAR MATRIX E IN LSQ SUBPROBLEM * */
364 /* * 6: SINGULAR MATRIX C IN LSQ SUBPROBLEM * */
365 /* * 7: RANK-DEFICIENT EQUALITY CONSTRAINT SUBPROBLEM HFTI* */
366 /* * 8: POSITIVE DIRECTIONAL DERIVATIVE FOR LINESEARCH * */
367 /* * 9: MORE THAN ITER ITERATIONS IN SQP * */
368 /* * >=10: WORKING SPACE W OR JW TOO SMALL, * */
369 /* * W SHOULD BE ENLARGED TO L_W=MODE/1000 * */
370 /* * JW SHOULD BE ENLARGED TO L_JW=MODE-1000*L_W * */
371 /* * * W(), L_W W() IS A ONE DIMENSIONAL WORKING SPACE, * */
372 /* * THE LENGTH L_W OF WHICH SHOULD BE AT LEAST * */
373 /* * (3*N1+M)*(N1+1) for LSQ * */
374 /* * +(N1-MEQ+1)*(MINEQ+2) + 2*MINEQ for LSI * */
375 /* * +(N1+MINEQ)*(N1-MEQ) + 2*MEQ + N1 for LSEI * */
376 /* * + N1*N/2 + 2*M + 3*N + 3*N1 + 1 for SLSQPB * */
377 /* * with MINEQ = M - MEQ + 2*N1 & N1 = N+1 * */
378 /* * NOTICE: FOR PROPER DIMENSIONING OF W IT IS RECOMMENDED TO * */
379 /* * COPY THE FOLLOWING STATEMENTS INTO THE HEAD OF * */
380 /* * THE CALLING PROGRAM (AND REMOVE THE COMMENT C) * */
381 /* ####################################################################### */
382 /* INT LEN_W, LEN_JW, M, N, N1, MEQ, MINEQ */
383 /* PARAMETER (M=... , MEQ=... , N=... ) */
384 /* PARAMETER (N1= N+1, MINEQ= M-MEQ+N1+N1) */
385 /* PARAMETER (LEN_W= */
386 /* $ (3*N1+M)*(N1+1) */
387 /* $ +(N1-MEQ+1)*(MINEQ+2) + 2*MINEQ */
388 /* $ +(N1+MINEQ)*(N1-MEQ) + 2*MEQ + N1 */
389 /* $ +(N+1)*N/2 + 2*M + 3*N + 3*N1 + 1, */
390 /* $ LEN_JW=MINEQ) */
391 /* DOUBLE PRECISION W(LEN_W) */
392 /* INT JW(LEN_JW) */
393 /* ####################################################################### */
394 /* * THE FIRST M+N+N*N1/2 ELEMENTS OF W MUST NOT BE * */
395 /* * CHANGED BETWEEN SUBSEQUENT CALLS OF SLSQP. * */
396 /* * ON RETURN W(1) ... W(M) CONTAIN THE MULTIPLIERS * */
397 /* * ASSOCIATED WITH THE GENERAL CONSTRAINTS, WHILE * */
398 /* * W(M+1) ... W(M+N(N+1)/2) STORE THE CHOLESKY FACTOR* */
399 /* * L*D*L(T) OF THE APPROXIMATE HESSIAN OF THE * */
400 /* * LAGRANGIAN COLUMNWISE DENSE AS LOWER TRIANGULAR * */
401 /* * UNIT MATRIX L WITH D IN ITS 'DIAGONAL' and * */
402 /* * W(M+N(N+1)/2+N+2 ... W(M+N(N+1)/2+N+2+M+2N) * */
403 /* * CONTAIN THE MULTIPLIERS ASSOCIATED WITH ALL * */
404 /* * ALL CONSTRAINTS OF THE QUADRATIC PROGRAM FINDING * */
405 /* * THE SEARCH DIRECTION TO THE SOLUTION X* * */
406 /* * * JW(), L_JW JW() IS A ONE DIMENSIONAL INT WORKING SPACE * */
407 /* * THE LENGTH L_JW OF WHICH SHOULD BE AT LEAST * */
408 /* * MINEQ * */
409 /* * with MINEQ = M - MEQ + 2*N1 & N1 = N+1 * */
410 /* * * */
411 /* * THE USER HAS TO PROVIDE THE FOLLOWING SUBROUTINES: * */
412 /* * LDL(N,A,Z,SIG,W) : UPDATE OF THE LDL'-FACTORIZATION. * */
413 /* * LINMIN(A,B,F,TOL) : LINESEARCH ALGORITHM IF EXACT = 1 * */
414 /* * LSQ(M,MEQ,LA,N,NC,C,D,A,B,XL,XU,X,LAMBDA,W,....) : * */
415 /* * * */
416 /* * SOLUTION OF THE QUADRATIC PROGRAM * */
417 /* * QPSOL IS RECOMMENDED: * */
418 /* * PE GILL, W MURRAY, MA SAUNDERS, MH WRIGHT: * */
419 /* * USER'S GUIDE FOR SOL/QPSOL: * */
420 /* * A FORTRAN PACKAGE FOR QUADRATIC PROGRAMMING, * */
421 /* * TECHNICAL REPORT SOL 83-7, JULY 1983 * */
422 /* * DEPARTMENT OF OPERATIONS RESEARCH, STANFORD UNIVERSITY * */
423 /* * STANFORD, CA 94305 * */
424 /* * QPSOL IS THE MOST ROBUST AND EFFICIENT QP-SOLVER * */
425 /* * AS IT ALLOWS WARM STARTS WITH PROPER WORKING SETS * */
426 /* * * */
427 /* * IF IT IS NOT AVAILABLE USE LSEI, A CONSTRAINT LINEAR LEAST * */
428 /* * SQUARES SOLVER IMPLEMENTED USING THE SOFTWARE HFTI, LDP, NNLS * */
429 /* * FROM C.L. LAWSON, R.J.HANSON: SOLVING LEAST SQUARES PROBLEMS, * */
430 /* * PRENTICE HALL, ENGLEWOOD CLIFFS, 1974. * */
431 /* * LSEI COMES WITH THIS PACKAGE, together with all necessary SR's. * */
432 /* * * */
433 /* * TOGETHER WITH A COUPLE OF SUBROUTINES FROM BLAS LEVEL 1 * */
434 /* * * */
435 /* * SQP IS HEAD SUBROUTINE FOR BODY SUBROUTINE SQPBDY * */
436 /* * IN WHICH THE ALGORITHM HAS BEEN IMPLEMENTED. * */
437 /* * * */
438 /* * IMPLEMENTED BY: DIETER KRAFT, DFVLR OBERPFAFFENHOFEN * */
439 /* * as described in Dieter Kraft: A Software Package for * */
440 /* * Sequential Quadratic Programming * */
441 /* * DFVLR-FB 88-28, 1988 * */
442 /* * which should be referenced if the user publishes results of SLSQP * */
443 /* * * */
444 /* * DATE: APRIL - OCTOBER, 1981. * */
445 /* * STATUS: DECEMBER, 31-ST, 1984. * */
446 /* * STATUS: MARCH , 21-ST, 1987, REVISED TO FORTAN 77 * */
447 /* * STATUS: MARCH , 20-th, 1989, REVISED TO MS-FORTRAN * */
448 /* * STATUS: APRIL , 14-th, 1989, HESSE in-line coded * */
449 /* * STATUS: FEBRUARY, 28-th, 1991, FORTRAN/2 Version 1.04 * */
450 /* * accepts Statement Functions * */
451 /* * STATUS: MARCH , 1-st, 1991, tested with SALFORD * */
452 /* * FTN77/386 COMPILER VERS 2.40* */
453 /* * in protected mode * */
454 /* * * */
455 /* *********************************************************************** */
456 /* * * */
457 /* * Copyright 1991: Dieter Kraft, FHM * */
458 /* * * */
459 /* *********************************************************************** */
460 /* dim(W) = N1*(N1+1) + MEQ*(N1+1) + MINEQ*(N1+1) for LSQ */
461 /* +(N1-MEQ+1)*(MINEQ+2) + 2*MINEQ for LSI */
462 /* +(N1+MINEQ)*(N1-MEQ) + 2*MEQ + N1 for LSEI */
463 /* + N1*N/2 + 2*M + 3*N +3*N1 + 1 for SLSQPB */
464 /* with MINEQ = M - MEQ + 2*N1 & N1 = N+1 */
465 /* CHECK LENGTH OF WORKING ARRAYS */
466 /* Parameter adjustments */
467 --c__;
468 a_dim1 = *la;
469 a_offset = 1 + a_dim1;
470 a -= a_offset;
471 --g;
472 --xu;
473 --xl;
474 --x;
475 --w;
476 --jw;
477
478 /* Function Body */
479 n1 = *n + 1;
480 mineq = *m - *meq + n1 + n1;
481 il = (n1 * 3 + *m) * (n1 + 1) + (n1 - *meq + 1) * (mineq + 2) + (mineq << 1) + (n1 + mineq) * (n1 - *meq) +
482 (*meq << 1) + n1 * *n / 2 + (*m << 1) + *n * 3 + (n1 << 2) + 1;
483 /* Computing MAX */
484 i__1 = mineq, i__2 = n1 - *meq;
485 im = std::max(i__1, i__2);
486 if (*l_w__ < il || *l_jw__ < im) {
487 *mode = std::max(10, il) * 1000;
488 *mode += std::max(10, im);
489 return 0;
490 }
491 /* PREPARE DATA FOR CALLING SQPBDY - INITIAL ADDRESSES IN W */
492 im = 1;
493 il = im + std::max(1, *m);
494 il = im + *la;
495 ix = il + n1 * *n / 2 + 1;
496 ir = ix + *n;
497 is = ir + *n + *n + std::max(1, *m);
498 is = ir + *n + *n + *la;
499 iu = is + n1;
500 iv = iu + n1;
501 iw = iv + n1;
502
503 slsqpb_(m, meq, la, n, &x[1], &xl[1], &xu[1], f, &c__[1], &g[1], &a[a_offset], acc, iter, mode, &w[ir], &w[il],
504 &w[ix], &w[im], &w[is], &w[iu], &w[iv], &w[iw], &jw[1]);
505 return 0;
506} /* slsqp_ */
507
508int slsqpb_(int *m, int *meq, int *la, int *n, double *x, double *xl, double *xu, const double *f, double *c__,
509 double *g, double *a, double *acc, int *iter, int *mode, double *r__, double *l, double *x0, double *mu,
510 double *s, double *u, double *v, double *w, int *iw) {
511 /* Initialized data */
512
513 static double zero = 0.;
514 static double one = 1.;
515 static double alfmin = .1;
516 static double hun = 100.;
517 static double ten = 10.;
518 static double two = 2.;
519
520 /* System generated locals */
521 int a_dim1, a_offset, i__1, i__2;
522 double d__1, d__2;
523
524 /* Local variables */
525 static int i__, j, k;
526 static double t, f0, h1, h2, h3, h4;
527 static int n1, n2, n3;
528 static double t0, gs;
529 static double tol;
530 static int line;
531 static double alpha;
532 static int iexact;
533 static int incons, ireset, itermx;
534
535 /* NONLINEAR PROGRAMMING BY SOLVING SEQUENTIALLY QUADRATIC PROGRAMS */
536 /* - L1 - LINE SEARCH, POSITIVE DEFINITE BFGS UPDATE - */
537 /* BODY SUBROUTINE FOR SLSQP */
538 /* dim(W) = N1*(N1+1) + MEQ*(N1+1) + MINEQ*(N1+1) for LSQ */
539 /* +(N1-MEQ+1)*(MINEQ+2) + 2*MINEQ */
540 /* +(N1+MINEQ)*(N1-MEQ) + 2*MEQ + N1 for LSEI */
541 /* with MINEQ = M - MEQ + 2*N1 & N1 = N+1 */
542 /* Parameter adjustments */
543 --mu;
544 --c__;
545 --v;
546 --u;
547 --s;
548 --x0;
549 --l;
550 --r__;
551 a_dim1 = *la;
552 a_offset = 1 + a_dim1;
553 a -= a_offset;
554 --g;
555 --xu;
556 --xl;
557 --x;
558 --w;
559 --iw;
560
561 /* Function Body */
562 if (*mode < 0) {
563 goto L260;
564 } else if (*mode == 0) {
565 goto L100;
566 } else {
567 goto L220;
568 }
569L100:
570 itermx = *iter;
571 if (*acc >= zero) {
572 iexact = 0;
573 } else {
574 iexact = 1;
575 }
576 *acc = std::abs(*acc);
577 tol = ten * *acc;
578 *iter = 0;
579 ireset = 0;
580 n1 = *n + 1;
581 n2 = n1 * *n / 2;
582 n3 = n2 + 1;
583 s[1] = zero;
584 mu[1] = zero;
585 dcopy___(n, &s[1], &c__0, &s[1], &c__1);
586 dcopy___(m, &mu[1], &c__0, &mu[1], &c__1);
587/* RESET BFGS MATRIX */
588L110:
589 ++ireset;
590 if (ireset > 5) {
591 goto L255;
592 }
593 l[1] = zero;
594 dcopy___(&n2, &l[1], &c__0, &l[1], &c__1);
595 j = 1;
596 i__1 = *n;
597 for (i__ = 1; i__ <= i__1; ++i__) {
598 l[j] = one;
599 j = j + n1 - i__;
600 /* L120: */
601 }
602/* MAIN ITERATION : SEARCH DIRECTION, STEPLENGTH, LDL'-UPDATE */
603L130:
604 ++(*iter);
605 *mode = 9;
606 if (*iter > itermx) {
607 goto L330;
608 }
609 /* SEARCH DIRECTION AS SOLUTION OF QP - SUBPROBLEM */
610 dcopy___(n, &xl[1], &c__1, &u[1], &c__1);
611 dcopy___(n, &xu[1], &c__1, &v[1], &c__1);
612 d__1 = -one;
613 daxpy_sl__(n, &d__1, &x[1], &c__1, &u[1], &c__1);
614 d__1 = -one;
615 daxpy_sl__(n, &d__1, &x[1], &c__1, &v[1], &c__1);
616 h4 = one;
617 lsq_(m, meq, n, &n3, la, &l[1], &g[1], &a[a_offset], &c__[1], &u[1], &v[1], &s[1], &r__[1], &w[1], &iw[1], mode);
618 /* AUGMENTED PROBLEM FOR INCONSISTENT LINEARIZATION */
619 if (*mode == 6) {
620 if (*n == *meq) {
621 *mode = 4;
622 }
623 }
624 if (*mode == 4) {
625 i__1 = *m;
626 for (j = 1; j <= i__1; ++j) {
627 if (j <= *meq) {
628 a[j + n1 * a_dim1] = -c__[j];
629 } else {
630 /* Computing MAX */
631 d__1 = -c__[j];
632 a[j + n1 * a_dim1] = std::max(d__1, zero);
633 }
634 /* L140: */
635 }
636 s[1] = zero;
637 dcopy___(n, &s[1], &c__0, &s[1], &c__1);
638 h3 = zero;
639 g[n1] = zero;
640 l[n3] = hun;
641 s[n1] = one;
642 u[n1] = zero;
643 v[n1] = one;
644 incons = 0;
645 L150:
646 lsq_(m, meq, &n1, &n3, la, &l[1], &g[1], &a[a_offset], &c__[1], &u[1], &v[1], &s[1], &r__[1], &w[1], &iw[1], mode);
647 h4 = one - s[n1];
648 if (*mode == 4) {
649 l[n3] = ten * l[n3];
650 ++incons;
651 if (incons > 5) {
652 goto L330;
653 }
654 goto L150;
655 } else if (*mode != 1) {
656 goto L330;
657 }
658 } else if (*mode != 1) {
659 goto L330;
660 }
661 /* UPDATE MULTIPLIERS FOR L1-TEST */
662 i__1 = *n;
663 for (i__ = 1; i__ <= i__1; ++i__) {
664 v[i__] = g[i__] - ddot_sl__(m, &a[i__ * a_dim1 + 1], &c__1, &r__[1], &c__1);
665 /* L160: */
666 }
667 f0 = *f;
668 dcopy___(n, &x[1], &c__1, &x0[1], &c__1);
669 gs = ddot_sl__(n, &g[1], &c__1, &s[1], &c__1);
670 h1 = std::abs(gs);
671 h2 = zero;
672 i__1 = *m;
673 for (j = 1; j <= i__1; ++j) {
674 if (j <= *meq) {
675 h3 = c__[j];
676 } else {
677 h3 = zero;
678 }
679 /* Computing MAX */
680 d__1 = -c__[j];
681 h2 += std::max(d__1, h3);
682 h3 = (d__1 = r__[j], std::abs(d__1));
683 /* Computing MAX */
684 d__1 = h3, d__2 = (mu[j] + h3) / two;
685 mu[j] = std::max(d__1, d__2);
686 h1 += h3 * (d__1 = c__[j], std::abs(d__1));
687 /* L170: */
688 }
689 /* CHECK CONVERGENCE */
690 *mode = 0;
691 if (h1 < *acc && h2 < *acc) {
692 goto L330;
693 }
694 h1 = zero;
695 i__1 = *m;
696 for (j = 1; j <= i__1; ++j) {
697 if (j <= *meq) {
698 h3 = c__[j];
699 } else {
700 h3 = zero;
701 }
702 /* Computing MAX */
703 d__1 = -c__[j];
704 h1 += mu[j] * std::max(d__1, h3);
705 /* L180: */
706 }
707 t0 = *f + h1;
708 h3 = gs - h1 * h4;
709 *mode = 8;
710 if (h3 >= zero) {
711 goto L110;
712 }
713 /* LINE SEARCH WITH AN L1-TESTFUNCTION */
714 line = 0;
715 alpha = one;
716 if (iexact == 1) {
717 goto L210;
718 }
719/* INEXACT LINESEARCH */
720L190:
721 ++line;
722 h3 = alpha * h3;
723 dscal_sl__(n, &alpha, &s[1], &c__1);
724 dcopy___(n, &x0[1], &c__1, &x[1], &c__1);
725 daxpy_sl__(n, &one, &s[1], &c__1, &x[1], &c__1);
726 *mode = 1;
727 goto L330;
728L200:
729 if (h1 <= h3 / ten || line > 10) {
730 goto L240;
731 }
732 /* Computing MAX */
733 d__1 = h3 / (two * (h3 - h1));
734 alpha = std::max(d__1, alfmin);
735 goto L190;
736/* EXACT LINESEARCH */
737L210:
738 if (line != 3) {
739 alpha = linmin_(&line, &alfmin, &one, &t, &tol);
740 dcopy___(n, &x0[1], &c__1, &x[1], &c__1);
741 daxpy_sl__(n, &alpha, &s[1], &c__1, &x[1], &c__1);
742 *mode = 1;
743 goto L330;
744 }
745 dscal_sl__(n, &alpha, &s[1], &c__1);
746 goto L240;
747/* CALL FUNCTIONS AT CURRENT X */
748L220:
749 t = *f;
750 i__1 = *m;
751 for (j = 1; j <= i__1; ++j) {
752 if (j <= *meq) {
753 h1 = c__[j];
754 } else {
755 h1 = zero;
756 }
757 /* Computing MAX */
758 d__1 = -c__[j];
759 t += mu[j] * std::max(d__1, h1);
760 /* L230: */
761 }
762 h1 = t - t0;
763 switch (iexact + 1) {
764 case 1:
765 goto L200;
766 case 2:
767 goto L210;
768 }
769/* CHECK CONVERGENCE */
770L240:
771 h3 = zero;
772 i__1 = *m;
773 for (j = 1; j <= i__1; ++j) {
774 if (j <= *meq) {
775 h1 = c__[j];
776 } else {
777 h1 = zero;
778 }
779 /* Computing MAX */
780 d__1 = -c__[j];
781 h3 += std::max(d__1, h1);
782 /* L250: */
783 }
784 if (((d__1 = *f - f0, std::abs(d__1)) < *acc || dnrm2___(n, &s[1], &c__1) < *acc) && h3 < *acc) {
785 *mode = 0;
786 } else {
787 *mode = -1;
788 }
789 goto L330;
790/* CHECK relaxed CONVERGENCE in case of positive directional derivative */
791L255:
792 if (((d__1 = *f - f0, std::abs(d__1)) < tol || dnrm2___(n, &s[1], &c__1) < tol) && h3 < tol) {
793 *mode = 0;
794 } else {
795 *mode = 8;
796 }
797 goto L330;
798/* CALL JACOBIAN AT CURRENT X */
799/* UPDATE CHOLESKY-FACTORS OF HESSIAN MATRIX BY MODIFIED BFGS FORMULA */
800L260:
801 i__1 = *n;
802 for (i__ = 1; i__ <= i__1; ++i__) {
803 u[i__] = g[i__] - ddot_sl__(m, &a[i__ * a_dim1 + 1], &c__1, &r__[1], &c__1) - v[i__];
804 /* L270: */
805 }
806 /* L'*S */
807 k = 0;
808 i__1 = *n;
809 for (i__ = 1; i__ <= i__1; ++i__) {
810 h1 = zero;
811 ++k;
812 i__2 = *n;
813 for (j = i__ + 1; j <= i__2; ++j) {
814 ++k;
815 h1 += l[k] * s[j];
816 /* L280: */
817 }
818 v[i__] = s[i__] + h1;
819 /* L290: */
820 }
821 /* D*L'*S */
822 k = 1;
823 i__1 = *n;
824 for (i__ = 1; i__ <= i__1; ++i__) {
825 v[i__] = l[k] * v[i__];
826 k = k + n1 - i__;
827 /* L300: */
828 }
829 /* L*D*L'*S */
830 for (i__ = *n; i__ >= 1; --i__) {
831 h1 = zero;
832 k = i__;
833 i__1 = i__ - 1;
834 for (j = 1; j <= i__1; ++j) {
835 h1 += l[k] * v[j];
836 k = k + *n - j;
837 /* L310: */
838 }
839 v[i__] += h1;
840 /* L320: */
841 }
842 h1 = ddot_sl__(n, &s[1], &c__1, &u[1], &c__1);
843 h2 = ddot_sl__(n, &s[1], &c__1, &v[1], &c__1);
844 h3 = h2 * .2;
845 if (h1 < h3) {
846 h4 = (h2 - h3) / (h2 - h1);
847 h1 = h3;
848 dscal_sl__(n, &h4, &u[1], &c__1);
849 d__1 = one - h4;
850 daxpy_sl__(n, &d__1, &v[1], &c__1, &u[1], &c__1);
851 }
852 d__1 = one / h1;
853 ldl_(n, &l[1], &u[1], &d__1, &v[1]);
854 d__1 = -one / h2;
855 ldl_(n, &l[1], &v[1], &d__1, &u[1]);
856 /* END OF MAIN ITERATION */
857 goto L130;
858/* END OF SLSQPB */
859L330:
860 return 0;
861} /* slsqpb_ */
862
863int lsq_(int *m, int *meq, int *n, const int *nl, int *la, double *l, double *g, double *a, double *b, double *xl,
864 double *xu, double *x, double *y, double *w, int *jw, int *mode) {
865 /* Initialized data */
866
867 static double zero = 0.;
868 static double one = 1.;
869
870 /* System generated locals */
871 int a_dim1, a_offset, i__1, i__2;
872 double d__1;
873
874 /* Local variables */
875 static int i__, i1, i2, i3, i4, m1, n1, n2, n3, ic, id, ie, if__, ig, ih, il, im, ip, iu, iw;
876 static double diag;
877 static int mineq;
878 static double xnorm;
879
880 /* MINIMIZE with respect to X */
881 /* ||E*X - F|| */
882 /* 1/2 T */
883 /* WITH UPPER TRIANGULAR MATRIX E = +D *L , */
884 /* -1/2 -1 */
885 /* AND VECTOR F = -D *L *G, */
886 /* WHERE THE UNIT LOWER TRIDIANGULAR MATRIX L IS STORED COLUMNWISE */
887 /* DENSE IN THE N*(N+1)/2 ARRAY L WITH VECTOR D STORED IN ITS */
888 /* 'DIAGONAL' THUS SUBSTITUTING THE ONE-ELEMENTS OF L */
889 /* SUBJECT TO */
890 /* A(J)*X - B(J) = 0 , J=1,...,MEQ, */
891 /* A(J)*X - B(J) >=0, J=MEQ+1,...,M, */
892 /* XL(I) <= X(I) <= XU(I), I=1,...,N, */
893 /* ON ENTRY, THE USER HAS TO PROVIDE THE ARRAYS L, G, A, B, XL, XU. */
894 /* WITH DIMENSIONS: L(N*(N+1)/2), G(N), A(LA,N), B(M), XL(N), XU(N) */
895 /* THE WORKING ARRAY W MUST HAVE AT LEAST THE FOLLOWING DIMENSION: */
896 /* DIM(W) = (3*N+M)*(N+1) for LSQ */
897 /* +(N-MEQ+1)*(MINEQ+2) + 2*MINEQ for LSI */
898 /* +(N+MINEQ)*(N-MEQ) + 2*MEQ + N for LSEI */
899 /* with MINEQ = M - MEQ + 2*N */
900 /* ON RETURN, NO ARRAY WILL BE CHANGED BY THE SUBROUTINE. */
901 /* X STORES THE N-DIMENSIONAL SOLUTION VECTOR */
902 /* Y STORES THE VECTOR OF LAGRANGE MULTIPLIERS OF DIMENSION */
903 /* M+N+N (CONSTRAINTS+LOWER+UPPER BOUNDS) */
904 /* MODE IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING MEANINGS: */
905 /* MODE=1: SUCCESSFUL COMPUTATION */
906 /* 2: ERROR RETURN BECAUSE OF WRONG DIMENSIONS (N<1) */
907 /* 3: ITERATION COUNT EXCEEDED BY NNLS */
908 /* 4: INEQUALITY CONSTRAINTS INCOMPATIBLE */
909 /* 5: MATRIX E IS NOT OF FULL RANK */
910 /* 6: MATRIX C IS NOT OF FULL RANK */
911 /* 7: RANK DEFECT IN HFTI */
912 /* coded Dieter Kraft, april 1987 */
913 /* revised march 1989 */
914 /* Parameter adjustments */
915 --y;
916 --x;
917 --xu;
918 --xl;
919 --g;
920 --l;
921 --b;
922 a_dim1 = *la;
923 a_offset = 1 + a_dim1;
924 a -= a_offset;
925 --w;
926 --jw;
927
928 /* Function Body */
929 n1 = *n + 1;
930 mineq = *m - *meq;
931 m1 = mineq + *n + *n;
932 /* determine whether to solve problem */
933 /* with inconsistent linerarization (n2=1) */
934 /* or not (n2=0) */
935 n2 = n1 * *n / 2 + 1;
936 if (n2 == *nl) {
937 n2 = 0;
938 } else {
939 n2 = 1;
940 }
941 n3 = *n - n2;
942 /* RECOVER MATRIX E AND VECTOR F FROM L AND G */
943 i2 = 1;
944 i3 = 1;
945 i4 = 1;
946 ie = 1;
947 if__ = *n * *n + 1;
948 i__1 = n3;
949 for (i__ = 1; i__ <= i__1; ++i__) {
950 i1 = n1 - i__;
951 diag = sqrt(l[i2]);
952 w[i3] = zero;
953 dcopy___(&i1, &w[i3], &c__0, &w[i3], &c__1);
954 i__2 = i1 - n2;
955 dcopy___(&i__2, &l[i2], &c__1, &w[i3], n);
956 i__2 = i1 - n2;
957 dscal_sl__(&i__2, &diag, &w[i3], n);
958 w[i3] = diag;
959 i__2 = i__ - 1;
960 w[if__ - 1 + i__] = (g[i__] - ddot_sl__(&i__2, &w[i4], &c__1, &w[if__], &c__1)) / diag;
961 i2 = i2 + i1 - n2;
962 i3 += n1;
963 i4 += *n;
964 /* L10: */
965 }
966 if (n2 == 1) {
967 w[i3] = l[*nl];
968 w[i4] = zero;
969 dcopy___(&n3, &w[i4], &c__0, &w[i4], &c__1);
970 w[if__ - 1 + *n] = zero;
971 }
972 d__1 = -one;
973 dscal_sl__(n, &d__1, &w[if__], &c__1);
974 ic = if__ + *n;
975 id = ic + *meq * *n;
976 if (*meq > 0) {
977 /* RECOVER MATRIX C FROM UPPER PART OF A */
978 i__1 = *meq;
979 for (i__ = 1; i__ <= i__1; ++i__) {
980 dcopy___(n, &a[i__ + a_dim1], la, &w[ic - 1 + i__], meq);
981 /* L20: */
982 }
983 /* RECOVER VECTOR D FROM UPPER PART OF B */
984 dcopy___(meq, &b[1], &c__1, &w[id], &c__1);
985 d__1 = -one;
986 dscal_sl__(meq, &d__1, &w[id], &c__1);
987 }
988 ig = id + *meq;
989 if (mineq > 0) {
990 /* RECOVER MATRIX G FROM LOWER PART OF A */
991 i__1 = mineq;
992 for (i__ = 1; i__ <= i__1; ++i__) {
993 dcopy___(n, &a[*meq + i__ + a_dim1], la, &w[ig - 1 + i__], &m1);
994 /* L30: */
995 }
996 }
997 /* AUGMENT MATRIX G BY +I AND -I */
998 ip = ig + mineq;
999 i__1 = *n;
1000 for (i__ = 1; i__ <= i__1; ++i__) {
1001 w[ip - 1 + i__] = zero;
1002 dcopy___(n, &w[ip - 1 + i__], &c__0, &w[ip - 1 + i__], &m1);
1003 /* L40: */
1004 }
1005 w[ip] = one;
1006 i__1 = m1 + 1;
1007 dcopy___(n, &w[ip], &c__0, &w[ip], &i__1);
1008 im = ip + *n;
1009 i__1 = *n;
1010 for (i__ = 1; i__ <= i__1; ++i__) {
1011 w[im - 1 + i__] = zero;
1012 dcopy___(n, &w[im - 1 + i__], &c__0, &w[im - 1 + i__], &m1);
1013 /* L50: */
1014 }
1015 w[im] = -one;
1016 i__1 = m1 + 1;
1017 dcopy___(n, &w[im], &c__0, &w[im], &i__1);
1018 ih = ig + m1 * *n;
1019 if (mineq > 0) {
1020 /* RECOVER H FROM LOWER PART OF B */
1021 dcopy___(&mineq, &b[*meq + 1], &c__1, &w[ih], &c__1);
1022 d__1 = -one;
1023 dscal_sl__(&mineq, &d__1, &w[ih], &c__1);
1024 }
1025 /* AUGMENT VECTOR H BY XL AND XU */
1026 il = ih + mineq;
1027 dcopy___(n, &xl[1], &c__1, &w[il], &c__1);
1028 iu = il + *n;
1029 dcopy___(n, &xu[1], &c__1, &w[iu], &c__1);
1030 d__1 = -one;
1031 dscal_sl__(n, &d__1, &w[iu], &c__1);
1032 iw = iu + *n;
1033 i__1 = std::max(1, *meq);
1034 lsei_(&w[ic], &w[id], &w[ie], &w[if__], &w[ig], &w[ih], &i__1, meq, n, n, &m1, &m1, n, &x[1], &xnorm, &w[iw], &jw[1],
1035 mode);
1036 if (*mode == 1) {
1037 /* restore Lagrange multipliers */
1038 dcopy___(m, &w[iw], &c__1, &y[1], &c__1);
1039 dcopy___(&n3, &w[iw + *m], &c__1, &y[*m + 1], &c__1);
1040 dcopy___(&n3, &w[iw + *m + *n], &c__1, &y[*m + n3 + 1], &c__1);
1041 }
1042 /* END OF SUBROUTINE LSQ */
1043 return 0;
1044} /* lsq_ */
1045
1046int lsei_(double *c__, double *d__, double *e, double *f, double *g, double *h__, int *lc, int *mc, int *le, int *me,
1047 int *lg, int *mg, int *n, double *x, double *xnrm, double *w, int *jw, int *mode) {
1048 /* Initialized data */
1049
1050 static double epmach = 2.22e-16;
1051 static double zero = 0.;
1052
1053 /* System generated locals */
1054 int c_dim1, c_offset, e_dim1, e_offset, g_dim1, g_offset, i__1, i__2, i__3;
1055 double d__1;
1056
1057 /* Local variables */
1058 static int i__, j, k, l;
1059 static double t;
1060 static int ie, if__, ig, iw, mc1;
1061 static int krank;
1062
1063 /* FOR MODE=1, THE SUBROUTINE RETURNS THE SOLUTION X OF */
1064 /* EQUALITY & INEQUALITY CONSTRAINED LEAST SQUARES PROBLEM LSEI : */
1065 /* MIN ||E*X - F|| */
1066 /* X */
1067 /* S.T. C*X = D, */
1068 /* G*X >= H. */
1069 /* USING QR DECOMPOSITION & ORTHOGONAL BASIS OF NULLSPACE OF C */
1070 /* CHAPTER 23.6 OF LAWSON & HANSON: SOLVING LEAST SQUARES PROBLEMS. */
1071 /* THE FOLLOWING DIMENSIONS OF THE ARRAYS DEFINING THE PROBLEM */
1072 /* ARE NECESSARY */
1073 /* DIM(E) : FORMAL (LE,N), ACTUAL (ME,N) */
1074 /* DIM(F) : FORMAL (LE ), ACTUAL (ME ) */
1075 /* DIM(C) : FORMAL (LC,N), ACTUAL (MC,N) */
1076 /* DIM(D) : FORMAL (LC ), ACTUAL (MC ) */
1077 /* DIM(G) : FORMAL (LG,N), ACTUAL (MG,N) */
1078 /* DIM(H) : FORMAL (LG ), ACTUAL (MG ) */
1079 /* DIM(X) : FORMAL (N ), ACTUAL (N ) */
1080 /* DIM(W) : 2*MC+ME+(ME+MG)*(N-MC) for LSEI */
1081 /* +(N-MC+1)*(MG+2)+2*MG for LSI */
1082 /* DIM(JW): MAX(MG,L) */
1083 /* ON ENTRY, THE USER HAS TO PROVIDE THE ARRAYS C, D, E, F, G, AND H. */
1084 /* ON RETURN, ALL ARRAYS WILL BE CHANGED BY THE SUBROUTINE. */
1085 /* X STORES THE SOLUTION VECTOR */
1086 /* XNORM STORES THE RESIDUUM OF THE SOLUTION IN EUCLIDIAN NORM */
1087 /* W STORES THE VECTOR OF LAGRANGE MULTIPLIERS IN ITS FIRST */
1088 /* MC+MG ELEMENTS */
1089 /* MODE IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING MEANINGS: */
1090 /* MODE=1: SUCCESSFUL COMPUTATION */
1091 /* 2: ERROR RETURN BECAUSE OF WRONG DIMENSIONS (N<1) */
1092 /* 3: ITERATION COUNT EXCEEDED BY NNLS */
1093 /* 4: INEQUALITY CONSTRAINTS INCOMPATIBLE */
1094 /* 5: MATRIX E IS NOT OF FULL RANK */
1095 /* 6: MATRIX C IS NOT OF FULL RANK */
1096 /* 7: RANK DEFECT IN HFTI */
1097 /* 18.5.1981, DIETER KRAFT, DFVLR OBERPFAFFENHOFEN */
1098 /* 20.3.1987, DIETER KRAFT, DFVLR OBERPFAFFENHOFEN */
1099 /* Parameter adjustments */
1100 --d__;
1101 --f;
1102 --h__;
1103 --x;
1104 g_dim1 = *lg;
1105 g_offset = 1 + g_dim1;
1106 g -= g_offset;
1107 e_dim1 = *le;
1108 e_offset = 1 + e_dim1;
1109 e -= e_offset;
1110 c_dim1 = *lc;
1111 c_offset = 1 + c_dim1;
1112 c__ -= c_offset;
1113 --w;
1114 --jw;
1115
1116 /* Function Body */
1117 *mode = 2;
1118 if (*mc > *n) {
1119 goto L75;
1120 }
1121 l = *n - *mc;
1122 mc1 = *mc + 1;
1123 iw = (l + 1) * (*mg + 2) + (*mg << 1) + *mc;
1124 ie = iw + *mc + 1;
1125 if__ = ie + *me * l;
1126 ig = if__ + *me;
1127 /* TRIANGULARIZE C AND APPLY FACTORS TO E AND G */
1128 i__1 = *mc;
1129 for (i__ = 1; i__ <= i__1; ++i__) {
1130 /* Computing MIN */
1131 i__2 = i__ + 1;
1132 j = std::min(i__2, *lc);
1133 i__2 = i__ + 1;
1134 i__3 = *mc - i__;
1135 h12_(&c__1, &i__, &i__2, n, &c__[i__ + c_dim1], lc, &w[iw + i__], &c__[j + c_dim1], lc, &c__1, &i__3);
1136 i__2 = i__ + 1;
1137 h12_(&c__2, &i__, &i__2, n, &c__[i__ + c_dim1], lc, &w[iw + i__], &e[e_offset], le, &c__1, me);
1138 /* L10: */
1139 i__2 = i__ + 1;
1140 h12_(&c__2, &i__, &i__2, n, &c__[i__ + c_dim1], lc, &w[iw + i__], &g[g_offset], lg, &c__1, mg);
1141 }
1142 /* SOLVE C*X=D AND MODIFY F */
1143 *mode = 6;
1144 i__2 = *mc;
1145 for (i__ = 1; i__ <= i__2; ++i__) {
1146 if ((d__1 = c__[i__ + i__ * c_dim1], std::abs(d__1)) < epmach) {
1147 goto L75;
1148 }
1149 i__1 = i__ - 1;
1150 x[i__] = (d__[i__] - ddot_sl__(&i__1, &c__[i__ + c_dim1], lc, &x[1], &c__1)) / c__[i__ + i__ * c_dim1];
1151 /* L15: */
1152 }
1153 *mode = 1;
1154 w[mc1] = zero;
1155 i__2 = *mg - *mc;
1156 dcopy___(&i__2, &w[mc1], &c__0, &w[mc1], &c__1);
1157 if (*mc == *n) {
1158 goto L50;
1159 }
1160 i__2 = *me;
1161 for (i__ = 1; i__ <= i__2; ++i__) {
1162 /* L20: */
1163 w[if__ - 1 + i__] = f[i__] - ddot_sl__(mc, &e[i__ + e_dim1], le, &x[1], &c__1);
1164 }
1165 /* STORE TRANSFORMED E & G */
1166 i__2 = *me;
1167 for (i__ = 1; i__ <= i__2; ++i__) {
1168 /* L25: */
1169 dcopy___(&l, &e[i__ + mc1 * e_dim1], le, &w[ie - 1 + i__], me);
1170 }
1171 i__2 = *mg;
1172 for (i__ = 1; i__ <= i__2; ++i__) {
1173 /* L30: */
1174 dcopy___(&l, &g[i__ + mc1 * g_dim1], lg, &w[ig - 1 + i__], mg);
1175 }
1176 if (*mg > 0) {
1177 goto L40;
1178 }
1179 /* SOLVE LS WITHOUT INEQUALITY CONSTRAINTS */
1180 *mode = 7;
1181 k = std::max(*le, *n);
1182 t = sqrt(epmach);
1183 hfti_(&w[ie], me, me, &l, &w[if__], &k, &c__1, &t, &krank, xnrm, &w[1], &w[l + 1], &jw[1]);
1184 dcopy___(&l, &w[if__], &c__1, &x[mc1], &c__1);
1185 if (krank != l) {
1186 goto L75;
1187 }
1188 *mode = 1;
1189 goto L50;
1190/* MODIFY H AND SOLVE INEQUALITY CONSTRAINED LS PROBLEM */
1191L40:
1192 i__2 = *mg;
1193 for (i__ = 1; i__ <= i__2; ++i__) {
1194 /* L45: */
1195 h__[i__] -= ddot_sl__(mc, &g[i__ + g_dim1], lg, &x[1], &c__1);
1196 }
1197 lsi_(&w[ie], &w[if__], &w[ig], &h__[1], me, me, mg, mg, &l, &x[mc1], xnrm, &w[mc1], &jw[1], mode);
1198 if (*mc == 0) {
1199 goto L75;
1200 }
1201 t = dnrm2___(mc, &x[1], &c__1);
1202 *xnrm = sqrt(*xnrm * *xnrm + t * t);
1203 if (*mode != 1) {
1204 goto L75;
1205 }
1206/* SOLUTION OF ORIGINAL PROBLEM AND LAGRANGE MULTIPLIERS */
1207L50:
1208 i__2 = *me;
1209 for (i__ = 1; i__ <= i__2; ++i__) {
1210 /* L55: */
1211 f[i__] = ddot_sl__(n, &e[i__ + e_dim1], le, &x[1], &c__1) - f[i__];
1212 }
1213 i__2 = *mc;
1214 for (i__ = 1; i__ <= i__2; ++i__) {
1215 /* L60: */
1216 d__[i__] = ddot_sl__(me, &e[i__ * e_dim1 + 1], &c__1, &f[1], &c__1) -
1217 ddot_sl__(mg, &g[i__ * g_dim1 + 1], &c__1, &w[mc1], &c__1);
1218 }
1219 for (i__ = *mc; i__ >= 1; --i__) {
1220 /* L65: */
1221 i__2 = i__ + 1;
1222 h12_(&c__2, &i__, &i__2, n, &c__[i__ + c_dim1], lc, &w[iw + i__], &x[1], &c__1, &c__1, &c__1);
1223 }
1224 for (i__ = *mc; i__ >= 1; --i__) {
1225 /* Computing MIN */
1226 i__2 = i__ + 1;
1227 j = std::min(i__2, *lc);
1228 i__2 = *mc - i__;
1229 w[i__] = (d__[i__] - ddot_sl__(&i__2, &c__[j + i__ * c_dim1], &c__1, &w[j], &c__1)) / c__[i__ + i__ * c_dim1];
1230 /* L70: */
1231 }
1232/* END OF SUBROUTINE LSEI */
1233L75:
1234 return 0;
1235} /* lsei_ */
1236
1237int lsi_(double *e, double *f, double *g, double *h__, int *le, int *me, int *lg, int *mg, int *n, double *x,
1238 double *xnorm, double *w, int *jw, int *mode) {
1239 /* Initialized data */
1240
1241 static double epmach = 2.22e-16;
1242 static double one = 1.;
1243
1244 /* System generated locals */
1245 int e_dim1, e_offset, g_dim1, g_offset, i__1, i__2, i__3;
1246 double d__1;
1247
1248 /* Local variables */
1249 static int i__, j;
1250 static double t;
1251
1252 /* FOR MODE=1, THE SUBROUTINE RETURNS THE SOLUTION X OF */
1253 /* INEQUALITY CONSTRAINED LINEAR LEAST SQUARES PROBLEM: */
1254 /* MIN ||E*X-F|| */
1255 /* X */
1256 /* S.T. G*X >= H */
1257 /* THE ALGORITHM IS BASED ON QR DECOMPOSITION AS DESCRIBED IN */
1258 /* CHAPTER 23.5 OF LAWSON & HANSON: SOLVING LEAST SQUARES PROBLEMS */
1259 /* THE FOLLOWING DIMENSIONS OF THE ARRAYS DEFINING THE PROBLEM */
1260 /* ARE NECESSARY */
1261 /* DIM(E) : FORMAL (LE,N), ACTUAL (ME,N) */
1262 /* DIM(F) : FORMAL (LE ), ACTUAL (ME ) */
1263 /* DIM(G) : FORMAL (LG,N), ACTUAL (MG,N) */
1264 /* DIM(H) : FORMAL (LG ), ACTUAL (MG ) */
1265 /* DIM(X) : N */
1266 /* DIM(W) : (N+1)*(MG+2) + 2*MG */
1267 /* DIM(JW): LG */
1268 /* ON ENTRY, THE USER HAS TO PROVIDE THE ARRAYS E, F, G, AND H. */
1269 /* ON RETURN, ALL ARRAYS WILL BE CHANGED BY THE SUBROUTINE. */
1270 /* X STORES THE SOLUTION VECTOR */
1271 /* XNORM STORES THE RESIDUUM OF THE SOLUTION IN EUCLIDIAN NORM */
1272 /* W STORES THE VECTOR OF LAGRANGE MULTIPLIERS IN ITS FIRST */
1273 /* MG ELEMENTS */
1274 /* MODE IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING MEANINGS: */
1275 /* MODE=1: SUCCESSFUL COMPUTATION */
1276 /* 2: ERROR RETURN BECAUSE OF WRONG DIMENSIONS (N<1) */
1277 /* 3: ITERATION COUNT EXCEEDED BY NNLS */
1278 /* 4: INEQUALITY CONSTRAINTS INCOMPATIBLE */
1279 /* 5: MATRIX E IS NOT OF FULL RANK */
1280 /* 03.01.1980, DIETER KRAFT: CODED */
1281 /* 20.03.1987, DIETER KRAFT: REVISED TO FORTRAN 77 */
1282 /* Parameter adjustments */
1283 --f;
1284 --jw;
1285 --h__;
1286 --x;
1287 g_dim1 = *lg;
1288 g_offset = 1 + g_dim1;
1289 g -= g_offset;
1290 e_dim1 = *le;
1291 e_offset = 1 + e_dim1;
1292 e -= e_offset;
1293 --w;
1294
1295 /* Function Body */
1296 /* QR-FACTORS OF E AND APPLICATION TO F */
1297 i__1 = *n;
1298 for (i__ = 1; i__ <= i__1; ++i__) {
1299 /* Computing MIN */
1300 i__2 = i__ + 1;
1301 j = std::min(i__2, *n);
1302 i__2 = i__ + 1;
1303 i__3 = *n - i__;
1304 h12_(&c__1, &i__, &i__2, me, &e[i__ * e_dim1 + 1], &c__1, &t, &e[j * e_dim1 + 1], &c__1, le, &i__3);
1305 /* L10: */
1306 i__2 = i__ + 1;
1307 h12_(&c__2, &i__, &i__2, me, &e[i__ * e_dim1 + 1], &c__1, &t, &f[1], &c__1, &c__1, &c__1);
1308 }
1309 /* TRANSFORM G AND H TO GET LEAST DISTANCE PROBLEM */
1310 *mode = 5;
1311 i__2 = *mg;
1312 for (i__ = 1; i__ <= i__2; ++i__) {
1313 i__1 = *n;
1314 for (j = 1; j <= i__1; ++j) {
1315 if ((d__1 = e[j + j * e_dim1], std::abs(d__1)) < epmach) {
1316 goto L50;
1317 }
1318 /* L20: */
1319 i__3 = j - 1;
1320 g[i__ + j * g_dim1] =
1321 (g[i__ + j * g_dim1] - ddot_sl__(&i__3, &g[i__ + g_dim1], lg, &e[j * e_dim1 + 1], &c__1)) / e[j + j * e_dim1];
1322 }
1323 /* L30: */
1324 h__[i__] -= ddot_sl__(n, &g[i__ + g_dim1], lg, &f[1], &c__1);
1325 }
1326 /* SOLVE LEAST DISTANCE PROBLEM */
1327 ldp_(&g[g_offset], lg, mg, n, &h__[1], &x[1], xnorm, &w[1], &jw[1], mode);
1328 if (*mode != 1) {
1329 goto L50;
1330 }
1331 /* SOLUTION OF ORIGINAL PROBLEM */
1332 daxpy_sl__(n, &one, &f[1], &c__1, &x[1], &c__1);
1333 for (i__ = *n; i__ >= 1; --i__) {
1334 /* Computing MIN */
1335 i__2 = i__ + 1;
1336 j = std::min(i__2, *n);
1337 /* L40: */
1338 i__2 = *n - i__;
1339 x[i__] = (x[i__] - ddot_sl__(&i__2, &e[i__ + j * e_dim1], le, &x[j], &c__1)) / e[i__ + i__ * e_dim1];
1340 }
1341 /* Computing MIN */
1342 i__2 = *n + 1;
1343 j = std::min(i__2, *me);
1344 i__2 = *me - *n;
1345 t = dnrm2___(&i__2, &f[j], &c__1);
1346 *xnorm = sqrt(*xnorm * *xnorm + t * t);
1347/* END OF SUBROUTINE LSI */
1348L50:
1349 return 0;
1350} /* lsi_ */
1351
1352int ldp_(double *g, const int *mg, int *m, int *n, double *h__, double *x, double *xnorm, double *w, int *index,
1353 int *mode) {
1354 /* Initialized data */
1355
1356 static double zero = 0.;
1357 static double one = 1.;
1358
1359 /* System generated locals */
1360 int g_dim1, g_offset, i__1, i__2;
1361 double d__1;
1362
1363 /* Local variables */
1364 static int i__, j, n1, if__, iw, iy, iz;
1365 static double fac;
1366 static double rnorm;
1367 static int iwdual;
1368
1369 /* T */
1370 /* MINIMIZE 1/2 X X SUBJECT TO G * X >= H. */
1371 /* C.L. LAWSON, R.J. HANSON: 'SOLVING LEAST SQUARES PROBLEMS' */
1372 /* PRENTICE HALL, ENGLEWOOD CLIFFS, NEW JERSEY, 1974. */
1373 /* PARAMETER DESCRIPTION: */
1374 /* G(),MG,M,N ON ENTRY G() STORES THE M BY N MATRIX OF */
1375 /* LINEAR INEQUALITY CONSTRAINTS. G() HAS FIRST */
1376 /* DIMENSIONING PARAMETER MG */
1377 /* H() ON ENTRY H() STORES THE M VECTOR H REPRESENTING */
1378 /* THE RIGHT SIDE OF THE INEQUALITY SYSTEM */
1379 /* REMARK: G(),H() WILL NOT BE CHANGED DURING CALCULATIONS BY LDP */
1380 /* X() ON ENTRY X() NEED NOT BE INITIALIZED. */
1381 /* ON EXIT X() STORES THE SOLUTION VECTOR X IF MODE=1. */
1382 /* XNORM ON EXIT XNORM STORES THE EUCLIDIAN NORM OF THE */
1383 /* SOLUTION VECTOR IF COMPUTATION IS SUCCESSFUL */
1384 /* W() W IS A ONE DIMENSIONAL WORKING SPACE, THE LENGTH */
1385 /* OF WHICH SHOULD BE AT LEAST (M+2)*(N+1) + 2*M */
1386 /* ON EXIT W() STORES THE LAGRANGE MULTIPLIERS */
1387 /* ASSOCIATED WITH THE CONSTRAINTS */
1388 /* AT THE SOLUTION OF PROBLEM LDP */
1389 /* INDEX() INDEX() IS A ONE DIMENSIONAL INT WORKING SPACE */
1390 /* OF LENGTH AT LEAST M */
1391 /* MODE MODE IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING */
1392 /* MEANINGS: */
1393 /* MODE=1: SUCCESSFUL COMPUTATION */
1394 /* 2: ERROR RETURN BECAUSE OF WRONG DIMENSIONS (N.LE.0) */
1395 /* 3: ITERATION COUNT EXCEEDED BY NNLS */
1396 /* 4: INEQUALITY CONSTRAINTS INCOMPATIBLE */
1397 /* Parameter adjustments */
1398 --index;
1399 --h__;
1400 --x;
1401 g_dim1 = *mg;
1402 g_offset = 1 + g_dim1;
1403 g -= g_offset;
1404 --w;
1405
1406 /* Function Body */
1407 *mode = 2;
1408 if (*n <= 0) {
1409 goto L50;
1410 }
1411 /* STATE DUAL PROBLEM */
1412 *mode = 1;
1413 x[1] = zero;
1414 dcopy___(n, &x[1], &c__0, &x[1], &c__1);
1415 *xnorm = zero;
1416 if (*m == 0) {
1417 goto L50;
1418 }
1419 iw = 0;
1420 i__1 = *m;
1421 for (j = 1; j <= i__1; ++j) {
1422 i__2 = *n;
1423 for (i__ = 1; i__ <= i__2; ++i__) {
1424 ++iw;
1425 /* L10: */
1426 w[iw] = g[j + i__ * g_dim1];
1427 }
1428 ++iw;
1429 /* L20: */
1430 w[iw] = h__[j];
1431 }
1432 if__ = iw + 1;
1433 i__1 = *n;
1434 for (i__ = 1; i__ <= i__1; ++i__) {
1435 ++iw;
1436 /* L30: */
1437 w[iw] = zero;
1438 }
1439 w[iw + 1] = one;
1440 n1 = *n + 1;
1441 iz = iw + 2;
1442 iy = iz + n1;
1443 iwdual = iy + *m;
1444 /* SOLVE DUAL PROBLEM */
1445 nnls_(&w[1], &n1, &n1, m, &w[if__], &w[iy], &rnorm, &w[iwdual], &w[iz], &index[1], mode);
1446 if (*mode != 1) {
1447 goto L50;
1448 }
1449 *mode = 4;
1450 if (rnorm <= zero) {
1451 goto L50;
1452 }
1453 /* COMPUTE SOLUTION OF PRIMAL PROBLEM */
1454 fac = one - ddot_sl__(m, &h__[1], &c__1, &w[iy], &c__1);
1455 d__1 = one + fac;
1456 if (d__1 - one <= zero) {
1457 goto L50;
1458 }
1459 *mode = 1;
1460 fac = one / fac;
1461 i__1 = *n;
1462 for (j = 1; j <= i__1; ++j) {
1463 /* L40: */
1464 x[j] = fac * ddot_sl__(m, &g[j * g_dim1 + 1], &c__1, &w[iy], &c__1);
1465 }
1466 *xnorm = dnrm2___(n, &x[1], &c__1);
1467 /* COMPUTE LAGRANGE MULTIPLIERS FOR PRIMAL PROBLEM */
1468 w[1] = zero;
1469 dcopy___(m, &w[1], &c__0, &w[1], &c__1);
1470 daxpy_sl__(m, &fac, &w[iy], &c__1, &w[1], &c__1);
1471/* END OF SUBROUTINE LDP */
1472L50:
1473 return 0;
1474} /* ldp_ */
1475
1476int nnls_(double *a, int *mda, int *m, int *n, double *b, double *x, double *rnorm, double *w, double *z__, int *index,
1477 int *mode) {
1478 /* Initialized data */
1479
1480 static double zero = 0.;
1481 static double one = 1.;
1482 static double factor = .01;
1483
1484 /* System generated locals */
1485 int a_dim1, a_offset, i__1, i__2;
1486 double d__1;
1487
1488 /* Local variables */
1489 static double c__;
1490 static int i__, j, k, l;
1491 static double s, t;
1492 static int ii, jj, ip, iz, jz;
1493 static double up;
1494 static int iz1, iz2, npp1, iter;
1495 static double wmax, alpha, asave;
1496 static int itmax, izmax, nsetp;
1497 static double unorm;
1498
1499 /* C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY: */
1500 /* 'SOLVING LEAST SQUARES PROBLEMS'. PRENTICE-HALL.1974 */
1501 /* ********** NONNEGATIVE LEAST SQUARES ********** */
1502 /* GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN */
1503 /* N-VECTOR, X, WHICH SOLVES THE LEAST SQUARES PROBLEM */
1504 /* A*X = B SUBJECT TO X >= 0 */
1505 /* A(),MDA,M,N */
1506 /* MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE ARRAY,A(). */
1507 /* ON ENTRY A() CONTAINS THE M BY N MATRIX,A. */
1508 /* ON EXIT A() CONTAINS THE PRODUCT Q*A, */
1509 /* WHERE Q IS AN M BY M ORTHOGONAL MATRIX GENERATED */
1510 /* IMPLICITLY BY THIS SUBROUTINE. */
1511 /* EITHER M>=N OR M<N IS PERMISSIBLE. */
1512 /* THERE IS NO RESTRICTION ON THE RANK OF A. */
1513 /* B() ON ENTRY B() CONTAINS THE M-VECTOR, B. */
1514 /* ON EXIT B() CONTAINS Q*B. */
1515 /* X() ON ENTRY X() NEED NOT BE INITIALIZED. */
1516 /* ON EXIT X() WILL CONTAIN THE SOLUTION VECTOR. */
1517 /* RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE */
1518 /* RESIDUAL VECTOR. */
1519 /* W() AN N-ARRAY OF WORKING SPACE. */
1520 /* ON EXIT W() WILL CONTAIN THE DUAL SOLUTION VECTOR. */
1521 /* W WILL SATISFY W(I)=0 FOR ALL I IN SET P */
1522 /* AND W(I)<=0 FOR ALL I IN SET Z */
1523 /* Z() AN M-ARRAY OF WORKING SPACE. */
1524 /* INDEX()AN INT WORKING ARRAY OF LENGTH AT LEAST N. */
1525 /* ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS */
1526 /* P AND Z AS FOLLOWS: */
1527 /* INDEX(1) THRU INDEX(NSETP) = SET P. */
1528 /* INDEX(IZ1) THRU INDEX (IZ2) = SET Z. */
1529 /* IZ1=NSETP + 1 = NPP1, IZ2=N. */
1530 /* MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING MEANING: */
1531 /* 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY. */
1532 /* 2 THE DIMENSIONS OF THE PROBLEM ARE WRONG, */
1533 /* EITHER M <= 0 OR N <= 0. */
1534 /* 3 ITERATION COUNT EXCEEDED, MORE THAN 3*N ITERATIONS. */
1535 /* Parameter adjustments */
1536 --z__;
1537 --b;
1538 --index;
1539 --w;
1540 --x;
1541 a_dim1 = *mda;
1542 a_offset = 1 + a_dim1;
1543 a -= a_offset;
1544
1545 /* Function Body */
1546 /* revised Dieter Kraft, March 1983 */
1547 *mode = 2;
1548 if (*m <= 0 || *n <= 0) {
1549 goto L290;
1550 }
1551 *mode = 1;
1552 iter = 0;
1553 itmax = *n * 3;
1554 /* STEP ONE (INITIALIZE) */
1555 i__1 = *n;
1556 for (i__ = 1; i__ <= i__1; ++i__) {
1557 /* L100: */
1558 index[i__] = i__;
1559 }
1560 iz1 = 1;
1561 iz2 = *n;
1562 nsetp = 0;
1563 npp1 = 1;
1564 x[1] = zero;
1565 dcopy___(n, &x[1], &c__0, &x[1], &c__1);
1566/* STEP TWO (COMPUTE DUAL VARIABLES) */
1567/* .....ENTRY LOOP A */
1568L110:
1569 if (iz1 > iz2 || nsetp >= *m) {
1570 goto L280;
1571 }
1572 i__1 = iz2;
1573 for (iz = iz1; iz <= i__1; ++iz) {
1574 j = index[iz];
1575 /* L120: */
1576 i__2 = *m - nsetp;
1577 w[j] = ddot_sl__(&i__2, &a[npp1 + j * a_dim1], &c__1, &b[npp1], &c__1);
1578 }
1579/* STEP THREE (TEST DUAL VARIABLES) */
1580L130:
1581 wmax = zero;
1582 i__2 = iz2;
1583 for (iz = iz1; iz <= i__2; ++iz) {
1584 j = index[iz];
1585 if (w[j] <= wmax) {
1586 goto L140;
1587 }
1588 wmax = w[j];
1589 izmax = iz;
1590 L140:;
1591 }
1592 /* .....EXIT LOOP A */
1593 if (wmax <= zero) {
1594 goto L280;
1595 }
1596 iz = izmax;
1597 j = index[iz];
1598 /* STEP FOUR (TEST INDEX J FOR LINEAR DEPENDENCY) */
1599 asave = a[npp1 + j * a_dim1];
1600 i__2 = npp1 + 1;
1601 h12_(&c__1, &npp1, &i__2, m, &a[j * a_dim1 + 1], &c__1, &up, &z__[1], &c__1, &c__1, &c__0);
1602 unorm = dnrm2___(&nsetp, &a[j * a_dim1 + 1], &c__1);
1603 t = factor * (d__1 = a[npp1 + j * a_dim1], std::abs(d__1));
1604 d__1 = unorm + t;
1605 if (d__1 - unorm <= zero) {
1606 goto L150;
1607 }
1608 dcopy___(m, &b[1], &c__1, &z__[1], &c__1);
1609 i__2 = npp1 + 1;
1610 h12_(&c__2, &npp1, &i__2, m, &a[j * a_dim1 + 1], &c__1, &up, &z__[1], &c__1, &c__1, &c__1);
1611 if (z__[npp1] / a[npp1 + j * a_dim1] > zero) {
1612 goto L160;
1613 }
1614L150:
1615 a[npp1 + j * a_dim1] = asave;
1616 w[j] = zero;
1617 goto L130;
1618/* STEP FIVE (ADD COLUMN) */
1619L160:
1620 dcopy___(m, &z__[1], &c__1, &b[1], &c__1);
1621 index[iz] = index[iz1];
1622 index[iz1] = j;
1623 ++iz1;
1624 nsetp = npp1;
1625 ++npp1;
1626 i__2 = iz2;
1627 for (jz = iz1; jz <= i__2; ++jz) {
1628 jj = index[jz];
1629 /* L170: */
1630 h12_(&c__2, &nsetp, &npp1, m, &a[j * a_dim1 + 1], &c__1, &up, &a[jj * a_dim1 + 1], &c__1, mda, &c__1);
1631 }
1632 k = std::min(npp1, *mda);
1633 w[j] = zero;
1634 i__2 = *m - nsetp;
1635 dcopy___(&i__2, &w[j], &c__0, &a[k + j * a_dim1], &c__1);
1636/* STEP SIX (SOLVE LEAST SQUARES SUB-PROBLEM) */
1637/* .....ENTRY LOOP B */
1638L180:
1639 for (ip = nsetp; ip >= 1; --ip) {
1640 if (ip == nsetp) {
1641 goto L190;
1642 }
1643 d__1 = -z__[ip + 1];
1644 daxpy_sl__(&ip, &d__1, &a[jj * a_dim1 + 1], &c__1, &z__[1], &c__1);
1645 L190:
1646 jj = index[ip];
1647 /* L200: */
1648 z__[ip] /= a[ip + jj * a_dim1];
1649 }
1650 ++iter;
1651 if (iter <= itmax) {
1652 goto L220;
1653 }
1654L210:
1655 *mode = 3;
1656 goto L280;
1657/* STEP SEVEN TO TEN (STEP LENGTH ALGORITHM) */
1658L220:
1659 alpha = one;
1660 jj = 0;
1661 i__2 = nsetp;
1662 for (ip = 1; ip <= i__2; ++ip) {
1663 if (z__[ip] > zero) {
1664 goto L230;
1665 }
1666 l = index[ip];
1667 t = -x[l] / (z__[ip] - x[l]);
1668 if (alpha < t) {
1669 goto L230;
1670 }
1671 alpha = t;
1672 jj = ip;
1673 L230:;
1674 }
1675 i__2 = nsetp;
1676 for (ip = 1; ip <= i__2; ++ip) {
1677 l = index[ip];
1678 /* L240: */
1679 x[l] = (one - alpha) * x[l] + alpha * z__[ip];
1680 }
1681 /* .....EXIT LOOP B */
1682 if (jj == 0) {
1683 goto L110;
1684 }
1685 /* STEP ELEVEN (DELETE COLUMN) */
1686 i__ = index[jj];
1687L250:
1688 x[i__] = zero;
1689 ++jj;
1690 i__2 = nsetp;
1691 for (j = jj; j <= i__2; ++j) {
1692 ii = index[j];
1693 index[j - 1] = ii;
1694 dsrotg_(&a[j - 1 + ii * a_dim1], &a[j + ii * a_dim1], &c__, &s);
1695 t = a[j - 1 + ii * a_dim1];
1696 dsrot_(n, &a[j - 1 + a_dim1], mda, &a[j + a_dim1], mda, &c__, &s);
1697 a[j - 1 + ii * a_dim1] = t;
1698 a[j + ii * a_dim1] = zero;
1699 /* L260: */
1700 dsrot_(&c__1, &b[j - 1], &c__1, &b[j], &c__1, &c__, &s);
1701 }
1702 npp1 = nsetp;
1703 --nsetp;
1704 --iz1;
1705 index[iz1] = i__;
1706 if (nsetp <= 0) {
1707 goto L210;
1708 }
1709 i__2 = nsetp;
1710 for (jj = 1; jj <= i__2; ++jj) {
1711 i__ = index[jj];
1712 if (x[i__] <= zero) {
1713 goto L250;
1714 }
1715 /* L270: */
1716 }
1717 dcopy___(m, &b[1], &c__1, &z__[1], &c__1);
1718 goto L180;
1719/* STEP TWELVE (SOLUTION) */
1720L280:
1721 k = std::min(npp1, *m);
1722 i__2 = *m - nsetp;
1723 *rnorm = dnrm2___(&i__2, &b[k], &c__1);
1724 if (npp1 > *m) {
1725 w[1] = zero;
1726 dcopy___(n, &w[1], &c__0, &w[1], &c__1);
1727 }
1728/* END OF SUBROUTINE NNLS */
1729L290:
1730 return 0;
1731} /* nnls_ */
1732
1733int hfti_(double *a, int *mda, int *m, int *n, double *b, int *mdb, int *nb, const double *tau, int *krank,
1734 double *rnorm, double *h__, double *g, int *ip) {
1735 /* Initialized data */
1736
1737 static double zero = 0.;
1738 static double factor = .001;
1739
1740 /* System generated locals */
1741 int a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
1742 double d__1;
1743
1744 /* Local variables */
1745 static int i__, j, k, l;
1746 static int jb, kp1;
1747 static double tmp, hmax;
1748 static int lmax, ldiag;
1749
1750 /* RANK-DEFICIENT LEAST SQUARES ALGORITHM AS DESCRIBED IN: */
1751 /* C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 */
1752 /* TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974 */
1753 /* A(*,*),MDA,M,N THE ARRAY A INITIALLY CONTAINS THE M x N MATRIX A */
1754 /* OF THE LEAST SQUARES PROBLEM AX = B. */
1755 /* THE FIRST DIMENSIONING PARAMETER MDA MUST SATISFY */
1756 /* MDA >= M. EITHER M >= N OR M < N IS PERMITTED. */
1757 /* THERE IS NO RESTRICTION ON THE RANK OF A. */
1758 /* THE MATRIX A WILL BE MODIFIED BY THE SUBROUTINE. */
1759 /* B(*,*),MDB,NB IF NB = 0 THE SUBROUTINE WILL MAKE NO REFERENCE */
1760 /* TO THE ARRAY B. IF NB > 0 THE ARRAY B() MUST */
1761 /* INITIALLY CONTAIN THE M x NB MATRIX B OF THE */
1762 /* THE LEAST SQUARES PROBLEM AX = B AND ON RETURN */
1763 /* THE ARRAY B() WILL CONTAIN THE N x NB SOLUTION X. */
1764 /* IF NB>1 THE ARRAY B() MUST BE DOUBLE SUBSCRIPTED */
1765 /* WITH FIRST DIMENSIONING PARAMETER MDB>=MAX(M,N), */
1766 /* IF NB=1 THE ARRAY B() MAY BE EITHER SINGLE OR */
1767 /* DOUBLE SUBSCRIPTED. */
1768 /* TAU ABSOLUTE TOLERANCE PARAMETER FOR PSEUDORANK */
1769 /* DETERMINATION, PROVIDED BY THE USER. */
1770 /* KRANK PSEUDORANK OF A, SET BY THE SUBROUTINE. */
1771 /* RNORM ON EXIT, RNORM(J) WILL CONTAIN THE EUCLIDIAN */
1772 /* NORM OF THE RESIDUAL VECTOR FOR THE PROBLEM */
1773 /* DEFINED BY THE J-TH COLUMN VECTOR OF THE ARRAY B. */
1774 /* H(), G() ARRAYS OF WORKING SPACE OF LENGTH >= N. */
1775 /* IP() INT ARRAY OF WORKING SPACE OF LENGTH >= N */
1776 /* RECORDING PERMUTATION INDICES OF COLUMN VECTORS */
1777 /* Parameter adjustments */
1778 --ip;
1779 --g;
1780 --h__;
1781 a_dim1 = *mda;
1782 a_offset = 1 + a_dim1;
1783 a -= a_offset;
1784 --rnorm;
1785 b_dim1 = *mdb;
1786 b_offset = 1 + b_dim1;
1787 b -= b_offset;
1788
1789 /* Function Body */
1790 k = 0;
1791 ldiag = std::min(*m, *n);
1792 if (ldiag <= 0) {
1793 goto L270;
1794 }
1795 /* COMPUTE LMAX */
1796 i__1 = ldiag;
1797 for (j = 1; j <= i__1; ++j) {
1798 if (j == 1) {
1799 goto L20;
1800 }
1801 lmax = j;
1802 i__2 = *n;
1803 for (l = j; l <= i__2; ++l) {
1804 /* Computing 2nd power */
1805 d__1 = a[j - 1 + l * a_dim1];
1806 h__[l] -= d__1 * d__1;
1807 /* L10: */
1808 if (h__[l] > h__[lmax]) {
1809 lmax = l;
1810 }
1811 }
1812 d__1 = hmax + factor * h__[lmax];
1813 if (d__1 - hmax > zero) {
1814 goto L50;
1815 }
1816 L20:
1817 lmax = j;
1818 i__2 = *n;
1819 for (l = j; l <= i__2; ++l) {
1820 h__[l] = zero;
1821 i__3 = *m;
1822 for (i__ = j; i__ <= i__3; ++i__) {
1823 /* L30: */
1824 /* Computing 2nd power */
1825 d__1 = a[i__ + l * a_dim1];
1826 h__[l] += d__1 * d__1;
1827 }
1828 /* L40: */
1829 if (h__[l] > h__[lmax]) {
1830 lmax = l;
1831 }
1832 }
1833 hmax = h__[lmax];
1834 /* COLUMN INTERCHANGES IF NEEDED */
1835 L50:
1836 ip[j] = lmax;
1837 if (ip[j] == j) {
1838 goto L70;
1839 }
1840 i__2 = *m;
1841 for (i__ = 1; i__ <= i__2; ++i__) {
1842 tmp = a[i__ + j * a_dim1];
1843 a[i__ + j * a_dim1] = a[i__ + lmax * a_dim1];
1844 /* L60: */
1845 a[i__ + lmax * a_dim1] = tmp;
1846 }
1847 h__[lmax] = h__[j];
1848 /* J-TH TRANSFORMATION AND APPLICATION TO A AND B */
1849 L70:
1850 /* Computing MIN */
1851 i__2 = j + 1;
1852 i__ = std::min(i__2, *n);
1853 i__2 = j + 1;
1854 i__3 = *n - j;
1855 h12_(&c__1, &j, &i__2, m, &a[j * a_dim1 + 1], &c__1, &h__[j], &a[i__ * a_dim1 + 1], &c__1, mda, &i__3);
1856 /* L80: */
1857 i__2 = j + 1;
1858 h12_(&c__2, &j, &i__2, m, &a[j * a_dim1 + 1], &c__1, &h__[j], &b[b_offset], &c__1, mdb, nb);
1859 }
1860 /* DETERMINE PSEUDORANK */
1861 i__2 = ldiag;
1862 for (j = 1; j <= i__2; ++j) {
1863 /* L90: */
1864 if ((d__1 = a[j + j * a_dim1], std::abs(d__1)) <= *tau) {
1865 goto L100;
1866 }
1867 }
1868 k = ldiag;
1869 goto L110;
1870L100:
1871 k = j - 1;
1872L110:
1873 kp1 = k + 1;
1874 /* NORM OF RESIDUALS */
1875 i__2 = *nb;
1876 for (jb = 1; jb <= i__2; ++jb) {
1877 /* L130: */
1878 i__1 = *m - k;
1879 rnorm[jb] = dnrm2___(&i__1, &b[kp1 + jb * b_dim1], &c__1);
1880 }
1881 if (k > 0) {
1882 goto L160;
1883 }
1884 i__1 = *nb;
1885 for (jb = 1; jb <= i__1; ++jb) {
1886 i__2 = *n;
1887 for (i__ = 1; i__ <= i__2; ++i__) {
1888 /* L150: */
1889 b[i__ + jb * b_dim1] = zero;
1890 }
1891 }
1892 goto L270;
1893L160:
1894 if (k == *n) {
1895 goto L180;
1896 }
1897 /* HOUSEHOLDER DECOMPOSITION OF FIRST K ROWS */
1898 for (i__ = k; i__ >= 1; --i__) {
1899 /* L170: */
1900 i__2 = i__ - 1;
1901 h12_(&c__1, &i__, &kp1, n, &a[i__ + a_dim1], mda, &g[i__], &a[a_offset], mda, &c__1, &i__2);
1902 }
1903L180:
1904 i__2 = *nb;
1905 for (jb = 1; jb <= i__2; ++jb) {
1906 /* SOLVE K*K TRIANGULAR SYSTEM */
1907 for (i__ = k; i__ >= 1; --i__) {
1908 /* Computing MIN */
1909 i__1 = i__ + 1;
1910 j = std::min(i__1, *n);
1911 /* L210: */
1912 i__1 = k - i__;
1913 b[i__ + jb * b_dim1] =
1914 (b[i__ + jb * b_dim1] - ddot_sl__(&i__1, &a[i__ + j * a_dim1], mda, &b[j + jb * b_dim1], &c__1)) /
1915 a[i__ + i__ * a_dim1];
1916 }
1917 /* COMPLETE SOLUTION VECTOR */
1918 if (k == *n) {
1919 goto L240;
1920 }
1921 i__1 = *n;
1922 for (j = kp1; j <= i__1; ++j) {
1923 /* L220: */
1924 b[j + jb * b_dim1] = zero;
1925 }
1926 i__1 = k;
1927 for (i__ = 1; i__ <= i__1; ++i__) {
1928 /* L230: */
1929 h12_(&c__2, &i__, &kp1, n, &a[i__ + a_dim1], mda, &g[i__], &b[jb * b_dim1 + 1], &c__1, mdb, &c__1);
1930 }
1931 /* REORDER SOLUTION ACCORDING TO PREVIOUS COLUMN INTERCHANGES */
1932 L240:
1933 for (j = ldiag; j >= 1; --j) {
1934 if (ip[j] == j) {
1935 goto L250;
1936 }
1937 l = ip[j];
1938 tmp = b[l + jb * b_dim1];
1939 b[l + jb * b_dim1] = b[j + jb * b_dim1];
1940 b[j + jb * b_dim1] = tmp;
1941 L250:;
1942 }
1943 }
1944L270:
1945 *krank = k;
1946 return 0;
1947} /* hfti_ */
1948
1949int h12_(const int *mode, const int *lpivot, const int *l1, const int *m, double *u, const int *iue, double *up,
1950 double *c__, const int *ice, const int *icv, const int *ncv) {
1951 /* Initialized data */
1952
1953 static double one = 1.;
1954 static double zero = 0.;
1955
1956 /* System generated locals */
1957 int u_dim1, u_offset, i__1, i__2;
1958 double d__1;
1959
1960 /* Local variables */
1961 static double b;
1962 static int i__, j, i2, i3, i4;
1963 static double cl, sm;
1964 static int incr;
1965 static double clinv;
1966
1967 /* C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 */
1968 /* TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974 */
1969 /* CONSTRUCTION AND/OR APPLICATION OF A SINGLE */
1970 /* HOUSEHOLDER TRANSFORMATION Q = I + U*(U**T)/B */
1971 /* MODE = 1 OR 2 TO SELECT ALGORITHM H1 OR H2 . */
1972 /* LPIVOT IS THE INDEX OF THE PIVOT ELEMENT. */
1973 /* L1,M IF L1 <= M THE TRANSFORMATION WILL BE CONSTRUCTED TO */
1974 /* ZERO ELEMENTS INDEXED FROM L1 THROUGH M. */
1975 /* IF L1 > M THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION. */
1976 /* U(),IUE,UP */
1977 /* ON ENTRY TO H1 U() STORES THE PIVOT VECTOR. */
1978 /* IUE IS THE STORAGE INCREMENT BETWEEN ELEMENTS. */
1979 /* ON EXIT FROM H1 U() AND UP STORE QUANTITIES DEFINING */
1980 /* THE VECTOR U OF THE HOUSEHOLDER TRANSFORMATION. */
1981 /* ON ENTRY TO H2 U() AND UP */
1982 /* SHOULD STORE QUANTITIES PREVIOUSLY COMPUTED BY H1. */
1983 /* THESE WILL NOT BE MODIFIED BY H2. */
1984 /* C() ON ENTRY TO H1 OR H2 C() STORES A MATRIX WHICH WILL BE */
1985 /* REGARDED AS A SET OF VECTORS TO WHICH THE HOUSEHOLDER */
1986 /* TRANSFORMATION IS TO BE APPLIED. */
1987 /* ON EXIT C() STORES THE SET OF TRANSFORMED VECTORS. */
1988 /* ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C(). */
1989 /* ICV STORAGE INCREMENT BETWEEN VECTORS IN C(). */
1990 /* NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. */
1991 /* IF NCV <= 0 NO OPERATIONS WILL BE DONE ON C(). */
1992 /* Parameter adjustments */
1993 u_dim1 = *iue;
1994 u_offset = 1 + u_dim1;
1995 u -= u_offset;
1996 --c__;
1997
1998 /* Function Body */
1999 if (0 >= *lpivot || *lpivot >= *l1 || *l1 > *m) {
2000 goto L80;
2001 }
2002 cl = (d__1 = u[*lpivot * u_dim1 + 1], std::abs(d__1));
2003 if (*mode == 2) {
2004 goto L30;
2005 }
2006 /* ****** CONSTRUCT THE TRANSFORMATION ****** */
2007 i__1 = *m;
2008 for (j = *l1; j <= i__1; ++j) {
2009 sm = (d__1 = u[j * u_dim1 + 1], std::abs(d__1));
2010 /* L10: */
2011 cl = std::max(sm, cl);
2012 }
2013 if (cl <= zero) {
2014 goto L80;
2015 }
2016 clinv = one / cl;
2017 /* Computing 2nd power */
2018 d__1 = u[*lpivot * u_dim1 + 1] * clinv;
2019 sm = d__1 * d__1;
2020 i__1 = *m;
2021 for (j = *l1; j <= i__1; ++j) {
2022 /* L20: */
2023 /* Computing 2nd power */
2024 d__1 = u[j * u_dim1 + 1] * clinv;
2025 sm += d__1 * d__1;
2026 }
2027 cl *= sqrt(sm);
2028 if (u[*lpivot * u_dim1 + 1] > zero) {
2029 cl = -cl;
2030 }
2031 *up = u[*lpivot * u_dim1 + 1] - cl;
2032 u[*lpivot * u_dim1 + 1] = cl;
2033 goto L40;
2034/* ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C ****** */
2035L30:
2036 if (cl <= zero) {
2037 goto L80;
2038 }
2039L40:
2040 if (*ncv <= 0) {
2041 goto L80;
2042 }
2043 b = *up * u[*lpivot * u_dim1 + 1];
2044 if (b >= zero) {
2045 goto L80;
2046 }
2047 b = one / b;
2048 i2 = 1 - *icv + *ice * (*lpivot - 1);
2049 incr = *ice * (*l1 - *lpivot);
2050 i__1 = *ncv;
2051 for (j = 1; j <= i__1; ++j) {
2052 i2 += *icv;
2053 i3 = i2 + incr;
2054 i4 = i3;
2055 sm = c__[i2] * *up;
2056 i__2 = *m;
2057 for (i__ = *l1; i__ <= i__2; ++i__) {
2058 sm += c__[i3] * u[i__ * u_dim1 + 1];
2059 /* L50: */
2060 i3 += *ice;
2061 }
2062 if (sm == zero) {
2063 goto L70;
2064 }
2065 sm *= b;
2066 c__[i2] += sm * *up;
2067 i__2 = *m;
2068 for (i__ = *l1; i__ <= i__2; ++i__) {
2069 c__[i4] += sm * u[i__ * u_dim1 + 1];
2070 /* L60: */
2071 i4 += *ice;
2072 }
2073 L70:;
2074 }
2075L80:
2076 return 0;
2077} /* h12_ */
2078
2079int ldl_(const int *n, double *a, double *z__, const double *sigma, double *w) {
2080 /* Initialized data */
2081
2082 static double zero = 0.;
2083 static double one = 1.;
2084 static double four = 4.;
2085 static double epmach = 2.22e-16;
2086
2087 /* System generated locals */
2088 int i__1, i__2;
2089
2090 /* Local variables */
2091 static int i__, j;
2092 static double t, u, v;
2093 static int ij;
2094 static double tp, beta, gamma, alpha, delta;
2095
2096 /* LDL LDL' - RANK-ONE - UPDATE */
2097 /* PURPOSE: */
2098 /* UPDATES THE LDL' FACTORS OF MATRIX A BY RANK-ONE MATRIX */
2099 /* SIGMA*Z*Z' */
2100 /* INPUT ARGUMENTS: (* MEANS PARAMETERS ARE CHANGED DURING EXECUTION) */
2101 /* N : ORDER OF THE COEFFICIENT MATRIX A */
2102 /* * A : POSITIVE DEFINITE MATRIX OF DIMENSION N; */
2103 /* ONLY THE LOWER TRIANGLE IS USED AND IS STORED COLUMN BY */
2104 /* COLUMN AS ONE DIMENSIONAL ARRAY OF DIMENSION N*(N+1)/2. */
2105 /* * Z : VECTOR OF DIMENSION N OF UPDATING ELEMENTS */
2106 /* SIGMA : SCALAR FACTOR BY WHICH THE MODIFYING DYADE Z*Z' IS */
2107 /* MULTIPLIED */
2108 /* OUTPUT ARGUMENTS: */
2109 /* A : UPDATED LDL' FACTORS */
2110 /* WORKING ARRAY: */
2111 /* W : VECTOR OP DIMENSION N (USED ONLY IF SIGMA .LT. ZERO) */
2112 /* METHOD: */
2113 /* THAT OF FLETCHER AND POWELL AS DESCRIBED IN : */
2114 /* FLETCHER,R.,(1974) ON THE MODIFICATION OF LDL' FACTORIZATION. */
2115 /* POWELL,M.J.D. MATH.COMPUTATION 28, 1067-1078. */
2116 /* IMPLEMENTED BY: */
2117 /* KRAFT,D., DFVLR - INSTITUT FUER DYNAMIK DER FLUGSYSTEME */
2118 /* D-8031 OBERPFAFFENHOFEN */
2119 /* STATUS: 15. JANUARY 1980 */
2120 /* SUBROUTINES REQUIRED: NONE */
2121 /* Parameter adjustments */
2122 --w;
2123 --z__;
2124 --a;
2125
2126 /* Function Body */
2127 if (*sigma == zero) {
2128 goto L280;
2129 }
2130 ij = 1;
2131 t = one / *sigma;
2132 if (*sigma > zero) {
2133 goto L220;
2134 }
2135 /* PREPARE NEGATIVE UPDATE */
2136 i__1 = *n;
2137 for (i__ = 1; i__ <= i__1; ++i__) {
2138 /* L150: */
2139 w[i__] = z__[i__];
2140 }
2141 i__1 = *n;
2142 for (i__ = 1; i__ <= i__1; ++i__) {
2143 v = w[i__];
2144 t += v * v / a[ij];
2145 i__2 = *n;
2146 for (j = i__ + 1; j <= i__2; ++j) {
2147 ++ij;
2148 /* L160: */
2149 w[j] -= v * a[ij];
2150 }
2151 /* L170: */
2152 ++ij;
2153 }
2154 if (t >= zero) {
2155 t = epmach / *sigma;
2156 }
2157 i__1 = *n;
2158 for (i__ = 1; i__ <= i__1; ++i__) {
2159 j = *n + 1 - i__;
2160 ij -= i__;
2161 u = w[j];
2162 w[j] = t;
2163 /* L210: */
2164 t -= u * u / a[ij];
2165 }
2166L220:
2167 /* HERE UPDATING BEGINS */
2168 i__1 = *n;
2169 for (i__ = 1; i__ <= i__1; ++i__) {
2170 v = z__[i__];
2171 delta = v / a[ij];
2172 if (*sigma < zero) {
2173 tp = w[i__];
2174 }
2175 if (*sigma > zero) {
2176 tp = t + delta * v;
2177 }
2178 alpha = tp / t;
2179 a[ij] = alpha * a[ij];
2180 if (i__ == *n) {
2181 goto L280;
2182 }
2183 beta = delta / tp;
2184 if (alpha > four) {
2185 goto L240;
2186 }
2187 i__2 = *n;
2188 for (j = i__ + 1; j <= i__2; ++j) {
2189 ++ij;
2190 z__[j] -= v * a[ij];
2191 /* L230: */
2192 a[ij] += beta * z__[j];
2193 }
2194 goto L260;
2195 L240:
2196 gamma = t / tp;
2197 i__2 = *n;
2198 for (j = i__ + 1; j <= i__2; ++j) {
2199 ++ij;
2200 u = a[ij];
2201 a[ij] = gamma * u + beta * z__[j];
2202 /* L250: */
2203 z__[j] -= v * u;
2204 }
2205 L260:
2206 ++ij;
2207 /* L270: */
2208 t = tp;
2209 }
2210L280:
2211 return 0;
2212 /* END OF LDL */
2213} /* ldl_ */
2214
2215double linmin_(int *mode, const double *ax, const double *bx, const double *f, const double *tol) {
2216 /* Initialized data */
2217
2218 static double c__ = .381966011;
2219 static double eps = 1.5e-8;
2220 static double zero = 0.;
2221
2222 /* System generated locals */
2223 double ret_val, d__1;
2224
2225 /* Local variables */
2226 static double a, b, d__, e, m, p, q, r__, u, v, w, x, fu, fv, fw, fx, tol1, tol2;
2227
2228 /* LINMIN LINESEARCH WITHOUT DERIVATIVES */
2229 /* PURPOSE: */
2230 /* TO FIND THE ARGUMENT LINMIN WHERE THE FUNCTION F TAKES IT'S MINIMUM */
2231 /* ON THE INTERVAL AX, BX. */
2232 /* COMBINATION OF GOLDEN SECTION AND SUCCESSIVE QUADRATIC INTERPOLATION. */
2233 /* INPUT ARGUMENTS: (* MEANS PARAMETERS ARE CHANGED DURING EXECUTION) */
2234 /* *MODE SEE OUTPUT ARGUMENTS */
2235 /* AX LEFT ENDPOINT OF INITIAL INTERVAL */
2236 /* BX RIGHT ENDPOINT OF INITIAL INTERVAL */
2237 /* F FUNCTION VALUE AT LINMIN WHICH IS TO BE BROUGHT IN BY */
2238 /* REVERSE COMMUNICATION CONTROLLED BY MODE */
2239 /* TOL DESIRED LENGTH OF INTERVAL OF UNCERTAINTY OF FINAL RESULT */
2240 /* OUTPUT ARGUMENTS: */
2241 /* LINMIN ABSCISSA APPROXIMATING THE POINT WHERE F ATTAINS A MINIMUM */
2242 /* MODE CONTROLS REVERSE COMMUNICATION */
2243 /* MUST BE SET TO 0 INITIALLY, RETURNS WITH INTERMEDIATE */
2244 /* VALUES 1 AND 2 WHICH MUST NOT BE CHANGED BY THE USER, */
2245 /* ENDS WITH CONVERGENCE WITH VALUE 3. */
2246 /* WORKING ARRAY: */
2247 /* NONE */
2248 /* METHOD: */
2249 /* THIS FUNCTION SUBPROGRAM IS A SLIGHTLY MODIFIED VERSION OF THE */
2250 /* ALGOL 60 PROCEDURE LOCALMIN GIVEN IN */
2251 /* R.P. BRENT: ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES, */
2252 /* PRENTICE-HALL (1973). */
2253 /* IMPLEMENTED BY: */
2254 /* KRAFT, D., DFVLR - INSTITUT FUER DYNAMIK DER FLUGSYSTEME */
2255 /* D-8031 OBERPFAFFENHOFEN */
2256 /* STATUS: 31. AUGUST 1984 */
2257 /* SUBROUTINES REQUIRED: NONE */
2258 /* EPS = SQUARE - ROOT OF MACHINE PRECISION */
2259 /* C = GOLDEN SECTION RATIO = (3-SQRT(5))/2 */
2260 switch (*mode) {
2261 case 1:
2262 goto L10;
2263 case 2:
2264 goto L55;
2265 }
2266 /* INITIALIZATION */
2267 a = *ax;
2268 b = *bx;
2269 e = zero;
2270 v = a + c__ * (b - a);
2271 w = v;
2272 x = w;
2273 ret_val = x;
2274 *mode = 1;
2275 goto L100;
2276/* MAIN LOOP STARTS HERE */
2277L10:
2278 fx = *f;
2279 fv = fx;
2280 fw = fv;
2281L20:
2282 m = (a + b) * .5;
2283 tol1 = eps * std::abs(x) + *tol;
2284 tol2 = tol1 + tol1;
2285 /* TEST CONVERGENCE */
2286 if ((d__1 = x - m, std::abs(d__1)) <= tol2 - (b - a) * .5) {
2287 goto L90;
2288 }
2289 r__ = zero;
2290 q = r__;
2291 p = q;
2292 if (std::abs(e) <= tol1) {
2293 goto L30;
2294 }
2295 /* FIT PARABOLA */
2296 r__ = (x - w) * (fx - fv);
2297 q = (x - v) * (fx - fw);
2298 p = (x - v) * q - (x - w) * r__;
2299 q -= r__;
2300 q += q;
2301 if (q > zero) {
2302 p = -p;
2303 }
2304 if (q < zero) {
2305 q = -q;
2306 }
2307 r__ = e;
2308 e = d__;
2309/* IS PARABOLA ACCEPTABLE */
2310L30:
2311 if (std::abs(p) >= (d__1 = q * r__, std::abs(d__1)) * .5 || p <= q * (a - x) || p >= q * (b - x)) {
2312 goto L40;
2313 }
2314 /* PARABOLIC INTERPOLATION STEP */
2315 d__ = p / q;
2316 /* F MUST NOT BE EVALUATED TOO CLOSE TO A OR B */
2317 if (u - a < tol2) {
2318 d__1 = m - x;
2319 d__ = d_sign(&tol1, &d__1);
2320 }
2321 if (b - u < tol2) {
2322 d__1 = m - x;
2323 d__ = d_sign(&tol1, &d__1);
2324 }
2325 goto L50;
2326/* GOLDEN SECTION STEP */
2327L40:
2328 if (x >= m) {
2329 e = a - x;
2330 }
2331 if (x < m) {
2332 e = b - x;
2333 }
2334 d__ = c__ * e;
2335/* F MUST NOT BE EVALUATED TOO CLOSE TO X */
2336L50:
2337 if (std::abs(d__) < tol1) {
2338 d__ = d_sign(&tol1, &d__);
2339 }
2340 u = x + d__;
2341 ret_val = u;
2342 *mode = 2;
2343 goto L100;
2344L55:
2345 fu = *f;
2346 /* UPDATE A, B, V, W, AND X */
2347 if (fu > fx) {
2348 goto L60;
2349 }
2350 if (u >= x) {
2351 a = x;
2352 }
2353 if (u < x) {
2354 b = x;
2355 }
2356 v = w;
2357 fv = fw;
2358 w = x;
2359 fw = fx;
2360 x = u;
2361 fx = fu;
2362 goto L85;
2363L60:
2364 if (u < x) {
2365 a = u;
2366 }
2367 if (u >= x) {
2368 b = u;
2369 }
2370 if (fu <= fw || w == x) {
2371 goto L70;
2372 }
2373 if (fu <= fv || v == x || v == w) {
2374 goto L80;
2375 }
2376 goto L85;
2377L70:
2378 v = w;
2379 fv = fw;
2380 w = u;
2381 fw = fu;
2382 goto L85;
2383L80:
2384 v = u;
2385 fv = fu;
2386L85:
2387 goto L20;
2388/* END OF MAIN LOOP */
2389L90:
2390 ret_val = x;
2391 *mode = 3;
2392L100:
2393 return ret_val;
2394 /* END OF LINMIN */
2395} /* linmin_ */
2396
2397/* ## Following a selection from BLAS Level 1 */
2398int daxpy_sl__(const int *n, const double *da, double *dx, const int *incx, double *dy, const int *incy) {
2399 /* System generated locals */
2400 int i__1;
2401
2402 /* Local variables */
2403 static int i__, m, ix, iy, mp1;
2404
2405 /* CONSTANT TIMES A VECTOR PLUS A VECTOR. */
2406 /* USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. */
2407 /* JACK DONGARRA, LINPACK, 3/11/78. */
2408 /* Parameter adjustments */
2409 --dy;
2410 --dx;
2411
2412 /* Function Body */
2413 if (*n <= 0) {
2414 return 0;
2415 }
2416 if (*da == 0.) {
2417 return 0;
2418 }
2419 if (*incx == 1 && *incy == 1) {
2420 goto L20;
2421 }
2422 /* CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS */
2423 /* NOT EQUAL TO 1 */
2424 ix = 1;
2425 iy = 1;
2426 if (*incx < 0) {
2427 ix = (-(*n) + 1) * *incx + 1;
2428 }
2429 if (*incy < 0) {
2430 iy = (-(*n) + 1) * *incy + 1;
2431 }
2432 i__1 = *n;
2433 for (i__ = 1; i__ <= i__1; ++i__) {
2434 dy[iy] += *da * dx[ix];
2435 ix += *incx;
2436 iy += *incy;
2437 /* L10: */
2438 }
2439 return 0;
2440/* CODE FOR BOTH INCREMENTS EQUAL TO 1 */
2441/* CLEAN-UP LOOP */
2442L20:
2443 m = *n % 4;
2444 if (m == 0) {
2445 goto L40;
2446 }
2447 i__1 = m;
2448 for (i__ = 1; i__ <= i__1; ++i__) {
2449 dy[i__] += *da * dx[i__];
2450 /* L30: */
2451 }
2452 if (*n < 4) {
2453 return 0;
2454 }
2455L40:
2456 mp1 = m + 1;
2457 i__1 = *n;
2458 for (i__ = mp1; i__ <= i__1; i__ += 4) {
2459 dy[i__] += *da * dx[i__];
2460 dy[i__ + 1] += *da * dx[i__ + 1];
2461 dy[i__ + 2] += *da * dx[i__ + 2];
2462 dy[i__ + 3] += *da * dx[i__ + 3];
2463 /* L50: */
2464 }
2465 return 0;
2466} /* daxpy_sl__ */
2467
2468int dcopy___(const int *n, double *dx, const int *incx, double *dy, const int *incy) {
2469 /* System generated locals */
2470 int i__1;
2471
2472 /* Local variables */
2473 static int i__, m, ix, iy, mp1;
2474
2475 /* COPIES A VECTOR, X, TO A VECTOR, Y. */
2476 /* USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. */
2477 /* JACK DONGARRA, LINPACK, 3/11/78. */
2478 /* Parameter adjustments */
2479 --dy;
2480 --dx;
2481
2482 /* Function Body */
2483 if (*n <= 0) {
2484 return 0;
2485 }
2486 if (*incx == 1 && *incy == 1) {
2487 goto L20;
2488 }
2489 /* CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS */
2490 /* NOT EQUAL TO 1 */
2491 ix = 1;
2492 iy = 1;
2493 if (*incx < 0) {
2494 ix = (-(*n) + 1) * *incx + 1;
2495 }
2496 if (*incy < 0) {
2497 iy = (-(*n) + 1) * *incy + 1;
2498 }
2499 i__1 = *n;
2500 for (i__ = 1; i__ <= i__1; ++i__) {
2501 dy[iy] = dx[ix];
2502 ix += *incx;
2503 iy += *incy;
2504 /* L10: */
2505 }
2506 return 0;
2507/* CODE FOR BOTH INCREMENTS EQUAL TO 1 */
2508/* CLEAN-UP LOOP */
2509L20:
2510 m = *n % 7;
2511 if (m == 0) {
2512 goto L40;
2513 }
2514 i__1 = m;
2515 for (i__ = 1; i__ <= i__1; ++i__) {
2516 dy[i__] = dx[i__];
2517 /* L30: */
2518 }
2519 if (*n < 7) {
2520 return 0;
2521 }
2522L40:
2523 mp1 = m + 1;
2524 i__1 = *n;
2525 for (i__ = mp1; i__ <= i__1; i__ += 7) {
2526 dy[i__] = dx[i__];
2527 dy[i__ + 1] = dx[i__ + 1];
2528 dy[i__ + 2] = dx[i__ + 2];
2529 dy[i__ + 3] = dx[i__ + 3];
2530 dy[i__ + 4] = dx[i__ + 4];
2531 dy[i__ + 5] = dx[i__ + 5];
2532 dy[i__ + 6] = dx[i__ + 6];
2533 /* L50: */
2534 }
2535 return 0;
2536} /* dcopy___ */
2537
2538double ddot_sl__(const int *n, double *dx, const int *incx, double *dy, const int *incy) {
2539 /* System generated locals */
2540 int i__1;
2541 double ret_val;
2542
2543 /* Local variables */
2544 static int i__, m, ix, iy, mp1;
2545 static double dtemp;
2546
2547 /* FORMS THE DOT PRODUCT OF TWO VECTORS. */
2548 /* USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. */
2549 /* JACK DONGARRA, LINPACK, 3/11/78. */
2550 /* Parameter adjustments */
2551 --dy;
2552 --dx;
2553
2554 /* Function Body */
2555 ret_val = 0.;
2556 dtemp = 0.;
2557 if (*n <= 0) {
2558 return ret_val;
2559 }
2560 if (*incx == 1 && *incy == 1) {
2561 goto L20;
2562 }
2563 /* CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS */
2564 /* NOT EQUAL TO 1 */
2565 ix = 1;
2566 iy = 1;
2567 if (*incx < 0) {
2568 ix = (-(*n) + 1) * *incx + 1;
2569 }
2570 if (*incy < 0) {
2571 iy = (-(*n) + 1) * *incy + 1;
2572 }
2573 i__1 = *n;
2574 for (i__ = 1; i__ <= i__1; ++i__) {
2575 dtemp += dx[ix] * dy[iy];
2576 ix += *incx;
2577 iy += *incy;
2578 /* L10: */
2579 }
2580 ret_val = dtemp;
2581 return ret_val;
2582/* CODE FOR BOTH INCREMENTS EQUAL TO 1 */
2583/* CLEAN-UP LOOP */
2584L20:
2585 m = *n % 5;
2586 if (m == 0) {
2587 goto L40;
2588 }
2589 i__1 = m;
2590 for (i__ = 1; i__ <= i__1; ++i__) {
2591 dtemp += dx[i__] * dy[i__];
2592 /* L30: */
2593 }
2594 if (*n < 5) {
2595 goto L60;
2596 }
2597L40:
2598 mp1 = m + 1;
2599 i__1 = *n;
2600 for (i__ = mp1; i__ <= i__1; i__ += 5) {
2601 dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[i__ + 2] * dy[i__ + 2] +
2602 dx[i__ + 3] * dy[i__ + 3] + dx[i__ + 4] * dy[i__ + 4];
2603 /* L50: */
2604 }
2605L60:
2606 ret_val = dtemp;
2607 return ret_val;
2608} /* ddot_sl__ */
2609
2610double dnrm2___(const int *n, double *dx, const int *incx) {
2611 /* Initialized data */
2612
2613 static double zero = 0.;
2614 static double one = 1.;
2615 static double cutlo = 8.232e-11;
2616 static double cuthi = 1.304e19;
2617
2618 /* Format strings */
2619 /* -- Unused -- */
2620 /* static char fmt_30[] = "";
2621 static char fmt_50[] = "";
2622 static char fmt_70[] = "";
2623 static char fmt_110[] = ""; */
2624
2625 /* System generated locals */
2626 int i__1, i__2;
2627 double ret_val, d__1;
2628
2629 /* Local variables */
2630 static int i__, j, nn;
2631 static double sum, xmax;
2632 static int next;
2633 static double hitest;
2634
2635 /* Assigned format variables */
2636 /* static char *next_fmt = fmt_30; */ /* Avoid compiler warning */
2637
2638 /* Parameter adjustments */
2639 --dx;
2640
2641 /* Function Body */
2642 /* EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE */
2643 /* INCREMENT INCX . */
2644 /* IF N .LE. 0 RETURN WITH RESULT = 0. */
2645 /* IF N .GE. 1 THEN INCX MUST BE .GE. 1 */
2646 /* C.L.LAWSON, 1978 JAN 08 */
2647 /* FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE */
2648 /* HOPEFULLY APPLICABLE TO ALL MACHINES. */
2649 /* CUTLO = MAXIMUM OF SQRT(U/EPS) OVER ALL KNOWN MACHINES. */
2650 /* CUTHI = MINIMUM OF SQRT(V) OVER ALL KNOWN MACHINES. */
2651 /* WHERE */
2652 /* EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1. */
2653 /* U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT) */
2654 /* V = LARGEST NO. (OVERFLOW LIMIT) */
2655 /* BRIEF OUTLINE OF ALGORITHM.. */
2656 /* PHASE 1 SCANS ZERO COMPONENTS. */
2657 /* MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO */
2658 /* MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO */
2659 /* MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M */
2660 /* WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX. */
2661 /* VALUES FOR CUTLO AND CUTHI.. */
2662 /* FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER */
2663 /* DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS.. */
2664 /* CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE */
2665 /* UNIVAC AND DEC AT 2**(-103) */
2666 /* THUS CUTLO = 2**(-51) = 4.44089E-16 */
2667 /* CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC. */
2668 /* THUS CUTHI = 2**(63.5) = 1.30438E19 */
2669 /* CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC. */
2670 /* THUS CUTLO = 2**(-33.5) = 8.23181D-11 */
2671 /* CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19 */
2672 /* DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / */
2673 /* DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / */
2674 if (*n > 0) {
2675 goto L10;
2676 }
2677 ret_val = zero;
2678 goto L300;
2679L10:
2680 next = 0;
2681 /* next_fmt = fmt_30; */
2682 sum = zero;
2683 nn = *n * *incx;
2684 /* BEGIN MAIN LOOP */
2685 i__ = 1;
2686L20:
2687 switch (next) {
2688 case 0:
2689 goto L30;
2690 case 1:
2691 goto L50;
2692 case 2:
2693 goto L70;
2694 case 3:
2695 goto L110;
2696 }
2697L30:
2698 if ((d__1 = dx[i__], std::abs(d__1)) > cutlo) {
2699 goto L85;
2700 }
2701 next = 1;
2702 /* next_fmt = fmt_50 */;
2703 xmax = zero;
2704/* PHASE 1. SUM IS ZERO */
2705L50:
2706 if (dx[i__] == zero) {
2707 goto L200;
2708 }
2709 if ((d__1 = dx[i__], std::abs(d__1)) > cutlo) {
2710 goto L85;
2711 }
2712 /* PREPARE FOR PHASE 2. */
2713 next = 2;
2714 /* next_fmt = fmt_70; */
2715 goto L105;
2716/* PREPARE FOR PHASE 4. */
2717L100:
2718 i__ = j;
2719 next = 3;
2720 /* next_fmt = fmt_110; */
2721 sum = sum / dx[i__] / dx[i__];
2722L105:
2723 xmax = (d__1 = dx[i__], std::abs(d__1));
2724 goto L115;
2725/* PHASE 2. SUM IS SMALL. */
2726/* SCALE TO AVOID DESTRUCTIVE UNDERFLOW. */
2727L70:
2728 if ((d__1 = dx[i__], std::abs(d__1)) > cutlo) {
2729 goto L75;
2730 }
2731/* COMMON CODE FOR PHASES 2 AND 4. */
2732/* IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. */
2733L110:
2734 if ((d__1 = dx[i__], std::abs(d__1)) <= xmax) {
2735 goto L115;
2736 }
2737 /* Computing 2nd power */
2738 d__1 = xmax / dx[i__];
2739 sum = one + sum * (d__1 * d__1);
2740 xmax = (d__1 = dx[i__], std::abs(d__1));
2741 goto L200;
2742L115:
2743 /* Computing 2nd power */
2744 d__1 = dx[i__] / xmax;
2745 sum += d__1 * d__1;
2746 goto L200;
2747/* PREPARE FOR PHASE 3. */
2748L75:
2749 sum = sum * xmax * xmax;
2750/* FOR REAL OR D.P. SET HITEST = CUTHI/N */
2751/* FOR COMPLEX SET HITEST = CUTHI/(2*N) */
2752L85:
2753 hitest = cuthi / static_cast<float>(*n);
2754 /* PHASE 3. SUM IS MID-RANGE. NO SCALING. */
2755 i__1 = nn;
2756 i__2 = *incx;
2757 for (j = i__; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
2758 if ((d__1 = dx[j], std::abs(d__1)) >= hitest) {
2759 goto L100;
2760 }
2761 /* L95: */
2762 /* Computing 2nd power */
2763 d__1 = dx[j];
2764 sum += d__1 * d__1;
2765 }
2766 ret_val = sqrt(sum);
2767 goto L300;
2768L200:
2769 i__ += *incx;
2770 if (i__ <= nn) {
2771 goto L20;
2772 }
2773 /* END OF MAIN LOOP. */
2774 /* COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. */
2775 ret_val = xmax * sqrt(sum);
2776L300:
2777 return ret_val;
2778} /* dnrm2___ */
2779
2780int dsrot_(const int *n, double *dx, const int *incx, double *dy, const int *incy, const double *c__, const double *s) {
2781 /* System generated locals */
2782 int i__1;
2783
2784 /* Local variables */
2785 static int i__, ix, iy;
2786 static double dtemp;
2787
2788 /* APPLIES A PLANE ROTATION. */
2789 /* JACK DONGARRA, LINPACK, 3/11/78. */
2790 /* Parameter adjustments */
2791 --dy;
2792 --dx;
2793
2794 /* Function Body */
2795 if (*n <= 0) {
2796 return 0;
2797 }
2798 if (*incx == 1 && *incy == 1) {
2799 goto L20;
2800 }
2801 /* CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL */
2802 /* TO 1 */
2803 ix = 1;
2804 iy = 1;
2805 if (*incx < 0) {
2806 ix = (-(*n) + 1) * *incx + 1;
2807 }
2808 if (*incy < 0) {
2809 iy = (-(*n) + 1) * *incy + 1;
2810 }
2811 i__1 = *n;
2812 for (i__ = 1; i__ <= i__1; ++i__) {
2813 dtemp = *c__ * dx[ix] + *s * dy[iy];
2814 dy[iy] = *c__ * dy[iy] - *s * dx[ix];
2815 dx[ix] = dtemp;
2816 ix += *incx;
2817 iy += *incy;
2818 /* L10: */
2819 }
2820 return 0;
2821/* CODE FOR BOTH INCREMENTS EQUAL TO 1 */
2822L20:
2823 i__1 = *n;
2824 for (i__ = 1; i__ <= i__1; ++i__) {
2825 dtemp = *c__ * dx[i__] + *s * dy[i__];
2826 dy[i__] = *c__ * dy[i__] - *s * dx[i__];
2827 dx[i__] = dtemp;
2828 /* L30: */
2829 }
2830 return 0;
2831} /* dsrot_ */
2832
2833int dsrotg_(double *da, double *db, double *c__, double *s) {
2834 /* Initialized data */
2835
2836 static double one = 1.;
2837 static double zero = 0.;
2838
2839 /* System generated locals */
2840 double d__1, d__2;
2841
2842 /* Local variables */
2843 static double r__, z__, roe, scale;
2844
2845 /* CONSTRUCT GIVENS PLANE ROTATION. */
2846 /* JACK DONGARRA, LINPACK, 3/11/78. */
2847 /* MODIFIED 9/27/86. */
2848 roe = *db;
2849 if (std::abs(*da) > std::abs(*db)) {
2850 roe = *da;
2851 }
2852 scale = std::abs(*da) + std::abs(*db);
2853 if (scale != zero) {
2854 goto L10;
2855 }
2856 *c__ = one;
2857 *s = zero;
2858 r__ = zero;
2859 goto L20;
2860L10:
2861 /* Computing 2nd power */
2862 d__1 = *da / scale;
2863 /* Computing 2nd power */
2864 d__2 = *db / scale;
2865 r__ = scale * sqrt(d__1 * d__1 + d__2 * d__2);
2866 r__ = d_sign(&one, &roe) * r__;
2867 *c__ = *da / r__;
2868 *s = *db / r__;
2869L20:
2870 z__ = *s;
2871 if (std::abs(*c__) > zero && std::abs(*c__) <= *s) {
2872 z__ = one / *c__;
2873 }
2874 *da = r__;
2875 *db = z__;
2876 return 0;
2877} /* dsrotg_ */
2878
2879int dscal_sl__(const int *n, const double *da, double *dx, const int *incx) {
2880 /* System generated locals */
2881 int i__1, i__2;
2882
2883 /* Local variables */
2884 static int i__, m, mp1, nincx;
2885
2886 /* SCALES A VECTOR BY A CONSTANT. */
2887 /* USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. */
2888 /* JACK DONGARRA, LINPACK, 3/11/78. */
2889 /* Parameter adjustments */
2890 --dx;
2891
2892 /* Function Body */
2893 if (*n <= 0) {
2894 return 0;
2895 }
2896 if (*incx == 1) {
2897 goto L20;
2898 }
2899 /* CODE FOR INCREMENT NOT EQUAL TO 1 */
2900 nincx = *n * *incx;
2901 i__1 = nincx;
2902 i__2 = *incx;
2903 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
2904 dx[i__] = *da * dx[i__];
2905 /* L10: */
2906 }
2907 return 0;
2908/* CODE FOR INCREMENT EQUAL TO 1 */
2909/* CLEAN-UP LOOP */
2910L20:
2911 m = *n % 5;
2912 if (m == 0) {
2913 goto L40;
2914 }
2915 i__2 = m;
2916 for (i__ = 1; i__ <= i__2; ++i__) {
2917 dx[i__] = *da * dx[i__];
2918 /* L30: */
2919 }
2920 if (*n < 5) {
2921 return 0;
2922 }
2923L40:
2924 mp1 = m + 1;
2925 i__2 = *n;
2926 for (i__ = mp1; i__ <= i__2; i__ += 5) {
2927 dx[i__] = *da * dx[i__];
2928 dx[i__ + 1] = *da * dx[i__ + 1];
2929 dx[i__ + 2] = *da * dx[i__ + 2];
2930 dx[i__ + 3] = *da * dx[i__ + 3];
2931 dx[i__ + 4] = *da * dx[i__ + 4];
2932 /* L50: */
2933 }
2934 return 0;
2935} /* dscal_sl__ */
2936
2938
2939} // namespace
2940
2941} // namespace Mantid::Kernel::Math
tmp
gsl_vector * tmp
Definition: AugmentedLagrangianOptimizer.cpp:49
tp
static std::unique_ptr< QThreadPool > tp
Definition: FindFilesThreadPoolManager.cpp:27
value
double value
The value of the point.
Definition: FitMW.cpp:51
sigma
double sigma
Definition: GetAllEi.cpp:156
index
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
Mantid::Kernel::Math::SLSQPMinimizer::initializeConstraints
void initializeConstraints(const DblMatrix &equality, const DblMatrix &inequality)
Create constraint array.
Definition: SLSQPMinimizer.cpp:129
Mantid::Kernel::Math::SLSQPMinimizer::minimize
std::vector< double > minimize(const std::vector< double > &x0) const
Perform the minimization.
Definition: SLSQPMinimizer.cpp:23
Mantid::Kernel::Math::SLSQPMinimizer::m_nparams
const size_t m_nparams
Number of parameters under minimization.
Definition: SLSQPMinimizer.h:145
Mantid::Kernel::Math::SLSQPMinimizer::evaluateConstraints
void evaluateConstraints(std::vector< double > &constrValues, const std::vector< double > &x) const
Compute values of constraints.
Definition: SLSQPMinimizer.cpp:106
Mantid::Kernel::Math::SLSQPMinimizer::fprime
void fprime(std::vector< double > &grad, const std::vector< double > &x) const
Compute derivative numerically.
Definition: SLSQPMinimizer.cpp:85
Mantid::Kernel::Math::SLSQPMinimizer::fvalue
double fvalue(const std::vector< double > &x) const
Compute the value of the objective function.
Definition: SLSQPMinimizer.h:104
Mantid::Kernel::Math::SLSQPMinimizer::m_constraintNorms
std::vector< double > m_constraintNorms
Holder for constraint normals.
Definition: SLSQPMinimizer.h:154
Mantid::Kernel::Matrix::numRows
size_t numRows() const
Return the number of rows in the matrix.
Definition: Matrix.h:144
Mantid::Kernel::Matrix::numCols
size_t numCols() const
Return the number of columns in the matrix.
Definition: Matrix.h:147
Generated on Tue Jan 3 2023 11:57:51 for Mantid by doxygen 1.9.5