0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024 function results = TEBD_NC_v2_2(par,prob)
0025
0026
0027 sigma = par.sigma;
0028
0029
0030
0031 x = prob.x0;
0032 u = prob.u0;
0033
0034
0035
0036 par.AlgName = '2EBD-v2.2-NoCell';
0037 fprintf('\n\n ****Begin 2EBD V2.2. iterations = %d ****\n', par.maxit);
0038 fprintf(par.AlgName);
0039 fprintf('\n');
0040 fprintf(TEBD_LegendAppend(par));
0041 fprintf('\n');
0042
0043 tstart = tic;
0044
0045
0046
0047
0048
0049 if isfield(par,'balancedinit') && par.balancedinit >0
0050 theta = TEBD_GetInitialScaling_NC_v2_2(prob,x,u,sigma);
0051 else
0052 theta = 1;
0053 end
0054
0055 lambda = sigma/sqrt(theta);
0056
0057
0058
0059
0060
0061 dynamicscaling = par.dynamicscaling;
0062
0063
0064 isdynamicscaling = 0;
0065
0066 if dynamicscaling == 1
0067
0068
0069 dyn_scale_update_it = par.dyn_scale_updateiteration;
0070
0071
0072
0073 scaleratio = par.scaleratio;
0074
0075
0076
0077
0078 scalecorrection = par.scalecorrection;
0079 end
0080
0081
0082 TEBD_PrintHeaders;
0083
0084
0085 runhist.preprocesstime = toc(tstart);
0086
0087
0088 stopstatus = 1;
0089
0090 tstart = tic;
0091
0092 Z_tilde.pos = 1;
0093
0094
0095 for iter = 1:par.maxit;
0096
0097
0098
0099 gradf_x = prob.gradf(x);
0100 gradf_x_plus_u = gradf_x+u;
0101 [x_tilde,Z_tilde] = prob.resolvent1(x - gradf_x_plus_u*(theta*lambda),Z_tilde,lambda*theta);
0102
0103
0104 u_plus_lambda_x = u +x_tilde*lambda;
0105 [resolv2] = prob.resolvent2(u_plus_lambda_x/lambda,lambda);
0106 u_tilde = u_plus_lambda_x - resolv2*lambda;
0107
0108
0109 v_1 = ((x-x_tilde)/lambda) + ((u_tilde-u)*theta);
0110 v_2 = (u-u_tilde)/lambda;
0111
0112 if par.adaptive == 0
0113 lambda_a = lambda;
0114 else
0115
0116 [lambda_a] = GetAdaptiveLambdaScaled(prob,lambda,sigma,x,x_tilde,v_1,u,u_tilde,v_2,gradf_x,theta);
0117 end
0118
0119
0120 x = x - v_1*lambda_a;
0121 u = u - v_2*lambda_a;
0122
0123
0124 secs = toc(tstart);
0125
0126
0127 dh2 = v_2+x_tilde;
0128 dh1 = v_1/theta - gradf_x - u_tilde;
0129 obj = prob.evalf2(x_tilde,u_tilde,dh2,dh1);
0130 residuals = prob.residuals2(x_tilde,u_tilde,v_1/theta,v_2,obj);
0131
0132
0133 runhist = TEBD_SaveHistory(obj,runhist,iter,residuals,secs,lambda_a/lambda,theta);
0134
0135
0136 [stopstatus,printstatus,statusstr,runhist] = TEBD_Check_StopCrit(par,residuals,runhist,iter,stopstatus);
0137 if printstatus
0138 TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0139 fprintf(statusstr);
0140 end
0141
0142 if residuals.resnorm < par.tol || iter >= par.maxit
0143
0144
0145 TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0146 fprintf(statusstr);
0147
0148
0149 [~,y_tilde] = prob.resolvent2(u_plus_lambda_x/lambda,lambda);
0150 obj_d = prob.evalf(x_tilde,y_tilde,Z_tilde);
0151 residuals_d = prob.residuals(x_tilde,u_tilde,v_1/theta,v_2,obj_d,y_tilde,Z_tilde,u_plus_lambda_x,lambda);
0152 results = TEBD_GetFinalValues_2(runhist,iter,par,prob,x,u,u_tilde,x_tilde,y_tilde,Z_tilde,residuals_d,secs);
0153 break
0154 end
0155
0156
0157 if iter == 1 || mod(iter,par.printit)==0
0158
0159
0160 TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0161 end
0162
0163
0164 if dynamicscaling == 1 && (max(residuals.res1,residuals.res2) > par.tol2||(residuals.res1>residuals.res2*5)||(residuals.res2>residuals.res1*5))
0165 if isdynamicscaling == 0
0166 fprintf(' ****Starting dynamic scaling at iteration = %d ****\n', iter);
0167 isdynamicscaling = 1;
0168 end
0169 if mod(iter,dyn_scale_update_it)==0
0170 wrel_norm2 = wrel_norm2 * residuals.res2^(1/dyn_scale_update_it);
0171 wrel_norm1 = wrel_norm1 * residuals.res1^(1/dyn_scale_update_it);
0172 if wrel_norm2/wrel_norm1 > scaleratio
0173 theta = theta/(scalecorrection^2);
0174 TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0175 elseif wrel_norm1/wrel_norm2 > scaleratio
0176 theta = theta*(scalecorrection^2);
0177 TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0178 end
0179 lambda = sigma/sqrt(theta);
0180 elseif mod(iter,dyn_scale_update_it) == 1
0181
0182 wrel_norm2 = residuals.res2^(1/dyn_scale_update_it);
0183 wrel_norm1 = residuals.res1^(1/dyn_scale_update_it);
0184 else
0185
0186 wrel_norm2 = wrel_norm2 * residuals.res2^(1/dyn_scale_update_it);
0187 wrel_norm1 = wrel_norm1 * residuals.res1^(1/dyn_scale_update_it);
0188 end
0189 elseif isdynamicscaling == 1 && dynamicscaling == 1
0190 fprintf(' ****Stopping dynamic scaling at iteration = %d ****\n', iter);
0191 isdynamicscaling = 0;
0192 end
0193
0194 end
0195
0196 fprintf('\n\n ****End 2EBD-HPE algorithm ****\n');
0197
0198
0199
0200
0201
0202 end
0203
0204 function [lambda_a] = GetAdaptiveLambdaScaled(prob,lambda,sigma,x,x_tilde,v_1,u,u_tilde,v_2,gradf_x,theta)
0205
0206 eps = prob.error1(x,x_tilde,gradf_x);
0207
0208 lambda_sq = lambda^2;
0209
0210
0211 v_x = v_1;
0212 v_u = v_2;
0213 Norm_v_x_sq = norm(v_x)^2;
0214 Norm_v_u_sq = norm(v_u)^2;
0215
0216
0217 Delta_x = (x_tilde-x);
0218 Delta_u = (u_tilde-u);
0219 Norm_Delta_x = norm(Delta_x);
0220 Norm_Delta_u = norm(Delta_u);
0221 Norm_Delta_x_sq = Norm_Delta_x^2;
0222 Norm_Delta_u_sq = Norm_Delta_u^2;
0223
0224 tr_v_Delta_x = v_x'*Delta_x;
0225 tr_v_Delta_u = v_u'*Delta_u;
0226
0227 alpha = sigma^2*(Norm_Delta_x_sq/theta+Norm_Delta_u_sq);
0228
0229 a = (Norm_v_x_sq)/theta + Norm_v_u_sq;
0230 b = (lambda*Norm_v_x_sq+tr_v_Delta_x)/theta + lambda*Norm_v_u_sq + tr_v_Delta_u + eps;
0231 c = (lambda_sq*Norm_v_x_sq+2*lambda*tr_v_Delta_x+Norm_Delta_x_sq)/theta+lambda_sq*Norm_v_u_sq+Norm_Delta_u_sq+2*lambda*(tr_v_Delta_u+eps)-alpha;
0232 discriminant = sqrt(b^2-a*c);
0233 if b <= 0
0234 beta = (-b+discriminant)/a;
0235 else
0236 beta = -c/(b+discriminant);
0237 end
0238 if beta <= 1e-13
0239 beta = 0;
0240 else
0241
0242 end
0243 lambda_a = lambda + beta;
0244 end