sp.c

Go to the documentation of this file.

/*
 * Copyright (c) 1994 by Lieven Vandenberghe and Stephen Boyd.
 * Permission to use, copy, modify, and distribute this software for 
 * any purpose without fee is hereby granted, provided that this entire 
 * notice is included in all copies of any software which is or includes
 * a copy or modification of this software and in all copies of the 
 * supporting documentation for such software.
 * This software is being provided "as is", without any express or 
 * implied warranty.  In particular, the authors do not make any
 * representation or warranty of any kind concerning the merchantability
 * of this software or its fitness for any particular purpose.
 */
 
#include <stdio.h> 
#include <math.h>
#include <string.h>
#include "spd.h"

#ifdef _MSC_VER 
extern void Scistring (char *str);
#endif

void cngrncb(itype,n,AP,B,CP,temp)
 int itype;
 int n;
 double *AP;
 double *B;
 double *CP;
 double *temp;

/* 
 * if itype = 1, computes C = B*A*B', otherwise, computes C = B'*A*B 
 * A and B are nxn with A symmetric.
 *
 * Arguments:
 * - itype  = 1: compute C = B*A*B'
 *          = any other integer: computes C = B'*A*B
 * - n      dimension of A and B
 * - AP     (input) double array of size n*(n+1)2;
 *          the lower triangle of A in packed storage
 * - B      (input) double array of size n*n;
 * - CP     (output) double array of size n*(n+1)/2; 
 *          the lower triangle of C in packed storage
 * - temp:  n-array, workspace
 */


{
 int j, pos, lngth = n*(n+1)/2;
 int int1=1;
 double dbl0=0.0, dbl1=1.0; 

 /* C := 0 */
 F2C(dscal)(&lngth, &dbl0, CP, &int1);

 if (itype == 1){

   for (j=0, pos=0;  j<n;  pos+=n-j, j++){

      /* temp = A*B(j,:)' */
      F2C(dspmv)("L", &n, &dbl1, AP, B+j, &n, &dbl0, temp, &int1);

      /* C(j:n,j) = B(j:n,:)*temp */
      lngth = n-j;
      F2C(dgemv)("N", &lngth, &n, &dbl1, B+j, &n, temp, &int1, &dbl0,
             CP+pos, &int1);

   }

 } else {
 
   for (j=0, pos=0;  j<n;  pos+=n-j, j++){

      /* temp = A*B(:,j) */
      F2C(dspmv)("L", &n, &dbl1, AP, B+j*n, &int1, &dbl0, temp, &int1);

      /* C(j:n,j) = B(:,j:n)'*temp */
      lngth = n-j;
      F2C(dgemv)("T", &n, &lngth, &dbl1, B+j*n, &n, temp, &int1, &dbl0,
             CP+pos, &int1);

   }
 } 
 
}


double inprd(X,Z,L,blck_szs)
 double *X;
 double *Z;
 int L;
 int *blck_szs;

/*
 * Computes Tr X*Z
 *
 * Arguments:
 * X,Z:       block diagonal matrices with L blocks X^0, ..., X^{L-1},
 *            and Z^0, ..., Z^{L-1}.  X^j and Z^j have size 
 *            blck_szs[j] times blck_szs[j].  Every block is stored 
 *            using packed storage of the lower triangle.
 * L:         number of blocks
 * blck_szs:  integer vector of length L 
 *            blck_szs[i], i=0,...,L-1 is the size of block i
 *
 */

{
 double result;
 int i, j, k, lngth, pos, sz, int1=1;
 
 /* sz = length of Z and X */  
 for (i=0, sz=0;  i<L;  i++)  sz += (blck_szs[i]*(blck_szs[i]+1))/2;

 /* result = Tr X Z + contributions of diagonal elements */
 result = 2.0*F2C(ddot)(&sz, X, &int1, Z, &int1);

 /* correct for diagonal elements 
  * loop over blocks, j=0,...,L-1  */
 for (j=0, pos=0;  j<L;  j++)

    /* loop over columns, k=0,...,blck_szs[j]-1 
     * pos is position of (k,k) element of block j 
     * lngth is length of column k */
    for (k=0, lngth=blck_szs[j];  k<blck_szs[j];  pos+=lngth, 
         lngth-=1, k++) 
    
       /* subtract Z^j_{kk}*X^j_{kk} from result */
       result -= Z[pos]*X[pos];

 return result;
}

