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_1(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.1-NoCell';
0037 fprintf('\n\n ****Begin 2EBD V2.1. 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_1(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,y_tilde] = 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 obj = prob.evalf(x_tilde,y_tilde,Z_tilde);
0128 residuals = prob.residuals(x_tilde,u_tilde,v_1/theta,v_2,obj,y_tilde,Z_tilde,u_plus_lambda_x,lambda);
0129
0130
0131 runhist = TEBD_SaveHistory(obj,runhist,iter,residuals,secs,lambda_a/lambda,theta);
0132
0133
0134 [stopstatus,printstatus,statusstr,runhist] = TEBD_Check_StopCrit(par,residuals,runhist,iter,stopstatus);
0135 if printstatus
0136 TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0137 fprintf(statusstr);
0138 end
0139
0140 if residuals.resnorm < par.tol || iter >= par.maxit
0141
0142
0143 TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0144 fprintf(statusstr);
0145
0146
0147 results = TEBD_GetFinalValues_1(runhist,iter,par,prob,x,u,u_tilde,x_tilde,y_tilde,Z_tilde,secs);
0148 break
0149 end
0150
0151
0152 if iter == 1 || mod(iter,par.printit)==0
0153
0154
0155 TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0156 end
0157
0158
0159 if dynamicscaling == 1 && (max(residuals.res1,residuals.res2) > par.tol2||(residuals.res1>residuals.res2*5)||(residuals.res2>residuals.res1*5))
0160 if isdynamicscaling == 0
0161 fprintf(' ****Starting dynamic scaling at iteration = %d ****\n', iter);
0162 isdynamicscaling = 1;
0163 end
0164 if mod(iter,dyn_scale_update_it)==0
0165 wrel_norm2 = wrel_norm2 * residuals.res2^(1/dyn_scale_update_it);
0166 wrel_norm1 = wrel_norm1 * residuals.res1^(1/dyn_scale_update_it);
0167 if wrel_norm2/wrel_norm1 > scaleratio
0168 theta = theta/(scalecorrection^2);
0169 TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0170 elseif wrel_norm1/wrel_norm2 > scaleratio
0171 theta = theta*(scalecorrection^2);
0172 TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0173 end
0174 lambda = sigma/sqrt(theta);
0175 elseif mod(iter,dyn_scale_update_it) == 1
0176
0177 wrel_norm2 = residuals.res2^(1/dyn_scale_update_it);
0178 wrel_norm1 = residuals.res1^(1/dyn_scale_update_it);
0179 else
0180
0181 wrel_norm2 = wrel_norm2 * residuals.res2^(1/dyn_scale_update_it);
0182 wrel_norm1 = wrel_norm1 * residuals.res1^(1/dyn_scale_update_it);
0183 end
0184 elseif isdynamicscaling == 1 && dynamicscaling == 1
0185 fprintf(' ****Stopping dynamic scaling at iteration = %d ****\n', iter);
0186 isdynamicscaling = 0;
0187 end
0188
0189 end
0190
0191 fprintf('\n\n ****End 2EBD-HPE algorithm ****\n');
0192
0193
0194
0195
0196
0197 end
0198
0199 function [lambda_a] = GetAdaptiveLambdaScaled(prob,lambda,sigma,x,x_tilde,v_1,u,u_tilde,v_2,gradf_x,theta)
0200
0201 eps = prob.error1(x,x_tilde,gradf_x);
0202
0203 lambda_sq = lambda^2;
0204
0205
0206 v_x = v_1;
0207 v_u = v_2;
0208 Norm_v_x_sq = norm(v_x)^2;
0209 Norm_v_u_sq = norm(v_u)^2;
0210
0211
0212 Delta_x = (x_tilde-x);
0213 Delta_u = (u_tilde-u);
0214 Norm_Delta_x = norm(Delta_x);
0215 Norm_Delta_u = norm(Delta_u);
0216 Norm_Delta_x_sq = Norm_Delta_x^2;
0217 Norm_Delta_u_sq = Norm_Delta_u^2;
0218
0219 tr_v_Delta_x = v_x'*Delta_x;
0220 tr_v_Delta_u = v_u'*Delta_u;
0221
0222 alpha = sigma^2*(Norm_Delta_x_sq/theta+Norm_Delta_u_sq);
0223
0224 a = (Norm_v_x_sq)/theta + Norm_v_u_sq;
0225 b = (lambda*Norm_v_x_sq+tr_v_Delta_x)/theta + lambda*Norm_v_u_sq + tr_v_Delta_u + eps;
0226 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;
0227 discriminant = sqrt(b^2-a*c);
0228 if b <= 0
0229 beta = (-b+discriminant)/a;
0230 else
0231 beta = -c/(b+discriminant);
0232 end
0233 if beta <= 1e-13
0234 beta = 0;
0235 else
0236
0237 end
0238 lambda_a = lambda + beta;
0239 end