Skip to content

Commit

Permalink
Add missing dgelsy to Lapack functions.
Browse files Browse the repository at this point in the history
- Add dgelsy to the runtime, the Lapack module and to the function
  evaluation in both the old and the new frontend.
- Improve the function evalution in the new frontend so it can handle
  variables that depend on each other better.

Belonging to [master]:
  - OpenModelica/OMCompiler#2717
  • Loading branch information
perost authored and OpenModelica-Hudson committed Oct 11, 2018
1 parent 2904807 commit e645dd7
Show file tree
Hide file tree
Showing 6 changed files with 149 additions and 1 deletion.
29 changes: 29 additions & 0 deletions Compiler/FrontEnd/CevalFunction.mo
Expand Up @@ -724,6 +724,35 @@ algorithm
then
(cache, env);

case("dgelsy", {arg_M, arg_N, arg_NRHS, arg_A, arg_LDA, arg_B, arg_LDB,
arg_JPVT, arg_RCOND, arg_RANK, arg_WORK, arg_LWORK, arg_INFO},
cache, env)
equation
(M, cache) = evaluateExtIntArg(arg_M, cache, env);
(N, cache) = evaluateExtIntArg(arg_N, cache, env);
(NRHS, cache) = evaluateExtIntArg(arg_NRHS, cache, env);
(A, cache) = evaluateExtRealMatrixArg(arg_A, cache, env);
(LDA, cache) = evaluateExtIntArg(arg_LDA, cache, env);
(B, cache) = evaluateExtRealMatrixArg(arg_B, cache, env);
(LDB, cache) = evaluateExtIntArg(arg_LDB, cache, env);
(JPVT, cache) = evaluateExtIntArrayArg(arg_JPVT, cache, env);
(RCOND, cache) = evaluateExtRealArg(arg_RCOND, cache, env);
(WORK, cache) = evaluateExtRealArrayArg(arg_WORK, cache, env);
(LWORK, cache) = evaluateExtIntArg(arg_LWORK, cache, env);
(A, B, JPVT, RANK, WORK, INFO) =
Lapack.dgelsy(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, WORK, LWORK);
val_A = ValuesUtil.makeRealMatrix(A);
val_B = ValuesUtil.makeRealMatrix(B);
val_JPVT = ValuesUtil.makeIntArray(JPVT);
val_RANK = ValuesUtil.makeInteger(RANK);
val_WORK = ValuesUtil.makeRealArray(WORK);
val_INFO = ValuesUtil.makeInteger(INFO);
arg_out = {arg_A, arg_B, arg_JPVT, arg_RANK, arg_WORK, arg_INFO};
val_out = {val_A, val_B, val_JPVT, val_RANK, val_WORK, val_INFO};
(cache, env) = assignExtOutputs(arg_out, val_out, cache, env);
then
(cache, env);

case("dgesv", {arg_N, arg_NRHS, arg_A, arg_LDA, arg_IPIV, arg_B, arg_LDB,
arg_INFO},
cache, env)
Expand Down
19 changes: 18 additions & 1 deletion Compiler/NFFrontEnd/NFEvalFunction.mo
Expand Up @@ -193,13 +193,21 @@ protected
algorithm
repl := ReplTree.new();

// Add inputs to the replacement tree. Since they can't be assigned to the
// replacements don't need to be mutable.
for i in fn.inputs loop
arg :: rest_args := rest_args;
repl := addInputReplacement(i, "", arg, repl);
end for;

// Add outputs and local variables to the replacement tree. These do need to
// be mutable to allow assigning to them.
repl := List.fold(fn.outputs, function addMutableReplacement(prefix = ""), repl);
repl := List.fold(fn.locals, function addMutableReplacement(prefix = ""), repl);

// Apply the replacements to the replacements themselves. This is done after
// building the tree to make sure all the replacements are available.
repl := ReplTree.map(repl, function applyBindingReplacement(repl = repl));
end createReplacements;