int C2F(spf)(m,L,F,blck_szs,c,x,Z,ul,nu,abstol,reltol,tv,iters, 
             work,lwork,iwork,info)
 int *m;                /* no of variables */
 int *L;                /* no of blocks in F */
 double *F;            /* F_i's in packed storage */
 int *blck_szs;        /* L-vector, dimensions of diagonal blocks */
 double *c;            /* m-vector */
 double *x;            /* m-vector */
 double *Z;            /* block diagonal matrix in packed storage */
 double *ul;           /* ul[0] = pr. obj, ul[1] = du. obj */ 
 double *nu;            /* >= 1.0 */
 double *abstol;        /* absolute accuracy */
 double *reltol;        /* relative accuracy */
 double *tv;            /* target value */ 
 int *iters;           /* on entry: the maximum number of iterations,
                        * on exit: the number of iterations taken */
 double *work;         /* work array */
 int *lwork;            /* size of work */
 int *iwork;           /* work array of m integers */
 int *info;            /* status on termination */
{
 return(sp(*m,*L,F,blck_szs,c,x,Z,ul,*nu,*abstol,*reltol,*tv,iters, work,
           *lwork,iwork,info));
}

int sp(m,L,F,blck_szs,c,x,Z,ul,nu,abstol,reltol,tv,iters,work,
       lwork,iwork,info)

 int m;                /* no of variables */
 int L;                /* no of blocks in F */
 double *F;            /* F_i's in packed storage */
 int *blck_szs;        /* L-vector, dimensions of diagonal blocks */
 double *c;            /* m-vector */
 double *x;            /* m-vector */
 double *Z;            /* block diagonal matrix in packed storage */
 double *ul;           /* ul[0] = pr. obj, ul[1] = du. obj */ 
 double nu;            /* >= 1.0 */
 double abstol;        /* absolute accuracy */
 double reltol;        /* relative accuracy */
 double tv;            /* target value */ 
 int *iters;           /* on entry: the maximum number of iterations,
                        * on exit: the number of iterations taken */
 double *work;         /* work array */
 int lwork;            /* size of work */
 int *iwork;           /* work array of m integers */
 int *info;             /* status on termination */

