Skip to content

Commit

Permalink
fix printSpasityPattern for daeMode
Browse files Browse the repository at this point in the history 
 - updated sparsity generation with leading 0
 - update all usage of this
  • Loading branch information
@OpenModelica-Hudson
Willi Braun authored and OpenModelica-Hudson committed Dec 7, 2016
1 parent 3f5a3c0 commit 0d56452
Show file tree
Hide file tree
Showing 18 changed files with 85 additions and 160 deletions.
44 changes: 15 additions & 29 deletions Compiler/Template/CodegenC.tpl
View file
Expand Up @@ -2692,9 +2692,9 @@ template generateStaticSparseData(String indexName, String systemType, list<tupl
inSysData->sparsePattern.maxColors = <%maxColor%>;
/* write lead index of compressed sparse column */
memcpy( inSysData->sparsePattern.leadindex, colPtrIndex, (<%sizeleadindex%>+1)*sizeof(int));
memcpy(inSysData->sparsePattern.leadindex, colPtrIndex, (<%sizeleadindex%>+1)*sizeof(int));
for(i=1;i<<%sizeleadindex%>+1;++i)
for(i=2;i<<%sizeleadindex%>+1;++i)
inSysData->sparsePattern.leadindex[i] += inSysData->sparsePattern.leadindex[i-1];
/* call sparse index */
Expand Down Expand Up @@ -4032,7 +4032,7 @@ template initializeDAEmodeData(Integer nResVars, Integer nAlgVars, list<tuple<In
/* write lead index of compressed sparse column */
memcpy(daeModeData->sparsePattern->leadindex, colPtrIndex, (1+<%sizeCols%>)*sizeof(int));
/* makek CRS compatible */
for(i=2;i<=<%sizeCols%>;++i)
for(i=2;i<<%sizeCols%>+1;++i)
daeModeData->sparsePattern->leadindex[i] += daeModeData->sparsePattern->leadindex[i-1];
/* call sparse index */
memcpy(daeModeData->sparsePattern->index, rowIndex, <%sizeNNZ%>*sizeof(int));
Expand Down Expand Up @@ -4585,24 +4585,10 @@ case _ then
let &eachCrefParts = buffer ""
let sp_size_index = lengthListElements(unzipSecond(sparsepattern))
let sizeleadindex = listLength(sparsepattern)
let leadindex = (sparsepattern |> (i, indexes) =>
<<
<%listLength(indexes)%>
>>
;separator=",")
let indexElems = ( sparsepattern |> (i, indexes) hasindex index0 =>
( indexes |> indexrow =>
<<
<%indexrow%>
>>
;separator=",")
;separator=",")
let colorArray = (colorList |> (indexes) hasindex index0 =>
let colorCol = ( indexes |> i_index =>
<<data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[<%i_index%>] = <%intAdd(index0,1)%>;>>
;separator="\n")
'<%colorCol%>'
;separator="\n")
let colPtr = genSPCRSPtr(listLength(sparsepattern), sparsepattern, "colPtrIndex")
let rowIndex = genSPCRSRows(lengthListElements(unzipSecond(sparsepattern)), sparsepattern, "rowIndex")
let colorString = genSPColors(colorList, "data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols")

let indexColumn = (jacobianColumn |> (eqs,vars,indxColumn) => indxColumn;separator="\n")
let tmpvarsSize = (jacobianColumn |> (_,vars,_) => listLength(vars);separator="\n")
let index_ = listLength(seedVars)
Expand All @@ -4613,8 +4599,8 @@ case _ then
TRACE_PUSH
DATA* data = ((DATA*)inData);
int index = <%symbolName(modelNamePrefix,"INDEX_JAC_")%><%matrixname%>;
const int tmp[<%sizeleadindex%>] = {<%leadindex%>};
const int tmpElem[<%sp_size_index%>] = {<%indexElems%>};
<%colPtr%>
<%rowIndex%>
int i = 0;
data->simulationInfo->analyticJacobians[index].sizeCols = <%index_%>;
Expand All @@ -4623,24 +4609,24 @@ case _ then
data->simulationInfo->analyticJacobians[index].seedVars = (modelica_real*) calloc(<%index_%>,sizeof(modelica_real));
data->simulationInfo->analyticJacobians[index].resultVars = (modelica_real*) calloc(<%indexColumn%>,sizeof(modelica_real));
data->simulationInfo->analyticJacobians[index].tmpVars = (modelica_real*) calloc(<%tmpvarsSize%>,sizeof(modelica_real));
data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex = (unsigned int*) malloc(<%sizeleadindex%>*sizeof(int));
data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex = (unsigned int*) malloc((<%sizeleadindex%>+1)*sizeof(int));
data->simulationInfo->analyticJacobians[index].sparsePattern.index = (unsigned int*) malloc(<%sp_size_index%>*sizeof(int));
data->simulationInfo->analyticJacobians[index].sparsePattern.numberOfNoneZeros = <%sp_size_index%>;
data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols = (unsigned int*) malloc(<%index_%>*sizeof(int));
data->simulationInfo->analyticJacobians[index].sparsePattern.maxColors = <%maxColor%>;
data->simulationInfo->analyticJacobians[index].jacobian = NULL;
/* write lead index of compressed sparse column */
memcpy(data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex, tmp, <%sizeleadindex%>*sizeof(int));
memcpy(data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex, colPtrIndex, (<%sizeleadindex%>+1)*sizeof(int));
for(i=1;i<<%sizeleadindex%>;++i)
data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[i] += data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[i-1];
for(i=2;i<<%sizeleadindex%>+1;++i)
data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[i] += data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[i-1];
/* call sparse index */
memcpy(data->simulationInfo->analyticJacobians[index].sparsePattern.index, tmpElem, <%sp_size_index%>*sizeof(int));
memcpy(data->simulationInfo->analyticJacobians[index].sparsePattern.index, rowIndex, <%sp_size_index%>*sizeof(int));
/* write color array */
<%colorArray%>
<%colorString%>
TRACE_POP
return 0;
}
Expand Down
4 changes: 0 additions & 4 deletions SimulationRuntime/c/optimization/DataManagement/DerStructure.c
View file
Expand Up @@ -230,11 +230,7 @@ static inline void local_jac_struct(DATA * data, OptDataDim * dim, OptDataStruct
lindex = (unsigned int*) data->simulationInfo->analyticJacobians[h_index].sparsePattern.leadindex;
pindex = data->simulationInfo->analyticJacobians[h_index].sparsePattern.index;

s->lindex[index] = (unsigned int*)calloc((sizeCols+1), sizeof(unsigned int));
memcpy(&s->lindex[index][1], lindex, sizeCols*sizeof(unsigned int));
lindex = s->lindex[index];
s->seedVec[index] = (modelica_real **)malloc((maxColors)*sizeof(modelica_real*));
free(data->simulationInfo->analyticJacobians[h_index].sparsePattern.leadindex);
/**********************/
if(sizeCols > 0){
for(ii = 1; ii < maxColors; ++ii){
Expand Down
4 changes: 2 additions & 2 deletions SimulationRuntime/c/optimization/DataManagement/MoveData.c
View file
Expand Up @@ -816,7 +816,7 @@ void diffSynColoredOptimizerSystem(OptData *optData, modelica_real **J, const in
const int h_index = optData->s.indexABCD[index];
const long double * scaldt = optData->bounds.scaldt[m];
const unsigned int * const cC = data->simulationInfo->analyticJacobians[h_index].sparsePattern.colorCols;
const unsigned int * const lindex = optData->s.lindex[index];
const unsigned int * const lindex = data->simulationInfo->analyticJacobians[h_index].sparsePattern.leadindex;
const int nx = data->simulationInfo->analyticJacobians[h_index].sizeCols;
const int Cmax = data->simulationInfo->analyticJacobians[h_index].sparsePattern.maxColors + 1;
const int dnx = optData->dim.nx;
Expand Down Expand Up @@ -869,7 +869,7 @@ void diffSynColoredOptimizerSystemF(OptData *optData, modelica_real **J){
const int index = 4;
const int h_index = optData->s.indexABCD[index];
const unsigned int * const cC = data->simulationInfo->analyticJacobians[h_index].sparsePattern.colorCols;
const unsigned int * const lindex = optData->s.lindex[index];
const unsigned int * const lindex = data->simulationInfo->analyticJacobians[h_index].sparsePattern.leadindex;
const int nx = data->simulationInfo->analyticJacobians[h_index].sizeCols;
const int Cmax = data->simulationInfo->analyticJacobians[h_index].sparsePattern.maxColors + 1;
const modelica_real * const resultVars = data->simulationInfo->analyticJacobians[h_index].resultVars;
Expand Down
14 changes: 4 additions & 10 deletions SimulationRuntime/c/simulation/solver/dassl.c
View file
Expand Up @@ -988,11 +988,8 @@ int functionJacAColored(DATA* data, threadData_t *threadData, double* jac)
{
if(data->simulationInfo->analyticJacobians[index].seedVars[j] == 1)
{
if(j==0)
ii = 0;
else
ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j-1];
while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j])
ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j];
while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j+1])
{
l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii];
k = j*data->simulationInfo->analyticJacobians[index].sizeRows + l;
Expand Down Expand Up @@ -1276,11 +1273,8 @@ int jacA_numColored(DATA* data, double *t, double *y, double *yprime, double *de
{
if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[ii]-1 == i)
{
if(ii==0)
j = 0;
else
j = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[ii-1];
while(j < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[ii])
j = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[ii];
while(j < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[ii+1])
{
l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[j];
k = l + ii*data->simulationInfo->analyticJacobians[index].sizeRows;
Expand Down
53 changes: 12 additions & 41 deletions SimulationRuntime/c/simulation/solver/ida_solver.c
View file
Expand Up @@ -1161,29 +1161,13 @@ int jacOwnNumColoredIDA(double tt, N_Vector yy, N_Vector yp, N_Vector rr, DlsMat
{
if(sparsePattern->colorCols[ii]-1 == i)
{
if (idaData->daeMode)
{
j = sparsePattern->leadindex[ii];
while(j < sparsePattern->leadindex[ii+1])
{
l = sparsePattern->index[j];
DENSE_ELEM(Jac, l, ii) = (newdelta[l] - delta[l]) * delta_hh[ii];
j++;
};
}
else
j = sparsePattern->leadindex[ii];
while(j < sparsePattern->leadindex[ii+1])
{
if(ii==0)
j = 0;
else
j = sparsePattern->leadindex[ii-1];
while(j < sparsePattern->leadindex[ii])
{
l = sparsePattern->index[j];
DENSE_ELEM(Jac, l, ii) = (newdelta[l] - delta[l]) * delta_hh[ii];
j++;
};
}
l = sparsePattern->index[j];
DENSE_ELEM(Jac, l, ii) = (newdelta[l] - delta[l]) * delta_hh[ii];
j++;
};
states[ii] = ysave[ii];
if (idaData->daeMode)
{
Expand Down Expand Up @@ -1448,26 +1432,13 @@ int jacobianSparseNumIDA(double tt, N_Vector yy, N_Vector yp, N_Vector rr, SlsMa
{
if(sparsePattern->colorCols[ii]-1 == i)
{
if (idaData->daeMode)
{
nth = sparsePattern->leadindex[ii];
while(nth < sparsePattern->leadindex[ii+1])
{
j = sparsePattern->index[nth];
setJacElementKluSparse(j, ii, (newdelta[j] - delta[j]) * delta_hh[ii], nth, Jac);
nth++;
};
}
else
nth = sparsePattern->leadindex[ii];
while(nth < sparsePattern->leadindex[ii+1])
{
nth = (ii == 0) ? 0 : sparsePattern->leadindex[ii-1];
while(nth < sparsePattern->leadindex[ii])
{
j = sparsePattern->index[nth];
setJacElementKluSparse(j, ii, (newdelta[j] - delta[j]) * delta_hh[ii], nth, Jac);
nth++;
};
}
j = sparsePattern->index[nth];
setJacElementKluSparse(j, ii, (newdelta[j] - delta[j]) * delta_hh[ii], nth, Jac);
nth++;
};
states[ii] = ysave[ii];
if (idaData->daeMode)
{
Expand Down
7 changes: 2 additions & 5 deletions SimulationRuntime/c/simulation/solver/linearSolverKlu.c
View file
Expand Up @@ -136,11 +136,8 @@ int getAnalyticalJacobian(DATA* data, threadData_t *threadData, int sysNumber)
{
if(data->simulationInfo->analyticJacobians[index].seedVars[j] == 1)
{
if(j==0)
ii = 0;
else
ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j-1];
while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j])
ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j];
while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j+1])
{
l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii];
systemData->setAElement(i, l, -data->simulationInfo->analyticJacobians[index].resultVars[l], nth, (void*) systemData, threadData);
Expand Down
7 changes: 2 additions & 5 deletions SimulationRuntime/c/simulation/solver/linearSolverLapack.c
View file
Expand Up @@ -119,11 +119,8 @@ int getAnalyticalJacobianLapack(DATA* data, threadData_t *threadData, double* ja
{
if(data->simulationInfo->analyticJacobians[index].seedVars[j] == 1)
{
if(j==0)
ii = 0;
else
ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j-1];
while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j])
ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j];
while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j+1])
{
l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii];
k = j*data->simulationInfo->analyticJacobians[index].sizeRows + l;
Expand Down
7 changes: 2 additions & 5 deletions SimulationRuntime/c/simulation/solver/linearSolverLis.c
View file
Expand Up @@ -156,11 +156,8 @@ int getAnalyticalJacobianLis(DATA* data, threadData_t *threadData, int sysNumber
{
if(data->simulationInfo->analyticJacobians[index].seedVars[j] == 1)
{
if(j==0)
ii = 0;
else
ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j-1];
while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j])
ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j-1];
while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j+1])
{
l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii];
/*infoStreamPrint(LOG_LS_V, 0, "set on Matrix A (%d, %d)(%d) = %f", i, l, nth, -data->simulationInfo->analyticJacobians[index].resultVars[l]); */
Expand Down
8 changes: 2 additions & 6 deletions SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.c
View file
Expand Up @@ -337,12 +337,8 @@ int getAnalyticalJacobianTotalPivot(DATA* data, threadData_t *threadData, double
{
if(data->simulationInfo->analyticJacobians[index].seedVars[j] == 1)
{
if(j==0) {
ii = 0;
} else {
ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j-1];
}
while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j]) {
ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j];
while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j+1]) {
l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii];
k = j*data->simulationInfo->analyticJacobians[index].sizeRows + l;
jac[k] = data->simulationInfo->analyticJacobians[index].resultVars[l];
Expand Down
7 changes: 2 additions & 5 deletions SimulationRuntime/c/simulation/solver/linearSolverUmfpack.c
View file
Expand Up @@ -150,11 +150,8 @@ int getAnalyticalJacobianUmfPack(DATA* data, threadData_t *threadData, int sysNu
{
if(data->simulationInfo->analyticJacobians[index].seedVars[j] == 1)
{
if(j==0)
ii = 0;
else
ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j-1];
while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j])
ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j];
while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j+1])
{
l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii];
/* infoStreamPrint(LOG_LS_V, 0, "set on Matrix A (%d, %d)(%d) = %f", i, l, nth, -data->simulationInfo->analyticJacobians[index].resultVars[l]); */
Expand Down
34 changes: 12 additions & 22 deletions SimulationRuntime/c/simulation/solver/model_help.c
View file
Expand Up @@ -350,44 +350,34 @@ void printParameters(DATA *data, int stream)
*
* prints sparse structure of jacobian A
*
* \param [in] [data]
* \param [in] [sparsePattern]
* \param [in] [sizeRow]
* \param [in] [sizeCol]
* \param [in] [stream]
*
* \author lochel
*/
void printSparseStructure(DATA *data, int stream)
void printSparseStructure(SPARSE_PATTERN *sparsePattern, int sizeRows, int sizeCols, int stream, const char* name)
{
const int index = data->callback->INDEX_JAC_A;
unsigned int row, col, i, j;
/* Will crash with a static size array */
char *buffer = NULL;

if (!ACTIVE_STREAM(stream))
return;

buffer = (char*)omc_alloc_interface.malloc(sizeof(char)* 2*data->simulationInfo->analyticJacobians[index].sizeCols + 4);

infoStreamPrint(stream, 1, "sparse structure of jacobian A [size: %ux%u]", data->simulationInfo->analyticJacobians[index].sizeRows, data->simulationInfo->analyticJacobians[index].sizeCols);
infoStreamPrint(stream, 0, "%u nonzero elements", data->simulationInfo->analyticJacobians[index].sparsePattern.numberOfNoneZeros);
/*
sprintf(buffer, "");
for(row=0; row < data->simulationInfo->analyticJacobians[index].sizeRows; row++)
sprintf(buffer, "%s%u ", buffer, data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[row]);
infoStreamPrint(stream, 0, "leadindex: %s", buffer);
sprintf(buffer, "");
for(i=0; i < data->simulationInfo->analyticJacobians[index].sparsePattern.numberOfNoneZeros; i++)
sprintf(buffer, "%s%u ", buffer, data->simulationInfo->analyticJacobians[index].sparsePattern.index[i]);
infoStreamPrint(stream, 0, "index: %s", buffer);
*/
buffer = (char*)omc_alloc_interface.malloc(sizeof(char)* 2*sizeCols + 4);

infoStreamPrint(stream, 1, "sparse structure of %s [size: %ux%u]", name, sizeRows, sizeCols);
infoStreamPrint(stream, 0, "%u nonzero elements", sparsePattern->numberOfNoneZeros);

infoStreamPrint(stream, 1, "transposed sparse structure (rows: states)");
i=0;
for(row=0; row < data->simulationInfo->analyticJacobians[index].sizeRows; row++)
for(row=0; row < sizeRows; row++)
{
j=0;
for(col=0; i < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[row]; col++)
for(col=0; i < sparsePattern->leadindex[row]+1; col++)
{
if(data->simulationInfo->analyticJacobians[index].sparsePattern.index[i] == col)
if(sparsePattern->index[i] == col)
{
buffer[j++] = '*';
++i;
Expand Down
2 changes: 1 addition & 1 deletion SimulationRuntime/c/simulation/solver/model_help.h
View file
Expand Up @@ -107,7 +107,7 @@ void printAllVars(DATA *data, int ringSegment, int stream);
void printRelations(DATA *data, int stream);
void printZeroCrossings(DATA *data, int stream);
void printParameters(DATA *data, int stream);
void printSparseStructure(DATA *data, int stream);
void printSparseStructure(SPARSE_PATTERN *sparsePattern, int sizeRows, int sizeCols, int stream, const char*);

void overwriteOldSimulationData(DATA *data);
void copyRingBufferSimulationData(DATA *data, threadData_t *threadData, SIMULATION_DATA **destData, RINGBUFFER* destRing);
Expand Down
7 changes: 2 additions & 5 deletions SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c
View file
Expand Up @@ -778,11 +778,8 @@ int getAnalyticalJacobianHomotopy(DATA_HOMOTOPY* solverData, double* jac)
{
if(data->simulationInfo->analyticJacobians[index].seedVars[j] == 1)
{
if(j==0)
ii = 0;
else
ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j-1];
while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j])
ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j];
while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j+1])
{
l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii];
k = j*data->simulationInfo->analyticJacobians[index].sizeRows + l;
Expand Down

0 comments on commit 0d56452

Please sign in to comment.