diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 3fcb600e56..8fd85ae7be 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -28,13 +28,12 @@ jobs:
strategy:
matrix:
build-option: [ Debug, Release ]
- sparse-solver: [ EIGEN, MKL, MKL_RUNTIME ]
steps:
- uses: actions/checkout@v3
- name: Run build script
- run: ./build.sh ${{ matrix.build-option }} ${{ matrix.sparse-solver }}
+ run: ./build.sh ${{ matrix.build-option }}
build-cpp-test-windows:
if: (github.event_name == 'push') || (github.event_name == 'workflow_dispatch') || (!startsWith(github.head_ref, 'release'))
@@ -42,10 +41,6 @@ jobs:
strategy:
matrix:
build-option: [ Debug, Release ]
- sparse-solver: [ EIGEN, MKL, MKL_RUNTIME ]
- env:
- MKL_INCLUDE: C:\conda_envs\cpp_pkgs\Library\include
- MKL_LIB: C:\conda_envs\cpp_pkgs\Library\lib
steps:
- uses: actions/checkout@v3
@@ -57,7 +52,7 @@ jobs:
- name: Install conda environment
# TODO removed pinned v2 of catch2, see https://github.com/alliander-opensource/power-grid-model/issues/73
run: |
- conda create --yes -p C:\conda_envs\cpp_pkgs -c conda-forge boost-cpp eigen nlohmann_json mkl mkl-devel mkl-include catch2==2.13.9
+ conda create --yes -p C:\conda_envs\cpp_pkgs -c conda-forge boost-cpp eigen nlohmann_json catch2==2.13.9
- name: Build and test
run: |
@@ -65,17 +60,16 @@ jobs:
Import-Module (Join-Path $vsPath 'Common7\Tools\Microsoft.VisualStudio.DevShell.dll')
Enter-VsDevShell -VsInstallPath $vsPath -SkipAutomaticLocation -DevCmdArguments '-arch=x64 -host_arch=x64'
$env:Path += ";C:\conda_envs\cpp_pkgs\Library\bin"
- mkdir cpp_build_${{ matrix.build-option }}_${{ matrix.sparse-solver }}
- cd cpp_build_${{ matrix.build-option }}_${{ matrix.sparse-solver }}
+ mkdir cpp_build_${{ matrix.build-option }}
+ cd cpp_build_${{ matrix.build-option }}
# generate cmake cache
cmake .. `
-G "Ninja" `
-DCMAKE_BUILD_TYPE=${{ matrix.build-option }} `
- -DPOWER_GRID_MODEL_SPARSE_SOLVER=${{ matrix.sparse-solver }} `
-DCMAKE_PREFIX_PATH=C:\conda_envs\cpp_pkgs\Library `
-DPOWER_GRID_MODEL_BUILD_BENCHMARK=1
# build
- cmake --build . --verbose
+ cmake --build . --verbose -j 1
# test
.\tests\cpp_unit_tests\power_grid_model_unit_tests.exe
@@ -85,12 +79,9 @@ jobs:
strategy:
matrix:
build-option: [ Debug, Release ]
- sparse-solver: [ EIGEN, MKL, MKL_RUNTIME ]
env:
CC: clang
CXX: clang++
- MKL_INCLUDE: /usr/local/include
- MKL_LIB: /usr/local/lib
CMAKE_PREFIX_PATH: /usr/local
steps:
@@ -101,14 +92,10 @@ jobs:
run: |
brew install ninja boost eigen nlohmann-json
curl https://raw.githubusercontent.com/Homebrew/homebrew-core/5e5abb11bf49787d01164c4066119365262c21ed/Formula/catch2.rb > $(find $(brew --repository) -name catch2.rb) && brew reinstall catch2
- sudo pip3 install mkl mkl-devel mkl-include
- name: Build and test
run: |
- # environment
- export LD_LIBRARY_PATH=${MKL_LIB}:${LD_LIBRARY_PATH}
- # bash
- ./build.sh ${{ matrix.build-option }} ${{ matrix.sparse-solver }}
+ ./build.sh ${{ matrix.build-option }}
build-and-test-python:
strategy:
@@ -122,9 +109,6 @@ jobs:
boost:
eigen:
cibw_build: "cp*-manylinux_x86_64 cp*-manylinux_aarch64"
- test_cmd: >
- pytest {package}/tests &&
- LD_LIBRARY_PATH= pytest {package}/tests
archs: "x86_64 aarch64"
- platform: macos
os: macos-latest
@@ -134,9 +118,6 @@ jobs:
eigen: /usr/local/include/eigen3
cibw_build: cp*-macosx_*
archs: "x86_64 arm64"
- test_cmd: >
- pytest {package}/tests &&
- LD_LIBRARY_PATH=/usr/local/lib:${LD_LIBRARY_PATH} pytest {package}/tests
- platform: windows
os: windows-latest
cc:
@@ -145,7 +126,6 @@ jobs:
eigen: C:\conda_envs\cpp_pkgs\Library\include\eigen3
cibw_build: cp*-win_amd64
archs: AMD64
- test_cmd: pytest {package}/tests
runs-on: ${{ matrix.os }}
needs: [build-cpp-test-linux, build-cpp-test-windows, build-cpp-test-macos]
@@ -160,7 +140,7 @@ jobs:
CIBW_BUILD: ${{ matrix.cibw_build }}
CIBW_ARCHS: ${{ matrix.archs }}
CIBW_TEST_REQUIRES: pytest pytest-cov
- CIBW_TEST_COMMAND: ${{ matrix.test_cmd }}
+ CIBW_TEST_COMMAND: pytest {package}/tests
# Skip trying to test arm64 builds on Intel Macs
CIBW_TEST_SKIP: "*-macosx_arm64 *-macosx_universal2:arm64"
CIBW_MANYLINUX_X86_64_IMAGE: ghcr.io/alliander-opensource/power-grid-model-build-env
@@ -169,11 +149,10 @@ jobs:
CIBW_BEFORE_ALL_MACOS: >
which clang++ &&
clang++ --version &&
- brew install boost eigen &&
- sudo pip3 install mkl mkl-devel mkl-include
+ brew install boost eigen
MACOSX_DEPLOYMENT_TARGET: 10.15
CIBW_BEFORE_ALL_WINDOWS : >
- conda create --yes -p C:\conda_envs\cpp_pkgs -c conda-forge boost-cpp eigen mkl
+ conda create --yes -p C:\conda_envs\cpp_pkgs -c conda-forge boost-cpp eigen
steps:
- uses: actions/checkout@v3
diff --git a/.github/workflows/sonar.yml b/.github/workflows/sonar.yml
index 62b34c6267..a8f8a03190 100644
--- a/.github/workflows/sonar.yml
+++ b/.github/workflows/sonar.yml
@@ -61,7 +61,7 @@ jobs:
- name: Run build-wrapper for C++
run: |
mkdir cpp_build
- cmake -S . -B cpp_build -GNinja -DCMAKE_BUILD_TYPE=Debug -DPOWER_GRID_MODEL_SPARSE_SOLVER=EIGEN -DCMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH}
+ cmake -S . -B cpp_build -GNinja -DCMAKE_BUILD_TYPE=Debug -DCMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH}
VERBOSE=1 build-wrapper-linux-x86-64 --out-dir ${{ env.BUILD_WRAPPER_OUT_DIR }} cmake --build cpp_build/
- name: C++ test and coverage
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7703b11ae5..7bf396cc21 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -31,52 +31,23 @@ else()
add_compile_options(-Wall -Wextra -pedantic -Werror)
add_compile_options(-Wno-deprecated-copy) # bug in boost
if(UNIX AND NOT APPLE)
- list(APPEND CMAKE_FIND_LIBRARY_SUFFIXES .so.1)
- list(APPEND CMAKE_FIND_LIBRARY_SUFFIXES .so.2)
- # test coverage for debug build in linux
- if (CMAKE_BUILD_TYPE STREQUAL "Debug")
+ # test coverage for debug build in linux, if specified
+ if ((CMAKE_BUILD_TYPE STREQUAL "Debug") AND (DEFINED POWER_GRID_MODEL_COVERAGE) AND (POWER_GRID_MODEL_COVERAGE EQUAL 1))
add_compile_options(-fprofile-arcs -ftest-coverage)
add_link_options(-fprofile-arcs)
endif()
else()
- list(APPEND CMAKE_FIND_LIBRARY_SUFFIXES .1.dylib)
- list(APPEND CMAKE_FIND_LIBRARY_SUFFIXES .2.dylib)
endif()
# thread
find_package(Threads REQUIRED)
target_link_libraries(power_grid_model INTERFACE Threads::Threads)
endif()
-
-# sparse solver, default eigen
-if (NOT DEFINED POWER_GRID_MODEL_SPARSE_SOLVER)
- set(POWER_GRID_MODEL_SPARSE_SOLVER "EIGEN")
-endif()
-
-# use mkl solver
-if((POWER_GRID_MODEL_SPARSE_SOLVER STREQUAL "MKL") OR (POWER_GRID_MODEL_SPARSE_SOLVER STREQUAL "MKL_RUNTIME"))
- # check load mkl at runtime or link-time
- if(POWER_GRID_MODEL_SPARSE_SOLVER STREQUAL "MKL_RUNTIME")
- target_compile_definitions(power_grid_model INTERFACE POWER_GRID_MODEL_USE_MKL_AT_RUNTIME=1)
- target_link_libraries(power_grid_model INTERFACE ${CMAKE_DL_LIBS})
- else()
- if(NOT DEFINED ENV{MKL_INCLUDE})
- message(FATAL_ERROR "MKL include dir not found!")
- endif()
- if(NOT DEFINED ENV{MKL_LIB})
- message(FATAL_ERROR "MKL lib dir not found!")
- endif()
- # mkl
- target_include_directories(power_grid_model INTERFACE $ENV{MKL_INCLUDE})
- list(APPEND CMAKE_PREFIX_PATH "$ENV{MKL_LIB}")
- find_library(MKL_RT mkl_rt)
- message("${MKL_RT}")
- target_compile_definitions(power_grid_model INTERFACE POWER_GRID_MODEL_USE_MKL=1)
- target_link_libraries(power_grid_model INTERFACE "${MKL_RT}")
- endif()
+if (CMAKE_CXX_COMPILER_ID STREQUAL "GNU")
+ # using GCC
+ add_compile_options(-Wno-maybe-uninitialized) # bug in boost
endif()
-
# get tests
add_subdirectory(tests/cpp_unit_tests)
@@ -85,3 +56,4 @@ if(DEFINED POWER_GRID_MODEL_BUILD_BENCHMARK AND POWER_GRID_MODEL_BUILD_BENCHMARK
add_subdirectory(tests/benchmark_cpp)
endif()
+
diff --git a/CMakeSettings.json b/CMakeSettings.json
index 5196a06947..cc7d611e8a 100644
--- a/CMakeSettings.json
+++ b/CMakeSettings.json
@@ -23,54 +23,6 @@
"ctestCommandArgs": "",
"inheritEnvironments": [ "msvc_x64_x64" ],
"cmakeToolchain": "${env.VCPKG_ROOT}\\scripts\\buildsystems\\vcpkg.cmake"
- },
- {
- "name": "x64-Debug-MKL",
- "generator": "Ninja",
- "configurationType": "Debug",
- "buildRoot": "${projectDir}\\cpp_build\\${name}",
- "installRoot": "${projectDir}\\install\\${name}",
- "cmakeCommandArgs": "-DPOWER_GRID_MODEL_BUILD_BENCHMARK=1 -DPOWER_GRID_MODEL_SPARSE_SOLVER=MKL",
- "buildCommandArgs": "-v",
- "ctestCommandArgs": "",
- "inheritEnvironments": [ "msvc_x64_x64" ],
- "cmakeToolchain": "${env.VCPKG_ROOT}\\scripts\\buildsystems\\vcpkg.cmake"
- },
- {
- "name": "x64-Release-MKL",
- "generator": "Ninja",
- "configurationType": "RelWithDebInfo",
- "buildRoot": "${projectDir}\\cpp_build\\${name}",
- "installRoot": "${projectDir}\\install\\${name}",
- "cmakeCommandArgs": "-DPOWER_GRID_MODEL_BUILD_BENCHMARK=1 -DPOWER_GRID_MODEL_SPARSE_SOLVER=MKL",
- "buildCommandArgs": "-v",
- "ctestCommandArgs": "",
- "inheritEnvironments": [ "msvc_x64_x64" ],
- "cmakeToolchain": "${env.VCPKG_ROOT}\\scripts\\buildsystems\\vcpkg.cmake"
- },
- {
- "name": "x64-Debug-MKL-at-runtime",
- "generator": "Ninja",
- "configurationType": "Debug",
- "buildRoot": "${projectDir}\\cpp_build\\${name}",
- "installRoot": "${projectDir}\\install\\${name}",
- "cmakeCommandArgs": "-DPOWER_GRID_MODEL_BUILD_BENCHMARK=1 -DPOWER_GRID_MODEL_SPARSE_SOLVER=MKL_RUNTIME",
- "buildCommandArgs": "-v",
- "ctestCommandArgs": "",
- "cmakeToolchain": "${env.VCPKG_ROOT}\\scripts\\buildsystems\\vcpkg.cmake",
- "inheritEnvironments": [ "msvc_x64_x64" ]
- },
- {
- "name": "x64-Release-MKL-at-runtime",
- "generator": "Ninja",
- "configurationType": "RelWithDebInfo",
- "buildRoot": "${projectDir}\\cpp_build\\${name}",
- "installRoot": "${projectDir}\\install\\${name}",
- "cmakeCommandArgs": "-DPOWER_GRID_MODEL_BUILD_BENCHMARK=1 -DPOWER_GRID_MODEL_SPARSE_SOLVER=MKL_RUNTIME",
- "buildCommandArgs": "-v",
- "ctestCommandArgs": "",
- "cmakeToolchain": "${env.VCPKG_ROOT}\\scripts\\buildsystems\\vcpkg.cmake",
- "inheritEnvironments": [ "msvc_x64_x64" ]
}
]
}
\ No newline at end of file
diff --git a/README.md b/README.md
index ac8c34f460..51c0528f86 100644
--- a/README.md
+++ b/README.md
@@ -29,39 +29,8 @@ Currently, it supports the following calculations:
# Installation
-## Runtime Dependencies
+## Install from PyPI
-The only Python runtime dependency is
-[numpy](https://numpy.org/). It will be automatically installed as the requirements.
-Moreover, the library optionally depends on
-[Intel Math Kernel Library (mkl)](https://software.intel.com/content/www/us/en/develop/tools/oneapi/components/onemkl.html),
-for its [PARDISO](https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-iface.html) sparse solver.
-It is recommended to install `mkl` because it gives huge performance boosts.
-
-The easiest way to install `mkl` is using `pip` or `conda`:
-
-```
-pip install mkl
-```
-
-or
-
-```
-conda install -c conda-forge mkl
-```
-
-You need to add the path to the `mkl` runtime file `libmkl_rt.so` or `mkl_rt.dll` the environment variable
-`LD_LIBRARY_PATH` in Linux or `Path` in Windows (`conda` does this automatically in the environment).
-If the library can find `mkl` runtime, it uses it as the sparse solver.
-It is recommended to set the environment variable `MKL_THREADING_LAYER` to `SEQUENTIAL`,
-as multi-threading is handled in a higher level.
-If the library cannot find `mkl` runtime, it will fall back to an internally built-in (and much slower)
-[Eigen SparseLU](https://eigen.tuxfamily.org/dox/classEigen_1_1SparseLU.html) solver.
-
-## Install from Pre-built Binary Package
-
-The `power-grid-model` python package is pre-built for Windows, Linux, and macOS (both Intel and Arm-based),
-for Python version 3.8, 3.9, and 3.10.
You can directly install the package from PyPI.
```
@@ -208,13 +177,6 @@ which are licensed under their own respective Open-Source licenses.
SPDX-License-Identifier headers are used to show which license is applicable.
The concerning license files can be found in the LICENSES directory.
-## Intel Math Kernel Library License
-
-The `power-grid-model` does not bundle or redistribute any MKL runtime library.
-It only detects if MKL library is installed in the target system.
-If so, it will use the library to accelerate the calculation.
-The user is responsible to acquire a suitable MKL license.
-
# Contributing
Please read [CODE_OF_CONDUCT](CODE_OF_CONDUCT.md) and [CONTRIBUTING](CONTRIBUTING.md) for details on the process
for submitting pull requests to us.
diff --git a/build.sh b/build.sh
index bfe104ec28..ed6b20e6fa 100755
--- a/build.sh
+++ b/build.sh
@@ -7,7 +7,7 @@
set -e
usage() {
- echo "$0 {Debug/Release} {EIGEN,MKL,MKL_RUNTIME}"
+ echo "$0 {Debug/Release}"
}
if [ ! "$1" = "Debug" ] && [ ! "$1" = "Release" ]; then
@@ -16,13 +16,13 @@ if [ ! "$1" = "Debug" ] && [ ! "$1" = "Release" ]; then
exit 1;
fi
-if [ ! "$2" = "EIGEN" ] && [ ! "$2" = "MKL" ] && [ ! "$2" = "MKL_RUNTIME" ]; then
- echo "Missing second argument"
- usage
- exit 1;
+if [[ $2 == "Coverage" ]]; then
+ BUILD_COVERAGE=-DPOWER_GRID_MODEL_COVERAGE=1
+else
+ BUILD_COVERAGE=
fi
-BUILD_DIR=cpp_build_$1_$2
+BUILD_DIR=cpp_build_$1
echo "Build dir: ${BUILD_DIR}"
if [[ ! -z "${VCPKG_ROOT}" ]]; then
@@ -41,27 +41,19 @@ cd ${BUILD_DIR}
# generate
cmake .. -GNinja \
-DCMAKE_BUILD_TYPE=$1 \
- -DPOWER_GRID_MODEL_SPARSE_SOLVER=$2 \
${PATH_FOR_CMAKE} \
- -DPOWER_GRID_MODEL_BUILD_BENCHMARK=1
+ -DPOWER_GRID_MODEL_BUILD_BENCHMARK=1 \
+ ${BUILD_COVERAGE}
# build
VERBOSE=1 cmake --build .
# test
+./tests/cpp_unit_tests/power_grid_model_unit_tests
-if [[ $2 == "MKL_RUNTIME" ]]; then
- LD_LIBRARY_PATH= ./tests/cpp_unit_tests/power_grid_model_unit_tests
- LD_LIBRARY_PATH= POWER_GRID_MODEL_SPARSE_SOLVER=MKL ./tests/cpp_unit_tests/power_grid_model_unit_tests
- LD_LIBRARY_PATH=${MKL_LIB} ./tests/cpp_unit_tests/power_grid_model_unit_tests
- LD_LIBRARY_PATH=${MKL_LIB} POWER_GRID_MODEL_SPARSE_SOLVER=MKL ./tests/cpp_unit_tests/power_grid_model_unit_tests
- POWER_GRID_MODEL_SPARSE_SOLVER=EIGEN ./tests/cpp_unit_tests/power_grid_model_unit_tests
-else
- LD_LIBRARY_PATH=${MKL_LIB} ./tests/cpp_unit_tests/power_grid_model_unit_tests
-fi
cd ..
# test coverage report for debug build and for linux
-if [[ "$1" = "Debug" ]] && [[ $3 == "Coverage" ]]; then
+if [[ "$1" = "Debug" ]] && [[ $2 == "Coverage" ]]; then
echo "Generating coverage report..."
if [[ ${CXX} == "clang++"* ]]; then
GCOV_TOOL="--gcov-tool llvm-gcov.sh"
diff --git a/docs/build-guide.md b/docs/build-guide.md
index ed11569ad7..062ce608e4 100644
--- a/docs/build-guide.md
+++ b/docs/build-guide.md
@@ -9,7 +9,7 @@ SPDX-License-Identifier: MPL-2.0
This document explains how you can build this library from source, including some examples of build environment. In this
repository there are two builds:
-* A `power-grid-model` [pip](https://pip.pypa.io/en/stable/) Python package with C++ extension as the calculation core.
+* A `power-grid-model` [pip](https://pypi.org/project/power-grid-model/) Python package with C++ extension as the calculation core.
* A [CMake](https://cmake.org/) C++ project to build native C++ unit tests and/or performance benchmark.
* In the future, the following two libraries might be provided:
* A C++ header-only library of the calculation core
@@ -23,9 +23,8 @@ section a list of general requirements are given. After this section there are e
## Architecture Support
-This library is written and tested on `x86_64` and `aarch64` architecture. Building the library in `x86_32` might be working, but is
-not tested. The MKL library is only
-available for `x86_64`. So only Eigen sparse solver is available for `aarch64`.
+This library is written and tested on `x86_64` and `arm64` architecture. Building the library in `x86_32` might be working, but is
+not tested.
The source code is written with the mindset of ISO standard C++ only, i.e. avoid compiler-extension or platform-specific
features as much as possible. In this way the effort to port the library to other platform/architecture might be
@@ -66,7 +65,6 @@ The table below shows the C++ build dependencies
| [eigen3](https://eigen.tuxfamily.org/) | Define environment variable `EIGEN_INCLUDE` to the include folder of `eigen3`\* | CMake needs to be able find `eigen3` | header-only | [Mozilla Public License, version 2.0](https://www.mozilla.org/en-US/MPL/2.0/) |
| [Catch2](https://github.com/catchorg/Catch2) | None | CMake needs to be able find `Catch2` | header-only | [Boost Software License 1.0](https://github.com/catchorg/Catch2/blob/devel/LICENSE.txt) |
| [nlohmann-json](https://github.com/nlohmann/json) | None | CMake needs to be able find `nlohmann_json` | header-only | [MIT](https://github.com/nlohmann/json/blob/develop/LICENSE.MIT) |
-| [MKL](https://software.intel.com/content/www/us/en/develop/tools/oneapi/components/onemkl.html) | Add path to MKL runtime library (`libmkl_rt.so` or `mkl_rt.dll`) into `PATH` (Windows) or `LD_LIBRARY_PATH` (Linux) |
- Define environment variable `MKL_INCLUDE` to the include folder
- Define environment variable `MKL_LIB` to the lib folder (`.lib` or `.so`)
- Add path to MKL runtime library (`libmkl_rt.so` or `mkl_rt.dll`) into `PATH` (Windows) or `LD_LIBRARY_PATH` (Linux)
| Optional for use of MKL PARDISO sparse solver. | [Intel Simplified Software License (proprietary license)](https://www.intel.com/content/www/us/en/developer/articles/license/onemkl-license-faq.html) |
\* The environment variables should point to the root include folder of the library, not a subfolder. For example in the
path `BOOST_INCLUDE` there should be a folder called `boost` which has all the `boost` header files.
@@ -117,12 +115,12 @@ library `power_grid_model`. There are two sub-project defined in the root cmake
In principle, you can use any C++ IDE with cmake and ninja support to develop the C++ project. When you
use `cmake build` for the root cmake file, the following additional options are available besides the standard cmake
-option. If no option is defined, the cmake project will build the unit tests with *Eigen sparse solver*.
+option. If no option is defined, the cmake project will build the unit tests.
-| Option | Description |
-|------------------------------------|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
-| `POWER_GRID_MODEL_SPARSE_SOLVER` | Specify which sparse solver to use during the build. There are three options:
`EIGEN`: use built-in `SparseLU` solver of `eigen3`.
`MKL`: use MKL PARDISO solver and link MKL library at link-time.
`MKL_RUNTIME`: build the project using both Eigen sparse solver and MKL PARDISO solver. The MKL library `mkl_rt` is loaded at runtime. If the MKL library cannot be found, it falls back to the built-in Eigen sparse solver. This is the build option for Python package. |
-| `POWER_GRID_MODEL_BUILD_BENCHMARK` | When set to `1`, build the both sub-projects: unit test and benchmarks. Otherwise only build the unit test. |
+| Option | Description |
+|------------------------------------|-------------------------------------------------------------------------------------------------------------|
+| `POWER_GRID_MODEL_BUILD_BENCHMARK` | When set to `1`, build the both sub-projects: unit test and benchmarks. Otherwise only build the unit test. |
+| `POWER_GRID_MODEL_COVERAGE` | When set to `1`, build with test coverage. This is only applicable for Linux. |
# Example Setup for Ubuntu 20.04 (in WSL or physical/virtual machine)
@@ -140,11 +138,7 @@ export VCPKG_FEATURE_FLAGS=-binarycaching
export VCPKG_ROOT=${HOME}/vcpkg
export EIGEN_INCLUDE=${VCPKG_ROOT}/installed/x64-linux/include/eigen3
export BOOST_INCLUDE=${VCPKG_ROOT}/installed/x64-linux/include
-export MKL_THREADING_LAYER=SEQUENTIAL
-export MKL_INTERFACE_LAYER=LP64
-export LD_LIBRARY_PATH=${HOME}/.local/lib:${LD_LIBRARY_PATH}
-export MKL_LIB=${HOME}/.local/lib
-export MKL_INCLUDE=${HOME}/.local/include
+export LLVM_COV=llvm-cov-11 # or llvm-cov-10 if you use clang 10
```
## Ubuntu Software Packages
@@ -154,7 +148,7 @@ compiler.
```shell
sudo apt update && sudo apt -y upgrade
-sudo apt install -y curl zip unzip tar git build-essential lcov gcc g++ gcc-10 g++-10 clang clang-11 make cmake gdb ninja-build pkg-config python3.9 python3.9-dev python3.9-venv python3-pip
+sudo apt install -y curl zip unzip tar git build-essential lcov gcovr gcc g++ gcc-10 g++-10 clang clang-11 make cmake gdb ninja-build pkg-config python3.9 python3.9-dev python3.9-venv python3-pip
```
## C++ package manager: vcpkg
@@ -171,17 +165,6 @@ cd vcpkg
**The installation of `boost` will take a long time, be patient**
-## MKL
-
-The easiest way to install `mkl` is through `pip`. Since you can have multiple Python environment which may depend on a
-single `mkl` library. It is better to install `mkl` without virtual environment. Make sure your current shell is **NOT**
-in a virtual environment (no parentheses before the prompt). Without `sudo`, it will install the `mkl` into the
-folder `${HOME}/.local`.
-
-```shell
-python3.9 -m pip install mkl mkl-include mkl-devel
-```
-
## Build Python Library from Source
It is recommended to create a virtual environment. Clone repository, create and activate virtual environment, and
@@ -205,48 +188,28 @@ pytest
## Build CMake Project
There is a convenient shell script to build the cmake project:
-[`build.sh`](../build.sh). You can study the file and write your own build script. Six configurations are pre-defined
+[`build.sh`](../build.sh). You can study the file and write your own build script. Four configurations are pre-defined
for two input arguments, which will be passed into `cmake`. It includes debug or release build, as well as the option to
-use MKL at link-time, or at runtime, or not at all.
+build test coverage or not.
* Option 1
* Debug
* Release
-* Option 2
- * EIGEN
- * MKL
- * MKL_RUNTIME
+* Option 2 (optional)
+ * Coverage (if specified, the test coverage will be run and a web report will be generated in `cpp_cov_html` folder.)
-As an example, go to the root folder of repo. Use the following command to build the project with MKL at link-time and
-benchmark:
+As an example, go to the root folder of repo. Use the following command to build the project in release mode:
```shell
-./build.sh Release MKL
+./build.sh Release
```
One can run the unit tests and benchmark by:
```shell
-./cpp_build_Relase_MKL/tests/cpp_unit_tests/power_grid_model_unit_tests
-
-./cpp_build_Relase_MKL/tests/benchmark_cpp/power_grid_model_benchmark_cpp
-```
-
-## Optional: to generate a C++ test coverage report
-
-If you want to generate a C++ test coverage report you should do a few extra steps.
-
-1. Install the following extra packages from Ubuntu.
-
-```shell
-sudo apt install lcov gcovr
-```
-
-2. Change your compiler in `${HOME}/.bashrc`
+./cpp_build_Release/tests/cpp_unit_tests/power_grid_model_unit_tests
-```shell
-export CXX=g++-10 # or clang++-10, or g++-9, or g++-10
-export CC=gcc-10 # or clang-10, or gcc-9, or gcc-10
+./cpp_build_Release/tests/benchmark_cpp/power_grid_model_benchmark_cpp
```
# Example Setup for Windows 10
@@ -263,12 +226,6 @@ Define the following environment variables in user wide.
| VCPKG_FEATURE_FLAGS | -binarycaching |
| EIGEN_INCLUDE | \installed\x64-windows\include\eigen3 |
| BOOST_INCLUDE | \installed\x64-windows\include |
-| MKL_LIB | \envs\mkl-base\Library\lib |
-| MKL_INCLUDE | \envs\mkl-base\Library\include |
-| MKL_THREADING_LAYER | SEQUENTIAL |
-| MKL_INTERFACE_LAYER | LP64 |
-
-Further, append `\envs\mkl-base\Library\bin` to environment variable `PATH` in user wide.
## Software Toolchains
@@ -296,18 +253,6 @@ cd vcpkg
./vcpkg install eigen3 catch2 boost nlohmann-json
```
-## MKL
-
-You can install `mkl` from `conda`. Since you might want to use `mkl` in different projects, it is recommended to (only)
-install `mkl` in a separate conda environment. In this example the environment is called `mkl-base`. Open a Miniconda
-PowerShell Prompt.
-
-```shell
-conda create -n mkl-base python=3.9
-conda activate mkl-base
-conda install mkl mkl-include mkl-devel
-```
-
## Build Python Library from Source
It is recommended to create a `conda` environment. Clone repository, create and activate `conda` environment, and
@@ -336,15 +281,10 @@ pytest
If you have installed Visual Studio 2019/2022 (not the build tools), you can open the repo folder as a cmake project.
The IDE should be able to automatically detect the Visual Studio cmake configuration file
-[`CMakeSettings.json`](../CMakeSettings.json). Six configurations are pre-defined. It includes debug or release build,
-as well as the option to use MKL at link-time, or at runtime, or not at all.
+[`CMakeSettings.json`](../CMakeSettings.json). Two configurations are pre-defined. It includes debug or release build.
* `x64-Debug`
-* `x64-Debug-MKL`
-* `x64-Debug-MKL-at-runtime`
* `x64-Release`
-* `x64-Release-MKL`
-* `x64-Release-MKL-at-runtime`
# Example Setup for macOS (Big Sur)
@@ -362,12 +302,6 @@ export VCPKG_FEATURE_FLAGS=-binarycaching
export VCPKG_ROOT=${HOME}/vcpkg
export EIGEN_INCLUDE=${VCPKG_ROOT}/installed/x64-osx/include/eigen3 # use arm64-osx in m1 Mac
export BOOST_INCLUDE=${VCPKG_ROOT}/installed/x64-osx/include # use arm64-osx in m1 Mac
-# Skip the following for Mac M1 (due to MKL)
-export MKL_THREADING_LAYER=SEQUENTIAL
-export MKL_INTERFACE_LAYER=LP64
-export LD_LIBRARY_PATH=/Library/Frameworks/Python.framework/Versions/3.9/lib:${LD_LIBRARY_PATH}
-export MKL_LIB=/Library/Frameworks/Python.framework/Versions/3.9/lib
-export MKL_INCLUDE=/Library/Frameworks/Python.framework/Versions/3.9/include
```
## macOS Software Packages
@@ -392,19 +326,6 @@ cd vcpkg
**The installation of `boost` will take a long time, be patient**
-## MKL
-
-Skip this step for Mac M1.
-
-The easiest way to install `mkl` is through `pip`. Since you can have multiple Python environment which may depend on a
-single `mkl` library. It is better to install `mkl` without virtual environment. Make sure your current shell is **NOT**
-in a virtual environment (no parentheses before the prompt). Without `sudo`, it will install the `mkl` into the
-folder `${HOME}/.local`.
-
-```shell
-python3.9 -m pip install mkl mkl-include mkl-devel
-```
-
## Build Python Library from Source
It is recommended to create a virtual environment. Clone repository, create and activate virtual environment, and
@@ -428,29 +349,26 @@ pytest
## Build CMake Project
There is a convenient shell script to build the cmake project:
-[`build.sh`](../build.sh). You can study the file and write your own build script. Six configurations are pre-defined
-for two input arguments, which will be passed into `cmake`. It includes debug or release build, as well as the option to
-use MKL at link-time, or at runtime, or not at all.
+[`build.sh`](../build.sh). You can study the file and write your own build script. Two configurations are pre-defined
+with one input argument, which will be passed into `cmake`. It includes debug or release build.
* Option 1
* Debug
* Release
-* Option 2
- * EIGEN
- * MKL
- * MKL_RUNTIME
-As an example, go to the root folder of repo. Use the following command to build the project with MKL at link-time and
-benchmark:
+**Note: the test coverage option is not supported in macOS.**
+
+As an example, go to the root folder of repo. Use the following command to build the project in release mode:
```shell
-./build.sh Release MKL
+./build.sh Release
```
One can run the unit tests and benchmark by:
```shell
-./cpp_build_Release_MKL/tests/cpp_unit_tests/power_grid_model_unit_tests
+./cpp_build_Release/tests/cpp_unit_tests/power_grid_model_unit_tests
+
+./cpp_build_Release/tests/benchmark_cpp/power_grid_model_benchmark_cpp
+
-./cpp_build_Release_MKL/tests/benchmark_cpp/power_grid_model_benchmark_cpp
-```
diff --git a/include/power_grid_model/calculation_parameters.hpp b/include/power_grid_model/calculation_parameters.hpp
index db588b31d0..b9a5b3e427 100644
--- a/include/power_grid_model/calculation_parameters.hpp
+++ b/include/power_grid_model/calculation_parameters.hpp
@@ -91,6 +91,7 @@ struct MathModelTopology {
Idx slack_bus_;
std::vector phase_shift;
std::vector branch_bus_idx;
+ std::vector fill_in;
IdxVector source_bus_indptr;
IdxVector shunt_bus_indptr;
IdxVector load_gen_bus_indptr;
diff --git a/include/power_grid_model/enum.hpp b/include/power_grid_model/enum.hpp
index b3298a0409..e822769c3a 100644
--- a/include/power_grid_model/enum.hpp
+++ b/include/power_grid_model/enum.hpp
@@ -54,7 +54,15 @@ enum class ComponentType : IntS {
// for 0b00 - 0b11
// 0bXY, where X, Y means the from(0)/to(1) side of branch
// i.e. 0b01 is the branch element for Yft
-enum class YBusElementType : IntS { bff = 0b00, bft = 0b01, btf = 0b10, btt = 0b11, shunt = 0b100 };
+enum class YBusElementType : IntS {
+ bff = 0b00,
+ bft = 0b01,
+ btf = 0b10,
+ btt = 0b11,
+ shunt = 0b100,
+ fill_in_ft = 0b101,
+ fill_in_tf = 0b110
+};
} // namespace power_grid_model
diff --git a/include/power_grid_model/main_model.hpp b/include/power_grid_model/main_model.hpp
index 99deccb749..f5f55d6025 100644
--- a/include/power_grid_model/main_model.hpp
+++ b/include/power_grid_model/main_model.hpp
@@ -1166,24 +1166,34 @@ class MainModelImpl, ComponentLis
template
void prepare_solvers() {
std::vector>& solvers = get_solvers();
+ // also get the vector of other solvers (sym -> asym, or asym -> sym)
+ std::vector>& other_solvers = get_solvers();
// rebuild topology if needed
if (!is_topology_up_to_date_) {
rebuild_topology();
}
// if solvers do not exist, build them
if (n_math_solvers_ != (Idx)solvers.size()) {
+ // check if other (sym/asym) solver exist
+ bool const other_solver_exist = (n_math_solvers_ == (Idx)other_solvers.size());
assert(solvers.size() == 0);
solvers.reserve(n_math_solvers_);
// get param, will be consumed
std::vector> math_params = get_math_param();
- // use transform to build
- std::transform(math_topology_.cbegin(), math_topology_.cend(), math_params.begin(),
- std::back_inserter(solvers),
- [](std::shared_ptr const& topo_ptr, MathModelParam& param) {
- return MathSolver{topo_ptr,
- // move parameter into a shared ownership for the math solver
- std::make_shared const>(std::move(param))};
- });
+ // loop to build
+ for (Idx i = 0; i != n_math_solvers_; ++i) {
+ // if other solver exists, construct from existing y bus struct
+ if (other_solver_exist) {
+ solvers.emplace_back(math_topology_[i],
+ std::make_shared const>(std::move(math_params[i])),
+ other_solvers[i].shared_y_bus_struct());
+ }
+ // else construct from scratch
+ else {
+ solvers.emplace_back(math_topology_[i],
+ std::make_shared const>(std::move(math_params[i])));
+ }
+ }
}
// if parameters are not up to date, update them
else if (!is_parameter_up_to_date()) {
diff --git a/include/power_grid_model/math_solver/block_matrix.hpp b/include/power_grid_model/math_solver/block_matrix.hpp
index d059e52f79..9f9bb9f692 100644
--- a/include/power_grid_model/math_solver/block_matrix.hpp
+++ b/include/power_grid_model/math_solver/block_matrix.hpp
@@ -17,48 +17,58 @@ namespace power_grid_model {
// hide implementation in inside namespace
namespace math_model_impl {
-template || std::is_same_v>>
-class BlockEntry {
+template >>
+struct block_trait {
+ static constexpr int n_row = sym ? n_sub_block : n_sub_block * 3;
+ static constexpr int n_col = is_tensor ? (sym ? n_sub_block : n_sub_block * 3) : 1;
+ using ArrayType = Eigen::Array;
+};
+
+template
+class Block : public block_trait::ArrayType {
public:
- static constexpr int scalar_size = sym ? 1 : 3;
- static constexpr int block_size = scalar_size * n_sub_block;
- static constexpr int size = block_size * block_size;
- static constexpr int size_in_double = (std::is_same_v ? 1 : 2) * size;
+ using ArrayType = typename block_trait::ArrayType;
+ using ArrayType::operator();
- using ArrayType = Eigen::Array;
- using GetterType = std::conditional_t>;
+ // default zero
+ Block() : ArrayType{ArrayType::Zero()} {};
+ // eigen expression
+ template
+ Block(Eigen::ArrayBase const& other) : ArrayType{other} {
+ }
+ template
+ Block& operator=(Eigen::ArrayBase const& other) {
+ this->ArrayType::operator=(other);
+ return *this;
+ }
- protected:
- // get position, can be 4 values
- // 0, 0: upper left
- // 0, 1: upper right
- // 1, 0: lower left
- // 1, 1: lower right
- template
- GetterType get_val() {
- if constexpr (sym) {
- return data_(row, col);
+ template
+ static auto get_asym_row_idx() {
+ return Eigen::seqN(Eigen::fix, Eigen::fix<3>);
+ }
+ template
+ static auto get_asym_col_idx() {
+ if constexpr (is_tensor) {
+ return Eigen::seqN(Eigen::fix, Eigen::fix<3>);
}
else {
- return data_.template block(row * scalar_size, col * scalar_size);
+ return Eigen::seqN(Eigen::fix<0>, Eigen::fix<1>);
}
}
- private:
- ArrayType data_{ArrayType::Zero()};
-};
+ template
+ using GetterType =
+ std::conditional_t()(get_asym_row_idx(), get_asym_col_idx()))>;
-template class BlockType>
-struct block_entry_trait {
- template
- struct internal_trait {
- static_assert(sizeof(BlockType) == sizeof(double[BlockType::size_in_double]));
- static_assert(alignof(BlockType) >= alignof(double[BlockType::size_in_double]));
- static_assert(std::is_standard_layout_v>);
- };
- static constexpr internal_trait sym_trait{};
- static constexpr internal_trait asym_trait{};
+ template
+ GetterType get_val() {
+ if constexpr (sym) {
+ return (*this)(r, c);
+ }
+ else {
+ return (*this)(get_asym_row_idx(), get_asym_col_idx());
+ }
+ }
};
} // namespace math_model_impl
diff --git a/include/power_grid_model/math_solver/bsr_solver.hpp b/include/power_grid_model/math_solver/bsr_solver.hpp
deleted file mode 100644
index 69e54bb332..0000000000
--- a/include/power_grid_model/math_solver/bsr_solver.hpp
+++ /dev/null
@@ -1,102 +0,0 @@
-// SPDX-FileCopyrightText: 2022 Contributors to the Power Grid Model project
-//
-// SPDX-License-Identifier: MPL-2.0
-
-#pragma once
-#ifndef POWER_GRID_MODEL_MATH_SOLVER_BSR_SOLVER_HPP
-#define POWER_GRID_MODEL_MATH_SOLVER_BSR_SOLVER_HPP
-// BSR solver adaptor based on which solver is selected
-
-// if use mkl at runtime, also set use mkl
-#ifdef POWER_GRID_MODEL_USE_MKL_AT_RUNTIME
-#ifndef POWER_GRID_MODEL_USE_MKL
-#define POWER_GRID_MODEL_USE_MKL
-#endif
-#endif
-
-#include "../exception.hpp"
-#include "../power_grid_model.hpp"
-// include both solver
-#include
-
-#include "eigen_superlu_solver.hpp"
-#include "mkl_pardiso_solver.hpp"
-
-namespace power_grid_model {
-
-#ifdef POWER_GRID_MODEL_USE_MKL_AT_RUNTIME
-
-// wrapper solver
-template
-class BSRSolver {
- using SolverType = std::variant, EigenSuperLUSolver>;
-
- public:
- template , BSRSolver> && ...)>>
- explicit BSRSolver(Args&&... args)
- : // template Args&&... args perfect forwarding
- solver_{get_pardiso_handle().has_pardiso ?
- // initialize with pardiso
- SolverType{std::in_place_type>, std::forward(args)...}
- :
- // initialize with eigen
- SolverType{std::in_place_type>, std::forward(args)...}} {
- }
-
- template
- void solve(Args&&... args) {
- // template Args&&... args perfect forwarding
- // call solver based on use mkl
- if (get_pardiso_handle().has_pardiso) {
- std::get>(solver_).solve(std::forward(args)...);
- }
- else {
- std::get>(solver_).solve(std::forward(args)...);
- }
- }
-
- template
- void prefactorize(Args&&... args) {
- // template Args&&... args perfect forwarding
- if (get_pardiso_handle().has_pardiso) {
- std::get>(solver_).prefactorize(std::forward(args)...);
- }
- else {
- std::get>(solver_).prefactorize(std::forward(args)...);
- }
- }
-
- void invalidate_prefactorization() {
- if (get_pardiso_handle().has_pardiso) {
- std::get>(solver_).invalidate_prefactorization();
- }
- else {
- std::get>(solver_).invalidate_prefactorization();
- }
- }
-
- private:
- SolverType solver_;
-};
-
-#else
-
-// select solver statically at compile time
-#ifdef POWER_GRID_MODEL_USE_MKL
-// use mkl solver
-template
-using BSRSolver = PARDISOSolver;
-#else // !POWER_GRID_MODEL_USE_MKL
-// use eigen solver
-template
-using BSRSolver = EigenSuperLUSolver;
-#endif // POWER_GRID_MODEL_USE_MKL
-
-#endif // POWER_GRID_MODEL_USE_MKL_AT_RUNTIME
-
-} // namespace power_grid_model
-
-#endif
diff --git a/include/power_grid_model/math_solver/eigen_superlu_solver.hpp b/include/power_grid_model/math_solver/eigen_superlu_solver.hpp
deleted file mode 100644
index 25ede61917..0000000000
--- a/include/power_grid_model/math_solver/eigen_superlu_solver.hpp
+++ /dev/null
@@ -1,165 +0,0 @@
-// SPDX-FileCopyrightText: 2022 Contributors to the Power Grid Model project
-//
-// SPDX-License-Identifier: MPL-2.0
-
-#pragma once
-#ifndef POWER_GRID_MODEL_MATH_SOLVER_EIGEN_SUPERLU_SOLVER_HPP
-#define POWER_GRID_MODEL_MATH_SOLVER_EIGEN_SUPERLU_SOLVER_HPP
-// BSR solver using Eigen SuperLU
-// https://eigen.tuxfamily.org/dox/group__Sparse__chapter.html
-
-// only enable if POWER_GRID_MODEL_USE_MKL_AT_RUNTIME is defined
-// or POWER_GRID_MODEL_USE_MKL is not defined
-#if defined(POWER_GRID_MODEL_USE_MKL_AT_RUNTIME) || !defined(POWER_GRID_MODEL_USE_MKL)
-
-#include
-#include
-
-#include "../exception.hpp"
-#include "../power_grid_model.hpp"
-
-namespace power_grid_model {
-
-template
-class EigenSuperLUSolver final {
- private:
- static_assert(std::is_same_v || std::is_same_v);
- using SparseMatrix = Eigen::SparseMatrix;
- using SparseIndexMatrix = Eigen::SparseMatrix;
- using SparseSolver = Eigen::SparseLU>;
- using BufferVector = Eigen::Map, Eigen::Aligned8>;
- using SparseIdxEntry = Eigen::Triplet;
-
- struct BSRHandle {
- Idx matrix_size_in_block;
- Idx block_size;
- Idx nnz_block;
- Idx block_capacity;
- Idx nnz;
- Idx matrix_size;
- std::shared_ptr data_mapping;
- };
-
- public:
- EigenSuperLUSolver(Idx matrix_size_in_block, // size of matrix in block form
- Idx block_size, // size of block
- std::shared_ptr const& ia, // indptr
- std::shared_ptr const& ja // column indices
- ) {
- bsr_handle_.matrix_size_in_block = matrix_size_in_block;
- bsr_handle_.block_size = block_size;
- bsr_handle_.nnz_block = ia->back();
- bsr_handle_.block_capacity = block_size * block_size;
- bsr_handle_.nnz = bsr_handle_.nnz_block * bsr_handle_.block_capacity;
- bsr_handle_.matrix_size = matrix_size_in_block * block_size;
- initialize_matrix(*ia, *ja);
- initialize_solver();
- }
-
- // copy construction
- EigenSuperLUSolver(EigenSuperLUSolver const& other)
- : bsr_handle_{other.bsr_handle_}, sparse_matrix_{other.sparse_matrix_} {
- initialize_solver();
- }
-
- // copy assignment
- EigenSuperLUSolver& operator=(EigenSuperLUSolver const& other) {
- // copy new handle
- bsr_handle_ = other.bsr_handle_;
- // copy new matrix
- sparse_matrix_ = other.sparse_matrix_;
- // re initialize
- initialize_solver();
- prefactorized_ = false;
- return *this;
- }
-
- void solve(void const* data, void* b, void* x, bool use_prefactorization = false) {
- // reset possible pre-factorization if we are not using prefactorization
- prefactorized_ = prefactorized_ && use_prefactorization;
- // run factorization
- if (!prefactorized_) {
- prefactorize(data);
- }
- // run solve
- BufferVector bm{std::launder(reinterpret_cast(b)), bsr_handle_.matrix_size};
- BufferVector xm{std::launder(reinterpret_cast(x)), bsr_handle_.matrix_size};
- xm = sparse_solver_.solve(bm);
- if (sparse_solver_.info() != Eigen::Success) {
- throw SparseMatrixError{sparse_solver_.info(), sparse_solver_.lastErrorMessage()};
- }
- }
-
- void prefactorize(void const* data) {
- // copy data
- copy_matrix_data(data);
- // factorize
- sparse_solver_.factorize(sparse_matrix_);
- if (sparse_solver_.info() != Eigen::Success) {
- throw SparseMatrixError{sparse_solver_.info(), sparse_solver_.lastErrorMessage()};
- }
- prefactorized_ = true;
- }
-
- void invalidate_prefactorization() {
- prefactorized_ = false;
- }
-
- private:
- BSRHandle bsr_handle_;
- SparseMatrix sparse_matrix_;
- SparseSolver sparse_solver_; // not copyable or movable
- bool prefactorized_{false};
-
- void initialize_matrix(IdxVector const& ia, // indptr
- IdxVector const& ja // column indices
- ) {
- // loop to add triplet
- std::vector idx_list;
- idx_list.reserve(bsr_handle_.nnz);
- // two loops for blocks in csr
- for (Idx bi = 0; bi != bsr_handle_.matrix_size_in_block; ++bi) {
- for (Idx block_ind = ia[bi]; block_ind != ia[bi + 1]; ++block_ind) {
- Idx const bj = ja[block_ind];
- // two loops for rows/columns inside block
- for (Idx ci = 0; ci != bsr_handle_.block_size; ++ci) {
- for (Idx cj = 0; cj != bsr_handle_.block_size; ++cj) {
- Idx const i = bi * bsr_handle_.block_size + ci;
- Idx const j = bj * bsr_handle_.block_size + cj;
- Idx const data_ind = block_ind * bsr_handle_.block_capacity + ci * bsr_handle_.block_size + cj;
- idx_list.emplace_back(i, j, data_ind);
- }
- }
- }
- }
- // idx matrix
- SparseIndexMatrix idx_matrix{bsr_handle_.matrix_size, bsr_handle_.matrix_size};
- idx_matrix.setFromTriplets(idx_list.begin(), idx_list.end());
- assert(idx_matrix.isCompressed());
- // copy idx data
- bsr_handle_.data_mapping =
- std::make_shared(idx_matrix.valuePtr(), idx_matrix.valuePtr() + bsr_handle_.nnz);
- // initialize sparse matrix
- sparse_matrix_ = idx_matrix.cast();
- assert(sparse_matrix_.isCompressed());
- }
-
- void initialize_solver() {
- sparse_solver_.analyzePattern(sparse_matrix_);
- }
-
- void copy_matrix_data(void const* data) {
- T const* const input_ptr = std::launder(reinterpret_cast(data));
- T* const data_ptr = sparse_matrix_.valuePtr();
- std::transform(bsr_handle_.data_mapping->begin(), bsr_handle_.data_mapping->end(), data_ptr,
- [&input_ptr](Idx j) {
- return input_ptr[j];
- });
- }
-};
-
-} // namespace power_grid_model
-
-#endif // defined(POWER_GRID_MODEL_USE_MKL_AT_RUNTIME) || !defined(POWER_GRID_MODEL_USE_MKL)
-
-#endif
diff --git a/include/power_grid_model/math_solver/iterative_current_pf_solver.hpp b/include/power_grid_model/math_solver/iterative_current_pf_solver.hpp
index f3b9c8b669..4da6a64af8 100644
--- a/include/power_grid_model/math_solver/iterative_current_pf_solver.hpp
+++ b/include/power_grid_model/math_solver/iterative_current_pf_solver.hpp
@@ -60,8 +60,8 @@ Iterative Power Flow
#include "../three_phase_tensor.hpp"
#include "../timer.hpp"
#include "block_matrix.hpp"
-#include "bsr_solver.hpp"
#include "iterative_pf_solver.hpp"
+#include "sparse_lu_solver.hpp"
#include "y_bus.hpp"
namespace power_grid_model {
@@ -72,29 +72,36 @@ namespace math_model_impl {
// solver
template
class IterativeCurrentPFSolver : public IterativePFSolver> {
- private:
- // block size 1 for symmetric, 3 for asym
- static constexpr Idx bsr_block_size_ = sym ? 1 : 3;
-
public:
+ using BlockPermArray =
+ typename SparseLUSolver, ComplexValue, ComplexValue>::BlockPermArray;
+
IterativeCurrentPFSolver(YBus const& y_bus, std::shared_ptr const& topo_ptr)
: IterativePFSolver{y_bus, topo_ptr},
- updated_u_(y_bus.size()),
- rhs_(y_bus.size()),
- mat_data_(y_bus.nnz()),
+ rhs_u_(y_bus.size()),
y_data_ptr_(nullptr),
- bsr_solver_{y_bus.size(), bsr_block_size_, y_bus.shared_indptr(), y_bus.shared_indices()} {
+ sparse_solver_{y_bus.shared_indptr_lu(), y_bus.shared_indices_lu(), y_bus.shared_diag_lu()} {
}
// Add source admittance to Y bus and set variable for prepared y bus to true
void initialize_derived_solver(YBus const& y_bus, MathOutput const&) {
IdxVector const& source_bus_indptr = *this->source_bus_indptr_;
ComplexTensorVector const& ydata = y_bus.admittance();
- IdxVector const& bus_entry = y_bus.bus_entry();
- // Build y bus data with source admittance
- // copy y bus data.
+ IdxVector const& bus_entry = y_bus.lu_diag();
+ // if Y bus is not up to date
+ // re-build matrix and prefactorize Build y bus data with source admittance
if (y_data_ptr_ != &y_bus.admittance()) {
- std::copy(ydata.begin(), ydata.end(), mat_data_.begin());
+ ComplexTensorVector mat_data(y_bus.nnz_lu());
+ // copy y bus data
+ std::transform(y_bus.map_lu_y_bus().cbegin(), y_bus.map_lu_y_bus().cend(), mat_data.begin(), [&](Idx k) {
+ if (k == -1) {
+ return ComplexTensor{};
+ }
+ else {
+ return ydata[k];
+ }
+ });
+
// loop bus
for (Idx bus_number = 0; bus_number != this->n_bus_; ++bus_number) {
Idx const data_sequence = bus_entry[bus_number];
@@ -102,10 +109,16 @@ class IterativeCurrentPFSolver : public IterativePFSolvern_bus_);
+ sparse_solver_.prefactorize(mat_data, perm);
+ // move pre-factorized version into shared ptr
+ mat_data_ = std::make_shared const>(std::move(mat_data));
+ perm_ = std::make_shared(std::move(perm));
+ // cache pointer
y_data_ptr_ = &y_bus.admittance();
}
}
@@ -118,7 +131,7 @@ class IterativeCurrentPFSolver : public IterativePFSolver const& load_gen_type = *this->load_gen_type_;
// set rhs to zero for iteration start
- std::fill(rhs_.begin(), rhs_.end(), ComplexValue{0.0});
+ std::fill(rhs_u_.begin(), rhs_u_.end(), ComplexValue{0.0});
// loop buses: i
for (Idx bus_number = 0; bus_number != this->n_bus_; ++bus_number) {
@@ -130,15 +143,16 @@ class IterativeCurrentPFSolver : public IterativePFSolver{input.source[source_number]});
+ rhs_u_[bus_number] += dot(y_bus.math_model_param().source_param[source_number],
+ ComplexValue{input.source[source_number]});
}
}
}
// Solve the linear equations I_inj = YU
+ // inplace
void solve_matrix() {
- bsr_solver_.solve(mat_data_.data(), rhs_.data(), updated_u_.data(), true);
+ sparse_solver_.solve_with_prefactorized_matrix(*mat_data_, *perm_, rhs_u_, rhs_u_);
}
// Find maximum deviation in voltage among all buses
@@ -165,21 +180,22 @@ class IterativeCurrentPFSolver : public IterativePFSolvern_bus_; ++bus_number) {
// Get maximum iteration for a bus
- double const dev = max_val(cabs(updated_u_[bus_number] - u[bus_number]));
+ double const dev = max_val(cabs(rhs_u_[bus_number] - u[bus_number]));
// Keep maximum deviation of all buses
max_dev = std::max(dev, max_dev);
// assign updated values
- u[bus_number] = updated_u_[bus_number];
+ u[bus_number] = rhs_u_[bus_number];
}
return max_dev;
}
private:
- ComplexValueVector updated_u_;
- ComplexValueVector rhs_;
- ComplexTensorVector mat_data_;
+ ComplexValueVector rhs_u_;
+ std::shared_ptr const> mat_data_;
ComplexTensorVector const* y_data_ptr_;
- BSRSolver bsr_solver_;
+ // sparse solver
+ SparseLUSolver, ComplexValue, ComplexValue> sparse_solver_;
+ std::shared_ptr perm_;
};
template class IterativeCurrentPFSolver;
diff --git a/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp b/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp
index 8240db8688..c86178b4ba 100644
--- a/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp
+++ b/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp
@@ -15,7 +15,6 @@ iterative linear state estimation solver
#include "../power_grid_model.hpp"
#include "../three_phase_tensor.hpp"
#include "../timer.hpp"
-#include "bsr_solver.hpp"
#include "y_bus.hpp"
namespace power_grid_model {
@@ -23,30 +22,34 @@ namespace power_grid_model {
// hide implementation in inside namespace
namespace math_model_impl {
-// block class for the unknown vector in state estimation equation
+// block class for the unknown vector and/or right hand side in state estimation equation
template
-struct SEUnknown {
- ComplexValue u; // real unknown
- ComplexValue phi; // artificial unknown
+struct SEUnknown : public Block {
+ template
+ using GetterType = typename Block::template GetterType;
+
+ // eigen expression
+ using Block::Block;
+ using Block::operator=;
+
+ GetterType<0, 0> u() {
+ return this->template get_val<0, 0>();
+ }
+ GetterType<1, 0> phi() {
+ return this->template get_val<1, 0>();
+ }
+
+ GetterType<0, 0> eta() {
+ return this->template get_val<0, 0>();
+ }
+ GetterType<1, 0> tau() {
+ return this->template get_val<1, 0>();
+ }
};
-static_assert(sizeof(SEUnknown) == sizeof(double[4]));
-static_assert(alignof(SEUnknown) == alignof(double[4]));
-static_assert(std::is_standard_layout_v>);
-static_assert(sizeof(SEUnknown) == sizeof(double[12]));
-static_assert(alignof(SEUnknown) >= alignof(double[12]));
-static_assert(std::is_standard_layout_v>);
+
// block class for the right hand side in state estimation equation
template
-struct SERhs {
- ComplexValue eta; // bus voltage, branch flow, shunt flow
- ComplexValue tau; // injection flow, zero injection constraint
-};
-static_assert(sizeof(SERhs) == sizeof(double[4]));
-static_assert(alignof(SERhs) == alignof(double[4]));
-static_assert(std::is_standard_layout_v>);
-static_assert(sizeof(SERhs) == sizeof(double[12]));
-static_assert(alignof(SERhs) >= alignof(double[12]));
-static_assert(std::is_standard_layout_v>);
+using SERhs = SEUnknown;
// class of 2*2 (6*6) se gain block
/*
@@ -56,23 +59,28 @@ static_assert(std::is_standard_layout_v>);
]
*/
template
-class SEGainBlock : public BlockEntry {
+class SEGainBlock : public Block {
public:
- using typename BlockEntry::GetterType;
- GetterType g() {
+ template
+ using GetterType = typename Block::template GetterType;
+
+ // eigen expression
+ using Block::Block;
+ using Block::operator=;
+
+ GetterType<0, 0> g() {
return this->template get_val<0, 0>();
}
- GetterType qh() {
+ GetterType<0, 1> qh() {
return this->template get_val<0, 1>();
}
- GetterType q() {
+ GetterType<1, 0> q() {
return this->template get_val<1, 0>();
}
- GetterType r() {
+ GetterType<1, 1> r() {
return this->template get_val<1, 1>();
}
};
-constexpr block_entry_trait se_gain_trait{};
// processed measurement struct
// combined all measurement of the same quantity
@@ -572,10 +580,10 @@ class IterativeLinearSESolver {
IterativeLinearSESolver(YBus const& y_bus, std::shared_ptr const& topo_ptr)
: n_bus_{y_bus.size()},
math_topo_{topo_ptr},
- data_gain_(y_bus.nnz()),
- x_(y_bus.size()),
- rhs_(y_bus.size()),
- bsr_solver_{y_bus.size(), bsr_block_size_, y_bus.shared_indptr(), y_bus.shared_indices()} {
+ data_gain_(y_bus.nnz_lu()),
+ x_rhs_(y_bus.size()),
+ sparse_solver_{y_bus.shared_indptr_lu(), y_bus.shared_indices_lu(), y_bus.shared_diag_lu()},
+ perm_(y_bus.size()) {
}
MathOutput run_state_estimation(YBus const& y_bus, StateEstimationInput const& input, double err_tol,
@@ -613,7 +621,7 @@ class IterativeLinearSESolver {
prepare_rhs(y_bus, measured_values, output.u);
// solve with prefactorization
sub_timer = Timer(calculation_info, 2225, "Solve sparse linear equation (pre-factorized)");
- bsr_solver_.solve(data_gain_.data(), rhs_.data(), x_.data(), true);
+ sparse_solver_.solve_with_prefactorized_matrix(data_gain_, perm_, x_rhs_, x_rhs_);
sub_timer = Timer(calculation_info, 2226, "Iterate unknown");
max_dev = iterate_unknown(output.u, measured_values.has_angle_measurement());
}
@@ -645,87 +653,101 @@ class IterativeLinearSESolver {
// data for gain matrix
std::vector> data_gain_;
// unknown and rhs
- std::vector> x_;
- std::vector> rhs_;
+ std::vector> x_rhs_;
// solver
- BSRSolver bsr_solver_;
+ SparseLUSolver, SERhs, SEUnknown> sparse_solver_;
+ typename SparseLUSolver, SERhs, SEUnknown>::BlockPermArray perm_;
void prepare_matrix(YBus const& y_bus, MeasuredValues const& measured_value) {
MathModelParam const& param = y_bus.math_model_param();
+ IdxVector const& row_indptr = y_bus.row_indptr_lu();
+ IdxVector const& col_indices = y_bus.col_indices_lu();
// loop data index, all rows and columns
- for (Idx data_idx = 0; data_idx != y_bus.nnz(); ++data_idx) {
- Idx const row = y_bus.row_indices()[data_idx];
- Idx const col = y_bus.col_indices()[data_idx];
- // get a reference and reset block to zero
- SEGainBlock& block = data_gain_[data_idx];
- block = {};
- // fill block with voltage measurement, only diagonal
- if ((row == col) && measured_value.has_voltage(row)) {
- // G += 1.0 / variance
- // for 3x3 tensor, fill diagonal
- block.g() += ComplexTensor{1.0 / measured_value.voltage_var(row)};
- }
- // fill block with branch, shunt measurement
- for (Idx element_idx = y_bus.y_bus_entry_indptr()[data_idx];
- element_idx != y_bus.y_bus_entry_indptr()[data_idx + 1]; ++element_idx) {
- Idx const obj = y_bus.y_bus_element()[element_idx].idx;
- YBusElementType const type = y_bus.y_bus_element()[element_idx].element_type;
- // shunt
- if (type == YBusElementType::shunt) {
- if (measured_value.has_shunt(obj)) {
- // G += Ys^H * Ys / variance
- block.g() += dot(hermitian_transpose(param.shunt_param[obj]), param.shunt_param[obj]) /
- measured_value.shunt_power(obj).variance;
- }
+ for (Idx row = 0; row != n_bus_; ++row) {
+ for (Idx data_idx_lu = row_indptr[row]; data_idx_lu != row_indptr[row + 1]; ++data_idx_lu) {
+ Idx const col = col_indices[data_idx_lu];
+ // get a reference and reset block to zero
+ SEGainBlock& block = data_gain_[data_idx_lu];
+ block = SEGainBlock{};
+ // get data idx of y bus,
+ // skip for a fill-in
+ Idx const data_idx = y_bus.map_lu_y_bus()[data_idx_lu];
+ if (data_idx == -1) {
+ continue;
}
- // branch
- else {
- // branch from- and to-side index at 0, and 1 position
- IntS const b0 = static_cast(type) / 2;
- IntS const b1 = static_cast(type) % 2;
- // measured at from-side: 0, to-side: 1
- for (IntS const measured_side : std::array{0, 1}) {
- // has measurement
- if (std::invoke(has_branch_[measured_side], measured_value, obj)) {
- // G += Y{side, b0}^H * Y{side, b1} / variance
- block.g() += dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b0]),
- param.branch_param[obj].value[measured_side * 2 + b1]) /
- std::invoke(branch_power_[measured_side], measured_value, obj).variance;
+ // fill block with voltage measurement, only diagonal
+ if ((row == col) && measured_value.has_voltage(row)) {
+ // G += 1.0 / variance
+ // for 3x3 tensor, fill diagonal
+ block.g() += ComplexTensor{1.0 / measured_value.voltage_var(row)};
+ }
+ // fill block with branch, shunt measurement
+ for (Idx element_idx = y_bus.y_bus_entry_indptr()[data_idx];
+ element_idx != y_bus.y_bus_entry_indptr()[data_idx + 1]; ++element_idx) {
+ Idx const obj = y_bus.y_bus_element()[element_idx].idx;
+ YBusElementType const type = y_bus.y_bus_element()[element_idx].element_type;
+ // shunt
+ if (type == YBusElementType::shunt) {
+ if (measured_value.has_shunt(obj)) {
+ // G += Ys^H * Ys / variance
+ block.g() += dot(hermitian_transpose(param.shunt_param[obj]), param.shunt_param[obj]) /
+ measured_value.shunt_power(obj).variance;
+ }
+ }
+ // branch
+ else {
+ // branch from- and to-side index at 0, and 1 position
+ IntS const b0 = static_cast(type) / 2;
+ IntS const b1 = static_cast(type) % 2;
+ // measured at from-side: 0, to-side: 1
+ for (IntS const measured_side : std::array{0, 1}) {
+ // has measurement
+ if (std::invoke(has_branch_[measured_side], measured_value, obj)) {
+ // G += Y{side, b0}^H * Y{side, b1} / variance
+ block.g() +=
+ dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b0]),
+ param.branch_param[obj].value[measured_side * 2 + b1]) /
+ std::invoke(branch_power_[measured_side], measured_value, obj).variance;
+ }
}
}
}
- }
- // fill block with injection measurement
- // injection measurement exist
- if (measured_value.has_bus_injection(row)) {
- // Q_ij = Y_bus_ij
- block.q() = y_bus.admittance()[data_idx];
- // R_ii = -variance, only diagonal
- if (row == col) {
- // assign variance to diagonal of 3x3 tensor, for asym
- block.r() = ComplexTensor{-measured_value.bus_injection_power(row).variance};
+ // fill block with injection measurement
+ // injection measurement exist
+ if (measured_value.has_bus_injection(row)) {
+ // Q_ij = Y_bus_ij
+ block.q() = y_bus.admittance()[data_idx];
+ // R_ii = -variance, only diagonal
+ if (row == col) {
+ // assign variance to diagonal of 3x3 tensor, for asym
+ block.r() = ComplexTensor{-measured_value.bus_injection_power(row).variance};
+ }
}
- }
- // injection measurement not exist
- else {
- // Q_ij = 0
- // R_ii = -1.0, only diagonal
- // assign -1.0 to diagonal of 3x3 tensor, for asym
- if (row == col) {
- block.r() = ComplexTensor{-1.0};
+ // injection measurement not exist
+ else {
+ // Q_ij = 0
+ // R_ii = -1.0, only diagonal
+ // assign -1.0 to diagonal of 3x3 tensor, for asym
+ if (row == col) {
+ block.r() = ComplexTensor{-1.0};
+ }
}
}
}
// loop all transpose entry for QH
// assign the hermitian transpose of the transpose entry of Q
- for (Idx data_idx = 0; data_idx != y_bus.nnz(); ++data_idx) {
- Idx const data_idx_tranpose = y_bus.transpose_entry()[data_idx];
- data_gain_[data_idx].qh() = hermitian_transpose(data_gain_[data_idx_tranpose].q());
+ for (Idx data_idx_lu = 0; data_idx_lu != y_bus.nnz_lu(); ++data_idx_lu) {
+ // skip for fill-in
+ if (y_bus.map_lu_y_bus()[data_idx_lu] == -1) {
+ continue;
+ }
+ Idx const data_idx_tranpose = y_bus.lu_transpose_entry()[data_idx_lu];
+ data_gain_[data_idx_lu].qh() = hermitian_transpose(data_gain_[data_idx_tranpose].q());
}
// prefactorize
- bsr_solver_.prefactorize(data_gain_.data());
+ sparse_solver_.prefactorize(data_gain_, perm_);
}
void prepare_rhs(YBus const& y_bus, MeasuredValues const& measured_value,
@@ -740,12 +762,12 @@ class IterativeLinearSESolver {
for (Idx bus = 0; bus != n_bus_; ++bus) {
Idx const data_idx = y_bus.bus_entry()[bus];
// reset rhs block to fill values
- SERhs& rhs_block = rhs_[bus];
- rhs_block = {};
+ SERhs& rhs_block = x_rhs_[bus];
+ rhs_block = SERhs{};
// fill block with voltage measurement
if (measured_value.has_voltage(bus)) {
// eta += u / variance
- rhs_block.eta += u[bus] / measured_value.voltage_var(bus);
+ rhs_block.eta() += u[bus] / measured_value.voltage_var(bus);
}
// fill block with branch, shunt measurement, need to convert to current
for (Idx element_idx = y_bus.y_bus_entry_indptr()[data_idx];
@@ -757,7 +779,7 @@ class IterativeLinearSESolver {
if (measured_value.has_shunt(obj)) {
SensorCalcParam const& m = measured_value.shunt_power(obj);
// eta -= Ys^H * i_shunt / variance
- rhs_block.eta -=
+ rhs_block.eta() -=
dot(hermitian_transpose(param.shunt_param[obj]), conj(m.value / u[bus])) / m.variance;
}
}
@@ -776,7 +798,7 @@ class IterativeLinearSESolver {
// NOTE: not the current bus!
Idx const measured_bus = branch_bus_idx[obj][measured_side];
// eta += Y{side, b}^H * i_branch_{f, t} / variance
- rhs_block.eta +=
+ rhs_block.eta() +=
dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b]),
conj(m.value / u[measured_bus])) /
m.variance;
@@ -786,7 +808,7 @@ class IterativeLinearSESolver {
}
// fill block with injection measurement, need to convert to current
if (measured_value.has_bus_injection(bus)) {
- rhs_block.tau = conj(measured_value.bus_injection_power(bus).value / u[bus]);
+ rhs_block.tau() = conj(measured_value.bus_injection_power(bus).value / u[bus]);
}
}
}
@@ -800,16 +822,16 @@ class IterativeLinearSESolver {
return 1.0;
}
if constexpr (sym) {
- return cabs(x_[math_topo_->slack_bus_].u) / x_[math_topo_->slack_bus_].u;
+ return cabs(x_rhs_[math_topo_->slack_bus_].u()) / x_rhs_[math_topo_->slack_bus_].u();
}
else {
- return cabs(x_[math_topo_->slack_bus_].u(0)) / x_[math_topo_->slack_bus_].u(0);
+ return cabs(x_rhs_[math_topo_->slack_bus_].u()(0)) / x_rhs_[math_topo_->slack_bus_].u()(0);
}
}();
for (Idx bus = 0; bus != n_bus_; ++bus) {
// phase offset to calculated voltage as normalized
- ComplexValue u_normalized = x_[bus].u * angle_offset;
+ ComplexValue u_normalized = x_rhs_[bus].u() * angle_offset;
// get dev of last iteration, get max
double const dev = max_val(cabs(u_normalized - u[bus]));
max_dev = std::max(dev, max_dev);
diff --git a/include/power_grid_model/math_solver/linear_pf_solver.hpp b/include/power_grid_model/math_solver/linear_pf_solver.hpp
index 2531105815..7905b09997 100644
--- a/include/power_grid_model/math_solver/linear_pf_solver.hpp
+++ b/include/power_grid_model/math_solver/linear_pf_solver.hpp
@@ -35,7 +35,7 @@ if there are sources
#include "../power_grid_model.hpp"
#include "../three_phase_tensor.hpp"
#include "../timer.hpp"
-#include "bsr_solver.hpp"
+#include "sparse_lu_solver.hpp"
#include "y_bus.hpp"
namespace power_grid_model {
@@ -51,16 +51,16 @@ class LinearPFSolver {
: n_bus_{y_bus.size()},
load_gen_bus_indptr_{topo_ptr, &topo_ptr->load_gen_bus_indptr},
source_bus_indptr_{topo_ptr, &topo_ptr->source_bus_indptr},
- mat_data_(y_bus.nnz()),
- rhs_(n_bus_),
- bsr_solver_{y_bus.size(), bsr_block_size_, y_bus.shared_indptr(), y_bus.shared_indices()} {
+ mat_data_(y_bus.nnz_lu()),
+ sparse_solver_{y_bus.shared_indptr_lu(), y_bus.shared_indices_lu(), y_bus.shared_diag_lu()},
+ perm_(n_bus_) {
}
MathOutput run_power_flow(YBus const& y_bus, PowerFlowInput const& input,
CalculationInfo& calculation_info) {
// getter
ComplexTensorVector const& ydata = y_bus.admittance();
- IdxVector const& bus_entry = y_bus.bus_entry();
+ IdxVector const& bus_entry = y_bus.lu_diag();
// output
MathOutput output;
output.u.resize(n_bus_);
@@ -71,9 +71,14 @@ class LinearPFSolver {
Timer sub_timer(calculation_info, 2221, "Prepare matrix");
// copy y bus data
- std::copy(ydata.begin(), ydata.end(), mat_data_.begin());
- // set rhs to zero
- std::fill(rhs_.begin(), rhs_.end(), ComplexValue{0.0});
+ std::transform(y_bus.map_lu_y_bus().cbegin(), y_bus.map_lu_y_bus().cend(), mat_data_.begin(), [&](Idx k) {
+ if (k == -1) {
+ return ComplexTensor{};
+ }
+ else {
+ return ydata[k];
+ }
+ });
// loop to all loads and sources, j as load number
IdxVector const& load_gen_bus_idxptr = *load_gen_bus_indptr_;
@@ -92,15 +97,15 @@ class LinearPFSolver {
// YBus_diag += Y_source
mat_data_[data_sequence] += y_bus.math_model_param().source_param[source_number];
// rhs += Y_source_j * U_ref_j
- rhs_[bus_number] += dot(y_bus.math_model_param().source_param[source_number],
- ComplexValue{input.source[source_number]});
+ output.u[bus_number] += dot(y_bus.math_model_param().source_param[source_number],
+ ComplexValue{input.source[source_number]});
}
}
// solve
// u vector will have I_injection for slack bus for now
sub_timer = Timer(calculation_info, 2222, "Solve sparse linear equation");
- bsr_solver_.solve(mat_data_.data(), rhs_.data(), output.u.data());
+ sparse_solver_.prefactorize_and_solve(mat_data_, perm_, output.u, output.u);
// calculate math result
sub_timer = Timer(calculation_info, 2223, "Calculate Math Result");
@@ -117,9 +122,9 @@ class LinearPFSolver {
std::shared_ptr source_bus_indptr_;
// sparse linear equation
ComplexTensorVector mat_data_;
- ComplexValueVector rhs_;
// sparse solver
- BSRSolver bsr_solver_;
+ SparseLUSolver, ComplexValue, ComplexValue> sparse_solver_;
+ typename SparseLUSolver, ComplexValue, ComplexValue>::BlockPermArray perm_;
void calculate_result(YBus const& y_bus, PowerFlowInput const& input, MathOutput& output) {
// call y bus
diff --git a/include/power_grid_model/math_solver/math_solver.hpp b/include/power_grid_model/math_solver/math_solver.hpp
index 0130a3c48c..5646abb2bf 100644
--- a/include/power_grid_model/math_solver/math_solver.hpp
+++ b/include/power_grid_model/math_solver/math_solver.hpp
@@ -27,9 +27,10 @@ template
class MathSolver {
public:
MathSolver(std::shared_ptr const& topo_ptr,
- std::shared_ptr const> const& param)
+ std::shared_ptr const> const& param,
+ std::shared_ptr const& y_bus_struct = {})
: topo_ptr_{topo_ptr},
- y_bus_{topo_ptr, param},
+ y_bus_{topo_ptr, param, y_bus_struct},
all_const_y_{std::all_of(topo_ptr->load_gen_type.cbegin(), topo_ptr->load_gen_type.cend(), [](LoadGenType x) {
return x == LoadGenType::const_y;
})} {
@@ -94,12 +95,17 @@ class MathSolver {
newton_pf_solver_.reset();
linear_pf_solver_.reset();
iterative_current_pf_solver_.reset();
+ iterative_linear_se_solver_.reset();
}
void update_value(std::shared_ptr const> const& math_model_param) {
y_bus_.update_admittance(math_model_param);
}
+ std::shared_ptr shared_y_bus_struct() const {
+ return y_bus_.shared_y_bus_struct();
+ }
+
private:
std::shared_ptr topo_ptr_;
YBus y_bus_;
diff --git a/include/power_grid_model/math_solver/mkl_pardiso_solver.hpp b/include/power_grid_model/math_solver/mkl_pardiso_solver.hpp
deleted file mode 100644
index 08085fade5..0000000000
--- a/include/power_grid_model/math_solver/mkl_pardiso_solver.hpp
+++ /dev/null
@@ -1,390 +0,0 @@
-// SPDX-FileCopyrightText: 2022 Contributors to the Power Grid Model project
-//
-// SPDX-License-Identifier: MPL-2.0
-
-#pragma once
-#ifndef POWER_GRID_MODEL_MATH_SOLVER_MKL_PARDISO_SOLVER_HPP
-#define POWER_GRID_MODEL_MATH_SOLVER_MKL_PARDISO_SOLVER_HPP
-// BSR solver using MKL PARDISO
-// https://software.intel.com/en-us/mkl-developer-reference-c-pardiso
-
-#ifdef POWER_GRID_MODEL_USE_MKL
-
-#include "../exception.hpp"
-#include "../power_grid_model.hpp"
-
-#ifdef POWER_GRID_MODEL_USE_MKL_AT_RUNTIME
-
-// add iostream header to report loading mkl status
-#include
-#include