//*************************************************************
//	resolvendo assignments com o método húngaro
//	gilberto miranda jr. - 2006
//*************************************************************
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <iomanip>
#include <fstream>
#include "lap.h"
#include "lap.cpp"
#include "splx_lap.cpp"

using namespace std;

typedef int boolean;

int main(int argc, char *argv[]) {

  int h, i, j, k, l, m, n, n2, n3, n4, matrix_size, le, *d, iter, **x, **xb;
  double **a, *f, *u, *v, *t, *w, lb, lbb, ub, ubb, cf, cv, gap, **c, *y, *r;
  double splx, splt, stpy, stpw, stpt, lambda, ynorma, ynorma2, wnorma, wnorma2, tnorma, tnorma2, wh;
  float lef;

  ifstream data(argv[1]);

  data>>n;;

  n2 = n*n;
  n3 = n2*n;
  n4 = n3*n;
  
//  alocacao dinamica : 
  a = new double* [n];
  for(i=0; i<n; i++) a[i] = new double [n]; 
  v = new double [n3];
  u = new double [n3]; 
  r = new double [n3]; 
  d = new int [n4];
  f = new double [n4];
  x = new int* [n];
  for(i=0; i<n; i++) x[i] = new int [n]; 
  x = new int* [n];
  for(i=0; i<n; i++) x[i] = new int [n]; 
  c = new double* [n];
  for(i=0; i<n; i++) c[i] = new double [n];

//  leitura de dados :
  for(i=0; i<n; i++){
   for(j=0; j<n; j++){
     data>>a[i][j];
   }
  }

for(k=0; k<n; k++){
 for(i=0; i<n; i++){
  for(l=0; l<n; l++){
   for(j=0; j<n; j++){
    data>>d[k*n3+i*n2+l*n+j];
   }
  }
 }
}

//  inicializando upperbound :
  ub = ubb = 1.e6;

  // solucao inicial:
  lb = splx_lap(n, a, x);
  lbb = lb;
  printf(" Limite Inferior Inicial: %10.5f \n\n", lbb);
  ubb = 0;
  for(k=0; k<n; k++){
   for(i=0; i<n; i++){
    ubb+= a[k][i]*x[k][i];
    for(l=(k+1); l<n; l++){
     for(j=0; j<n; j++){
      ubb += (d[k*n3+i*n2+l*n+j]+d[l*n3+j*n2+k*n+i])*x[l][j]*x[k][i];
     }
    }
   }
  }

  if (ub<ubb) ubb = ub;
  printf(" Limite Superior Inicial: %10.5f \n\n", ubb);
  gap = 100*(ubb-lbb)/ubb;
  printf(" Gap-otimalidade Inicial: %10.5f \n\n", gap);

  // dimensionando subgradientes:

   w = new double [n3]; // u
   y = new double [n3]; // v
   t = new double [n3]; // r

   // inicializando vetor de precos:
    for (i=0; i<n3; i++) w[i] = 0.0;
    for (i=0; i<n3; i++) y[i] = 0.0;
    for (i=0; i<n3; i++) t[i] = 0.0;
    for (i=0; i<n3; i++) v[i] = 0.0;
    for (i=0; i<n3; i++) u[i] = 0.0;
    for (i=0; i<n3; i++) r[i] = 0.0;
    for (i=0; i<n4; i++) f[i] = 0.0;
    lambda = atof(argv[2]);

    // disparando método do subgradiente :
    iter = 0;
    int iternmelhora = 0;
    while(1){
    //while(iter<=8850){

     //inicializando matriz c :
      for(k=0; k<n; k++){
       for(i=0; i<n; i++){
        c[k][i] = a[k][i];
        for(j=0; j<n; j++) if(i!=j)c[k][i] -= r[k*n2+i*n+j];
       }
      }

     // subproblema lagrangeano em x :
     // montando matriz c:
	for(k=0;k<n;k++){
	    for(l=(k+1);l<n;l++){
                for(j=0;j<n;j++){
	         c[l][j] -= v[k*n2+l*n+j];
	        }
	        for(i=0;i<n;i++){
	         c[k][i] -= u[k*n2+l*n+i];
	        }
            }
	}

     // tratando custos negativos:
     double Min = 1.e8;
     int flag1 = 0;
      for(k=0; k<n; k++){
       for(i=0; i<n; i++){
        if(c[k][i]<Min) Min = c[k][i];
       }
      }
      if(Min < 0.0){
	      for(k=0; k<n; k++){
	       for(i=0; i<n; i++){
	        c[k][i]-=Min;
	       }
	      }
	flag1 = 1;
      }
 
     splx = splx_lap(n, c, x);
     if(flag1) splx+=n*Min;

       // subproblema lagrangeano em f :
         splt = 0.0;
           for(k=0; k<n; k++){
             for(l=(k+1); l<n; l++){
                 for(i=0; i<n; i++){
                    for(j=0; j<n; j++){
			if(i!=j){
			double del = d[k*n3+i*n2+l*n+j]+d[l*n3+j*n2+k*n+i]+v[k*n2+l*n+j]+u[k*n2+l*n+i]+r[k*n2+i*n+j]+r[l*n2+j*n+i];
			if(del>1.e-12) {splt+=0.0; f[k*n3+i*n2+l*n+j]=0;}
			else {splt+= del; f[k*n3+i*n2+l*n+j]=1;}
			}
                    }
                 }
             }
           }
       // calculando lowerbound lagrangeano
       lb = splx + splt; 

       int imprimir = 0; 
       if (lb>lbb) 
	{ 
	   lbb = lb; 
	   imprimir = 1;
	   iternmelhora = 0;
           gap = 100*(ubb-lbb)/ubb;
	}
	else iternmelhora++;

       // determinando subgradiente y:
        for(k=0; k<n; k++){
          for(l=(k+1); l<n; l++){
	   for(j=0; j<n; j++){
            y[k*n2+l*n+j] = 0.0;
            for(i=0; i<n; i++) y[k*n2+l*n+j] += f[k*n3+i*n2+l*n+j];
            y[k*n2+l*n+j] -= x[l][j];
           }
          }
        }

        // determinando norma-2 do subgradiente :
        ynorma2 = 0.0;
        for(k=0; k<n; k++){
          for(l=(k+1); l<n; l++){
	   for(j=0; j<n; j++){
            ynorma2 += y[k*n2+l*n+j] * y[k*n2+l*n+j];
          }
         }
        }
        ynorma = sqrt(ynorma2);

       // determinando subgradiente w:
        for(k=0; k<n; k++){
          for(l=(k+1); l<n; l++){
	   for(i=0; i<n; i++){
            w[k*n2+l*n+i] = 0.0;
            for(j=0; j<n; j++) w[k*n2+l*n+i] += f[k*n3+i*n2+l*n+j];
            w[k*n2+l*n+i] -= x[k][i];
           }
          }
        }

        // determinando norma-2 do subgradiente :
        wnorma2 = 0.0;
        for(k=0; k<n; k++){
          for(l=(k+1); l<n; l++){
	   for(i=0; i<n; i++){
            wnorma2 += w[k*n2+l*n+i] * w[k*n2+l*n+i];
          }
         }
        }
        wnorma = sqrt(wnorma2);

       // determinando subgradiente t:
        for(k=0; k<n; k++){
          for(i=0; i<n; i++){
	   for(j=0; j<n; j++){
            if(i!=j){
            t[k*n2+i*n+j] = 0.0;
            for(l=(k+1); l<n; l++) t[k*n2+i*n+j] += f[k*n3+i*n2+l*n+j];
                for(l=0; l<k; l++) t[k*n2+i*n+j] += f[l*n3+j*n2+k*n+i];
            t[k*n2+i*n+j] -= x[k][i];
            }
           }
          }
        }

        // determinando norma-2 do subgradiente :
        tnorma2 = 0.0;
        for(k=0; k<n; k++){
          for(i=0; i<n; i++){
	   for(j=0; j<n; j++){
            tnorma2 += t[k*n2+i*n+j] * t[k*n2+i*n+j];
          }
         }
        }
        tnorma = sqrt(wnorma2);

        // determinando passo do subgradiente :
        stpy = lambda*(ubb - lb)/ynorma2;
        stpw = lambda*(ubb - lb)/wnorma2;
        stpt = lambda*(ubb - lb)/tnorma2;
        //atualizando variáveis duais :
        for(k=0; k<n; k++){
           for(l=(k+1); l<n; l++){
             for(j=0; j<n; j++){
             v[k*n2+l*n+j] = v[k*n2+l*n+j] + stpy*y[k*n2+l*n+j];
             }
           }
        }
        for(k=0; k<n; k++){
           for(l=(k+1); l<n; l++){
             for(i=0; i<n; i++){
             u[k*n2+l*n+i] = u[k*n2+l*n+i] + stpw*w[k*n2+l*n+i];
             }
           }
        }
        for(k=0; k<n; k++){
           for(i=0; i<n; i++){
             for(j=0; j<n; j++){
             r[k*n2+i*n+j] = r[k*n2+i*n+j] + stpt*t[k*n2+i*n+j];
             }
           }
        }

        iter ++;
      ub = 0;
      for(k=0; k<n; k++){
       for(i=0; i<n; i++){
        ub+= a[k][i]*x[k][i];
        for(l=(k+1); l<n; l++){
         for(j=0; j<n; j++){
          ub += (d[k*n3+i*n2+l*n+j]+d[l*n3+j*n2+k*n+i])*x[l][j]*x[k][i];
         }
        }
       }
      }

      if (ub<ubb){ ubb = ub;
      gap = 100*(ubb-lbb)/ubb;
      }

       // impressao intermediária :
       if(imprimir) printf("  Iteracao: %d --- Lowerbound: %10.5f --- Upperbound: %10.5f --- Gap: %10.5f \n", iter, lbb, ubb, gap);
       //printf("  Iteracao: %d --- Lowerbound: %10.5f --- NormaSG: %10.5f \n", iter, lbb, ynorma);
        int iterlim = 500 + iter/100;
      if (iternmelhora==iterlim) {
           lambda = 0.5*lambda;
           cout<<"diminuindo passo : "<<iter<<" "<<lambda<<"\n";
           iternmelhora = 0;
	   printf("  Iteracao: %d --- Lowerbound: %10.5f --- Upperbound: %10.5f --- Gap: %10.5f \n", iter, lbb, ubb, gap);
       }
      if(lambda < 1.e-8) break;

    } // end while : fim do método do subgradiente

 

  return 0;
}

