152 void Thomas(
double *zMax,
double *z0,
double * Av,
int nn){
185 for (i=1;i < nn; i++){
187 z0[i]=tt - tt / (i+2);
198 z_max= fabs(z0[nn-1]);
200 for (i=nn-2; i>=0; i--){
202 z0[i]+= z0[i+1] - z0[i+1]/ (i+2);
223 void Rose(
double *zMax,
double *z0,
double * Av,
int nn){
265 z0[nn-1]=Av[nn-1]- s;
266 for (i=nn-2;i >=0; i--){
267 z0[i]=Av[i] + z0[i+1];
274 for (i=0; i<nn; i++){
306 int supportSet(
double *x,
double *v,
double *z,
double *g,
int * S,
double lambda,
int nn){
308 int i, j, n=nn+1, numS=0;
318 if ( ( (z[i]==lambda) && (g[i] <
delta) ) || ( (z[i]==-lambda) && (g[i] >
delta) )){
350 for (i=0;i<=S[0]; i++)
353 temp=( temp + z[ S[0] ] ) / (S[0] +1);
354 for (i=0;i<=S[0]; i++)
362 for (j=1; j < numS; j++){
364 for (i= S[j-1] +1; i<= S[j]; i++){
370 temp=(temp - z[ S[j-1] ] + z[ S[j] ])/ (S[j]- S[j-1]);
372 for (i= S[j-1] +1; i<= S[j]; i++){
381 for (i=S[numS-1] +1 ;i< n; i++)
385 temp=( temp - z[ S[numS-1] ] ) / (nn - S[numS-1]);
387 for (i=S[numS-1] +1 ;i< n; i++)
414 void dualityGap(
double *gap,
double *z,
double *g,
double *s,
double *Av,
double lambda,
int nn){
420 g[0]=z[0] + z[0] - z[1] - Av[0];
421 for (i=1;i<nn-1;i++){
422 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
424 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
443 temp=temp + s[i] *g[i]
458 void dualityGap2(
double *gap,
double *z,
double *g,
double *s,
double *Av,
double lambda,
int nn){
489 temp=temp + s[i] *g[i]
507 double *v,
double *Av,
508 double *g,
double *s,
int *S,
509 double lambda,
int nn){
511 int i, m, numS, n=nn+1;
512 double temp, funVal1, funVal2;
532 temp+=z[i]*(g[i] + Av[i]);
535 temp=temp + z[i] *(g[i] + Av[i])
536 + z[i+1]*(g[i+1] + Av[i+1])
537 + z[i+2]*(g[i+2] + Av[i+2])
538 + z[i+3]*(g[i+3] + Av[i+3])
539 + z[i+4]*(g[i+4] + Av[i+4]);
549 temp=temp + fabs(g[i])
554 funVal1=funVal1+ temp*lambda;
571 temp+=(x[i]-v[i]) * (x[i]-v[i]);
574 temp=temp + (x[i]- v[i]) * ( x[i]- v[i])
575 + (x[i+1]-v[i+1]) * (x[i+1]-v[i+1])
576 + (x[i+2]-v[i+2]) * (x[i+2]-v[i+2])
577 + (x[i+3]-v[i+3]) * (x[i+3]-v[i+3])
578 + (x[i+4]-v[i+4]) * (x[i+4]-v[i+4]);
585 temp+=fabs( x[i+1]-x[i] );
588 temp=temp + fabs( x[i+1]-x[i] )
589 + fabs( x[i+2]-x[i+1] )
590 + fabs( x[i+3]-x[i+2] )
591 + fabs( x[i+4]-x[i+3] )
592 + fabs( x[i+5]-x[i+4] );
593 funVal2=funVal2 + lambda * temp;
603 if (funVal2 > funVal1){
608 x[i]= v[i] - z[i-1] + z[i];
609 x[n-1]=v[n-1] - z[n-2];
619 *gap=*gap - (funVal1- funVal2);
629 double lambda,
int nn)
634 int* S=(
int *) malloc(
sizeof(
int)*nn);
635 double *x=(
double *)malloc(
sizeof(
double)*n);
636 double *s=(
double *)malloc(
sizeof(
double)*nn);
637 double *Av=(
double *)malloc(
sizeof(
double)*nn);
655 g[0]=z[0] + z[0] - z[1] - Av[0];
656 for (i=1;i<nn-1;i++){
657 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
659 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
674 s[i]=Av[i] - x[i+1] + x[i];
846 int sfa(
double *x,
double *gap,
int * activeS,
847 double *z,
double *z0,
double * v,
double * Av,
848 double lambda,
int nn,
int maxStep,
849 double *s,
double *g,
850 double tol,
int tau,
int flag){
852 int i, iterStep, m, tFlag=0, n=nn+1;
853 double alphap=0, alpha=1, beta=0, temp;
854 int* S=(
int *) malloc(
sizeof(
int)*nn);
855 double gapp=-1, gappp=-1;
856 int numS=-1, numSp=-2, numSpp=-3;;
872 for (iterStep=1; iterStep<=maxStep; iterStep++){
877 beta=(alphap -1 ) / alpha;
888 s[i]=z[i]+ beta* z0[i];
891 s[i] =z[i] + beta* z0[i];
892 s[i+1] =z[i+1] + beta* z0[i+1];
893 s[i+2] =z[i+2] + beta* z0[i+2];
894 s[i+3] =z[i+3] + beta* z0[i+3];
895 s[i+4] =z[i+4] + beta* z0[i+4];
914 g[0]=s[0] + s[0] - s[1] - Av[0];
915 for (i=1;i<nn-1;i++){
916 g[i]= - s[i-1] + s[i] + s[i] - s[i+1] - Av[i];
918 g[nn-1]= -s[nn-2] + s[nn-1] + s[nn-1] - Av[nn-1];
929 for (i=m; i <nn; i+=7){
948 for (i=m;i<nn; i+=5){
949 z[i] = s[i] - g[i] /4;
950 z[i+1] = s[i+1] - g[i+1]/4;
951 z[i+2] = s[i+2] - g[i+2]/4;
952 z[i+3] = s[i+3] - g[i+3]/4;
953 z[i+4] = s[i+4] - g[i+4]/4;
991 alpha=(1+sqrt(4*alpha*alpha+1))/2;
996 if (iterStep%tau==0){
1057 numS =
supportSet(x, v, z, g, S, lambda, nn);
1076 if ( abs(numS-numSp) < m){
1079 g, s, S, lambda, nn);
1090 if ( (*gap ==gappp) && (numS==numSpp) ){
1107 s[i]=Av[i] - x[i+1] + x[i];
1141 for (i=m; i<nn; i+=7){
1177 temp=temp + fabs(z0[i])
1214 g[0]=z[0] + z[0] - z[1] - Av[0];
1215 for (i=1;i<nn-1;i++){
1216 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1219 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1225 s[i]=-lambda + z[i];
1234 temp=temp + s[i] *g[i]
1285 g[0]=z[0] + z[0] - z[1] - Av[0];
1286 for (i=1;i<nn-1;i++){
1287 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1290 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1296 s[i]=-lambda + z[i];
1305 temp=temp + s[i] *g[i]
1322 temp+=z[i] * (g[i] - Av[i]);
1325 temp=temp + z[i] * (g[i] - Av[i])
1326 + z[i+1]* (g[i+1]- Av[i+1])
1327 + z[i+2]* (g[i+2]- Av[i+2])
1328 + z[i+3]* (g[i+3]- Av[i+3])
1329 + z[i+4]* (g[i+4]- Av[i+4]);
1362 temp=temp + fabs(z0[i])
1391 if ( (flag !=1) || (tFlag==0) ){
1394 x[i]= v[i] - z[i-1] + z[i];
1395 x[n-1]=v[n-1] - z[n-2];
1398 if ( (flag==1) && (tFlag==1)){
1420 g[0]=z[0] + z[0] - z[1] - Av[0];
1421 for (i=1;i<nn-1;i++){
1422 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1424 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1435 for (i=m;i<nn; i+=5){
1436 z[i] = z[i] - g[i] /4;
1437 z[i+1] = z[i+1] - g[i+1]/4;
1438 z[i+2] = z[i+2] - g[i+2]/4;
1439 z[i+3] = z[i+3] - g[i+3]/4;
1440 z[i+4] = z[i+4] - g[i+4]/4;
1448 for (i=0;i<nn; i++){
1463 g[0]=z[0] + z[0] - z[1] - Av[0];
1464 for (i=1;i<nn-1;i++){
1465 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1467 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1471 g, s, S, lambda, nn);
1499 double *z,
double * v,
double * Av,
1500 double lambda,
int nn,
int maxStep,
1501 double *s,
double *g,
1502 double tol,
int tau){
1508 int* S=(
int *) malloc(
sizeof(
int)*nn);
1510 int numS=-1, numSp=-1;
1518 for (iterStep=1; iterStep<=maxStep; iterStep++){
1521 g[0]=z[0] + z[0] - z[1] - Av[0];
1522 for (i=1;i<nn-1;i++){
1523 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1525 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1528 numS =
supportSet(x, v, z, g, S, lambda, nn);
1540 s[i]=Av[i] - x[i+1] + x[i];
1566 if (iterStep%tau==0){
1581 if ( (*gap <1) && (numS==numSp) && (gapp == *gap) ){
1608 int sfa_one(
double *x,
double *gap,
int * activeS,
1609 double *z,
double * v,
double * Av,
1610 double lambda,
int nn,
int maxStep,
1611 double *s,
double *g,
1612 double tol,
int tau){
1618 int* S=(
int *) malloc(
sizeof(
int)*nn);
1619 double gapp=-1, gappp=-2;
1620 int numS=-100, numSp=-200, numSpp=-300;
1646 g[0]=z[0] + z[0] - z[1] - Av[0];
1647 for (i=1;i<nn-1;i++){
1648 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1650 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1661 for (i=m;i<nn; i+=5){
1662 z[i] = z[i] - g[i] /4;
1663 z[i+1] = z[i+1] - g[i+1]/4;
1664 z[i+2] = z[i+2] - g[i+2]/4;
1665 z[i+3] = z[i+3] - g[i+3]/4;
1666 z[i+4] = z[i+4] - g[i+4]/4;
1674 for (i=0;i<nn; i++){
1690 g[0]=z[0] + z[0] - z[1] - Av[0];
1691 for (i=1;i<nn-1;i++){
1692 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1694 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1696 for (iterStep=1; iterStep<=maxStep; iterStep++){
1706 numS =
supportSet(x, v, z, g, S, lambda, nn);
1718 s[i]=Av[i] - x[i+1] + x[i];
1755 g[0]=z[0] + z[0] - z[1] - Av[0];
1756 for (i=1;i<nn-1;i++){
1757 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1759 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1770 for (i=m;i<nn; i+=5){
1771 z[i] = z[i] - g[i] /4;
1772 z[i+1] = z[i+1] - g[i+1]/4;
1773 z[i+2] = z[i+2] - g[i+2]/4;
1774 z[i+3] = z[i+3] - g[i+3]/4;
1775 z[i+4] = z[i+4] - g[i+4]/4;
1783 for (i=0;i<nn; i++){
1798 g[0]=z[0] + z[0] - z[1] - Av[0];
1799 for (i=1;i<nn-1;i++){
1800 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1802 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1805 if (iterStep % tau==0){
1845 if ( abs( numS-numSp) <m ){
1852 g, s, S, lambda, nn);
1863 if ( (*gap ==gappp) && (numS==numSpp) ){
int sfa_one(double *x, double *gap, int *activeS, double *z, double *v, double *Av, double lambda, int nn, int maxStep, double *s, double *g, double tol, int tau)
void Rose(double *zMax, double *z0, double *Av, int nn)
int generateSolution(double *x, double *z, double *gap, double *v, double *Av, double *g, double *s, int *S, double lambda, int nn)
void restartMapping(double *g, double *z, double *v, double lambda, int nn)
void Thomas(double *zMax, double *z0, double *Av, int nn)
int sfa_special(double *x, double *gap, int *activeS, double *z, double *v, double *Av, double lambda, int nn, int maxStep, double *s, double *g, double tol, int tau)
void dualityGap2(double *gap, double *z, double *g, double *s, double *Av, double lambda, int nn)
void dualityGap(double *gap, double *z, double *g, double *s, double *Av, double lambda, int nn)
int sfa(double *x, double *gap, int *activeS, double *z, double *z0, double *v, double *Av, double lambda, int nn, int maxStep, double *s, double *g, double tol, int tau, int flag)
int supportSet(double *x, double *v, double *z, double *g, int *S, double lambda, int nn)