Rdonlp2 - an R interface to DONLP2

Ryuichi Tamura(ry.tamura @ gmail.com)


Date: June 7, 2007(Version 0.3-1)


Contents

1 Abount This Package

Rdonlp2 is an R package to use Peter Spellucci's DONLP2 from R. DONLP2 is the copyrighted software written by Peter Spellucci for solving nonlinear programming problems. DONLP2 is available from: Rdonlp2 is a wrapper for ANSI-C version of DONLP2(also called DONLP3): Current Rdonlp2 package is available from: Since Rdonlp2 is simply wrapper program, user is required to refer PDF manual included in donlp2_intv_dyn.tar.gz for the detail of the algorithm. As for condition of use, Please refer Section8.

2 Problem Definition

R Package Rdonlp2 solves following nonlinear minimization problem with linear, nonlinear constraints as well as parameter bounds:

$\displaystyle \min_x f(x)$    subject to $\displaystyle \quad x\in \mathcal{S}$    
$\displaystyle \mathcal{S} \in$ $\displaystyle \quad\{ \quad x\in R^n,$    
  $\displaystyle \qquad x_l\le x\le x_u,$    
  $\displaystyle \qquad b_l \le Ax \le b_u,$    
  $\displaystyle \qquad c_l \le c(x) \le c_u \quad\},$    

where, $ f(x)$ is a continuous function, $ x_l,x_u$ are bounds for parameters($ x$), $ b_l, b_u$ are bounds for linear combinations $ Ax$(linear constraints), and $ c_l,c_u$ are bounds for nonlinear function $ c(x)$(nonlinear constraints). To describe equality constraint or parameter constancy, let lower and upper bounds for constraint be equal.

2.1 R function donlp2()

Rdonlp provides single function for this problem: 
 donlp2 <- function(par, fn,
                   par.upper=rep(+Inf, length(par)),
                   par.lower=rep(-Inf, length(par)),

                   A = NULL,
                   lin.upper=rep(+Inf, length(par)),
                   lin.lower=rep(-Inf, length(par)),

                   nlin = list(),
                   nlin.upper=rep(+Inf, length(nlin)),
                   nlin.lower=rep(-Inf, length(nlin)),

                   control=donlp2.control(),
                   control.fun=function(lst){return(TRUE)},
                   env=.GlobalEnv, name="Rdonlp2")
where,
fn - the objective function to be minimized. Currently, fn must take only one argument, and the parameter vector(par) will be passed to fn during the optimization. The first element of return value must be the evaluated value.
par - parameter vector(vector object)
par.upper, par.lower - upper and lower bounds for parameter vector, respectively. Their length must equal to length(par). If some elements are unbounded, specify +Inf or -Inf explicitly.
A - the matrix object that represents linear constraints. Its columns must be equal to length(par), and its rows must be equal to the number of linear constraints.
lin.upper, lin.lower - upper and lower bounds for linear constraints, respectively. Their length must equal to the number of linear constraints. If some elements are unbounded, specify +Inf or -Inf explicitly.
nlin - list object whose elements are functions that represents nonlinear constraints. rule for argument and return value is the same as fn, i.e., these functions take only one arugument(par), and return a vector object whose first element is the evaluated value.
lin.upper, lin.lower - upper and lower bounds for nonlinear constraints, respectively. Their length must equal to length(nlin). If some elements are unbounded, specify +Inf or -Inf explicitly.
control - control parameters that defines the behavior of DONLP2. See below for details.
control.fun - donlp2() reports a group of optimization parameters in every iteration(see below for details). This (read-only) information can be available within control.fun(), in which user can decide whether the optimization should be iterrupted. Set its return value to FALSE to interrupt donlp2().
env - the environment in which objective, constraint, control functions are evaluated.
name - an character object that specify file name(without extension) of output file. DONLP2 can output following 2 files(name.pro and name.mes) in working directory which contain detailed information during the optimization.

3 Details

3.1 Initial Values

Although intial values(given to par) should be carefully determined by user, DONLP2 has the feature to correct initial values automatically that violate given constraints.


3.2 Objective Function and its Gradients

Objective function should be implemented so that parameter vector is the only argument and the first element of return value is numeric. Typically, it looks like: 
objective.fun <- function(par){
  # calculation using par, and results are stored to ans
  ans # return value
}

