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