lswitchrk_matrix(ig,n,b1,(long)1); 
      i=lin_dependent[1];  
      lin_dependent[1]=lin_dependent[b1];
      lin_dependent[b1]=i;
    }
    else {
      lswitchrk_matrix(ig,n,b2,(long)1); 
      i=lin_dependent[1];  
      lin_dependent[1]=lin_dependent[b2];
      lin_dependent[b2]=i;
    }
    lswitchrk_matrix(ig,n,b1,(long)0); 
    i=lin_dependent[0];  
    lin_dependent[0]=lin_dependent[b1];
    lin_dependent[b1]=i;
  }
  /* lprint_matrix(d,n); */
  /* lprint_matrix(ig,n); */

  lcopy_matrix(g,n,g_new);   /* restore g_new matrix */
  if(add>0)
    for(j=0;j<n;j++) {
      for(k=0;k<n;k++) {
	g_new[j*n+k]+=add*ce[j]*ce[k];
      }
    }

  for(i=0;i<n;i++) {  /* fix linear dependent vectors */
    g0_new[i]=g0[i]+add*ce0[0]*ce[i];
  }
  if(m>0) ce0_new[0]=-ce0[0];
  for(i=0;i<n;i++) {  /* fix linear dependent vectors */
    if(lin_dependent[i]) {
      for(j=0;j<n;j++) {
	if(!lin_dependent[j]) {
	  g0_new[j]+=start[i]*g_new[i*n+j];
	}
      }
      if(m>0) ce0_new[0]-=(start[i]*ce[i]);
    }
  }
  from=0;   /* remove linear dependent vectors */
  to=0;
  n_indep=0;
  for(i=0;i<n;i++) {
    if(!lin_dependent[i]) {
      g0_new[n_indep]=g0_new[i];
      ce_new[n_indep]=ce[i]; 
      low_new[n_indep]=low[i];
      up_new[n_indep]=up[i];
      primal[n_indep]=start[i];
      n_indep++;
    }
    for(j=0;j<n;j++) {
      if((!lin_dependent[i]) && (!lin_dependent[j])) {
        ig[to]=ig[from];
        g_new[to]=g_new[from];
	to++;
      }
      from++;
    }
  }

  if(verbosity>=3) {
    printf("real_qp_size(%ld)...",n_indep);
  }
  
  /* cannot optimize with only one variable */
  if((n_indep<=1) && (m>0) && (!changed)) { 
    for(i=n-1;i>=0;i--) {
      primal[i]=init[i];
    }
    return((int)ONLY_ONE_VARIABLE);
  }

  if((!changed) || (n_indep>1)) { 
    result=solve_dual(n_indep,m,precision,epsilon_crit,maxiter,g_new,g0_new,
		      ce_new,ce0_new,low_new,up_new,primal,d,d0,ig,
		      dual,dual_old,temp,goal);
  }
  else {
    result=PRIMAL_OPTIMAL;
  }
  
  j=n_indep;
  for(i=n-1;i>=0;i--) {
    if(!lin_dependent[i]) {
      j--;
      primal[i]=primal[j];
    }
    else {
      primal[i]=start[i];  /* leave as is */
    }
    temp[i]=primal[i];
  }
   
  obj_before=calculate_qp_objective(n,g,g0,init);
  obj_after=calculate_qp_objective(n,g,g0,primal);
  (*progress)=obj_before-obj_after;
  if(verbosity>=3) {
    printf("before(%.30f)...after(%.30f)...result_sd(%d)...",
	   obj_before,obj_after,result); 
  }

  return((int)result);
}


int solve_dual(n,m,precision,epsilon_crit,maxiter,g,g0,ce,ce0,low,up,primal,
	       d,d0,ig,dual,dual_old,temp,goal)
     /* Solves the dual using the method of Hildreth and D'Espo. */
     /* Can only handle problems with zero or exactly one */
     /* equality constraints. */

     long   n;            /* number of variables */
     long   m;            /* number of linear equality constraints */
     double precision;    /* solve at least to this dual precision */
     double epsilon_crit; /* stop, if KT-Conditions approx fulfilled */
     long   maxiter;      /* stop after that many iterations */
     double *g;
     double *g0;          /* linear part of objective */
     double *ce,*ce0;     /* linear equality constraints */
     double *low,*up;     /* box constraints */
     double *primal;      /* variables (with initial values) */
     double *d,*d0,*ig,*dual,*dual_old,*temp;       /* buffer  */
     long goal;
{
  long i,j,k,iter;
  double sum,w,maxviol,viol,temp1,temp2,isnantest;
  double model_b,dist;
  long retrain,maxfaktor,primal_optimal=0,at_bound,scalemaxiter;
  double epsilon_a=1E-15,epsilon_hideo;
  double eq; 

  if((m<0) || (m>1)) 
    perror("SOLVE DUAL: inappropriate number of eq-constrains!");

  /*  
  printf("\n");
  for(i=0;i<n;i++) {
    printf("%f: ",g0[i]);
    for(j=0;j<n;j++) {
      printf("%f ",g[i*n+j]);
    }
    printf(": a=%.30f",primal[i]);
    printf(": y=%f\n",ce[i]);
  }
  */

  for(i=0;i<2*(n+m);i++) {
    dual[i]=0;
    dual_old[i]=0;
  }
  for(i=0;i<n;i++) {   
    for(j=0;j<n;j++) {   /* dual hessian for box constraints */
      d[i*2*(n+m)+j]=ig[i*n+j];
      d[(i+n)*2*(n+m)+j]=-ig[i*n+j];
      d[i*2*(n+m)+j+n]=-ig[i*n+j];
      d[(i+n)*2*(n+m)+j+n]=ig[i*n+j];
    }
    if(m>0) {
      sum=0;              /* dual hessian for eq constraints */
      for(j=0;j<n;j++) {
	sum+=(ce[j]*ig[i*n+j]);
      }
      d[i*2*(n+m)+2*n]=sum;
      d[i*2*(n+m)+2*n+1]=-sum;
      d[(n+i)*2*(n+m)+2*n]=-sum;
      d[(n+i)*2*(n+m)+2*n+1]=sum;
      d[(n+n)*2*(n+m)+i]=sum;
      d[(n+n+1)*2*(n+m)+i]=-sum;
      d[(n+n)*2*(n+m)+(n+i)]=-sum;
      d[(n+n+1)*2*(n+m)+(n+i)]=sum;
      
      sum=0;
      for(j=0;j<n;j++) {
	for(k=0;k<n;k++) {
	  sum+=(ce[k]*ce[j]*ig[j*n+k]);
	}
      }
      d[(n+n)*2*(n+m)+2*n]=sum;
      d[(n+n)*2*(n+m)+2*n+1]=-sum;
      d[(n+n+1)*2*(n+m)+2*n]=-sum;
      d[(n+n+1)*2*(n+m)+2*n+1]=sum;
    } 
  }

  for(i=0;i<n;i++) {   /* dual linear component for the box constraints */
    w=0;
    for(j=0;j<n;j++) {
      w+=(ig[i*n+j]*g0[j]); 
    }
    d0[i]=up[i]+w;
    d0[i+n]=-low[i]-w;
  }

  if(m>0) {  
    sum=0;             /* dual linear component for eq constraints */
    for(j=0;j<n;j++) {
      for(k=0;k<n;k++) {
	sum+=(ce[k]*ig[k*n+j]*g0[j]); 
      }
    }
    d0[2*n]=ce0[0]+sum;
    d0[2*n+1]=-ce0[0]-sum;
  }

  maxviol=999999;
  iter=0;
  retrain=1;
  maxfaktor=1;
  scalemaxiter=maxiter/5;
  while((retrain) && (maxviol > 0) && (iter < (scalemaxiter*maxfaktor))) {
    iter++;
    
    while((maxviol > precision) && (iter < (scalemaxiter*maxfaktor))) {
      iter++;
      maxviol=0;
      for(i=0;i<2*(n+m);i++) {
	sum=d0[i];
	for(j=0;j<2*(n+m);j++) {
	  sum+=d[i*2*(n+m)+j]*dual_old[j];
	}
	sum-=d[i*2*(n+m)+i]*dual_old[i];
	dual[i]=-sum/d[i*2*(n+m)+i];
	if(dual[i]<0) dual[i]=0;
	
	viol=fabs(dual[i]-dual_old[i]);
	if(viol>maxviol) 
	  maxviol=viol;
	dual_old[i]=dual[i];
      }
      /*
      printf("%d) maxviol=%20f precision=%f\n",iter,maxviol,precision); 
      */
    }
  
    if(m>0) {
      for(i=0;i<n;i++) {
	temp[i]=dual[i]-dual[i+n]+ce[i]*(dual[n+n]-dual[n+n+1])+g0[i];
      }
    } 
    else {
      for(i=0;i<n;i++) {
	temp[i]=dual[i]-dual[i+n]+g0[i];
      }
    }
    for(i=0;i<n;i++) {
      primal[i]=0;             /* calc value of primal variables */
      for(j=0;j<n;j++) {
	primal[i]+=ig[i*n+j]*temp[j];
      }
      primal[i]*=-1.0;
      if(primal[i]<=(low[i])) {  /* clip conservatively */
	primal[i]=low[i];
      }
      else if(primal[i]>=(up[i])) {
	primal[i]=up[i];
      }
    }

    if(m>0) 
      model_b=dual[n+n+1]-dual[n+n];
    else
      model_b=0;

    epsilon_hideo=EPSILON_HIDEO;
    for(i=0;i<n;i++) {           /* check precision of alphas */
      dist=-model_b*ce[i]; 
      dist+=(g0[i]+1.0);
      for(j=0;j<i;j++) {
	dist+=(primal[j]*g[j*n+i]);
      }
      for(j=i;j<n;j++) {
	dist+=(primal[j]*g[i*n+j]);
      }
      if((primal[i]<(up[i]-epsilon_hideo)) && (dist < (1.0-epsilon_crit))) {
	epsilon_hideo=(up[i]-primal[i])*2.0;
      }
      else if((primal[i]>(low[i]+epsilon_hideo)) &&(dist>(1.0+epsilon_crit))) {
	epsilon_hideo=(primal[i]-low[i])*2.0;
      }
    }
    /* printf("\nEPSILON_HIDEO=%.30f\n",epsilon_hideo); */

    for(i=0;i<n;i++) {           /* clip alphas to bounds */
      if(primal[i]<=(low[i]+epsilon_hideo)) {
	primal[i]=low[i];
      }
      else if(primal[i]>=(up[i]-epsilon_hideo)) {
	primal[i]=up[i];
      }
    }

    retrain=0;
    primal_optimal=1;
    at_bound=0;
    for(i=0;(i<n);i++) {  /* check primal KT-Conditions */
      dist=-model_b*ce[i]; 
      dist+=(g0[i]+1.0);
      for(j=0;j<i;j++) {
	dist+=(primal[j]*g[j*n+i]);
      }
      for(j=i;j<n;j++) {
	dist+=(primal[j]*g[i*n+j]);
      }
      if((primal[i]<(up[i]-epsilon_a)) && (dist < (1.0-epsilon_crit))) {
	retrain=1;
	primal_optimal=0;
      }
      else if((primal[i]>(low[i]+epsilon_a)) && (dist > (1.0+epsilon_crit))) {
	retrain=1;
	primal_optimal=0;
      }
      if((primal[i]<=(low[i]+epsilon_a)) || (primal[i]>=(up[i]-epsilon_a))) {
	at_bound++;
      }
      /*    printf("HIDEOtemp: a[%ld]=%.30f, dist=%.6f, b=%f, at_bound=%ld\n",i,primal[i],dist,model_b,at_bound);  */
    }
    if(m>0) {
      eq=-ce0[0];               /* check precision of eq-constraint */
      for(i=0;i<n;i++) { 
	eq+=(ce[i]*primal[i]);
      }
      if((EPSILON_EQ < fabs(eq)) 
	 /*
	 && !((goal==PRIMAL_OPTIMAL) 
	       && (at_bound==n)) */
	 ) {
	retrain=1;
	primal_optimal=0;
      }
      /* printf("\n eq=%.30f ce0=%f at-bound=%ld\n",eq,ce0[0],at_bound);  */
    }

    if(retrain) {
      precision/=10;
      if(((goal == PRIMAL_OPTIMAL) && (maxfaktor < 50000))
	 || (maxfaktor < 5)) {
	maxfaktor++;
      }
    }
  }

  if(!primal_optimal) {
    for(i=0;i<n;i++) {
      primal[i]=0;             /* calc value of primal variables */
      for(j=0;j<n;j++) {
	primal[i]+=ig[i*n+j]*temp[j];
      }
      primal[i]*=-1.0;
      if(primal[i]<=(low[i]+epsilon_a)) {  /* clip conservatively */
	primal[i]=low[i];
      }
      else if(primal[i]>=(up[i]-epsilon_a)) {
	primal[i]=up[i];
      }
    }
  }

  isnantest=0;
  for(i=0;i<n;i++) {           /* check for isnan */
    isnantest+=primal[i];
  }

  if(m>0) {
    temp1=dual[n+n+1];   /* copy the dual variables for the eq */
    temp2=dual[n+n];     /* constraints to a handier location */
    for(i=n+n+1;i>=2;i--) {
      dual[i]=dual[i-2];
    }
    dual[0]=temp2;
    dual[1]=temp1;
    isnantest+=temp1+temp2;
  }

  if(isnan(isnantest)) {
    return((int)NAN_SOLUTION);
  }
  else if(primal_optimal) {
    return((int)PRIMAL_OPTIMAL);
  }
  else if(maxviol == 0.0) {
    return((int)DUAL_OPTIMAL);
  }
  else {
    return((int)MAXITER_EXCEEDED);
  }
}


void linvert_matrix(matrix,depth,inverse,lindep_sensitivity,lin_dependent)
double *matrix;
long depth;
double *inverse,lindep_sensitivity;
long *lin_dependent;  /* indicates the active parts of matrix on 
			 input and output*/
{
  long i,j,k;
  double factor;

  for(i=0;i<depth;i++) {
    /*    lin_dependent[i]=0; */
    for(j=0;j<depth;j++) {
      inverse[i*depth+j]=0.0;
    }
    inverse[i*depth+i]=1.0;
  }
  for(i=0;i<depth;i++) {
    if(lin_dependent[i] || (fabs(matrix[i*depth+i])<lindep_sensitivity)) {
      lin_dependent[i]=1;
    }
    else {
      for(j=i+1;j<depth;j++) {
	factor=matrix[j*depth+i]/matrix[i*depth+i];
	for(k=i;k<depth;k++) {
	  matrix[j*depth+k]-=(factor*matrix[i*depth+k]);
	}
	for(k=0;k<depth;k++) {
	  inverse[j*depth+k]-=(factor*inverse[i*depth+k]);
	}
      }
    }
  }
  for(i=depth-1;i>=0;i--) {
    if(!lin_dependent[i]) {
      factor=1/matrix[i*depth+i];
      for(k=0;k<depth;k++) {
	inverse[i*depth+k]*=factor;
      }
      matrix[i*depth+i]=1;
      for(j=i-1;j>=0;j--) {
	factor=matrix[j*depth+i];
	matrix[j*depth+i]=0;
	for(k=0;k<depth;k++) {
	  inverse[j*depth+k]-=(factor*inverse[i*depth+k]);
	}
      }
    }
  }
}