/*
 * Solves semidefinite program 
 *
 *  minimize    c'*x 
 *  subject to  F_0 + x_1*F_1 + ... + x_m*F_m  >= 0
 * 
 * and its dual
 * 
 *  maximize    -Tr F_0*Z 
 *  subject to  Z >= 0
 *              Tr F_i*Z = c_i, i=1,...,m
 *
 *
 * Convergence criteria: 
 * (1) maxiters is exceeded
 * (2) duality gap is less than abstol 
 * (3) primal and dual objective are both positive and 
 *     duality gap is less than reltol * dual objective         
 *     or primal and dual objective are both negative and
 *     duality gap is less than reltol * minus the primal objective
 * (4) reltol is negative and primal objective is less than tv 
 * (5) reltol is negative and dual objective is greater than tv
 * 
 * Arguments:
 * - m:        number of variables x_i. m >= 1.
 * - L:        number of diagonal blocks in F_i. L >= 1.
 * - F:        the block diagonal matrices F_i, i=0,...,m. 
 *             it is assumed that the matrices F_i are linearly 
 *             independent. 
 *             let F_i^j, i=0,..,m, j=0,...,L-1 denote the jth 
 *             diagonal block of F_i, 
 *             the array F contains F_0^0, ..., F_0^{L-1}, F_1^0, ..., 
 *             F_1^{L-1}, ..., F_m^0, ..., F_m^{L-1}, in this order, 
 *             using packed storage for the lower triangular part of 
 *             F_i^j.
 * - blck_szs: an integer L-vector. blck_szs[j], j=0,....L-1 gives the 
 *             size of block j, ie, F_i^j has size blck_szs[j] 
 *             times blck_szs[j].
 * - c:        m-vector, primal objective.
 * - x:        m-vector.  On entry, a strictly primal feasible point. 
 *             On exit, the last iterate for x.
 * - Z:        block diagonal matrix with L blocks Z^0, ..., Z^{L-1}.
 *             Z^j has size blck_szs[j] times blck_szs[j].
 *             Every block is stored using packed storage of the lower 
 *             triangular part.
 *             On entry, a strictly dual feasible point.  On exit, the 
 *             last dual iterate.
 * - ul:       two-vector.  On exit, ul[0] is the primal objective value
 *             c'*x;  ul[1] is the dual objective value -Tr F_0*Z.
 * - nu:       >= 1.0. Controls the rate of convergence.
 * - abstol:   absolute tolerance, >= MINABSTOL.
 * - reltol:   relative tolerance.  Has a special meaning when negative.
 * - tv:       target value, only referenced if reltol < 0.
 * - iters:    on entry: maximum number of iterations >= 0,
 *             on exit: the number of iterations taken.
 * - work:     work array of size lwork.
 * - lwork:    size of work, must be at least:
 *             (m+2)*sz + up_sz + 2*n + ltemp, with 
 *             ltemp = max( m+sz*nb, 3max_n + max_n*(max_n+1), 3*m )
 *             (sz: space needed to store one matrix F_i in packed
 *             storage, ie, 
 *                sum_{j=0}^{L-1} blck_szs[j]*(blck_szs[j]+1)/2;
 *             up_sz: space needed to store one matrix F_i in 
 *             unpacked storage, ie, 
 *                sum_{j=0}^{L-1} blck_szs[j]*blck_szs[j];
 *             max_n: max block size;
 *             n: sum of the block sizes.
 *             nb >= 1, for best performance, nb should be at least
 *             equal to the optimal block size for dgels.
 * - iwork:    work array of m integers
 * - info:     returns 1 if maxiters exceeded,  2 if absolute accuracy
 *             is reached, 3 if relative accuracy is reached,
 *             4 if target value is reached, 5 if target value is
 *             not achievable; 
 *             negative values indicate errors: -i means argument i 
 *             has an illegal value, -18 stands for all other errors.
 *
 *         
 * Returns 0 for normal exit, 1 if an error occurred.
 *
 */
 

