Skip to content

Commit

Permalink
workaround fix for discrete relations in cpp tempalte
Browse files Browse the repository at this point in the history 
reactivated jacobians in cpp template


git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@19063 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
  • Loading branch information
@niklwors
niklwors committed Feb 12, 2014
1 parent 16aeaf9 commit a6179af
Show file tree
Hide file tree
Showing 4 changed files with 195 additions and 54 deletions.
108 changes: 59 additions & 49 deletions Compiler/Template/CodegenCpp.tpl
View file
Expand Up @@ -2776,10 +2776,10 @@ case SIMCODE(modelInfo = MODELINFO(__)) then
bool isConsistent();
/*
void initialAnalyticJacobian();
void calcJacobianColumn();
*/
//Variables:
EventHandling _event_handling;
Expand Down Expand Up @@ -8023,7 +8023,10 @@ template giveZeroFunc3(Integer index1, Exp relation, Text &varDecls /*BUFP*/,Tex
f[<%index1%>]=(<%e1%> - EPSILON - <%e2%>);
>>
else
error(sourceInfo(), 'Unknown relation: <%printExpStr(rel)%> for <%index1%>')
<<
f[<%index1%>]=-1;
/*error(sourceInfo(), 'Unknown relation: <%printExpStr(rel)%> for <%index1%>')*/
>>
end match
case CALL(path=IDENT(name="sample"), expLst={_, start, interval}) then
//error(sourceInfo(), ' sample not supported for <%index1%> ')
Expand Down Expand Up @@ -8423,21 +8426,6 @@ case "A" then
}
>>
case _ then
/* wbraun: not used anywhere
let sp_size_index = lengthListElements(sparsepattern)
let leadindex = ( sparsepattern |> (indexes) hasindex index0 =>
if index0 then
<<data->simulationInfo.analyticJacobians[index].sparsePattern.leadindex[<%index0%>] = data->simulationInfo.analyticJacobians[index].sparsePattern.leadindex[<%intSub(index0,1)%>] + <%listLength(indexes)%>;>>
else
<<data->simulationInfo.analyticJacobians[index].sparsePattern.leadindex[<%index0%>] = <%listLength(indexes)%>;>>
;separator="\n")
let indexElems = ( flatten(sparsepattern) |> (indexes) hasindex index0 =>
<<data->simulationInfo.analyticJacobians[index].sparsePattern.index[<%index0%>] = <%indexes%>; >>
;separator="\n")
let colorArray = ( colorList |> (indexes) hasindex index0 =>
<<data->simulationInfo.analyticJacobians[index].sparsePattern.colorCols[<%index0%>] = <%indexes%>; >>
;separator="\n")
*/
let sp_size_index = lengthListElements(splitTuple212List(sparsepattern))
let indexColumn = (jacobianColumn |> (eqs,vars,indxColumn) => indxColumn;separator="\n")
let tmpvarsSize = (jacobianColumn |> (_,vars,_) => listLength(vars);separator="\n")
Expand All @@ -8461,8 +8449,7 @@ end initialAnalyticJacobians;
template functionAnalyticJacobians(list<JacobianMatrix> JacobianMatrixes, SimCode simCode)
"Generates Matrixes for Linear Model."
::=
/*
let initialjacMats = (JacobianMatrixes |> (mat, vars, name, (sparsepattern,(_,_)), colorList, _) hasindex index0 =>
let initialjacMats = (JacobianMatrixes |> (mat, vars, name, (sparsepattern,(_,_)), colorList, _) hasindex index0 =>
initialAnalyticJacobians(index0, mat, vars, name, sparsepattern, colorList,simCode)
;separator="\n\n";empty)
let jacMats = (JacobianMatrixes |> (mat, vars, name, (sparsepattern,(_,_)), colorList, maxColor) hasindex index0 =>
Expand All @@ -8472,15 +8459,7 @@ template functionAnalyticJacobians(list<JacobianMatrix> JacobianMatrixes, SimCod
<%initialjacMats%>
<%jacMats%>
>>
*/
match simCode
case SIMCODE(modelInfo = MODELINFO(__)) then
let classname = lastIdentOfPath(modelInfo.name)
<<
void <%classname%>::getJacobian(SparseMatrix& matrix)
{
}
>>

end functionAnalyticJacobians;


