📄 aprecon.c
字号:
/*
* MATLAB Compiler: 2.0.1
* Date: Tue May 08 21:28:16 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 "aprecon.h"
/*
* The function "Maprecon" is the implementation version of the "aprecon"
* M-function from file "C:\MATLABR11\toolbox\optim\aprecon.m" (lines 1-93). 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[RPCMTX,ppvec] = aprecon(DM,DG,A,upperbandw);
*/
static mxArray * Maprecon(mxArray * * ppvec,
int nargout_,
mxArray * DM,
mxArray * DG,
mxArray * A,
mxArray * upperbandw) {
mxArray * RPCMTX = mclGetUninitializedArray();
mxArray * DDG = mclGetUninitializedArray();
mxArray * M = mclGetUninitializedArray();
mxArray * TM = mclGetUninitializedArray();
mxArray * d = mclGetUninitializedArray();
mxArray * dgg = mclGetUninitializedArray();
mxArray * dnrms = mclGetUninitializedArray();
mxArray * epsi = mclGetUninitializedArray();
mxArray * info = mclGetUninitializedArray();
mxArray * lambda = mclGetUninitializedArray();
mxArray * m = mclGetUninitializedArray();
mxArray * mdiag = mclGetUninitializedArray();
mxArray * n = mclGetUninitializedArray();
mxArray * nargin_ = mclGetUninitializedArray();
mxArray * p = mclGetUninitializedArray();
mlfAssign(&nargin_, mlfNargin(0, DM, DG, A, upperbandw, NULL));
mclValidateInputs("aprecon", 4, &DM, &DG, &A, &upperbandw);
mclCopyArray(&upperbandw);
/*
* %APRECON Banded preconditioner function for least-squares problems.
* %
* % [RPCMTX,PPVEC] = APRECON(DM,DG,A,UPPERBW) produces the sparse
* % nonsingular upper triangular matrix RPCMTX such that
* % RPCMTX'*RPCMTX is a preconditioner for the
* % matrix M = DM*(A'*A)*DM + DG, where DM is a positive
* % diagonal matrix, DG is a non-negative diagonal matrix,
* % and A is sparse rectangular matrix with more rows than columns.
* % PPVEC is the associated permutation (row) vector and UPPERBW
* % specifies the upperbandwidth of RPCMTX.
*
* % Copyright (c) 1990-98 by The MathWorks, Inc.
* % $Revision: 1.2 $ $Date: 1998/07/21 22:50:21 $
*
* % Initialization
* if nargin < 3,
*/
if (mlfTobool(mlfLt(nargin_, mlfScalar(3.0)))) {
/*
* error('aprecon requires at least 3 input parameters.'),
*/
mlfError(
mxCreateString("aprecon requires at least 3 input parameters."));
/*
* end
*/
}
/*
* [m,n] = size(A);
*/
mlfSize(mlfVarargout(&m, &n, NULL), A, NULL);
/*
* if m < n
*/
if (mlfTobool(mlfLt(m, n))) {
/*
* error('A matrix must have at least as many rows as columns.')
*/
mlfError(
mxCreateString(
"A matrix must have at least as many rows as columns."));
/*
* end
*/
}
/*
* n = length(DM);
*/
mlfAssign(&n, mlfLength(DM));
/*
* RPCMTX = speye(n);
*/
mlfAssign(&RPCMTX, mlfSpeye(n, NULL));
/*
* ppvec = (1:n);
*/
mlfAssign(ppvec, mlfColon(mlfScalar(1.0), n, NULL));
/*
* if isempty(A)
*/
if (mlfTobool(mlfIsempty(A))) {
/*
* RPCMTX = sparse([]);
*/
mlfAssign(
&RPCMTX,
mlfSparse(mclCreateEmptyArray(), NULL, NULL, NULL, NULL, NULL));
/*
* ppvec=[];
*/
mlfAssign(ppvec, mclCreateEmptyArray());
/*
* return
*/
goto return_;
/*
* end
*/
}
/*
* epsi = sqrt(eps);
*/
mlfAssign(&epsi, mlfSqrt(mlfEps()));
/*
* if nargin < 4
*/
if (mlfTobool(mlfLt(nargin_, mlfScalar(4.0)))) {
/*
* upperbandw = 0;
*/
mlfAssign(&upperbandw, mlfScalar(0.0));
/*
* elseif isempty(upperbandw)
*/
} else if (mlfTobool(mlfIsempty(upperbandw))) {
/*
* upperbandw = 0;
*/
mlfAssign(&upperbandw, mlfScalar(0.0));
/*
* end
*/
}
/*
*
* % Form matrix M
* TM = A*DM;
*/
mlfAssign(&TM, mlfMtimes(A, DM));
/*
*
* % Determine factor of preconditioner.
* if upperbandw == 0 % Diag. preconditioner based on column 2-norms
*/
if (mlfTobool(mlfEq(upperbandw, mlfScalar(0.0)))) {
/*
* M = TM'*TM + DG;
*/
mlfAssign(&M, mlfPlus(mlfMtimes(mlfCtranspose(TM), TM), DG));
/*
* % M cannot be a row vector so no problem with sum.
* dnrms = sqrt(sum(M.*M))';
*/
mlfAssign(&dnrms, mlfCtranspose(mlfSqrt(mlfSum(mlfTimes(M, M), NULL))));
/*
* d = max(sqrt(dnrms),epsi);
*/
mlfAssign(&d, mlfMax(NULL, mlfSqrt(dnrms), epsi, NULL));
/*
* RPCMTX = sparse(1:n,1:n,full(d));
*/
mlfAssign(
&RPCMTX,
mlfSparse(
mlfColon(mlfScalar(1.0), n, NULL),
mlfColon(mlfScalar(1.0), n, NULL),
mlfFull(d),
NULL,
NULL,
NULL));
/*
* ppvec = (1:n);
*/
mlfAssign(ppvec, mlfColon(mlfScalar(1.0), n, NULL));
/*
* % But what if RPCMTX is singular?
* elseif upperbandw >= n-1 % Attempt sparse QR
*/
} else if (mlfTobool(mlfGe(upperbandw, mlfMinus(n, mlfScalar(1.0))))) {
/*
* dgg = sqrt(diag(DG));
*/
mlfAssign(&dgg, mlfSqrt(mlfDiag(DG, NULL)));
/*
* dgg = full(dgg);
*/
mlfAssign(&dgg, mlfFull(dgg));
/*
* DDG = sparse(1:n,1:n,dgg);
*/
mlfAssign(
&DDG,
mlfSparse(
mlfColon(mlfScalar(1.0), n, NULL),
mlfColon(mlfScalar(1.0), n, NULL),
dgg,
NULL,
NULL,
NULL));
/*
* TM = [TM;DDG];
*/
mlfAssign(
&TM, mlfVertcat(mlfHorzcat(TM, NULL), mlfHorzcat(DDG, NULL), NULL));
/*
* p = colmmd(TM);
*/
mlfAssign(&p, mlfColmmd(TM));
/*
* RPCMTX = qr(TM(:,p));
*/
mlfAssign(
&RPCMTX,
mlfQr(
NULL,
NULL,
mlfIndexRef(TM, "(?,?)", mlfCreateColonIndex(), p),
NULL,
NULL));
/*
* RPCMTX = RPCMTX(1:n,1:n);
*/
mlfAssign(
&RPCMTX,
mlfIndexRef(
RPCMTX,
"(?,?)",
mlfColon(mlfScalar(1.0), n, NULL),
mlfColon(mlfScalar(1.0), n, NULL)));
/*
* ppvec = p;
*/
mlfAssign(ppvec, p);
/*
*
* % Modify for singularity?
* mdiag = min(abs(diag(RPCMTX)));
*/
mlfAssign(
&mdiag, mlfMin(NULL, mlfAbs(mlfDiag(RPCMTX, NULL)), NULL, NULL));
/*
* lambda = 1;
*/
mlfAssign(&lambda, mlfScalar(1.0));
/*
* while mdiag < sqrt(eps);
*/
while (mlfTobool(mlfLt(mdiag, mlfSqrt(mlfEps())))) {
/*
* TM = [A*DM; DDG + lambda*speye(n)];
*/
mlfAssign(
&TM,
mlfVertcat(
mlfHorzcat(mlfMtimes(A, DM), NULL),
mlfHorzcat(
mlfPlus(DDG, mlfMtimes(lambda, mlfSpeye(n, NULL))), NULL),
NULL));
/*
* lambda = 4*lambda;
*/
mlfAssign(&lambda, mlfMtimes(mlfScalar(4.0), lambda));
/*
* p = colmmd(TM);
*/
mlfAssign(&p, mlfColmmd(TM));
/*
* RPCMTX = qr(TM(:,p));
*/
mlfAssign(
&RPCMTX,
mlfQr(
NULL,
NULL,
mlfIndexRef(TM, "(?,?)", mlfCreateColonIndex(), p),
NULL,
NULL));
/*
* RPCMTX = RPCMTX(1:n,1:n);
*/
mlfAssign(
&RPCMTX,
mlfIndexRef(
RPCMTX,
"(?,?)",
mlfColon(mlfScalar(1.0), n, NULL),
mlfColon(mlfScalar(1.0), n, NULL)));
/*
* ppvec = p;
*/
mlfAssign(ppvec, p);
/*
* mdiag = min(abs(diag(RPCMTX)));
*/
mlfAssign(
&mdiag, mlfMin(NULL, mlfAbs(mlfDiag(RPCMTX, NULL)), NULL, NULL));
/*
* end
*/
}
/*
*
* elseif (upperbandw > 0) & (upperbandw < n-1) % Band approximation.
*/
} else {
mxArray * a_ = mclInitialize(mlfGt(upperbandw, mlfScalar(0.0)));
if (mlfTobool(a_)
&& mlfTobool(
mlfAnd(a_, mlfLt(upperbandw, mlfMinus(n, mlfScalar(1.0)))))) {
mxDestroyArray(a_);
/*
* M = TM'*TM + DG;
*/
mlfAssign(&M, mlfPlus(mlfMtimes(mlfCtranspose(TM), TM), DG));
/*
* p = (1:n);
*/
mlfAssign(&p, mlfColon(mlfScalar(1.0), n, NULL));
/*
* M = tril(triu(M(p,p),-upperbandw),upperbandw);
*/
mlfAssign(
&M,
mlfTril(
mlfTriu(mlfIndexRef(M, "(?,?)", p, p), mlfUminus(upperbandw)),
upperbandw));
/*
* RPCMTX = sparse(n,n);
*/
mlfAssign(&RPCMTX, mlfSparse(n, n, NULL, NULL, NULL, NULL));
/*
* [RPCMTX,info] = chol(M);
*/
mlfAssign(&RPCMTX, mlfChol(&info, M));
/*
* lambda = 1;
*/
mlfAssign(&lambda, mlfScalar(1.0));
/*
* while info > 0
*/
while (mlfTobool(mlfGt(info, mlfScalar(0.0)))) {
/*
* M = M + lambda*speye(n);
*/
mlfAssign(&M, mlfPlus(M, mlfMtimes(lambda, mlfSpeye(n, NULL))));
/*
* RPCMTX = sparse(n,n);
*/
mlfAssign(&RPCMTX, mlfSparse(n, n, NULL, NULL, NULL, NULL));
/*
* [RPCMTX,info] = chol(M);
*/
mlfAssign(&RPCMTX, mlfChol(&info, M));
/*
* lambda = 4*lambda;
*/
mlfAssign(&lambda, mlfMtimes(mlfScalar(4.0), lambda));
/*
* end
*/
}
/*
* ppvec = p;
*/
mlfAssign(ppvec, p);
/*
* else
*/
} else {
mxDestroyArray(a_);
/*
* error('upperbandw must be >= 0.');
*/
mlfError(mxCreateString("upperbandw must be >= 0."));
}
/*
* end
*/
}
/*
*
*
*/
return_:
mclValidateOutputs("aprecon", 2, nargout_, &RPCMTX, ppvec);
mxDestroyArray(DDG);
mxDestroyArray(M);
mxDestroyArray(TM);
mxDestroyArray(d);
mxDestroyArray(dgg);
mxDestroyArray(dnrms);
mxDestroyArray(epsi);
mxDestroyArray(info);
mxDestroyArray(lambda);
mxDestroyArray(m);
mxDestroyArray(mdiag);
mxDestroyArray(n);
mxDestroyArray(nargin_);
mxDestroyArray(p);
mxDestroyArray(upperbandw);
return RPCMTX;
}
/*
* The function "mlfAprecon" contains the normal interface for the "aprecon"
* M-function from file "C:\MATLABR11\toolbox\optim\aprecon.m" (lines 1-93).
* This function processes any input arguments and passes them to the
* implementation version of the function, appearing above.
*/
mxArray * mlfAprecon(mxArray * * ppvec,
mxArray * DM,
mxArray * DG,
mxArray * A,
mxArray * upperbandw) {
int nargout = 1;
mxArray * RPCMTX = mclGetUninitializedArray();
mxArray * ppvec__ = mclGetUninitializedArray();
mlfEnterNewContext(1, 4, ppvec, DM, DG, A, upperbandw);
if (ppvec != NULL) {
++nargout;
}
RPCMTX = Maprecon(&ppvec__, nargout, DM, DG, A, upperbandw);
mlfRestorePreviousContext(1, 4, ppvec, DM, DG, A, upperbandw);
if (ppvec != NULL) {
mclCopyOutputArg(ppvec, ppvec__);
} else {
mxDestroyArray(ppvec__);
}
return mlfReturnValue(RPCMTX);
}
/*
* The function "mlxAprecon" contains the feval interface for the "aprecon"
* M-function from file "C:\MATLABR11\toolbox\optim\aprecon.m" (lines 1-93).
* The feval function calls the implementation version of aprecon through this
* function. This function processes any input arguments and passes them to the
* implementation version of the function, appearing above.
*/
void mlxAprecon(int nlhs, mxArray * plhs[], int nrhs, mxArray * prhs[]) {
mxArray * mprhs[4];
mxArray * mplhs[2];
int i;
if (nlhs > 2) {
mlfError(
mxCreateString(
"Run-time Error: File: aprecon Line: 1 Column: "
"0 The function \"aprecon\" was called with mor"
"e than the declared number of outputs (2)"));
}
if (nrhs > 4) {
mlfError(
mxCreateString(
"Run-time Error: File: aprecon Line: 1 Column:"
" 0 The function \"aprecon\" was called with m"
"ore than the declared number of inputs (4)"));
}
for (i = 0; i < 2; ++i) {
mplhs[i] = NULL;
}
for (i = 0; i < 4 && i < nrhs; ++i) {
mprhs[i] = prhs[i];
}
for (; i < 4; ++i) {
mprhs[i] = NULL;
}
mlfEnterNewContext(0, 4, mprhs[0], mprhs[1], mprhs[2], mprhs[3]);
mplhs[0]
= Maprecon(&mplhs[1], nlhs, mprhs[0], mprhs[1], mprhs[2], mprhs[3]);
mlfRestorePreviousContext(0, 4, mprhs[0], mprhs[1], mprhs[2], mprhs[3]);
plhs[0] = mplhs[0];
for (i = 1; i < 2 && i < nlhs; ++i) {
plhs[i] = mplhs[i];
}
for (; i < 2; ++i) {
mxDestroyArray(mplhs[i]);
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -