Fixed convergence issue

Author: Leily Rabbani

src/KRBalancing.cpp

@@ -34,26 +34,15 @@ class kr_balancing{
         SparseMatrixCol I;
         I.resize(A.rows(),A.cols());
         I.setIdentity();
-        I = I*0.00001;
-        A = A + I;
-        x = e.sparseView();
-        // Eigen::MatrixXd test;
-        // test.resize(A.rows(),A.cols());
-        // SparseMatrixCol temp = test.sparseView();
-        // A = A + temp;
-        // A = A*1.00001;
-        // for (int k=0; k<A.outerSize(); ++k){
-        //     for(SparseMatrixCol::InnerIterator it(A,k); it; ++it){
-        //       std::cout << it.value() << " , ";
-        //     }
-        // }
+  //      I = I*0.00001;
+  //      A = A + I;
       }
 
       void outter_loop(){
+        x = e.sparseView();
         size_t outter_loop_count = 0;
         double stop_tol = tol*0.5;
         double eta = etamax;
-  //      x = x0;
         assert(x.isVector());
         double rt = std::pow(tol,2);
         v = x.cwiseProduct(A*x);
@@ -66,26 +55,22 @@ class kr_balancing{
         double rold = rout;
         if(fl == 1) std::cout<< 'intermediate convergence statistics is off'<<std::endl;
         size_t nochange = 0;
-        while(rout > rt && nochange < 100){//outer itteration
+        while(rout > rt){//outer itteration
           outter_loop_count ++;
+//          if(outter_loop_count == 30) break;
           i=i+1; k=0; y=e.sparseView();
           innertol = std::max(std::pow(eta,2)*rout,rt);
           inner_loop();
-          assert(rho_km1.coeff(0,0) > 0);
+  //        assert(rho_km1.coeff(0,0) > 0);
           assert(rho_km1.coeff(0,0) != std::numeric_limits<double>::infinity());
           x=x.cwiseProduct(y);
-//          std::cout << "x" << x << std::endl;
           assert(x.isVector());
           v=x.cwiseProduct((A*x));
-  //        std::cout << "v "<< v << std::endl;
           rk = Eigen::VectorXd::Constant(v.rows(),1) - v;
-          rho_km1 = rk.transpose()*rk;
-          assert(rho_km1.coeff(0,0) > 0);
+          rho_km1 = rk.conjugate().transpose()*rk;
+  //        assert(rho_km1.coeff(0,0) > 0);
           assert(rho_km1.coeff(0,0) != std::numeric_limits<double>::infinity());
-          if(std::abs(rout - rho_km1.coeff(0,0)) < 0.0000001 ||rho_km1.coeff(0,0) ==std::numeric_limits<double>::infinity()){
-            nochange++;
-  //          std::cout << "nochange " << rho_km1.coeff(0,0)<<std::endl;
-          }
+
           rout = rho_km1.coeff(0,0);
           MVP = MVP + k + 1;
           //Update inner iteration stopping criterion.
@@ -99,7 +84,13 @@ class kr_balancing{
             res.push_back(res_norm);
           }
           if(outter_loop_count %50 ==0) std::cout << "outer loop number "<< outter_loop_count <<std::endl;
-         }
+        }//End of outer loop
+    //    if (nochange >= 100) {
+    //       std::cout << "nochange >=100" <<std::endl;
+    //       return 0;
+    //    }else{
+    //    return &x;
+    //    }
       }
       void inner_loop(){
         size_t inner_loop_count = 0;
@@ -108,77 +99,65 @@ class kr_balancing{
           k++;
           if(k == 1){
             Z = rk.cwiseQuotient(v);
-  //          std::cout << "Z "<< Z << std::endl;
-  //          std::cout << "rk " << rk << std::endl;
-  //          std::cout << "v "<< v <<std::endl;
             p=Z;
-            rho_km1 = rk.transpose()*Z;
+            rho_km1 = rk.conjugate().transpose()*Z;
           }else{
-  //          std::cout << "rho_km2 " << rho_km2 << std::endl;
             Eigen::VectorXd beta=rho_km1.cwiseQuotient(rho_km2);
             p = Z + (beta(0)*p);
 
           }
-//          std::cout << "p "<< p <<std::endl;
           //Update search direction efficiently.
           SparseMatrixCol w=x.cwiseProduct(A*(x.cwiseProduct(p))) + v.cwiseProduct(p);
           SparseMatrixCol alpha_mat = rho_km1.cwiseQuotient(p.conjugate().transpose()*w);
           double alpha = alpha_mat.coeff(0,0);
           Eigen::VectorXd ap = alpha * p;
           //Test distance to boundary of cone.
           Eigen::VectorXd ynew = y + ap;
-  //        std::cout << "ynew.minCoeff() " << ynew.minCoeff() <<std::endl;
           if(ynew.minCoeff() <= delta){
             if(delta == 0) break;
             Eigen::Matrix<bool, Eigen::Dynamic , Eigen::Dynamic> ind_helper = (ap.array()<0);
             Eigen::VectorXi ind = Eigen::VectorXi::LinSpaced(ind_helper.size(),0,ind_helper.size()-1);
             ind.conservativeResize(std::stable_partition(
               ind.data(), ind.data()+ind.size(), [&ind_helper](int i){return ind_helper(i);})-ind.data());
-            Eigen::VectorXd y_copy = y;
-            for(size_t i = 0; i < ind.size(); i++){ //TODO proper masking? //ind_vec.unaryExpr(x);
-               y_copy(ind(i)) = (delta - y_copy(ind(i)))/ap(i);
+            Eigen::VectorXd y_copy;
+            y_copy.resize(ind.size());
+            for(size_t i = 0; i < ind.size(); i++){
+               y_copy(i) = (delta - y(ind(i)))/ap(ind(i));
             }
             double gamma = y_copy.minCoeff();
             y = y + gamma*ap;
-  //          std::cout << y << std::endl;
             break;
           }
-          if(ynew.minCoeff() >= Delta){
+          if(ynew.maxCoeff() >= Delta){
             Eigen::Matrix<bool, Eigen::Dynamic , Eigen::Dynamic> ind_helper = (ynew.array() > Delta);
             Eigen::VectorXi ind = Eigen::VectorXi::LinSpaced(ind_helper.size(),0,ind_helper.size()-1);
             ind.conservativeResize(std::stable_partition(
               ind.data(), ind.data()+ind.size(), [&ind_helper](int i){return ind_helper(i);})-ind.data());
-            Eigen::VectorXd y_copy = y;
-            for(size_t i = 0; i < ind.size(); i++){ //TODO proper masking? //ind_vec.unaryExpr(x);
-                 y_copy(ind(i)) = (Delta - y_copy(ind(i)))/ap(i);
+            Eigen::VectorXd y_copy;
+            y_copy.resize(ind.size());
+            for(size_t i = 0; i < ind.size(); i++){
+                 y_copy(i) = (Delta - y(ind(i)))/ap(ind(i));
             }
             double gamma = y_copy.minCoeff();
-            y = y + gamma*ap;
+            y = y + (gamma*ap);
             break;
           }
-
           y = ynew;
           rk = rk - (alpha*w); rho_km2 = rho_km1;
           Z = rk.cwiseQuotient(v); rho_km1 = rk.conjugate().transpose()*Z;
-  //        std::cout << "rho_km1 "<< rho_km1 <<std::endl;
-          if(inner_loop_count %100 ==0) std::cout << "inner loop number "<< inner_loop_count <<std::endl;
-
 
         }//End of the inner 'while'
       }
       const SparseMatrixCol* get_output(){
         assert(A.rows() == A.cols());
         A = SparseMatrixCol(A.triangularView<Eigen::Lower>());
-  //      std::cout << "x is "<< std::endl;
-  //      std::cout << x << std::endl;
+
         #pragma omp parallel for num_threads(num_threads) schedule(dynamic)
         for (int k=0; k<A.outerSize(); ++k){
             for(SparseMatrixCol::InnerIterator it(A,k); it; ++it){
               it.valueRef() = it.value()*x.coeff(it.row(),0)*x.coeff(it.col(),0);
-  //            std::cout << it.value() << " ";
             }
           }
-  //        std::cout << " "<<std::endl;
         return &A;
       }
      private: