{                    if ( iter > m ) bound = m;                    ys = ddot ( n , w , iypt + npt , 1 , w , ispt + npt , 1 );                    if ( ! diagco )                    {                        yy = ddot ( n , w , iypt + npt , 1 , w , iypt + npt , 1 );                        for ( i = 1 ; i <= n ; i += 1 )                        {                            diag [ i -1] = ys / yy;                        }                    }                    else                    {                        iflag[0]=2;                        return;                    }                }            }            if ( execute_entire_while_loop || iflag[0] == 2 )            {                if ( iter != 1 )                {                    if ( diagco )                    {                        for ( i = 1 ; i <= n ; i += 1 )                        {                            if ( diag [ i -1] <= 0 )                            {                                iflag[0]=-2;                                throw new ExceptionWithIflag( iflag[0], "The "+i+"-th diagonal element of the inverse hessian approximation is not positive." );                            }                        }                    }                    cp= point;                    if ( point == 0 ) cp = m;                    w [ n + cp -1] = 1 / ys;                    for ( i = 1 ; i <= n ; i += 1 )                    {                        w [ i -1] = - g [ i -1];                    }                    cp= point;                    for ( i = 1 ; i <= bound ; i += 1 )                    {                        cp=cp-1;                        if ( cp == - 1 ) cp = m - 1;                        sq = ddot ( n , w , ispt + cp * n , 1 , w , 0 , 1 );                        inmc=n+m+cp+1;                        iycn=iypt+cp*n;                        w [ inmc -1] = w [ n + cp + 1 -1] * sq;                        daxpy ( n , - w [ inmc -1] , w , iycn , 1 , w , 0 , 1 );                    }                    for ( i = 1 ; i <= n ; i += 1 )                    {                        w [ i -1] = diag [ i -1] * w [ i -1];                    }                    for ( i = 1 ; i <= bound ; i += 1 )                    {                        yr = ddot ( n , w , iypt + cp * n , 1 , w , 0 , 1 );                        beta = w [ n + cp + 1 -1] * yr;                        inmc=n+m+cp+1;                        beta = w [ inmc -1] - beta;                        iscn=ispt+cp*n;                        daxpy ( n , beta , w , iscn , 1 , w , 0 , 1 );                        cp=cp+1;                        if ( cp == m ) cp = 0;                    }                    for ( i = 1 ; i <= n ; i += 1 )                    {                        w [ ispt + point * n + i -1] = w [ i -1];                    }                }                nfev[0]=0;                stp[0]=1;                if ( iter == 1 ) stp[0] = stp1;                for ( i = 1 ; i <= n ; i += 1 )                {                    w [ i -1] = g [ i -1];                }            }            Mcsrch.mcsrch ( n , x , f , g , w , ispt + point * n , stp , ftol , xtol , maxfev , info , nfev , diag );            if ( info[0] == - 1 )            {                iflag[0]=1;                return;            }            if ( info[0] != 1 )            {                iflag[0]=-1;                throw new ExceptionWithIflag( iflag[0], "Line search failed. See documentation of routine mcsrch. Error return of line search: info = "+info[0]+" Possible causes: function or gradient are incorrect, or incorrect tolerances." );            }            nfun= nfun + nfev[0];            npt=point*n;            for ( i = 1 ; i <= n ; i += 1 )            {                w [ ispt + npt + i -1] = stp[0] * w [ ispt + npt + i -1];                w [ iypt + npt + i -1] = g [ i -1] - w [ i -1];            }            point=point+1;            if ( point == m ) point = 0;            gnorm = Math.sqrt ( ddot ( n , g , 0 , 1 , g , 0 , 1 ) );            xnorm = Math.sqrt ( ddot ( n , x , 0 , 1 , x , 0 , 1 ) );            xnorm = Math.max ( 1.0 , xnorm );            if ( gnorm / xnorm <= eps ) finish = true;            if ( iprint [ 1 -1] >= 0 ) lb1 ( iprint , iter , nfun , gnorm , n , m , x , f , g , stp , finish );            // Cache the current solution vector. Due to the spaghetti-like            // nature of this code, it's not possible to quit here and return;            // we need to go back to the top of the loop, and eventually call            // mcsrch one more time -- but that will modify the solution vector.            // So we need to keep a copy of the solution vector as it was at            // the completion (info[0]==1) of the most recent line search.            System.arraycopy( x, 0, solution_cache, 0, n );            if ( finish )            {                iflag[0]=0;                    return;            }            execute_entire_while_loop = true;		// from now on, execute whole loop        }    }    /** Print debugging and status messages for <code>lbfgs</code>.      * Depending on the parameter <code>iprint</code>, this can include      * number of function evaluations, current function value, etc.      * The messages are output to <code>System.err</code>.      *      * @param iprint Specifies output generated by <code>lbfgs</code>.<p>      *		<code>iprint[0]</code> specifies the frequency of the output:      *		<ul>      *		<li> <code>iprint[0] &lt; 0</code>: no output is generated,      *		<li> <code>iprint[0] = 0</code>: output only at first and last iteration,      *		<li> <code>iprint[0] &gt; 0</code>: output every <code>iprint[0]</code> iterations.      *		</ul><p>      *      *		<code>iprint[1]</code> specifies the type of output generated:      *		<ul>      *		<li> <code>iprint[1] = 0</code>: iteration count, number of function      *			evaluations, function value, norm of the gradient, and steplength,      *		<li> <code>iprint[1] = 1</code>: same as <code>iprint[1]=0</code>, plus vector of      *			variables and  gradient vector at the initial point,      *		<li> <code>iprint[1] = 2</code>: same as <code>iprint[1]=1</code>, plus vector of      *			variables,      *		<li> <code>iprint[1] = 3</code>: same as <code>iprint[1]=2</code>, plus gradient vector.      *		</ul>      * @param iter Number of iterations so far.      * @param nfun Number of function evaluations so far.      * @param gnorm Norm of gradient at current solution <code>x</code>.      * @param n Number of free parameters.      * @param m Number of corrections kept.      * @param x Current solution.      * @param f Function value at current solution.      * @param g Gradient at current solution <code>x</code>.      * @param stp Current stepsize.      * @param finish Whether this method should print the ``we're done'' message.      */    public static void lb1 ( int[] iprint , int iter , int nfun , double gnorm , int n , int m , double[] x , double f , double[] g , double[] stp , boolean finish )    {        int i;        if ( iter == 0 )        {            System.err.println( "*************************************************" );            System.err.println( "  n = " + n + "   number of corrections = " + m + "\n       initial values" );            System.err.println( " f =  " + f + "   gnorm =  " + gnorm );            if ( iprint [ 2 -1] >= 1 )            {                System.err.print( " vector x =" );                for ( i = 1; i <= n; i++ )                    System.err.print( "  "+x[i-1] );                System.err.println( "" );                System.err.print( " gradient vector g =" );                for ( i = 1; i <= n; i++ )                    System.err.print( "  "+g[i-1] );                System.err.println( "" );            }            System.err.println( "*************************************************" );            System.err.println( "\ti\tnfn\tfunc\tgnorm\tsteplength" );        }        else        {            if ( ( iprint [ 1 -1] == 0 ) && ( iter != 1 && ! finish ) ) return;            if ( iprint [ 1 -1] != 0 )            {                if ( (iter - 1) % iprint [ 1 -1] == 0 || finish )                {                    if ( iprint [ 2 -1] > 1 && iter > 1 )                        System.err.println( "\ti\tnfn\tfunc\tgnorm\tsteplength" );                    System.err.println( "\t"+iter+"\t"+nfun+"\t"+f+"\t"+gnorm+"\t"+stp[0] );                }                else                {                    return;                }            }            else            {                if ( iprint [ 2 -1] > 1 && finish )                        System.err.println( "\ti\tnfn\tfunc\tgnorm\tsteplength" );                System.err.println( "\t"+iter+"\t"+nfun+"\t"+f+"\t"+gnorm+"\t"+stp[0] );            }            if ( iprint [ 2 -1] == 2 || iprint [ 2 -1] == 3 )            {                if ( finish )                {                    System.err.print( " final point x =" );                }                else                {                    System.err.print( " vector x =  " );                }                for ( i = 1; i <= n; i++ )                    System.err.print( "  "+x[i-1] );                System.err.println( "" );                if ( iprint [ 2 -1] == 3 )                {                    System.err.print( " gradient vector g =" );                    for ( i = 1; i <= n; i++ )                        System.err.print( "  "+g[i-1] );                    System.err.println( "" );                }            }            if ( finish )                System.err.println( " The minimization terminated without detecting errors. iflag = 0" );        }        return;    }    /** Compute the sum of a vector times a scalara plus another vector.      * Adapted from the subroutine <code>daxpy</code> in <code>lbfgs.f</code>.      * There could well be faster ways to carry out this operation; this      * code is a straight translation from the Fortran.      */    public static void daxpy ( int n , double da , double[] dx , int ix0, int incx , double[] dy , int iy0, int incy )    {        int i, ix, iy, m, mp1;        if ( n <= 0 ) return;        if ( da == 0 ) return;        if  ( ! ( incx == 1 && incy == 1 ) )        {            ix = 1;            iy = 1;            if ( incx < 0 ) ix = ( - n + 1 ) * incx + 1;            if ( incy < 0 ) iy = ( - n + 1 ) * incy + 1;            for ( i = 1 ; i <= n ; i += 1 )            {                dy [ iy0+iy -1] = dy [ iy0+iy -1] + da * dx [ ix0+ix -1];                ix = ix + incx;                iy = iy + incy;            }            return;        }        m = n % 4;        if ( m != 0 )        {            for ( i = 1 ; i <= m ; i += 1 )            {                dy [ iy0+i -1] = dy [ iy0+i -1] + da * dx [ ix0+i -1];            }            if ( n < 4 ) return;        }        mp1 = m + 1;        for ( i = mp1 ; i <= n ; i += 4 )        {            dy [ iy0+i -1] = dy [ iy0+i -1] + da * dx [ ix0+i -1];            dy [ iy0+i + 1 -1] = dy [ iy0+i + 1 -1] + da * dx [ ix0+i + 1 -1];            dy [ iy0+i + 2 -1] = dy [ iy0+i + 2 -1] + da * dx [ ix0+i + 2 -1];            dy [ iy0+i + 3 -1] = dy [ iy0+i + 3 -1] + da * dx [ ix0+i + 3 -1];        }        return;    }    /** Compute the dot product of two vectors.      * Adapted from the subroutine <code>ddot</code> in <code>lbfgs.f</code>.      * There could well be faster ways to carry out this operation; this      * code is a straight translation from the Fortran.      */    public static double ddot ( int n, double[] dx, int ix0, int incx, double[] dy, int iy0, int incy )    {        double dtemp;        int i, ix, iy, m, mp1;        dtemp = 0;        if ( n <= 0 ) return 0;        if ( !( incx == 1 && incy == 1 ) )        {            ix = 1;            iy = 1;            if ( incx < 0 ) ix = ( - n + 1 ) * incx + 1;            if ( incy < 0 ) iy = ( - n + 1 ) * incy + 1;            for ( i = 1 ; i <= n ; i += 1 )            {                dtemp = dtemp + dx [ ix0+ix -1] * dy [ iy0+iy -1];                ix = ix + incx;                iy = iy + incy;            }            return dtemp;        }        m = n % 5;        if ( m != 0 )        {            for ( i = 1 ; i <= m ; i += 1 )            {                dtemp = dtemp + dx [ ix0+i -1] * dy [ iy0+i -1];            }            if ( n < 5 ) return dtemp;        }        mp1 = m + 1;        for ( i = mp1 ; i <= n ; i += 5 )        {            dtemp = dtemp + dx [ ix0+i -1] * dy [ iy0+i -1] + dx [ ix0+i + 1 -1] * dy [ iy0+i + 1 -1] + dx [ ix0+i + 2 -1] * dy [ iy0+i + 2 -1] + dx [ ix0+i + 3 -1] * dy [ iy0+i + 3 -1] + dx [ ix0+i + 4 -1] * dy [ iy0+i + 4 -1];        }        return dtemp;    }}