Login
[x]
Log in using an account from:
Fedora Account System
Red Hat Associate
Red Hat Customer
Or login using a Red Hat Bugzilla account
Forgot Password
Login:
Hide Forgot
Create an Account
Red Hat Bugzilla – Attachment 267611 Details for
Bug 396921
gcc 4.1.2 produces bad code
[?]
New
Simple Search
Advanced Search
My Links
Browse
Requests
Reports
Current State
Search
Tabular reports
Graphical reports
Duplicates
Other Reports
User Changes
Plotly Reports
Bug Status
Bug Severity
Non-Defaults
|
Product Dashboard
Help
Page Help!
Bug Writing Guidelines
What's new
Browser Support Policy
5.0.4.rh83 Release notes
FAQ
Guides index
User guide
Web Services
Contact
Legal
This site requires JavaScript to be enabled to function correctly, please enable it.
First demonstration file
f1.c (text/x-csrc), 110.21 KB, created by
Antti Honkela
on 2007-11-23 15:36:58 UTC
(
hide
)
Description:
First demonstration file
Filename:
MIME Type:
Creator:
Antti Honkela
Created:
2007-11-23 15:36:58 UTC
Size:
110.21 KB
patch
obsolete
>/* Conditions of use: */ >/* 1. donlp2 is under the exclusive copyright of P. Spellucci */ >/* (e-mail:spellucci@mathematik.tu-darmstadt.de) */ >/* "donlp2" is a reserved name */ >/* 2. donlp2 and its constituent parts come with no warranty, whether ex- */ >/* pressed or implied, that it is free of errors or suitable for any */ >/* specific purpose. */ >/* It must not be used to solve any problem, whose incorrect solution */ >/* could result in injury to a person , institution or property. */ >/* It is at the users own risk to use donlp2 or parts of it and the */ >/* author disclaims all liability for such use. */ >/* 3. donlp2 is distributed "as is". In particular, no maintenance, support */ >/* or trouble-shooting or subsequent upgrade is implied. */ >/* 4. The use of donlp2 must be acknowledged, in any publication which */ >/* contains */ >/* results obtained with it or parts of it. Citation of the authors name */ >/* and netlib-source is suitable. */ >/* 5. The free use of donlp2 and parts of it is restricted for research */ >/* purposes */ >/* commercial uses require permission and licensing from P. Spellucci. */ > >/* d o n l p 2 */ > >/* version 28/11/2001 (*) */ >/* tauqp dependent on scalres only */ >/* weights computed in a modified version in the singular case */ >/* some comparisons relative to zero changed */ >/* error-return for function evaluation as added feature */ >/* termination of QP-solver changed (d not set to zero) */ >/* new version of BFGS: if nr = 0 take Powell's update */ >/* no suppression of update beforehand (with exception dg = 0) */ >/* plus some minor corrections */ > >/* for consistency reasons variable names cf and cgf changed into */ >/* icf and icgf */ >/* added feature numerical differentiation of order 1,2,6 */ >/* requires new parameters epsfcn, taubnd, difftype */ >/* added feature of external blockwise evaluation of functions */ >/* rather than individual evaluation */ >/* requires new parameter bloc */ >/* new feature of user-supplied scaling of donlp2_x */ > >/* ************************************************************************* */ >/* (*) using a c-version of donlp2 as of 1998 made by */ >/* by S. Schoeffert, ASB, Bourges, France */ >/* dynamic memory allocation added by k.s. Cover */ >/* copyrights and right of commercial exploitation remaining by contract */ >/* with P. Spellucci */ >/* ************************************************************************ */ >/* new function format */ >/* developed in 2002/2003 by P. Spellucci copyrigth P. Spellucci */ >/* The problem format here is */ >/* min_x f(donlp2_x) */ >/* */ >/* subject to */ >/* low(i) <= donlp2_x(i) <= up(i) , i=1,..,n */ >/* low(i+n) <= (Ax)_i <= up(i+n) , i=1,nlin nlin=0 is allowed */ >/*low(i+nlin+n) <= c_i(donlp2_x) <= up(i+n+nlin) i=1,nonlin, nonlin=0 allowe */ >/* */ >/* */ > >/* low(i)=up(i) is allowed . In this case only one constraint is used */ >/* low(i)=-big, up(i)=big is allowed, big is userdefined */ >/* these are a total of 2*(n+nlin+nonlin) constraints, numbered */ >/* consecutively */ >/* such that odd numbers correspond to the lower and even numbered to the */ >/* upper */ >/* bound. A is to be stored in the nlin first columns of gres. */ >/* the user evaluation routines for the constraints only deal with the */ >/* nonlinear part and the objective function. */ >/* gradients of the bound constraints are never explicitly stored . */ >/* the relevant parts of the code have been changed to reflect this. */ >/* the parameters */ >/* n,nlin,nonlin,big */ >/* low */ >/* up */ >/* A = gres[][i], i=1,...,nlin */ >/* are to be set (or read from a file) in the basic routine user_init */ >/* the parameters tau0, del0, delmin etc keep there meaning and usage. */ >/* */ >/* A new form of evaluation interface is used for the gradients: */ >/* only the gradients of the nonlinear constraints must be evaluated, */ >/* for a given list of indizes, that is grad_c_i , i in list. */ >/* Also the nonlinear functions c_i(donlp2_x) are evaluated blockwise */ >/* bug fix: mes-file written before opened (correction of initial value) */ >/* bug fix: dead loop in o8qpdu in case of degenerate and illconditioned */ >/* primal solution */ >/* bug fix: log(b2n) now only for b2n>0 */ >/* bug fix in user_eval.c: nres replaced by nonlin */ >/* fix of a problem with glibc : bind and bind0 renamed to o8bind and o8bind0*/ >/*****************************************************************************/ > >/* **************************************************************************** */ >/* o8gene.h */ >/* **************************************************************************** */ > >/* **************************************************************************** */ >/* o8para.h */ >/* **************************************************************************** */ > >#include <stdio.h> >#include <stdlib.h> >#include <string.h> >#include <time.h> >#include <math.h> > >/* Macros */ > >#define max(A,B) ((A) > (B) ? (A) : (B)) >#define min(A,B) ((A) < (B) ? (A) : (B)) > >/* Types */ > >typedef int IINTEGER; >typedef int LLOGICAL; >enum {FFALSE,TTRUE}; >typedef float RREAL; >typedef double DDOUBLE; > > >#define X > >/* **************************************************************************** */ >/* o8comm.h */ >/* **************************************************************************** */ > >/* **************************************************************************** */ >/* o8fuco.h */ >/* **************************************************************************** */ > >X LLOGICAL *val,*llow,*lup; > >X IINTEGER n,nr,nres,nlin,nonlin, nstep, ndualm, mdualm; > >X DDOUBLE epsmac,tolmac,deldif; > >X char name[41]; > >X DDOUBLE epsdif; > >X LLOGICAL intakt,te0,te1,te2,te3,singul; >X LLOGICAL ident,eqres,silent,analyt,cold; > >X IINTEGER icf,icgf,cfincr,*cres,*cgres; > >X LLOGICAL ffuerr,*confuerr; > >/* special variables for the interface to old fashioned function */ >/* specification */ >/* can be removed is only problems in the new formulation are to be used */ >X IINTEGER nh,ng; > >X IINTEGER *nonlinlist; > >X IINTEGER **gunit; > >X LLOGICAL *gconst; > >X LLOGICAL *cfuerr; >/* this is necessary because of the different index positions used */ >/* in the old and new versions */ > > >X RREAL runtim; >X RREAL optite; >X DDOUBLE **accinf; >X IINTEGER itstep,phase; > >X DDOUBLE upsi,upsi0,upsi1,upsist,psi,psi0, > psi1,psist,psimin, > phi,phi0,phi1,phimin,fx,fx0,fx1, > fxst,donlp2_fmin,b2n,b2n0,xnorm,x0norm,sig0,dscal,ddnorm,d0norm; >X DDOUBLE sig,sigmin,dirder,cosphi,upsim; >X DDOUBLE *donlp2_x,*x0,*x1,*xmin,*d,*d0, > *dd,*difx,*resmin; > >X DDOUBLE *gradf,gfn,*qgf,*gphi0,*gphi1, > **gres,*gresn; > >X IINTEGER *perm,*perm1,*colno,rank; >X DDOUBLE **qr,*betaq,*diag, > *cscal, > *colle; > >/* colno also used o8qpso with double length ! */ > >X DDOUBLE **a,scalm,scalm2,*diag0,matsc; > >X IINTEGER *violis,*alist,*o8bind, > *o8bind0, *aalist,*clist; > >X DDOUBLE *u,*u0, > *w,*w1,*res, > *res0,*res1, > *resst,scf,scf0, > *yu,*slack,infeas,*work; > >X IINTEGER iterma; >X DDOUBLE del,del0,del01,delmin,tau0,tau,ny; >X DDOUBLE smalld,smallw,rho,rho1,eta,epsx,c1d, > scfmax,updmy0,tauqp,taufac,taumax; > >X DDOUBLE alpha,bbeta,theta,sigsm,sigla,delta,stptrm; >X DDOUBLE delta1,stmaxl; > >X DDOUBLE level; >X IINTEGER clow,lastdw,lastup,lastch; > >X FILE *prou,*meu; > >X DDOUBLE *ug,*og; > >X DDOUBLE *low,*up,big; > >X IINTEGER nreset; > >X DDOUBLE *xst; > > > >/* **************************************************************************** */ >/* o8fint.h */ >/* **************************************************************************** */ > >/* if bloc = TRUE then it is assumed that functionevaluation takes place */ >/* at once in an external routine and */ >/* that user_eval has been called before calling for evaluation of functions */ >/* the latter then simply consists in copying data */ >/* to donlp2's own data */ >/* user_eval must set valid = TRUE, if functionvalues are valid for the */ >/* current xtr */ >/* corr is set to true by donlp2, if the initial x does not satisfy */ >/* the bound constraints. x is modified in this case */ >/* difftype = 1,2,3 numerical differentiation by the ordinary forward */ >/* differences, by central differences or by Richardson-extrapolation */ >/* of order 6, requiring n, 2n , 6n additional function evaluations */ >/* respectively */ >/* epsfcn is the assumed precision of the function evaluation, to be */ >/* set by the user */ >/* taubnd: amount by which bound constraints may be violated during */ >/* finite differencing, set by the user */ > >X LLOGICAL bloc,valid,corr; >X IINTEGER difftype; >X DDOUBLE *xtr,*xsc,*fu,**fugrad, > **fud,epsfcn,taubnd; > >#undef X > >const DDOUBLE zero = 0.e0,one = 1.e0,two = 2.e0,three = 3.e0,four = 4.e0, >five = 5.e0,six = 6.e0,seven = 7.e0, >twom2 = .25e0,twop4 = 16.e0,twop11 = 2048.e0,onep3 = 1.3e0, >onep5 = 1.5e0,p2 = .2e0,p4 = .4e0,p5 = .5e0,p7 = .7e0,p8 = .8e0,p9 = .9e0, >c45 = 45.e0,tm12 = 1.e-12, >tm10 = 1.e-10,tm9 = 1.e-9,tm8 = 1.e-8,tm7 = 1.e-7,tm6 = 1.e-6,tm5 = 1.e-5, >tm4 = 1.e-4,tm3 = 1.e-3,tm2 = 1.e-2,tm1 = 1.e-1,tp1 = 1.e1,tp2 = 1.e2, >tp3 = 1.e3,tp4 = 1.e4 ; > >static IINTEGER qpterm,fcount,qpitma; > >static IINTEGER ndual,mi,me,iq; >static DDOUBLE *np,rnorm,rlow,**xj, > *ddual, **r,*ud, > *ud1; > >static IINTEGER iptr,iqtr,*aitr; >static DDOUBLE sstr,riitr; > > >#define ABC 0 > > >/* Global variables only used locally */ > >DDOUBLE *donlp2_yy; >DDOUBLE *o8bfgs_dg,*o8bfgs_adx,*o8bfgs_ltdx, > *o8bfgs_gtdx,*o8bfgs_updx,*o8bfgs_updz; >DDOUBLE *o8opti_qtx; >DDOUBLE *o8opti_yy,*o8opti_yx,*o8opti_trvec,*o8opti_con; >IINTEGER *o8opti_delist,*o8opti_bindba; >DDOUBLE *o8eval_con; >DDOUBLE *o8unim_step; >DDOUBLE *o8dec_qri,*o8dec_qri0; >DDOUBLE *o8elim_qri,*o8elim_y,*o8elim_rhsscal; >IINTEGER *o8elim_col; >DDOUBLE *o8sol_xl; >DDOUBLE *o8upd_sdiag,*o8upd_rn1,*o8upd_w; >DDOUBLE *o8qpdu_g0,*o8qpdu_ci0,*o8qpdu_cii, > *o8qpdu_cei, > *o8qpdu_xd,*o8qpdu_s,*o8qpdu_z,*o8qpdu_vr; >DDOUBLE *o8qpdu_xdold,*o8qpdu_udold; >IINTEGER *o8qpdu_ai,*o8qpdu_iai,*o8qpdu_iaexcl,*o8qpdu_aiold; >IINTEGER *o8qpdu_eqlist,*o8qpdu_iqlist; >DDOUBLE *o8qpdu_y,*o8qpdu_mult; >IINTEGER *o8qpdu_qpdel; >LLOGICAL *escongrad_errloc; >DDOUBLE *escongrad_fhelp1,*escongrad_fhelp2, > *escongrad_fhelp3,*escongrad_fhelp4, > *escongrad_fhelp5,*escongrad_fhelp6; > >typedef void (*func_void_void_t) (void); >typedef void (*func_void_type_liste_donlp_x_err_t)(IINTEGER type, IINTEGER liste[], > DDOUBLE donlp2_x[], DDOUBLE con[], LLOGICAL err[]); >typedef void (*func_void_liste_shift_donlp_x_grad_t)(IINTEGER liste[], IINTEGER shift , > DDOUBLE donlp2_x[], DDOUBLE **grad); >typedef void (*func_void_donlp_x_fx_t)(DDOUBLE donlp2_x[],DDOUBLE *fx); >typedef void (*func_void_donlp_x_gradf_t)(DDOUBLE donlp2_x[],DDOUBLE gradf[]); >typedef void (*func_void_mode_t)(IINTEGER mode); >typedef void (*func_void_t)(); > >//func_void_void_t user_init; >func_void_type_liste_donlp_x_err_t econ; >func_void_liste_shift_donlp_x_grad_t econgrad; >func_void_donlp_x_fx_t ef; >func_void_donlp_x_gradf_t egradf; >func_void_mode_t eval_extern; >func_void_t freemem; >func_void_t initialparams; >//func_void_void_t setup; >func_void_void_t solchk; >//func_void_void_t user_init; >//func_void_void_t user_init_size; >func_void_t allocatemem; > > >void donlp2(void) { > > void o8st (void); > void o8fin (void); > void o8opti (void); > DDOUBLE o8vecn (IINTEGER nl,IINTEGER nm,DDOUBLE donlp2_x[]); > void esgradf(DDOUBLE donlp2_x[],DDOUBLE gradf[]); > void global_mem_malloc(); > void global_mem_free(); > > void setup0 (void); > void setup (void); > void user_init_size (void); > void user_init (void); > void o8msg(IINTEGER num); > void o8elim(void); > static DDOUBLE term ; > static IINTEGER i,j,k,iz; > static char fil[13],xxx[9] = "xxxxxxxx",name1[9]; > >/* ---------------> step one of problem initialization by user */ > /* user_init_size must initialize n, nlin, nonlin, iterma, nstep. */ > /* These values may not be changed. */ > > user_init_size(); > ndualm = 2*n+nlin+nonlin; > mdualm = 2*(n+nlin+nonlin); > > /* allocate the global memory */ > > global_mem_malloc(); > > /* default settings of new parameters */ > > bloc = FFALSE; > analyt = TTRUE; > valid = FFALSE; > epsfcn = 1.e-16; > difftype = 3; > taubnd = 1.; > for (i = 1 ; i <= n ; i++) { > xsc[i] = one; > xtr[i] = zero; > } > epsdif = tm8; > for (i = 0 ; i <= nlin+nonlin ; i++) { > val[i] = FFALSE; > if ( i > 0 ) gresn[i] = one; > } > ffuerr = FFALSE; > fcount = 0; > > /* some standard initialization which may be overwritten by */ > /* user_init or setup */ > > silent = FFALSE; > > /* the interactive input feature is no longer supported here. for the sake */ > /* of easy revision the variable is supplied here however */ > /* if intakt is set TTRUE, output to protocol file is copied to stdout in */ > /* addition */ > > intakt = FFALSE; > te0 = TTRUE; /* control running optimizer on stdout */ > te1 = FFALSE; /* information about iteration in compressed form */ > te2 = FFALSE; /* detailed output in case of trouble */ > te3 = FFALSE; /* print gradients and Hessian */ > cold = TTRUE; > big = 1.e20 ; > >/* -> here user initialization continued */ > > > /* user_init must initialize analyt, epsdif, del0, tau0 , */ > /* low, up, and the gradients of the linear constraints, if any */ > /* (stored in gres[1:n][j], j=1,nlin) */ > /* bloc */ > /* analyt and if analyt = FFALSE then also epsfcn , taubnd , viobnd , */ > /* difftype */ > /* and the initial value for donlp2_x */ > /* may also change settings of all variables initialized above */ > > user_init(); > >/* open files for output if wanted */ > > j = 0; > while ( name[j] == ' ' ) { > j = j+1; > } > if ( name[j] == '\0' ) { > strcpy(name1,xxx); > } else { > k = j+1; > while ( name[k] != ' ' && name[k] != '\0' && k-j < 8 ) { > k = k+1; > } > strncpy(name1,&name[j],k-j); > name1[k-j] = '\0'; > > for (i = 0 ; i <= k-j-1 ; i++) { > iz = name1[i]; > if ( iz < 48 || ( iz > 57 && iz < 65 ) > || ( iz > 90 && iz < 97 ) || iz > 122 ) name1[i] = 'x'; > } > if ( k - j < 8 ) strncat(name1,xxx,8-k+j); > } > if ( ! silent ) > { > strcpy(fil,name1); > meu = fopen(strcat(fil,".mes"),"w"); > strcpy(fil,name1); > prou = fopen(strcat(fil,".pro"),"w"); > } > > > > >/* ************************************************************************** */ >/* ***** automatic correction of del0 */ >/******************************************************************************/ > > for ( i = 1 ; i <= n ; i++ ) > { > if ( xsc[i] == zero ) > { > fprintf(stderr,"scaling variable %i is zero\n",i); > exit(1); > } > } > nres = n+nlin+nonlin ; > for ( i = 1 ; i <= nres ; i++ ) > { > if ( i <= n ) > { > term = xsc[i] ; > } > else > { > term = one ; > } > if (! (low[i] == up[i]) ) > { > del0 = min ( del0 , (up[i]-low[i])/(four*term) ); > } > } > delmin = tm6; > > for (i = 1 ; i <= n ; i++) > { > > ug[i] = low[i]; > if ( ug[i] <= -big ) > { > llow[i] = FFALSE ; > } > else > { > llow[i] = TTRUE ; > } > og[i] = up[i]; > if ( og[i] >= big ) > { > lup[i] = FFALSE ; > } > else > { > lup[i] = TTRUE ; > } > > } > > for (i=1; i<=n; i++) > { > if ( donlp2_x[i] < low[i] || donlp2_x[i] > up[i] ) > { > corr = TTRUE ; > if ( llow[i] && lup[i] ) donlp2_x[i] = (up[i]+low[i])/2.0; > if ( llow[i] && !lup[i] ) donlp2_x[i] = low[i]+2.0*del0; > if ( lup[i] && !llow[i]) donlp2_x[i] = up[i]-2.0*del0; > } > } > > if ( corr && ! silent ) o8msg(13); > > for (i = 1 ; i <= n ; i++) { > xst[i] = donlp2_x[i]; > donlp2_x[i] = donlp2_x[i]/xsc[i]; > o8opti_yy[i] = xsc[i]; > > } > for (i = 1 ; i <= n ; i++) { > > /* ug and og have been evaluted for the original variables */ > /* here we use them internally for the scaled ones */ > /* ug = low up to scaling etc */ > ug[i] = ug[i]/xsc[i]; > og[i] = og[i]/xsc[i]; > } > for ( i = 1 ; i <= nlin ; i++ ) > { > /* rescale the gradients of the linear constraints */ > gresn[i] = zero ; > for ( j = 1 ; j <= n ; j++ ) > { > gres[j][i] = xsc[j]*gres[j][i] ; > gresn[i] += pow(gres[j][i],2); > } > gresn[i] = max ( one , zero ) ; > cgres[i] = 1 ; > val[i] = TTRUE ; > } >/* the sign of the gradient is determined from the current constraint */ >/* and stored in gres[0][i]. done afterwards */ > nreset = n; > > o8st(); > > /* setup may change standard settings of parameters */ > /* and add some computations in the user environment */ > qpitma = nres ; > setup(); > > if ( taubnd <= 0 ) { > fprintf(stderr,"taubnd le zero is not allowed "); > exit(1); > } > for ( i = 1 ; i<= n ; i++ ) { > if ( o8opti_yy[i] != xsc[i] ) { > fprintf(stderr,"setup has changed xsc, not allowed"); > exit(1); > } > } > /* preevaluation of gradients of linear functions */ > /* done now in user_init */ > > > runtim = clock(); > /* check for redundant linear equality constraints */ > o8elim(); > > > > /* call the optimizer */ > > o8opti(); > > runtim = (clock()-runtim)/CLOCKS_PER_SEC; > > /* do final solution check and output */ > > o8fin(); > > /* free up memory */ > > global_mem_free(); > > return; >} > >/* **************************************************************************** */ >/* initialization program , standard parameter settings done here */ >/* **************************************************************************** */ >void o8st(void) { > return; >} > >/* **************************************************************************** */ >/* do final solution check and output */ >/* **************************************************************************** */ >void o8fin(void) { > return; >} >/* **************************************************************************** */ >/* prints informations if te2 = TTRUE */ >/* **************************************************************************** */ >void o8info(IINTEGER icase) { >} > >/* **************************************************************************** */ >/* computation of new scaling factors for L1-penalty-function */ >/* **************************************************************************** */ >void o8sce(void) { > return; >} > >/* **************************************************************************** */ >/* computation of the Pantoja-Mayne BFGS-update of hessian */ >/* **************************************************************************** */ >void o8bfgs(void) { > return; >} > >/* **************************************************************************** */ >/* write short information on standard out */ >/* **************************************************************************** */ >void o8shms(void) { > return; >} > >/* **************************************************************************** */ >/* write messages on "special events" */ >/* **************************************************************************** */ >void o8msg(IINTEGER num) >{ >} > >/* **************************************************************************** */ >/* equality constrained recursive quadratic programming with */ >/* multiple inactivation and superlinearly convergent projected */ >/* BFGS-update (version 12/93 Spellucci, modified subsequently ) */ >/* this version from 31.10.2001 */ >/* **************************************************************************** */ >void o8opti(void) >{ > void o8info (IINTEGER icase); > void o8sce (void); > void o8bfgs (void); > void o8shms (void); > void o8msg (IINTEGER num); > void o8inim (void); > void o8dird (void); > void o8cutd (void); > void o8smax (void); > void o8unim (DDOUBLE sig1th); > void o8egph (DDOUBLE gphi[]); > void o8dec (IINTEGER nlow,IINTEGER nrl); > void o8ht (IINTEGER id,IINTEGER incr,IINTEGER is1, > IINTEGER is2,IINTEGER m, > DDOUBLE **a,DDOUBLE bbeta[],DDOUBLE b[],DDOUBLE c[]); > void o8sol (IINTEGER nlow,IINTEGER nup,DDOUBLE b[],DDOUBLE donlp2_x[]); > void o8solt (IINTEGER nlow,IINTEGER nup,DDOUBLE b[],DDOUBLE donlp2_x[]); > void o8rght (DDOUBLE **a,DDOUBLE b[],DDOUBLE y[], > DDOUBLE *yl,IINTEGER n); > void o8left (DDOUBLE **a,DDOUBLE b[],DDOUBLE y[], > DDOUBLE *yl,IINTEGER n); > DDOUBLE o8vecn (IINTEGER nl,IINTEGER nm,DDOUBLE donlp2_x[]); > DDOUBLE o8sc3 ( IINTEGER i , IINTEGER j , IINTEGER k , > DDOUBLE **mat, DDOUBLE donlp2_x[] ); > DDOUBLE o8sc4 ( IINTEGER i , IINTEGER j , IINTEGER k , > DDOUBLE **mat) ; > void o8qpdu (void); > void esf (DDOUBLE donlp2_x[],DDOUBLE *fx); > void esgradf(DDOUBLE donlp2_x[],DDOUBLE gradf[]); > void escon (IINTEGER type, IINTEGER liste[], DDOUBLE donlp2_x[], > DDOUBLE constr[], > LLOGICAL errlist[]); > void escongrad( IINTEGER liste[], IINTEGER shift , > DDOUBLE donlp2_x[], DDOUBLE **grad_constr); > > void user_eval(DDOUBLE xvar[],IINTEGER mode); > void newx ( DDOUBLE donlp2_x[], DDOUBLE u[], IINTEGER itstep, DDOUBLE **accinf, > LLOGICAL *cont ) ; > > static IINTEGER l0,i,j,k,csssig,csirup,csreg,cschgx; > static IINTEGER csmdph; > static DDOUBLE delsig,delx,sum,term; > static DDOUBLE umin,term1,scfh,unorm; > static DDOUBLE compl,del1; > static IINTEGER iumin,rank0,nr0,csdifx,clwold; > static IINTEGER nrbas; > static DDOUBLE eps,delold,uminsc,fac,slackn,tauqp0,denom; > static LLOGICAL nperm,qpnew,etaini; > /* added feature 25.01.2000 */ > static LLOGICAL viobnd ; > /* end added feature */ > /* added feature user interrupt */ > static LLOGICAL cont ; > /* initialization */ > > > for (i = 1 ; i <= n ; i++) > { > d[i] = zero; > d0[i] = zero; > } > itstep = 0; > alist[0] = 0; /* active general constraints */ > aalist[0] = 0 ; /* all active constraints in 1,...,2*nres */ > o8opti_delist[0] = 0; > violis[0] = 0; > clist[0] = 0; /* the nonlinear active constraints */ > upsi = zero; > psi = zero; > psi0 = zero; > sig0 = zero; > d0norm = one; > unorm = one; > fx = zero ; > fx0 = zero ; > /* in order to have cosphi well defined for itstep = 1 */ > > ddnorm = one; > del = del0; > > /* count successive regularization steps */ > > csreg = 0; > > /* count successive small changes in donlp2_x */ > > cschgx = 0; > > /* count small differences of fx */ > > csdifx = 0; > > /* count irregular quasi-Newton-updates */ > > csirup = 0; > > /* count successive small stepsizes */ > > csssig = 0; > > /* count successive small differences of penalty-function */ > > csmdph = 0; > /* moved to o8st : matsc = one; */ > > /* formerly tauqp = tp2 is worse */ > > tauqp = one; > nperm = FFALSE; > ident = FFALSE; > etaini = FFALSE; > if ( n > 100 || nres > 100 ) te3 = FFALSE; > for (i = 1 ; i <= n ; i ++) > { > perm[i] = i; > perm1[i] = i; > } > /* compute first list of binding constraints */ > for ( i = 1 ; i <= nres ; i++ ) > { > /* equality constraint is always active . use the odd numbered formal > constraint only */ > if ( low[i] == up[i] ) > { > aalist[0] += 1 ; > aalist[aalist[0]] = 2*i-1 ; > if ( i > n ) > /* a general constraint */ > { > alist[0] += 1; > alist[alist[0]] = i-n ; > gres[0][i-n] = one ; > } > if ( i > n+nlin ) > { > /* a nonlinear equality constraint */ > clist[0] += 1; > clist[clist[0]] = i-n-nlin ; > } > > o8bind[2*i-1] = 1 ; > o8bind[2*i] = 0 ; > o8bind0[2*i-1] = 1 ; > o8bind0[2*i] = 0 ; > } > } > if ( analyt ) > { > eps = min(epsx,epsx); > } > else > { > eps = epsdif; > if ( epsx < pow(epsdif,2) ) epsx = pow(epsdif,2); > } > eps = max(epsmac*tp3,min(tm3,eps)); > > /* calling for external function evaluation necessary only if corr = TTRUE */ > > /* function and gradient values, from xtr = donlp2_x*xsc */ > > /* remember that that es_routines consider bloc/non bloc call */ > /* in case of bloc == TTRUE the evaluation has been done in o8st already */ > > for ( i = 1 ; i <= n ; i++ ) > { res[2*i-1] = donlp2_x[i]-ug[i]; /* the scaled variables */ > res[2*i] = og[i]-donlp2_x[i] ; > } > for ( i = 1 ; i<=nlin ; i++ ) > { > term = o8sc3(1,n,i,gres,donlp2_x); /* the gradients of the linear */ > /* constraints are scaled accordingly */ > /* in o8st already */ > res[2*(n+i)-1] = term - low[i+n] ; > res[2*(n+i)] = up[i+n] - term ; > cres[i] += 1 ; > } > /* ------- > evaluate _all_ nonlinear constraints */ > for ( i = 1 ; i <= nonlin ; i++ ) > { > confuerr[i] = FFALSE ; > } > escon(1,clist,donlp2_x,o8opti_con,confuerr); > /* 1. parameter = 1 makes second parameter meaningless */ > /* */ > > /* this evaluates the nonlinear constraints at donlp2_x */ > for ( i = 1 ; i <= nonlin ; i++ ) > { > res[2*(i+n+nlin)-1] = o8opti_con[i]-low[i+n+nlin] ; > res[2*(i+n+nlin)] = up[i+n+nlin]-o8opti_con[i] ; > } > for ( i =1 ; i <= nres ; i++ ) >/* check for active inequality constraints and compute the scaled and */ >/* the unscaled penalty function */ > { > if ( low[i] == up[i] ) > { > term = fabs(res[2*i-1]); > term1 = zero ; > } > else > { > term = - min(zero,res[2*i-1]); > term1= - min(zero,res[2*i]); > if ( res[2*i-1] <= delmin ) > { > aalist[0] += 1 ; > aalist[aalist[0]] = 2*i-1 ; > o8bind[2*i-1] = 1 ; > } > if ( res[2*i] <= delmin ) > { > aalist[0] += 1 ; > aalist[aalist[0]] = 2*i ; > o8bind[2*i] = 1 ; > } > if ( o8bind[2*i-1]+o8bind[2*i] == 2 ) > { > fprintf(stderr," donlp2: lower and upper bound are binding "); > fprintf(stdout," decrease delmin !\n "); > exit(1); > } > if ( i > n && (o8bind[2*i-1]+o8bind[2*i] == 1 ) ) > { > alist[0] += 1 ; > alist[alist[0]] = i-n ; > if ( o8bind[2*i-1] == 1 ) > { > gres[0][i-n] = one; > } > else > { > gres[0][i-n] = - one ; > } > } > if ( i > n+nlin && (o8bind[2*i-1]+o8bind[2*i] == 1 ) ) > { > clist[0] += 1 ; > clist[clist[0]] = i-n-nlin ; > } > > > } >/* because of the definition of delmin the lower and the upper bound cannot */ >/* be active simultaneously */ > upsi += term+term1 ; > psi += w[2*i-1]*term+w[2*i]*term1 ; > } > /* only the gradients of the nonlinear constraints need to be evaluated */ > /* they are stored in gres[][nlin+clist[]] */ > escongrad(clist,nlin,donlp2_x,gres); > for ( i = 1 ; i <= clist[0] ; i++ ) > { > cgres[clist[i]+nlin] += 1 ; > val[clist[i]+nlin] = TTRUE ; > gresn[clist[i]+nlin] = max ( one , one ); > } > > /* end first evaluation of all contraints */ > L100: > > /* ******************************************** */ > /* obtaining a point feasible within tau0 first */ > /* ******************************************** */ > > if ( upsi >= tau0 ) > { > scf = zero; > phase = -1; > } > else > { > ffuerr = FFALSE; > > esf(donlp2_x,&fx); > icf += 1 ; > if ( ffuerr ) > { > if ( ! silent ) o8msg(22); > optite = -9; > > return; > } > if ( ! val[0] ) > { > > /* we assume that the gradient evaluation can be done whenever */ > /* the function has been evaluated */ > > esgradf(donlp2_x,gradf); > icgf += 1 ; > val[0] = TTRUE; > } > scf = one; > phase = 0; > fxst = fx; > psist = psi; > upsist = upsi; > for (j = 1 ; j <= 2*nres ; j++) > { > resst[j] = res[j]; > } > eta = zero; > } > L200: > > > > for (i = 1 ; i <= 2*nres ; i++) > { > u[i] = zero; > } > iumin = 0; > uminsc = zero; > for (i = 1 ; i <= rank ; i++) > { > unorm = max(unorm,fabs(yu[i])); > k = aalist[colno[i]]; > u[k] = -yu[i]; > term = one ; > if ( k > 2*n ) term=gresn[(k-2*n+1)/2]; > if ( !(low[(k+1)/2] == up[(k+1)/2]) ) > { > if ( -yu[i]/term < uminsc ) > { > iumin = k; > uminsc = -yu[i]/term; > } > } > } > if ( scf != zero ) > { > for (i = 1 ; i <= n ; i++) > { > o8opti_yx[i] = scf*gradf[i]; > for (j = 1 ; j <= nlin+nonlin ; j++) > { > o8opti_yx[i] = o8opti_yx[i]-gres[i][j]*(u[2*j-1+2*n]-u[2*j+2*n]); > } > } > for (i = 1 ; i <= n ; i++ ) > { > o8opti_yx[i] -= xsc[i]*(u[2*i-1]-u[2*i]) ; > } > /* the l2-norm of the original unscaled gradient of the lagrangian */ > } > > if ( clist[0] != 0 ) > { > escongrad(clist,nlin,donlp2_x,gres); > for ( i = 1 ; i <= clist[0] ; i++ ) > { > cgres[clist[i]+nlin] += 1 ; > val[clist[i]+nlin] = TTRUE ; > gresn[clist[i]+nlin] = max(one,one); > } > } > > umin = zero; > unorm = zero; > for (i = 1 ; i <= 2*nres ; i++) > { > u[i] = zero; > } > iumin = 0; > uminsc = zero; > for (i = 1 ; i <= rank ; i++) > { > unorm = max(unorm,fabs(yu[i])); > k = aalist[colno[i]]; > u[k] = -yu[i]; > term = one ; > term = ( k > 2*n )? gresn[(k-2*n+1)/2] : xsc[(k+1)/2]; > if ( !(low[(k+1)/2] == up[(k+1)/2]) ) > { > if ( -yu[i]/term < uminsc ) > { > iumin = k; > uminsc = -yu[i]/term; > } > } > > } > if ( scf != zero ) > { > > for (i = 1 ; i <= n ; i++) > { > o8opti_yx[i] = scf*gradf[i]; > for (j = 1 ; j <= nlin+nonlin ; j++) > { > o8opti_yx[i] -= gres[i][j]*(u[2*j-1+2*n]-u[2*j+2*n]); > } > } > for (i = 1 ; i <= n ; i++ ) > { > o8opti_yx[i] -= xsc[i]*(u[2*i-1]-u[2*i]) ; > } > b2n = o8vecn(1,n,o8opti_yx)/scf; > } > > eqres = TTRUE; > for (i = 1 ; i <= 2*nres ; i ++) > { > eqres = eqres && ( o8bind[i] == o8bind0[i] ); > } > /* compute condition number estimators of diag-r and diag of */ > /* Cholesky-decomposition of b */ > > if ( nr > 1 ) > { > term = zero; > term1 = one; > for (i = 1 ; i <= nr ; i++) > { > term = max(term, fabs(diag[i])); > term1 = min(term1,fabs(diag[i])); > } > accinf[itstep][13] = term/term1; > } > else if ( nr == 1 ) > { > accinf[itstep][13] = 1.; > } > else > { > accinf[itstep][13] = -1.; > } > term = fabs(a[1][1]); > term1 = fabs(a[1][1]); > i = 2; > while ( i <= n ) > { > term = max(term, fabs(a[i][i])); > term1 = min(term1,fabs(a[i][i])); > i = i+1; > } > accinf[itstep][14] = pow(term/term1,2); > > if ( ! silent ) o8info(5); > > /* since a represents the Cholesky-factor, this square */ > slackn = zero; > for (i = 1 ; i <= nres ; i++) > { > term = res[2*i-1]; > term1 = res[2*i]; > denom = (i > n ) ? gresn[i-n] : max(one,xsc[i]); > slackn += > max(zero,term/denom-delmin)*max(zero,u[2*i-1]-smallw)/denom > +max(zero,term1/denom-delmin)*max(zero,u[2*i]-smallw)/denom; > } > > if ( umin >= -smallw && > slackn <= delmin*smallw*nres && > upsi <= nres*delmin && upsi0 <= nres*delmin > && fabs(fx-fx0) <= eps*(fabs(fx)+one) && > b2n != -one && > b2n <= tp2*epsx*(gfn+one) ) > { > > csdifx = csdifx+1; > } > else > { > csdifx = 0; > } > /* compute damping factor for tangential component if upsi>tau0/2 */ > > scf0 = one; > if ( phase >= 0 && upsi > tau0*p5 ) > scf0 = max(one/scfmax,(two*(tau0-upsi)/tau0)*upsi*tm1/max(one,gfn) )/scf; > accinf[itstep][15] = scf0; > > /* compute new penalty weights : regular case */ > > clwold = clow; > if ( phase >= 0 ) o8sce(); > if ( clow > clwold ) > { > > /* tau_qp depends on the (new) weights */ > > term = w[1]; > for (i = 1 ; i <= 2*nres ; i++) > { > term = max(term,w[i]); > } > tauqp = max(one,min(tauqp,term)); > } > > /* terminate if correction is small */ > > if ( ddnorm <= epsx*(xnorm+epsx) && upsi <= delmin > && b2n != -one && uminsc >= -smallw && b2n <= epsx*(gfn+one) ) > { > > optite = one; > > return; > } > L350: > > /* reenter from the singular case: dirder has been computed already */ > > accinf[itstep][16] = xnorm; > accinf[itstep][17] = ddnorm; > accinf[itstep][18] = phase; > > /* compute stepsize */ > > cfincr = icf; > > /* if no descent direction is obtained, check whether restarting the method */ > /* might help */ > > > if ( phase == 2 && ddnorm > xnorm*epsmac && ! singul > && ! viobnd && nres > 0 ) > { > /* compute the Maratos correction using a difference approximation */ > if ( bloc ) > { > valid = FFALSE ; > /* user_eval must reset valid */ > user_eval(x1,-1); > } > > if ( term > p5*ddnorm ) > { > > /* second order correction almost as large as first order one: */ > /* not useful */ > > for (i = 1 ; i <= n ; i++) > { > dd[i] = zero; > } > if ( ! silent ) o8msg(6); > } > } > L355: > > if ( ! silent ) o8info(7); > > sig = min(one,stmaxl); >// printf("3\n"); > o8unim(sig); >// printf("4\n"); > > L360: > > if ( stptrm < zero ) > { > > /* stepsize selection failed */ > > if ( ! ident ) > { > > /* try restart with a = identity scaled */ > > if ( ! silent ) o8msg(7); > ident = TTRUE; > o8opti_delist[0] = 0; > violis[0] = 0; > csreg = 0; > csssig = 0; > csirup = 0; > > o8inim(); > > aalist[0] = nrbas; > for (i = 1 ; i <= 2*nres ; i++) > { > o8bind[i] = o8opti_bindba[i]; > } > if ( upsi >= tau0 ) > { > > goto L100; > > } > else > { > > goto L200; > } > } > if ( ! singul && ident && accinf[itstep][13] > tp4 && nres > 0 ) > { > > /* try the full SQP-direction */ > /* this may be the third try for this point */ > > singul = TTRUE; > ident = TTRUE; > > goto L400; > } > if ( stptrm == -two ) > { > optite = -four; > > return; > } > if ( sig == zero && optite == -one ) return; > > /* unidimensional search unsuccessfully terminated */ > > optite = -two; > > return; > } > > L400: > > if ( phase == -1 ) > { > > goto L100; > > } > else > { > > goto L200; > } >} > >void o8inim(void) >{ > return; >} >void o8dird(void) >{ > return; >} > >/* **************************************************************************** */ >/* cut d if appropriate and rescale */ >/* **************************************************************************** */ >void o8cutd(void) >{ > return; >} > >/* **************************************************************************** */ >/* compute maximum stepsize stmaxl such that projection on the box of lower */ >/* and upper bounds changes for sig in [0,stmaxl],if such exists */ >/* **************************************************************************** */ >void o8smax(void) { > return; >} > >/* **************************************************************************** */ >/* restore the best point found so far to be the current new point */ >/* **************************************************************************** */ >void o8rest(void) { > return; >} > >/* **************************************************************************** */ >/* save the best point found so far in ...min variables */ >/* **************************************************************************** */ >void o8save(void) { > return; >} > >/* **************************************************************************** */ >/* evaluate the functions at the new point */ >/* **************************************************************************** */ >void o8eval(DDOUBLE sigact,DDOUBLE *sigres,LLOGICAL *reject,LLOGICAL *error) >{ > return; >} > >/* **************************************************************************** */ >/* determination of stepsize by an Armijo-like test for descent */ >/* **************************************************************************** */ >void o8unim(DDOUBLE sig1th) { > return; >} > >/* **************************************************************************** */ >/* scalar product of two vectors or parts of vectors */ >/* **************************************************************************** */ >DDOUBLE o8sc1(IINTEGER i,IINTEGER j,DDOUBLE a[],DDOUBLE b[]) >{ > return zero; >} > >/* **************************************************************************** */ >/* multiply row j of matrix a with vector b */ >/* **************************************************************************** */ >DDOUBLE o8sc2(IINTEGER n,IINTEGER m,IINTEGER j,DDOUBLE **a,DDOUBLE b[]) { > return zero; >} > >/* **************************************************************************** */ >/* multiply column j section (n to m) of matrix a with vector b */ >/* **************************************************************************** */ >DDOUBLE o8sc3(IINTEGER n,IINTEGER m,IINTEGER j,DDOUBLE **a,DDOUBLE b[]) >{ > return zero; >} >/*************************************************************************/ >DDOUBLE o8sc3b(IINTEGER n,IINTEGER m,IINTEGER j,DDOUBLE **a,DDOUBLE b[]) >{ > return zero; >} > >DDOUBLE o8sc4(IINTEGER n,IINTEGER m,IINTEGER j,DDOUBLE **a) >{ > return zero; >} > >DDOUBLE o8sc3_(IINTEGER n,IINTEGER m,IINTEGER j,DDOUBLE **a,DDOUBLE b[]) { > return zero; >} > >/* **************************************************************************** */ >/* subprogram for structured output of a submatrix a[ma][na] */ >/* on channel lognum in fix or float format with heading "head". */ >/* uses a fixed format string with 70 print columns */ >/* **************************************************************************** */ >void o8mdru(DDOUBLE **a,IINTEGER ma,IINTEGER na,char head[], > FILE *lognum,LLOGICAL fix) { > return; >} >void o8mdru_(DDOUBLE **a,IINTEGER ma,IINTEGER na,char head[], > FILE *lognum,LLOGICAL fix) { > return; >} > >/* **************************************************************************** */ >/* compute gradient of lagrangian */ >/* **************************************************************************** */ >void o8egph(DDOUBLE gphi[]) { > return; >} > >/* detection and elimination of redundant equality constraints */ >/* this is done making them always inactive inequality constraints */ >/* by redefining low and up */ >void o8elim(void) { > return; >} > > >/* **************************************************************************** */ >/* QR-decomposition of matrix of gradients of binding constraints */ >/* this set may be expanded using multiple calls to o8dec. */ >/* No exit on singular r-factor here. Information on */ >/* the decompostion is stored in betaq and in and below the */ >/* diagonal of qr. r-factor is stored in diag (diagonal ) and */ >/* above the diagonal of qr. cscal is the column scaling of the */ >/* original matrix. Column pivoting is done here and is stored in colno */ >/* **************************************************************************** */ >void o8dec(IINTEGER nlow,IINTEGER nrl) { > return; >} > >/* **************************************************************************** */ >/* application of Householder transformations */ >/* **************************************************************************** */ >void o8ht(IINTEGER id,IINTEGER incr,IINTEGER is1,IINTEGER is2,IINTEGER m, > DDOUBLE **a,DDOUBLE bbeta[],DDOUBLE b[],DDOUBLE c[]) { > return; >} > >/* **************************************************************************** */ >/* solve triangular system r*donlp2_x = b, r defined by Householder-QR- */ >/* decomposition decomp (with column scaling) */ >/* **************************************************************************** */ >void o8sol(IINTEGER nlow,IINTEGER nup,DDOUBLE b[],DDOUBLE donlp2_x[]) { > return; >} > >/* **************************************************************************** */ >/* solve triangular system r(transpose)*donlp2_x = b, r defined by */ >/* Householder-QR-decomposition decomp (with column scaling] */ >/* **************************************************************************** */ >void o8solt(IINTEGER nlow,IINTEGER nup,DDOUBLE b[],DDOUBLE donlp2_x[]) { > return; >} > >/* **************************************************************************** */ >/* lenght of vector (a,b). numerically stable version with */ >/* overflow / underflow saveguard */ >/* **************************************************************************** */ >DDOUBLE o8dsq1(DDOUBLE a,DDOUBLE b) { > return zero; >} > >/* **************************************************************************** */ >/* computes the upper triangular Cholesky-factor */ >/* **************************************************************************** */ >void o8upd(DDOUBLE **r,DDOUBLE z[],DDOUBLE y[],IINTEGER n,LLOGICAL *fail) { > return; >} > >/* **************************************************************************** */ >/* o8rght */ >/* **************************************************************************** */ >void o8rght(DDOUBLE **a,DDOUBLE b[],DDOUBLE y[],DDOUBLE *yl,IINTEGER n) { > return; >} > >/* **************************************************************************** */ >/* o8left */ >/* **************************************************************************** */ >void o8left(DDOUBLE **a,DDOUBLE b[],DDOUBLE y[],DDOUBLE *yl,IINTEGER n) { > return; >} > >/* **************************************************************************** */ >/* euclidean norm of donlp2_x , avoid overflow */ >/* **************************************************************************** */ >DDOUBLE o8vecn(IINTEGER nl,IINTEGER nm,DDOUBLE donlp2_x[]) { > return (zero); >} > >/* ************************************************************************* */ >/* solution of extended quadratic program */ >/* */ >/* scf*gradf(donlp2_x)*dir+(1/2)*dir*a*dir+summe(tauqp*sl[i]+( my/2)*pow(sl[i],2) )*/ >/* minimal subject to */ >/* sl[i] >= 0, i = 1,...,nr (the slacks) */ >/* (gres[.][j]*dir+res[j])+vz*sl[j] = 0, j = 1,...,nh vz = -sign(res[j]) */ >/* (gres[.][aalist[j]])*dir+res[aalist[j]])+sl[j] >= 0, j = nh+1,....,nr */ >/* the weight tauqp is adapted during solution */ >/* the quasi-Newton-matrix a is taken from o8comm.h */ >/* a is regularized if not sufficiently well conditioned */ >/* the resulting dir=xd[1+nr],...,xd[n+nr] is a direction of descent for */ >/* the Zangwill function of the corresponding nonlinear */ >/* optimization problem */ >/* f(donlp2_x) = min, res[j] = 0, j = 1,..nh, res[j] >= 0 , j = nh+1,nres */ >/* at the currrent point donlp2_x if the weight tauqp is chosen appropriately */ >/* the quadratic programming problem is solved using the method */ >/* of Goldfarb and Idnani */ >/* variables are stored in xd (solution) and ud (multipliers) */ >/* in the following order xd = ( sl[j], j = 1,nr; dir = direct. of desc.) */ >/* ud = ( multipliers for sl[i] >= 0 , i = 1,..,nr ; */ >/* multipliers for the equality constraints , */ >/* multipliers for the general inequality constraints ) */ >/* ***************************************************************************/ >void o8qpdu(void) >{ > > void o8info(IINTEGER icase); > void o8msg (IINTEGER num); > void o8dird(void); > void o8cutd(void); > DDOUBLE o8sc1 (IINTEGER i,IINTEGER j,DDOUBLE a[],DDOUBLE b[]); > DDOUBLE o8vecn(IINTEGER nl,IINTEGER nm,DDOUBLE donlp2_x[]); > void o8zup (DDOUBLE z[]); > void o8rup (DDOUBLE rv[]); > void o8dlcd(IINTEGER ai[],IINTEGER l); > void o8adcd(void); > void o8rinv(IINTEGER n,DDOUBLE **a,IINTEGER ndual,DDOUBLE **donlp2_x); > > static DDOUBLE infe1,s1,s2,tiny, > my,zz,ss,su,t,t1,t2,f,fmax,psid,c1,c2,cdiag,term, > su1,su2,condr,infiny,term1, > diff0; > static IINTEGER i,j,k,ip,l,incr,nosucc; > static LLOGICAL wlow; > > ndual = n+nr; /* nr=aalist[0] */ > > /* number of equality constraints in QP-problem */ > incr = nr ; > mi = 0 ; > me = 0 ; > /* compute the index lists of the QP's equality constraints > eqlist and inequality constraints iqlist from aalist and > low and up */ > for ( i = 1; i <= aalist[0] ; i++ ) > { > l = (aalist[i]+1)/2; > if ( low[l] == up[l] ) > { > me += 1 ; > o8qpdu_eqlist[me] = aalist[i] ; > } > else > { > mi += 1 ; > o8qpdu_iqlist[mi+incr] = aalist[i]+nr ; > } > } > for ( i = 1 ; i <= nr ; i++ ) > { > o8qpdu_iqlist[i] = i ; /* the slacks of the QP model */ > } > mi += nr ; > > /* QP inequality constraints = active constraints - */ > /* equality-constraints + slack's */ > > infiny = epsmac/tolmac; > if ( analyt ) > { > tiny = 2*nr*epsmac*tp3; > } > else > { > tiny = 2*nr*max(epsdif,epsmac*tp3); > } > > qpterm = 0; > for (i = 1 ; i <= nr ; i++) > { > > /* check gradients of active constraints against zero */ > o8qpdu_mult [i] = one ; > /* for a zero gradient the slack cannot be diminished:*/ > /* set mult[i] = 0 */ > l = aalist[i] ; > if ( l > 2*n ) > { /* not a bound */ > for (j = 1 ; j <= n ; j++) > { > o8qpdu_y[j] = gres[j][(l-2*n+1)/2]; > } > if ( o8vecn(1,n,o8qpdu_y) == zero ) > { > o8qpdu_mult[i] = zero; > if ( ! silent ) o8msg(8); > } > } > } > /* restart point in case of increase of tauqp */ > > L10: > nosucc = 0 ; > /* initialize matrices j and r */ > > for (i = 1 ; i <= ndual ; i++) > { > ddual[i] = zero; > for (j = 1 ; j <= ndual ; j++) > { > r[j][i] = zero; > xj[j][i] = zero; > } > } > rnorm = one; > rlow = one; > term1 = zero; > for (i = 1 ; i <= 2*nres ; i++) > { > u[i] = zero; > if ( w[i] > term1 ) term1 = w[i]; > } > accinf[itstep][19] = clow; > accinf[itstep][20] = term1; > accinf[itstep][31] = tauqp; > for (i = 1 ; i <= 2*nr ; i++) > { /* the multipliers in the QP */ > ud[i] = zero; > } > c1 = fabs(a[1][1]); > for (i = 1 ; i <= n ; i++) > { > c1 = max(c1,fabs(a[i][i])); > } > c1 = c1*tp1; > > /* we require much more regularity of a in the singular case */ > > for (i = 1 ; i <= n ; i++) > { > if ( fabs(a[i][i]) < rho1*c1 ) a[i][i] = rho1*c1; > } > /* invert the Cholesky-factor and store in > the right lower block of dimension n of xj ( Idnanis j-matrix) */ > > o8rinv(n,a,ndual,xj); > > c1 = fabs(a[1][1]); > incr = nr; > c2 = fabs(xj[1+incr][1+incr]); > for (i = 1 ; i <= n ; i++) > { > c1 = max(c1,fabs(a[i][i])); > c2 = max(c2,fabs(xj[i+incr][i+incr])); > } > my = zero; > for (i = 1 ; i <= n ; i++) > { > for (j = i ; j <= n ; j++) > { > my = my+pow(a[i][j],2); > } > } > my = my/n; > cdiag = one/my; > for ( i = 1 ; i <= incr ; i++) > { > xj[i][i] = cdiag; > /* the quadratic part for the slacks */ > /* has same order of magnitude as the part */ > /* coming from the original Lagrangian */ > } > for (i = 1 ; i <= ndual ; i++) > { > if ( i > incr ) o8qpdu_g0[i] = gradf[i-incr]*scf; > if ( i <= incr ) o8qpdu_g0[i] = tauqp; > /* the linear part for the slacks */ > } > /* compute unconstrained solution */ > /* the Cholesky-factor of a is stored in the upper triangle */ > > for (i = 1 ; i <= n ; i++) > { > su = zero; > for (j = 1 ; j <= i-1 ; j++) > { > su = su+a[j][i]*o8qpdu_y[j+incr]; > } > o8qpdu_y[i+incr] = (o8qpdu_g0[i+incr]-su)/a[i][i]; > } > for (i = n ; i >= 1 ; i--) > { > su = zero; > for (j = i+1 ; j <= n ; j++) > { > su = su+a[i][j]*o8qpdu_y[j+incr]; > } > o8qpdu_y[i+incr] = (o8qpdu_y[i+incr]-su)/a[i][i]; > } > for (i = 1 ; i <= incr ; i++) > { > > /* initially assume the slacks being zero */ > > o8qpdu_y[i] = zero; > } > for (i = 1 ; i <= ndual ; i++) > { > o8qpdu_xd[i] = -o8qpdu_y[i]; > } > /* unconstrained minimizer of the QP: slacks come first */ > > f = p5*o8sc1(1,ndual,o8qpdu_g0,o8qpdu_xd); > fmax = f ; > /* define the initial working set: all slacks are at their */ > /* lower bounds. iq = dimension of the QP's working set */ > > iq = nr; > > for (i = 1 ; i <= iq ; i++) > { > o8qpdu_ai[i] = i; > r[i][i] = one; > ud[i] = tauqp; > } > /* slacks are at zero, multipliers for the slacks at tauqp */ > rnorm = one; > rlow = one; > > /* introduction of equality constraints */ > /* these are characterized by a negative index in ai and never*/ > /* considered for inactivation. they also play no role for the*/ > /* dual step size t1 */ > > > for (i = 1 ; i <= me ; i++) > { > for ( j = 1 ; j <= iq ; j ++) > { > ud1[j] = ud[j]; > } > L20: > > ud1[iq+1] = zero; > > /* an equality constraint is indicated by the negative index */ > > o8qpdu_ai[iq+1] = -o8qpdu_eqlist[i] ; > l = (o8qpdu_eqlist[i]+1)/2 ; > /* store the gradient of this constraint in np */ > for (j = 1 ; j <= n ; j++) > { > if ( l > n ) > { /* sign is 1 for equality constraints */ > o8qpdu_cei[j+incr] = gres[j][l-n]; > } > else > { > o8qpdu_cei[j+incr] = zero ; > } > } > for (j = 1 ; j <= incr ; j++) > { > o8qpdu_cei[j] = zero; > } > o8qpdu_cei[i] = one; > if ( res[o8qpdu_eqlist[i]] > zero ) o8qpdu_cei[i] = -one; > if ( l <= n ) > { > o8qpdu_cei[l+incr] = xsc[l] ; /* a fixed primal variable */ > } > for (j = 1 ; j <= ndual ; j++) > { > np[j] = o8qpdu_cei[j]; > } > > o8zup(o8qpdu_z); > > if ( iq != 0 ) o8rup(o8qpdu_vr); > > /* |z| = 0? */ > > zz = o8vecn(1,ndual,o8qpdu_z); > term = o8sc1(1,ndual,o8qpdu_z,np); > > if ( zz >= tiny*rnorm && term > zero ) > { > t2 = (-o8sc1(1,ndual,np,o8qpdu_xd)-res[o8qpdu_eqlist[i]])/term; > } > else if ( (-o8sc1(1,ndual,np,o8qpdu_xd)-res[o8qpdu_eqlist[i]]) >= zero ) > { > t2 = infiny; > } > else > { > t2 = -infiny; > } > /* for an equality constraint t2 may be positive or negative*/ > > if ( iq != 0 ) o8rup(o8qpdu_vr); > l = 0; > if ( t2 > zero ) > { > t1 = infiny; > for (k = 1 ; k <= iq ; k++) > { > if( o8qpdu_vr[k] > zero && o8qpdu_ai[k] > 0 ) > { > if( ud1[k]/o8qpdu_vr[k] < t1 ) > { > t1 = ud1[k]/o8qpdu_vr[k]; > } > } > } > t = min(t1,t2); > } > else > { > t1 = infiny; > for ( k = 1 ; k <= iq ; k++) > { > if( o8qpdu_vr[k] < zero && o8qpdu_ai[k] > 0 ) > { > if( ud1[k]/fabs(o8qpdu_vr[k]) < t1 ) > { > t1 = ud1[k]/fabs(o8qpdu_vr[k]); > } > } > } > t1 = -t1; > t = max(t1,t2); > > /* t now negative */ > } > /* add constraint , otherwise we must first delete some */ > /* inequality constraint with zero multiplier */ > /* first delete then add! */ > > if ( fabs(t) >= infiny ) > { > qpterm = -2 ; > if ( ! silent ) o8msg(20); > goto L2000; > } > if ( fabs(t2) >= infiny ) > { > > /* purely dual step */ > > for (k = 1 ; k <= iq ; k++) > { > ud1[k] = ud1[k]+t*(-o8qpdu_vr[k]); > if ( ud1[k] < zero && o8qpdu_ai[k] > 0 ) ud1[k] = zero; > /* ud1[k] < 0 is a roundoff effect */ > } > ud1[iq+1] = ud1[iq+1]+t; > o8qpdu_qpdel[0] = 0; > for (j = 1 ; j <= iq ; j++) > { > if ( ud1[j] <= tiny && o8qpdu_ai[j] > 0 ) > { > o8qpdu_qpdel[0] = o8qpdu_qpdel[0]+1; > o8qpdu_qpdel[o8qpdu_qpdel[0]] = o8qpdu_ai[j]; > } > } > for (k = 1 ; k <= o8qpdu_qpdel[0] ; k++) > { > l = o8qpdu_qpdel[k]; > o8qpdu_iai[l] = l; > > o8dlcd(o8qpdu_ai,l); > } > goto L20; > } > /* primal and dual step */ > for (k = 1 ; k <= ndual ; k++) > { > o8qpdu_xd[k] = o8qpdu_xd[k]+t*o8qpdu_z[k]; > } > for (k = 1 ; k <= iq ; k++) > { > ud1[k] = ud1[k]+t*(-o8qpdu_vr[k]); > if ( ud1[k] < zero && o8qpdu_ai[k] > 0 ) ud1[k] = zero; > } > ud1[iq+1] = ud1[iq+1]+t; > > f = f+t*o8sc1(1,ndual,o8qpdu_z,np)*(p5*t+ud1[iq+1]); > if ( f <= fmax*(one+epsmac)+epsmac ) > { > nosucc += 1 ; > if ( nosucc > qpitma ) > { > qpterm = - 3; > if ( ! silent ) o8msg(26); > goto L2000; > } > } > else > { > nosucc = 0 ; > } > fmax = max( f , fmax ) ; > > if ( fabs(t2-t1) <= tiny ) > { > o8qpdu_qpdel[0] = 0; > for (j = 1 ; j <= iq ; j++) > { > if ( ud1[j] <= tiny && o8qpdu_ai[j] > 0 ) > { > o8qpdu_qpdel[0] = o8qpdu_qpdel[0]+1; > o8qpdu_qpdel[o8qpdu_qpdel[0]] = o8qpdu_ai[j]; > } > } > for (k = 1 ; k <= o8qpdu_qpdel[0] ; k++) > { > l = o8qpdu_qpdel[k]; > o8qpdu_iai[l] = l; > > o8dlcd(o8qpdu_ai,l); > } > o8qpdu_ai[iq+1] = -o8qpdu_eqlist[i] ; > > o8adcd(); > > } > else if ( t == t2 ) > { > o8qpdu_ai[iq+1] = -o8qpdu_eqlist[i] ; > > o8adcd(); > > } > else > { > o8qpdu_qpdel[0] = 0; > for (j = 1 ; j <= iq ; j++) > { > if ( ud1[j] <= tiny && o8qpdu_ai[j] > 0 ) > { > o8qpdu_qpdel[0] = o8qpdu_qpdel[0]+1; > o8qpdu_qpdel[o8qpdu_qpdel[0]] = o8qpdu_ai[j]; > } > } > for (k = 1 ; k <= o8qpdu_qpdel[0] ; k++) > { > l = o8qpdu_qpdel[k]; > o8qpdu_iai[l] = l; > o8dlcd(o8qpdu_ai,l); > } > goto L20; > } > for (j = 1 ; j <= iq ; j++) > { > ud[j] = ud1[j]; > } > /* end of loop over equality constraints */ > } > > /* set iai = k\ai */ > > for (i = 1 ; i <= mi ; i++) > { > o8qpdu_iai[i] = i; > } > /* step 1 */ > > L50: > > /* ai = QP - working set , iai[i] = 0 if i in ai */ > > for (i = 1 ; i <= iq ; i++) > { > ip = o8qpdu_ai[i]; > if ( ip > 0 ) o8qpdu_iai[ip] = 0; > } > /* s[o8qpdu_xd] = ci(trans)*o8qpdu_xd+ci0 >= 0 ? */ > > psid = zero; > > /* psid : the measure of infeasibility */ > > for (i = 1 ; i <= mi ; i++) > { > > /* iaexcl: if = 0, exclude from addition in this cycle */ > > o8qpdu_iaexcl[i] = 1; > su = zero; > k = 0 ; > /* numbers of inequality constraints: */ > /* i = 1,...,nr corresponds to the constraints v >= 0, u_a >= 0 */ > /* i = nr+1,....,mi to the regularized general inequalities */ > > if ( i > nr ) > { > /* an original inequality constraint */ > k = (o8qpdu_iqlist[i]-nr+1)/2 ; > if ( k > n ) > { /* not a bound, gradient stored in gres, sign is in gres[0][] */ > for (j = 1 ; j <= n ; j++) > { > o8qpdu_cii[j+incr] = gres[j][k-n]*gres[0][k-n]; > } > for (j = 1 ; j <= incr ; j++) > { > o8qpdu_cii[j] = zero; > } > o8qpdu_cii[me+i-incr] = one; > o8qpdu_ci0[i] = res[o8qpdu_iqlist[i]-nr]; > } > else > { > for ( j = 1 ; j <= ndual ; j++ ) > { > o8qpdu_cii[j] = zero ; > } > o8qpdu_cii[me+i-incr] = one ; > if ( (o8qpdu_iqlist[i]-nr) % 2 == 0 ) > { > o8qpdu_cii[k+incr] = -xsc[k]; > o8qpdu_ci0[i] = res[o8qpdu_iqlist[i]-nr]; > } > else > { > o8qpdu_cii[k+incr] = xsc[k]; > o8qpdu_ci0[i] = res[o8qpdu_iqlist[i]-nr]; > } > } > > } > else > { /* a QP slack variable */ > for (j = 1 ; j <= ndual ; j++) > { > o8qpdu_cii[j] = zero; > } > o8qpdu_ci0[i] = zero; > o8qpdu_cii[i] = one; > } >/* end computation of data of mi-th QP constraint data */ > > su = o8sc1(1,ndual,o8qpdu_cii,o8qpdu_xd)+o8qpdu_ci0[i]; > o8qpdu_s[i] = su; > psid = psid+min(zero,su); > } > > for (i = 1 ; i <= iq ; i++) { > o8qpdu_udold[i] = ud[i]; > o8qpdu_aiold[i] = o8qpdu_ai[i]; > } > for (i = 1 ; i <= ndual ; i++) { > o8qpdu_xdold[i] = o8qpdu_xd[i]; > } > L60: > ss = zero; > ip = 0; > > /* introduce most violated inequality constraint */ > > for (i = 1 ; i <= mi ; i++) { > if( o8qpdu_s[i] < ss && o8qpdu_iai[i] != 0 && o8qpdu_iaexcl[i] != 0 ) { > ss = o8qpdu_s[i]; > ip = i; > } > } > if ( iq > 1 ) { > condr = rnorm/rlow; > } else { > condr = one; > } > > if ( fabs(psid) <= tiny*(c1*c2+condr) || ip == 0) { > > /* successful termination of QP-solver for current tauqp */ > > if ( fabs(psid) > tiny*(c1*c2+condr) && ! silent ) o8msg(10); > qpterm = 1; > accinf[itstep][30] = one; > accinf[itstep][13] = condr; > accinf[itstep][14] = c1*c2; > for (i = 1 ; i <= n ; i++) { > d[i] = o8qpdu_xd[i+incr]; > } > /* new : ddnorm added */ > > ddnorm = o8vecn(1,n,d); > infeas = zero; > > for (i = 1 ; i <= incr ; i++) { > infeas = infeas+fabs(o8qpdu_xd[i]); > } > /* L1-norm of slack variables */ > > accinf[itstep][31] = tauqp; > accinf[itstep][32] = infeas; > wlow = FFALSE; > su1 = zero; > su2 = zero; > for (i = 1 ; i <= iq ; i++) > { > if ( o8qpdu_ai[i] < 0 ) > { > u[-o8qpdu_ai[i]] = ud[i]; > } > else > { > if ( o8qpdu_ai[i] > nr ) u[o8qpdu_iqlist[o8qpdu_ai[i]]-nr] = ud[i]; > } > } > > term1 = zero; > for (j = 1 ; j <= n ; j++) { > np[j] = gradf[j]*scf; > } > for (i = 1 ; i <= nres ; i++) > { > if ( i > n ) > { > for (j = 1 ; j <= n ; j++) > { > np[j] = np[j]-gres[j][i-n]*(u[2*i-1]-u[2*i]); > } > } > else > { > np[i] = np[i]-xsc[i]*(u[2*i-1]-u[2*i]); > } > } > > b2n = o8vecn(1,n,np); > > if ( scf != zero ) b2n = b2n/scf; > > /* correction in the original variables */ > > infe1 = zero; > for (i = 1 ; i <= nr ; i++) { > infe1 = infe1+fabs(o8qpdu_xd[i])*o8qpdu_mult[i]; > } > if ( upsi <= delmin*nres > && b2n <= (gfn+one)*epsx*tp2 && phase >= 0 > && infeas <= delmin*nres ) { > /* dual feasibility is automatic here */ > /* since multipliers may be incorrect for infeas != zero be careful */ > /* we consider the problem as successfully solved with reduced */ > /* requirements. we terminate here */ > > for (i = 1 ; i <= n ; i++) > { > d[i] = zero; > } > ddnorm = zero; > optite = three; > dirder = zero ; > qpterm = 1 ; > return; > } > /* there may be an additional increase of tauqp necessary again */ > if ( infe1 > (one-delta1/tauqp)*upsi && > ( o8vecn(1,n,d) <= min(infe1,pow(infe1,2))*tp1 > || upsi > tau0*p5 ) ) > { > > /* further increase tauqp ! */ > > for (i = 1 ; i <= 2*nres ; i++) > { > u[i] = zero; > slack[i] = zero; > } > if ( ! silent ) o8msg(17); > if ( tauqp*taufac > taumax ) { > if ( ! silent ) o8msg(5); > qpterm = -1; > accinf[itstep][30] = qpterm; > accinf[itstep][31] = tauqp; > accinf[itstep][32] = infeas; > for (i = 1 ; i <= n ; i++) { > d[i] = zero; > } > ddnorm = zero; > > return; > > } else { > tauqp = tauqp*taufac; > > goto L10; > } > } > /* compute new weights for the penalty-function */ > > for (i = 1 ; i <= 2*nres ; i++) > { > slack[i] = zero; > } > for (i = 1 ; i <= nr ; i++) > { > slack[aalist[i]] = o8qpdu_xd[i]; > } > wlow = FFALSE; > for (i = 1 ; i <= nres ; i++) > { > w1[2*i-1] = w[2*i-1]; > w1[2*i] = w[2*i] ; > if ( low[i] == up[i] ) > { > if ( fabs(slack[2*i-1]) > fabs(res[2*i-1])+tiny ) > { > w1[2*i-1] = fabs(u[2*i-1]); > } > else > { > w1[2*i-1] = ny*fabs(u[2*i-1])+tau; > } > } > else > { > if ( o8bind[2*i-1] == 0 ) > { > w1[2*i-1] = max(w[2*i-1]*p8,tau); > } > else > { > if ( res[2*i-1] >= zero && slack[2*i-1] <= tiny ) > w1[2*i-1] = > max(ny*fabs(u[2*i-1])+tau,(fabs(u[2*i-1])+w1[2*i-1])*p5); > if ( res[2*i-1] >= zero && slack[2*i-1] > tiny ) > w1[2*i-1] = fabs(u[2*i-1]); > if ( res[2*i-1] < zero && slack[2*i-1] <= -res[2*i-1]+tiny ) > w1[2*i-1] = > max(ny*fabs(u[2*i-1])+tau,(fabs(u[2*i-1])+w1[2*i-1])*p5); > if ( res[2*i-1] < zero && slack[2*i-1] > -res[2*i-1]+tiny ) > w1[2*i-1] = fabs(u[2*i-1]); > } > if ( o8bind[2*i] == 0 ) > { > w1[2*i] = max(w[2*i]*p8,tau); > } > else > { > if ( res[2*i] >= zero && slack[2*i] <= tiny ) > w1[2*i] = > max(ny*fabs(u[2*i])+tau,(fabs(u[2*i])+w1[2*i])*p5); > if ( res[2*i] >= zero && slack[2*i] > tiny ) > w1[2*i] = fabs(u[2*i]); > if ( res[2*i] < zero && slack[2*i] <= -res[2*i]+tiny ) > w1[2*i] = > max(ny*fabs(u[2*i])+tau,(fabs(u[2*i])+w1[2*i])*p5); > if ( res[2*i] < zero && slack[2*i] > -res[2*i]+tiny ) > w1[2*i] = fabs(u[2*i]); > } > } > if ( w1[2*i-1] < w[2*i-1] || w1[2*i] < w[2*i] ) wlow = TTRUE; > } > if ( wlow ) > { > s1 = zero; > s2 = zero; > for (i = 1 ; i <= nres ; i++) > { > if ( low[i] == up[i] ) > { > s1 += w1[2*i-1]*fabs(resst[2*i-1]); > s2 += w1[2*i-1]*fabs(res[2*i-1]); > } > else > { > s1 -= min(zero,resst[2*i-1])*w1[2*i-1]; > s2 -= min(zero,res[2*i-1] )*w1[2*i-1]; > s1 -= min(zero,resst[2*i])*w1[2*i]; > s2 -= min(zero,res[2*i] )*w1[2*i]; > > } > } > diff0 = (fxst-fx)*scf+(s1-s2); > if ( diff0 >= eta*clow && itstep-lastdw >= max(5,min(n/10,20)) ) { > > /* accept new (diminished ) weights */ > > lastdw = itstep; > lastch = itstep; > level = diff0/iterma; > psist = s1; > psi = s2; > for (i = 1 ; i <= 2*nres ; i++) { > if ( w1[i] != w[i] ) lastch = itstep; > w[i] = w1[i]; > } > clow = clow+one; > if ( clow > itstep/10 ) { > > /* additional increase of eta */ > > eta = eta*onep3; > if ( ! silent ) o8info(11); > } > if ( ! silent ) o8info(12); > > goto L1000; > } > } > /* we cannot accept new weights */ > /* reset weights */ > > for (i = 1 ; i <= nres ; i++ ) > { > w1[2*i-1] = w[2*i-1]; > w1[2*i ] = w[2*i] ; > if ( low[i] == up[i] ) > { > if ( slack[2*i-1] > fabs(res[2*i-1]) ) > w1[2*i-1] = fabs(u[2*i-1]); > if ( slack[2*i-1] <= fabs(res[2*i-1]) ) > { > if ( w[2*i-1] <= fabs(u[2*i-1]) > && fabs(u[2*i-1]) <= w[2*i-1]+tau ) > { > w1[2*i-1] = w[2*i-1]+two*tau; > } > else > { > w1[2*i-1] = max(w[2*i-1],ny*fabs(u[2*i-1])+tau); > } > } > } > else > { > if ( slack[2*i-1] > -min(-tiny,res[2*i-1]) && o8bind[2*i-1] == 1) > { > w1[2*i-1] = fabs(u[2*i-1]); > } > else if( o8bind[2*i-1] == 1 && > slack[2*i-1] <= -min(-tiny,res[2*i-1]) > && u[2*i-1] <= w[2*i-1]+tau && w[2*i-1] >= u[2*i-1] ) > { > > w1[2*i-1] = w[2*i-1]+two*tau; > } > else if ( o8bind[2*i-1] == 1 ) > { > w1[2*i-1] = max(w[2*i-1],ny*fabs(u[2*i-1])+tau); > } > if ( slack[2*i] > -min(-tiny,res[2*i]) && o8bind[2*i] == 1 ) > { > w1[2*i] = fabs(u[2*i]); > } > else if ( o8bind[2*i] == 1 && slack[2*i] <= -min(-tiny,res[2*i]) > && u[2*i] <= w[2*i]+tau && w[2*i] >= u[2*i] ) > { > > w1[2*i] = w[2*i]+two*tau; > } > else if ( o8bind[2*i] == 1 ) > { > w1[2*i] = max(w[2*i],ny*fabs(u[2*i])+tau); > } > > } > } > term1 = zero; > for (i = 1 ; i <= 2*nres ; i++) { > if ( w1[i] > w[i] || w1[i] < w[i] ) lastch = itstep; > if ( w1[i] > w[i] ) lastup = itstep; > if ( w1[i] < w[i] ) lastdw = itstep; > w[i] = w1[i]; > term1 = max(term1,w[i]); > } > s1 = zero; > s2 = zero; > for (i = 1 ; i <= nres ; i++) { > if ( low[i] == up[i] ) > { > s1 += w[2*i-1]*fabs(resst[2*i-1]); > s2 += w[2*i-1]*fabs(res[2*i-1]); > } else { > s1 -= w[2*i-1]*min(zero,resst[2*i-1]); > s2 -= w[2*i-1]*min(zero,res[2*i-1]); > s1 -= w[2*i]*min(zero,resst[2*i]); > s2 -= w[2*i]*min(zero,res[2*i]); > > } > } > psist = s1; > psi = s2; > if ( ! silent ) o8info(12); > accinf[itstep][20] = term1; > accinf[itstep][19] = clow; > > goto L1000; > } > /************************************************************/ > /** continue QP solver **/ > /************************************************************/ > > if ( ip > nr ) > { > /* compute gradient of current QP inequality constraint */ > /* this is a linearized original constraint */ > k = (o8qpdu_iqlist[ip]-nr+1)/2 ; > /* k is in {1,...,nres} */ > if ( k > n ) > { /* a general constraint */ > for (j = 1 ; j <= n ; j++) > { > o8qpdu_cii[j+incr] = gres[j][k-n]*gres[0][k-n]; > } > for (j = 1 ; j <= incr ; j++) > { > o8qpdu_cii[j] = zero; > } > o8qpdu_cii[me+ip-nr] = one; > o8qpdu_ci0[ip] = res[o8qpdu_iqlist[ip]-nr]; > } > else > { /* an active bound constraint */ > for ( j = 1 ; j <= ndual ; j++ ) > { > o8qpdu_cii[j] = zero ; > } > o8qpdu_cii[me+ip-nr] = one ; > o8qpdu_ci0[ip ] = res[o8qpdu_iqlist[ip]-nr]; > o8qpdu_cii[k+incr ] = ( (o8qpdu_iqlist[ip]-nr) %2 ==0 ) ? -xsc[k] : xsc[k] ; > } > > } > else > { > for (j = 1 ; j <= ndual ; j++) > { > o8qpdu_cii[j] = zero; > } > o8qpdu_ci0[ip] = zero; > o8qpdu_cii[ip] = one; > } > for (i = 1 ; i <= ndual ; i++) { > np[i] = o8qpdu_cii[i]; > } > > for (i = 1 ; i <= iq ; i++) { > ud1[i] = ud[i]; > } > ud1[iq+1] = zero; > o8qpdu_ai[iq+1] = ip; /* this is in {1,..,mi} */ > > L100: > > /* step 2a */ > > o8zup(o8qpdu_z); > > if ( iq != 0 ) o8rup(o8qpdu_vr); > > l = 0; > t1 = infiny; > for (k = 1 ; k <= iq ; k++) { > if(o8qpdu_ai[k] > 0 && o8qpdu_vr[k] > zero) { > if ( ud1[k]/o8qpdu_vr[k] < t1 ) { > t1 = ud1[k]/o8qpdu_vr[k]; > } > } > } > /* |z| = 0? */ > > /* old zz = o8sc1(1,ndual,o8qpdu_z,o8qpdu_z) */ > > zz = o8vecn(1,ndual,o8qpdu_z); > term = o8sc1(1,ndual,o8qpdu_z,np); > > /* old if (zz != zero && term > zero ) { */ > > if ( zz >= tiny*rnorm && term > zero ) { > t2 = -o8qpdu_s[ip]/term; > } else { > t2 = infiny; > } > t = min(t1,t2); > > if ( t >= infiny ) > { > qpterm = -2 ; > if ( ! silent ) o8msg(20); > goto L2000; > } > > > if( t2 >= infiny ) { > for (k = 1 ; k <= iq ; k++) { > ud1[k] = ud1[k]+t*(-o8qpdu_vr[k]); > if ( ud1[k] < zero && o8qpdu_ai[k] > 0 ) ud1[k] = zero; > } > ud1[iq+1] = ud1[iq+1]+t; > o8qpdu_qpdel[0] = 0; > for (i = 1 ; i <= iq ; i++) { > if ( ud1[i] <= tiny && o8qpdu_ai[i] > 0 ) { > o8qpdu_qpdel[0] = o8qpdu_qpdel[0]+1; > o8qpdu_qpdel[o8qpdu_qpdel[0]] = o8qpdu_ai[i]; > } > } > for (k = 1 ; k <= o8qpdu_qpdel[0] ; k++) { > l = o8qpdu_qpdel[k]; > o8qpdu_iai[l] = l; > > o8dlcd(o8qpdu_ai,l); > } > goto L100; > } > for (k = 1 ; k <= ndual ; k++) { > o8qpdu_xd[k] = o8qpdu_xd[k]+t*o8qpdu_z[k]; > } > for (k = 1 ; k <= iq ; k++) { > ud1[k] = ud1[k]+t*(-o8qpdu_vr[k]); > if ( ud1[k] < zero && o8qpdu_ai[k] > 0 ) ud1[k] = zero; > } > ud1[iq+1] = ud1[iq+1]+t; > > f = f+t*o8sc1(1,ndual,o8qpdu_z,np)*(p5*t+ud1[iq+1]); > if ( f <= fmax*(one+epsmac)+epsmac ) > { > nosucc += 1 ; > if ( nosucc > qpitma ) > { > qpterm = - 3; > if ( ! silent ) o8msg(26); > goto L2000; > } > } > else > { > nosucc = 0 ; > } > fmax = max( f , fmax ) ; > > if ( t2 <= t1-tiny ) { > > /* ddual is computed by o8zup */ > > if ( o8vecn(iq+1,ndual,ddual) < epsmac*rnorm ) { > > /* degeneracy: adding this constraint gives a singular working set */ > /* theoretically impossible, but due to roundoff this may occur. */ > /* mark this constraint and try to add another one */ > > iptr = ip; > iqtr = iq; > for (i = 1 ; i <= iq ; i++) { > aitr[i] = o8qpdu_ai[i]; > } > sstr = ss; > riitr = o8vecn(iq+1,ndual,ddual); > > if ( ! silent ) o8msg(19); > o8qpdu_iaexcl[ip] = 0; > for (i = 1 ; i <= mi ; i++) { > o8qpdu_iai[i] = i; > } > for (i = 1 ; i <= iq ; i++) { > o8qpdu_ai[i] = o8qpdu_aiold[i]; > if (o8qpdu_ai[i] > 0 ) o8qpdu_iai[o8qpdu_ai[i]] = 0; > ud1[i] = o8qpdu_udold[i]; > } > for (i = 1 ; i <= ndual ; i++) { > o8qpdu_xd[i] = o8qpdu_xdold[i]; > } > goto L60; > } > /* add constraint, l-pair */ > > o8adcd(); > > o8qpdu_iai[ip] = 0; > for (i = 1 ; i <= iq ; i++) { > ud[i] = ud1[i]; > } > goto L50; > } > su = zero; > if ( ip > nr ) > { > > /* an original linearized inequality constraint */ > k = (o8qpdu_iqlist[ip]-nr+1)/2 ; > /* k is in {1,...,nres} */ > if ( k > n ) > { /* a general constraint */ > for (j = 1 ; j <= n ; j++) > { > o8qpdu_cii[j+incr] = gres[j][k-n]*gres[0][k-n]; > } > for (j = 1 ; j <= incr ; j++) > { > o8qpdu_cii[j] = zero; > } > o8qpdu_cii[me+ip-nr] = one; > o8qpdu_ci0[ip] = res[o8qpdu_iqlist[ip]-nr]; > } > else > { /* an active bound constraint */ > for ( j = 1 ; j <= ndual ; j++ ) > { > o8qpdu_cii[j] = zero ; > } > o8qpdu_cii[me+ip-nr] = one ; > o8qpdu_ci0[ip ] = res[o8qpdu_iqlist[ip]-nr]; > o8qpdu_cii[k+incr ] = ( (o8qpdu_iqlist[ip]-nr) %2 ==0 ) ? -xsc[k] : xsc[k] ; > } > o8qpdu_s[ip] = o8sc1(1,ndual,o8qpdu_cii,o8qpdu_xd)+o8qpdu_ci0[ip]; > } > else > { > /* a slack constraint */ > > o8qpdu_s[ip] = o8qpdu_xd[ip]; > } > /* now t = t1 */ > > o8qpdu_qpdel[0] = 0; > for (i = 1 ; i <= iq ; i++) { > if ( ud1[i] <= tiny && o8qpdu_ai[i] > 0 ) { > o8qpdu_qpdel[0] = o8qpdu_qpdel[0]+1; > o8qpdu_qpdel[o8qpdu_qpdel[0]] = o8qpdu_ai[i]; > } > } > for (k = 1 ; k <= o8qpdu_qpdel[0] ; k++) { > l = o8qpdu_qpdel[k]; > o8qpdu_iai[l] = l; > > o8dlcd(o8qpdu_ai,l); > } > if ( t2 <= t1+tiny ) { > if ( o8vecn(iq+1,ndual,ddual) < epsmac*rnorm ) { > > /* degeneracy */ > > iptr = ip; > iqtr = iq; > for (i = 1 ; i <= iq ; i++) { > aitr[i] = o8qpdu_ai[i]; > } > sstr = ss; > riitr = o8vecn(iq+1,ndual,ddual); > > if ( ! silent ) o8msg(19); > o8qpdu_iaexcl[ip] = 0; > for (i = 1 ; i <= mi ; i++) { > o8qpdu_iai[i] = i; > } > for (i = 1 ; i <= iq ; i++) { > o8qpdu_ai[i] = o8qpdu_aiold[i]; > if ( o8qpdu_ai[i] > 0 ) o8qpdu_iai[o8qpdu_ai[i]] = 0; > ud1[i] = o8qpdu_udold[i]; > } > for (i = 1 ; i <= ndual ; i++) { > o8qpdu_xd[i] = o8qpdu_xdold[i]; > } > goto L60; > } > /* add constraint, l-pair */ > > o8adcd(); > > o8qpdu_iai[ip] = 0; > for (i = 1 ; i <= iq ; i++) { > ud[i] = ud1[i]; > } > goto L50; > > } else { > > goto L100; > } > /* this is the exit point of o8qpdu */ > /* we either may have successful or unsuccessful termination here */ > /* the latter with qpterm = -2 or -3, in which case it may nevertheless */ > /* be possible to use the computed d. -2 or -3 exit is theoretically */ > /* impossible but may occur due to roundoff effects. */ > /* we check the directional derivative of the penalty-function now */ > > L1000: > > /* cut and rescale d if appropriate */ > > o8cutd(); > > /* compute the directional derivative dirder */ > > o8dird(); > > if ( dirder >= zero || ( -dirder <= epsmac*tp2*(scf*fabs(fx)+psi+one) && > infeas > max(upsi,nres*delmin) ) ) { > if ( ! silent ) o8msg(18); > if ( tauqp <= taumax/taufac ) { > tauqp = tauqp*taufac; > > goto L10; > > } else { > if ( ! silent ) o8msg(5); > qpterm = -1; > accinf[itstep][30] = qpterm; > accinf[itstep][31] = tauqp; > accinf[itstep][32] = infeas; > for (i = 1 ; i <= n ; i++) { > d[i] = zero; > } > ddnorm = zero; > } > } > return; > > L2000: > > /* QP infeasible ( in this application impossible , theoretically) */ > > > accinf[itstep][30] = -two; > accinf[itstep][13] = condr; > accinf[itstep][14] = c1*c2; > for (i = 1 ; i <= n ; i++) { > d[i] = o8qpdu_xd[i+incr]; > } > ddnorm = o8vecn(1,n,d); > su1 = zero; > for (i = 1 ; i <= incr ; i++) { > su1 = su1+fabs(o8qpdu_xd[i]); > } > /* L1-norm of slack variables */ > > accinf[itstep][32] = su1; > wlow = FFALSE; > su1 = zero; > su2 = zero; > for (i = 1 ; i <= iq ; i++) { > if ( o8qpdu_ai[i] < 0 ) { > u[-o8qpdu_ai[i]] = ud[i]; > } else { > if ( o8qpdu_ai[i] > nr ) u[o8qpdu_iqlist[o8qpdu_ai[i]]-nr] = ud[i]; > /* o8qpdu_iqlist[ai[i]] = aalist[k]+nr for some k */ > } > } > /* compute Lagrangian violation */ > term1 = zero; > for (j = 1 ; j <= n ; j++) { > np[j] = gradf[j]*scf; > } > for (i = 1 ; i <= nres ; i++) > { > if ( i > n ) > { > for (j = 1 ; j <= n ; j++) > { > np[j] = np[j]-gres[j][i-n]*(u[2*i-1]-u[2*i]); > } > } > else > { > np[i] -= xsc[i]*(u[2*i-1]-u[2*i]) ; > } > w1[2*i-1] = w[2*i-1]; > w1[2*i ] = w[2*i] ; > if ( low[i] == up[i] ) > { > if ( slack[2*i-1] > fabs(res[2*i-1]) ) w1[2*i-1] = fabs(u[2*i-1]); > if ( slack[2*i-1] <= fabs(res[2*i-1]) ) > { > if ( w[2*i-1] <= fabs(u[2*i-1]) > && fabs(u[2*i-1]) <= w[2*i-1]+tau ) > { > w1[2*i-1] = w[2*i-1]+two*tau; > } > else > { > w1[2*i-1] = max(w[2*i-1],ny*fabs(u[2*i-1])+tau); > } > } > su1 += fabs(res[2*i-1] )*w1[2*i-1]; > su2 += fabs(resst[2*i-1])*w1[2*i-1]; > } > else > { > if ( slack[2*i-1] > -min(-tiny,res[2*i-1]) && o8bind[2*i-1] == 1 ) > { > w1[2*i-1] = fabs(u[2*i-1]); > } > else if ( o8bind[2*i-1] == 1 && slack[2*i-1] <= -min(-tiny,res[2*i-1]) > && u[2*i-1] <= w[2*i-1]+tau && w[2*i-1] >= u[2*i-1] ) > { > w1[2*i-1] = w[2*i-1]+two*tau; > } > else if ( o8bind[2*i-1] == 1 ) > { > w1[2*i-1] = max(w[2*i-1],ny*fabs(u[2*i-1])+tau); > } > su1 -= w1[2*i-1]*min(zero,res[2*i-1]); > su2 -= w1[2*i-1]*min(zero,resst[2*i-1]); > > if ( slack[2*i] > -min(-tiny,res[2*i]) && o8bind[2*i] == 1 ) > { > w1[2*i] = fabs(u[2*i]); > } > else if ( o8bind[2*i] == 1 && slack[2*i] <= -min(-tiny,res[2*i]) > && u[2*i] <= w[2*i]+tau && w[2*i] >= u[2*i] ) > { > w1[2*i] = w[2*i]+two*tau; > } > else if ( o8bind[2*i] == 1 ) > { > w1[2*i] = max(w[2*i],ny*fabs(u[2*i])+tau); > } > su1 -= w1[2*i]*min(zero,res[2*i]); > su2 -= w1[2*i]*min(zero,resst[2*i]); > > } > if ( w[i] != w1[i] ) lastch = itstep; > w[i] = w1[i]; > term1 = max(term1,w[i]); > } > psist = su2; > psi = su1; > > b2n = one; > > if ( scf != zero ) b2n = b2n/scf; > if ( wlow ) { > clow = clow+one; > lastch = itstep; > lastdw = itstep; > } > if ( ! silent ) o8info(12); > accinf[itstep][19] = clow; > accinf[itstep][20] = term1; > accinf[itstep][31] = tauqp; > if ( upsi <= delmin*nres > && b2n <= (gfn+one)*epsx*tp2 && phase >= 0 > && infeas <= delmin*nres ) > { > > /* since multipliers may be incorrect for infeas != zero be careful */ > /* we consider the problem as successfully solved with reduced */ > /* requirements */ > > for (i = 1 ; i <= n ; i++) { > d[i] = zero; > } > ddnorm = zero; > optite = three; > qpterm = 1 ; /* new */ > return; > } > > goto L1000; >} > >/* **************************************************************************** */ >/* compute updated projected gradient (primal) */ >/* **************************************************************************** */ >void o8zup(DDOUBLE z[]) { >} > >/* **************************************************************************** */ >/* compute correction of dual multipliers */ >/* **************************************************************************** */ >void o8rup(DDOUBLE rv[]) { >} > >/* **************************************************************************** */ >/* delete constraint nr. l */ >/* **************************************************************************** */ >void o8dlcd(IINTEGER ai[],IINTEGER l) { > return; >} > >/* **************************************************************************** */ >/* add constraint whose gradient is given by np */ >/* **************************************************************************** */ >void o8adcd(void) { > return; >} > >/* **************************************************************************** */ >/* computes the inverse of the upper triangular matrix part */ >/* **************************************************************************** */ >void o8rinv(IINTEGER n,DDOUBLE **a,IINTEGER ndual,DDOUBLE **donlp2_x) { > return; >} > >/* **************************************************************************** */ >/* suite of function interfaces, for performing scaling and */ >/* external evaluations */ >/* if bloc = TTRUE then it is assumed that prior to calls to es<xyz> */ >/* valid function information is computed (via a call of user_eval) */ >/* and stored in fu and fugrad , setting valid = TTRUE afterwards */ >/* **************************************************************************** */ > >/* **************************************************************************** */ >/* objective function */ >/* **************************************************************************** */ >void esf(DDOUBLE donlp2_x[],DDOUBLE *fx) >{ > return; >} > >/* **************************************************************************** */ >/* gradient of objective function */ >/* **************************************************************************** */ >void esgradf(DDOUBLE donlp2_x[],DDOUBLE gradf[]) >{ > return; >} > >/********************************************************************/ > >void escon( IINTEGER call_type , IINTEGER liste[], DDOUBLE donlp2_x[], DDOUBLE constr[], > LLOGICAL errliste[]) >{ > return; >} >void escongrad(IINTEGER liste[], IINTEGER shift , DDOUBLE donlp2_x[], > DDOUBLE **grad_constr) >/********************************************************************/ >/* evaluate gradients of nonlinear constraints from liste for the */ >/* internally scaled donlp2_x */ >/********************************************************************/ >{ > return; >} >/* **************************************************************************** */ >/* suite of functions to allocate and free memory for dynamic memory allocation */ >/* **************************************************************************** */ > > >/**********************************************************************/ >/*** allocate memory for integer 1D array *******/ >/**********************************************************************/ > >IINTEGER* i1_malloc(IINTEGER size1, IINTEGER init) >{ > > /* Allocate IINTEGER precision array with length "size1" */ > /* assuming zero origin. If initialize is non zero */ > /* then initize the allocated array to zero. */ > > IINTEGER* array; > IINTEGER i; > > array = (IINTEGER*) malloc((size_t) (size1+2*ABC)*sizeof(IINTEGER)); > if (!array) { > fprintf(stderr, "ERROR: i1_malloc: memory error: malloc failed"); > exit(-1); > } > > /* initialize the array to 0 */ > if (init) { > for (i=0; i<size1; i++) > array[i] = 0; > } > > return array; >} > > >/**********************************************************************/ >/*** free memory for IINTEGER 1D array *******/ >/**********************************************************************/ >void i1_free(IINTEGER* array) >{ > /* Check for null pointer */ > if (!array) { > fprintf(stderr, "ERROR: i1_free: memory error: pointer is null"); > exit(-1); > } > > /* free the memory for the IINTEGER 1D array */ > free(array-ABC); >} > >/**********************************************************************/ >/*** allocate memory for IINTEGER 2D array *******/ >/**********************************************************************/ >IINTEGER** i2_malloc(IINTEGER size1, IINTEGER size2, IINTEGER init) >{ > > /* Allocate IINTEGER precision array with lengths "size1" and */ > /* size2 assuming zero origin. If initialize is non zero */ > /* then initize the allocated array to zero. */ > > IINTEGER** array; > IINTEGER* arraytemp; > IINTEGER i,j; > > array = (IINTEGER**) malloc((size_t) size1*sizeof(IINTEGER*)); > if (!array) { > fprintf(stderr, "ERROR: d2_malloc: memory error: malloc failed"); > exit(-1); > } > for (i=0; i<size1; i++) { > arraytemp = (IINTEGER*) malloc((size_t) (size2+2*ABC)*sizeof(IINTEGER)); > if (!arraytemp) { > fprintf(stderr, "ERROR: d2_malloc: memory error: malloc failed"); > exit(-1); > } > array[i] = arraytemp + ABC; > } > > if (init) { > for (i=0; i<size1; i++) > for (j=0; j<size2; j++) > array[i][j] = 0.0; > } > > return array; >} > >/**********************************************************************/ >/*** free memory for IINTEGER 2D array *******/ >/**********************************************************************/ >void i2_free(IINTEGER** array, IINTEGER size1) >{ > IINTEGER i; > > /* Check for null pointer */ > if (!array) { > fprintf(stderr, "ERROR: d2_free: memory error: pointer is null"); > exit(-1); > } > > /* free the memory for the IINTEGER 2D array piece by piece */ > for (i=0; i<size1; i++) > free(array[i]-ABC); > free(array); >} > > >/**********************************************************************/ >/*** allocate memory for double 1D array *******/ >/**********************************************************************/ >DDOUBLE* d1_malloc(IINTEGER size1, IINTEGER init) >{ > > /* Allocate DDOUBLE precision array with length "size1" */ > /* assuming zero origin. If initialize is non zero */ > /* then initize the allocated array to zero. */ > > DDOUBLE* array; > IINTEGER i; > > array = (DDOUBLE*) malloc((size_t) (size1+2*ABC)*sizeof(DDOUBLE)); > if (!array) { > fprintf(stderr, "ERROR: d1_malloc: memory error: malloc failed"); > exit(-1); > } > > if (init) { > for (i=0; i<size1; i++) > array[i] = 0.0; > } > > return array; >} > >/**********************************************************************/ >/*** free memory for DDOUBLE 1D array *******/ >/**********************************************************************/ >void d1_free(DDOUBLE* array) >{ > /* Check for null pointer */ > if (!array) { > fprintf(stderr, "ERROR: d1_free: memory error: pointer is null"); > exit(-1); > } > > /* free the memory for the DDOUBLE 1D array */ > free(array-ABC); >} > >/**********************************************************************/ >/*** allocate memory for DDOUBLE 2D array *******/ >/**********************************************************************/ >DDOUBLE** d2_malloc(IINTEGER size1, IINTEGER size2, IINTEGER init) >{ > > /* Allocate DDOUBLE precision array with lengths "size1" and */ > /* size2 assuming zero origin. If initialize is non zero */ > /* then initize the allocated array to zero. */ > > DDOUBLE** array; > DDOUBLE* arraytemp; > IINTEGER i,j; > > array = (DDOUBLE**) malloc((size_t) size1*sizeof(DDOUBLE*)); > if (!array) { > fprintf(stderr, "ERROR: d2_malloc: memory error: malloc failed"); > exit(-1); > } > for (i=0; i<size1; i++) { > > /* Allocate memory */ > arraytemp = (DDOUBLE*) malloc((size_t) (size2+2*ABC)*sizeof(DDOUBLE)); > if (!arraytemp) { > fprintf(stderr, "ERROR: d2_malloc: memory error: malloc failed"); > exit(-1); > } > > array[i] = arraytemp + ABC; > } > > if (init) { > for (i=0; i<size1; i++) > for (j=0; j<size2; j++) > array[i][j] = 0.0; > } > > return array; >} > >/**********************************************************************/ >/*** free memory for DDOUBLE 2D array *******/ >/**********************************************************************/ >void d2_free(DDOUBLE** array, IINTEGER size1) >{ > IINTEGER i; > > /* Check for null pointer */ > if (!array) { > fprintf(stderr, "ERROR: d2_free: memory error: pointer is null"); > exit(-1); > } > > /* free the memory for the DDOUBLE 2D array piece by piece */ > for (i=0; i<size1; i++) > free(array[i]); > free(array-ABC); >} > >/**********************************************************************/ >/*** allocate memory for LLOGICAL 1D array *******/ >/**********************************************************************/ >LLOGICAL* l1_malloc(IINTEGER size1, IINTEGER init) >{ > > /* Allocate LLOGICAL precision array with length "size1" */ > /* assuming zero origin. If initialize is non zero */ > /* then initize the allocated array to zero. */ > > LLOGICAL* array; > IINTEGER i; > > array = (LLOGICAL*) malloc((size_t) (size1+2*ABC)*sizeof(LLOGICAL)); > if (!array) { > fprintf(stderr, "ERROR: l1_malloc: memory error: malloc failed"); > exit(-1); > } > > if (init) { > for (i=0; i<size1; i++) > array[i] = 0.0; > } > > return array; >} > >/**********************************************************************/ >/*** free memory for LLOGICAL 1D array *******/ >/**********************************************************************/ >void l1_free(LLOGICAL* array) >{ > /* Check for null pointer */ > if (!array) { > fprintf(stderr, "ERROR: l1_free: memory error: pointer is null"); > exit(-1); > } > > /* free the memory for the LLOGICAL 1D array */ > free(array-ABC); >} > >/**********************************************************************/ >/*** allocate memory for LLOGICAL 2D array *******/ >/**********************************************************************/ >LLOGICAL** l2_malloc(IINTEGER size1, IINTEGER size2, IINTEGER init) >{ > > /* Allocate LLOGICAL precision array with lengths "size1" and */ > /* size2 assuming zero origin. If initialize is non zero */ > /* then initize the allocated array to zero. */ > > LLOGICAL** array; > LLOGICAL* arraytemp; > IINTEGER i,j; > > array = (LLOGICAL**) malloc((size_t) size1*sizeof(LLOGICAL*)); > if (!array) { > fprintf(stderr, "ERROR: l2_malloc: memory error: malloc failed"); > exit(-1); > } > for (i=0; i<size1; i++) { > arraytemp = (LLOGICAL*) malloc((size_t) (size2+2*ABC)*sizeof(LLOGICAL)); > if (!arraytemp) { > fprintf(stderr, "ERROR: l2_malloc: memory error: malloc failed"); > exit(-1); > } > > array[i] = arraytemp + ABC; > } > > if (init) { > for (i=0; i<size1; i++) > for (j=0; j<size2; j++) > array[i][j] = 0.0; > } > > return array; >} > >/**********************************************************************/ >/*** free memory for LLOGICAL 2D array *******/ >/**********************************************************************/ >void l2_free(LLOGICAL** array, IINTEGER size1) >{ > IINTEGER i; > > /* Check for null pointer */ > if (!array) { > fprintf(stderr, "ERROR: l2_free: memory error: pointer is null"); > exit(-1); > } > > /* free the memory for the LLOGICAL 2D array piece by piece */ > for (i=0; i<size1; i++) > free(array[i]-ABC); > free(array); >} > > >/**********************************************************************/ >/*** allocate memory for global arrays *******/ >/**********************************************************************/ >void global_mem_malloc() { > > /* o8comm.h */ > accinf = d2_malloc(iterma+1, 33, 1); > donlp2_x = d1_malloc(n+1, 1); > x0 = d1_malloc(n+1, 1); > x1 = d1_malloc(n+1, 1); > xmin = d1_malloc(n+1, 1); > d = d1_malloc(n+1, 1); > d0 = d1_malloc(n+1, 1); > dd = d1_malloc(n+1, 1); > difx = d1_malloc(n+1, 1); > resmin = d1_malloc(2*(n+nlin+nonlin+1), 1); > gradf = d1_malloc(n+1, 1); > qgf = d1_malloc(n+1, 1); > gphi0 = d1_malloc(n+1, 1); > gphi1 = d1_malloc(n+1, 1); > gres = d2_malloc(n+1, nlin+nonlin+1, 1); > gresn = d1_malloc(nlin+nonlin+1, 1); > perm = i1_malloc(n+1, 1); > perm1 = i1_malloc(n+1, 1); > colno = i1_malloc(2*(n+nlin+nonlin)+1, 1); > qr = d2_malloc(n+1, n+nlin+nonlin+1, 1); > betaq = d1_malloc(n+nlin+nonlin+1, 1); > diag = d1_malloc(n+nlin+nonlin+1, 1); > cscal = d1_malloc(n+nlin+nonlin+1, 1); > colle = d1_malloc(n+nlin+nonlin+1, 1); > > a = d2_malloc(n+1, n+1, 1); > diag0 = d1_malloc(n+1, 1); > violis = i1_malloc(2*nstep*(n+nlin+nonlin)+1, 1); > alist = i1_malloc(n+nlin+nonlin+1, 1); > o8bind = i1_malloc(2*(n+nlin+nonlin)+1, 1); > o8bind0 = i1_malloc(2*(n+nlin+nonlin)+1, 1); > aalist = i1_malloc(n+nlin+nonlin+1, 1); > clist = i1_malloc(nlin+nonlin+1, 1); > u = d1_malloc(2*(n+nlin+nonlin)+1, 1); > u0 = d1_malloc(2*(n+nlin+nonlin)+1, 1); > w = d1_malloc(2*(n+nlin+nonlin)+1, 1); > w1 = d1_malloc(2*(n+nlin+nonlin)+1, 1); > res = d1_malloc(2*(n+nlin+nonlin)+1, 1); > res0 = d1_malloc(2*(n+nlin+nonlin)+1, 1); > res1 = d1_malloc(2*(n+nlin+nonlin)+1, 1); > resst = d1_malloc(2*(n+nlin+nonlin)+1, 1); > yu = d1_malloc(n+nlin+nonlin+1, 1); > slack = d1_malloc(2*(n+nlin+nonlin)+1, 1); > work = d1_malloc(2*(n+nlin+nonlin)+1, 1); > ug = d1_malloc(n+1, 1); > og = d1_malloc(n+1, 1); > low = d1_malloc(n+nlin+nonlin+1, 1); > up = d1_malloc(n+nlin+nonlin+1, 1); > xst = d1_malloc(n+1, 1); > > > /* o8fint.h */ > xtr = d1_malloc(n+1, 1); > xsc = d1_malloc(n+1, 1); > fu = d1_malloc(nlin+nonlin+1, 1); > fugrad = d2_malloc(n+1, nlin+nonlin+1, 1); > fud = d2_malloc(nlin+nonlin+1, 7, 1); > > /* o8fuco.h */ > val = l1_malloc(nlin+nonlin+1, 1); > llow = l1_malloc(n+1, 1); > lup = l1_malloc(n+1, 1); > cres = i1_malloc(nlin+nonlin+1, 1); > cgres = i1_malloc(nlin+nonlin+1, 1); > confuerr = l1_malloc(nlin+nonlin+1, 1); > nonlinlist = i1_malloc(nlin+nonlin+1, 1); > gunit = i2_malloc(4, nlin+nonlin+1, 1); > gconst = l1_malloc(nlin+nonlin+1, 1); > cfuerr = l1_malloc(nlin+nonlin+1, 1); > > /* o8gene.h */ > np = d1_malloc(ndualm+1, 1); > xj = d2_malloc(ndualm+1, ndualm+1, 1); > ddual = d1_malloc(ndualm+1, 1); > r = d2_malloc(ndualm+1, ndualm+1, 1); > ud = d1_malloc(mdualm+1, 1); > ud1 = d1_malloc(mdualm+1, 1); > aitr = i1_malloc(mdualm+1, 1); > > > /* GLOBAL VARIABLES USED LOCALLY */ > > donlp2_yy = d1_malloc(n+1, 1); > o8bfgs_dg = d1_malloc(n+1, 1); > o8bfgs_adx = d1_malloc(n+1, 1); > o8bfgs_ltdx = d1_malloc(n+1, 1); >/* printf("d_borders_count: %d\n", d_borders_count); */ > o8bfgs_gtdx = d1_malloc(n+nlin+nonlin+1, 1); > o8bfgs_updx = d1_malloc(n+1, 1); > o8bfgs_updz = d1_malloc(n+1, 1); > o8opti_qtx = d1_malloc(n+1, 1); > o8opti_yy = d1_malloc(n+1, 1); > o8opti_yx = d1_malloc(n+1, 1); > o8opti_trvec = d1_malloc(n+1, 1); > o8opti_con = d1_malloc(nlin+nonlin+1, 1); > o8opti_delist = i1_malloc(n+nlin+nonlin+1, 1); > o8opti_bindba = i1_malloc(2*(n+nlin+nonlin)+1, 1); > o8eval_con = d1_malloc(nlin+nonlin+1, 1); > o8unim_step = d1_malloc(nstep+1, 1); > o8dec_qri = d1_malloc(n+1, 1); > o8dec_qri0 = d1_malloc(n+1, 1); > o8elim_qri = d1_malloc(n+1, 1); > o8elim_y = d1_malloc(n+1, 1); > o8elim_rhsscal= d1_malloc(n+nlin+1, 1); > o8elim_col = i1_malloc(n+nlin+1, 1); > o8sol_xl = d1_malloc(n+1, 1); > o8upd_sdiag = d1_malloc(n+1, 1); > o8upd_rn1 = d1_malloc(n+1, 1); > o8upd_w = d1_malloc(n+1, 1); > o8qpdu_g0 = d1_malloc(ndualm+1, 1); > o8qpdu_ci0 = d1_malloc(mdualm+1, 1); > o8qpdu_cii = d1_malloc(mdualm+1, 1); > o8qpdu_cei = d1_malloc(ndualm+1, 1); > o8qpdu_xd = d1_malloc(ndualm+1, 1); > o8qpdu_s = d1_malloc(mdualm+1, 1); > o8qpdu_z = d1_malloc(ndualm+1, 1); > o8qpdu_vr = d1_malloc(mdualm+1, 1); > o8qpdu_xdold = d1_malloc(ndualm+1, 1); > o8qpdu_udold = d1_malloc(mdualm+1, 1); > o8qpdu_ai = i1_malloc(mdualm+1, 1); > o8qpdu_iai = i1_malloc(mdualm+1, 1); > o8qpdu_iaexcl = i1_malloc(mdualm+1, 1); > o8qpdu_aiold = i1_malloc(mdualm+1, 1); > o8qpdu_eqlist = i1_malloc(n+1, 1); > o8qpdu_iqlist = i1_malloc(mdualm+1, 1); > o8qpdu_y = d1_malloc(ndualm+1, 1); > o8qpdu_mult = d1_malloc(n+nlin+nonlin+1, 1); > o8qpdu_qpdel = i1_malloc(mdualm+1, 1); > escongrad_errloc = i1_malloc(nlin+nonlin+1, 1); > escongrad_fhelp1 = d1_malloc(nlin+nonlin+1, 1); > escongrad_fhelp2 = d1_malloc(nlin+nonlin+1, 1); > escongrad_fhelp3 = d1_malloc(nlin+nonlin+1, 1); > escongrad_fhelp4 = d1_malloc(nlin+nonlin+1, 1); > escongrad_fhelp5 = d1_malloc(nlin+nonlin+1, 1); > escongrad_fhelp6 = d1_malloc(nlin+nonlin+1, 1); > >} > >/**********************************************************************/ >/*** free memory from global arrays *******/ >/**********************************************************************/ >void global_mem_free() { > > /* GLOBAL VARIABLES USED GLOBALLY */ > > /* o8comm.h */ > d2_free(accinf, iterma+1); > d1_free(donlp2_x); > d1_free(x0); > d1_free(x1); > d1_free(xmin); > d1_free(d); > d1_free(d0); > d1_free(dd); > d1_free(difx); > d1_free(resmin); > d1_free(gradf); > d1_free(qgf); > d1_free(gphi0); > d1_free(gphi1); > d2_free(gres, n+1); > d1_free(gresn); > i1_free(perm); > i1_free(perm1); > i1_free(colno); > d2_free(qr, n+1); > d1_free(betaq); > d1_free(diag); > d1_free(cscal); > d1_free(colle); > > d2_free(a, n+1); > d1_free(diag0); > i1_free(violis); > i1_free(alist); > i1_free(o8bind); > i1_free(o8bind0); > i1_free(aalist); > i1_free(clist); > d1_free(u); > d1_free(u0); > d1_free(w); > d1_free(w1); > d1_free(res); > d1_free(res0); > d1_free(res1); > d1_free(resst); > d1_free(yu); > d1_free(slack); > d1_free(work); > d1_free(ug); > d1_free(og); > d1_free(low); > d1_free(up); > d1_free(xst); > > > /* o8fint.h */ > d1_free(xtr); > d1_free(xsc); > d1_free(fu); > d2_free(fugrad, n+1); > d2_free(fud, nlin+nonlin+1); > > /* o8fuco.h */ > l1_free(val); > l1_free(llow); > l1_free(lup); > i1_free(cres); > i1_free(cgres); > l1_free(confuerr); > i1_free(nonlinlist); > i2_free(gunit, 4); > l1_free(gconst); > l1_free(cfuerr); > > /* o8gene.h */ > d1_free(np); > d2_free(xj, ndualm+1); > d1_free(ddual); > d2_free(r, ndualm+1); > d1_free(ud); > d1_free(ud1); > i1_free(aitr); > > /* GLOBAL VARIABLES USED LOCALLY */ > d1_free(donlp2_yy); > d1_free(o8bfgs_dg); > d1_free(o8bfgs_adx); > d1_free(o8bfgs_ltdx); > d1_free(o8bfgs_gtdx); > d1_free(o8bfgs_updx); > d1_free(o8bfgs_updz); > d1_free(o8opti_qtx); > d1_free(o8opti_yy); > d1_free(o8opti_yx); > d1_free(o8opti_trvec); > d1_free(o8opti_con); > i1_free(o8opti_delist); > i1_free(o8opti_bindba); > d1_free(o8eval_con); > d1_free(o8unim_step); > d1_free(o8dec_qri); > d1_free(o8dec_qri0); > d1_free(o8elim_qri); > d1_free(o8elim_y); > d1_free(o8elim_rhsscal); > i1_free(o8elim_col); > d1_free(o8sol_xl); > d1_free(o8upd_sdiag); > d1_free(o8upd_rn1); > d1_free(o8upd_w); > d1_free(o8qpdu_g0); > d1_free(o8qpdu_ci0); > d1_free(o8qpdu_cii); > d1_free(o8qpdu_cei); > d1_free(o8qpdu_xd); > d1_free(o8qpdu_s); > d1_free(o8qpdu_z); > d1_free(o8qpdu_vr); > d1_free(o8qpdu_xdold); > d1_free(o8qpdu_udold); > i1_free(o8qpdu_ai); > i1_free(o8qpdu_iai); > i1_free(o8qpdu_iaexcl); > i1_free(o8qpdu_aiold); > i1_free(o8qpdu_eqlist); > i1_free(o8qpdu_iqlist); > d1_free(o8qpdu_y); > d1_free(o8qpdu_mult); > i1_free(o8qpdu_qpdel); > i1_free(escongrad_errloc); > d1_free(escongrad_fhelp1); > d1_free(escongrad_fhelp2); > d1_free(escongrad_fhelp3); > d1_free(escongrad_fhelp4); > d1_free(escongrad_fhelp5); > d1_free(escongrad_fhelp6); >} > >
<b>You cannot view the attachment while viewing its details because your browser does not support IFRAMEs. <a href="attachment.cgi?id=267611">View the attachment on a separate page</a>.</b>
View Attachment As Raw
Actions:
View
Attachments on
bug 396921
: 267611 |
267621