Skip to content

Commit 7ef2b83

Browse files
author
Leily Rabbani
committed
loops are limited to nochange < 100
1 parent 3c38a7c commit 7ef2b83

File tree

1 file changed

+41
-36
lines changed

1 file changed

+41
-36
lines changed

‎src/KRBalancing.cpp

Lines changed: 41 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -36,36 +36,56 @@ class kr_balancing{
3636
I.setIdentity();
3737
I = I*0.00001;
3838
A = A + I;
39-
x0 = e.sparseView();
39+
x = e.sparseView();
40+
// Eigen::MatrixXd test;
41+
// test.resize(A.rows(),A.cols());
42+
// SparseMatrixCol temp = test.sparseView();
43+
// A = A + temp;
44+
// A = A*1.00001;
45+
// for (int k=0; k<A.outerSize(); ++k){
46+
// for(SparseMatrixCol::InnerIterator it(A,k); it; ++it){
47+
// std::cout << it.value() << " , ";
48+
// }
49+
// }
4050
}
4151

4252
void outter_loop(){
4353
size_t outter_loop_count = 0;
4454
double stop_tol = tol*0.5;
4555
double eta = etamax;
46-
x = x0;
56+
// x = x0;
4757
assert(x.isVector());
4858
double rt = std::pow(tol,2);
4959
v = x.cwiseProduct(A*x);
5060
rk = Eigen::VectorXd::Constant(v.rows(),1) - v; //Vecotr of 1 - v
5161
assert(rk.isVector());
5262
rho_km1 = rk.conjugate().transpose()*rk;
63+
assert(rho_km1.coeff(0,0)>=0);
5364
rho_km2 = rho_km1;
5465
double rout = rho_km1.coeff(0,0);
5566
double rold = rout;
5667
if(fl == 1) std::cout<< 'intermediate convergence statistics is off'<<std::endl;
57-
while(rout > rt){//outer itteration
68+
size_t nochange = 0;
69+
while(rout > rt && nochange < 100){//outer itteration
5870
outter_loop_count ++;
5971
i=i+1; k=0; y=e.sparseView();
6072
innertol = std::max(std::pow(eta,2)*rout,rt);
6173
inner_loop();
62-
74+
assert(rho_km1.coeff(0,0) > 0);
75+
assert(rho_km1.coeff(0,0) != std::numeric_limits<double>::infinity());
6376
x=x.cwiseProduct(y);
77+
// std::cout << "x" << x << std::endl;
6478
assert(x.isVector());
65-
66-
v=x.cwiseProduct(A*x);
79+
v=x.cwiseProduct((A*x));
80+
// std::cout << "v "<< v << std::endl;
6781
rk = Eigen::VectorXd::Constant(v.rows(),1) - v;
6882
rho_km1 = rk.transpose()*rk;
83+
assert(rho_km1.coeff(0,0) > 0);
84+
assert(rho_km1.coeff(0,0) != std::numeric_limits<double>::infinity());
85+
if(std::abs(rout - rho_km1.coeff(0,0)) < 0.0000001 ||rho_km1.coeff(0,0) ==std::numeric_limits<double>::infinity()){
86+
nochange++;
87+
// std::cout << "nochange " << rho_km1.coeff(0,0)<<std::endl;
88+
}
6989
rout = rho_km1.coeff(0,0);
7090
MVP = MVP + k + 1;
7191
//Update inner iteration stopping criterion.
@@ -88,33 +108,39 @@ class kr_balancing{
88108
k++;
89109
if(k == 1){
90110
Z = rk.cwiseQuotient(v);
111+
// std::cout << "Z "<< Z << std::endl;
112+
// std::cout << "rk " << rk << std::endl;
113+
// std::cout << "v "<< v <<std::endl;
91114
p=Z;
92115
rho_km1 = rk.transpose()*Z;
93116
}else{
117+
// std::cout << "rho_km2 " << rho_km2 << std::endl;
94118
Eigen::VectorXd beta=rho_km1.cwiseQuotient(rho_km2);
95119
p = Z + (beta(0)*p);
96120

97121
}
122+
// std::cout << "p "<< p <<std::endl;
98123
//Update search direction efficiently.
99124
SparseMatrixCol w=x.cwiseProduct(A*(x.cwiseProduct(p))) + v.cwiseProduct(p);
100125
SparseMatrixCol alpha_mat = rho_km1.cwiseQuotient(p.conjugate().transpose()*w);
101126
double alpha = alpha_mat.coeff(0,0);
102127
Eigen::VectorXd ap = alpha * p;
103128
//Test distance to boundary of cone.
104129
Eigen::VectorXd ynew = y + ap;
105-
130+
// std::cout << "ynew.minCoeff() " << ynew.minCoeff() <<std::endl;
106131
if(ynew.minCoeff() <= delta){
107132
if(delta == 0) break;
108133
Eigen::Matrix<bool, Eigen::Dynamic , Eigen::Dynamic> ind_helper = (ap.array()<0);
109134
Eigen::VectorXi ind = Eigen::VectorXi::LinSpaced(ind_helper.size(),0,ind_helper.size()-1);
110135
ind.conservativeResize(std::stable_partition(
111136
ind.data(), ind.data()+ind.size(), [&ind_helper](int i){return ind_helper(i);})-ind.data());
112137
Eigen::VectorXd y_copy = y;
113-
for(i = 0; i < ind.size(); i++){ //TODO proper masking? //ind_vec.unaryExpr(x);
138+
for(size_t i = 0; i < ind.size(); i++){ //TODO proper masking? //ind_vec.unaryExpr(x);
114139
y_copy(ind(i)) = (delta - y_copy(ind(i)))/ap(i);
115140
}
116141
double gamma = y_copy.minCoeff();
117142
y = y + gamma*ap;
143+
// std::cout << y << std::endl;
118144
break;
119145
}
120146
if(ynew.minCoeff() >= Delta){
@@ -123,7 +149,7 @@ class kr_balancing{
123149
ind.conservativeResize(std::stable_partition(
124150
ind.data(), ind.data()+ind.size(), [&ind_helper](int i){return ind_helper(i);})-ind.data());
125151
Eigen::VectorXd y_copy = y;
126-
for(i = 0; i < ind.size(); i++){ //TODO proper masking? //ind_vec.unaryExpr(x);
152+
for(size_t i = 0; i < ind.size(); i++){ //TODO proper masking? //ind_vec.unaryExpr(x);
127153
y_copy(ind(i)) = (Delta - y_copy(ind(i)))/ap(i);
128154
}
129155
double gamma = y_copy.minCoeff();
@@ -134,46 +160,25 @@ class kr_balancing{
134160
y = ynew;
135161
rk = rk - (alpha*w); rho_km2 = rho_km1;
136162
Z = rk.cwiseQuotient(v); rho_km1 = rk.conjugate().transpose()*Z;
163+
// std::cout << "rho_km1 "<< rho_km1 <<std::endl;
137164
if(inner_loop_count %100 ==0) std::cout << "inner loop number "<< inner_loop_count <<std::endl;
138165

139166

140167
}//End of the inner 'while'
141168
}
142169
const SparseMatrixCol* get_output(){
143-
// std::cout << "export output of shape "<< A.cwiseProduct(x*x.transpose()).rows() <<
144-
// " by " << A.cwiseProduct(x*x.transpose()).rows()<<std::endl;
145-
// std::cout << (x*x.transpose()).triangularView<Eigen::Lower>()<<std::endl;
146-
// SparseMatrixCol D = x*x.transpose();
147-
// Eigen::SparseTriangularView<Eigen::SparseMatrix<double>, Eigen::Lower> MView = D.triangularView<Eigen::Lower>();
148-
// std::cout << D.triangularView<Eigen::Lower>() <<std::endl;
149-
// std::cout << "export output"<<std::endl;
150-
// output = (A.triangularView<Eigen::Lower>()).cwiseProduct((x*x.transpose()).triangularView<Eigen::Lower>());
151-
// std::cout << A.cwiseProduct(x*x.transpose()) << std::endl;
152-
// SparseMatrixCol diag_x;
153-
// diag_x.resize(x.rows(),x.rows());
154-
// diag_x = Eigen::MatrixXd(Eigen::MatrixXd(x).asDiagonal()).sparseView();
155-
// diag_x.makeCompressed();
156-
// std::cout << "diag x size " << diag_x.rows() << " x " << diag_x.cols()<<std::endl;
157-
158-
// Eigen::MatrixXd dense_mat_x = dense_x.replicate(1,x.rows());
159-
// Eigen::MatrixXd dense_mat_x = dense_x*dense_x.transpose();
160-
// std::cout << diag_x <<std::endl;
161-
// SparseMatrixCol test = dense_x.sparseView();
162-
163-
// std::cout << test <<std::endl;
164-
// SparseMatrixCol temp = (dense_x.asDiagonal()*A).sparseView();
165-
// SparseMatrixCol temp1 = temp.cwiseProduct(x);
166-
// std::cout << "temp" <<std::endl;
167170
assert(A.rows() == A.cols());
168171
A = SparseMatrixCol(A.triangularView<Eigen::Lower>());
172+
// std::cout << "x is "<< std::endl;
173+
// std::cout << x << std::endl;
169174
#pragma omp parallel for num_threads(num_threads) schedule(dynamic)
170175
for (int k=0; k<A.outerSize(); ++k){
171176
for(SparseMatrixCol::InnerIterator it(A,k); it; ++it){
172177
it.valueRef() = it.value()*x.coeff(it.row(),0)*x.coeff(it.col(),0);
173-
178+
// std::cout << it.value() << " ";
174179
}
175180
}
176-
181+
// std::cout << " "<<std::endl;
177182
return &A;
178183
}
179184
private:
@@ -184,7 +189,7 @@ class kr_balancing{
184189
double tol = 1e-6;
185190
double g = 0.9;
186191
double etamax = 0.1;
187-
SparseMatrixCol x0;
192+
// SparseMatrixCol x0;
188193
Eigen::MatrixXd e;
189194
SparseMatrixCol A;
190195
SparseMatrixCol rho_km1;

0 commit comments

Comments
 (0)