Skip to content

Commit ea1e2c1

Browse files
Andreas Heuermannlochel
Andreas Heuermann
authored and
lochel
committed
Linear solvers: Check residuals and fixes for Lis
 - Fixed wrong index in analytic Jacobian lead index for Lis solver - Switched column and rows in getAnalytical Jacobian for Lis - Added check for all linear solver using jacobian matrices for residual vector equals null - Added missing linear solver for linear solver tests
1 parent 4dc44f9 commit ea1e2c1

File tree

10 files changed

+189
-28
lines changed

10 files changed

+189
-28
lines changed

‎OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverKlu.c‎

Lines changed: 15 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,7 @@
4141
#include "simulation_data.h"
4242
#include "simulation/simulation_info_json.h"
4343
#include "util/omc_error.h"
44+
#include "omc_math.h"
4445
#include "util/varinfo.h"
4546
#include "model_help.h"
4647

@@ -182,6 +183,7 @@ solveKlu(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
182183
void *dataAndThreadData[2] = {data, threadData};
183184
LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]);
184185
DATA_KLU* solverData = (DATA_KLU*)systemData->solverData[0];
186+
_omc_scalar residualNorm = 0;
185187

186188
int i, j, status = 0, success = 0, n = systemData->size, eqSystemNumber = systemData->equationIndex, indexes[2] = {1,eqSystemNumber};
187189
double tmpJacEvalTime;
@@ -299,14 +301,26 @@ solveKlu(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
299301

300302
/* update inner equations */
301303
residual_wrapper(aux_x, solverData->work, dataAndThreadData, sysNumber);
304+
residualNorm = _omc_gen_euclideanVectorNorm(solverData->work, solverData->n_row);
305+
306+
if ((isnan(residualNorm)) || (residualNorm>1e-4)){
307+
warningStreamPrint(LOG_LS, 0,
308+
"Failed to solve linear system of equations (no. %d) at time %f. Residual norm is %.15g.",
309+
(int)systemData->equationIndex, data->localData[0]->timeValue, residualNorm);
310+
success = 0;
311+
}
302312
} else {
303313
/* the solution is automatically in x */
304314
memcpy(aux_x, systemData->b, sizeof(double)*systemData->size);
305315
}
306316

307317
if (ACTIVE_STREAM(LOG_LS_V))
308318
{
309-
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
319+
if (1 == systemData->method) {
320+
infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
321+
} else {
322+
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
323+
}
310324
infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);
311325

312326
for(i = 0; i < systemData->size; ++i)

‎OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLapack.c‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -304,7 +304,11 @@ int solveLapack(DATA *data, threadData_t *threadData, int sysNumber, double* aux
304304
}
305305

306306
if (ACTIVE_STREAM(LOG_LS_V)){
307-
infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
307+
if (1 == systemData->method) {
308+
infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
309+
} else {
310+
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
311+
}
308312
infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);
309313

