lbfgsp.h

/*
* All of the documentation and software included in the
* Alchemy Software is copyrighted by Stanley Kok, Parag
* Singla, Matthew Richardson, Pedro Domingos, Marc
* Sumner and Hoifung Poon.
* 
* Copyright [2004-06] Stanley Kok, Parag Singla, Matthew
* Richardson, Pedro Domingos, Marc Sumner and Hoifung
* Poon. All rights reserved.
* 
* Contact: Pedro Domingos, University of Washington
* (pedrod@cs.washington.edu).
* 
* 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. All advertising materials mentioning features or use
* of this software must display the following
* acknowledgment: "This product includes software
* developed by Stanley Kok, Parag Singla, Matthew
* Richardson, Pedro Domingos, Marc Sumner and Hoifung
* Poon in the Department of Computer Science and
* Engineering at the University of Washington".
* 
* 4. Your publications acknowledge the use or
* contribution made by the Software to your research
* using the following citation(s): 
* Stanley Kok, Parag Singla, Matthew Richardson and
* Pedro Domingos (2005). "The Alchemy System for
* Statistical Relational AI", Technical Report,
* Department of Computer Science and Engineering,
* University of Washington, Seattle, WA.
* http://www.cs.washington.edu/ai/alchemy.
* 
* 5. Neither the name of the University of Washington nor
* the names of its contributors may be used to endorse or
* promote products derived from this software without
* specific prior written permission.
* 
* THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF WASHINGTON
* AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
* WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
* PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
* OF WASHINGTON OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
* INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
* PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
* ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN
* IF ADVISED OF THE POSS IBILITY OF SUCH DAMAGE.
* 
*/
#ifndef LBFGSP_H_SEP_17_2005
#define LBFGSP_H_SEP_17_2005

#include <cmath>
#include <cstdlib>
#include <fstream>
#include "timer.h"
#include "Polynomial.h"

const int LBFGSP_PRINT = -1; //-1: no output; 1: output for every iteration

// The code below from setulb() onwards is ported from: 
// C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a
// limited memory FORTRAN code for solving bound constrained
// optimization problems'', Tech. Report, NAM-11, EECS Department,
// Northwestern University.

class LBFGSP
{
public:
        LBFGSP(const int& maxIter, const double& ftol, 
                const int& numWts) 
                : maxIter_(maxIter), ftol_(ftol), 
                nbd_(NULL), iwa_(NULL), l_(NULL), u_(NULL), g_(NULL), wa_(NULL)
        { init(numWts); }

        LBFGSP() 
                : maxIter_(-1), ftol_(-1), 
                nbd_(NULL), iwa_(NULL), l_(NULL), u_(NULL), g_(NULL), wa_(NULL)
        { init(0); }


        ~LBFGSP() { destroy(); }


        void setMaxIter(const int& iter) { maxIter_ = iter; }

        void setFtol(const double& tol)  { ftol_ = tol; }


        void init(const int& numWts)
        {
                destroy();
                numWts_ = numWts;
                int mmax = 17;
                nbd_ = new int[numWts_+1];
                iwa_ = new int[3*numWts_+1];
                l_   = new double[numWts_+1];
                u_   = new double[numWts_+1];
                g_   = new double[numWts_+1];
                wa_  = new double[2*mmax*numWts_+4*numWts_+12*mmax*mmax+12*mmax+1];    
        }


        void reInit(const int& numWts)
        {
                if (numWts_ == numWts) return;
                init(numWts);
        }

        void destroy() 
        {
                if (nbd_ != NULL) delete [] nbd_;
                if (iwa_ != NULL) delete [] iwa_;
                if (l_   != NULL) delete [] l_;
                if (u_   != NULL) delete [] u_;
                if (g_   != NULL) delete [] g_;
                if (wa_  != NULL) delete [] wa_;
        }

        double minimize(const int& numWts, double* const & wts, int& iter,bool& error, PolyNomial& pl)
        {
                reInit(numWts);
                return minimize(iter, error, pl);
        }

        double startLBFGS(PolyNomial& pl)
        {
                int m = 5; //max number of limited memory corrections
                double f; // value of function to be optimized
                double factr = 0;
                double pgtol = 0;
                // -1: silent (-1); 1: out at every iteration
                ofstream* itfile = NULL;
                int iprint = LBFGSP_PRINT; 
                if (iprint >= 1) itfile = new ofstream("iterate.dat");
                //indicate that the elements of x[] are unbounded
                //for (int i = 0; i <= numWts_; i++) nbd_[i] = 0;

                //indicate that the elements of x[] have both upper bounds and lower bounds
                for (int i = 0; i <= numWts_; i++) nbd_[i] = 2;

                strcpy(task_,"START");
                for (int i = 5; i <= 60; i++) task_[i]=' ';


                double* wts = new double[numWts_ + 1];
                // set the initial value of pl here
                pl.ClearVarValue();
                for(int i = 0; i < numWts_; i++) 
                {
                        wts[i+1] = 0;
                        pl.AppendVarValue(wts[i+1]);
                }

                setulb(numWts_,m,wts,l_,u_,nbd_,f,g_,factr,pgtol,wa_,iwa_,task_,
                        iprint,csave_,lsave_,isave_,dsave_,itfile);
                pl.ClearVarValue();
                for(int i = 0; i < numWts_; i++) pl.AppendVarValue(wts[i+1]);
                delete[] wts;
                return f;
        }
        // assume the lower bound and upper bound are set
        // assume the initial value of variables are already in pl
        double proceedOneStep(PolyNomial& pl)
        {
        int m = 5; //max number of limited memory corrections
                double f = -999; // value of function to be optimized
                double factr = 0;
                double pgtol = 0;
                // -1: silent (-1); 1: out at every iteration
                ofstream* itfile = NULL;
                int iprint = LBFGSP_PRINT; 
                if (iprint >= 1) itfile = new ofstream("iterate.dat");
                double* wts = new double[numWts_ + 1];
                for(int i = 0; i < numWts_; i++) {
                        wts[i+1] = pl.GetVarValue(i);
                }
                //iter = 0;
                //while routine returns "FG" or "NEW_X" in task, keep calling it
                if (strncmp(task_,"FG",2)==0)
                {
                        f = getValueAndGradient(g_, wts, pl);
                        setulb(numWts_,m,wts,l_,u_,nbd_,f,g_,factr,pgtol,wa_,iwa_,task_,
                                iprint,csave_,lsave_,isave_,dsave_,itfile);
                        pl.ClearVarValue();
                        for(int i = 0; i < numWts_; i++) pl.AppendVarValue(wts[i+1]);
                }
                else
                {
                        //the minimization routine has returned with a new iterate,
                        //and we have opted to continue the iteration
                        setulb(numWts_,m,wts,l_,u_,nbd_,f,g_,factr,pgtol,wa_,iwa_,task_,
                                iprint,csave_,lsave_,isave_,dsave_,itfile);

                        pl.ClearVarValue();
                        for(int i = 0; i < numWts_; i++) pl.AppendVarValue(wts[i+1]);
                }
                //If task is neither FG nor NEW_X we terminate execution.
                //the minimization routine has returned with one of
                //{CONV, ABNO, ERROR}
                delete[] wts;
                return pl.ComputePlValue();
        }


        // we assume the current variable values are stored in the pl itself
        double minimize(int& iter, bool& error, PolyNomial& pl)
        {
//              cout << "LBFGSP::minimize() begins" << endl;
                error = false;
                int m = 5; //max number of limited memory corrections
                double f; // value of function to be optimized
                double factr = 0;
                double pgtol = 0;
                // -1: silent (-1); 1: out at every iteration
                ofstream* itfile = NULL;
                int iprint = LBFGSP_PRINT; 
                if (iprint >= 1) itfile = new ofstream("iterate.dat");
                iter = 0;
                //indicate that the elements of x[] are unbounded
                //for (int i = 0; i <= numWts_; i++) nbd_[i] = 0;

                //indicate that the elements of x[] have both upper bounds and lower bounds
                for (int i = 0; i <= numWts_; i++) nbd_[i] = 2;

                strcpy(task_,"START");
                for (int i = 5; i <= 60; i++) task_[i]=' ';

                // set the initial value of pl here

                double* wts = new double[numWts_ + 1];
                pl.ClearVarValue();
                for(int i = 0; i < numWts_; i++) 
                {
                        wts[i + 1] = 0;
                        pl.AppendVarValue(wts[i+1]);
                }

                setulb(numWts_,m,wts,l_,u_,nbd_,f,g_,factr,pgtol,wa_,iwa_,task_,
                        iprint,csave_,lsave_,isave_,dsave_,itfile);

                pl.ClearVarValue();
                for(int i = 0; i < numWts_; i++) pl.AppendVarValue(wts[i+1]);

                double initialValue = 0, prevValue = 0, newValue;
                bool firstIter = true;

                //while routine returns "FG" or "NEW_X" in task, keep calling it
                while (strncmp(task_,"FG",2)==0 || strncmp(task_,"NEW_X",5)==0)
                {
                        if (strncmp(task_,"FG",2)==0)
                        {
                                f = getValueAndGradient(g_, wts, pl);

                                setulb(numWts_,m,wts,l_,u_,nbd_,f,g_,factr,pgtol,wa_,iwa_,task_,
                                        iprint,csave_,lsave_,isave_,dsave_,itfile);

                                pl.ClearVarValue();
                                for(int i = 0; i < numWts_; i++) pl.AppendVarValue(wts[i+1]);

                                if (firstIter) { firstIter = false; prevValue = f; initialValue = f; }
                        }
                        else
                        {
                                //the minimization routine has returned with a new iterate,
                                //and we have opted to continue the iteration
                                if (iter+1 > maxIter_) break;
                                ++iter;
                                newValue = f;

                                if (fabs(newValue-prevValue) < ftol_*fabs(prevValue)) break;
                                prevValue = newValue;
                                setulb(numWts_,m,wts,l_,u_,nbd_,f,g_,factr,pgtol,wa_,iwa_,task_,
                                        iprint,csave_,lsave_,isave_,dsave_,itfile);
                                pl.ClearVarValue();
                                for(int i = 0; i < numWts_; i++) pl.AppendVarValue(wts[i+1]);
                        }
                }

                //If task is neither FG nor NEW_X we terminate execution.
                //the minimization routine has returned with one of
                //{CONV, ABNO, ERROR}

                if (strncmp(task_,"ABNO",4) == 0)
                {
                        cout << "ERROR: LBFGSP failed. Returned ABNO" << endl;
                        error = true;
                        delete[] wts;
                        return initialValue;
                }

                if (strncmp(task_,"ERROR",5)==0)
                {
                        cout << "ERROR: LBFGSP failed. Returned ERROR" << endl;
                        error = true;
                        delete[] wts;
                        return initialValue;
                }

                if (strncmp(task_,"CONV",4)==0)
                {
                        cout << "LBFGSP converged!" << endl;
                }

                delete[] wts;
//              cout << "LBFGSP::minimize() ends" << endl;
                return f;
        }


        void setUpperBounds(const double* u)
        {
                if (numWts_ > 0)
                {
                        memcpy(u_+1, u, numWts_*sizeof(double));
                }       
        }

        void setLowerBounds(const double* l)
        {
                if(numWts_ > 0)
                {
                        memcpy(l_+1, l, numWts_*sizeof(double));
                }
        }

private:
        int maxIter_;
        double ftol_;
        int numWts_;

        // params of lbfgsb algorithm
        int*    nbd_;
        int*    iwa_;
        double* l_;
        double* u_;
        double* g_;
        double* wa_;

        char task_[61];
        char csave_[61];
        bool lsave_[5];
        int  isave_[45];
        double dsave_[30];


private :
        // g for storing gradient, wts for storing current assignment to variables
        double getValueAndGradient(double* const & g, const double* const & wts, PolyNomial& pl)
        {
                g[0] = 0;
                //set pl var values
                //get Pl gradient 
                Array<double> gr = pl.GetGradient();
                for(int i = 0; i < gr.size(); i++)
                {
                        g[i+1] = gr[i];
                }

                return pl.ComputePlValue();
        }
        /*double getValueAndGradient(double* const & g, const double* const & wts)
        {
                return (pll_->getValueAndGradient(g+1, wts+1, numWts_));
        }*/

        //function
        int getIdx(const int& i, const int& j, const int& idim)
        {
                return (j-1)*idim + i;
        }


        double max(const double& a, const double& b) { if (a>=b) return a; return b; }

        double min(const double& a, const double& b) { if (a<=b) return a; return b; }

        int min(const int& a, const int& b) { if (a<=b) return a; return b; }

        double max(const double& a, const double& b, const double& c)
        {
                if (a >= b && a >= c) return a;
                if (b >= a && b >= c) return b;
                assert(c >= a && c >= b);
                return c;
        }


        // The code below is ported from: 
        // C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a
        // limited memory FORTRAN code for solving bound constrained
        // optimization problems'', Tech. Report, NAM-11, EECS Department,
        // Northwestern University.

        // The comments before each function are from the original fortran code.

        //   ************
        //
        //   Subroutine setulb
        //
        //   This subroutine partitions the working arrays wa and iwa, and 
        //     then uses the limited memory BFGS method to solve the bound
        //     constrained optimization problem by calling mainlb.
        //     (The direct method will be used in the subspace minimization.)
        //
        //   n is an integer variable.
        //     On entry n is the dimension of the problem.
        //     On exit n is unchanged.
        //
        //   m is an integer variable.
        //     On entry m is the maximum number of variable metric corrections
        //       used to define the limited memory matrix.
        //     On exit m is unchanged.
        //
        //   x is a double precision array of dimension n.
        //     On entry x is an approximation to the solution.
        //     On exit x is the current approximation.
        //
        //   l is a double precision array of dimension n.
        //     On entry l is the lower bound on x.
        //     On exit l is unchanged.
        //
        //   u is a double precision array of dimension n.
        //     On entry u is the upper bound on x.
        //     On exit u is unchanged.
        //
        //   nbd is an integer array of dimension n.
        //     On entry nbd represents the type of bounds imposed on the
        //       variables, and must be specified as follows:
        //       nbd(i)=0 if x(i) is unbounded,
        //              1 if x(i) has only a lower bound,
        //              2 if x(i) has both lower and upper bounds, and
        //              3 if x(i) has only an upper bound.
        //     On exit nbd is unchanged.
        //
        //   f is a double precision variable.
        //     On first entry f is unspecified.
        //     On final exit f is the value of the function at x.
        //
        //   g is a double precision array of dimension n.
        //     On first entry g is unspecified.
        //     On final exit g is the value of the gradient at x.
        //
        //   factr is a double precision variable.
        //     On entry factr >= 0 is specified by the user.  The iteration
        //       will stop when
        //
        //       (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch
        //
        //       where epsmch is the machine precision, which is automatically
        //       generated by the code. Typical values for factr: 1.d+12 for
        //       low accuracy; 1.d+7 for moderate accuracy; 1.d+1 for extremely
        //       high accuracy.
        //     On exit factr is unchanged.
        //
        //   pgtol is a double precision variable.
        //     On entry pgtol >= 0 is specified by the user.  The iteration
        //       will stop when
        //
        //               max{|proj g_i | i = 1, ..., n} <= pgtol
        //
        //       where pg_i is the ith component of the projected gradient.   
        //     On exit pgtol is unchanged.
        //
        //   wa is a double precision working array of length 
        //     (2mmax + 4)nmax + 12mmax^2 + 12mmax.
        //
        //   iwa is an integer working array of length 3nmax.
        //
        //   task is a working string of characters of length 60 indicating
        //     the current job when entering and quitting this subroutine.
        //
        //   iprint is an integer variable that must be set by the user.
        //     It controls the frequency and type of output generated:
        //      iprint<0    no output is generated;
        //      iprint=0    print only one line at the last iteration;
        //      0<iprint<99 print also f and |proj g| every iprint iterations;
        //      iprint=99   print details of every iteration except n-vectors;
        //      iprint=100  print also the changes of active set and final x;
        //      iprint>100  print details of every iteration including x and g;
        //     When iprint > 0, the file iterate.dat will be created to
        //                      summarize the iteration.
        //
        //   csave is a working string of characters of length 60.
        //
        //   lsave is a logical working array of dimension 4.
        //     On exit with 'task' = NEW_X, the following information is 
        //                                                           available:
        //       If lsave(1) = true  then  the initial X has been replaced by
        //                                   its projection in the feasible set;
        //       If lsave(2) = true  then  the problem is constrained;
        //       If lsave(3) = true  then  each variable has upper and lower
        //                                   bounds;
        //
        //   isave is an integer working array of dimension 44.
        //     On exit with 'task' = NEW_X, the following information is 
        //                                                           available:
        //       isave(22) = the total number of intervals explored in the 
        //                       search of Cauchy points;
        //       isave(26) = the total number of skipped BFGS updates before 
        //                       the current iteration;
        //       isave(30) = the number of current iteration;
        //       isave(31) = the total number of BFGS updates prior the current
        //                       iteration;
        //       isave(33) = the number of intervals explored in the search of
        //                       Cauchy point in the current iteration;
        //       isave(34) = the total number of function and gradient 
        //                       evaluations;
        //       isave(36) = the number of function value or gradient
        //                                evaluations in the current iteration;
        //       if isave(37) = 0  then the subspace argmin is within the box;
        //       if isave(37) = 1  then the subspace argmin is beyond the box;
        //       isave(38) = the number of free variables in the current
        //                       iteration;
        //       isave(39) = the number of active constraints in the current
        //                       iteration;
        //       n + 1 - isave(40) = the number of variables leaving the set of
        //                         active constraints in the current iteration;
        //       isave(41) = the number of variables entering the set of active
        //                       constraints in the current iteration.
        //
        //   dsave is a double precision working array of dimension 29.
        //     On exit with 'task' = NEW_X, the following information is
        //                                                           available:
        //       dsave(1) = current 'theta' in the BFGS matrix;
        //       dsave(2) = f(x) in the previous iteration;
        //       dsave(3) = factr*epsmch;
        //       dsave(4) = 2-norm of the line search direction vector;
        //       dsave(5) = the machine precision epsmch generated by the code;
        //       dsave(7) = the accumulated time spent on searching for
        //                                                       Cauchy points;
        //       dsave(8) = the accumulated time spent on
        //                                               subspace minimization;
        //       dsave(9) = the accumulated time spent on line search;
        //       dsave(11) = the slope of the line search function at
        //                                the current point of line search;
        //       dsave(12) = the maximum relative step length imposed in
        //                                                         line search;
        //       dsave(13) = the infinity norm of the projected gradient;
        //       dsave(14) = the relative step length in the line search;
        //       dsave(15) = the slope of the line search function at
        //                               the starting point of the line search;
        //       dsave(16) = the square of the 2-norm of the line search
        //                                                    direction vector.
        //
        //   Subprograms called:
        //
        //     L-BFGS-B Library ... mainlb.    
        //
        //
        //   References:
        //
        //     [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
        //     memory algorithm for bound constrained optimization'',
        //     SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
        //
        //     [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a
        //     limited memory FORTRAN code for solving bound constrained
        //     optimization problems'', Tech. Report, NAM-11, EECS Department,
        //     Northwestern University, 1994.
        //
        //     (Postscript files of these papers are available via anonymous
        //      ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        // ************
        //function
        void setulb(const int& n, const int& m, double* const & x, 
                const double* const & l, const double* const & u, 
                const int* const & nbd, double& f, double* const & g, 
                const double& factr, const double& pgtol, 
                double* const & wa, int* const & iwa,
                char* const & task, const int& iprint,  
                char* const & csave, bool* const & lsave, 
                int* const & isave, double* const & dsave,
                ofstream* itfile)
        { 
                int l1,l2,l3,lws,lr,lz,lt,ld,lsg,lwa,lyg,
                        lsgo,lwy,lsy,lss,lyy,lwt,lwn,lsnd,lygo;

                if (strncmp(task,"START",5)==0)
                {
                        isave[1]  = m*n;
                        isave[2]  = m*m;
                        isave[3]  = 4*m*m;
                        isave[4]  = 1;
                        isave[5]  = isave[4]  + isave[1];
                        isave[6]  = isave[5]  + isave[1];
                        isave[7]  = isave[6]  + isave[2];
                        isave[8]  = isave[7]  + isave[2];
                        isave[9]  = isave[8]  + isave[2];
                        isave[10] = isave[9]  + isave[2];
                        isave[11] = isave[10] + isave[3];
                        isave[12] = isave[11] + isave[3];
                        isave[13] = isave[12] + n;
                        isave[14] = isave[13] + n;
                        isave[15] = isave[14] + n;
                        isave[16] = isave[15] + n;
                        isave[17] = isave[16] + 8*m;
                        isave[18] = isave[17] + m;
                        isave[19] = isave[18] + m;
                        isave[20] = isave[19] + m;
                }
                l1   = isave[1];
                l2   = isave[2];
                l3   = isave[3];
                lws  = isave[4];
                lwy  = isave[5];
                lsy  = isave[6];
                lss  = isave[7];
                lyy  = isave[8];
                lwt  = isave[9];
                lwn  = isave[10];
                lsnd = isave[11];
                lz   = isave[12];
                lr   = isave[13];
                ld   = isave[14];
                lt   = isave[15];
                lwa  = isave[16];
                lsg  = isave[17];
                lsgo = isave[18];
                lyg  = isave[19];
                lygo = isave[20];

                timer_.reset();

                mainlb(n,m,x,l,u,nbd,f,g,factr,pgtol,
                        &(wa[lws-1]),&(wa[lwy-1]),&(wa[lsy-1]),&(wa[lss-1]),&(wa[lyy-1]),
                        &(wa[lwt-1]),&(wa[lwn-1]),&(wa[lsnd-1]),&(wa[lz-1]),&(wa[lr-1]),
                        &(wa[ld-1]),&(wa[lt-1]),&(wa[lwa-1]),&(wa[lsg-1]),&(wa[lsgo-1]),
                        &(wa[lyg-1]),&(wa[lygo-1]),&(iwa[1-1]),&(iwa[n+1-1]),&(iwa[2*n+1-1]),
                        task,iprint,csave,lsave,&(isave[22-1]),dsave, itfile);
        } //setulb()


        //   ************
        //
        //   Subroutine mainlb
        //
        //   This subroutine solves bound constrained optimization problems by
        //     using the compact formula of the limited memory BFGS updates.
        //     
        //   n is an integer variable.
        //     On entry n is the number of variables.
        //     On exit n is unchanged.
        //
        //   m is an integer variable.
        //     On entry m is the maximum number of variable metri//
        //     corrections allowed in the limited memory matrix.
        //     On exit m is unchanged.
        //
        //   x is a double precision array of dimension n.
        //     On entry x is an approximation to the solution.
        //     On exit x is the current approximation.
        //
        //   l is a double precision array of dimension n.
        //     On entry l is the lower bound of x.
        //     On exit l is unchanged.
        //
        //   u is a double precision array of dimension n.
        //     On entry u is the upper bound of x.
        //     On exit u is unchanged.
        //
        //   nbd is an integer array of dimension n.
        //     On entry nbd represents the type of bounds imposed on the
        //       variables, and must be specified as follows:
        //       nbd(i)=0 if x(i) is unbounded,
        //              1 if x(i) has only a lower bound,
        //              2 if x(i) has both lower and upper bounds,
        //              3 if x(i) has only an upper bound.
        //     On exit nbd is unchanged.
        //
        //   f is a double precision variable.
        //     On first entry f is unspecified.
        //     On final exit f is the value of the function at x.
        //
        //   g is a double precision array of dimension n.
        //     On first entry g is unspecified.
        //     On final exit g is the value of the gradient at x.
        //
        //   factr is a double precision variable.
        //     On entry factr >= 0 is specified by the user.  The iteration
        //       will stop when
        //
        //       (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch
        //
        //       where epsmch is the machine precision, which is automatically
        //       generated by the code.
        //     On exit factr is unchanged.
        //
        //   pgtol is a double precision variable.
        //     On entry pgtol >= 0 is specified by the user.  The iteration
        //       will stop when
        //
        //               max{|proj g_i | i = 1, ..., n} <= pgtol
        //
        //       where pg_i is the ith component of the projected gradient.
        //     On exit pgtol is unchanged.
        //
        //   ws, wy, sy, and wt are double precision working arrays used to
        //     store the following information defining the limited memory
        //        BFGS matrix:
        //        ws, of dimension n x m, stores S, the matrix of s-vectors;
        //        wy, of dimension n x m, stores Y, the matrix of y-vectors;
        //        sy, of dimension m x m, stores S'Y;
        //        ss, of dimension m x m, stores S'S;
        //          yy, of dimension m x m, stores Y'Y;
        //        wt, of dimension m x m, stores the Cholesky factorization
        //                                of (theta*S'S+LD^(-1)L'); see eq.
        //                                (2.26) in [3].
        //
        //   wn is a double precision working array of dimension 2m x 2m
        //     used to store the LEL^T factorization of the indefinite matrix
        //               K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
        //                   [L_a -R_z           theta*S'AA'S ]
        //
        //     where     E = [-I  0]
        //                   [ 0  I]
        //
        //   snd is a double precision working array of dimension 2m x 2m
        //     used to store the lower triangular part of
        //               N = [Y' ZZ'Y   L_a'+R_z']
        //                   [L_a +R_z  S'AA'S   ]
        //           
        //   z(n),r(n),d(n),t(n),wa(8*m) are double precision working arrays.
        //     z is used at different times to store the Cauchy point and
        //     the Newton point.
        //
        //   sg(m),sgo(m),yg(m),ygo(m) are double precision working arrays. 
        //
        //   index is an integer working array of dimension n.
        //     In subroutine freev, index is used to store the free and fixed
        //        variables at the Generalized Cauchy Point (GCP).
        //
        //   iwhere is an integer working array of dimension n used to record
        //     the status of the vector x for GCP computation.
        //     iwhere(i)=0 or -3 if x(i) is free and has bounds,
        //               1       if x(i) is fixed at l(i), and l(i) != u(i)
        //               2       if x(i) is fixed at u(i), and u(i) != l(i)
        //               3       if x(i) is always fixed, i.e.,  u(i)=x(i)=l(i)
        //              -1       if x(i) is always free, i.e., no bounds on it.
        //
        //   indx2 is an integer working array of dimension n.
        //     Within subroutine cauchy, indx2 corresponds to the array iorder.
        //     In subroutine freev, a list of variables entering and leaving
        //     the free set is stored in indx2, and it is passed on to
        //     subroutine formk with this information.
        //
        //   task is a working string of characters of length 60 indicating
        //     the current job when entering and leaving this subroutine.
        //
        //   iprint is an INTEGER variable that must be set by the user.
        //     It controls the frequency and type of output generated:
        //      iprint<0    no output is generated;
        //      iprint=0    print only one line at the last iteration;
        //      0<iprint<99 print also f and |proj g| every iprint iterations;
        //      iprint=99   print details of every iteration except n-vectors;
        //      iprint=100  print also the changes of active set and final x;
        //      iprint>100  print details of every iteration including x and g;
        //     When iprint > 0, the file iterate.dat will be created to
        //                      summarize the iteration.
        //
        //   csave is a working string of characters of length 60.
        //
        //   lsave is a logical working array of dimension 4.
        //
        //   isave is an integer working array of dimension 23.
        //
        //   dsave is a double precision working array of dimension 29.
        //
        //
        //   Subprograms called
        //
        //     L-BFGS-B Library ... cauchy, subsm, lnsrlb, formk, 
        //
        //      errclb, prn1lb, prn2lb, prn3lb, active, projgr,
        //
        //      freev, cmprlb, matupd, formt.
        //
        //     Minpack2 Library ... timer, dpmeps.
        //
        //     Linpack Library ... dcopy, ddot.
        //
        //
        //   References:
        //
        //     [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
        //     memory algorithm for bound constrained optimization'',
        //     SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
        //
        //     [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN
        //     Subroutines for Large Scale Bound Constrained Optimization''
        //     Tech. Report, NAM-11, EECS Department, Northwestern University,
        //     1994.
        // 
        //     [3] R. Byrd, J. Nocedal and R. Schnabel "Representations of
        //     Quasi-Newton Matrices and their use in Limited Memory Methods'',
        //     Mathematical Programming 63 (1994), no. 4, pp. 129-156.
        //
        //     (Postscript files of these papers are available via anonymous
        //      ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        //   ************
        //function
        void mainlb(const int& n, const int& m, double* const & x, 
                const double* const & l, const double* const & u, 
                const int* const & nbd, double& f,  double* const & g, 
                const double& factr, const double& pgtol, 
                double* const & ws, double* const & wy, double* const & sy,
                double* const & ss, double* const & yy, double* const & wt,
                double* const & wn, double* const & snd, double* const & z,
                double* const & r, double* const & d, double* const & t, 
                double* const & wa, double* const & sg, double* const& sgo,
                double* const & yg, double* const & ygo, int* const& index,
                int* const & iwhere, int* const & indx2, char* const& task,
                const int& iprint, char* const& csave, bool * const& lsave,
                int* const & isave, double* const & dsave, ofstream* itfile)
        {
                bool   prjctd,cnstnd,boxed,updatd,wrk;
                char   word[4];
                int    k,nintol,iback,nskip,
                        head,col,iter,itail,iupdat,
                        nint,nfgv,info,ifun=0,
                        iword,nfree,nact,ileave,nenter;
                double theta,fold=0,dr,rr,tol,
                        xstep,sbgnrm,ddum,dnorm=0,dtd,epsmch,
                        cpu1=0,cpu2,cachyt,sbtime,lnscht,time1,time2,
                        gd,gdold=0,stp,stpmx=0,time;
                double one=1.0,zero=0.0;

                if (iprint > 1) { assert(itfile); }

                if (strncmp(task,"START",5)==0)
                {
                        time1 = timer_.time();

                        //Generate the current machine precision.
                        epsmch = dpmeps();
                        //epsmch = 1.08420217E-19;
                        //cout << "L-BFGS-B computed machine precision = " << epsmch << endl;

                        //Initialize counters and scalars when task='START'.

                        //for the limited memory BFGS matrices:
                        col    = 0;
                        head   = 1;
                        theta  = one;
                        iupdat = 0;
                        updatd = false;

                        //for operation counts:
                        iter   = 0;
                        nfgv   = 0;
                        nint   = 0;
                        nintol = 0;
                        nskip  = 0;
                        nfree  = n;

                        //for stopping tolerance:
                        tol = factr*epsmch;

                        //for measuring running time:
                        cachyt = 0;
                        sbtime = 0;
                        lnscht = 0;


                        //'word' records the status of subspace solutions.
                        strcpy(word, "---");

                        //'info' records the termination information.
                        info = 0;

                        //commented out: file opened at beginning of function
                        //if (iprint >= 1)
                        //{
                        //  //open a summary file 'iterate.dat'
                        //  ofstream itfile("lbfgsb-iterate.dat");
                        //}

                        //Check the input arguments for errors.
                        errclb(n,m,factr,l,u,nbd,task,info,k);
                        if (strncmp(task,"ERROR",5)==0)
                        {
                                prn3lb(n,x,f,task,iprint,info,itfile,
                                        iter,nfgv,nintol,nskip,nact,sbgnrm,
                                        zero,nint,word,iback,stp,xstep,k,
                                        cachyt,sbtime,lnscht);
                                return;
                        }

                        prn1lb(n,m,l,u,x,iprint,itfile,epsmch);

                        //Initialize iwhere & project x onto the feasible set.
                        active(n,l,u,nbd,x,iwhere,iprint,prjctd,cnstnd,boxed) ;

                        //The end of the initialization.
                } //if (strncmp(task,"START",5)==0)
                else
                {
                        //restore local variables.
                        prjctd = lsave[1];
                        cnstnd = lsave[2];
                        boxed  = lsave[3];
                        updatd = lsave[4];

                        nintol = isave[1];
                        //itfile = isave[3];
                        iback  = isave[4];
                        nskip  = isave[5];
                        head   = isave[6];
                        col    = isave[7];
                        itail  = isave[8];
                        iter   = isave[9];
                        iupdat = isave[10];
                        nint   = isave[12];
                        nfgv   = isave[13];
                        info   = isave[14];
                        ifun   = isave[15];
                        iword  = isave[16];
                        nfree  = isave[17];
                        nact   = isave[18];
                        ileave = isave[19];
                        nenter = isave[20];

                        theta  = dsave[1];
                        fold   = dsave[2];
                        tol    = dsave[3];
                        dnorm  = dsave[4];
                        epsmch = dsave[5];
                        cpu1   = dsave[6];
                        cachyt = dsave[7];
                        sbtime = dsave[8];
                        lnscht = dsave[9];
                        time1  = dsave[10];
                        gd     = dsave[11];
                        stpmx  = dsave[12];
                        sbgnrm = dsave[13];
                        stp    = dsave[14];
                        gdold  = dsave[15];
                        dtd    = dsave[16];


                        //After returning from the driver go to the point where execution
                        //is to resume.
                        if (strncmp(task,"FG_LN",5)==0) goto goto66;
                        if (strncmp(task,"NEW_X",5)==0) goto goto777;
                        if (strncmp(task,"FG_ST",5)==0) goto goto111;
                        if (strncmp(task,"STOP",4)==0)  
                        {
                                if (strncmp(&(task[6]),"CPU",3)==0)
                                {
                                        //restore the previous iterate.
                                        dcopy(n,t,1,x,1);
                                        dcopy(n,r,1,g,1);
                                        f = fold;
                                }
                                goto goto999;
                        }
                }

                //Compute f0 and g0.

                strcpy(task, "FG_START");
                //return to the driver to calculate f and g; reenter at goto111.
                goto goto1000;

goto111:

                nfgv = 1;

                //Compute the infinity norm of the (-) projected gradient.

                projgr(n,l,u,nbd,x,g,sbgnrm);

                if (iprint >= 1)
                {
                        cout << "At iterate " << iter << ": f=" <<f<<",  |proj g|="<<sbgnrm<<endl;
                        *itfile << "At iterate " << iter << ": nfgv=" << nfgv 
                                << ", sbgnrm=" << sbgnrm << ", f=" << f << endl;
                }

                if (sbgnrm <= pgtol) 
                {
                        //terminate the algorithm.
                        strcpy(task,"CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL");
                        goto goto999;
                }

                //----------------- the beginning of the loop ---------------------

goto222:
                if (iprint >= 99) cout << "ITERATION " << (iter + 1) << endl;
                iword = -1;

                if (!cnstnd && col > 0) 
                {
                        //skip the search for GCP.
                        dcopy(n,x,1,z,1);
                        wrk = updatd;
                        nint = 0;
                        goto goto333;
                }

                /*********************************************************************
                *
                *     Compute the Generalized Cauchy Point (GCP).
                *
                *********************************************************************/

                cpu1 = timer_.time();

                cauchy(n,x,l,u,nbd,g,indx2,iwhere,t,d,z,
                        m,wy,ws,sy,wt,theta,col,head,
                        &(wa[1-1]),&(wa[2*m+1-1]),&(wa[4*m+1-1]),&(wa[6*m+1-1]),nint,
                        sg,yg,iprint,sbgnrm,info,epsmch);

                if (info > 0)
                { 
                        //singular triangular system detected; refresh the lbfgs memory.
                        if (iprint >= 1) 
                                cout << "Singular triangular system detected; " 
                                << "refresh the lbfgs memory and restart the iteration." <<endl;
                        info   = 0;
                        col    = 0;
                        head   = 1;
                        theta  = one;
                        iupdat = 0;
                        updatd = false;
                        cpu2 = timer_.time();
                        cachyt = cachyt + cpu2 - cpu1;
                        goto goto222;
                }
                cpu2 = timer_.time(); 
                cachyt = cachyt + cpu2 - cpu1;
                nintol = nintol + nint;

                //Count the entering and leaving variables for iter > 0; 
                //find the index set of free and active variables at the GCP.

                freev(n,nfree,index,nenter,ileave,indx2,
                        iwhere,wrk,updatd,cnstnd,iprint,iter);

                nact = n - nfree;

goto333:

                //If there are no free variables or B=theta*I, then
                //skip the subspace minimization.

                if (nfree == 0 || col == 0) goto goto555;

                /**********************************************************************
                *
                *     Subspace minimization.
                *
                **********************************************************************/

                cpu1 = timer_.time();

                //Form  the LEL^T factorization of the indefinite
                //matrix    K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
                //              [L_a -R_z           theta*S'AA'S ]
                //where     E = [-I  0]
                //              [ 0  I]

                if (wrk) formk(n,nfree,index,nenter,ileave,indx2,iupdat,
                        updatd,wn,snd,m,ws,wy,sy,theta,col,head,info);

                if (info != 0) 
                {
                        //nonpositive definiteness in Cholesky factorization;
                        //refresh the lbfgs memory and restart the iteration.
                        if(iprint >= 1) 
                                cout << "Nonpositive definiteness in Cholesky factorization in formk;"
                                << "refresh the lbfgs memory and restart the iteration." << endl;
                        info   = 0;
                        col    = 0;
                        head   = 1;
                        theta  = one;
                        iupdat = 0;
                        updatd = false;
                        cpu2 = timer_.time() ;
                        sbtime = sbtime + cpu2 - cpu1;
                        goto goto222;
                }

                //compute r=-Z'B(xcp-xk)-Z'g (using wa(2m+1)=W'(xcp-x)
                //from 'cauchy').
                cmprlb(n,m,x,g,ws,wy,sy,wt,z,r,wa,index,
                        theta,col,head,nfree,cnstnd,info);

                if (info != 0) goto goto444;
                //call the direct method.
                subsm(n,m,nfree,index,l,u,nbd,z,r,ws,wy,theta,
                        col,head,iword,wa,wn,iprint,info);


goto444:
                if (info != 0) 
                {
                        //singular triangular system detected;
                        //refresh the lbfgs memory and restart the iteration.
                        if(iprint >= 1)         
                                cout << "Singular triangular system detected; " 
                                << "refresh the lbfgs memory and restart the iteration." <<endl;
                        info   = 0;
                        col    = 0;
                        head   = 1;
                        theta  = one;
                        iupdat = 0;
                        updatd = false;
                        cpu2   = timer_.time(); 
                        sbtime = sbtime + cpu2 - cpu1;
                        goto goto222;
                }

                cpu2 = timer_.time(); 
                sbtime = sbtime + cpu2 - cpu1 ;
goto555:  

                /*********************************************************************
                *
                *     Line search and optimality tests.
                *
                *********************************************************************/

                //Generate the search direction d:=z-x.

                for (int i = 1; i <= n; i++)
                        d[i] = z[i] - x[i];

                cpu1 = timer_.time();
goto66: 

                lnsrlb(n,l,u,nbd,x,f,fold,gd,gdold,g,d,r,t,z,stp,dnorm,
                        dtd,xstep,stpmx,iter,ifun,iback,nfgv,info,task,
                        boxed,cnstnd,csave,&(isave[22-1]),&(dsave[17-1]));

                if (info != 0 || iback >= 20)
                {
                        //restore the previous iterate.
                        dcopy(n,t,1,x,1);
                        dcopy(n,r,1,g,1);
                        f = fold;
                        if (col == 0) 
                        {
                                //abnormal termination.
                                if (info == 0) 
                                {
                                        info = -9;
                                        //restore the actual number of f and g evaluations etc.
                                        nfgv  -=  1;
                                        ifun  -=  1;
                                        iback -= 1;
                                }
                                strcpy(task, "ABNORMAL_TERMINATION_IN_LNSRCH");
                                iter += 1;
                                goto goto999;
                        }
                        else
                        {
                                //refresh the lbfgs memory and restart the iteration.
                                if(iprint >= 1)
                                        cout << "Bad direction in the line search; " 
                                        << "the lbfgs memory and restart the iteration" << endl;
                                if (info == 0) nfgv = nfgv - 1;
                                info   = 0;
                                col    = 0;
                                head   = 1;
                                theta  = one;
                                iupdat = 0;
                                updatd = false;
                                strcpy(task, "RESTART_FROM_LNSRCH");
                                cpu2 = timer_.time();
                                lnscht += cpu2 - cpu1;
                                goto goto222;
                        }
                }
                else 
                        if (strncmp(task,"FG_LN",5)==0) 
                        {
                                //return to the driver for calculating f and g; reenter at goto66.
                                goto goto1000;
                        }
                        else 
                        {
                                //calculate and print out the quantities related to the new X.
                                cpu2 = timer_.time() ;
                                lnscht += cpu2 - cpu1;
                                iter += 1;

                                //Compute the infinity norm of the projected (-)gradient.

                                projgr(n,l,u,nbd,x,g,sbgnrm);

                                //Print iteration information.

                                prn2lb(n,x,f,g,iprint,itfile,iter,nfgv,nact,
                                        sbgnrm,nint,word,iword,iback,stp,xstep);
                                goto goto1000;
                        }
goto777:

                        //Test for termination.

                        if (sbgnrm <= pgtol) 
                        {
                                //terminate the algorithm.
                                strcpy(task, "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL");
                                goto goto999;
                        } 

                        ddum = max(fabs(fold), fabs(f), one);
                        if ((fold - f) <= tol*ddum) 
                        {
                                //terminate the algorithm.
                                strcpy(task, "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH");
                                if (iback >= 10) info = -5;
                                //i.e., to issue a warning if iback>10 in the line search.
                                goto goto999;
                        } 

                        //Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's.

                        for (int i = 1; i <= n; i++)
                                r[i] = g[i] - r[i];
                        rr = ddot(n,r,1,r,1);
                        if (stp == one) 
                        {  
                                dr = gd - gdold;
                                ddum = -gdold;
                        }
                        else
                        {
                                dr = (gd - gdold)*stp;
                                dscal(n,stp,d,1);
                                ddum = -gdold*stp;
                        }

                        if (dr <= epsmch*ddum) 
                        {
                                //skip the L-BFGS update.
                                nskip = nskip + 1;
                                updatd = false;
                                if (iprint >= 1) 
                                        cout << "  ys=" << dr << "   -gs=" <<ddum<<" BFSG update SKIPPED"<<endl;
                                goto goto888;
                        } 

                        /*********************************************************************
                        *
                        *     Update the L-BFGS matrix.
                        *
                        *********************************************************************/

                        updatd = true;
                        iupdat += 1;

                        //Update matrices WS and WY and form the middle matrix in B.

                        matupd(n,m,ws,wy,sy,ss,d,r,itail,
                                iupdat,col,head,theta,rr,dr,stp,dtd);

                        //Form the upper half of the pds T = theta*SS + L*D^(-1)*L';
                        //Store T in the upper triangular of the array wt;
                        //Cholesky factorize T to J*J' with
                        //J' stored in the upper triangular of wt.

                        formt(m,wt,sy,ss,col,theta,info);

                        if (info != 0) 
                        {
                                //nonpositive definiteness in Cholesky factorization;
                                //refresh the lbfgs memory and restart the iteration.
                                if(iprint >= 1) 
                                        cout << "Nonpositive definiteness in Cholesky factorization in formt; "
                                        << "refresh the lbfgs memory and restart the iteration." << endl;
                                info = 0;
                                col  = 0;
                                head = 1;
                                theta = one;
                                iupdat = 0;
                                updatd = false;
                                goto goto222;
                        }

                        //Now the inverse of the middle matrix in B is

                        //[  D^(1/2)      O ] [ -D^(1/2)  D^(-1/2)*L' ]
                        //[ -L*D^(-1/2)   J ] [  0        J'          ]

goto888:

                        //-------------------- the end of the loop -----------------------------

                        goto goto222;

goto999:
                        time2 = timer_.time();
                        time = time2 - time1;
                        prn3lb(n,x,f,task,iprint,info,itfile,
                                iter,nfgv,nintol,nskip,nact,sbgnrm,
                                time,nint,word,iback,stp,xstep,k,
                                cachyt,sbtime,lnscht);
goto1000:

                        //Save local variables.
                        lsave[1]  = prjctd;
                        lsave[2]  = cnstnd;
                        lsave[3]  = boxed;
                        lsave[4]  = updatd;

                        isave[1]  = nintol; 
                        //isave[3]  = itfile;
                        isave[4]  = iback;
                        isave[5]  = nskip; 
                        isave[6]  = head;
                        isave[7]  = col;
                        isave[8]  = itail;
                        isave[9]  = iter; 
                        isave[10] = iupdat; 
                        isave[12] = nint; 
                        isave[13] = nfgv; 
                        isave[14] = info; 
                        isave[15] = ifun; 
                        isave[16] = iword; 
                        isave[17] = nfree; 
                        isave[18] = nact; 
                        isave[19] = ileave; 
                        isave[20] = nenter; 

                        dsave[1]  = theta; 
                        dsave[2]  = fold; 
                        dsave[3]  = tol; 
                        dsave[4]  = dnorm; 
                        dsave[5]  = epsmch; 
                        dsave[6]  = cpu1; 
                        dsave[7]  = cachyt; 
                        dsave[8]  = sbtime; 
                        dsave[9]  = lnscht; 
                        dsave[10] = time1;
                        dsave[11] = gd;
                        dsave[12] = stpmx; 
                        dsave[13] = sbgnrm;
                        dsave[14] = stp;
                        dsave[15] = gdold;
                        dsave[16] = dtd;

        }//mainlb()