function addMutableReplacement
Expand All @@ -211,7 +219,6 @@ protected
Expression repl_exp;
algorithm
repl_exp := getBindingExp(node, repl);
repl_exp := Expression.map(repl_exp, function applyReplacements2(repl = repl));
repl_exp := Expression.makeMutable(repl_exp);
repl := ReplTree.add(repl, prefix + InstNode.name(node), repl_exp);
end addMutableReplacement;
Expand Down Expand Up @@ -303,6 +310,15 @@ algorithm
repl := ReplTree.add(repl, prefix + InstNode.name(node), argument);
end addInputReplacement;

function applyBindingReplacement
input String name;
input Expression exp;
input ReplTree.Tree repl;
output Expression outExp;
algorithm
outExp := Expression.map(exp, function applyReplacements2(repl = repl));
end applyBindingReplacement;

function applyReplacements
input ReplTree.Tree repl;
input output list<Statement> fnBody;
Expand Down Expand Up @@ -987,6 +1003,7 @@ algorithm
case "dgegv" algorithm EvalFunctionExt.Lapack_dgegv(args); then ();
case "dgels" algorithm EvalFunctionExt.Lapack_dgels(args); then ();
case "dgelsx" algorithm EvalFunctionExt.Lapack_dgelsx(args); then ();
case "dgelsy" algorithm EvalFunctionExt.Lapack_dgelsy(args); then ();
case "dgesv" algorithm EvalFunctionExt.Lapack_dgesv(args); then ();
case "dgglse" algorithm EvalFunctionExt.Lapack_dgglse(args); then ();
case "dgtsv" algorithm EvalFunctionExt.Lapack_dgtsv(args); then ();
Expand Down
35 changes: 35 additions & 0 deletions Compiler/NFFrontEnd/NFEvalFunctionExt.mo
Expand Up @@ -181,6 +181,41 @@ algorithm
assignVariable(info, Expression.makeInteger(INFO));
end Lapack_dgelsx;

function Lapack_dgelsy
input list<Expression> args;
protected
Expression m, n, nrhs, a, lda, b, ldb, jpvt, rcond, rank, work, lwork, info;
Integer M, N, NRHS, LDA, LDB, RANK, LWORK, INFO;
list<list<Real>> A, B;
list<Integer> JPVT;
Real RCOND;
list<Real> WORK;
algorithm
{m, n, nrhs, a, lda, b, ldb, jpvt, rcond, rank, work, lwork, info} := args;

M := evaluateExtIntArg(m);
N := evaluateExtIntArg(n);
NRHS := evaluateExtIntArg(nrhs);
A := evaluateExtRealMatrixArg(a);
LDA := evaluateExtIntArg(lda);
B := evaluateExtRealMatrixArg(b);
LDB := evaluateExtIntArg(ldb);
JPVT := evaluateExtIntArrayArg(jpvt);
RCOND := evaluateExtRealArg(rcond);
WORK := evaluateExtRealArrayArg(work);
LWORK := evaluateExtIntArg(lwork);

(A, B, JPVT, RANK, WORK, INFO) :=
Lapack.dgelsy(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, WORK, LWORK);

assignVariableExt(a, Expression.makeRealMatrix(A));
assignVariableExt(b, Expression.makeRealMatrix(B));
assignVariable(jpvt, Expression.makeIntegerArray(JPVT));
assignVariable(rank, Expression.makeInteger(RANK));
assignVariable(work, Expression.makeRealArray(WORK));
assignVariable(info, Expression.makeInteger(INFO));
end Lapack_dgelsy;