:

ret <- donlp2(par=par, fn=objective.fun, ....)
User can implement gradient function to improve accuracy and efficiency. Gradient function should be implement so that parameter vector is the only argument and return value is vector of length equal to length(par). Then, assign it to the attibute gr of objective function: 
# par is vector of length n
grad.fun <- function(par){
  
  c(v1, v2, ..., vn)
}
# assign grad.fun to objective.fun's gr attribute
attr(objective.fn, "gr") <- grad.fun

3.3 Bounds and Equality Constraints

Bounds for parameter, linear and nonlinear constraints are given as vector of appropriate length(par.upper,par.lower,lin.upper,lin.lower,nlin.upper,nlin.lower). If some parameter or constraints are bounded from below (above), then specify +Inf(-Inf), respectively. Set upper and lower bounds to be equal if equality constraint is needed. 
# par[1]<0, 0<par[2]<1, par[3]>1
par.lower <- c(-Inf, 0,    1)
par.upper <- c(   0, 1, +Inf)

# two linear constraints on two parameters
# (1) par[1]+par[2]=0
# (2) par[1]-2*par[2]+10>0
lin.lower <- c(0,  -10)
lin.upper <- c(0, +Inf)

3.4 Linear Constraints

Bounds arguments and the coefficient matrix A represents the linear constraints. Every row of A stands for linear combination of parameters: 
# two linear constraints on two parameters
# (1) par[1]+par[2]=0
# (2) par[1]-2*par[2]+10>0
lin.lower <- c(0,  -10)
lin.upper <- c(0, +Inf)

A = rbind( c(1, 1),   # 1st linear compination
           c(1,-2) )  # 2nd linear combination

3.5 Nonlinear Constraints and its Gradients

Nonlinear constraints are represented as their bounds given to nlin.upper and nlin.lower, and user-defined function(and gradients). The way to implement function and gradients are the same as objective function and its gradients(see Section 3.2): 
# Nonlinear constraints: 1 constraint on 2 parameters
# par[1]*par[2] = 1
nlcon1 <- function(par){
  par[1]*par[2]
}
nlcon1.gr <- function(par){
  c(par[2], par[1])
}
attr(nlcon1, "gr")<-nlcon1.gr

nlin.upper = c(1)
nlin.lower = c(1)
:
ret <- donlp2(par, fn, 
              nlin=list(nlcon1), 
              nlin.upper=nlin.upper,nlin.lower=nlin.lower,.....)
All the nonlinear constraint function are collected in a list(nlin).


3.6 Numerical Gradients

User can omit the implementation of gradients. In this case one of 3 algorithm fro numerical differentiation provided by DONLP2 will be performed. If there are n parameters,
  1. the forward difference: requires n additional evaluations of function(difftype=1).
  2. the central difference: requires 2n additional evaluations of function(difftype=2).
  3. a sixth order approximation computing a Richardson extrapolation of the three symmetric difference: requires 6n additional evaluations of function(difftype=3).
By default, Rdonlp2 uses 3rd algorithm(most acculate, but quite costly). User can change this by the control variable(below) difftype. Currently, if user want to use analytical gradients, he must implement all of the gradients for both objective function and nonlinear constraint functions. If one of them are not implemented, Rdonlp2 gives up using any gradient function and uses numerical method instead.

3.7 Control Variables

User can control the behavior of donlp2() by donlp2.control(). donlp2.control()returns a group of default control parameters as list object, so user change some of them by giving tag=value pairs as arguments. Currently following tags(control variables) are available (values in () are the defaults). Some of parameteters listed above are not well documented in this tutorial. Please refer the original PDF manual included in DONLP distribution for details.


4 Information during the Optimization

User may want to know what happens during the optimation, and if some parameters are 'undesirable' he may stop the execution. Rdonlp2 provides the way to access the information via control.fun(lst). On each iterantion, 35 parameters that tell us how the optimization is working are passed to control.fun(lst). The last four parameters are not reported until the optimization has finished. Parameters are collected into single list object, so user can easily access some of them by specifying the tags. Complete tag list is given Section 4.1. control.fun() should return TRUE if user want to continue and FALSE if user want to interrupt. 
## Example 1
## keep track of current lagrange multiplier values
mycontrol <- function(lst){
 print(lst$u)
 TRUE # return TRUE to continue execution
}