        //   ************
        //
        //   Subroutine active
        //
        //   This subroutine initializes iwhere and projects the initial x to
        //     the feasible set if necessary.
        //
        //   iwhere is an integer array of dimension n.
        //     On entry iwhere is unspecified.
        //     On exit iwhere(i)=-1  if x(i) has no bounds
        //                       3   if l(i)=u(i)
        //                       0   otherwise.
        //     In cauchy, iwhere is given finer gradations.
        //
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        //   ************
        //function
        void active(const int& n, const double* const & l, 
                const double* const & u, const int* const & nbd, 
                double* const & x, int* const & iwhere, const int& iprint,
                bool& prjctd, bool& cnstnd, bool& boxed)
        {
                //int i;
                int nbdd;
                double zero=0.0;

                //Initialize nbdd, prjctd, cnstnd and boxed.
                nbdd   = 0;
                prjctd = false;
                cnstnd = false;
                boxed  = true;

                //Project the initial x to the easible set if necessary.

                for (int i = 1; i <= n; i++)
                {
                        if (nbd[i] > 0)
                        {
                                if (nbd[i] <= 2 && x[i] <= l[i])
                                {
                                        if (x[i] < l[i])
                                        {
                                                prjctd = true;
                                                x[i] = l[i];
                                        }
                                        nbdd += 1;
                                }
                                else 
                                        if (nbd[i] >= 2 && x[i] >= u[i])
                                        {
                                                if (x[i] > u[i])
                                                {
                                                        prjctd = true;
                                                        x[i] = u[i];
                                                }
                                                nbdd += 1;
                                        }
                        }
                }

                //Initialize iwhere and assign values to cnstnd and boxed.

                for (int i = 1; i <= n; i++)
                {
                        if (nbd[i] != 2) boxed = false;
                        if (nbd[i] == 0)
                        {
                                //this variable is always free
                                iwhere[i] = -1;
                                //otherwise set x(i)=mid(x(i), u(i), l(i)).
                        }
                        else
                        {
                                cnstnd = true;
                                if (nbd[i] == 2 && (u[i]-l[i]) <= zero)
                                {
                                        //this variable is always fixed
                                        iwhere[i] = 3;
                                }
                                else 
                                        iwhere[i] = 0;
                        }
                }

                if (iprint >= 0)
                {
                        if (prjctd) 
                                cout << "The initial X is infeasible.  Restart with its projection."
                                << endl;
                        if (!cnstnd) cout << "This problem is unconstrained." << endl;
                }

                if (iprint > 0) 
                        cout << "At X0 " << nbdd << " variables are exactly at the bounds" <<endl;
        }//active()


        //   ************
        //
        //   Subroutine bmv
        //
        //   This subroutine computes the product of the 2m x 2m middle matrix 
        //       in the compact L-BFGS formula of B and a 2m vector v;  
        //       it returns the product in p.
        //      
        //   m is an integer variable.
        //     On entry m is the maximum number of variable metric corrections
        //       used to define the limited memory matrix.
        //     On exit m is unchanged.
        //
        //   sy is a double precision array of dimension m x m.
        //     On entry sy specifies the matrix S'Y.
        //     On exit sy is unchanged.
        //
        //   wt is a double precision array of dimension m x m.
        //     On entry wt specifies the upper triangular matrix J' which is 
        //       the Cholesky factor of (thetaS'S+LD^(-1)L').
        //     On exit wt is unchanged.
        //
        //   col is an integer variable.
        //     On entry col specifies the number of s-vectors (or y-vectors)
        //       stored in the compact L-BFGS formula.
        //     On exit col is unchanged.
        //
        //   v is a double precision array of dimension 2col.
        //     On entry v specifies vector v.
        //     On exit v is unchanged.
        //
        //   p is a double precision array of dimension 2col.
        //     On entry p is unspecified.
        //     On exit p is the product Mv.
        //
        //   info is an integer variable.
        //     On entry info is unspecified.
        //     On exit info = 0       for normal return,
        //                  = nonzero for abnormal return when the system
        //                              to be solved by dtrsl is singular.
        //
        //   Subprograms called:
        //
        //     Linpack ... dtrsl.
        //
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        //   ************
        //function
        void bmv(const int& m, double* const & sy, double* const & wt, 
                const int& col, const double* const & v, double* const & p, 
                int& info)
        {
                //int i,k;
                int i2;
                double sum;

                if (col == 0) return;

                //PART I: solve [  D^(1/2)      O ] [ p1 ] = [ v1 ]
                //              [ -L*D^(-1/2)   J ] [ p2 ]   [ v2 ].

                //solve Jp2=v2+LD^(-1)v1.
                p[col+1] = v[col+1];

                for (int i = 2; i <= col; i++)
                {
                        i2 = col + i;
                        sum = 0;
                        for (int k = 1; k <= i-1; k++)
                                sum += sy[getIdx(i,k,m)]*v[k]/sy[getIdx(k,k,m)];

                        p[i2] = v[i2] + sum;
                }

                //Solve the triangular system
                dtrsl(wt,m,col,&(p[col+1-1]),11,info);
                if (info != 0) return; 

                //solve D^(1/2)p1=v1.
                for (int i = 1; i <= col; i++)
                        p[i] = v[i]/sqrt(sy[getIdx(i,i,m)]);

                //PART II: solve [ -D^(1/2)   D^(-1/2)*L'  ] [ p1 ] = [ p1 ]
                //               [  0         J'           ] [ p2 ]   [ p2 ]. 

                //solve J^Tp2=p2. 
                dtrsl(wt,m,col,&(p[col+1-1]),1,info);
                if (info != 0) return;

                //compute p1=-D^(-1/2)(p1-D^(-1/2)L'p2)
                //          =-D^(-1/2)p1+D^(-1)L'p2.  
                for (int i = 1; i <= col; i++)
                        p[i] = -p[i]/sqrt(sy[getIdx(i,i,m)]);

                for (int i = 1; i <= col; i++)
                {
                        sum = 0;
                        for (int k = i+1; i <= col; i++)
                                sum += sy[getIdx(k,i,m)]*p[col+k]/sy[getIdx(i,i,m)];
                        p[i] += sum;
                }
        }//bmv()