Expand All @@ -8504,7 +8483,7 @@ template functionJac(list<SimEqSystem> jacEquations, list<SimVar> tmpVars, Strin
<<
void <%classname%>::calcJacobianColumn()
{
<%varDecls%>
<%eqns_%>
}
Expand Down Expand Up @@ -8559,7 +8538,7 @@ case _ then

let jaccol = ( indexes |> i_index =>
//' _jacobian(<%index0%>,<%intSub(cref(i_index),1)%>) = _jac_y(<%intSub(cref(i_index),1)%>);'
' _jacobian(<%index0%>,<%cref(i_index)%><%matrixname%>$indexdiff) = _jac_y(<%cref(i_index)%><%matrixname%>$indexdiff);'
' _jacobian(<%index0%>,<%crefWithoutIndexOperator(i_index)%><%matrixname%>$indexdiff) = _jac_y(<%crefWithoutIndexOperator(i_index)%><%matrixname%>$indexdiff);'
;separator="\n" )
' _jac_x(<%index0%>)=1;
calcJacobianColumn();
Expand All @@ -8568,6 +8547,7 @@ _jac_x.clear();
;separator="\n")

<<
/*test*/
<%jacMats%>
void <%classname%>::getJacobian(SparseMatrix& matrix)
{
Expand All @@ -8580,7 +8560,7 @@ end generateMatrix;
template variableDefinitionsJacobians(list<JacobianMatrix> JacobianMatrixes)
"Generates defines for jacobian vars."
::=
/*

let analyticVars = (JacobianMatrixes |> (jacColumn, seedVars, name, (_,(diffVars,diffedVars)), _, _) hasindex index0 =>
let varsDef = variableDefinitionsJacobians2(index0, jacColumn, seedVars, name)
let sparseDef = defineSparseIndexes(diffVars, diffedVars, name)
Expand All @@ -8591,11 +8571,10 @@ template variableDefinitionsJacobians(list<JacobianMatrix> JacobianMatrixes)
;separator="\n";empty)

<<
/* Jacobian Variables */
<%analyticVars%>
/* Jacobian Variables */
<%analyticVars%>
>>
*/
""

end variableDefinitionsJacobians;

template variableDefinitionsJacobians2(Integer indexJacobian, list<JacobianColumn> jacobianColumn, list<SimVar> seedVars, String name)
Expand Down Expand Up @@ -8628,11 +8607,11 @@ case "jacobianVars" then
match index
case -1 then
<<
#define <%cref(name)%> _jac_tmp(<%index0%>)
#define <%crefWithoutIndexOperator(name)%> _jac_tmp(<%index0%>)
>>
case _ then
<<
#define <%cref(name)%> _jac_y(<%index%>)
#define <%crefWithoutIndexOperator(name)%> _jac_y(<%index%>)
>>
end match
end match
Expand All @@ -8641,20 +8620,12 @@ case "jacobianVarsSeed" then
case SIMVAR(aliasvar=NOALIAS()) then
let tmp = System.tmpTick()
<<
#define _<%jacobianVarsSeedDefine(name)%>$pDER<%matrixName%>$P<%jacobianVarsSeedDefine(name)%> _jac_x(<%index0%>)
#define <%crefWithoutIndexOperator(name)%>$pDER<%matrixName%><%crefWithoutIndexOperator(name)%> _jac_x(<%index0%>)
>>
end match
end jacobianVarDefine;

template jacobianVarsSeedDefine(ComponentRef cr)
::=
match cr
case CREF_IDENT(__) then '<%ident%>'
case CREF_QUAL(__) then '<%ident%><%subscriptsToCStrForArray(subscriptLst)%>$P<%crefToCStr1(componentRef)%>'

case WILD(__) then ' '
else "CREF_NOT_IDENT_OR_QUAL"
end jacobianVarsSeedDefine;


template defineSparseIndexes(list<SimVar> diffVars, list<SimVar> diffedVars, String matrixName) "template variableDefinitionsJacobians2
Expand All @@ -8663,12 +8634,12 @@ template defineSparseIndexes(list<SimVar> diffVars, list<SimVar> diffedVars, Str
match matrixName
case "A" then
let diffVarsResult = (diffVars |> var as SIMVAR(name=name) hasindex index0 =>
'#define <%cref(name)%><%matrixName%>$indexdiff <%index0%>'
'#define <%crefWithoutIndexOperator(name)%><%matrixName%>$indexdiff <%index0%>'
;separator="\n")
/* generate at least one print command to have the same index and avoid the strange side effect */
<<
/* <%matrixName%> sparse indexes */
/*<%diffVarsResult%>*/
<%diffVarsResult%>
>>
else " "
end defineSparseIndexes;
Expand Down Expand Up @@ -9262,5 +9233,44 @@ match simVar
end setVariablesDefault;



template crefWithoutIndexOperator(ComponentRef cr)
"Generates C equivalent name for component reference."
::=
match cr
case CREF_IDENT(ident = "xloc") then crefStr(cr)
case CREF_IDENT(ident = "time") then "time"
case WILD(__) then ''
else "$P" + crefToCStrfWithoutIndexOperator(cr)
end crefWithoutIndexOperator;

template crefToCStrfWithoutIndexOperator(ComponentRef cr)
"Helper function to cref."
::=
match cr
case CREF_IDENT(__) then '<%unquoteIdentifier(ident)%><%subscriptsToCStrWithoutIndexOperator(subscriptLst)%>'
case CREF_QUAL(__) then '<%unquoteIdentifier(ident)%><%subscriptsToCStrWithoutIndexOperator(subscriptLst)%>$P<%crefToCStrWithoutIndexOperator(componentRef)%>'
case WILD(__) then ''
else "CREF_NOT_IDENT_OR_QUAL"
end crefToCStrfWithoutIndexOperator;

template subscriptsToCStrWithoutIndexOperator(list<Subscript> subscripts)
::=
if subscripts then
'$lB<%subscripts |> s => subscriptToCStr(s) ;separator="$c"%>$rB'
end subscriptsToCStrWithoutIndexOperator;


template crefToCStrWithoutIndexOperator(ComponentRef cr)
"Helper function to cref."
::=
match cr
case CREF_IDENT(__) then '<%unquoteIdentifier(ident)%><%subscriptsToCStrWithoutIndexOperator(subscriptLst)%>'
case CREF_QUAL(__) then '<%unquoteIdentifier(ident)%><%subscriptsToCStrWithoutIndexOperator(subscriptLst)%>$P<%crefToCStrWithoutIndexOperator(componentRef)%>'
case WILD(__) then ''
else "CREF_NOT_IDENT_OR_QUAL"
end crefToCStrWithoutIndexOperator;


end CodegenCpp;

132 changes: 132 additions & 0 deletions SimulationRuntime/cpp/Core/Math/Functions.cpp
View file
Expand Up @@ -2,6 +2,27 @@
#include <Math/Functions.h>
#include <stdexcept>

/* Matrixes using column major order (as in Fortran) */
#ifndef set_matrix_elt
#define set_matrix_elt(A,r,c,n_rows,value) A[r + n_rows * c] = value
#endif

#ifndef get_matrix_elt
#define get_matrix_elt(A,r,c,n_rows) A[r + n_rows * c]
#endif


/* Matrixes using column major order (as in Fortran) */
/* colInd, rowInd, n_rows is added implicitly, makes code easier to read but may be considered bad programming style! */
#define set_pivot_matrix_elt(A,r,c,value) set_matrix_elt(A,rowInd[r],colInd[c],n_rows,value)
/* #define set_pivot_matrix_elt(A,r,c,value) set_matrix_elt(A,colInd[c],rowInd[r],n_cols,value) */
#define get_pivot_matrix_elt(A,r,c) get_matrix_elt(A,rowInd[r],colInd[c],n_rows)
/* #define get_pivot_matrix_elt(A,r,c) get_matrix_elt(A,colInd[c],rowInd[r],n_cols) */
#define swap(a,b) { int _swap=a; a=b; b=_swap; }

#ifndef min
#define min(a,b) ((a > b) ? (b) : (a))
#endif

double division (const double &a,const double &b,std::string text)
{
Expand All @@ -13,3 +34,114 @@ double division (const double &a,const double &b,std::string text)
throw std::invalid_argument(error_msg+text);
}
}

/*
find the maximum element below (and including) line/row start
*/
int maxsearch( double *A, int start, int n_rows, int n_cols, int *rowInd, int *colInd, int *maxrow, int *maxcol, double *maxabsval)
{
/* temporary variables */
int row;
int col;

/* Initialization */
int mrow = -1;
int mcol = -1;
double mabsval = 0.0;

/* go through all rows and columns */
for(row=start; row<n_rows; row++)
{
for(col=start; col<n_cols; col++)
{
double tmp = fabs(get_pivot_matrix_elt(A,row,col));
/* Compare element to current maximum */
if (tmp > mabsval)
{
mrow = row;
mcol = col;
mabsval = tmp;
}
}
}

/* assert that the matrix is not identical to zero */
if ((mrow < 0) || (mcol < 0)) return -1;

/* return result */
*maxrow = mrow;
*maxcol = mcol;
*maxabsval = mabsval;
return 0;
}



/*
pivot performs a full pivotization of a rectangular matrix A of dimension n_cols x n_rows
rowInd and colInd are vectors of length nrwos and n_cols respectively.
They hold the old (and new) pivoting information, such that
A_pivoted[i,j] = A[rowInd[i], colInd[j]]
*/
int pivot( double *A, int n_rows, int n_cols, int *rowInd, int *colInd )
{
/* parameter, determines how much larger an element should be before rows and columns are interchanged */
const double fac = 1.125; /* approved by dymola ;) */

/* temporary variables */
int row;
int i,j;
int maxrow;
int maxcol;
double maxabsval;
double pivot;

/* go over all pivot elements */
for(row=0; row<min(n_rows,n_cols); row++)
{
/* get current pivot */
pivot = fabs(get_pivot_matrix_elt(A,row,row));

/* find the maximum element in matrix
result is stored in maxrow, maxcol and maxabsval */
if (maxsearch(A, row, n_rows, n_cols, rowInd, colInd, &maxrow, &maxcol, &maxabsval) != 0) return -1;


/* compare max element and pivot (scaled by fac) */
if (maxabsval > (fac*pivot))
{
/* row interchange */
swap(rowInd[row], rowInd[maxrow]);
/* column interchange */
swap(colInd[row], colInd[maxcol]);
}

/* get pivot (without abs, may have changed because of row/column interchange */
pivot = get_pivot_matrix_elt(A,row,row);
/* internal error, pivot element should never be zero if maxsearch succeeded */
if(pivot != 0)
throw std::invalid_argument("pivot element is zero ");

/* perform one step of Gaussian Elimination */
for(i=row+1;i<n_rows;i++)
{
double leader = get_pivot_matrix_elt(A,i,row);
if (leader != 0.0)
{
double scale = -leader/pivot;
/* set leader to zero */
set_pivot_matrix_elt(A,i,row, 0.0);
/* subtract scaled equation from pivot row from current row */
for(j=row+1;j<n_cols;j++)
{
double t1 = get_pivot_matrix_elt(A,i,j);
double t2 = get_pivot_matrix_elt(A,row,j);
double tmp = t1 + scale*t2;
set_pivot_matrix_elt(A,i,j, tmp);
}
}
}
}
/* all fine */
return 0;
}
2 changes: 1 addition & 1 deletion SimulationRuntime/cpp/Include/Core/Math/Functions.h
View file
Expand Up @@ -236,7 +236,7 @@ inline bool in_range(T i,T start,T stop)




int pivot( double *A, int n_rows, int n_cols, int *rowInd, int *colInd );


// (C) Copyright Gennadiy Rozental 2001-2002.
Expand Down
7 changes: 3 additions & 4 deletions SimulationRuntime/cpp/Solver/CVode/CVode.cpp
View file
Expand Up @@ -177,14 +177,13 @@ void Cvode::initialize()
{
_idid = CVodeRootInit(_cvodeMem,_dimZeroFunc, CV_ZerofCallback);

for(int i=0;i<_dimZeroFunc;i++)
_zeroSign[i] = 1;


_idid = CVodeSetRootDirection(_cvodeMem, _zeroSign);
if(_idid < 0)
throw std::invalid_argument(/*_idid,_tCurrent,*/"CVode::initialize()");
memset(_zeroSign,0,_dimZeroFunc*sizeof(int));
memset(_zeroVal,0,_dimZeroFunc*sizeof(int));
memset(_zeroSign,-1,_dimZeroFunc*sizeof(int));
memset(_zeroVal,-1,_dimZeroFunc*sizeof(int));

}
_cvode_initialized = true;
Expand Down

0 comments on commit a6179af

Please sign in to comment.