Permalink
Browse files

New initial solution routine.

  • Loading branch information...
brianborchers committed Jun 10, 2018
1 parent f7d51dd commit 9ab14dca290997c3e37897b4472e19b02eb2b77c
Showing with 273 additions and 46 deletions.
  1. +273 −46 lib/initsoln.c
@@ -17,16 +17,18 @@ void initsoln(n,k,C,a,constraints,pX0,py0,pZ0)
double **py0;
struct blockmatrix *pZ0;
{
int i,j;
double alpha,beta;
double maxnrmA,nrmA;
double nrmC;
double scale1, scale2;
int i,j,blk;
double xi,eta;
double minxi,mineta;
double *normsofA;
double normofC;
double rescale;
struct sparseblock *ptr;

/*
* First, allocate storage for X0, y0, and Z0.
*/

alloc_mat(C,pX0);
alloc_mat(C,pZ0);
*py0=(double *)malloc(sizeof(double)*(k+1));
@@ -38,74 +40,299 @@ void initsoln(n,k,C,a,constraints,pX0,py0,pZ0)
};

/*
* next, pick alpha=n*max_i((1+|a_i|)/(1+||A_i||)).
* Go ahead and set y to 0.
*/

maxnrmA=0.0;
alpha=0.0;
for (i=1; i<=k; i++)
(*py0)[i]=0.0;

/*
* Allocate space to store the norms of the constraint matrices.
*/

normsofA=(double *)malloc(sizeof(double)*(k+1));
if (normsofA == NULL)
{
nrmA=0.0;
ptr=constraints[i].blocks;
printf("Storage allocation failed!\n");
exit(205);
};

/*
* Keep track of the minimum values of xi and eta.
*/

minxi=1.0e300;
mineta=1.0e300;

/*
* Next, work through the blocks one by one.
*/

while (ptr != NULL)
for (blk=1; blk<=C.nblocks; blk++)
{
/*
* Work out the norms of the objective and constraint blocks.
*/

switch (C.blocks[blk].blockcategory)
{
for (j=1; j<=ptr->numentries; j++)
case DIAG:
/*
* Get the norm of this block of C.
*/

normofC=0;
for (i=1; i<=C.blocks[blk].blocksize; i++)
normofC=normofC+C.blocks[blk].data.vec[i]*C.blocks[blk].data.vec[i];
normofC=sqrt(normofC);

/*
* Get the norms of this block in constraints 1, 2, ...
*/

for (i=1; i<=k; i++)
{
normsofA[i]=0.0;
ptr=constraints[i].blocks;
while (ptr != NULL)
{
if (ptr->blocknum == blk)
{
for (j=1; j<=ptr->numentries; j++)
{
normsofA[i]=normsofA[i]+(ptr->entries[j])*(ptr->entries[j]);
};
normsofA[i]=sqrt(normsofA[i]);
break;
};
ptr=ptr->next;
};
};
break;
case MATRIX:
/*
* Get the norm of this block of C.
*/

normofC=0;
for (j=1; j<=C.blocks[blk].blocksize; j++)
for (i=1; i<=C.blocks[blk].blocksize; i++)
normofC=normofC+C.blocks[blk].data.mat[ijtok(i,j,C.blocks[blk].blocksize)]*C.blocks[blk].data.mat[ijtok(i,j,C.blocks[blk].blocksize)];
normofC=sqrt(normofC);

/*
* Get the norms of this block in constraints 1, 2, ...
*/

for (i=1; i<=k; i++)
{
nrmA += (ptr->entries[j])*(ptr->entries[j]);
if (ptr->iindices[j] != ptr->jindices[j])
nrmA += (ptr->entries[j])*(ptr->entries[j]);
normsofA[i]=0.0;
ptr=constraints[i].blocks;
while (ptr != NULL)
{
if (ptr->blocknum == blk)
{
for (j=1; j<=ptr->numentries; j++)
{
if (ptr->iindices[j] == ptr->jindices[j])
{
normsofA[i]=normsofA[i]+(ptr->entries[j])*(ptr->entries[j]);
}
else
{
normsofA[i]=normsofA[i]+2*(ptr->entries[j])*(ptr->entries[j]);
};
};
normsofA[i]=sqrt(normsofA[i]);
break;
};
ptr=ptr->next;
};
};
ptr=ptr->next;
break;
case PACKEDMATRIX:
default:
printf("initsoln illegal block type \n");
exit(206);
};

/*
* Now, setup X0 and Z0 for this block.
*/

switch (C.blocks[blk].blockcategory)
{
case DIAG:
/*
* First, deal with the X0 block.
*/

xi=10;

if (sqrt(C.blocks[blk].blocksize) > xi)
xi=sqrt(C.blocks[blk].blocksize);

for (i=1; i<=k; i++)
{
if (sqrt(C.blocks[blk].blocksize)*(1.0+fabs(a[i]))/(1+normsofA[i]) > xi)
xi=sqrt(C.blocks[blk].blocksize)*(1.0+fabs(a[i]))/(1+normsofA[i]);
};

nrmA=sqrt(nrmA);
if (nrmA > maxnrmA)
maxnrmA=nrmA;
if (0==1)
printf("block=%d, xi=%e\n",blk,xi);

for (i=1; i<=C.blocks[blk].blocksize; i++)
(*pX0).blocks[blk].data.vec[i]=xi;

if (xi < minxi)
minxi=xi;

/*
* Next, deal with the Z0 block.
*/

eta=10;

if (sqrt(C.blocks[blk].blocksize) > eta)
eta=sqrt(C.blocks[blk].blocksize);

if (normofC > eta)
eta=normofC;

for (i=1; i<=k; i++)
{
if (normsofA[i] > eta)
eta=normsofA[i];
};

if ((1+fabs(a[i]))/(1+nrmA) > alpha)
alpha=(1+fabs(a[i]))/(1+nrmA);
if (0==1)
printf("block=%d, eta=%e\n",blk,eta);

for (i=1; i<=C.blocks[blk].blocksize; i++)
(*pZ0).blocks[blk].data.vec[i]=eta;

if (eta < mineta)
mineta=eta;

/*
* Done with current block.
*/
break;
case MATRIX:
/*
* First, deal with the X0 block.
*/

xi=10;

if (sqrt(C.blocks[blk].blocksize) > xi)
xi=sqrt(C.blocks[blk].blocksize);

for (i=1; i<=k; i++)
{
if (C.blocks[blk].blocksize*(1.0+fabs(a[i]))/(1+normsofA[i]) > xi)
xi=C.blocks[blk].blocksize*(1.0+fabs(a[i]))/(1+normsofA[i]);
};

if (0==1)
printf("block=%d, xi=%e\n",blk,xi);

for (j=1; j<=C.blocks[blk].blocksize; j++)
for (i=1; i<=C.blocks[blk].blocksize; i++)
(*pX0).blocks[blk].data.mat[ijtok(i,j,C.blocks[blk].blocksize)]=0.0;

for (i=1; i<=C.blocks[blk].blocksize; i++)
(*pX0).blocks[blk].data.mat[ijtok(i,i,C.blocks[blk].blocksize)]=xi;

if (xi < minxi)
minxi=xi;

/*
* Next, deal with the Z0 block.
*/

eta=10;

if (sqrt(C.blocks[blk].blocksize) > eta)
eta=sqrt(C.blocks[blk].blocksize);

if (normofC > eta)
eta=normofC;

for (i=1; i<=k; i++)
{
if (normsofA[i] > eta)
eta=normsofA[i];
};

if (0==1)
printf("block=%d, eta=%e\n",blk,eta);

for (j=1; j<=C.blocks[blk].blocksize; j++)
for (i=1; i<=C.blocks[blk].blocksize; i++)
(*pZ0).blocks[blk].data.mat[ijtok(i,j,C.blocks[blk].blocksize)]=0.0;

for (i=1; i<=C.blocks[blk].blocksize; i++)
(*pZ0).blocks[blk].data.mat[ijtok(i,i,C.blocks[blk].blocksize)]=eta;

if (eta < mineta)
mineta=eta;

/*
* Done with current block.
*/
break;
case PACKEDMATRIX:
default:
printf("initsoln illegal block type \n");
exit(206);
};
};

alpha=n*alpha;
/*
* Look for situations where X0 or Z0 is far too large and scale down
* as needed.
*/

/*
* Next, calculate the F norm of C.
* First, look at A(X0) versus the right hand side. We'll reuse
* normsofA as a work vector. If norm(A(X0) is 10 times bigger than
* norm(a), then scale it back down, but keep all the diagonal elements
* at least 10.
*/

nrmC=Fnorm(C);
op_a(k,constraints,*pX0,normsofA);

if (nrmC > maxnrmA)
{
beta=(1+nrmC)/sqrt(n*1.0);
}
else
if (norm2(k,normsofA+1) > 10.0*norm2(k,a+1))
{
beta=(1+maxnrmA)/sqrt(n*1.0);
rescale=10*norm2(k,a+1)/norm2(k,normsofA+1);
if (rescale*minxi < 10)
rescale=10/minxi;
if (0==1)
printf("Rescaling X0 by factor of %e\n",rescale);
scalemat(rescale,*pX0,*pX0);
};

/*
* Now that we have alpha and beta, make X=100*alpha*I, Z=100*beta*I, y=0.
* Next, look at norm(C) versus norm(Z0) and rescale Z0 as with X0.
*/

make_i(*pX0);
scale1=0.0;
scale2=10*alpha;
mat_mult(scale1,scale2,*pX0,*pX0,*pX0);

make_i(*pZ0);
scale1=0.0;
scale2=10*beta;
mat_mult(scale1,scale2,*pZ0,*pZ0,*pZ0);

if (Fnorm(*pZ0) > 10*Fnorm(C))
{
rescale=10*Fnorm(C)/Fnorm(*pZ0);
if (rescale*mineta < 10)
rescale=10/mineta;
if (0==1)
printf("rescaling Z0 by %e\n",rescale);
scalemat(rescale,*pZ0,*pZ0);
};


/*
* Set y to 0.
* Free up working storage.
*/
for (i=1; i<=k; i++)
(*py0)[i]=0;

free(normsofA);

}



0 comments on commit 9ab14dc

Please sign in to comment.