/*
Copyright (C) 2010 Richard Gomes
This file is part of JQuantLib, a free-software/open-source library
for financial quantitative analysts and developers - http://jquantlib.org/
JQuantLib is free software: you can redistribute it and/or modify it
under the terms of the QuantLib license. You should have received a
copy of the license along with this program; if not, please email
. The license is also available online at
.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the license for more details.
JQuantLib is based on QuantLib. http://quantlib.org/
When applicable, the original copyright notice follows this notice.
*/
/*
The original Fortran version is Copyright (C) 1999 University of Chicago.
All rights reserved.
Redistribution and use in source and binary forms, with or
without modification, are permitted provided that the
following conditions are met:
1. Redistributions of source code must retain the above
copyright notice, this list of conditions and the following
disclaimer.
2. Redistributions in binary form must reproduce the above
copyright notice, this list of conditions and the following
disclaimer in the documentation and/or other materials
provided with the distribution.
3. The end-user documentation included with the
redistribution, if any, must include the following
acknowledgment:
This product includes software developed by the
University of Chicago, as Operator of Argonne National
Laboratory.
Alternately, this acknowledgment may appear in the software
itself, if and wherever such third-party acknowledgments
normally appear.
4. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
BE CORRECTED.
5. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
(INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
POSSIBILITY OF SUCH LOSS OR DAMAGES.
C translation Copyright (C) Steve Moshier
What you see here may be used freely but it comes with no support
or guarantee.
*/
package org.jquantlib.math.optimization;
import org.jquantlib.math.matrixutilities.Array;
import org.jquantlib.math.matrixutilities.Matrix;
import org.jquantlib.math.matrixutilities.internal.Address;
/**
* This class provides linear algebra optimization algorithms.
*
* The original Fortran version is Copyright (C) 1999 University of Chicago. All rights reserved.
* C translation :: Copyright (C) Steve Moshier
* Adapted to JQuantLib :: Copyright (C) 2010 Richard Gomes
*
* Brief history:
* Original sources were developed in FORTRAN language.
* Later sources were translated to C language.
* Finally we made this very crude (and certainly insufficient) adaptation to work with JQuantLib/Java.
*
* Ideally, these sources should be ready to smoothly work with classes {@link Matrix} and {@link Array}.
* At the moment, we require that underlying data structures are accessed directly, which is clearly very
* far from ideal because it breaks the paradigm offered by interface {@link Address}.
*
* @author Richard Gomes
*/
public class Minpack {
public static void qrfac(
final int m,
final int n,
final Matrix a,
final boolean pivot,
final int[] ipvt,
final Array rdiag,
final Array acnorm,
final Array wa) {
MinpackC.qrfac(m, n, a.$, pivot?1:0, ipvt, rdiag.$, acnorm.$, wa.$);
}
private static class MinpackC {
private static final double MACHEP = 1.2e-16; // resolution of arithmetic
private static final double DWARF = 1.0e-38; // smallest nonzero number
private static final double rdwarf = 3.834e-20;
private static final double rgiant = 1.304e19;
private static final double zero = 0.0;
private static final double one = 1.0;
private static final double p5 = 0.5;
private static final double p05 = 0.05;
private static final double p25 = 0.25;
private static double enorm(final int n, final double[] x, final int offset)
{
/*
* **********
*
* function enorm
*
* given an n-vector x, this function calculates the
* euclidean norm of x.
*
* the euclidean norm is computed by accumulating the sum of
* squares in three different sums. the sums of squares for the
* small and large components are scaled so that no overflows
* occur. non-destructive underflows are permitted. underflows
* and overflows do not occur in the computation of the unscaled
* sum of squares for the intermediate components.
* the definitions of small, intermediate and large components
* depend on two constants, rdwarf and rgiant. the main
* restrictions on these constants are that rdwarf**2 not
* underflow and rgiant**2 not overflow. the constants
* given here are suitable for every known computer.
*
* the function statement is
*
* double precision function enorm(n,x)
*
* where
*
* n is a positive integer input variable.
*
* x is an input array of length n.
*
* subprograms called
*
* fortran-supplied ... dabs,dsqrt
*
* argonne national laboratory. minpack project. march 1980.
* burton s. garbow, kenneth e. hillstrom, jorge j. more
*
* **********
*/
int i;
double agiant,floatn,s1,s2,s3;
double xabs;
double x1max, x3max;
double ans, temp;
s1 = zero;
s2 = zero;
s3 = zero;
x1max = zero;
x3max = zero;
floatn = n;
agiant = rgiant/floatn;
for( i=0; i rdwarf) && (xabs < agiant) )
{
/*
* sum for intermediate components.
*/
s2 += xabs*xabs;
continue;
}
if(xabs > rdwarf)
{
/*
* sum for large components.
*/
if(xabs > x1max)
{
temp = x1max/xabs;
s1 = one + s1*temp*temp;
x1max = xabs;
}
else
{
temp = xabs/x1max;
s1 += temp*temp;
}
continue;
}
/*
* sum for small components.
*/
if(xabs > x3max)
{
temp = x3max/xabs;
s3 = one + s3*temp*temp;
x3max = xabs;
}
else
{
if(xabs != zero)
{
temp = xabs/x3max;
s3 += temp*temp;
}
}
}
/*
* calculation of norm.
*/
if(s1 != zero)
{
temp = s1 + (s2/x1max)/x1max;
ans = x1max*Math.sqrt(temp);
return(ans);
}
if(s2 != zero)
{
if(s2 >= x3max) {
temp = s2*(one+(x3max/s2)*(x3max*s3));
} else {
temp = x3max*((s2/x3max)+(x3max*s3));
}
ans = Math.sqrt(temp);
}
else
{
ans = x3max*Math.sqrt(s3);
}
return(ans);
}
// void
// fdjac2(final int m,final int n,double* x,double* fvec,double* fjac,int,
// int* iflag,double epsfcn,double* wa,
// const QuantLib::MINPACK::LmdifCostFunction& fcn)
// {
// /*
// * **********
// *
// * subroutine fdjac2
// *
// * this subroutine computes a forward-difference approximation
// * to the m by n jacobian matrix associated with a specified
// * problem of m functions in n variables.
// *
// * the subroutine statement is
// *
// * subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
// *
// * where
// *
// * fcn is the name of the user-supplied subroutine which
// * calculates the functions. fcn must be declared
// * in an external statement in the user calling
// * program, and should be written as follows.
// *
// * subroutine fcn(m,n,x,fvec,iflag)
// * integer m,n,iflag
// * double precision x(n),fvec(m)
// * ----------
// * calculate the functions at x and
// * return this vector in fvec.
// * ----------
// * return
// * end
// *
// * the value of iflag should not be changed by fcn unless
// * the user wants to terminate execution of fdjac2.
// * in this case set iflag to a negative integer.
// *
// * m is a positive integer input variable set to the number
// * of functions.
// *
// * n is a positive integer input variable set to the number
// * of variables. n must not exceed m.
// *
// * x is an input array of length n.
// *
// * fvec is an input array of length m which must contain the
// * functions evaluated at x.
// *
// * fjac is an output m by n array which contains the
// * approximation to the jacobian matrix evaluated at x.
// *
// * ldfjac is a positive integer input variable not less than m
// * which specifies the leading dimension of the array fjac.
// *
// * iflag is an integer variable which can be used to terminate
// * the execution of fdjac2. see description of fcn.
// *
// * epsfcn is an input variable used in determining a suitable
// * step length for the forward-difference approximation. this
// * approximation assumes that the relative errors in the
// * functions are of the order of epsfcn. if epsfcn is less
// * than the machine precision, it is assumed that the relative
// * errors in the functions are of the order of the machine
// * precision.
// *
// * wa is a work array of length m.
// *
// * subprograms called
// *
// * user-supplied ...... fcn
// *
// * minpack-supplied ... dpmpar
// *
// * fortran-supplied ... dabs,dmax1,dsqrt
// *
// * argonne national laboratory. minpack project. march 1980.
// * burton s. garbow, kenneth e. hillstrom, jorge j. more
// *
// **********
// */
// int i,j,ij;
// double eps,h,temp;
// static double zero = 0.0;
//
//
// temp = dmax1(final epsfcn,MACHEP);
// eps = Math.sqrt(temp);
// ij = 0;
// for( j=0; j rdiag[kmax]) {
kmax = k;
}
}
if(kmax != j) {
ij = m * j;
jj = m * kmax;
for( i=0; i kp1 )
{
ik = kk + 1;
for( i=kp1; i= 1) {
for( k=0; k jp1)
{
ij = jp1 + ldr * j;
for( i=jp1; i= 1)
// {
// for( k=0; k= 0)
// {
// ij = ldr * j;
// for( i=0; i<=jm1; i++ )
// {
// wa1[i] -= r[ij]*temp;
// ij += 1;
// }
// }
// }
// }
//
// for( j=0; j= n)
// {
// for( j=0; j= 0)
// {
// ij = jj;
// for( i=0; i<=jm1; i++ )
// {
// sum += r[ij]*wa1[i];
// ij += 1;
// }
// }
// wa1[j] = (wa1[j] - sum)/r[j+ldr*j];
// jj += ldr; /* [i+ldr*j] */
// }
// temp = enorm(n,wa1);
// parl = ((fp/delta)/temp)/temp;
// }
// /*
// * calculate an upper bound, paru, for the zero of the function.
// */
// jj = 0;
// for( j=0; j zero)
// parl = dmax1(parl, *par);
// if(fp < zero)
// paru = dmin1(paru, *par);
// /*
// * compute an improved estimate for par.
// */
// *par = dmax1(parl, *par + parc);
// /*
// * end of an iteration.
// */
// goto L150;
//
// L220:
// /*
// * termination.
// */
// if(iter == 0)
// *par = zero;
// /*
// * last card of subroutine lmpar.
// */
// }
//
//
//
//
//
// void lmdif(final int m,final int n,double* x,double* fvec,double ftol,
// double xtol,double gtol,int maxfev,double epsfcn,
// double* diag, int mode, double factor,
// int nprint, int* info,int* nfev,double* fjac,
// int ldfjac,int* ipvt,double* qtf,
// double* wa1,double* wa2,double* wa3,double* wa4,
// const QuantLib::MINPACK::LmdifCostFunction& fcn)
// {
// /*
// * **********
// *
// * subroutine lmdif
// *
// * the purpose of lmdif is to minimize the sum of the squares of
// * m nonlinear functions in n variables by a modification of
// * the levenberg-marquardt algorithm. the user must provide a
// * subroutine which calculates the functions. the jacobian is
// * then calculated by a forward-difference approximation.
// *
// * the subroutine statement is
// *
// * subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn,
// * diag,mode,factor,nprint,info,nfev,fjac,
// * ldfjac,ipvt,qtf,wa1,wa2,wa3,wa4)
// *
// * where
// *
// * fcn is the name of the user-supplied subroutine which
// * calculates the functions. fcn must be declared
// * in an external statement in the user calling
// * program, and should be written as follows.
// *
// * subroutine fcn(m,n,x,fvec,iflag)
// * integer m,n,iflag
// * double precision x(n),fvec(m)
// * ----------
// * calculate the functions at x and
// * return this vector in fvec.
// * ----------
// * return
// * end
// *
// * the value of iflag should not be changed by fcn unless
// * the user wants to terminate execution of lmdif.
// * in this case set iflag to a negative integer.
// *
// * m is a positive integer input variable set to the number
// * of functions.
// *
// * n is a positive integer input variable set to the number
// * of variables. n must not exceed m.
// *
// * x is an array of length n. on input x must contain
// * an initial estimate of the solution vector. on output x
// * contains the final estimate of the solution vector.
// *
// * fvec is an output array of length m which contains
// * the functions evaluated at the output x.
// *
// * ftol is a nonnegative input variable. termination
// * occurs when both the actual and predicted relative
// * reductions in the sum of squares are at most ftol.
// * therefore, ftol measures the relative error desired
// * in the sum of squares.
// *
// * xtol is a nonnegative input variable. termination
// * occurs when the relative error between two consecutive
// * iterates is at most xtol. therefore, xtol measures the
// * relative error desired in the approximate solution.
// *
// * gtol is a nonnegative input variable. termination
// * occurs when the cosine of the angle between fvec and
// * any column of the jacobian is at most gtol in absolute
// * value. therefore, gtol measures the orthogonality
// * desired between the function vector and the columns
// * of the jacobian.
// *
// * maxfev is a positive integer input variable. termination
// * occurs when the number of calls to fcn is at least
// * maxfev by the end of an iteration.
// *
// * epsfcn is an input variable used in determining a suitable
// * step length for the forward-difference approximation. this
// * approximation assumes that the relative errors in the
// * functions are of the order of epsfcn. if epsfcn is less
// * than the machine precision, it is assumed that the relative
// * errors in the functions are of the order of the machine
// * precision.
// *
// * diag is an array of length n. if mode = 1 (see
// * below), diag is internally set. if mode = 2, diag
// * must contain positive entries that serve as
// * multiplicative scale factors for the variables.
// *
// * mode is an integer input variable. if mode = 1, the
// * variables will be scaled internally. if mode = 2,
// * the scaling is specified by the input diag. other
// * values of mode are equivalent to mode = 1.
// *
// * factor is a positive input variable used in determining the
// * initial step bound. this bound is set to the product of
// * factor and the euclidean norm of diag*x if nonzero, or else
// * to factor itself. in most cases factor should lie in the
// * interval (.1,100.). 100. is a generally recommended value.
// *
// * nprint is an integer input variable that enables controlled
// * printing of iterates if it is positive. in this case,
// * fcn is called with iflag = 0 at the beginning of the first
// * iteration and every nprint iterations thereafter and
// * immediately prior to return, with x and fvec available
// * for printing. if nprint is not positive, no special calls
// * of fcn with iflag = 0 are made.
// *
// * info is an integer output variable. if the user has
// * terminated execution, info is set to the (negative)
// * value of iflag. see description of fcn. otherwise,
// * info is set as follows.
// *
// * info = 0 improper input parameters.
// *
// * info = 1 both actual and predicted relative reductions
// * in the sum of squares are at most ftol.
// *
// * info = 2 relative error between two consecutive iterates
// * is at most xtol.
// *
// * info = 3 conditions for info = 1 and info = 2 both hold.
// *
// * info = 4 the cosine of the angle between fvec and any
// * column of the jacobian is at most gtol in
// * absolute value.
// *
// * info = 5 number of calls to fcn has reached or
// * exceeded maxfev.
// *
// * info = 6 ftol is too small. no further reduction in
// * the sum of squares is possible.
// *
// * info = 7 xtol is too small. no further improvement in
// * the approximate solution x is possible.
// *
// * info = 8 gtol is too small. fvec is orthogonal to the
// * columns of the jacobian to machine precision.
// *
// * nfev is an integer output variable set to the number of
// * calls to fcn.
// *
// * fjac is an output m by n array. the upper n by n submatrix
// * of fjac contains an upper triangular matrix r with
// * diagonal elements of nonincreasing magnitude such that
// *
// * t t t
// * p *(jac *jac)*p = r *r,
// *
// * where p is a permutation matrix and jac is the final
// * calculated jacobian. column j of p is column ipvt(j)
// * (see below) of the identity matrix. the lower trapezoidal
// * part of fjac contains information generated during
// * the computation of r.
// *
// * ldfjac is a positive integer input variable not less than m
// * which specifies the leading dimension of the array fjac.
// *
// * ipvt is an integer output array of length n. ipvt
// * defines a permutation matrix p such that jac*p = q*r,
// * where jac is the final calculated jacobian, q is
// * orthogonal (not stored), and r is upper triangular
// * with diagonal elements of nonincreasing magnitude.
// * column j of p is column ipvt(j) of the identity matrix.
// *
// * qtf is an output array of length n which contains
// * the first n elements of the vector (q transpose)*fvec.
// *
// * wa1, wa2, and wa3 are work arrays of length n.
// *
// * wa4 is a work array of length m.
// *
// * subprograms called
// *
// * user-supplied ...... fcn
// *
// * minpack-supplied ... dpmpar,enorm,fdjac2,lmpar,qrfac
// *
// * fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod
// *
// * argonne national laboratory. minpack project. march 1980.
// * burton s. garbow, kenneth e. hillstrom, jorge j. more
// *
// * **********
// */
// int i,iflag,ij,jj,iter,j,l;
// double actred,delta=0,dirder,fnorm,fnorm1,gnorm;
// double par,pnorm,prered,ratio;
// double sum,temp,temp1,temp2,temp3,xnorm=0;
// static double one = 1.0;
// static double p1 = 0.1;
// static double p5 = 0.5;
// static double p25 = 0.25;
// static double p75 = 0.75;
// static double p0001 = 1.0e-4;
// static double zero = 0.0;
//
// *info = 0;
// iflag = 0;
// *nfev = 0;
// /*
// * check the input parameters for errors.
// */
// if( (n <= 0) || (m < n) || (ldfjac < m) || (ftol < zero)
// || (xtol < zero) || (gtol < zero) || (maxfev <= 0)
// || (factor <= zero) )
// goto L300;
//
// if( mode == 2 )
// { /* scaling by diag[] */
// for( j=0; j 0 )
// {
// iflag = 0;
// if(mod(iter-1,nprint) == 0)
// {
// fcn(m,n,x,fvec,&iflag);
// if(iflag < 0)
// goto L300;
// }
// }
// /*
// * compute the qr factorization of the jacobian.
// */
// qrfac(m,n,fjac,ldfjac,1,ipvt,n,wa1,wa2,wa3);
// /*
// * on the first iteration and if mode is 1, scale according
// * to the norms of the columns of the initial jacobian.
// */
// if(iter == 1)
// {
// if(mode != 2)
// {
// for( j=0; j= zero)
// temp = p5;
// else
// temp = p5*dirder/(dirder + p5*actred);
// if( ((p1*fnorm1) >= fnorm)
// || (temp < p1) )
// temp = p1;
// delta = temp*dmin1(delta,pnorm/p1);
// par = par/temp;
// }
// else
// {
// if( (par == zero) || (ratio >= p75) )
// {
// delta = pnorm/p5;
// par = p5*par;
// }
// }
// /*
// * test for successful iteration.
// */
// if(ratio >= p0001)
// {
// /*
// * successful iteration. update x, fvec, and their norms.
// */
// for( j=0; j= maxfev)
// *info = 5;
// if( (Math.fabs(actred) <= MACHEP)
// && (prered <= MACHEP)
// && (p5*ratio <= one) )
// *info = 6;
// if(delta <= MACHEP*xnorm)
// *info = 7;
// if(gnorm <= MACHEP)
// *info = 8;
// if( *info != 0)
// goto L300;
// /*
// * end of the inner loop. repeat if iteration unsuccessful.
// */
// if(ratio < p0001)
// goto L200;
// /*
// * end of the outer loop.
// */
// goto L30;
//
// L300:
// /*
// * termination, either normal or user imposed.
// */
// if(iflag < 0)
// *info = iflag;
// iflag = 0;
// if(nprint > 0)
// fcn(m,n,x,fvec,&iflag);
// /*
// last card of subroutine lmdif.
// */
// }
private static double dmax1(final double a,final double b)
{
if( a >= b )
return(a);
else
return(b);
}
private double dmin1(final double a,final double b)
{
if( a <= b )
return(a);
else
return(b);
}
private static int min0(final int a,final int b)
{
if( a <= b )
return(a);
else
return(b);
}
private int mod( final int k, final int m )
{
return( k % m );
}
/***********Sample of user supplied function****************
* m = number of functions
* n = number of variables
* x = vector of function arguments
* fvec = vector of function values
* iflag = error return variable
*/
//void fcn(int m,int n, double* x, double* fvec,int *iflag)
//{
// QuantLib::LevenbergMarquardt::fcn(m, n, x, fvec, iflag);
//}
} // private static class MinpackC
}
Ads help maintain this website.