# tell donlp2 to use mycontrol()
donlp2(.....,control.fun=mycontrol,....)

## Example 2(useless example)
## force to terminate optimization after 10 iterations
mycontrol2 <- function(lst){
  lst$step.nr <= 10 # return FALSE when step.nr>11
}
# tell donlp2 to use mycontrol2()
donlp2(.....,control.fun=mycontrol2,....)


4.1 Tag list

Some of tags listed here are not well documented in this tutorial. Please refer the original PDF manual included in DONLP distribution for details.

5 Value from donlp2()

The return value from donlp2() is the list object with 38(35+3; if hessian=FALSE(default)) or 39(35+4; if hessian=TRUE) elements with specified tags. The 35 pareameters are identical to those listed in Section 4.1, and the rest parameters are:


5.1 The termination Reason

When the optimization finishes, DONLP2 returns one of 19 messages listed below. They are classified to following 3 groups, last of which user need to decide the result is 'reasonable' and accestable.
  1. "constraint evaluation returns error with current point"
  2. "objective evaluation returns error with current point"
  3. "qpsolver: working set singular in dual extended qp "
  4. "qpsolver: extended qp-problem seemingly infeasible "
  5. "qpsolver: no descent for infeas from qp for tau=tau_max"
  6. "qpsolver: on exit correction small, infeasible point"
  7. "stepsizeselection: computed d from qp not a dir. of descent"
  8. "more than maxit iteration steps"
  9. "stepsizeselection: no acceptable stepsize in [sigsm,sigla]"
  10. "small correction from qp, infeasible point"
  11. "kt-conditions satisfied, no further correction computed"
  12. "computed correction small, regular case "
  13. "stepsizeselection: x almost feasible, dir. deriv. very small"
  14. "kt-conditions (relaxed) satisfied, singular point"
  15. "very slow primal progress, singular or illconditoned problem"
  16. "more than nreset small corrections in x "
  17. "correction from qp very small, almost feasible, singular "
  18. "numsm small differences in penalty function,terminate"
  19. "user required termination"
Some of messages listed above are not well documented in this tutorial. Please refer the original PDF manual included in DONLP distribution for details.

6 Examples

Example C source files are included in the original DONLP2 distribution. In this section, we rewrite some of them in R and show you how to code constrained optimization problem with Rdonlp2.

6.1 examples/simple.c: linear constraint

$\displaystyle \min_{x,y}x^2+y^2$   subject to$\displaystyle \quad 0< x,y< 100,\quad x+y=1$ (1)

with intial value: $ (x,y)=(-10,10)$. Note that initial value is infeasible because $ x=-10\notin (0,100)$. This problem has 2 parameters, 2 parameter bounds, 1 linear equality constraint. R script looks like: 
p <- c(-10,10)
par.l <- c(0,0); par.u <- c(100,100)

lin.u <- 1; lin.l <- 1
A <- t(c(1,1))

fn <- function(x){
  x[1]^2+x[2]^2
}
ret <- donlp2(p, fn, par.lower=par.l, par.upper=par.u,
              A=A, lin.u=lin.u, lin.l=lin.l, name="simple")
Note that A must be represented as row vector(1x2) with single linear constraint. Also we use numerical gradients. Since control variables are all default values(specifically te0=TRUE and silent=FALSE), running the script outputs detailed information on console: 
    1 fx=   0.000000e+00 upsi=  5.9e+01 b2n= -1.0e+00 umi=  0.0e+00 nr   1 si-1
    2 fx=   0.000000e+00 upsi=  2.0e+01 b2n= -1.0e+00 umi=  0.0e+00 nr   2 si-1
    3 fx=   1.000000e+00 upsi=  0.0e+00 b2n=  4.4e-16 umi=  0.0e+00 nr   2 si-1
simple

     n=         2    nlin=         1    nonlin=         0

  epsx= 1.000e-05 sigsm= 3.293e-10

startvalue
  5.0000000e+01   1.0000000e+01 

  eps=  1.08e-19  tol=  0.00e+00 del0=  1.00e+00 delm=  1.00e-06 tau0=  1.00e+00
  tau=  1.00e-01   sd=  1.00e-01   sw=  2.27e-13  rho=  1.00e-06 rho1=  1.00e-10
 scfm=  1.00e+04  c1d=  1.00e-02 epdi=  1.00e-16
  nre=         4 anal=         0