        //   ************
        //
        //   Subroutine cauchy
        //
        //   For given x, l, u, g (with sbgnrm > 0), and a limited memory
        //     BFGS matrix B defined in terms of matrices WY, WS, WT, and
        //     scalars head, col, and theta, this subroutine computes the
        //     generalized Cauchy point (GCP), defined as the first local
        //     minimizer of the quadrati//
        //
        //                Q(x + s) = g's + 1/2 s'Bs
        //
        //     along the projected gradient direction P(x-tg,l,u).
        //     The routine returns the GCP in xcp. 
        //     
        //   n is an integer variable.
        //     On entry n is the dimension of the problem.
        //     On exit n is unchanged.
        //
        //   x is a double precision array of dimension n.
        //     On entry x is the starting point for the GCP computation.
        //     On exit x is unchanged.
        //
        //   l is a double precision array of dimension n.
        //     On entry l is the lower bound of x.
        //     On exit l is unchanged.
        //
        //   u is a double precision array of dimension n.
        //     On entry u is the upper bound of x.
        //     On exit u is unchanged.
        //
        //   nbd is an integer array of dimension n.
        //     On entry nbd represents the type of bounds imposed on the
        //       variables, and must be specified as follows:
        //       nbd(i)=0 if x(i) is unbounded,
        //              1 if x(i) has only a lower bound,
        //              2 if x(i) has both lower and upper bounds, and
        //              3 if x(i) has only an upper bound. 
        //     On exit nbd is unchanged.
        //
        //   g is a double precision array of dimension n.
        //     On entry g is the gradient of f(x).  g must be a nonzero vector.
        //     On exit g is unchanged.
        //
        //   iorder is an integer working array of dimension n.
        //     iorder will be used to store the breakpoints in the piecewise
        //     linear path and free variables encountered. On exit,
        //       iorder(1),...,iorder(nleft) are indices of breakpoints
        //                              which have not been encountered; 
        //       iorder(nleft+1),...,iorder(nbreak) are indices of
        //                                   encountered breakpoints; and
        //       iorder(nfree),...,iorder(n) are indices of variables which
        //               have no bound constraits along the search direction.
        //
        //   iwhere is an integer array of dimension n.
        //     On entry iwhere indicates only the permanently fixed (iwhere=3)
        //     or free (iwhere= -1) components of x.
        //     On exit iwhere records the status of the current x variables.
        //     iwhere(i)=-3  if x(i) is free and has bounds, but is not moved
        //               0   if x(i) is free and has bounds, and is moved
        //               1   if x(i) is fixed at l(i), and l(i) != u(i)
        //               2   if x(i) is fixed at u(i), and u(i) != l(i)
        //               3   if x(i) is always fixed, i.e.,  u(i)=x(i)=l(i)
        //               -1  if x(i) is always free, i.e., it has no bounds.
        //
        //   t is a double precision working array of dimension n. 
        //     t will be used to store the break points.
        //
        //   d is a double precision array of dimension n used to store
        //     the Cauchy direction P(x-tg)-x.
        //
        //   xcp is a double precision array of dimension n used to return the
        //     GCP on exit.
        //
        //   m is an integer variable.
        //     On entry m is the maximum number of variable metric corrections 
        //       used to define the limited memory matrix.
        //     On exit m is unchanged.
        //
        //   ws, wy, sy, and wt are double precision arrays.
        //     On entry they store information that defines the
        //                           limited memory BFGS matrix:
        //       ws(n,m) stores S, a set of s-vectors;
        //       wy(n,m) stores Y, a set of y-vectors;
        //       sy(m,m) stores S'Y;
        //       wt(m,m) stores the
        //               Cholesky factorization of (theta*S'S+LD^(-1)L').
        //     On exit these arrays are unchanged.
        //
        //   theta is a double precision variable.
        //     On entry theta is the scaling factor specifying B_0 = theta I.
        //     On exit theta is unchanged.
        //
        //   col is an integer variable.
        //     On entry col is the actual number of variable metri//
        //       corrections stored so far.
        //     On exit col is unchanged.
        //
        //   head is an integer variable.
        //     On entry head is the location of the first s-vector (or y-vector)
        //       in S (or Y).
        //     On exit col is unchanged.
        //
        //   p is a double precision working array of dimension 2m.
        //     p will be used to store the vector p = W^(T)d.
        //
        //   c is a double precision working array of dimension 2m.
        //     c will be used to store the vector c = W^(T)(xcp-x).
        //
        //   wbp is a double precision working array of dimension 2m.
        //     wbp will be used to store the row of W corresponding
        //       to a breakpoint.
        //
        //   v is a double precision working array of dimension 2m.
        //
        //   nint is an integer variable.
        //     On exit nint records the number of quadratic segments explored
        //       in searching for the GCP.
        //
        //   sg and yg are double precision arrays of dimension m.
        //     On entry sg  and yg store S'g and Y'g correspondingly.
        //     On exit they are unchanged. 
        // 
        //   iprint is an INTEGER variable that must be set by the user.
        //     It controls the frequency and type of output generated:
        //      iprint<0    no output is generated;
        //      iprint=0    print only one line at the last iteration;
        //      0<iprint<99 print also f and |proj g| every iprint iterations;
        //      iprint=99   print details of every iteration except n-vectors;
        //      iprint=100  print also the changes of active set and final x;
        //      iprint>100  print details of every iteration including x and g;
        //     When iprint > 0, the file iterate.dat will be created to
        //                      summarize the iteration.
        //
        //   sbgnrm is a double precision variable.
        //     On entry sbgnrm is the norm of the projected gradient at x.
        //     On exit sbgnrm is unchanged.
        //
        //   info is an integer variable.
        //     On entry info is 0.
        //     On exit info = 0       for normal return,
        //                  = nonzero for abnormal return when the the system
        //                            used in routine bmv is singular.
        //
        //   Subprograms called:
        // 
        //     L-BFGS-B Library ... hpsolb, bmv.
        //
        //     Linpack ... dscal dcopy, daxpy.
        //
        //
        //   References:
        //
        //     [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
        //     memory algorithm for bound constrained optimization'',
        //     SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
        //
        //     [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN
        //     Subroutines for Large Scale Bound Constrained Optimization''
        //     Tech. Report, NAM-11, EECS Department, Northwestern University,
        //     1994.
        //
        //     (Postscript files of these papers are available via anonymous
        //      ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        //   ************  
        //function
        void cauchy(const int& n, const double* const & x, 
                const double* const & l, const double* const & u, 
                const int* const & nbd, const double* const & g, 
                int* const & iorder, int* const & iwhere, 
                double* const & t, double* const & d, double* const & xcp, 
                const int& m, double* const & wy, 
                double* const & ws, double* const & sy, 
                double* const & wt, double& theta, 
                const int& col, const int& head, double* const & p, 
                double* const & c, double* const & wbp, 
                double* const & v, int& nint, const double* const & sg, 
                const double* const & yg, const int& iprint, 
                const double& sbgnrm, int& info, const double& epsmch)
        {
                bool xlower,xupper,bnded;
                //int i,j;
                int col2,nfree,nbreak,pointr,
                        ibp,nleft,ibkmin,iter;
                double f1,f2,dt,dtm,tsum,dibp,zibp,dibp2,bkmin,
                        tu=0,tl=0,wmc,wmp,wmw,tj,tj0,neggi,
                        f2_org;
                double one=1,zero=0;

                //Check the status of the variables, reset iwhere(i) if necessary;
                //compute the Cauchy direction d and the breakpoints t; initialize
                //the derivative f1 and the vector p = W'd (for theta = 1).

                if (sbgnrm <= zero)
                {
                        if (iprint >= 0) cout << "Subgnorm = 0.  GCP = X." << endl;
                        dcopy(n,x,1,xcp,1);
                        return;
                }
                bnded = true;
                nfree = n + 1;
                nbreak = 0;
                ibkmin = 0;
                bkmin = zero;
                col2 = 2*col;
                f1 = zero;
                if (iprint >= 99) 
                        cout << "---------------- CAUCHY entered-------------------" << endl;

                //We set p to zero and build it up as we determine d.
                for (int i = 1; i <= col2; i++)
                        p[i] = zero;

                //In the following loop we determine for each variable its bound
                //status and its breakpoint, and update p accordingly.
                //Smallest breakpoint is identified.

                for (int i = 1; i <= n; i++)      
                {
                        neggi = -g[i];  
                        if (iwhere[i] != 3 && iwhere[i] != -1)
                        {
                                //if x(i) is not a constant and has bounds,
                                //compute the difference between x(i) and its bounds.
                                if (nbd[i] <= 2) tl = x[i] - l[i];
                                if (nbd[i] >= 2) tu = u[i] - x[i];

                                //If a variable is close enough to a bound
                                //we treat it as at bound.
                                xlower = (nbd[i] <= 2 && tl <= zero);
                                xupper = (nbd[i] >= 2 && tu <= zero);

                                //reset iwhere(i).
                                iwhere[i] = 0;
                                if (xlower)
                                { 
                                        if (neggi <= zero) iwhere[i] = 1;
                                }
                                else 
                                        if (xupper)
                                        {
                                                if (neggi >= zero) iwhere[i] = 2;
                                        }
                                        else
                                        {
                                                if (fabs(neggi) <= zero) iwhere[i] = -3;
                                        }
                        }
                        pointr = head;
                        if (iwhere[i] != 0 && iwhere[i] != -1)
                                d[i] = zero;
                        else
                        {
                                d[i] = neggi;
                                f1 -= neggi*neggi;
                                //calculate p := p - W'e_i* (g_i).
                                for (int j = 1; j <= col; j++)
                                {
                                        p[j] += wy[getIdx(i,pointr,n)] * neggi;
                                        p[col + j] += ws[getIdx(i,pointr,n)] * neggi;
                                        pointr = pointr%m + 1;
                                }
                                if (nbd[i] <= 2 && nbd[i] != 0 && neggi < zero)
                                {
                                        //x[i] + d[i] is bounded; compute t[i].
                                        nbreak += 1;
                                        iorder[nbreak] = i;
                                        t[nbreak] = tl/(-neggi);
                                        if (nbreak == 1 || t[nbreak] < bkmin)
                                        {
                                                bkmin = t[nbreak];
                                                ibkmin = nbreak;
                                        }
                                }
                                else
                                        if (nbd[i] >= 2 && neggi > zero)
                                        {
                                                //x(i) + d(i) is bounded; compute t(i).
                                                nbreak += 1;
                                                iorder[nbreak] = i;
                                                t[nbreak] = tu/neggi;
                                                if (nbreak == 1 || t[nbreak] < bkmin)
                                                {
                                                        bkmin = t[nbreak];
                                                        ibkmin = nbreak;
                                                }          
                                        }
                                        else
                                        {
                                                //x(i) + d(i) is not bounded.
                                                nfree -= 1;
                                                iorder[nfree] = i;
                                                if (fabs(neggi) > zero) bnded = false;
                                        }
                        }
                } //for (int i = 1; i <= n; i++)      


                //The indices of the nonzero components of d are now stored
                //in iorder(1),...,iorder(nbreak) and iorder(nfree),...,iorder(n).
                //The smallest of the nbreak breakpoints is in t(ibkmin)=bkmin.

                if (theta != one)
                {
                        //complete the initialization of p for theta not= one.
                        dscal(col,theta,&(p[col+1-1]),1);
                }

                //Initialize GCP xcp = x.

                dcopy(n,x,1,xcp,1);

                if (nbreak == 0 && nfree == n + 1)
                {
                        //is a zero vector, return with the initial xcp as GCP.
                        if (iprint > 100) 
                        {
                                cout << "Cauchy X = ";
                                for (int i = 1; i <= n; i++)
                                        cout << xcp[i] << " ";
                                cout << endl;
                        }
                        return;
                }    

                //Initialize c = W'(xcp - x) = 0.

                for (int j = 1; j <= col2; j++)
                        c[j] = zero;

                //Initialize derivative f2.

                f2 = -theta*f1;
                f2_org = f2;
                if (col > 0) 
                {
                        bmv(m,sy,wt,col,p,v,info);
                        if (info != 0) return;
                        f2 -= ddot(col2,v,1,p,1);
                }
                dtm = -f1/f2;
                tsum = zero;
                nint = 1;
                if (iprint >= 99)
                        cout << "There are " << nbreak << " breakpoints" << endl;

                //If there are no breakpoints, locate the GCP and return. 

                if (nbreak == 0) goto goto888;

                nleft = nbreak;
                iter = 1;

                tj = zero;

                //------------------- the beginning of the loop -------------------------

goto777:

                //Find the next smallest breakpoint;
                //compute dt = t(nleft) - t(nleft + 1).

                tj0 = tj;
                if (iter == 1) 
                {
                        //Since we already have the smallest breakpoint we need not do
                        //heapsort yet. Often only one breakpoint is used and the
                        //cost of heapsort is avoided.
                        tj = bkmin;
                        ibp = iorder[ibkmin];
                }
                else
                {
                        if (iter == 2)
                        {
                                //Replace the already used smallest breakpoint with the
                                //breakpoint numbered nbreak > nlast, before heapsort call.
                                if (ibkmin != nbreak)
                                {
                                        t[ibkmin] = t[nbreak];
                                        iorder[ibkmin] = iorder[nbreak];
                                }
                                //Update heap structure of breakpoints
                                //(if iter=2, initialize heap).
                        }

                        hpsolb(nleft,t,iorder,iter-2);
                        tj = t[nleft];
                        ibp = iorder[nleft];    
                }        

                dt = tj - tj0;

                if (dt != zero && iprint >= 100)
                {
                        cout << "Piece    " << nint << " --f1, f2 at start point " 
                                << f1 << "," << f2 << endl;
                        cout << "Distance to the next break point = " << dt << endl;
                        cout << "Distance to the stationary point = " << dtm << endl;
                }

                //If a minimizer is within this interval, locate the GCP and return. 

                if (dtm < dt) goto goto888;

                //Otherwise fix one variable and
                //reset the corresponding component of d to zero.

                tsum += dt;
                nleft -= 1;
                iter += 1;
                dibp = d[ibp];
                d[ibp] = zero;
                if (dibp > zero) 
                {
                        zibp = u[ibp] - x[ibp];
                        xcp[ibp] = u[ibp];
                        iwhere[ibp] = 2;
                }
                else
                {
                        zibp = l[ibp] - x[ibp];
                        xcp[ibp] = l[ibp];
                        iwhere[ibp] = 1;
                }
                if (iprint >= 100) cout << "Variable  " << ibp << " is fixed." << endl;
                if (nleft == 0 && nbreak == n) 
                {
                        //all n variables are fixed, return with xcp as GCP.
                        dtm = dt;
                        goto goto999;
                }

                //Update the derivative information.

                nint += 1;
                dibp2 = dibp*dibp;

                //Update f1 and f2.

                //temporarily set f1 and f2 for col=0.
                f1 += dt*f2 + dibp2 - theta*dibp*zibp;
                f2 -= theta*dibp2;

                if (col > 0)
                {
                        //update c = c + dt*p.
                        daxpy(col2,dt,p,1,c,1);

                        //choose wbp,
                        //the row of W corresponding to the breakpoint encountered.
                        pointr = head;
                        for (int j = 1; j <= col; j++)
                        {
                                wbp[j] = wy[getIdx(ibp,pointr,n)];
                                wbp[col + j] = theta*ws[getIdx(ibp,pointr,n)];
                                pointr = pointr%m + 1;
                        }
                        //compute (wbp)Mc, (wbp)Mp, and (wbp)M(wbp)'.
                        bmv(m,sy,wt,col,wbp,v,info);
                        if (info != 0) return;
                        wmc = ddot(col2,c,1,v,1);
                        wmp = ddot(col2,p,1,v,1); 
                        wmw = ddot(col2,wbp,1,v,1);

                        //update p = p - dibp*wbp. 
                        daxpy(col2,-dibp,wbp,1,p,1);

                        //complete updating f1 and f2 while col > 0.
                        f1 += dibp*wmc;
                        f2 += 2.0*dibp*wmp - dibp2*wmw;
                }

                f2 = max(epsmch*f2_org,f2);
                if (nleft > 0)
                {
                        dtm = -f1/f2;
                        goto goto777;
                        //to repeat the loop for unsearched intervals. 
                }
                else
                        if(bnded)
                        {
                                f1 = zero;
                                f2 = zero;
                                dtm = zero;
                        }
                        else
                                dtm = -f1/f2;

                //------------------- the end of the loop -------------------------------

goto888:
                if (iprint >= 99)
                {
                        cout << endl << "GCP found in this segment" << endl
                                << "Piece    " << nint << " --f1, f2 at start point " 
                                << f1 << ","  << f2 << endl
                                << "Distance to the stationary point = " << dtm << endl;
                }

                if (dtm <= zero) dtm = zero;
                tsum += dtm;

                //Move free variables (i.e., the ones w/o breakpoints) and 
                //the variables whose breakpoints haven't been reached.

                daxpy(n,tsum,d,1,xcp,1);

goto999:

                //Update c = c + dtm*p = W'(x^c - x) 
                //which will be used in computing r = Z'(B(x^c - x) + g).

                if (col > 0) daxpy(col2,dtm,p,1,c,1);
                if (iprint > 100)
                {
                        cout<< "Cauchy X = ";
                        for (int i = 1; i <=n; i++)
                                cout << xcp[i] << " ";
                        cout << endl;
                }
                if (iprint >= 99) 
                        cout << "---------------- exit CAUCHY----------------------" << endl;
        }//cauchy()


        //   ************
        //
        //   Subroutine cmprlb 
        //
        //     This subroutine computes r=-Z'B(xcp-xk)-Z'g by using 
        //       wa(2m+1)=W'(xcp-x) from subroutine cauchy.
        //
        //   Subprograms called:
        //
        //     L-BFGS-B Library ... bmv.
        //
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        //   ************
        //function
        void cmprlb(const int& n, const int& m, const double* const & x, 
                const double* const & g, double* const & ws, 
                double* const & wy, double* const & sy, 
                double* const & wt, const double* const & z, 
                double* const & r, double* const & wa, 
                const int* const & index, const double& theta, 
                const int& col, const int& head, const int& nfree, 
                const bool& cnstnd, int& info)
        {
                //int i,j;
                int k,pointr,idx;
                double a1,a2;

                if (!cnstnd && col > 0)
                {
                        for (int i = 1; i <= n; i++)
                                r[i] = -g[i];
                }
                else
                {
                        for (int i = 1; i <= nfree; i++)
                        {
                                k = index[i];
                                r[i] = -theta*(z[k] - x[k]) - g[k];
                        }
                        bmv(m,sy,wt,col,&(wa[2*m+1-1]),&(wa[1-1]),info);
                        if (info != 0)
                        {
                                info = -8;
                                return;
                        }
                        pointr = head;
                        for (int j = 1; j <= col; j++)
                        {
                                a1 = wa[j];
                                a2 = theta*wa[col + j];
                                for (int i = 1; i <= nfree; i++)
                                {
                                        k = index[i];
                                        idx = getIdx(k,pointr,n);
                                        r[i] += wy[idx]*a1 + ws[idx]*a2;
                                }

                                pointr = pointr%m + 1;
                        }
                }
        }//cmprlb()


        //   ************
        //
        //   Subroutine errclb
        //
        //   This subroutine checks the validity of the input data.
        //
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        //   ************
        //function
        void errclb(const int& n, const int& m, const double& factr, 
                const double* const & l, const double* const & u, 
                const int* const & nbd, char* const & task, int& info, 
                int& k)
        {
                //int i;
                //double one=1.0;
                double zero=0.0;

                //Check the input arguments for errors.

                if (n <= 0) strcpy(task,"ERROR: N <= 0");
                if (m <= 0) strcpy(task,"ERROR: M <= 0");
                if (factr < zero) strcpy(task,"ERROR: FACTR < 0");

                //Check the validity of the arrays nbd[i], u[i], and l[i].
                for (int i = 1; i <= n; i++)
                {
                        if (nbd[i] < 0 || nbd[i] > 3)
                        {
                                //return
                                strcpy(task,"ERROR: INVALID NBD");
                                info = -6;
                                k = i;
                        }
                        if (nbd[i] == 2)
                        {
                                if (l[i] > u[i])
                                {
                                        //return
                                        strcpy(task, "ERROR: NO FEASIBLE SOLUTION");
                                        info = -7;
                                        k = i;
                                }
                        }
                }
        }//errclb