void lprint_matrix(matrix,depth)
double *matrix;
long depth;
{
  long i,j;
  for(i=0;i<depth;i++) {
    for(j=0;j<depth;j++) {
      printf("%5.2f ",(double)(matrix[i*depth+j]));
    }
    printf("\n");
  }
  printf("\n");
}

void ladd_matrix(matrix,depth,scalar)
double *matrix;
long depth;
double scalar;
{
  long i,j;
  for(i=0;i<depth;i++) {
    for(j=0;j<depth;j++) {
      matrix[i*depth+j]+=scalar;
    }
  }
}

void lcopy_matrix(matrix,depth,matrix2) 
double *matrix;
long depth;
double *matrix2;
{
  long i;
  
  for(i=0;i<(depth)*(depth);i++) {
    matrix2[i]=matrix[i];
  }
}

void lswitch_rows_matrix(matrix,depth,r1,r2) 
double *matrix;
long depth,r1,r2;
{
  long i;
  double temp;

  for(i=0;i<depth;i++) {
    temp=matrix[r1*depth+i];
    matrix[r1*depth+i]=matrix[r2*depth+i];
    matrix[r2*depth+i]=temp;
  }
}

void lswitchrk_matrix(matrix,depth,rk1,rk2) 
double *matrix;
long depth,rk1,rk2;
{
  long i;
  double temp;

  for(i=0;i<depth;i++) {
    temp=matrix[rk1*depth+i];
    matrix[rk1*depth+i]=matrix[rk2*depth+i];
    matrix[rk2*depth+i]=temp;
  }
  for(i=0;i<depth;i++) {
    temp=matrix[i*depth+rk1];
    matrix[i*depth+rk1]=matrix[i*depth+rk2];
    matrix[i*depth+rk2]=temp;
  }
}

double calculate_qp_objective(opt_n,opt_g,opt_g0,alpha)
long opt_n;
double *opt_g,*opt_g0,*alpha;
{
  double obj;
  long i,j;
  obj=0;  /* calculate objective  */
  for(i=0;i<opt_n;i++) {
    obj+=(opt_g0[i]*alpha[i]);
    obj+=(0.5*alpha[i]*alpha[i]*opt_g[i*opt_n+i]);
    for(j=0;j<i;j++) {
      obj+=(alpha[j]*alpha[i]*opt_g[j*opt_n+i]);
    }
  }
  return(obj);
}