{
 int i, j, k, n, sz, up_sz, max_n, lngth, pos, pos2, pos3, pos4, ltemp, 
     maxiters, info2, minlwork; 
 double q, *rhs, *Fsc, *R, *X, rho, *dx, *sigx, *sigz, *dZ, *temp, scal,
        scal2, XdZ, ZdX, alphax, alphaz, lambda_ls, gradx, hessx,
        gradz, hessz, dalphax, dalphaz, gap, newgap=0.0, newu=0.0, 
        newl=0.0, maxpossigx, minnegsigx, maxpossigz, minnegsigz, nrmc,
        nrmx, nrmz, nrmmax, rcond; 
 int int2=2, int1=1;  
 double dbl1=1.0, dbl0=0.0, sqrt2=sqrt(2.0);
 char str[100];
 double dbl_epsilon;

 if (m < 1){
    sprintf(str, "m must be at least one. \n");
    Scistring(str);
    *info = -1;
    return 1;
 }
 if (L < 1){
    sprintf(str, "L must be at least one. \n");
    Scistring(str);
    *info = -2;
    return 1;
 }
 for (i=0; i<L; i++) if (blck_szs[i] < 1){
    sprintf(str, "blck_szs[%d] must be at least one.\n", i);
    Scistring(str);
    *info = -4;
    return 1;
 }
 if (nu < 1.0){
    sprintf(str, "nu must be at least 1.0. \n");
    Scistring(str);
    *info = -9;
    return 1;
 }
 

 /*
  * calculate dimensions:
  * n:      total size of semidefinite program
  * sz:     length of one block-diagonal matrix in packed storage
  * up_sz:  length of one block-diagonal matrix in unpacked storage
  * max_n:  size of biggest block
  */

 for (i=0, n=0, sz=0, up_sz=0, max_n=0;  i<L;  i++){
    n     += blck_szs[i];
    sz    += blck_szs[i]*(blck_szs[i]+1)/2;
    up_sz += blck_szs[i]*blck_szs[i];
    max_n  = MAX(max_n, blck_szs[i]);
 } 
 if (m > sz){
     sprintf(str, "The matrices Fi, i=1,...,m are linearly dependent.\n");
     Scistring(str);
    *info = -3;  return 1;
 }

 q = (double)n + nu*sqrt((double)n); 


 /*
  * check if Tr Fi*Z = c_i, i=1,...,m
  */

 nrmc = F2C(dnrm2)(&m, c, &int1);
 for (i=0; i<m; i++) 
 if (fabs(inprd(F+(i+1)*sz, Z, L, blck_szs) - c[i]) > nrmc*TOLC){
     sprintf(str, "Z0 does not satisfy equality conditions\
 for dual feasibility.\n");
     Scistring(str);
    *info = -7;
    return 1;
 }


 /*
  * organize workspace
  *
  * work:  (m+2)*sz + up_sz + 2*n + ltemp
  * minimum ltemp: the maximum of 
  *         m+sz*nb, 3*max_n + max_n*(max_n+1), and 3*m  
  *         (nb is at least one)
  * 
  * for dgels:        m + sz*nb, nb at least 1 
  * for dspev("N"):   3*max_n + max_n*(max_n+1)
  * for dspgv("N"):   3*max_n + max_n*(max_n+1)
  * for dspgv("V"):   3*max_n + max_n*(max_n+1)/2  
  * for cngrncb:      max_n
  * for dtrcon:       3*m
  * 
  * rhs  (sz):        work[0 ... sz-1] 
  * Fsc  (m*sz):      work[sz ... (m+1)*sz-1]
  * R    (up_sz):     work[(m+1)*sz ... (m+1)*sz+up_sz-1]
  * X    (sz):        work[(m+1)*sz+up_sz ... (m+2)*sz+up_sz-1]
  * sigx (n):         work[(m+2)*sz+up_sz ... (m+2)*sz+up_sz+n-1]
  * sigz (n):         work[(m+2)*sz+up_sz+n ... (m+2)*sz+up_sz+2*n-1]
  * temp (remainder): work[(m+2)*sz+up_sz+2*n ... lwork-1]
  */

 /* check lwork */
 minlwork = (m+2)*sz + up_sz + 2*n + 
            MAX( MAX( m+sz, 3*max_n + max_n*(max_n+1) ), 3*m ); 
 if (lwork < minlwork){
    sprintf(str, "Work space is too small.  Need at least\
 %d*sizeof(double).\n", minlwork);
    Scistring(str);
    *info = -15;
    return 1;
 } 

 rhs   = work;        /* rhs for ls problem */
 dx    = work;        /* solution of ls system; overlaps with rhs  */
 Fsc   = rhs + sz;    /* scaled matrices */
 dZ    = rhs + sz;    /* overlaps with first column of Fsc */
 R     = Fsc + m*sz;  /* eigenvectors of Z*F */
 X     = R + up_sz;   /* F(x) */
 sigx  = X + sz;      /* generalized eigenvalues of (dX,X) */
 sigz  = sigx + n;    /* generalized eigenvalues of (dZ,Z) */
 temp  = sigz + n;
 ltemp = lwork - (m+2)*sz - up_sz - 2*n; 


 maxiters = (*iters >= 0) ? *iters : MAXITERS;
 for (*iters=0; *iters <= maxiters; (*iters)++){


    /* compute F(x) = F_0 + x_1*F_1 + ... + x_m*F_m, store in X */
    F2C(dcopy)(&sz, F, &int1, X, &int1);
    F2C(dgemv)("N", &sz, &m, &dbl1, F+sz, &sz, x, &int1, &dbl1, X, &int1);


    /* 
     * compute generalized eigendecomp  Z*F*x = lambda*x
     * loop over blocks, i=0,...,L-1 
     * pos:  position of (0,0) element of block i in packed storage
     * pos2: position of (0,0) element of block i in unpacked
     *       storage
     * pos3: position of first eigenvalue of block i in sigx
     */

    for (i=0, pos=0, pos2=0, pos3=0, gap=0.0;  i<L; 
         pos += blck_szs[i]*(blck_szs[i]+1)/2, 
         pos2 += blck_szs[i]*blck_szs[i], 
         pos3 += blck_szs[i], i++){

       lngth = blck_szs[i]*(blck_szs[i]+1)/2;

       /* copy block i of Z in temp (need max_n*(max_n+1)/2) */
       F2C(dcopy)(&lngth, Z+pos, &int1, temp, &int1); 

       /* generalized eigenvalue decomposition Z*F*x = lambda*x
        * - eigenvectors V are normalized s.t. V^T*F*V = I
        * - store block i of V in R+pos2
        * - store eigenvalues of block i in sigx+pos3
        * - dspgv replaces X+pos by cholesky factor L of ith 
        *   block of F (F = L*L^T) 
        * use temp+lngth as workspace (need at least 3*max_n) */
       F2C(dspgv)(&int2, "V", "L", blck_szs+i, temp, X+pos, sigx+pos3, 
              R+pos2, blck_szs+i, temp+lngth, &info2);
       if (info2){
          sprintf(str,"Error in dspgv, info = %d.\n", info2);
          Scistring(str);
          if (*iters == 0 && info2 > blck_szs[i]){
             sprintf(str, "x0 is not strictly primal feasible.\n");
             Scistring(str);
             *info = -6;
          } else *info = -18;  
          return 1;
       }

       /* - replace sigx+pos3 by lambda^(1/2)
        * - normalize block i of V (stored in R+pos2) s.t. 
        *   V^T*F*V = Lambda^(1/2) */
       for (k=0; k<blck_szs[i]; k++){
          scal = sigx[pos3+k];
          if (scal < 0.0){       
             if (*iters == 0){ 
                sprintf(str, "Z0 is not positive definite.\n");
                Scistring(str);
                *info = 7;
             } else {
                sprintf(str, "F(x)*Z has a negative eigenvalue.\n");
                Scistring(str);
                *info = -18;
             }
             return 1;
          }
          gap += scal;    /* duality gap is sum of eigenvalues of ZF */
          scal2 = sqrt(scal);
          scal = sqrt(scal2);
          sigx[pos3+k] = scal2; 
          F2C(dscal)(blck_szs+i, &scal, R+pos2+k*blck_szs[i], &int1);
       }

    }


    /* 
     * check convergence 
     */

    ul[1] = -inprd(F,Z,L,blck_szs);         /* -Tr F_0 Z */
    ul[0] = F2C(ddot)(&m, c, &int1, x, &int1);  /* c^T x */
    if (*iters == 0){
        sprintf(str,"\n    primal obj.  dual obj.  dual. gap  \n");
        Scistring(str);
    }
    sprintf(str,"% 13.2e % 12.2e %10.2e \n", ul[0], ul[1], gap);
    Scistring(str);
    if (gap <= MAX(abstol, MINABSTOL))  *info = 2;
    else if ( (ul[1] > 0.0 && gap <= reltol*ul[1]) ||
              (ul[0] < 0.0 && gap <= reltol*(-ul[0])) ) *info = 3;
    else if ( reltol < 0.0 && ul[0] <= tv ) *info = 4;
    else if ( reltol < 0.0 && ul[1] >= tv ) *info = 5;
    else if ( *iters == maxiters ) *info = 1;
    else *info = 0;
    if (*info) return 0; 



    /* 
     * compute scaled matrices F 
     */

    for (j=0, pos=0;  j<m;  j++) for (i=0, pos2=0;  i<L; 
         pos += blck_szs[i]*(blck_szs[i]+1)/2, 
         pos2 += blck_szs[i]*blck_szs[i], i++) {

       /* compute R' * Fj(i) * R, store in Fsc+pos */
       cngrncb(2, blck_szs[i], F+sz+pos, R+pos2, Fsc+pos, temp);

       /* correct diagonal elements */
       for (k=0, pos4=pos;  k<blck_szs[i];  pos4 += blck_szs[i]-k, k++)
          Fsc[pos4] /= sqrt2;

    }


    /* 
     * form rhs = Lambda^(-1/2) - (q/gap) * Lambda^(1/2) 
     */

    F2C(dscal)(&sz, &dbl0, rhs, &int1);    /* rhs := 0 */
    rho = -q/gap;
    for (i=0, pos=0, pos3=0;  i<L;  
         pos += blck_szs[i]*(blck_szs[i]+1)/2, 
         pos3 += blck_szs[i], i++)
       for (k=0, pos4=pos;  k<blck_szs[i];  pos4+=blck_szs[i]-k, k++){
          scal = sigx[pos3+k];
          rhs[pos4] = (1.0/scal + rho*scal)/sqrt2; 
    }


    /*
     * solve least-squares problem; need workspace of size m + nb*sz
     * - rhs is overwritten by dx
     * - in first iteration, estimate condition number of Fsc
     */
  
    F2C(dgels)("N", &sz, &m, &int1, Fsc, &sz, rhs, &sz, temp, &ltemp, 
           &info2);
    if (info2){
       sprintf(str,"Error in dgels, info = %d.\n", info2);
       Scistring(str);
       *info = -18; return 1;
    }

    if (*iters == 0){
       
       /* estimate the condition number in 1-norm of the R-factor of 
        * the qr-decomposition of Fsc (is stored in Fsc) 
        * need work space of size 3*m */
       F2C(dtrcon)("1", "U", "N", &m, Fsc, &sz, &rcond, temp, iwork, 
                &info2);
       if (info2 < 0){
          sprintf(str,"Error in dtrcon, info = %d.\n", info2);
          Scistring(str);
          *info = -18; return 1;
       }
       if (rcond < MINRCOND) {
          sprintf(str,"The matrices F_i, i=1,...,m are linearly\
 dependent (or the initial points are very badly conditioned).\n");
          Scistring(str);
          *info = -3; return 1;
       }

    }
    


    /*
     * - compute dZ = 
     *   R*((q/gap)*Lambda^(1/2) - Lambda^(-1/2) + R^T*dF*R )*R^T
     * - compute generalized eigenvalues of (dF, F), store in sigx
     * - compute generalized eigenvalues of (dZ, Z), store in sigz
     * 
     * loop over blocks i=0,...,L-1
     * pos:  position of (0,0) element of block i in packed storage
     * pos2: position of (0,0) element of block i in unpacked storage 
     * pos3: position of first eigenvalue of in sigx and sigz
     */

    for (i=0, pos=0, pos2=0, pos3=0;  i<L; 
         pos  += blck_szs[i]*(blck_szs[i]+1)/2, 
         pos2 += blck_szs[i]*blck_szs[i], 
         pos3 += blck_szs[i], i++){

       lngth = blck_szs[i]*(blck_szs[i]+1)/2;

       /* compute ith block of dF = \sum \delta x_i F_i, 
        * store in temp */
       F2C(dgemv)("N", &lngth, &m, &dbl1, F+sz+pos, &sz, dx, &int1, 
              &dbl0, temp, &int1);

       /* scale dF as R'*dF*R, store in temp + lngth */
       cngrncb(2, blck_szs[i], temp, R+pos2, temp+lngth, temp+2*lngth);

       /* add (q/gap)*Lambda^(1/2) - Lambda^(-1/2) */
       for (k=0, pos4=lngth;  k<blck_szs[i];  pos4+=blck_szs[i]-k, k++)
          temp[pos4] -= rho*sigx[pos3+k] + 1.0/sigx[pos3+k];

       /* replace dF in temp by L^{-1}*dF*L^{-T},
        * (L: cholesky factor of F, stored in X)
        * and compute eigenvalues of L^{-1}*dF*L^{-T}  */
       F2C(dspgst)(&int1, "L", blck_szs+i, temp, X+pos, &info2);
       if (info2){ 
          sprintf(str,"Error in dspst, info = %d.\n", info2);
          Scistring(str);
          *info = -18;  return 1; 
       }
       /* temp has to be of size max_n*(max_n+1)+3*max_n */
       F2C(dspev)("N", "L", blck_szs+i, temp, sigx+pos3, NULL, &int1,
              temp+2*lngth, &info2);
       if (info2){
           sprintf(str,"Error in dspev, info = %d.\n", info2);
           Scistring(str);
          *info = -18;  return 1;
       }

       /* dZ := R*((q/gap)*Lambda^(1/2) - Lambda^(-1/2) + R'*dF*R)*R' */
       cngrncb(1, blck_szs[i], temp+lngth, R+pos2, dZ+pos, 
               temp+2*lngth);

       /* copy ith block of dZ to temp */
       F2C(dcopy)(&lngth, dZ+pos, &int1, temp, &int1);

       /* copy ith block of Z to temp + lngth */
       F2C(dcopy)(&lngth, Z+pos, &int1, temp+lngth, &int1);

       /* sigz: generalized eigenvalues of (dZ,Z)
        * required size of temp: 3*max_n + max_n*(max_n+1) */
       F2C(dspgv)(&int1, "N", "L", blck_szs+i, temp, temp+lngth, sigz+pos3,
              NULL, &int1, temp+2*lngth, &info2);
       if (info2){
           sprintf(str,"Error in dspgv, info = %d.\n", info2);Scistring(str);
          *info = -18;  return 1; 
       }

    }
    

    /* 
     * compute feasible rectangle for plane search
     */ 
 
    maxpossigx = 0.0;  minnegsigx = 0.0;
    maxpossigz = 0.0;  minnegsigz = 0.0;
    for (i=0; i<n; i++) {
       if ( sigx[i] > maxpossigx ) 
          maxpossigx = sigx[i];  /* max pos eigenvalue in sigx */
       else if ( sigx[i] < minnegsigx ) 
          minnegsigx = sigx[i];  /* min neg eigenvalue in sigx */
       if ( sigz[i] > maxpossigz ) 
          maxpossigz = sigz[i];  /* max pos eigenvalue in sigz */
       else if ( sigz[i] < minnegsigz ) 
          minnegsigz = sigz[i];  /* min neg eigenvalue in sigz */
    }
    nrmx = F2C(dnrm2)(&n, sigx, &int1);        /* norm of scaled dx */ 
    nrmz = F2C(dnrm2)(&n, sigz, &int1);        /* norm of scaled dZ */
    nrmmax = MAX( nrmx, nrmz);

    XdZ = inprd(F,dZ,L,blck_szs);          /* Tr F0*dZ */ 
    ZdX = F2C(ddot)(&m, c, &int1, dx, &int1);  /* c^T*dx */ 


    /*
     * check corners of feasible rectangle
     */

   dbl_epsilon = F2C(dlamch)("e"); 
   if (nrmx > SIGTOL*nrmmax)
      if (ZdX < 0.0) 
          alphax = (minnegsigx < -dbl_epsilon) ? -1.0/minnegsigx : 0.0;
      else 
          alphax = (maxpossigx >  dbl_epsilon) ? -1.0/maxpossigx : 0.0;
    else alphax = 0.0;
    
    if (nrmz > SIGTOL*nrmmax)
       if (XdZ < 0.0)
          alphaz = (minnegsigz < -dbl_epsilon) ? -1.0/minnegsigz : 0.0;
       else 
          alphaz = (maxpossigz >  dbl_epsilon) ? -1.0/maxpossigz : 0.0;
    else alphaz = 0.0;

    newgap = gap + alphax*ZdX + alphaz*XdZ;
    newu = ul[0] + alphax*ZdX;
    newl = ul[1] - alphaz*XdZ;

    if (newgap <= MAX(abstol, MINABSTOL))  *info = 2;
    else if ( (newl > 0.0 && newgap <= reltol*newl) ||
              (newu < 0.0 && newgap <= -reltol*newu) ) *info = 3;
    else if ( reltol < 0.0 && newu <= tv ) *info = 4;
    else if ( reltol < 0.0 && newl >= tv ) *info = 5;
    else if ( *iters == maxiters ) *info = 1;
    else *info = 0;
    
    if (*info) {   
       F2C(daxpy)(&m, &alphax, dx, &int1, x, &int1);
       F2C(daxpy)(&sz, &alphaz, dZ, &int1, Z, &int1);
       gap = newgap;  ul[0] = newu;   ul[1] = newl;
       sprintf(str,"% 13.2e % 12.2e %10.2e \n", ul[0], ul[1], gap);
       Scistring(str);
       (*iters)++;
       return 0;
    }


    /*
     * plane search 
     *  minimize   phi(alphax,alphaz) = 
     *    q*log(dual_gap + alphax*c^T*dx + alphaz* Tr F_0 dZ)
     *  - sum log (1+alphax*sigx_i) - sum log (1+alphaz*sigz)
     */

    alphax = 0.0;  alphaz = 0.0;  lambda_ls = 1.0;

    if (nrmx > SIGTOL*nrmmax)
       if (nrmz > SIGTOL*nrmmax)    /* compute primal and dual steps */
          while ( lambda_ls > 1e-4 ) {

             /* compute 1st and 2nd derivatives of phi */
             rho = q/(gap + alphax*ZdX + alphaz*XdZ);
             gradx = rho*ZdX;  hessx = 0.0;
             gradz = rho*XdZ;  hessz = 0.0;
             for (i=0; i<n; i++){
                gradx -= sigx[i] / (1.0+alphax*sigx[i]);
                hessx += SQR( sigx[i] / (1.0+alphax*sigx[i]) );
                gradz -= sigz[i] / (1.0+alphaz*sigz[i]);
                hessz += SQR( sigz[i] / (1.0+alphaz*sigz[i]) );
             }

             /* newton step */
             dalphax = -gradx/hessx;  dalphaz = -gradz/hessz;
             lambda_ls = sqrt( SQR(gradx)/hessx + SQR(gradz)/hessz ); 
             alphax += (lambda_ls > 0.25) ? 
                       dalphax/(1.0+lambda_ls) : dalphax;
             alphaz += (lambda_ls > 0.25) ? 
                       dalphaz/(1.0+lambda_ls) : dalphaz;

         }
         
       else while ( lambda_ls > 1e-4 ) {  /* primal step only */

             /* compute 1st and 2nd derivatives of phi */
             rho = q/(gap + alphax*ZdX);
             gradx = rho*ZdX;  hessx = 0.0;
             for (i=0; i<n; i++){
                gradx -= sigx[i] / (1.0+alphax*sigx[i]);
                hessx += SQR( sigx[i] / (1.0+alphax*sigx[i]) );
             }

             /* newton step */
             dalphax = -gradx/hessx;
             lambda_ls = fabs(gradx)/sqrt(hessx);
             alphax += (lambda_ls > 0.25) ? 
                       dalphax/(1.0+lambda_ls) : dalphax;

       }

       else if (nrmz > SIGTOL*nrmmax)        /* dual step only */
          while ( lambda_ls > 1e-4 ) {

             /* compute 1st and 2nd derivatives of phi */
             rho = q/(gap + alphaz*XdZ);
             gradz = rho*XdZ;  hessz = 0.0;
             for (i=0; i<n; i++){
                gradz -= sigz[i] / (1.0+alphaz*sigz[i]);
                hessz += SQR( sigz[i] / (1.0+alphaz*sigz[i]) );
             }
       
             /* newton step */
             dalphaz = -gradz/hessz;
             lambda_ls = fabs(gradz)/sqrt(hessz);
             alphaz += (lambda_ls > 0.25) ? 
                       dalphaz/(1.0+lambda_ls) : dalphaz;
        }
    


    /* update x and Z */
    F2C(daxpy)(&m, &alphax, dx, &int1, x, &int1);
    F2C(daxpy)(&sz, &alphaz, dZ, &int1, Z, &int1);

 }
 
 return -1;   /* should never happen */
}

---