        //   ************
        //
        //   Subroutine formk 
        //
        //   This subroutine forms  the LEL^T factorization of the indefinite
        //
        //     matrix    K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
        //                   [L_a -R_z           theta*S'AA'S ]
        //                                                  where E = [-I  0]
        //                                                            [ 0  I]
        //   The matrix K can be shown to be equal to the matrix M^[-1]N
        //     occurring in section 5.1 of [1], as well as to the matrix
        //     Mbar^[-1] Nbar in section 5.3.
        //
        //   n is an integer variable.
        //     On entry n is the dimension of the problem.
        //     On exit n is unchanged.
        //
        //   nsub is an integer variable
        //     On entry nsub is the number of subspace variables in free set.
        //     On exit nsub is not changed.
        //
        //   ind is an integer array of dimension nsub.
        //     On entry ind specifies the indices of subspace variables.
        //     On exit ind is unchanged. 
        //
        //   nenter is an integer variable.
        //     On entry nenter is the number of variables entering the 
        //       free set.
        //     On exit nenter is unchanged. 
        //
        //   ileave is an integer variable.
        //     On entry indx2(ileave),...,indx2(n) are the variables leaving
        //       the free set.
        //     On exit ileave is unchanged. 
        //
        //   indx2 is an integer array of dimension n.
        //     On entry indx2(1),...,indx2(nenter) are the variables entering
        //       the free set, while indx2(ileave),...,indx2(n) are the
        //       variables leaving the free set.
        //     On exit indx2 is unchanged. 
        //
        //   iupdat is an integer variable.
        //     On entry iupdat is the total number of BFGS updates made so far.
        //     On exit iupdat is unchanged. 
        //
        //   updatd is a logical variable.
        //     On entry 'updatd' is true if the L-BFGS matrix is updatd.
        //     On exit 'updatd' is unchanged. 
        //
        //   wn is a double precision array of dimension 2m x 2m.
        //     On entry wn is unspecified.
        //     On exit the upper triangle of wn stores the LEL^T factorization
        //       of the 2*col x 2*col indefinite matrix
        //                   [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
        //                   [L_a -R_z           theta*S'AA'S ]
        //
        //   wn1 is a double precision array of dimension 2m x 2m.
        //     On entry wn1 stores the lower triangular part of 
        //                   [Y' ZZ'Y   L_a'+R_z']
        //                   [L_a+R_z   S'AA'S   ]
        //       in the previous iteration.
        //     On exit wn1 stores the corresponding updated matrices.
        //     The purpose of wn1 is just to store these inner products
        //     so they can be easily updated and inserted into wn.
        //
        //   m is an integer variable.
        //     On entry m is the maximum number of variable metric corrections
        //       used to define the limited memory matrix.
        //     On exit m is unchanged.
        //
        //   ws, wy, sy, and wtyy are double precision arrays;
        //   theta is a double precision variable;
        //   col is an integer variable;
        //   head is an integer variable.
        //     On entry they store the information defining the
        //                                        limited memory BFGS matrix:
        //       ws(n,m) stores S, a set of s-vectors;
        //       wy(n,m) stores Y, a set of y-vectors;
        //       sy(m,m) stores S'Y;
        //       wtyy(m,m) stores the Cholesky factorization
        //                                 of (theta*S'S+LD^(-1)L')
        //       theta is the scaling factor specifying B_0 = theta I;
        //       col is the number of variable metric corrections stored;
        //       head is the location of the 1st s- (or y-) vector in S (or Y).
        //     On exit they are unchanged.
        //
        //   info is an integer variable.
        //     On entry info is unspecified.
        //     On exit info =  0 for normal return;
        //                  = -1 when the 1st Cholesky factorization failed;
        //                  = -2 when the 2st Cholesky factorization failed.
        //
        //   Subprograms called:
        //
        //     Linpack ... dcopy, dpofa, dtrsl.
        //
        //
        //   References:
        //     [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
        //     memory algorithm for bound constrained optimization'',
        //     SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
        //
        //     [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a
        //     limited memory FORTRAN code for solving bound constrained
        //     optimization problems'', Tech. Report, NAM-11, EECS Department,
        //     Northwestern University, 1994.
        //
        //     (Postscript files of these papers are available via anonymous
        //      ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        //   ************
        //function
        void formk(const int& n, const int& nsub, const int* const & ind, 
                const int& nenter, const int& ileave, 
                const int* const & indx2, const int& iupdat, 
                const bool& updatd, double* const & wn, 
                double* const & wn1, const int& m, 
                double* const & ws, double* const & wy, 
                double* const & sy, const double& theta, 
                const int& col, const int& head, int& info)
        {
                //int i,k;
                int m2,ipntr,jpntr,iy,is,jy,js,is1,js1,k1,
                        col2,pbegin,pend,dbegin,dend,upcl;
                int idx1, idx2;
                double temp1,temp2,temp3,temp4;
                //double one=1.0;
                double zero=0.0;

                m2 = 2*m;

                //Form the lower triangular part of
                //WN1 = [Y' ZZ'Y   L_a'+R_z'] 
                //      [L_a+R_z   S'AA'S   ]
                //where L_a is the strictly lower triangular part of S'AA'Y
                //      R_z is the upper triangular part of S'ZZ'Y.

                if (updatd)
                {
                        if (iupdat > m)
                        { 
                                //shift old part of WN1.
                                for (int jy = 1; jy <= m-1; jy++)
                                { 
                                        js = m + jy;
                                        dcopy(m-jy,&(wn1[getIdx(jy+1,jy+1,m2)-1]),1,
                                                &(wn1[getIdx(jy,jy,m2)-1]),1);
                                        dcopy(m-jy,&(wn1[getIdx(js+1,js+1,m2)-1]),1,
                                                &(wn1[getIdx(js,js,m2)-1]), 1);
                                        dcopy(m-1, &(wn1[getIdx(m+2,jy+1,m2)-1]), 1,
                                                &(wn1[getIdx(m+1,jy,m2)-1]),1);
                                }
                        }

                        //put new rows in blocks (1,1), (2,1) and (2,2).
                        pbegin = 1;
                        pend = nsub;
                        dbegin = nsub + 1;
                        dend = n;
                        iy = col;
                        is = m + col;
                        ipntr = head + col - 1;
                        if (ipntr > m) ipntr -= m ;
                        jpntr = head;

                        for (jy = 1; jy <= col; jy++)
                        {
                                js = m + jy;
                                temp1 = zero;
                                temp2 = zero;
                                temp3 = zero;
                                //compute element jy of row 'col' of Y'ZZ'Y
                                for (int k = pbegin; k <= pend; k++)
                                {
                                        k1 = ind[k];
                                        temp1 += wy[getIdx(k1,ipntr,n)]*wy[getIdx(k1,jpntr,n)];
                                }
                                //compute elements jy of row 'col' of L_a and S'AA'S
                                for (int k = dbegin; k <= dend; k++)
                                {
                                        k1 = ind[k];
                                        idx1 = getIdx(k1,ipntr,n);
                                        idx2 = getIdx(k1,jpntr,n);
                                        temp2 += ws[idx1] * ws[idx2];
                                        temp3 += ws[idx1] * wy[idx2];
                                }
                                wn1[getIdx(iy,jy,m2)] = temp1;
                                wn1[getIdx(is,js,m2)] = temp2;
                                wn1[getIdx(is,jy,m2)] = temp3;
                                jpntr = jpntr%m + 1;
                        }

                        //put new column in block (2,1).
                        jy = col;
                        jpntr = head + col - 1;
                        if (jpntr > m) jpntr -= m;
                        ipntr = head;
                        for (int i = 1; i <= col; i++)
                        {
                                is = m + i;
                                temp3 = zero;
                                //compute element i of column 'col' of R_z
                                for (int k = pbegin; k <= pend; k++)
                                {
                                        k1 = ind[k];
                                        temp3 += ws[getIdx(k1,ipntr,n)]*wy[getIdx(k1,jpntr,n)];
                                }
                                ipntr = ipntr%m + 1;
                                wn1[getIdx(is,jy,m2)] = temp3;
                        }
                        upcl = col - 1;
                }
                else
                        upcl = col;

                //modify the old parts in blocks (1,1) and (2,2) due to changes
                //in the set of free variables.
                ipntr = head;

                for (int iy = 1; iy <= upcl; iy++)
                {
                        is = m + iy;
                        jpntr = head;

                        for (int jy = 1; jy <= iy; jy++)
                        {
                                js = m + jy;
                                temp1 = zero;
                                temp2 = zero;
                                temp3 = zero;
                                temp4 = zero;
                                for (int k = 1; k <= nenter; k++)
                                {
                                        k1 = indx2[k];
                                        idx1 = getIdx(k1,ipntr,n);
                                        idx2 = getIdx(k1,jpntr,n);
                                        temp1 += wy[idx1]*wy[idx2];
                                        temp2 += ws[idx1]*ws[idx2];
                                }
                                for (int k = ileave; k <= n; k++)
                                {
                                        k1 = indx2[k];
                                        idx1 = getIdx(k1,ipntr,n);
                                        idx2 = getIdx(k1,jpntr,n);
                                        temp3 += wy[idx1]*wy[idx2];
                                        temp4 += ws[idx1]*ws[idx2];
                                }
                                wn1[getIdx(iy,jy,m2)] += temp1 - temp3;
                                wn1[getIdx(is,js,m2)] += temp4 - temp2;
                                jpntr = jpntr%m + 1;
                        }
                        ipntr = ipntr%m + 1;
                } 
                //modify the old parts in block (2,1).
                ipntr = head;

                for (int is = m+1; is <= m + upcl; is++)
                {
                        jpntr = head;
                        for (int jy = 1; jy <= upcl; jy++)
                        {
                                temp1 = zero;
                                temp3 = zero;

                                for (int k = 1; k <= nenter; k++)
                                {
                                        k1 = indx2[k];
                                        temp1 += ws[getIdx(k1,ipntr,n)]*wy[getIdx(k1,jpntr,n)];
                                }
                                for (int k = ileave; k <= n; k++)
                                {
                                        k1 = indx2[k];
                                        temp3 += ws[getIdx(k1,ipntr,n)]*wy[getIdx(k1,jpntr,n)];
                                }
                                if (is <= jy + m)
                                        wn1[getIdx(is,jy,m2)] += temp1 - temp3;
                                else
                                        wn1[getIdx(is,jy,m2)] += temp3 - temp1;

                                jpntr = jpntr%m + 1;
                        }
                        ipntr = ipntr%m + 1;
                } 
                //Form the upper triangle of WN = [D+Y' ZZ'Y/theta   -L_a'+R_z' ] 
                //                                [-L_a +R_z        S'AA'S*theta]

                for (int iy = 1; iy <= col; iy++)
                {
                        is = col + iy;
                        is1 = m + iy;
                        for (int jy = 1; jy <= iy; jy++)
                        {
                                js = col + jy;
                                js1 = m + jy;
                                wn[getIdx(jy,iy,m2)] = wn1[getIdx(iy,jy,m2)]/theta;
                                wn[getIdx(js,is,m2)] = wn1[getIdx(is1,js1,m2)]*theta;
                        }

                        for (int jy = 1; jy <= iy-1; jy++)
                                wn[getIdx(jy,is,m2)] = -wn1[getIdx(is1,jy,m2)];

                        for (int jy = iy; jy <= col; jy++)
                                wn[getIdx(jy,is,m2)] = wn1[getIdx(is1,jy,m2)];
                        wn[getIdx(iy,iy,m2)] += sy[getIdx(iy,iy,m)];
                }

                //Form the upper triangle of WN= [  LL'            L^-1(-L_a'+R_z')] 
                //                               [(-L_a +R_z)L'^-1   S'AA'S*theta  ]

                //first Cholesky factor (1,1) block of wn to get LL'
                //with L' stored in the upper triangle of wn.
                dpofa(wn,m2,col,info);
                if (info != 0) 
                {
                        info = -1;
                        return;
                }

                //then form L^-1(-L_a'+R_z') in the (1,2) block.
                col2 = 2*col;
                for (int js = col+1; js <= col2; js++)
                        dtrsl(wn,m2,col,&(wn[getIdx(1,js,m2)-1]),11,info);



                //Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the
                //upper triangle of (2,2) block of wn.
                //HH                      
                for (int is = col+1; is <= col2; is++)
                        for (int js = is; js <= col2; js++)
                                wn[getIdx(is,js,m2)] += ddot(col,&(wn[getIdx(1,is,m2)-1]),1,
                                &(wn[getIdx(1,js,m2)-1]),1);

                //Cholesky factorization of (2,2) block of wn.
                dpofa(&(wn[getIdx(col+1,col+1,m2)-1]),m2,col,info);

                if (info != 0)
                {
                        info = -2;
                        return;
                }
        }//formk()


        //   ************
        //
        //   Subroutine formt
        //
        //     This subroutine forms the upper half of the pos. def. and symm.
        //       T = theta*SS + L*D^(-1)*L', stores T in the upper triangle
        //       of the array wt, and performs the Cholesky factorization of T
        //       to produce J*J', with J' stored in the upper triangle of wt.
        //
        //   Subprograms called:
        //
        //     Linpack ... dpofa.
        //
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        //   ***********
        //function
        void formt(const int& m, double* const & wt, double* const& sy,
                double* const &  ss, const int& col, const double& theta, 
                int& info)
        {
                //int i,j,k;
                int k1, idx;
                double ddum;
                double zero=0.0;

                //Form the upper half of  T = theta*SS + L*D^(-1)*L',
                //store T in the upper triangle of the array wt.
                for (int j = 1; j <= col; j++)
                {
                        idx = getIdx(1,j,m);
                        wt[idx] = theta*ss[idx];
                }

                for (int i = 2; i <= col; i++) 
                        for (int j = i; j <= col; j++) 
                        {
                                k1 = min(i,j) - 1;
                                ddum = zero;
                                for (int k = 1; k <= k1; k++) 
                                        ddum  += sy[getIdx(i,k,m)]*sy[getIdx(j,k,m)]/sy[getIdx(k,k,m)];
                                wt[getIdx(i,j,m)] = ddum + theta*ss[getIdx(i,j,m)];
                        }

                        //Cholesky factorize T to J*J' with 
                        //J' stored in the upper triangle of wt.

                        dpofa(wt,m,col,info);
                        if (info != 0) info = -3;
        }//formt()


        //   ************
        //
        //   Subroutine freev 
        //
        //   This subroutine counts the entering and leaving variables when
        //     iter > 0, and finds the index set of free and active variables
        //     at the GCP.
        //
        //   cnstnd is a logical variable indicating whether bounds are present
        //
        //   index is an integer array of dimension n
        //     for i=1,...,nfree, index(i) are the indices of free variables
        //     for i=nfree+1,...,n, index(i) are the indices of bound variables
        //     On entry after the first iteration, index gives 
        //       the free variables at the previous iteration.
        //     On exit it gives the free variables based on the determination
        //       in cauchy using the array iwhere.
        //
        //   indx2 is an integer array of dimension n
        //     On entry indx2 is unspecified.
        //     On exit with iter>0, indx2 indicates which variables
        //        have changed status since the previous iteration.
        //     For i= 1,...,nenter, indx2(i) have changed from bound to free.
        //     For i= ileave+1,...,n, indx2(i) have changed from free to bound.
        // 
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        //   ************
        //function
        void freev(const int& n, int& nfree, int* const & index, int& nenter, 
                int& ileave, int* const & indx2, int* const & iwhere, 
                bool& wrk, const bool& updatd, const bool& cnstnd, 
                const int& iprint, const int& iter)
        {
                //int i;
                int iact,k;

                nenter = 0;
                ileave = n + 1;
                if (iter > 0 && cnstnd) 
                {
                        //count the entering and leaving variables.
                        for (int i = 1; i <= nfree; i++) 
                        {
                                k = index[i];
                                if (iwhere[k] > 0)
                                {
                                        ileave -= 1;
                                        indx2[ileave] = k;
                                        if (iprint >= 100) 
                                                cout << "Variable " << k<<" leaves the set of free variables"<<endl;
                                }
                        }

                        for (int i = 1+nfree; i <= n; i++) 
                        {
                                k = index[i];
                                if (iwhere[k] <= 0)
                                {
                                        nenter += 1;
                                        indx2[nenter] = k;
                                        if (iprint >= 100) 
                                                cout << "Variable " << k<<" enters the set of free variables"<<endl;
                                }
                        }
                        if (iprint >= 99) 
                                cout << n+1-ileave << " variables leave; " 
                                << nenter << " variables enter" << endl;
                }
                wrk = (ileave < n+1) || (nenter > 0) || updatd;

                //Find the index set of free and active variables at the GCP.

                nfree = 0;
                iact = n + 1;
                for (int i = 1; i <= n; i++) 
                {
                        if (iwhere[i] <= 0)
                        {
                                nfree += 1;
                                index[nfree] = i;
                        }
                        else
                        {
                                iact -= 1;
                                index[iact] = i;
                        }
                }

                if (iprint >= 99) 
                        cout << nfree << " variables are free at GCP " << iter + 1 << endl;
        }//freev()


        //   ************
        //
        //   Subroutine hpsolb 
        //
        //   This subroutine sorts out the least element of t, and puts the
        //     remaining elements of t in a heap.
        // 
        //   n is an integer variable.
        //     On entry n is the dimension of the arrays t and iorder.
        //     On exit n is unchanged.
        //
        //   t is a double precision array of dimension n.
        //     On entry t stores the elements to be sorted,
        //     On exit t(n) stores the least elements of t, and t(1) to t(n-1)
        //       stores the remaining elements in the form of a heap.
        //
        //   iorder is an integer array of dimension n.
        //     On entry iorder(i) is the index of t(i).
        //     On exit iorder(i) is still the index of t(i), but iorder may be
        //       permuted in accordance with t.
        //
        //   iheap is an integer variable specifying the task.
        //     On entry iheap should be set as follows:
        //       iheap == 0 if t(1) to t(n) is not in the form of a heap,
        //       iheap != 0 if otherwise.
        //     On exit iheap is unchanged.
        //
        //
        //   References:
        //     Algorithm 232 of CACM (J. W. J. Williams): HEAPSORT.
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //   ************
        //function
        void hpsolb(const int& n, double* const & t, int* const & iorder, 
                const int& iheap)
        {
                //int k;
                int i,j,indxin,indxou;
                double ddum,out;

                if (iheap == 0)
                {
                        //Rearrange the elements t[1] to t[n] to form a heap.
                        for (int k = 2; k <= n; k++) 
                        {
                                ddum  = t[k];
                                indxin = iorder[k];

                                //Add ddum to the heap.
                                i = k;
goto10:
                                if (i>1)
                                {
                                        j = i/2;
                                        if (ddum < t[j])
                                        {
                                                t[i] = t[j];
                                                iorder[i] = iorder[j];
                                                i = j;
                                                goto goto10;
                                        }

                                }  
                                t[i] = ddum;
                                iorder[i] = indxin;
                        }
                }

                //Assign to 'out' the value of t(1), the least member of the heap,
                //and rearrange the remaining members to form a heap as
                //elements 1 to n-1 of t.

                if (n > 1)
                {
                        i = 1;
                        out = t[1];
                        indxou = iorder[1];
                        ddum = t[n];
                        indxin = iorder[n];

                        //Restore the heap 
goto30:
                        j = i+i;
                        if (j <= n-1) 
                        {
                                if (t[j+1] < t[j]) j = j+1;
                                if (t[j] < ddum )
                                {
                                        t[i] = t[j];
                                        iorder[i] = iorder[j];
                                        i = j;
                                        goto goto30;
                                } 
                        } 
                        t[i] = ddum;
                        iorder[i] = indxin;

                        //Put the least member in t(n). 

                        t[n] = out;
                        iorder[n] = indxou;
                }
        }//hpsolb()


        //   **********
        //
        //   Subroutine lnsrlb
        //
        //   This subroutine calls subroutine dcsrch from the Minpack2 library
        //     to perform the line search.  Subroutine dscrch is safeguarded so
        //     that all trial points lie within the feasible region.
        //
        //   Subprograms called:
        //
        //     Minpack2 Library ... dcsrch.
        //
        //     Linpack ... dtrsl, ddot.
        //
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        //   **********
        //function
        void lnsrlb(const int& n, const double* const & l, const double* const & u, 
                const int* const& nbd, double* const & x, double& f, 
                double& fold, double& gd, double& gdold, double* const& g, 
                double* const& d, double* const& r, double* const& t,
                double* const& z, double& stp, double& dnorm, double& dtd, 
                double& xstep, double& stpmx, const int& iter, int& ifun,
                int& iback, int& nfgv, int& info, char* const & task, 
                const bool& boxed, const bool& cnstnd, char* const & csave,
                int* const & isave, double* const & dsave)
        {
                //int i;
                double a1,a2;
                double one=1.0,zero=0.0,big=1e10;
                double ftol=1.0e-3,gtol=0.9,xtol=0.1;


                if (strncmp(task,"FG_LN",5)==0) goto goto556;

                dtd = ddot(n,d,1,d,1);
                dnorm = sqrt(dtd);

                //Determine the maximum step length.

                stpmx = big;
                if (cnstnd)
                {
                        if (iter == 0)
                                stpmx = one;
                        else
                        {
                                for (int i = 1; i <= n; i++) 
                                {
                                        a1 = d[i];
                                        if (nbd[i] != 0)
                                        {
                                                if (a1 < zero && nbd[i] <= 2)
                                                {
                                                        a2 = l[i] - x[i];
                                                        if (a2 >= zero)
                                                                stpmx = zero;
                                                        else 
                                                                if (a1*stpmx < a2)
                                                                        stpmx = a2/a1;
                                                }
                                                else 
                                                        if (a1 > zero && nbd[i] >= 2)
                                                        {
                                                                a2 = u[i] - x[i];
                                                                if (a2 <= zero)
                                                                        stpmx = zero;
                                                                else 
                                                                        if (a1*stpmx > a2)
                                                                                stpmx = a2/a1;
                                                        }
                                        }
                                } // for (int i = 1; i <= n; i++) 
                        }
                }

                if (iter == 0 && !boxed)
                        stp = min(one/dnorm, stpmx);
                else
                        stp = one;

                dcopy(n,x,1,t,1);
                dcopy(n,g,1,r,1);
                fold = f;
                ifun = 0;
                iback = 0;
                strcpy(csave,"START");
goto556:
                gd = ddot(n,g,1,d,1);
                if (ifun == 0)
                {
                        gdold=gd;
                        if (gd >= zero)
                        {
                                //the directional derivative >=0.
                                //Line search is impossible.
                                info = -4;
                                return;
                        }
                }

                dcsrch(f,gd,stp,ftol,gtol,xtol,zero,stpmx,csave,isave,dsave);

                xstep = stp*dnorm;
                if (strncmp(csave,"CONV",4) != 0 && strncmp(csave,"WARN",4) != 0) 
                {
                        strcpy(task,"FG_LNSRCH");
                        ifun += 1;
                        nfgv += 1;
                        iback = ifun - 1;
                        if (stp == one)
                                dcopy(n,z,1,x,1);
                        else
                        {
                                for (int i = 1; i <= n; i++) 
                                        x[i] = stp*d[i] + t[i];
                        }
                }
                else
                        strcpy(task,"NEW_X");    
        } //lnsrlb()


        //   ************
        //
        //   Subroutine matupd
        //
        //     This subroutine updates matrices WS and WY, and forms the
        //       middle matrix in B.
        //
        //   Subprograms called:
        //
        //     Linpack ... dcopy, ddot.
        //
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        //   ************
        //function
        void matupd(const int& n, const int& m, double* const & ws, 
                double* const & wy, double* const & sy, double* const & ss, 
                double* const & d, double* const & r, int& itail, int& iupdat, 
                int& col, int& head, double& theta, double& rr, double& dr, 
                double& stp, double& dtd)
        {
                //int j;
                int pointr;
                double one=1.0;
                int idx;

                //Set pointers for matrices WS and WY.
                if (iupdat <= m)
                {
                        col = iupdat;
                        itail = (head+iupdat-2)%m + 1;
                }
                else
                {
                        itail = itail%m + 1;
                        head = head%m + 1;
                }

                //Update matrices WS and WY.
                idx = getIdx(1,itail,n)-1;
                dcopy(n,d,1,&(ws[idx]),1);
                dcopy(n,r,1,&(wy[idx]),1);

                //Set theta=yy/ys.
                theta = rr/dr;

                //Form the middle matrix in B.
                //update the upper triangle of SS,
                //and the lower triangle of SY:
                if (iupdat > m)
                {
                        //move old information
                        for (int j = 1; j <= col-1; j++)
                        {
                                dcopy(j,&(ss[getIdx(2,j+1,m)-1]),1,&(ss[getIdx(1,j,m)-1]),1);
                                dcopy(col-j,&(sy[getIdx(j+1,j+1,m)-1]),1,&(sy[getIdx(j,j,m)-1]),1);
                        }
                }

                //add new information: the last row of SY
                //and the last column of SS:
                pointr = head;
                for (int j = 1; j <= col-1; j++)
                {
                        idx = getIdx(1,pointr,n)-1;
                        sy[getIdx(col,j,m)] = ddot(n,d,1,&(wy[idx]),1);
                        ss[getIdx(j,col,m)] = ddot(n,&(ws[idx]),1,d,1);   
                        pointr = pointr%m + 1;
                }

                idx = getIdx(col,col,m);
                if (stp == one)
                        ss[idx] = dtd;
                else
                        ss[idx] = stp*stp*dtd;

                sy[idx] = dr;
        }//matupd() 


        //   ************
        //
        //   Subroutine prn1lb
        //
        //   This subroutine prints the input data, initial point, upper and
        //     lower bounds of each variable, machine precision, as well as 
        //     the headings of the output.
        //
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        //   ************
        //function
        void prn1lb(const int& n, const int& m, const double* const & l, 
                const double* const & u, const double* const & x, 
                const int& iprint, ofstream* itfile, const double& epsmch)
        {
                //int i;

                if (iprint >= 0)
                {
                        cout << "RUNNING THE L-BFGS-B CODE, " << endl
                                << "epsmch = machine precision" << endl
                                << "it    = iteration number" << endl
                                << "nf    = number of function evaluations" << endl
                                << "nint  = number of segments explored during the Cauchy search" 
                                << endl
                                << "nact  = number of active bounds at the generalized Cauchy point"
                                << endl
                                << "sub   = manner in which the subspace minimization terminated:"
                                << endl
                                <<"con = converged, bnd = a bound was reached" << endl
                                << "itls  = number of iterations performed in the line search" 
                                << endl
                                << "stepl = step length used" << endl
                                << "tstep = norm of the displacement (total step)" << endl
                                << "projg = norm of the projected gradient" << endl
                                << "f     = function value" << endl
                                << "Machine precision = " << epsmch << endl;
                        cout << "N = " << n << ",    M = " << m << endl;
                        if (iprint >= 1)
                        {
                                *itfile << "RUNNING THE L-BFGS-B CODE, epsmch = " 
                                        << "it    = iteration number" << endl
                                        << "nf    = number of function evaluations" << endl
                                        <<"nint  = number of segments explored during the Cauchy search"
                                        << endl
                                        <<"nact  = number of active bounds at generalized Cauchy point"
                                        << endl
                                        <<"sub   =manner in which the subspace minimization terminated:"
                                        << endl
                                        <<"con = converged, bnd = a bound was reached" << endl
                                        << "itls  = number of iterations performed in the line search" 
                                        << endl
                                        << "stepl = step length used" << endl
                                        << "tstep = norm of the displacement (total step)" << endl
                                        << "projg = norm of the projected gradient" << endl
                                        << "f     = function value" << endl
                                        << "Machine precision = " << epsmch << endl;
                                *itfile << "N = " << n << ",    M = " << m << endl;
                                *itfile << "   it   nf  nint  nact  sub  itls  stepl    tstep     projg"
                                        << "        f" << endl;
                                if (iprint > 100)
                                {
                                        cout << "L = ";
                                        for (int i = 1; i <= n; i++) cout << l[i] << " "; 
                                        cout << endl;
                                        cout << "X0 = ";
                                        for (int i = 1; i <= n; i++) cout << x[i] << " "; 
                                        cout << endl;
                                        cout << "U = ";
                                        for (int i = 1; i <= n; i++) cout << u[i] << " ";
                                        cout << endl;
                                }
                        }
                }
        } //prn2lb()


        //   ************
        //
        //   Subroutine prn2lb
        //
        //   This subroutine prints out new information after a successful
        //     line search. 
        //
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        //   ************
        //function
        void prn2lb(const int& n, const double* const & x, const double& f, 
                const double* const & g, const int& iprint, 
                ofstream* itfile, const int& iter, const int& nfgv, 
                const int& nact, const double& sbgnrm, const int& nint, 
                char* const & word, const int& iword, const int& iback, 
                const double& stp, const double& xstep)
        {
                //int i;
                int imod;

                //'word' records the status of subspace solutions.
                if (iword == 0)
                {
                        //the subspace minimization converged.
                        strcpy(word, "con");
                }
                else 
                        if (iword == 1)
                        {
                                //the subspace minimization stopped at a bound.
                                strcpy(word, "bnd");
                        }
                        else 
                                if (iword == 5)
                                {
                                        //the truncated Newton step has been used.
                                        strcpy(word, "TNT");

                                }      
                                else
                                        strcpy(word, "---");

                if (iprint >= 99)
                {
                        cout << "LINE SEARCH " << iback <<" times; norm of step = " <<xstep<<endl;
                        cout << "At iterate " << iter << " f=" << f <<"  |proj g|="<<sbgnrm<<endl;
                        if (iprint > 100)
                        {       
                                cout << "X = " << endl;
                                for (int i = 1; i <= n; i++) cout << x[i] << " ";
                                cout << endl;
                                cout << "G = " << endl;
                                for (int i = 1; i <= n; i++) cout << g[i] << " ";
                                cout << endl;
                        }
                }
                else 
                        if (iprint > 0)
                        {
                                imod = iter%iprint;
                                if (imod == 0) 
                                        cout << "At iterate " << iter << " f="<<f<< "  |proj g|="<<sbgnrm<<endl;
                        }
                        if (iprint >= 1)
                                *itfile << " " << iter << " " << nfgv << " " << nint << " " << nact << " "
                                << word << " " << iback << " " << stp << " " << xstep << " " 
                                << sbgnrm << " " << f << endl;
        }//prn2lb()


        //   ************
        //
        //   Subroutine prn3lb
        //
        //   This subroutine prints out information when either a built-in
        //     convergence test is satisfied or when an error message is
        //     generated.
        //     
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        //   ************
        //function
        void prn3lb(const int& n, const double* const & x, double& f, 
                const char* const & task, const int& iprint, 
                const int& info, ofstream* itfile, const int& iter, 
                const int& nfgv, const int& nintol, const int& nskip, 
                const int& nact, const double& sbgnrm, const double& time, 
                const int& nint, const char* const & word, const int& iback, 
                const double& stp,const double& xstep, const int& k, 
                const double& cachyt, const double& sbtime, const double& lnscht)
        {
                //int i;
                if (strncmp(task,"ERROR",5)==0) goto goto999;

                if (iprint >= 0)
                {
                        cout << "           * * *" 
                                << "Tit   = total number of iterations" << endl
                                << "Tnf   = total number of function evaluations" << endl
                                << "Tnint = total number of segments explored during"
                                << " Cauchy searches" << endl
                                << "Skip  = number of BFGS updates skipped" << endl
                                << "Nact  = number of active bounds at final generalized"
                                << " Cauchy point" << endl
                                << "Projg = norm of the final projected gradient" << endl
                                << "F     = final function value" << endl
                                << "           * * *" << endl;

                        cout << "   N   Tit  Tnf  Tnint  Skip  Nact     Projg        F" << endl;
                        cout << " " << n << " " << iter << " " << nfgv << " " << nintol << " " 
                                << nskip << " " << nact << " " << sbgnrm << " " << f << endl;
                        if (iprint >= 100)
                        {
                                cout << "X = ";
                                for (int i = 1; i <= n; i++)
                                        cout << x[i] << " ";
                                cout << endl;
                        }
                        if (iprint >= 1) cout << " F = " << f << endl;
                }

goto999:

                if (iprint >= 0)
                {
                        cout << task << endl;
                        if (info != 0)
                        {
                                if (info == -1)
                                        cout << "Matrix in 1st Cholesky factorization in formk is not " 
                                        << "Pos. Def." << endl;
                                if (info == -2)
                                        cout << "Matrix in 2st Cholesky factorization in formk is not "
                                        << "Pos. Def." << endl;
                                if (info == -3) 
                                        cout << "Matrix in the Cholesky factorization in formt is not "
                                        << "Pos. Def." << endl;
                                if (info == -4)
                                        cout << "Derivative >= 0, backtracking line search impossible."
                                        << endl << "Previous x, f and g restored." << endl
                                        << "Possible causes: 1 error in function or gradient "
                                        <<"evaluation;" << endl
                                        << "2 rounding errors dominate computation." << endl;
                                if (info == -5) 
                                        cout << "Warning:  more than 10 function and gradient" << endl
                                        << "evaluations in the last line search.  Termination"
                                        << "may possibly be caused by a bad search direction."<<endl;
                                if (info == -6) 
                                        cout << " Input nbd(" << k << ") is invalid."<<endl;
                                if (info == -7)
                                        cout << " l(" << k << ") > u(" << k << ").  No feasible solution."
                                        << endl;
                                if (info == -8) 
                                        cout << "The triangular system is singular" << endl;
                                if (info == -9)
                                        cout << " Line search cannot locate an adequate point after"<<endl
                                        << " 20 function and gradient evaluations. Previous x, f and"
                                        << " g restored." << endl
                                        << "Possible causes: 1 error in function or gradient "
                                        << "evaluation; 2 rounding error dominate computation" 
                                        << endl;
                        }

                        if (iprint >= 1) 
                                cout << "Cauchy                time " << cachyt << " seconds."<<endl
                                << "Subspace minimization time " << sbtime << " seconds."<<endl
                                << "Line search           time " << lnscht << " seconds."<<endl;

                        cout << "Total User time " << time << " seconds." << endl;

                        if (iprint >= 1)
                        {
                                if (info == -4 || info == -9)
                                {              
                                        *itfile << " " << n << " " << iter << " " << nfgv << " " << nintol 
                                                << " " << nskip << " " << nact << " " << sbgnrm << " " << f
                                                << endl;
                                }
                                *itfile << task << endl;
                                if (info != 0)
                                {
                                        if (info == -1)
                                                *itfile << "Matrix in 1st Cholesky factorization in formk is not"
                                                << " Pos. Def." << endl;
                                        if (info == -2)
                                                *itfile << "Matrix in 2st Cholesky factorization in formk is not"
                                                << " Pos. Def." << endl;
                                        if (info == -3) 
                                                *itfile << "Matrix in the Cholesky factorization in formt is not"
                                                << " Pos. Def." << endl;
                                        if (info == -4)
                                                *itfile << "Derivative>=0, backtracking line search impossible."
                                                << endl << "Previous x, f and g restored." << endl
                                                << "Possible causes: 1 error in function or gradient "
                                                <<"evaluation;" << endl
                                                << "2 rounding errors dominate computation." << endl;
                                        if (info == -5) 
                                                *itfile << "Warning:  more than 10 function and gradient" << endl
                                                << "evaluations in the last line search.  Termination"
                                                << "may possibly be caused by a bad search direction."
                                                << endl;
                                        if (info == -8) 
                                                cout << "The triangular system is singular" << endl;
                                        if (info == -9)
                                                cout << " Line search cannot locate an adequate point after"
                                                << endl
                                                << " 20 function and gradient evaluations. Previous x, f "
                                                << "and g restored." << endl
                                                << "Possible causes: 1 error in function or gradient "
                                                << "evaluation; 2 rounding error dominate computation" 
                                                << endl;
                                }
                                *itfile << "Total User time " << time << " seconds." << endl;
                        }
                }    
        }//prn3lb()


        //   ************
        //
        //   Subroutine projgr
        //
        //   This subroutine computes the infinity norm of the projected
        //     gradient.
        //
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        //   ************
        //functions
        void projgr(const int& n, const double* const & l, const double* const & u, 
                const int* const & nbd, const double* const & x, 
                const double* const & g, double& sbgnrm)
        {
                //int i;
                double gi;
                //double one=1.0;
                double zero=0.0;

                sbgnrm = zero;
                for (int i = 1; i <= n; i++)
                {
                        gi = g[i];
                        if (nbd[i] != 0)
                        {
                                if (gi < zero)
                                {
                                        if (nbd[i] >= 2) gi = max((x[i]-u[i]),gi);
                                }
                                else
                                {
                                        if (nbd[i] <= 2) gi = min((x[i]-l[i]),gi);
                                }
                        }
                        sbgnrm = max(sbgnrm,fabs(gi));
                }
        }//projgr()


        //   ************
        //
        //   Subroutine subsm
        //
        //   Given xcp, l, u, r, an index set that specifies
        //       the active set at xcp, and an l-BFGS matrix B 
        //       (in terms of WY, WS, SY, WT, head, col, and theta), 
        //       this subroutine computes an approximate solution
        //       of the subspace problem
        //
        //      (P)   min Q(x) = r'(x-xcp) + 1/2 (x-xcp)' B (x-xcp)
        //
        //           subject to l<=x<=u
        //                      x_i=xcp_i for all i in A(xcp)
        //                    
        //      along the subspace unconstrained Newton direction 
        //      
        //         d = -(Z'BZ)^(-1) r.
        //
        //     The formula for the Newton direction, given the L-BFGS matrix
        //     and the Sherman-Morrison formula, is
        //
        //         d = (1/theta)r + (1/theta*2) Z'WK^(-1)W'Z r.
        //  
        //     where
        //               K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
        //                   [L_a -R_z           theta*S'AA'S ]
        //
        //   Note that this procedure for computing d differs 
        //   from that described in [1]. One can show that the matrix K is
        //   equal to the matrix M^[-1]N in that paper.
        //
        //   n is an integer variable.
        //     On entry n is the dimension of the problem.
        //     On exit n is unchanged.
        //
        //   m is an integer variable.
        //     On entry m is the maximum number of variable metric corrections
        //       used to define the limited memory matrix.
        //     On exit m is unchanged.
        //
        //   nsub is an integer variable.
        //     On entry nsub is the number of free variables.
        //     On exit nsub is unchanged.
        //
        //   ind is an integer array of dimension nsub.
        //     On entry ind specifies the coordinate indices of free variables.
        //     On exit ind is unchanged.
        //
        //   l is a double precision array of dimension n.
        //     On entry l is the lower bound of x.
        //     On exit l is unchanged.
        //
        //   u is a double precision array of dimension n.
        //     On entry u is the upper bound of x.
        //     On exit u is unchanged.
        //
        //   nbd is a integer array of dimension n.
        //     On entry nbd represents the type of bounds imposed on the
        //       variables, and must be specified as follows:
        //       nbd(i)=0 if x(i) is unbounded,
        //              1 if x(i) has only a lower bound,
        //              2 if x(i) has both lower and upper bounds, and
        //              3 if x(i) has only an upper bound.
        //     On exit nbd is unchanged.
        //
        //   x is a double precision array of dimension n.
        //     On entry x specifies the Cauchy point xcp. 
        //     On exit x(i) is the minimizer of Q over the subspace of
        //                                                      free variables. 
        //
        //   d is a double precision array of dimension n.
        //     On entry d is the reduced gradient of Q at xcp.
        //     On exit d is the Newton direction of Q. 
        //
        //   ws and wy are double precision arrays;
        //   theta is a double precision variable;
        //   col is an integer variable;
        //   head is an integer variable.
        //     On entry they store the information defining the
        //                                        limited memory BFGS matrix:
        //       ws(n,m) stores S, a set of s-vectors;
        //       wy(n,m) stores Y, a set of y-vectors;
        //       theta is the scaling factor specifying B_0 = theta I;
        //       col is the number of variable metric corrections stored;
        //       head is the location of the 1st s- (or y-) vector in S (or Y).
        //     On exit they are unchanged.
        //
        //   iword is an integer variable.
        //     On entry iword is unspecified.
        //     On exit iword specifies the status of the subspace solution.
        //       iword = 0 if the solution is in the box,
        //               1 if some bound is encountered.
        //
        //   wv is a double precision working array of dimension 2m.
        //
        //   wn is a double precision array of dimension 2m x 2m.
        //     On entry the upper triangle of wn stores the LEL^T factorization
        //       of the indefinite matrix
        //
        //            k = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
        //                [L_a -R_z           theta*S'AA'S ]
        //                                                  where E = [-I  0]
        //                                                            [ 0  I]
        //     On exit wn is unchanged.
        //
        //   iprint is an INTEGER variable that must be set by the user.
        //     It controls the frequency and type of output generated:
        //      iprint<0    no output is generated;
        //      iprint=0    print only one line at the last iteration;
        //      0<iprint<99 print also f and |proj g| every iprint iterations;
        //      iprint=99   print details of every iteration except n-vectors;
        //      iprint=100  print also the changes of active set and final x;
        //      iprint>100  print details of every iteration including x and g;
        //     When iprint > 0, the file iterate.dat will be created to
        //                      summarize the iteration.
        //
        //   info is an integer variable.
        //     On entry info is unspecified.
        //     On exit info = 0       for normal return,
        //                  = nonzero for abnormal return 
        //                                when the matrix K is ill-conditioned.
        //
        //   Subprograms called:
        //
        //     Linpack dtrsl.
        //
        //
        //   References:
        //
        //     [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
        //     memory algorithm for bound constrained optimization'',
        //     SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
        //
        //
        //
        //                         *  *  *
        //
        //   NEOS, November 1994. (Latest revision June 1996.)
        //   Optimization Technology Center.
        //   Argonne National Laboratory and Northwestern University.
        //   Written by
        //                      Ciyou Zhu
        //   in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
        //
        //
        //   ************
        //function
        void subsm(const int& n, const int& m, const int& nsub, const int* const& ind,
                const double* const & l, const double* const & u, 
                const int* const & nbd, double* const & x, double* const & d, 
                double* const & ws, double* const & wy, const double& theta, 
                const int& col, const int& head, int& iword, double* const & wv, 
                double* const & wn, const int& iprint, int& info) 
        {
                //int jy,i,j;
                int pointr,m2,col2,ibd=0,js,k;
                double alpha,dk,temp1,temp2;
                double one=1.0,zero=0.0;
                int idx;

                if (nsub <= 0) return;
                if (iprint >= 99) 
                        cout << "----------------SUBSM entered-----------------" << endl;

                //Compute wv = W'Zd.

                pointr = head;
                for (int i = 1; i <= col; i++)
                {
                        temp1 = zero;
                        temp2 = zero;
                        for (int j = 1; j <= nsub; j++)
                        {
                                k = ind[j];
                                idx = getIdx(k,pointr,n);
                                temp1 += wy[idx]*d[j];
                                temp2 += ws[idx]*d[j];
                        }
                        wv[i] = temp1;
                        wv[col + i] = theta*temp2;
                        pointr = pointr%m + 1;
                }


                //Compute wv:=K^(-1)wv.

                m2 = 2*m;
                col2 = 2*col;
                dtrsl(wn,m2,col2,wv,11,info);

                if (info != 0) return;
                for (int i = 1; i <= col; i++)
                        wv[i] = -wv[i];
                dtrsl(wn,m2,col2,wv,01,info);

                if (info != 0) return;

                //Compute d = (1/theta)d + (1/theta**2)Z'W wv.

                pointr = head;
                for (int jy = 1; jy <= col; jy++)
                {
                        js = col + jy;
                        for (int i = 1; i <= nsub; i++)        
                        {
                                k = ind[i];
                                idx = getIdx(k,pointr,n);
                                d[i] += wy[idx]*wv[jy]/theta + ws[idx]*wv[js];
                        }
                        pointr = pointr%m + 1;
                }

                for (int i = 1; i <= nsub; i++)
                        d[i] /= theta;

                //Backtrack to the feasible region.

                alpha = one;
                temp1 = alpha;

                for (int i = 1; i <= nsub; i++)
                {
                        k = ind[i];
                        dk = d[i];
                        if (nbd[k] != 0)
                        {
                                if (dk < zero && nbd[k] <= 2)
                                {
                                        temp2 = l[k] - x[k];
                                        if (temp2 >= zero)
                                                temp1 = zero;
                                        else 
                                                if (dk*alpha < temp2)
                                                        temp1 = temp2/dk;               
                                }
                                else 
                                        if (dk > zero && nbd[k] >= 2)
                                        {
                                                temp2 = u[k] - x[k];
                                                if (temp2 <= zero)
                                                        temp1 = zero;
                                                else
                                                        if (dk*alpha > temp2)
                                                                temp1 = temp2/dk;
                                        }
                                        if (temp1 < alpha)
                                        {
                                                alpha = temp1;
                                                ibd = i;
                                        }
                        }
                }


                if (alpha < one)
                {
                        dk = d[ibd];
                        k = ind[ibd];
                        if (dk > zero)
                        {
                                x[k] = u[k];
                                d[ibd] = zero;
                        }
                        else 
                                if (dk < zero)
                                {
                                        x[k] = l[k];
                                        d[ibd] = zero;
                                }
                }

                for (int i = 1; i <= nsub; i++)
                {
                        k = ind[i];
                        x[k] += alpha*d[i];
                } 

                if (iprint >= 99)
                {
                        if (alpha < one)
                                cout << "ALPHA = " << alpha << " backtrack to the BOX" << endl;
                        else
                                cout  << "SM solution inside the box" << endl;

                        if (iprint >100)
                        {
                                cout << "Subspace solution X = ";
                                for (int i = 1; i <= n; i++) cout << x[i] << " ";
                                cout << endl;
                        }
                }

                if (alpha < one)
                        iword = 1;
                else
                        iword = 0;

                if (iprint >= 99) 
                        cout << "----------------exit SUBSM --------------------" << endl;    
        }//subsm()