function Lapack_dgesv
input list<Expression> args;
protected
Expand Down
1 change: 1 addition & 0 deletions Compiler/Script/CevalScript.mo
Expand Up @@ -2610,6 +2610,7 @@ algorithm
case "dgegv" then ();
case "dgels" then ();
case "dgelsx" then ();
case "dgelsy" then ();
case "dgeqpf" then ();
case "dgesv" then ();
case "dgesvd" then ();
Expand Down
23 changes: 23 additions & 0 deletions Compiler/Util/Lapack.mo
Expand Up @@ -128,6 +128,29 @@ public function dgelsx
annotation(Library = {"omcruntime", "Lapack"});
end dgelsx;

public function dgelsy
input Integer inM;
input Integer inN;
input Integer inNRHS;
input list<list<Real>> inA;
input Integer inLDA;
input list<list<Real>> inB;
input Integer inLDB;
input list<Integer> inJPVT;
input Real inRCOND;
input list<Real> inWORK;
input Integer inLWORK;
output list<list<Real>> outA;
output list<list<Real>> outB;
output list<Integer> outJPVT;
output Integer outRANK;
output list<Real> outWORK;
output Integer outINFO;
external "C" LapackImpl__dgelsy(inM, inN, inNRHS, inA, inLDA, inB, inLDB,
inJPVT, inRCOND, inWORK, inLWORK, outA, outB, outJPVT, outRANK, outWORK, outINFO)
annotation(Library = {"omcruntime", "Lapack"});
end dgelsy;

public function dgesv
input Integer inN;
input Integer inNRHS;
Expand Down
43 changes: 43 additions & 0 deletions Compiler/runtime/lapackimpl.c
Expand Up @@ -100,6 +100,10 @@ extern int dgelsx_(integer *m, integer *n, integer *nrhs, doublereal *a,
integer *lda, doublereal *b, integer *ldb, integer *jpvt, doublereal *rcond,
integer *rank, doublereal *work, integer *info);

extern int dgelsy_(integer *m, integer *n, integer *nrhs, doublereal *a,
integer *lda, doublereal *b, integer *ldb, integer *jpvt, doublereal *rcond,
integer *rank, doublereal *work, integer *lwork, integer *info);

extern int dgesv_(integer *n, integer *nrhs, doublereal *a, integer *lda,
integer *ipiv, doublereal *b, integer *ldb, integer *info);

Expand Down Expand Up @@ -467,6 +471,45 @@ void LapackImpl__dgelsx(int M, int N, int NRHS, void *inA, int LDA,
#endif
}

void LapackImpl__dgelsy(int M, int N, int NRHS, void *inA, int LDA,
void *inB, int LDB, void *inJPVT, double rcond, void *inWORK, int LWORK,
void **outA, void **outB, void **outJPVT, int *RANK, void **outWORK, int *INFO)
{
#ifndef NO_LAPACK
integer m, n, nrhs, lda, ldb, rank = 0, info = 0, lwork;
double *a, *b, *work;
integer *jpvt;

m = M;
n = N;
nrhs = NRHS;
lda = LDA;
ldb = LDB;
lwork = LWORK;

a = alloc_real_matrix(lda, n, inA);
b = alloc_real_matrix(ldb, nrhs, inB);
work = alloc_real_vector(lwork, inWORK);
jpvt = alloc_int_vector(n, inJPVT);

dgelsy_(&m, &n, &nrhs, a, &lda, b, &ldb, jpvt, &rcond, &rank, work, &lwork, &info);

*outA = mk_rml_real_matrix(lda, n, a);
*outB = mk_rml_real_matrix(lda, nrhs, b);
*outJPVT = mk_rml_int_vector(n, jpvt);
*RANK = rank;
*outWORK = mk_rml_real_vector(lwork, work);
*INFO = info;

free(a);
free(b);
free(work);
free(jpvt);
#else
MMC_THROW();
#endif
}

void LapackImpl__dgesv(int N, int NRHS, void *inA, int LDA, void *inB,
int LDB, void **outA, void **IPIV, void **outB, int *INFO)
{
Expand Down

0 comments on commit e645dd7

Please sign in to comment.