taubnd=  1.00e+00 epsfcn=  1.00e-16 difftype=3

 termination reason:
 KT-conditions satisfied, no further correction computed
 evaluations of f                            3
 evaluations of grad f                       2
 evaluations of constraints                  6
 evaluations of grads of constraints         0
 final scaling of objective           1.000000e+00
 norm of grad(f)                      1.414214e+00
 lagrangian violation                 4.718134e-14
 feasibility violation                0.000000e+00
 dual feasibility violation           0.000000e+00
 optimizer runtime sec's              0.000000e+00


 optimal value of f =   5.00000000000000e-01

 optimal solution  x =
  4.99999999999991e-01  5.00000000000009e-01

  multipliers are relativ to scf=1
  nr.    constraint     multiplier norm(grad) or 1 
    1   5.0000000e-01    0.0000000e+00 
    2   9.9500000e+01    0.0000000e+00 
    3   5.0000000e-01    0.0000000e+00 
    4   9.9500000e+01    0.0000000e+00 
    5   0.0000000e+00    1.0000000e+00   1.3906925e-309
    6   0.0000000e+00    0.0000000e+00 

 evaluations of restrictions and their gradients
 (     6,     0)
 last estimate of cond.nr. of active gradients   1.414e+00

 last estimate of cond.nr. of approx.  hessian   1.000e+00
iterative steps total               3
# of restarts                       0
# of full regular updates           1
# of updates                        1
# of regularized full SQP-steps     0
If user want to access parameter values, simply 
> ret$par
[1] 0.5 0.5

6.2 examples/simple2.c: nonlinear constraint

Second example shows optimization problem with parameter bounds and nonlinear constraint.

$\displaystyle \min_{x,y}x^2+y^2$   subject to$\displaystyle \quad -100< x,y< 100,\quad xy=2$ (2)

with intial value: $ (x,y)=(10,10)$. We give gradient functions for objective and nonlinear constraint(dfn() and dnlcon(), respectively). 
p <- c(10,10)
par.l <- c(-100,-100); par.u <- c(100,100)
nlin.l <- nlin.u <- 2

fn <- function(x){
  x[1]^2+x[2]^2
}
dfn <- function(x){
  c(2*x[1], 2*x[2])
}
attr(fn, "gr") <- dfn

nlcon <- function(x){
  x[1]*x[2]
}
dnlcon <- function(x){
  c(x[2], x[1])
}
attr(nlcon, "gr") <- dnlcon

ret<-donlp2(p, fn, par.u=par.u, par.l=par.l,
            nlin=list(nlcon), nlin.u=nlin.u, nlin.l=nlin.l)

6.3 example/hs211.c

This problem uses all types of constraints:

$\displaystyle \min_{x_i, i=1\ldots 10}$ $\displaystyle 5.04x_1+0.035x_2+10x_3+3.36x_5-0.063x_4x_7$    
subject to:    
$\displaystyle h_1(x)$ $\displaystyle = 1.22x_4 - x_1 - x_5 = 0$    
$\displaystyle h_2(x)$ $\displaystyle = 98000x_3/(x_4x_9+1000x_3)-x_6=0$    
$\displaystyle h_3(x)$ $\displaystyle = (x_2+x_5)/x_1-x_8 = 0$    
$\displaystyle g_1(x)$ $\displaystyle = 35.82-0.222x_{10}-bx_9\ge 0, b=0.9$    
$\displaystyle g_2(x)$ $\displaystyle = -133+3x_7 -ax_10 \ge 0, a=0.99$    
$\displaystyle g_3(x)$ $\displaystyle = -g_1(x) + x_9(1/b-b)\ge 0$    
$\displaystyle g_4(x)$ $\displaystyle = -g_2(x) +(1/a-a)x_{10} \ge 0$    
$\displaystyle g_5(x)$ $\displaystyle = 1.12x_1+0.13167x_1x_8 - 0.00667x_1x_8^2-ax_4 \ge 0$    
$\displaystyle g_6(x)$ $\displaystyle = 57.425 + 1.098x_8 - 0.038x_8^2 + 0.325x_6 - ax_7 \ge 0$    
$\displaystyle g_7(x)$ $\displaystyle = -g_5(x) + (1/a-a)x_4 \ge 0$    
$\displaystyle g_8(x)$ $\displaystyle = -g_6(x) + (1/a-a)x_7 \ge 0$    
$\displaystyle 0.00001$ $\displaystyle \le x_1 \le 2000$    
$\displaystyle 0.00001$ $\displaystyle \le x_2 \le 16000$    
$\displaystyle 0.00001$ $\displaystyle \le x_3 \le 120$    
$\displaystyle 0.00001$ $\displaystyle \le x_4 \le 5000$    
$\displaystyle 0.00001$ $\displaystyle \le x_5 \le 2000$    
$\displaystyle 85$ $\displaystyle \le x_6 \le 93$    
$\displaystyle 90$ $\displaystyle \le x_7 \le 95$    
$\displaystyle 3$ $\displaystyle \le x_8 \le 12$    
$\displaystyle 1.2$ $\displaystyle \le x_9 \le 4$    
$\displaystyle 145$ $\displaystyle \le x_{10} \le 162$    

We have 10 parameter bounds, 5 linear constraints(arranging terms):

$\displaystyle h_1$ $\displaystyle \rightarrow 1.22x_4 - x_1 - x_5 = 0$    
$\displaystyle g_1$ $\displaystyle \rightarrow -0.222x_10-bx_9\ge -35.82, b=0.9$    
$\displaystyle g_2$ $\displaystyle \rightarrow 3x_7 -ax_10 \ge 133, a=0.99$    
$\displaystyle g_3$ $\displaystyle \rightarrow x_9(1/b-b+b)+0.222x_{10}\ge 35.82$    
$\displaystyle g_4$ $\displaystyle \rightarrow -3x_7+(1/a-a+a)x_{10} \ge -133,$    

1 of which($ h_1$) is equality constraint, and 6 nonlinear constraints:

$\displaystyle h_2$ $\displaystyle \rightarrow 98000x_3/(x_4x_9+1000x_3)-x_6=0$    
$\displaystyle h_3$ $\displaystyle \rightarrow (x_2+x_5)/x_1-x_8 = 0$    
$\displaystyle g_5$ $\displaystyle \rightarrow 1.12x_1+0.13167x_1x_8 - 0.00667x_1x_8^2-ax_4 \ge 0$    
$\displaystyle g_6$ $\displaystyle \rightarrow 57.425 + 1.098x_8 - 0.038x_8^2 + 0.325x_6 - ax_7 \ge 0$    
$\displaystyle g_7$ $\displaystyle \rightarrow -g_5(x) + (1/a-a)x_4 \ge 0$    
$\displaystyle g_8$ $\displaystyle \rightarrow -g_6(x) + (1/a-a)x_7 \ge 0,$    

2 of which($ h_2$, $ h_3$) are equality constraint. R code will be follows: 
a <- 0.99; b <- 0.9
##
## Objective Function
##
fn <- function(par){
  x1 <- par[1]; x2 <- par[2]; x3 <- par[3]; x4 <- par[4]; x5 <- par[5]
  x6 <- par[6]; x7 <- par[7]; x8 <- par[8]; x9 <- par[9]; x10 <- par[10]

  5.04*x1 + 0.035*x2 + 10*x3 +3.36*x5 - 0.063*x4*x7
}
##
## Parameter Bounds
##
par.l <- c(rep(1e-5, 5), 85, 90, 3, 1.2, 145)
par.u <- c(2000, 16000, 120, 5000, 2000, 93, 95, 12, 4, 162)

##
## Constraints
##
linbd  <- matrix(0, nr=5, nc=2)
nlinbd <- matrix(0, nr=6, nc=2)

## linear equality 
linbd[1,] <- c(0,0)          # h1
linbd[2,] <- c(-35.82, Inf)  # g1
linbd[3,] <- c(133, Inf)     # g2
linbd[4,] <- c(35.82,Inf)    # g3
linbd[5,] <- c(-133, Inf)    # g4

## nonlinear equality
h2 <- function(par){
  x1 <- par[1]; x2 <- par[2]; x3 <- par[3]; x4 <- par[4]; x5 <- par[5]
  x6 <- par[6]; x7 <- par[7]; x8 <- par[8]; x9 <- par[9]; x10 <- par[10]

  98000*x3/(x4*x9+1000*x3)-x6
}
nlinbd[1,] <- c(0,0)

## nonlinear equality
h3 <- function(par){
  x1 <- par[1]; x2 <- par[2]; x3 <- par[3]; x4 <- par[4]; x5 <- par[5]
  x6 <- par[6]; x7 <- par[7]; x8 <- par[8]; x9 <- par[9]; x10 <- par[10]

  (x2+x5)/x1 - x8
}
nlinbd[2,] <- c(0,0)

## nonlinear inequality
g5 <- function(par){
  x1 <- par[1]; x2 <- par[2]; x3 <- par[3]; x4 <- par[4]; x5 <- par[5]
  x6 <- par[6]; x7 <- par[7]; x8 <- par[8]; x9 <- par[9]; x10 <- par[10]

  1.12*x1 + 0.13167*x1*x8 - 0.00667*x1*x8^2 - a*x4
}
nlinbd[3,] <- c(0,Inf)

## nonlinear inequality
g6 <- function(par){
  x1 <- par[1]; x2 <- par[2]; x3 <- par[3]; x4 <- par[4]; x5 <- par[5]
  x6 <- par[6]; x7 <- par[7]; x8 <- par[8]; x9 <- par[9]; x10 <- par[10]

  57.425 + 1.098*x8 - 0.038*x8^2 + 0.325*x6 - a*x7
}
nlinbd[4,] <- c(0,Inf)

## nonlinear inequality
g7 <- function(par){
  x1 <- par[1]; x2 <- par[2]; x3 <- par[3]; x4 <- par[4]; x5 <- par[5]
  x6 <- par[6]; x7 <- par[7]; x8 <- par[8]; x9 <- par[9]; x10 <- par[10]

  -g5(par) + (1/a-a)*x4
}
nlinbd[5,] <- c(0,Inf)

## nonlinear inequality
g8 <- function(par){
  x1 <- par[1]; x2 <- par[2]; x3 <- par[3]; x4 <- par[4]; x5 <- par[5]
  x6 <- par[6]; x7 <- par[7]; x8 <- par[8]; x9 <- par[9]; x10 <- par[10]

  -g6(par) + (1/a-a)*x7
}
nlinbd[6,] <- c(0,Inf)

## A is 5(linear constraints) x 10(params) matrix
A <- rbind(c(-1, 0, 0, 1.22,-1, 0,  0, 0,        0,         0), #h1
           c( 0, 0, 0,    0, 0, 0,  0, 0,       -b,    -0.222), #g1
           c( 0, 0, 0,    0, 0, 0,  3, 0,        0,        -a), #g2
           c( 0, 0, 0,    0, 0, 0,  0, 0,(1/b-b)+b,     0.222), #g3
           c( 0, 0, 0,    0, 0, 0, -3, 0,        0,(1/a-a)+a))  #g4

## initial values
p0 <- c(1745, 12e3, 11e1, 3048, 1974, 89.2, 92.8, 8, 3.6, 145)

## control variables
cntl <- donlp2.control(del0=0.2, tau0=1.0, tau=0.1, taubnd=5e-6)

## start constrained optimization
ret <- donlp2(par=p0, fn=fn,
              par.u=par.u, par.l=par.l,
              A=A,
              lin.u=linbd[,2], lin.l=linbd[,1],
              nlin=list(h2,h3,g5,g6,g7,g8),
              nlin.upper=nlinbd[,2], nlin.lower=nlinbd[,1], name="alkylation",
              control=cntl)
Since gradient functions are not implemented in the program and the numerical differentiation algorithm is difftype=3(default), it takes about 0.90 sec to finish the computation on my machine(Intel CoreDuo 2G, Memory 2G), while original C version does within 0.1 sec.

7 Bugs

DONLP2 provides much more 'fine-tuning' parameters than those exported to Rdonlp2. So if you want other parameters that should be exported, please e-mail me. Also, any comments of suggestions are highly welcome.


8 Copyright

Original DONLP2: 

/* 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.    */

Rdonlp2

Copyright (C) 2007 Ryuichi Tamura(ry.tamura at gmail.com). You may re-distribute or modify this library under GNU LGPL ver.2.
TAMURA Ryuichi
2007-06-07