Skip to content

Commit

Permalink
- linear solver "lis" added
Browse files Browse the repository at this point in the history
git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@16871 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
  • Loading branch information
Marcus Walther committed Aug 20, 2013
1 parent 8945a66 commit 1d12ed1
Show file tree
Hide file tree
Showing 8 changed files with 265 additions and 17 deletions.
10 changes: 6 additions & 4 deletions Compiler/Template/CodegenC.tpl
Expand Up @@ -948,9 +948,11 @@ template functionInitialLinearSystemsTemp(list<SimEqSystem> allEquations)
case eq as SES_MIXED(__) then functionInitialLinearSystemsTemp(fill(eq.cont,1))
case eq as SES_LINEAR(__) then
let size = listLength(vars)
let nnz = listLength(simJac)
<<
linearSystemData[<%indexLinearSystem%>].equationIndex = <%index%>;
linearSystemData[<%indexLinearSystem%>].size = <%size%>;
linearSystemData[<%indexLinearSystem%>].nnz = <%nnz%>;
linearSystemData[<%indexLinearSystem%>].setA = setLinearMatrixA<%index%>;
linearSystemData[<%indexLinearSystem%>].setb = setLinearVectorb<%index%>;
>>
Expand Down Expand Up @@ -987,10 +989,10 @@ template functionSetupLinearSystemsTemp(list<SimEqSystem> allEquations)
case eq as SES_MIXED(__) then functionSetupLinearSystemsTemp(fill(eq.cont,1))
case eq as SES_LINEAR(__) then
let &varDecls = buffer "" /*BUFD*/
let MatrixA = (simJac |> (row, col, eq as SES_RESIDUAL(__)) =>
let MatrixA = (simJac |> (row, col, eq as SES_RESIDUAL(__)) hasindex i0 =>
let &preExp = buffer "" /*BUFD*/
let expPart = daeExp(eq.exp, contextSimulationDiscrete, &preExp /*BUFC*/, &varDecls /*BUFD*/)
'<%preExp%>linearSystemData->A[<%row%> + <%col%> * linearSystemData->size] = <%expPart%>;'
'<%preExp%>linearSystemData->setAElement(<%row%>, <%col%>, <%expPart%>, <%i0%>, linearSystemData);'
;separator="\n")
let &varDecls2 = buffer "" /*BUFD*/
let vectorb = (beqs |> exp hasindex i0 =>
Expand Down Expand Up @@ -3231,7 +3233,7 @@ case SIMCODE(modelInfo=MODELINFO(__), makefileParams=MAKEFILE_PARAMS(__), simula
# /MD - link with MSVCRT.LIB
# /link - [linker options and libraries]
# /LIBPATH: - Directories where libs can be found
LDFLAGS=/MD /link /STACK:0x2000000 /pdb:"<%fileNamePrefix%>.pdb" /LIBPATH:"<%makefileParams.omhome%>/lib/omc/msvc/" /LIBPATH:"<%makefileParams.omhome%>/lib/omc/msvc/release/" <%dirExtra%> <%libsPos1%> <%libsPos2%> f2c.lib initialization.lib libexpat.lib math-support.lib meta.lib results.lib simulation.lib solver.lib sundials_kinsol.lib sundials_nvecserial.lib util.lib lapack_win32_MT.lib
LDFLAGS=/MD /link /STACK:0x2000000 /pdb:"<%fileNamePrefix%>.pdb" /LIBPATH:"<%makefileParams.omhome%>/lib/omc/msvc/" /LIBPATH:"<%makefileParams.omhome%>/lib/omc/msvc/release/" <%dirExtra%> <%libsPos1%> <%libsPos2%> f2c.lib initialization.lib libexpat.lib math-support.lib meta.lib results.lib simulation.lib solver.lib sundials_kinsol.lib sundials_nvecserial.lib util.lib lapack_win32_MT.lib lis.lib

# /MDd link with MSVCRTD.LIB debug lib
# lib names should not be appended with a d just switch to lib/omc/msvc/debug
Expand Down Expand Up @@ -3274,7 +3276,7 @@ case SIMCODE(modelInfo=MODELINFO(__), makefileParams=MAKEFILE_PARAMS(__), simula
CFLAGS=$(CFLAGS_BASED_ON_INIT_FILE) <%makefileParams.cflags%> <%match sopt case SOME(s as SIMULATION_SETTINGS(__)) then s.cflags /* From the simulate() command */%>
CPPFLAGS=-I"<%makefileParams.omhome%>/include/omc" -I. <%dirExtra%> <%makefileParams.includes ; separator=" "%> -DOPENMODELICA_XML_FROM_FILE_AT_RUNTIME
LIBSIMULATIONRUNTIMEC=<% if boolAnd(boolNot(stringEq(os(), "OSX")), boolOr(acceptMetaModelicaGrammar(), Flags.isSet(Flags.GEN_DEBUG_SYMBOLS))) then "-Wl,-whole-archive "%>-lSimulationRuntimeC <% if boolAnd(boolNot(stringEq(os(), "OSX")), boolOr(acceptMetaModelicaGrammar(), Flags.isSet(Flags.GEN_DEBUG_SYMBOLS))) then " -Wl,-no-whole-archive"%> <% if stringEq(makefileParams.platform, "win32") then "" else " -ldl"%>
LDFLAGS=-L"<%makefileParams.omhome%>/lib/omc" -Wl,<% if stringEq(makefileParams.platform, "win32") then "--stack,0x2000000,"%>-rpath,'<%makefileParams.omhome%>/lib/omc' $(LIBSIMULATIONRUNTIMEC) -linteractive <%ParModelicaLibs%> <%makefileParams.ldflags%> <%makefileParams.runtimelibs%> <%match System.os() case "OSX" then "-lf2c" else "-Wl,-Bstatic -lf2c -Wl,-Bdynamic"%>
LDFLAGS=-L"<%makefileParams.omhome%>/lib/omc" -Wl,<% if stringEq(makefileParams.platform, "win32") then "--stack,0x2000000,"%>-rpath,'<%makefileParams.omhome%>/lib/omc' $(LIBSIMULATIONRUNTIMEC) -linteractive <%ParModelicaLibs%> <%makefileParams.ldflags%> <%makefileParams.runtimelibs%> <%match System.os() case "OSX" then "-lf2c" else "-Wl,-Bstatic -lf2c -Wl,-Bdynamic -llis"%>
PERL=perl
FILEPREFIX=<%fileNamePrefix%>
MAINFILE=$(FILEPREFIX)<% if boolOr(acceptMetaModelicaGrammar(), Flags.isSet(Flags.GEN_DEBUG_SYMBOLS)) then ".conv"%>.c
Expand Down
4 changes: 4 additions & 0 deletions SimulationRuntime/c/simulation/simulation_runtime.cpp
Expand Up @@ -295,10 +295,14 @@ int getlinearSolverMethod(int argc, char**argv)
if(*method == string("lapack"))
return LS_LAPACK;

if(*method == string("lis"))
return LS_LIS;

WARNING1(LOG_STDOUT, "unrecognized option -ls %s", method->c_str());
WARNING(LOG_STDOUT, "current options are:");
INDENT(LOG_STDOUT);
WARNING2(LOG_STDOUT, "%-18s [%s]", "lapack", "default method");
WARNING2(LOG_STDOUT, "%-18s [%s]", "lis", "Lis");
THROW("see last warning");
return LS_NONE;
}
Expand Down
2 changes: 2 additions & 0 deletions SimulationRuntime/c/simulation/solver/CMakeLists.txt
Expand Up @@ -13,6 +13,7 @@ SET(solver_sources dassl.c
kinsolSolver.c
linearSystem.c
linearSolverLapack.c
linearSolverLis.c
nonlinearSystem.c
nonlinearSolverHybrd.c
nonlinearSolverNewton.c
Expand All @@ -28,6 +29,7 @@ SET(solver_headers dassl.h
kinsolSolver.h
linearSystem.h
linearSolverLapack.h
linearSolverLis.h
model_help.h
nonlinearSolverHybrd.h
nonlinearSolverNewton.h
Expand Down
142 changes: 142 additions & 0 deletions SimulationRuntime/c/simulation/solver/linearSolverLis.c
@@ -0,0 +1,142 @@
/*
* This file is part of OpenModelica.
*
* Copyright (c) 1998-CurrentYear, Linköping University,
* Department of Computer and Information Science,
* SE-58183 Linköping, Sweden.
*
* All rights reserved.
*
* THIS PROGRAM IS PROVIDED UNDER THE TERMS OF GPL VERSION 3
* AND THIS OSMC PUBLIC LICENSE (OSMC-PL).
* ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES RECIPIENT'S
* ACCEPTANCE OF THE OSMC PUBLIC LICENSE.
*
* The OpenModelica software and the Open Source Modelica
* Consortium (OSMC) Public License (OSMC-PL) are obtained
* from Linköping University, either from the above address,
* from the URLs: http://www.ida.liu.se/projects/OpenModelica or
* http://www.openmodelica.org, and in the OpenModelica distribution.
* GNU version 3 is obtained from: http://www.gnu.org/copyleft/gpl.html.
*
* This program is distributed WITHOUT ANY WARRANTY; without
* even the implied warranty of MERCHANTABILITY or FITNESS
* FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY SET FORTH
* IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS
* OF OSMC-PL.
*
* See the full OSMC Public License conditions for more details.
*
*/

/*! \file linearSolverLis.c
*/

#include <math.h>
#include <stdlib.h>
#include <string.h>

#include "simulation_data.h"
#include "omc_error.h"
#include "varinfo.h"
#include "model_help.h"

#include "linearSystem.h"
#include "linearSolverLis.h"

/*! \fn allocate memory for linear system solver Lis
*
*/
int
allocateLisData(int n_row, int n_col, int nz, void** voiddata)
{
DATA_LIS* data = (DATA_LIS*) malloc(sizeof(DATA_LIS));
ASSERT(data, "Could not allocate data for linear solver Lis.");

data->n_col = n_col;
data->n_row = n_row;
data->nnz = nz;

lis_vector_create(0, &(data->b));
lis_vector_set_size(data->b, data->n_row, 0);

lis_solver_create(&(data->solver));
lis_solver_set_option("-i fgmres ",data->solver);

*voiddata = (void*)data;
return 0;
}

/*! \fn free memory for linear system solver Lis
*
*/
int
freeLisData(void **voiddata)
{
DATA_LIS* data = (DATA_LIS*) *voiddata;

lis_solver_destroy(data->solver);
lis_vector_destroy(data->b);
lis_vector_destroy(data->x);

return 0;
}

/*! \fn solve linear system with Lis method
*
* \param [in] [data]
* [sysNumber] index of the corresponing non-linear system
*
* \author swagner
*/
int
solveLis(DATA *data, int sysNumber)
{
int i, ret, success = 1;
LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo.linearSystemData[sysNumber]);
DATA_LIS* sData = (DATA_LIS*)systemData->solverData;

/* Destroy the old matrix, create a new one and fill it with values */
lis_matrix_destroy(sData->A);
lis_matrix_create(0, &(sData->A));
lis_matrix_set_size(sData->A, sData->n_row, 0);
systemData->setA(data, systemData);
lis_matrix_set_type(sData->A, LIS_MATRIX_CSR);
lis_matrix_assemble(sData->A);

/* fill b with values */
systemData->setb(data, systemData);
for (i=0; i<sData->n_row; i++){
lis_vector_set_value(LIS_INS_VALUE, i, systemData->b[i], sData->b);
}

/* Create a new Vector for the solution */
lis_vector_destroy(sData->x);
lis_vector_duplicate(sData->A,&(sData->x));

/* solve */
ret = lis_solve(sData->A,sData->b,sData->x,sData->solver);

/* handle return status */
switch(ret){
case LIS_SUCCESS:
success = 1;
break;
case LIS_ILL_OPTION:
case LIS_BREAKDOWN:
case LIS_OUT_OF_MEMORY:
case LIS_MAXITER:
case LIS_ERR_NOT_IMPLEMENTED:
case LIS_ERR_FILE_IO:
default:
success = 0;
break;
}

/* write solution */
lis_vector_get_values(sData->x, 0, sData->n_col, systemData->x);

lis_matrix_destroy(sData->A);

return success;
}
70 changes: 70 additions & 0 deletions SimulationRuntime/c/simulation/solver/linearSolverLis.h
@@ -0,0 +1,70 @@
/*
* This file is part of OpenModelica.
*
* Copyright (c) 1998-CurrentYear, Linköping University,
* Department of Computer and Information Science,
* SE-58183 Linköping, Sweden.
*
* All rights reserved.
*
* THIS PROGRAM IS PROVIDED UNDER THE TERMS OF GPL VERSION 3
* AND THIS OSMC PUBLIC LICENSE (OSMC-PL).
* ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES RECIPIENT'S
* ACCEPTANCE OF THE OSMC PUBLIC LICENSE.
*
* The OpenModelica software and the Open Source Modelica
* Consortium (OSMC) Public License (OSMC-PL) are obtained
* from Linköping University, either from the above address,
* from the URLs: http://www.ida.liu.se/projects/OpenModelica or
* http://www.openmodelica.org, and in the OpenModelica distribution.
* GNU version 3 is obtained from: http://www.gnu.org/copyleft/gpl.html.
*
* This program is distributed WITHOUT ANY WARRANTY; without
* even the implied warranty of MERCHANTABILITY or FITNESS
* FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY SET FORTH
* IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS
* OF OSMC-PL.
*
* See the full OSMC Public License conditions for more details.
*
*/

/*! \file linearSolverLis.h
*/

#ifndef _LINEARSOLVERLIS_H_
#define _LINEARSOLVERLIS_H_

#include "simulation_data.h"
#include <lis.h>

#ifdef __cplusplus
extern "C" {
#endif

#ifdef VOID
#undef VOID
#endif

#ifdef __cplusplus
}
#endif

typedef struct DATA_LIS
{
LIS_MATRIX A;
LIS_VECTOR x,b;
LIS_SOLVER solver;
LIS_INT err;
int n_col;
int n_row;
int nnz;
} DATA_LIS;


int allocateLisData(int n_row, int n_col, int nz, void **data);
int freeLisData(void **data);
int solveLis(DATA *data, int sysNumber);

#endif

48 changes: 35 additions & 13 deletions SimulationRuntime/c/simulation/solver/linearSystem.c
Expand Up @@ -37,10 +37,12 @@
#include "omc_error.h"
#include "linearSystem.h"
#include "linearSolverLapack.h"
#include "linearSolverLis.h"
#include "simulation_info_xml.h"
#include "blaswrap.h"
#include "f2c.h"

int globalCounterIterations = 0;
/*! \fn int allocatelinearSystem(DATA *data)
*
* This function allocates memory for all linear systems.
Expand All @@ -49,25 +51,31 @@
*/
int allocatelinearSystem(DATA *data)
{
int i;
int i,nnz;
int size;
LINEAR_SYSTEM_DATA *linsys = data->simulationInfo.linearSystemData;

for(i=0; i<data->modelData.nLinearSystems; ++i)
{
size = linsys[i].size;
nnz = linsys[i].nnz;

/* allocate system data */
linsys[i].x = (double*) malloc(size*sizeof(double));
linsys[i].b = (double*) malloc(size*sizeof(double));
linsys[i].A = (double*) malloc(size*size*sizeof(double));

/* allocate solver data */
switch(data->simulationInfo.lsMethod)
{
/* the implementation of matrix A is solver-specific */
switch(data->simulationInfo.lsMethod){
case LS_LAPACK:
linsys[i].A = (double*) malloc(size*size*sizeof(double));
linsys[i].setAElement = setAElementLAPACK;
allocateLapackData(size, &linsys[i].solverData);
break;
case LS_LIS:
linsys[i].setAElement = setAElementLis;
allocateLisData(size, size, nnz, &linsys[i].solverData);
break;
default:
THROW("unrecognized linear solver");
}
Expand All @@ -88,15 +96,17 @@ int freelinearSystem(DATA *data)

for(i=0;i<data->modelData.nLinearSystems;++i)
{
/* free system and solver data */
free(linsys[i].x);
free(linsys[i].b);
free(linsys[i].A);

/* allocate solver data */
switch(data->simulationInfo.lsMethod)
{
switch(data->simulationInfo.lsMethod){
case LS_LAPACK:
freeLapackData(&linsys[i].solverData);
free(linsys[i].A);
break;
case LS_LIS:
freeLisData(&linsys[i].solverData);
break;
default:
THROW("unrecognized linear solver");
Expand All @@ -117,17 +127,16 @@ int freelinearSystem(DATA *data)
*/
int solve_linear_system(DATA *data, int sysNumber)
{
/* NONLINEAR_SYSTEM_DATA
* system = &(data->simulationInfo.nonlinearSystemData[sysNumber]); */
int success;
LINEAR_SYSTEM_DATA* linsys = data->simulationInfo.linearSystemData;

/* for now just use lapack solver as before */
switch(data->simulationInfo.lsMethod)
{
switch(data->simulationInfo.lsMethod){
case LS_LAPACK:
success = solveLapack(data, sysNumber);
break;
case LS_LIS:
success = solveLis(data, sysNumber);
break;
default:
THROW("unrecognized linear solver");
}
Expand Down Expand Up @@ -166,3 +175,16 @@ int check_linear_solutions(DATA *data, int printFailingSystems)

return retVal;
}

void setAElementLAPACK(int row, int col, double value, int nth, void *data )
{
LINEAR_SYSTEM_DATA* linsys = (LINEAR_SYSTEM_DATA*) data;
linsys->A[row + col * linsys->size] = value;
}

void setAElementLis(int row, int col, double value, int nth, void *data )
{
LINEAR_SYSTEM_DATA* linSys = (LINEAR_SYSTEM_DATA*) data;
DATA_LIS* sData = (DATA_LIS*) linSys->solverData;
lis_matrix_set_value(LIS_INS_VALUE, row, col, value, sData->A);
}

0 comments on commit 1d12ed1

Please sign in to comment.