        //   **********
        //
        //   Subroutine dcsrch
        //
        //   This subroutine finds a step that satisfies a sufficient
        //   decrease condition and a curvature condition.
        //
        //   Each call of the subroutine updates an interval with 
        //   endpoints stx and sty. The interval is initially chosen 
        //   so that it contains a minimizer of the modified function
        //
        //         psi(stp) = f(stp) - f(0) - ftol*stp*f'(0).
        //
        //   If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the
        //   interval is chosen so that it contains a minimizer of f. 
        //
        //   The algorithm is designed to find a step that satisfies 
        //   the sufficient decrease condition 
        //
        //         f(stp) <= f(0) + ftol*stp*f'(0),
        //
        //   and the curvature condition
        //
        //         abs(f'(stp)) <= gtol*abs(f'(0)).
        //
        //   If ftol is less than gtol and if, for example, the function
        //   is bounded below, then there is always a step which satisfies
        //   both conditions. 
        //
        //   If no step can be found that satisfies both conditions, then 
        //   the algorithm stops with a warning. In this case stp only 
        //   satisfies the sufficient decrease condition.
        //
        //   A typical invocation of dcsrch has the following outline:
        //
        //   task = 'START'
        //  10 continue
        //      call dcsrch( ... )
        //      if (task == 'FG') then
        //         Evaluate the function and the gradient at stp 
        //         goto 10
        //         end if
        //
        //   NOTE: The user must no alter work arrays between calls.
        //
        //   The subroutine statement is
        //
        //      subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax,
        //                        task,isave,dsave)
        //   where
        //
        //     f is a double precision variable.
        //       On initial entry f is the value of the function at 0.
        //          On subsequent entries f is the value of the 
        //          function at stp.
        //       On exit f is the value of the function at stp.
        //
        //    g is a double precision variable.
        //       On initial entry g is the derivative of the function at 0.
        //          On subsequent entries g is the derivative of the 
        //          function at stp.
        //       On exit g is the derivative of the function at stp.
        //
        //      stp is a double precision variable. 
        //       On entry stp is the current estimate of a satisfactory 
        //          step. On initial entry, a positive initial estimate 
        //          must be provided. 
        //       On exit stp is the current estimate of a satisfactory step
        //          if task = 'FG'. If task = 'CONV' then stp satisfies
        //          the sufficient decrease and curvature condition.
        //
        //     ftol is a double precision variable.
        //       On entry ftol specifies a nonnegative tolerance for the 
        //          sufficient decrease condition.
        //       On exit ftol is unchanged.
        //
        //     gtol is a double precision variable.
        //       On entry gtol specifies a nonnegative tolerance for the 
        //          curvature condition. 
        //       On exit gtol is unchanged.
        //
        //      xtol is a double precision variable.
        //       On entry xtol specifies a nonnegative relative tolerance
        //          for an acceptable step. The subroutine exits with a
        //          warning if the relative difference between sty and stx
        //          is less than xtol.
        //       On exit xtol is unchanged.
        //
        //      stpmin is a double precision variable.
        //       On entry stpmin is a nonnegative lower bound for the step.
        //       On exit stpmin is unchanged.
        //
        //    stpmax is a double precision variable.
        //       On entry stpmax is a nonnegative upper bound for the step.
        //       On exit stpmax is unchanged.
        //
        //     task is a character variable of length at least 60.
        //       On initial entry task must be set to 'START'.
        //       On exit task indicates the required action:
        //
        //          If task(1:2) = 'FG' then evaluate the function and 
        //          derivative at stp and call dcsrch again.
        //
        //          If task(1:4) = 'CONV' then the search is successful.
        //
        //          If task(1:4) = 'WARN' then the subroutine is not able
        //          to satisfy the convergence conditions. The exit value of
        //          stp contains the best point found during the search.
        //
        //          If task(1:5) = 'ERROR' then there is an error in the
        //          input arguments.
        //
        //       On exit with convergence, a warning or an error, the
        //          variable task contains additional information.
        //
        //     isave is an integer work array of dimension 2.
        //       
        //     dsave is a double precision work array of dimension 13.
        //
        //   Subprograms called
        //
        //       MINPACK-2 ... dcstep
        //
        //   MINPACK-1 Project. June 1983.
        //   Argonne National Laboratory. 
        //   Jorge J. More' and David J. Thuente.
        //
        //   MINPACK-2 Project. October 1993.
        //   Argonne National Laboratory and University of Minnesota. 
        //   Brett M. Averick, Richard G. Carter, and Jorge J. More'. 
        //
        //   **********
        //function
        void dcsrch(double& f, double& g, double& stp, const double& ftol, 
                const double& gtol, const double& xtol, const double& stpmin, 
                const double& stpmax, char* const & task, int* const & isave, 
                double* const & dsave)
        {
                double zero=0.0,p5=0.5,p66=0.66;
                double xtrapl=1.1,xtrapu=4.0;

                bool brackt;
                int stage;
                double finit,ftest,fm,fx,fxm,fy,fym,ginit,gtest,
                        gm,gx,gxm,gy,gym,stx,sty,stmin,stmax,width,width1;

                //Initialization block.

                if (strncmp(task,"START",5)==0)
                {
                        //Check the input arguments for errors.

                        if (stp < stpmin)    strcpy(task, "ERROR: STP < STPMIN");
                        if (stp > stpmax)    strcpy(task, "ERROR: STP > STPMAX");
                        if (g >zero)         strcpy(task, "ERROR: INITIAL G >ZERO");
                        if (ftol < zero)     strcpy(task, "ERROR: FTOL < ZERO");
                        if (gtol < zero)     strcpy(task, "ERROR: GTOL < ZERO");
                        if (xtol < zero)     strcpy(task, "ERROR: XTOL < ZERO");
                        if (stpmin < zero)   strcpy(task, "ERROR: STPMIN < ZERO");
                        if (stpmax < stpmin) strcpy(task, "ERROR: STPMAX < STPMIN");

                        //Exit if there are errors on input.
                        if (strncmp(task,"ERROR",5)==0) return;

                        //Initialize local variables.

                        brackt = false;
                        stage = 1;
                        finit = f;
                        ginit = g;
                        gtest = ftol*ginit;
                        width = stpmax - stpmin;
                        width1 = width/p5;

                        //The variables stx, fx, gx contain the values of the step, 
                        //function, and derivative at the best step. 
                        //The variables sty, fy, gy contain the value of the step, 
                        //function, and derivative at sty.
                        //The variables stp, f, g contain the values of the step, 
                        //function, and derivative at stp.

                        stx = zero;
                        fx = finit;
                        gx = ginit;
                        sty = zero;
                        fy = finit;
                        gy = ginit;
                        stmin = zero;
                        stmax = stp + xtrapu*stp;
                        strcpy(task, "FG");

                        goto goto1000;
                }
                else
                {
                        //Restore local variables.

                        if (isave[1] == 1)
                                brackt = true;
                        else
                                brackt = false;

                        stage  = isave[2]; 
                        ginit  = dsave[1]; 
                        gtest  = dsave[2]; 
                        gx     = dsave[3]; 
                        gy     = dsave[4]; 
                        finit  = dsave[5]; 
                        fx     = dsave[6]; 
                        fy     = dsave[7]; 
                        stx    = dsave[8]; 
                        sty    = dsave[9]; 
                        stmin  = dsave[10]; 
                        stmax  = dsave[11]; 
                        width  = dsave[12]; 
                        width1 = dsave[13]; 
                }

                //If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the
                //algorithm enters the second stage.

                ftest = finit + stp*gtest;
                if (stage == 1 && f <= ftest && g >= zero) stage = 2;

                //Test for warnings.

                if (brackt && (stp <= stmin || stp >= stmax))
                        strcpy(task, "WARNING: ROUNDING ERRORS PREVENT PROGRESS");
                if (brackt && stmax - stmin <= xtol*stmax) 
                        strcpy(task, "WARNING: XTOL TEST SATISFIED");
                if (stp == stpmax && f <= ftest && g <= gtest) 
                        strcpy(task, "WARNING: STP = STPMAX");
                if (stp == stpmin && (f > ftest || g >= gtest)) 
                        strcpy(task, "WARNING: STP = STPMIN");

                //Test for convergence.

                if (f <= ftest && fabs(g) <= gtol*(-ginit)) 
                        strcpy(task, "CONVERGENCE");

                //Test for termination.

                if (strncmp(task,"WARN",4)==0 || strncmp(task,"CONV",4)==0) goto goto1000;

                //A modified function is used to predict the step during the
                //first stage if a lower function value has been obtained but 
                //the decrease is not sufficient.

                if (stage == 1 && f <= fx && f > ftest)
                {
                        //Define the modified function and derivative values.

                        fm = f - stp*gtest;
                        fxm = fx - stx*gtest;
                        fym = fy - sty*gtest;
                        gm = g - gtest;
                        gxm = gx - gtest;
                        gym = gy - gtest;

                        //Call dcstep to update stx, sty, and to compute the new step.

                        dcstep(stx,fxm,gxm,sty,fym,gym,stp,fm,gm,
                                brackt,stmin,stmax);

                        //Reset the function and derivative values for f.

                        fx = fxm + stx*gtest;
                        fy = fym + sty*gtest;
                        gx = gxm + gtest;
                        gy = gym + gtest;
                }
                else
                {
                        //Call dcstep to update stx, sty, and to compute the new step.

                        dcstep(stx,fx,gx,sty,fy,gy,stp,f,g,
                                brackt,stmin,stmax);
                }

                //Decide if a bisection step is needed.

                if (brackt)
                {
                        if (fabs(sty-stx) >= p66*width1) stp = stx + p5*(sty - stx);
                        width1 = width;
                        width = fabs(sty-stx);
                }

                //Set the minimum and maximum steps allowed for stp.

                if (brackt)
                {
                        stmin = min(stx,sty);
                        stmax = max(stx,sty);
                }
                else
                {
                        stmin = stp + xtrapl*(stp - stx);
                        stmax = stp + xtrapu*(stp - stx);
                }

                //Force the step to be within the bounds stpmax and stpmin.

                stp = max(stp,stpmin);
                stp = min(stp,stpmax);

                //If further progress is not possible, let stp be the best
                //point obtained during the search.

                if (brackt && (stp <= stmin || stp >= stmax)
                        || (brackt && stmax-stmin <= xtol*stmax)) stp = stx;

                //Obtain another function and derivative.
                strcpy(task,"FG");

goto1000:

                //Save local variables.

                if (brackt)
                        isave[1]  = 1;
                else    
                        isave[1]  = 0;

                isave[2]  = stage;
                dsave[1]  = ginit;
                dsave[2]  = gtest;
                dsave[3]  = gx;
                dsave[4]  = gy;
                dsave[5]  = finit;
                dsave[6]  = fx;
                dsave[7]  = fy;
                dsave[8]  = stx;
                dsave[9]  = sty;
                dsave[10] = stmin;
                dsave[11] = stmax;
                dsave[12] = width;
                dsave[13] = width1;
        }//dcsrch()


        //   **********
        //
        //   Subroutine dcstep
        //
        //   This subroutine computes a safeguarded step for a search
        //   procedure and updates an interval that contains a step that
        //   satisfies a sufficient decrease and a curvature condition.
        //
        //   The parameter stx contains the step with the least function
        //   value. If brackt is set to true then a minimizer has
        //   been bracketed in an interval with endpoints stx and sty.
        //   The parameter stp contains the current step. 
        //   The subroutine assumes that if brackt is set to true then
        //
        //         min(stx,sty) < stp < max(stx,sty),
        //
        //   and that the derivative at stx is negative in the direction 
        //   of the step.
        //
        //   The subroutine statement is
        //
        //     subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt,
        //                       stpmin,stpmax)
        //
        //   where
        //
        //     stx is a double precision variable.
        //       On entry stx is the best step obtained so far and is an
        //          endpoint of the interval that contains the minimizer. 
        //       On exit stx is the updated best step.
        //
        //     fx is a double precision variable.
        //       On entry fx is the function at stx.
        //       On exit fx is the function at stx.
        //
        //     dx is a double precision variable.
        //       On entry dx is the derivative of the function at 
        //          stx. The derivative must be negative in the direction of 
        //          the step, that is, dx and stp - stx must have opposite 
        //          signs.
        //       On exit dx is the derivative of the function at stx.
        //
        //     sty is a double precision variable.
        //       On entry sty is the second endpoint of the interval that 
        //          contains the minimizer.
        //       On exit sty is the updated endpoint of the interval that 
        //          contains the minimizer.
        //
        //     fy is a double precision variable.
        //       On entry fy is the function at sty.
        //       On exit fy is the function at sty.
        //
        //     dy is a double precision variable.
        //       On entry dy is the derivative of the function at sty.
        //       On exit dy is the derivative of the function at the exit sty.
        //
        //     stp is a double precision variable.
        //       On entry stp is the current step. If brackt is set to true
        //          then on input stp must be between stx and sty. 
        //       On exit stp is a new trial step.
        //
        //     fp is a double precision variable.
        //       On entry fp is the function at stp
        //       On exit fp is unchanged.
        //
        //     dp is a double precision variable.
        //       On entry dp is the the derivative of the function at stp.
        //       On exit dp is unchanged.
        //
        //     brackt is an logical variable.
        //       On entry brackt specifies if a minimizer has been bracketed.
        //          Initially brackt must be set to false
        //       On exit brackt specifies if a minimizer has been bracketed.
        //          When a minimizer is bracketed brackt is set to true
        //
        //     stpmin is a double precision variable.
        //       On entry stpmin is a lower bound for the step.
        //       On exit stpmin is unchanged.
        //
        //     stpmax is a double precision variable.
        //       On entry stpmax is an upper bound for the step.
        //       On exit stpmax is unchanged.
        //
        //   MINPACK-1 Project. June 1983
        //   Argonne National Laboratory. 
        //   Jorge J. More' and David J. Thuente.
        //
        //   MINPACK-2 Project. October 1993.
        //   Argonne National Laboratory and University of Minnesota. 
        //   Brett M. Averick and Jorge J. More'.
        //
        //   **********
        //function
        void dcstep(double& stx, double& fx, double& dx, double& sty, double& fy, 
                double& dy, double& stp, const double& fp, const double& dp, 
                bool& brackt, const double& stpmin, const double& stpmax)
        {
                double zero=0.0,p66=0.66,two=2.0,three=3.0;
                double gamma,p,q,r,s,sgnd,stpc,stpf,stpq,theta;

                sgnd = dp*(dx/fabs(dx));

                //First case: A higher function value. The minimum is bracketed. 
                //If the cubic step is closer to stx than the quadratic step, the 
                //cubic step is taken, otherwise the average of the cubic and 
                //quadratic steps is taken.

                if (fp > fx)
                {
                        theta = three*(fx - fp)/(stp - stx) + dx + dp;
                        s = max(fabs(theta),fabs(dx),fabs(dp));
                        gamma = s*sqrt( (theta/s)*(theta/s) - (dx/s)*(dp/s) );
                        if (stp < stx) gamma = -gamma;
                        p = (gamma - dx) + theta;
                        q = ((gamma - dx) + gamma) + dp;
                        r = p/q;
                        stpc = stx + r*(stp - stx);
                        stpq = stx + ((dx/((fx - fp)/(stp - stx) + dx))/two)*(stp - stx);
                        if (fabs(stpc-stx) < fabs(stpq-stx))
                                stpf = stpc;
                        else
                                stpf = stpc + (stpq - stpc)/two;
                        brackt = true;
                }

                //Second case: A lower function value and derivatives of opposite 
                //sign. The minimum is bracketed. If the cubic step is farther from 
                //stp than the secant step, the cubic step is taken, otherwise the 
                //secant step is taken.

                else 
                        if (sgnd < zero)
                        {
                                theta = three*(fx - fp)/(stp - stx) + dx + dp;
                                s = max(fabs(theta),fabs(dx),fabs(dp));
                                gamma = s*sqrt( (theta/s)*(theta/s) - (dx/s)*(dp/s) );
                                if (stp > stx) gamma = -gamma;
                                p = (gamma - dp) + theta;
                                q = ((gamma - dp) + gamma) + dx;
                                r = p/q;
                                stpc = stp + r*(stx - stp);
                                stpq = stp + (dp/(dp - dx))*(stx - stp);
                                if (fabs(stpc-stp) > fabs(stpq-stp))
                                        stpf = stpc;
                                else
                                        stpf = stpq;
                                brackt = true;
                        }

                        //Third case: A lower function value, derivatives of the same sign,
                        //and the magnitude of the derivative decreases.

                        else 
                                if (fabs(dp) < fabs(dx))
                                {
                                        //The cubic step is computed only if the cubic tends to infinity 
                                        //in the direction of the step or if the minimum of the cubic
                                        //is beyond stp. Otherwise the cubic step is defined to be the 
                                        //secant step.

                                        theta = three*(fx - fp)/(stp - stx) + dx + dp;
                                        s = max(fabs(theta),fabs(dx),fabs(dp));

                                        //The case gamma = 0 only arises if the cubic does not tend
                                        //to infinity in the direction of the step.

                                        gamma = s*sqrt(max(zero, (theta/s)*(theta/s)-(dx/s)*(dp/s)));
                                        if (stp > stx) gamma = -gamma;
                                        p = (gamma - dp) + theta;
                                        q = (gamma + (dx - dp)) + gamma;
                                        r = p/q;
                                        if (r < zero && gamma != zero)
                                                stpc = stp + r*(stx - stp);
                                        else 
                                                if (stp > stx)
                                                        stpc = stpmax;
                                                else
                                                        stpc = stpmin;

                                        stpq = stp + (dp/(dp - dx))*(stx - stp);

                                        if (brackt)
                                        {

                                                //A minimizer has been bracketed. If the cubic step is 
                                                //closer to stp than the secant step, the cubic step is 
                                                //taken, otherwise the secant step is taken.

                                                if (fabs(stpc-stp) < fabs(stpq-stp))
                                                        stpf = stpc;
                                                else
                                                        stpf = stpq;

                                                if (stp > stx) 
                                                        stpf = min(stp+p66*(sty-stp),stpf);
                                                else
                                                        stpf = max(stp+p66*(sty-stp),stpf);
                                        }
                                        else
                                        {
                                                //A minimizer has not been bracketed. If the cubic step is 
                                                //farther from stp than the secant step, the cubic step is 
                                                //taken, otherwise the secant step is taken.

                                                if (fabs(stpc-stp) > fabs(stpq-stp))
                                                        stpf = stpc;
                                                else
                                                        stpf = stpq;

                                                stpf = min(stpmax,stpf);
                                                stpf = max(stpmin,stpf);
                                        }
                                }
                                //Fourth case: A lower function value, derivatives of the same sign, 
                                //and the magnitude of the derivative does not decrease. If the 
                                //minimum is not bracketed, the step is either stpmin or stpmax, 
                                //otherwise the cubic step is taken.
                                else
                                {
                                        if (brackt)
                                        {
                                                theta = three*(fp - fy)/(sty - stp) + dy + dp;
                                                s = max(fabs(theta),fabs(dy),fabs(dp));
                                                gamma = s*sqrt((theta/s)*(theta/s) - (dy/s)*(dp/s));
                                                if (stp > sty) gamma = -gamma;
                                                p = (gamma - dp) + theta;
                                                q = ((gamma - dp) + gamma) + dy;
                                                r = p/q;
                                                stpc = stp + r*(sty - stp);
                                                stpf = stpc;
                                        }
                                        else 
                                                if (stp > stx)
                                                        stpf = stpmax;
                                                else
                                                        stpf = stpmin;
                                }

                                //Update the interval which contains a minimizer.

                                if (fp > fx)
                                {
                                        sty = stp;
                                        fy = fp;
                                        dy = dp;
                                }
                                else
                                {
                                        if (sgnd < zero)
                                        {
                                                sty = stx;
                                                fy = fx;
                                                dy = dx;
                                        }
                                        stx = stp;
                                        fx = fp;
                                        dx = dp;
                                }

                                //Compute the new step.

                                stp = stpf;    
        }//dcstep()


        //   **********
        //
        //   Function dnrm2
        //
        //   Given a vector x of length n, this function calculates the
        //   Euclidean norm of x with stride incx.
        //
        //   The function statement is
        //
        //     double precision function dnrm2(n,x,incx)
        //
        //   where
        //
        //     n is a positive integer input variable.
        //
        //     x is an input array of length n.
        //
        //     incx is a positive integer variable that specifies the 
        //       stride of the vector.
        //
        //   Subprograms called
        //
        //     FORTRAN-supplied ... abs, max, sqrt
        //
        //   MINPACK-2 Project. February 1991.
        //   Argonne National Laboratory.
        //   Brett M. Averick.
        //
        //   **********
        //function
        double dnrm2(int& n,double* const & x, int& incx)
        {
                //int i;
                double scale = 0.0;
                double ddnrm2 = 0.0;

                for (int i = 1; i <= n; i+=incx)
                        scale = max(scale, fabs(x[i]));

                if (scale == 0.0) return ddnrm2;

                for (int i = 1; i <= n; i+=incx)  
                        ddnrm2 += (x[i]/scale)*(x[i]/scale);

                ddnrm2 = scale*sqrt(ddnrm2);
                return ddnrm2;
        }


        //   **********
        //
        //   Subroutine dpeps
        //
        //   This subroutine computes the machine precision parameter
        //   dpmeps as the smallest floating point number such that
        //   1 + dpmeps differs from 1.
        //
        //   This subroutine is based on the subroutine machar described in
        //
        //   W. J. Cody,
        //   MACHAR: A subroutine to dynamically determine machine parameters,
        //   ACM Transactions on Mathematical Software, 14, 1988, pages 303-311.
        //
        //   The subroutine statement is:
        //
        //     subroutine dpeps(dpmeps)
        //
        //   where
        //
        //     dpmeps is a double precision variable.
        //       On entry dpmeps need not be specified.
        //       On exit dpmeps is the machine precision.
        //
        //   MINPACK-2 Project. February 1991.
        //   Argonne National Laboratory and University of Minnesota.
        //   Brett M. Averick.
        //
        //   *******
        //function
        double dpmeps()
        {
                //int i ;
                long int ibeta,irnd,it,itemp,negep;
                long double a,b,beta,betain,betah,temp,tempa,temp1,
                        zero=0.0,one=1.0,two=2.0;
                long double ddpmeps;

                //determine ibeta, beta ala malcolm.

                a = one;
                b = one;
goto10:
                a = a + a;
                temp = a + one;
                temp1 = temp - a;
                if (temp1 - one == zero) goto goto10;

goto20:
                b = b + b;
                temp = a + b;
                itemp = int(temp - a);
                if (itemp == 0) goto goto20;
                ibeta = itemp;
                beta = double(ibeta);

                //determine it, irnd.

                it = 0;
                b = one;
goto30:
                it += 1;
                b *= beta;
                temp = b + one;
                temp1 = temp - b;
                if (temp1 - one == zero) goto goto30;
                irnd = 0;
                betah = beta/two;
                temp = a + betah;
                if (temp - a != zero) irnd = 1;
                tempa = a + beta;
                temp = tempa + betah;
                if ((irnd == 0) && (temp - tempa != zero)) irnd = 2;

                //determine ddpmeps.

                negep = it + 3;
                betain = one/beta;
                a = one;
                for (int i = 1; i <= negep; i++)
                        a *= betain;

goto50:
                temp = one + a;
                if (temp - one != zero) goto goto60;
                a *= beta;
                goto goto50;

goto60:
                ddpmeps = a;
                if ((ibeta == 2) || (irnd == 0)) goto goto70;
                a = (a*(one + a))/two;
                temp = one + a;
                if (temp - one != zero) ddpmeps = a;

goto70:
                return ddpmeps;
        }//dpmeps()


        //   constant times a vector plus a vector.
        //   uses unrolled loops for increments equal to one.
        //   jack dongarra, linpack, 3/11/78.
        //function
        void daxpy(const int& n, const double& da, const double* const & dx, 
                const int& incx, double* const & dy, const int& incy)
        {
                //int i;
                int ix,iy,m,mp1;
                if (n <= 0) return;
                if (da == 0.0) return;
                if(incx==1 && incy==1) goto goto20;

                //code for unequal increments or equal increments
                //not equal to 1

                ix = 1;
                iy = 1;
                if(incx < 0) ix = (-n+1)*incx + 1;
                if(incy < 0) iy = (-n+1)*incy + 1;

                for (int i = 1; i <= n; i++)
                {    
                        dy[iy] += da*dx[ix];
                        ix += incx;
                        iy += incy;
                }
                return;

                //code for both increments equal to 1
                //clean-up loop

goto20:
                m = n%4;
                if (m == 0) goto goto40;

                for (int i = 1; i <= m; i++)
                        dy[i] += da*dx[i];

                if (n < 4) return;

goto40:
                mp1 = m + 1;

                for (int i = mp1; i <= n; i+=4)
                {
                        dy[i] += da*dx[i];
                        dy[i + 1] += da*dx[i + 1];
                        dy[i + 2] += da*dx[i + 2];
                        dy[i + 3] += da*dx[i + 3];
                }
        }//daxpy()


        //   copies a vector, x, to a vector, y.
        //   uses unrolled loops for increments equal to one.
        //   jack dongarra, linpack, 3/11/78.
        //function  
        void dcopy(const int& n, const double* const & dx, const int& incx, 
                double* const & dy,const int& incy)
        {
                //int i;
                int ix,iy,m,mp1;

                if (n <= 0) return;
                if(incx==1 && incy==1) goto goto20;

                //code for unequal increments or equal increments
                //not equal to 1

                ix = 1;
                iy = 1;
                if (incx < 0) ix = (-n+1)*incx + 1;
                if (incy < 0) iy = (-n+1)*incy + 1;
                for (int i = 1; i <= n; i++)
                {
                        dy[iy] = dx[ix];
                        ix += incx;
                        iy += incy;
                }
                return;

                //code for both increments equal to 1

                //clean-up loop

goto20:
                m = n%7;
                if (m == 0) goto goto40;
                for (int i = 1; i <= m; i++)
                        dy[i] = dx[i];

                if (n < 7) return;

goto40:
                mp1 = m + 1;
                for (int i = mp1; i <= n; i+=7)
                {
                        dy[i] = dx[i];
                        dy[i + 1] = dx[i + 1];
                        dy[i + 2] = dx[i + 2];
                        dy[i + 3] = dx[i + 3];
                        dy[i + 4] = dx[i + 4];
                        dy[i + 5] = dx[i + 5];
                        dy[i + 6] = dx[i + 6];
                }      
        }//dcopy()


        //forms the dot product of two vectors.
        //uses unrolled loops for increments equal to one.
        //jack dongarra, linpack, 3/11/78.
        //function
        double ddot(const int& n, const double* const & dx, const int& incx, 
                const double* const & dy, const int& incy)
        {
                double dtemp;
                //int i;
                int ix,iy,m,mp1;

                double dddot = 0.0;
                dtemp = 0.0;
                if (n <= 0) return dddot;
                if (incx==1 && incy==1) goto goto20;

                //code for unequal increments or equal increments
                //not equal to 1

                ix = 1;
                iy = 1;
                if (incx < 0) ix = (-n+1)*incx + 1;
                if (incy < 0) iy = (-n+1)*incy + 1;

                for (int i = 1; i <= n; i++)
                {
                        dtemp += dx[ix]*dy[iy];
                        ix += incx;
                        iy += incy;
                }
                dddot = dtemp;
                return dddot;

                //code for both increments equal to 1

                //clean-up loop

goto20:
                m = n%5;
                if (m == 0) goto goto40;

                for (int i = 1; i <= m; i++)
                        dtemp += dx[i]*dy[i];

                if( n < 5 ) goto goto60;

goto40:
                mp1 = m + 1;


                for (int i = mp1; i <= n; i+=5)
                {
                        dtemp += dx[i]*dy[i] + dx[i + 1]*dy[i + 1] +
                                dx[i + 2]*dy[i + 2] + dx[i + 3]*dy[i + 3] + dx[i + 4]*dy[i + 4];
                }

goto60:
                dddot = dtemp;
                return dddot;
        }//ddot()


        //   dpofa factors a double precision symmetric positive definite
        //   matrix.
        //
        //   dpofa is usually called by dpoco, but it can be called
        //   directly with a saving in time if  rcond  is not needed.
        //   (time for dpoco) = (1 + 18/n)*(time for dpofa) .
        //
        //   on entry
        //
        //      a       double precision(lda, n)
        //              the symmetric matrix to be factored.  only the
        //              diagonal and upper triangle are used.
        //
        //      lda     integer
        //              the leading dimension of the array  a .
        //
        //      n       integer
        //              the order of the matrix  a .
        //
        //   on return
        //
        //      a       an upper triangular matrix  r  so that  a = trans(r)*r
        //              where  trans(r)  is the transpose.
        //              the strict lower triangle is unaltered.
        //              if  info != 0 , the factorization is not complete.
        //
        //      info    integer
        //              = 0  for normal return.
        //              = k  signals an error condition.  the leading minor
        //                   of order  k  is not positive definite.
        //
        //   linpack.  this version dated 08/14/78 .
        //   cleve moler, university of new mexico, argonne national lab.
        //
        //   subroutines and functions
        //
        //   blas ddot
        //   fortran sqrt
        //
        //   internal variables
        //  
        //function
        void dpofa(double* const & a, const int& lda, const int& n, int& info,
                bool smallerMatrix=false)
        {
                double t;
                double s;
                //int j,k
                int jm1;
                int idx;
                //begin block with ...exits to 40

                for (int j = 1; j <= n; j++)
                {
                        info = j;
                        s = 0.0;
                        jm1 = j - 1;
                        if (jm1 < 1) goto goto20;
                        for (int k = 1; k <= jm1; k++)
                        {
                                idx = getIdx(k,j,lda);
                                t = a[idx] - ddot(k-1,&(a[getIdx(1,k,lda)-1]),1,
                                        &(a[getIdx(1,j,lda)-1]),1);
                                t = t/a[getIdx(k,k,lda)];
                                a[idx] = t;
                                s += t*t;
                        }

goto20:
                        idx = getIdx(j,j,lda);
                        s = a[idx] - s;

                        //......exit
                        if (s <= 0.0) goto goto40;
                        a[idx] = sqrt(s);
                }
                info = 0;
goto40:
                t = 0;
        }//dpofa()


        //scales a vector by a constant.
        //uses unrolled loops for increment equal to one.
        //jack dongarra, linpack, 3/11/78.
        //modified 3/93 to return if incx <= 0.
        //function    
        void dscal(const int& n, const double& da, double* const & dx, 
                const int& incx)
        {
                //int i;
                int m,mp1,nincx;

                if (n <= 0 || incx <= 0) return;
                if (incx == 1) goto goto20;

                //code for increment not equal to 1

                nincx = n*incx;

                for (int i = 1; i <= nincx; i+=incx)
                        dx[i] = da*dx[i];
                return;

                //code for increment equal to 1

                //clean-up loop

goto20:
                m = n%5;
                if (m == 0) goto goto40;

                for (int i = 1; i <= m; i++)
                        dx[i] = da*dx[i];

                if (n < 5) return;

goto40:
                mp1 = m + 1;
                for (int i = mp1; i <= n; i+=5)
                {
                        dx[i] = da*dx[i];
                        dx[i + 1] = da*dx[i + 1];
                        dx[i + 2] = da*dx[i + 2];
                        dx[i + 3] = da*dx[i + 3];
                        dx[i + 4] = da*dx[i + 4];
                }
        }//dscal()

        //
        //
        //   dtrsl solves systems of the form
        //
        //                 t * x = b
        //   or
        //                 trans(t) * x = b
        //
        //   where t is a triangular matrix of order n. here trans(t)
        //   denotes the transpose of the matrix t.
        //
        //   on entry
        //
        //       t         double precision(ldt,n)
        //                 t contains the matrix of the system. the zero
        //                 elements of the matrix are not referenced, and
        //                 the corresponding elements of the array can be
        //                 used to store other information.
        //
        //       ldt       integer
        //                 ldt is the leading dimension of the array t.
        //
        //       n         integer
        //                 n is the order of the system.
        //
        //       b         double precision(n).
        //                 b contains the right hand side of the system.
        //
        //       job       integer
        //                 job specifies what kind of system is to be solved.
        //                 if job is
        //
        //                      00   solve t*x=b, t lower triangular,
        //                      01   solve t*x=b, t upper triangular,
        //                      10   solve trans(t)*x=b, t lower triangular,
        //                      11   solve trans(t)*x=b, t upper triangular.
        //
        //   on return
        //
        //       b         b contains the solution, if info == 0.
        //                 otherwise b is unaltered.
        //
        //       info      integer
        //                 info contains zero if the system is nonsingular.
        //                 otherwise info contains the index of
        //                 the first zero diagonal element of t.
        //
        //   linpack. this version dated 08/14/78 .
        //   g. w. stewart, university of maryland, argonne national lab.
        //
        //   subroutines and functions
        //
        //   blas daxpy,ddot
        //   fortran mod
        //function
        void dtrsl(double* const & t, const int& ldt, const int& n, 
                double* const &b, const int& job, int& info)
        {
                double temp;
                int ccase,j;
                //int jj;

                //begin block permitting ...exits to 150

                //check for zero diagonal elements.

                for (info = 1; info <= n; info++)
                {
                        //......exit
                        if (t[getIdx(info,info,ldt)] == 0.0) goto goto150;
                }
                info = 0;

                //determine the task and go to it.

                ccase = 1;
                if (job%10 != 0) ccase = 2;
                if ((job%100)/10 != 0) ccase += 2;
                if (ccase==1) goto goto20;
                if (ccase==2) goto goto50;
                if (ccase==3) goto goto80;
                if (ccase==4) goto goto110;

                //solve t*x=b for t lower triangular

goto20:
                b[1] /= t[getIdx(1,1,ldt)];
                if (n < 2) goto goto40;
                for (int j = 2; j <= n; j++)
                {
                        temp = -b[j-1];
                        daxpy(n-j+1,temp,&(t[getIdx(j,j-1,ldt)-1]),1,&(b[j-1]),1);
                        b[j] /= t[getIdx(j,j,ldt)];
                }
goto40:
                goto goto140;

                //solve t*x=b for t upper triangular.

goto50:    
                b[n] /= t[getIdx(n,n,ldt)];
                if (n < 2) goto goto70;
                for (int jj = 2; jj <= n; jj++)
                {
                        j = n - jj + 1;
                        temp = -b[j+1];
                        daxpy(j,temp,&(t[getIdx(1,j+1,ldt)-1]),1,&(b[1-1]),1);
                        b[j] /= t[getIdx(j,j,ldt)];
                }

goto70:
                goto goto140;

                //solve trans(t)*x=b for t lower triangular.

goto80:
                b[n] /= t[getIdx(n,n,ldt)];
                if (n < 2) goto goto100;
                for (int jj = 2; jj <= n; jj++)
                {
                        j = n - jj + 1;
                        b[j] = b[j] - ddot(jj-1,&(t[getIdx(j+1,j,ldt)-1]),1,&(b[j+1-1]),1);
                        b[j] /= t[getIdx(j,j,ldt)];
                }
goto100:
                goto goto140;

                //solve trans(t)*x=b for t upper triangular.

goto110:
                b[1] /= t[getIdx(1,1,ldt)];
                if (n < 2) goto goto130;
                for (int j = 2; j <= n; j++)
                {
                        b[j] = b[j] - ddot(j-1,&(t[getIdx(1,j,ldt)-1]),1,&(b[1-1]),1);
                        b[j] /= t[getIdx(j,j,ldt)];
                }

goto130:

                if (false);

goto140:

                if (false);

goto150:

                if (false);

        }//dtrsl()

private:
        Timer timer_;

};
#endif

---