310314
for(i = 0; i < systemData->size; ++i) {

‎OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLis.c‎

Lines changed: 36 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,7 @@
3838
#include "simulation_data.h"
3939
#include "simulation/simulation_info_json.h"
4040
#include "util/omc_error.h"
41+
#include "omc_math.h"
4142
#include "util/varinfo.h"
4243
#include "model_help.h"
4344

@@ -91,8 +92,7 @@ allocateLisData(int n_row, int n_col, int nz, void** voiddata)
9192
/*! \fn free memory for linear system solver Lis
9293
*
9394
*/
94-
int
95-
freeLisData(void **voiddata)
95+
int freeLisData(void **voiddata)
9696
{
9797
DATA_LIS* data = (DATA_LIS*) *voiddata;
9898

@@ -107,18 +107,26 @@ freeLisData(void **voiddata)
107107
}
108108

109109

110+
/*!
111+
* Print LIS_MATRIX provided in column sparse row format
112+
*
113+
* \param [in] A Matrix A of linear problem A*x=b in CSR format
114+
* \param [in] n Dimension of A
115+
*/
110116
void printLisMatrixCSR(LIS_MATRIX A, int n)
111117
{
112118
int i, j;
113119
/* A matrix */
114-
infoStreamPrint(LOG_LS_V, 1, "A matrix [%dx%d] nnz = %d ", n, n, A->nnz);
120+
infoStreamPrint(LOG_LS_V, 1, "A matrix [%dx%d] nnz = %d", n, n, A->nnz);
121+
infoStreamPrint(LOG_LS_V, 0, "Column Sparse Row format. Print tuple (index,value) for each row:");
115122
for(i=0; i<n; i++)
116123
{
117124
char *buffer = (char*)malloc(sizeof(char)*A->ptr[i+1]*50);
118125
buffer[0] = 0;
119-
for(j=A->ptr[i]; j<A->ptr[i+1]; j++)
126+
sprintf(buffer, "column %d: ", i);
127+
for(j = A->ptr[i]; j < A->ptr[i+1]; j++)
120128
{
121-
sprintf(buffer, "%s(%d,%d,%g) ", buffer, i, A->index[j], A->value[j]);
129+
sprintf(buffer, "%s(%d,%g) ", buffer, A->index[j], A->value[j]);
122130
}
123131
infoStreamPrint(LOG_LS_V, 0, "%s", buffer);
124132
free(buffer);
@@ -158,12 +166,12 @@ int getAnalyticalJacobianLis(DATA* data, threadData_t *threadData, int sysNumber
158166
{
159167
if(jacobian->seedVars[j] == 1)
160168
{
161-
ii = jacobian->sparsePattern.leadindex[j-1];
169+
ii = jacobian->sparsePattern.leadindex[j];
162170
while(ii < jacobian->sparsePattern.leadindex[j+1])
163171
{
164172
l = jacobian->sparsePattern.index[ii];
165-
/*infoStreamPrint(LOG_LS_V, 0, "set on Matrix A (%d, %d)(%d) = %f", i, l, nth, -jacobian->resultVars[l]); */
166-
systemData->setAElement(i, l, -jacobian->resultVars[l], nth, (void*) systemData, threadData);
173+
/* infoStreamPrint(LOG_LS_V, 0, "set on Matrix A (%d, %d)(%d) = %f", l, i, nth, -jacobian->resultVars[l]); */
174+
systemData->setAElement(l, i, -jacobian->resultVars[l], nth, (void*) systemData, threadData);
167175
nth++;
168176
ii++;
169177
};
@@ -175,7 +183,7 @@ int getAnalyticalJacobianLis(DATA* data, threadData_t *threadData, int sysNumber
175183
return 0;
176184
}
177185

178-
/*! \fn wrapper_fvec_umfpack for the residual function
186+
/*! \fn wrapper_fvec_lis for the residual function
179187
*
180188
*/
181189
static int wrapper_fvec_lis(double* x, double* f, void** data, int sysNumber)
@@ -193,16 +201,15 @@ static int wrapper_fvec_lis(double* x, double* f, void** data, int sysNumber)
193201
* [sysNumber] index of the corresponding linear system
194202
*
195203
*/
196-
197-
int
198-
solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
204+
int solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
199205
{
200206
void *dataAndThreadData[2] = {data, threadData};
201207
LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]);
202208
DATA_LIS* solverData = (DATA_LIS*)systemData->solverData[0];
203209
int i, ret, success = 1, ni, iflag = 1, n = systemData->size, eqSystemNumber = systemData->equationIndex;
204210
char *lis_returncode[] = {"LIS_SUCCESS", "LIS_ILL_OPTION", "LIS_BREAKDOWN", "LIS_OUT_OF_MEMORY", "LIS_MAXITER", "LIS_NOT_IMPLEMENTED", "LIS_ERR_FILE_IO"};
205211
LIS_INT err;
212+
_omc_scalar residualNorm = 0;
206213

207214
int indexes[2] = {1,eqSystemNumber};
208215
double tmpJacEvalTime;
@@ -216,10 +223,10 @@ solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
216223
}
217224

218225
rt_ext_tp_tick(&(solverData->timeClock));
226+
227+
lis_matrix_set_size(solverData->A, solverData->n_row, 0);
219228
if (0 == systemData->method)
220229
{
221-
222-
lis_matrix_set_size(solverData->A, solverData->n_row, 0);
223230
/* set A matrix */
224231
systemData->setA(data, threadData, systemData);
225232
lis_matrix_assemble(solverData->A);
@@ -228,8 +235,6 @@ solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
228235
systemData->setb(data, threadData, systemData);
229236

230237
} else {
231-
232-
lis_matrix_set_size(solverData->A, solverData->n_row, 0);
233238
/* calculate jacobian -> matrix A*/
234239
if(systemData->jacobianIndex != -1){
235240
getAnalyticalJacobianLis(data, threadData, sysNumber);
@@ -280,25 +285,37 @@ solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
280285
free(buffer);
281286
}
282287

283-
/* print solution */
288+
/* Log solution */
284289
if (1 == success){
285290

286-
if (1 == systemData->method){
291+
if (1 == systemData->method){ /* Case calculate jacobian -> matrix A*/
287292
/* take the solution */
288293
lis_vector_get_values(solverData->x, 0, solverData->n_row, aux_x);
289294
for(i = 0; i < solverData->n_row; ++i)
290295
aux_x[i] += solverData->work[i];
291296

292297
/* update inner equations */
293298
wrapper_fvec_lis(aux_x, solverData->work, dataAndThreadData, sysNumber);
299+
residualNorm = _omc_gen_euclideanVectorNorm(solverData->work, solverData->n_row);
300+
301+
if ((isnan(residualNorm)) || (residualNorm>1e-4)){
302+
warningStreamPrint(LOG_LS, 0,
303+
"Failed to solve linear system of equations (no. %d) at time %f. Residual norm is %.15g.",
304+
(int)systemData->equationIndex, data->localData[0]->timeValue, residualNorm);
305+
success = 0;
306+
}
294307
} else {
295308
/* write solution */
296309
lis_vector_get_values(solverData->x, 0, solverData->n_row, aux_x);
297310
}
298311

299312
if (ACTIVE_STREAM(LOG_LS_V))
300313
{
301-
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
314+
if (1 == systemData->method) {
315+
infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
316+
} else {
317+
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
318+
}
302319
infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);
303320

304321
for(i = 0; i < systemData->size; ++i)

‎OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLis.h‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -57,6 +57,5 @@ typedef struct DATA_LIS
5757
int allocateLisData(int n_row, int n_col, int nz, void **data);
5858
int freeLisData(void **data);
5959
int solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x);
60-
6160
#endif
6261

‎OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.c‎

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,7 @@
3838
#include "../../simulation_data.h"
3939
#include "../simulation_info_json.h"
4040
#include "../../util/omc_error.h"
41+
#include "omc_math.h"
4142
#include "../../util/varinfo.h"
4243
#include "model_help.h"
4344

@@ -320,8 +321,8 @@ int freeTotalPivotData(void** voiddata)
320321
*/
321322
int getAnalyticalJacobianTotalPivot(DATA* data, threadData_t *threadData, double* jac, int sysNumber)
322323
{
323-
int i,j,k,l,ii,currentSys = sysNumber;
324-
LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[currentSys]);
324+
int i,j,k,l,ii;
325+
LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[sysNumber]);
325326

326327
const int index = systemData->jacobianIndex;
327328
ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]);
@@ -394,6 +395,7 @@ int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber, double*
394395
int eqSystemNumber = systemData->equationIndex;
395396
int indexes[2] = {1,eqSystemNumber};
396397
int rank;
398+
_omc_scalar residualNorm = 0;
397399

398400
/* We are given the number of the linear system.
399401
* We want to look it up among all equations. */
@@ -459,7 +461,11 @@ int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber, double*
459461
}
460462

461463
if (ACTIVE_STREAM(LOG_LS_V)){
462-
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
464+
if (1 == systemData->method) {
465+
infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
466+
} else {
467+
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
468+
}
463469
infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);
464470
for(i=0; i<systemData->size; ++i)
465471
{

‎OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverUmfpack.c‎

Lines changed: 15 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,7 @@
4141
#include "simulation_data.h"
4242
#include "simulation/simulation_info_json.h"
4343
#include "util/omc_error.h"
44+
#include "omc_math.h"
4445
#include "util/varinfo.h"
4546
#include "model_help.h"
4647

@@ -196,6 +197,7 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
196197
void *dataAndThreadData[2] = {data, threadData};
197198
LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]);
198199
DATA_UMFPACK* solverData = (DATA_UMFPACK*)systemData->solverData[0];
200+
_omc_scalar residualNorm = 0;
199201

200202
int i, j, status = UMFPACK_OK, success = 0, ni=0, n = systemData->size, eqSystemNumber = systemData->equationIndex, indexes[2] = {1,eqSystemNumber};
201203
int casualTearingSet = systemData->strictTearingFunctionCall != NULL;
@@ -305,13 +307,25 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
305307

306308
/* update inner equations */
307309
wrapper_fvec_umfpack(aux_x, solverData->work, dataAndThreadData, sysNumber);
310+
residualNorm = _omc_gen_euclideanVectorNorm(solverData->work, solverData->n_row);
311+
312+
if ((isnan(residualNorm)) || (residualNorm>1e-4)){
313+
warningStreamPrint(LOG_LS, 0,
314+
"Failed to solve linear system of equations (no. %d) at time %f. Residual norm is %.15g.",
315+
(int)systemData->equationIndex, data->localData[0]->timeValue, residualNorm);
316+
success = 0;
317+
}
308318
} else {
309319
/* the solution is automatically in x */
310320
}
311321

312322
if (ACTIVE_STREAM(LOG_LS_V))
313323
{
314-
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
324+
if (1 == systemData->method) {
325+
infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
326+
} else {
327+
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
328+
}
315329
infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);
316330

317331
for(i = 0; i < systemData->size; ++i)

‎OMCompiler/SimulationRuntime/c/simulation/solver/omc_math.c‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -804,6 +804,25 @@ _omc_scalar _omc_euclideanVectorNorm(const _omc_vector* vec)
804804
return sqrt(result);
805805
}
806806

807+
/*! \fn _omc_scalar _omc_gen_euclideanVectorNorm(_omc_vector* vec)
808+
*
809+
* calculates the euclidean vector norm
810+
*
811+
* \param [in] [vec] !TODO: DESCRIBE ME!
812+
*/
813+
_omc_scalar _omc_gen_euclideanVectorNorm(const _omc_scalar* vec_data, const _omc_size vec_size)
814+
{
815+
_omc_size i;
816+
_omc_scalar result = 0;
817+
assertStreamPrint(NULL, vec_size > 0, "Vector size is greater than zero");
818+
assertStreamPrint(NULL, NULL != vec_data, "Vector data is NULL pointer");
819+
for (i = 0; i < vec_size; ++i) {
820+
result += pow(fabs(vec_data[i]),2.0);
821+
}
822+
823+
return sqrt(result);
824+
}
825+
807826
/*! \fn _omc_scalar _omc_maximumVectorNorm(_omc_vector* vec)
808827
*
809828
* calculates the maximum vector norm

‎OMCompiler/SimulationRuntime/c/simulation/solver/omc_math.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -119,6 +119,7 @@ void _omc_printMatrix(_omc_matrix* mat, const char* name, const int logLevel);
119119

120120
/* norm functions */
121121
_omc_scalar _omc_euclideanVectorNorm(const _omc_vector* vec);
122+
_omc_scalar _omc_gen_euclideanVectorNorm(const _omc_scalar* vec_data, const _omc_size vec_size);
122123
_omc_scalar _omc_maximumVectorNorm(const _omc_vector* vec);
123124

124125
#endif

0 commit comments

Comments
 (0)