📄 optim_private_trdog.c
字号:
/*
* MATLAB Compiler: 2.0.1
* Date: Tue May 08 21:28:22 2001
* Arguments: "-B" "sgl" "-m" "-W" "mainhg" "-L" "C" "asufit.m" "absfun.m"
* "absfunfree.m" "absfunshift.m" "absfunwidth.m" "dispfit.m" "dispfitdisp.m"
* "dispfitfree.m" "dispfitwidth.m" "seqmodfree.m" "spcfun.m" "spcfun1.m"
* "atamult.m" "aprecon.m"
*/
#include "optim_private_trdog.h"
#include "optim_private_pcgr.h"
#include "optim_private_quad1d.h"
#include "optim_private_trust.h"
/*
* The function "Moptim_private_trdog" is the implementation version of the
* "optim/private/trdog" M-function from file
* "C:\MATLABR11\toolbox\optim\private\trdog.m" (lines 1-226). It contains the
* actual compiled code for that M-function. It is a static function and must
* only be called from one of the interface functions, appearing below.
*/
/*
* function[s,snod,qpval,posdef,pcgit,Z] = trdog(x,g,H,fdata,D,delta,dv,...
*/
static mxArray * Moptim_private_trdog(mxArray * * snod,
mxArray * * qpval,
mxArray * * posdef,
mxArray * * pcgit,
mxArray * * Z,
int nargout_,
mxArray * x,
mxArray * g,
mxArray * H,
mxArray * fdata,
mxArray * D,
mxArray * delta,
mxArray * dv,
mxArray * mtxmpy,
mxArray * pcmtx,
mxArray * pcoptions,
mxArray * tol,
mxArray * kmax,
mxArray * theta,
mxArray * l,
mxArray * u,
mxArray * Z_,
mxArray * dnewt,
mxArray * preconflag) {
mxArray * s = mclGetUninitializedArray();
mxArray * DG = mclGetUninitializedArray();
mxArray * DM = mclGetUninitializedArray();
mxArray * HH = mclGetUninitializedArray();
mxArray * MM = mclGetUninitializedArray();
mxArray * R = mclGetUninitializedArray();
mxArray * W = mclGetUninitializedArray();
mxArray * WW = mclGetUninitializedArray();
mxArray * ZZ = mclGetUninitializedArray();
mxArray * alpha = mclGetUninitializedArray();
mxArray * arg = mclGetUninitializedArray();
mxArray * dis = mclGetUninitializedArray();
mxArray * fcnt = mclGetUninitializedArray();
mxArray * grad = mclGetUninitializedArray();
mxArray * ipt = mclGetUninitializedArray();
mxArray * lambda = mclGetUninitializedArray();
mxArray * mdis = mclGetUninitializedArray();
mxArray * mmdis = mclGetUninitializedArray();
mxArray * n = mclGetUninitializedArray();
mxArray * ng = mclGetUninitializedArray();
mxArray * ngrad = mclGetUninitializedArray();
mxArray * nrhs = mclGetUninitializedArray();
mxArray * nrmv2 = mclGetUninitializedArray();
mxArray * ns = mclGetUninitializedArray();
mxArray * nss = mclGetUninitializedArray();
mxArray * nst = mclGetUninitializedArray();
mxArray * nx = mclGetUninitializedArray();
mxArray * permR = mclGetUninitializedArray();
mxArray * po = mclGetUninitializedArray();
mxArray * qpval0 = mclGetUninitializedArray();
mxArray * qpval1 = mclGetUninitializedArray();
mxArray * qpval2 = mclGetUninitializedArray();
mxArray * qpval3 = mclGetUninitializedArray();
mxArray * r = mclGetUninitializedArray();
mxArray * rhs = mclGetUninitializedArray();
mxArray * sg = mclGetUninitializedArray();
mxArray * ss = mclGetUninitializedArray();
mxArray * ssave = mclGetUninitializedArray();
mxArray * ssg = mclGetUninitializedArray();
mxArray * sssave = mclGetUninitializedArray();
mxArray * ssssave = mclGetUninitializedArray();
mxArray * st = mclGetUninitializedArray();
mxArray * stsave = mclGetUninitializedArray();
mxArray * tau = mclGetUninitializedArray();
mxArray * tol2 = mclGetUninitializedArray();
mxArray * v1 = mclGetUninitializedArray();
mxArray * v2 = mclGetUninitializedArray();
mclValidateInputs(
"optim/private/trdog",
18,
&x,
&g,
&H,
&fdata,
&D,
&delta,
&dv,
&mtxmpy,
&pcmtx,
&pcoptions,
&tol,
&kmax,
&theta,
&l,
&u,
&Z_,
&dnewt,
&preconflag);
mclCopyArray(&tol);
mclCopyInputArg(Z, Z_);
/*
* mtxmpy,pcmtx,pcoptions,tol,kmax,theta,l,u,Z,dnewt,preconflag);
* %TRDOG Reflected (2-D) trust region trial step (box constraints)
* %
* % [s,snod,qpval,posdef,pcgit,Z] = TRDOG(x,g,H,fdata,D,delta,dv,...
* % mtxmpy,pcmtx,pcoptions,tol,theta,l,u,Z,dnewt,preconflag);
* %
* % Determine the trial step `s', an approx. trust region solution.
* % `s' is chosen as the best of 3 steps: the scaled gradient
* % (truncated to maintain strict feasibility),
* % a 2-D trust region solution (truncated to remain strictly feas.),
* % and the reflection of the 2-D trust region solution,
* % (truncated to remain strictly feasible).
* %
* % The 2-D subspace (defining the trust region
* % problem) is defined by the scaled gradient
* % direction and a CG process (returning
* % either an approximate Newton step of a dirn of negative curvature.
* % See ?? for more detail.
* % Driver function is SFMIN
*
* % Copyright (c) 1990-98 by The MathWorks, Inc.
* % $Revision: 1.5 $ $Date: 1998/09/15 22:43:40 $
*
* % Initialization
* n = length(g);
*/
mlfAssign(&n, mlfLength(g));
/*
* pcgit = 0;
*/
mlfAssign(pcgit, mlfScalar(0.0));
/*
* grad = D*g;
*/
mlfAssign(&grad, mlfMtimes(D, g));
/*
* DM = D;
*/
mlfAssign(&DM, D);
/*
* DG = sparse(1:n,1:n,full(abs(g).*dv));
*/
mlfAssign(
&DG,
mlfSparse(
mlfColon(mlfScalar(1.0), n, NULL),
mlfColon(mlfScalar(1.0), n, NULL),
mlfFull(mlfTimes(mlfAbs(g), dv)),
NULL,
NULL,
NULL));
/*
* posdef = 1;
*/
mlfAssign(posdef, mlfScalar(1.0));
/*
* pcgit = 0;
*/
mlfAssign(pcgit, mlfScalar(0.0));
/*
* tol2 = sqrt(eps);
*/
mlfAssign(&tol2, mlfSqrt(mlfEps()));
/*
* v1 = dnewt;
*/
mlfAssign(&v1, dnewt);
/*
* qpval1 = inf;
*/
mlfAssign(&qpval1, mlfInf());
/*
* qpval2 = inf;
*/
mlfAssign(&qpval2, mlfInf());
/*
* qpval3 = inf;
*/
mlfAssign(&qpval3, mlfInf());
/*
*
* % DETERMINE A 2-DIMENSIONAL SUBSPACE
* if isempty(Z)
*/
if (mlfTobool(mlfIsempty(*Z))) {
/*
* if isempty(v1)
*/
if (mlfTobool(mlfIsempty(v1))) {
/*
* switch preconflag
* case 'hessprecon'
*/
if (mclSwitchCompare(preconflag, mxCreateString("hessprecon"))) {
/*
* HH = DM*H*DM + DG;
*/
mlfAssign(&HH, mlfPlus(mlfMtimes(mlfMtimes(DM, H), DM), DG));
/*
* [R,permR] = feval(pcmtx,HH,pcoptions);
*/
mlfFeval(
mlfVarargout(&R, &permR, NULL),
mclFevalLookup(pcmtx, 0, NULL),
HH,
pcoptions,
NULL);
/*
* case 'jacobprecon'
*/
} else if (mclSwitchCompare(
preconflag, mxCreateString("jacobprecon"))) {
/*
* [R,permR] = feval(pcmtx,DM,DG,H,pcoptions);
*/
mlfFeval(
mlfVarargout(&R, &permR, NULL),
mclFevalLookup(pcmtx, 0, NULL),
DM,
DG,
H,
pcoptions,
NULL);
/*
* otherwise
*/
} else {
/*
* error('Invalid string used for PRECONFLAG argument to TRDOG');
*/
mlfError(
mxCreateString(
"Invalid string used for PRECONFLAG argument to TRDOG"));
/*
* end
*/
}
/*
* % We now pass kmax in from calling function
* %kmax = max(1,floor(n/2));
* if tol <= 0,
*/
if (mlfTobool(mlfLe(tol, mlfScalar(0.0)))) {
/*
* tol = .1;
*/
mlfAssign(&tol, mlfScalar(.1));
/*
* end
*/
}
/*
* [v1,posdef,pcgit] = pcgr(DM,DG,grad,kmax,tol,mtxmpy,fdata,H,R,permR);
*/
mlfAssign(
&v1,
mlfOptim_private_pcgr(
posdef,
pcgit,
DM,
DG,
grad,
kmax,
tol,
mtxmpy,
fdata,
H,
R,
permR));
/*
* end
*/
}
/*
* if norm(v1) > 0
*/
if (mlfTobool(mlfGt(mlfNorm(v1, NULL), mlfScalar(0.0)))) {
/*
* v1 = v1/norm(v1);
*/
mlfAssign(&v1, mlfMrdivide(v1, mlfNorm(v1, NULL)));
/*
* end
*/
}
/*
* Z(:,1) = v1;
*/
mlfIndexAssign(Z, "(?,?)", mlfCreateColonIndex(), mlfScalar(1.0), v1);
/*
* if n > 1
*/
if (mlfTobool(mlfGt(n, mlfScalar(1.0)))) {
/*
* if (posdef < 1)
*/
if (mlfTobool(mlfLt(*posdef, mlfScalar(1.0)))) {
/*
* v2 = D*sign(grad);
*/
mlfAssign(&v2, mlfMtimes(D, mlfSign(grad)));
/*
* if norm(v2) > 0
*/
if (mlfTobool(mlfGt(mlfNorm(v2, NULL), mlfScalar(0.0)))) {
/*
* v2 = v2/norm(v2);
*/
mlfAssign(&v2, mlfMrdivide(v2, mlfNorm(v2, NULL)));
/*
* end
*/
}
/*
* v2 = v2 - v1*(v1'*v2);
*/
mlfAssign(
&v2,
mlfMinus(
v2, mlfMtimes(v1, mlfMtimes(mlfCtranspose(v1), v2))));
/*
* nrmv2 = norm(v2);
*/
mlfAssign(&nrmv2, mlfNorm(v2, NULL));
/*
* if nrmv2 > tol2
*/
if (mlfTobool(mlfGt(nrmv2, tol2))) {
/*
* v2 = v2/nrmv2;
*/
mlfAssign(&v2, mlfMrdivide(v2, nrmv2));
/*
* Z(:,2) = v2;
*/
mlfIndexAssign(
Z, "(?,?)", mlfCreateColonIndex(), mlfScalar(2.0), v2);
/*
* end
*/
}
/*
* else
*/
} else {
/*
* if norm(grad) > 0
*/
if (mlfTobool(mlfGt(mlfNorm(grad, NULL), mlfScalar(0.0)))) {
/*
* v2 = grad/norm(grad);
*/
mlfAssign(&v2, mlfMrdivide(grad, mlfNorm(grad, NULL)));
/*
* else
*/
} else {
/*
* v2 = grad;
*/
mlfAssign(&v2, grad);
/*
* end
*/
}
/*
* v2 = v2 - v1*(v1'*v2);
*/
mlfAssign(
&v2,
mlfMinus(
v2, mlfMtimes(v1, mlfMtimes(mlfCtranspose(v1), v2))));
/*
* nrmv2 = norm(v2);
*/
mlfAssign(&nrmv2, mlfNorm(v2, NULL));
/*
* if nrmv2 > tol2
*/
if (mlfTobool(mlfGt(nrmv2, tol2))) {
/*
* v2 = v2/nrmv2;
*/
mlfAssign(&v2, mlfMrdivide(v2, nrmv2));
/*
* Z(:,2) = v2;
*/
mlfIndexAssign(
Z, "(?,?)", mlfCreateColonIndex(), mlfScalar(2.0), v2);
/*
* end
*/
}
/*
* end
*/
}
/*
* end
*/
}
/*
* end
*/
}
/*
*
* % REDUCE TO THE CHOSEN SUBSPACE
* W = DM*Z;
*/
mlfAssign(&W, mlfMtimes(DM, *Z));
/*
* WW = feval(mtxmpy,W,H,fdata);
*/
mlfAssign(
&WW,
mlfFeval(
mclValueVarargout(),
mclFevalLookup(mtxmpy, 0, NULL),
W,
H,
fdata,
NULL));
/*
* W = DM*WW;
*/
mlfAssign(&W, mlfMtimes(DM, WW));
/*
* MM = full(Z'*W + Z'*DG*Z);
*/
mlfAssign(
&MM,
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -