if(Double.isNaN(m_f)) throw new Exception("Objective function value is NaN!");
    public double[] lnsrch(double[] xold, double[] gradient, 			   double[] direct, double stpmax,			   boolean[] isFixed, double[][] nwsBounds,			   FastVector wsBdsIndx)	throws Exception {		int i, j, k,len=xold.length, 	    fixedOne=-1; // idx of variable to be fixed	double alam, alamin; // lambda to be found, and its lower bound		// For convergence and bound test	double temp,test,alpha=Double.POSITIVE_INFINITY,fold=m_f,sum; 		// For cubic interpolation	double a,alam2=0,b,disc=0,maxalam=1.0,rhs1,rhs2,tmplam;		double[] x = new double[len]; // New variable values		// Scale the step 	for (sum=0.0,i=0;i<len;i++){	    if(!isFixed[i]) // For fixed variables, direction = 0		sum += direct[i]*direct[i];	}		sum = Math.sqrt(sum);		if (m_Debug)	    System.err.println("fold:  "+Utils.doubleToString(fold,10,7)+"\n"+			       "sum:  "+Utils.doubleToString(sum,10,7)+"\n"+			       "stpmax:  "+Utils.doubleToString(stpmax,10,7));	if (sum > stpmax){	    for (i=0;i<len;i++)		if(!isFixed[i])		    direct[i] *= stpmax/sum;			}	else	    maxalam = stpmax/sum;		// Compute initial rate of decrease, g'*d 	m_Slope=0.0;	for (i=0;i<len;i++){	    x[i] = xold[i];	    if(!isFixed[i])		m_Slope += gradient[i]*direct[i];	}		if (m_Debug)	    System.err.print("slope:  " + Utils.doubleToString(m_Slope,10,7)+ "\n");		// Slope too small	if(Math.abs(m_Slope)<=m_Zero){	    if (m_Debug)		System.err.println("Gradient and direction orthogonal -- "+				   "Min. found with current fixed variables"+				   " (or all variables fixed). Try to release"+				   " some variables now.");	    return x;	}		// Err: slope > 0	if(m_Slope > m_Zero){	    if(m_Debug)		for(int h=0; h<x.length; h++)		    System.err.println(h+": isFixed="+isFixed[h]+", x="+				       x[h]+", grad="+gradient[h]+", direct="+				       direct[h]);	    throw new Exception("g'*p positive! -- Try to debug from here: line 263.");	}		// Compute LAMBDAmin and upper bound of lambda--alpha	test=0.0;	for(i=0;i<len;i++){	    	    if(!isFixed[i]){// No need for fixed variables		temp=Math.abs(direct[i])/Math.max(Math.abs(x[i]),1.0);		if (temp > test) test=temp;	    }	}		if(test>m_Zero) // Not converge	    alamin = m_TOLX/test;	else{	    if (m_Debug)		System.err.println("Zero directions for all free variables -- "+				   "Min. found with current fixed variables"+				   " (or all variables fixed). Try to release"+				   " some variables now.");	    return x;	}			// Check whether any non-working-set bounds are "binding"	for(i=0;i<len;i++){	    if(!isFixed[i]){// No need for fixed variables		double alpi;		if((direct[i]<-m_Epsilon) && !Double.isNaN(nwsBounds[0][i])){//Not feasible		    alpi = (nwsBounds[0][i]-xold[i])/direct[i];		    if(alpi <= m_Zero){ // Zero			if (m_Debug)			    System.err.println("Fix variable "+i+					       " to lower bound "+ nwsBounds[0][i]+					       " from value "+ xold[i]);			x[i] = nwsBounds[0][i];			isFixed[i]=true; // Fix this variable			alpha = 0.0;			nwsBounds[0][i]=Double.NaN; //Add cons. to working set			wsBdsIndx.addElement(new Integer(i));		    }		    else if(alpha > alpi){ // Fix one variable in one iteration			alpha = alpi;			fixedOne = i;		    }					}		else if((direct[i]>m_Epsilon) && !Double.isNaN(nwsBounds[1][i])){//Not feasible		    alpi = (nwsBounds[1][i]-xold[i])/direct[i];		    if(alpi <= m_Zero){ // Zero			if (m_Debug)			    System.err.println("Fix variable "+i+					       " to upper bound "+ nwsBounds[1][i]+					       " from value "+ xold[i]);			x[i] = nwsBounds[1][i];			isFixed[i]=true; // Fix this variable			alpha = 0.0;			nwsBounds[1][i]=Double.NaN; //Add cons. to working set			wsBdsIndx.addElement(new Integer(i));		    }		    else if(alpha > alpi){			alpha = alpi;			fixedOne = i;		    }					}					    }	}			if (m_Debug){	    System.err.println("alamin: " + Utils.doubleToString(alamin,10,7));	    System.err.println("alpha: " + Utils.doubleToString(alpha,10,7));	}		if(alpha <= m_Zero){ // Zero	   	    m_IsZeroStep = true;	    if (m_Debug)		System.err.println("Alpha too small, try again");	    return x;	}		alam = alpha; // Always try full feasible newton step 	if(alam > 1.0)	    alam = 1.0;		// Iteration of one newton step, if necessary, backtracking is done	double initF=fold, // Initial function value	    hi=alam, lo=alam, newSlope=0, fhi=m_f, flo=m_f;// Variables used for beta condition	double[] newGrad;  // Gradient on the new variable values		kloop:	for (k=0;;k++) {	    if(m_Debug)		System.err.println("\nIteration: " + k);	    	    for (i=0;i<len;i++){		if(!isFixed[i]){		    x[i] = xold[i]+alam*direct[i];  // Compute xnew		    if(!Double.isNaN(nwsBounds[0][i]) && (x[i]<nwsBounds[0][i])){    			x[i] = nwsBounds[0][i]; //Rounding error			    }		    else if(!Double.isNaN(nwsBounds[1][i]) && (x[i]>nwsBounds[1][i])){					x[i] = nwsBounds[1][i]; //Rounding error			    }		}	    }	    	    m_f = objectiveFunction(x);    // Compute fnew	    	    while(Double.isInfinite(m_f)){ // Avoid infinity		if(m_Debug)		    System.err.println("Too large m_f.  Shrink step by half.");		alam *= 0.5; // Shrink by half		if(alam <= m_Epsilon){		    if(m_Debug)			System.err.println("Wrong starting points, change them!");		    return x;		}				for (i=0;i<len;i++)		    if(!isFixed[i])			x[i] = xold[i]+alam*direct[i]; 				m_f = objectiveFunction(x); 		initF = Double.POSITIVE_INFINITY;	    }	    	    if(m_Debug) {		System.err.println("obj. function: " + 				 Utils.doubleToString(m_f, 10, 7));		System.err.println("threshold: " + 				 Utils.doubleToString(fold+m_ALF*alam*m_Slope,10,7));	    }	    	    if(m_f<=fold+m_ALF*alam*m_Slope){// Alpha condition: sufficient function decrease		if(m_Debug)				    System.err.println("Sufficient function decrease (alpha condition): "); 		newGrad = evaluateGradient(x);		for(newSlope=0.0,i=0; i<len; i++)		    if(!isFixed[i])			newSlope += newGrad[i]*direct[i];		if(newSlope >= m_BETA*m_Slope){ // Beta condition: ensure pos. defnty.			    if(m_Debug)					System.err.println("Increasing derivatives (beta condition): "); 			    if((fixedOne!=-1) && (alam>=alpha)){ // Has bounds and over			if(direct[fixedOne] > 0){			    x[fixedOne] = nwsBounds[1][fixedOne]; // Avoid rounding error			    nwsBounds[1][fixedOne]=Double.NaN; //Add cons. to working set			}			else{			    x[fixedOne] = nwsBounds[0][fixedOne]; // Avoid rounding error			    nwsBounds[0][fixedOne]=Double.NaN; //Add cons. to working set			}						if(m_Debug)			    System.err.println("Fix variable "					       +fixedOne+" to bound "+ x[fixedOne]+					       " from value "+ xold[fixedOne]);			isFixed[fixedOne]=true; // Fix the variable			wsBdsIndx.addElement(new Integer(fixedOne));		    }				    return x;		}		else if(k==0){ // First time: increase alam 		    // Search for the smallest value not complying with alpha condition		    double upper = Math.min(alpha,maxalam); 		    if(m_Debug)			System.err.println("Alpha condition holds, increase alpha... ");		    while(!((alam>=upper) || (m_f>fold+m_ALF*alam*m_Slope))){			lo = alam;			flo = m_f;			alam *= 2.0;			if(alam>=upper)			    alam=upper;			for (i=0;i<len;i++)			    if(!isFixed[i])				x[i] = xold[i]+alam*direct[i];			m_f = objectiveFunction(x);			newGrad = evaluateGradient(x);			for(newSlope=0.0,i=0; i<len; i++)			    if(!isFixed[i])				newSlope += newGrad[i]*direct[i];			if(newSlope >= m_BETA*m_Slope){			    if (m_Debug)						System.err.println("Increasing derivatives (beta condition): ");			    if((fixedOne!=-1) && (alam>=alpha)){ // Has bounds and over				if(direct[fixedOne] > 0){				    x[fixedOne] = nwsBounds[1][fixedOne]; // Avoid rounding error				    nwsBounds[1][fixedOne]=Double.NaN; //Add cons. to working set				}				else{				    x[fixedOne] = nwsBounds[0][fixedOne]; // Avoid rounding error				    nwsBounds[0][fixedOne]=Double.NaN; //Add cons. to working set				}								if(m_Debug)				    System.err.println("Fix variable "						       +fixedOne+" to bound "+ x[fixedOne]+						       " from value "+ xold[fixedOne]);				isFixed[fixedOne]=true; // Fix the variable				wsBdsIndx.addElement(new Integer(fixedOne));			    }		 				    			    return x;			}		    }		    hi = alam;		    fhi = m_f;					    break kloop;		}		else{		    if(m_Debug)			System.err.println("Alpha condition holds.");		    hi = alam2; lo = alam; flo = m_f;		    break kloop;		}		    	    }        	    else if (alam < alamin) { // No feasible lambda found       		if(initF<fold){ 		    alam = Math.min(1.0,alpha);		    for (i=0;i<len;i++)			if(!isFixed[i])			    x[i] = xold[i]+alam*direct[i]; //Still take Alpha		    		    if (m_Debug)			System.err.println("No feasible lambda: still take"+					   " alpha="+alam);		    		    if((fixedOne!=-1) && (alam>=alpha)){ // Has bounds and over			if(direct[fixedOne] > 0){			    x[fixedOne] = nwsBounds[1][fixedOne]; // Avoid rounding error			    nwsBounds[1][fixedOne]=Double.NaN; //Add cons. to working set			}			else{			    x[fixedOne] = nwsBounds[0][fixedOne]; // Avoid rounding error			    nwsBounds[0][fixedOne]=Double.NaN; //Add cons. to working set			}						if(m_Debug)			    System.err.println("Fix variable "					       +fixedOne+" to bound "+ x[fixedOne]+					       " from value "+ xold[fixedOne]);			isFixed[fixedOne]=true; // Fix the variable			wsBdsIndx.addElement(new Integer(fixedOne));		    }		 		    		}		else{   // Convergence on delta(x)		    for(i=0;i<len;i++) 			x[i]=xold[i];		    m_f=fold;		    if (m_Debug)			System.err.println("Cannot find feasible lambda"); 		}				return x; 	    }	    else { // Backtracking by polynomial interpolation		if(k==0){ // First time backtrack: quadratic interpolation		    if(!Double.isInfinite(initF))			initF = m_f;		    		    // lambda = -g'(0)/(2*g''(0))		    tmplam = -0.5*alam*m_Slope/((m_f-fold)/alam-m_Slope);		}		else {    // Subsequent backtrack: cubic interpolation 		    rhs1 = m_f-fold-alam*m_Slope;		    rhs2 = fhi-fold-alam2*m_Slope;		    a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);		    b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);		    if (a == 0.0) tmplam = -m_Slope/(2.0*b);		    else {			disc=b*b-3.0*a*m_Slope;			if (disc < 0.0) disc = 0.0;			tmplam=(-b+Math.sqrt(disc))/(3.0*a);		    }		    if (m_Debug)			System.err.print("newstuff: \n" + 					 "a:   " + Utils.doubleToString(a,10,7)+ "\n" +					 "b:   " + Utils.doubleToString(b,10,7)+ "\n" +					 "disc:   " + Utils.doubleToString(disc,10,7)+ "\n" +					 "tmplam:   " + tmplam + "\n" +					 "alam:   " + Utils.doubleToString(alam,10,7)+ "\n");			    if (tmplam>0.5*alam)			tmplam=0.5*alam;             // lambda <= 0.5*lambda_old		}	    }	    alam2=alam;	    fhi=m_f;	    alam=Math.max(tmplam,0.1*alam);          // lambda >= 0.1*lambda_old	    	    if(alam>alpha){		throw new Exception("Sth. wrong in lnsrch:"+				    "Lambda infeasible!(lambda="+alam+				    ", alpha="+alpha+", upper="+tmplam+				    "|"+(-alpha*m_Slope/(2.0*((m_f-fold)/alpha-m_Slope)))+				    ", m_f="+m_f+", fold="+fold+				    ", slope="+m_Slope);	    }	    	} // Endfor(k=0;;k++)		// Quadratic interpolation between lamda values between lo and hi.	// If cannot find a value satisfying beta condition, use lo.	double ldiff = hi-lo, lincr;	if(m_Debug)	    System.err.println("Last stage of searching for beta condition (alam between "			       +Utils.doubleToString(lo,10,7)+" and "			       +Utils.doubleToString(hi,10,7)+")...");	while((newSlope<m_BETA*m_Slope) && (ldiff>=alamin)){	    lincr = -0.5*newSlope*ldiff*ldiff/(fhi-flo-newSlope*ldiff);	    if(lincr<0.2*ldiff) lincr=0.2*ldiff;	    alam = lo+lincr;	    for (i=0;i<len;i++)		if(!isFixed[i])		    x[i] = xold[i]+alam*direct[i];	    m_f = objectiveFunction(x);	    if(m_f>fold+m_ALF*alam*m_Slope){ 		// Alpha condition fails, shrink lambda_upper		ldiff = lincr;		fhi = m_f;	    }	    	    else{ // Alpha condition holds	    		newGrad = evaluateGradient(x);		for(newSlope=0.0,i=0; i<len; i++)		    if(!isFixed[i])			newSlope += newGrad[i]*direct[i];				if(newSlope < m_BETA*m_Slope){		    // Beta condition fails, shrink lambda_lower		    lo = alam;		    ldiff -= lincr;		    flo = m_f;		}	    }	} 		if(newSlope < m_BETA*m_Slope){ // Cannot satisfy beta condition, take lo	    if(m_Debug)		System.err.println("Beta condition cannot be satisfied, take alpha condition");	    alam=lo;	    for (i=0;i<len;i++)		if(!isFixed[i])		    x[i] = xold[i]+alam*direct[i];	    m_f = flo;	}	else if(m_Debug)	    System.err.println("Both alpha and beta conditions are satisfied. alam="			       +Utils.doubleToString(alam,10,7));		if((fixedOne!=-1) && (alam>=alpha)){ // Has bounds and over	    if(direct[fixedOne] > 0){		x[fixedOne] = nwsBounds[1][fixedOne]; // Avoid rounding error		nwsBounds[1][fixedOne]=Double.NaN; //Add cons. to working set	    }	    else{		x[fixedOne] = nwsBounds[0][fixedOne]; // Avoid rounding error		nwsBounds[0][fixedOne]=Double.NaN; //Add cons. to working set	    }	    	    if(m_Debug)		System.err.println("Fix variable "				   +fixedOne+" to bound "+ x[fixedOne]+				   " from value "+ xold[fixedOne]);	    isFixed[fixedOne]=true; // Fix the variable	    wsBdsIndx.addElement(new Integer(fixedOne));	}		return x;    }