Merge ~bullwinkle-team/ubuntu/+source/rocalution:bullwinkle/llvm-21/ubuntu/devel into ubuntu/+source/rocalution:ubuntu/devel

[Bruno Bernardo de Moura]

Update to new upstream version 7.1.0

Link for build reference: https://launchpad.net/~bruno-bdmoura/+archive/ubuntu/lp-2138877

[Bruno Bernardo de Moura (bruno-bdmoura)]

Uploaded package to this ppa: https://launchpad.net/~bruno-bdmoura/+archive/ubuntu/lp-2138877

(-proposed and all archs enabled)

[Talha Can Havadar (tchavadar)]

Hello Bruno, I see we have segfaults for arm64 in the ppa you shared, do you happen to know the reason?

[Bruno Bernardo de Moura (bruno-bdmoura)]

Unfortunately I do not know the reason why this happens. However, the test that's causing the segmentation fault is the same one who failed whenever I was mistakenly trying to build this package using clang instead of the default upstream compiler.

I have asked Cory on matrix regarding this topic, let's see what he has to say.

[Talha Can Havadar (tchavadar)]

Can you point me the source code where it fails? is there like an obvious null pointer dereference some where?

[Bruno Bernardo de Moura (bruno-bdmoura)]

The error occurs on this test function: https://git.launchpad.net/~bullwinkle-team/ubuntu/+source/rocalution/tree/clients/tests/test_local_vector.cpp?h=bullwinkle/llvm-21/ubuntu/devel#n439.

This snippet then calls the testing_extract_coarse_mapping function, located here: https://git.launchpad.net/~bullwinkle-team/ubuntu/+source/rocalution/tree/clients/include/testing_local_vector.hpp?h=bullwinkle/llvm-21/ubuntu/devel#n1403, which is where the test fails.

The method invoked by the custom vector can be seen here: https://git.launchpad.net/~bullwinkle-team/ubuntu/+source/rocalution/tree/src/base/local_vector.cpp?h=bullwinkle/llvm-21/ubuntu/devel#n1210

I honestly don't know exactly how different compilers treat c++ template usage, but I would guess that's where the error might be coming from between the two platforms.

Besides that I don't see any obvious errors in that function.

[Bruno Bernardo de Moura (bruno-bdmoura)]

And also, I just have to say I'm by no means a c++ expert, so there might still be fact some intricate details on the above implementations that I'm missing.

[Bruno Bernardo de Moura (bruno-bdmoura)]

Talha, just as an update, I got access to the hinyari machine, and even though I had to skip the symbols file verification using the -c0 flag, the package builds and runs the tests without breaking, even the test that's breaking on launchpad.

[Talha Can Havadar (tchavadar)]

ok that is suspicious can you check the what is the revision of arm in hinyari? maybe it is different than the builders we have in launchpad and maybe the extensions supported by the arm cpu in hinyari so we can ask launchpad team

[Talha Can Havadar (tchavadar)]

it is strange but maybe some vector instructions are not there for builders we have in LP best to be sure

[Talha Can Havadar (tchavadar)]

That said, I think we need more information about the failure, I see some log functions called on the entrance of each function but cant see the output anywhere in build log maybe there is a way to get them to see at which step it fails.

[Andreas Hasenack (ahasenack)]

Don't delete the previous changelog entry from resolute:
diff --git a/debian/changelog b/debian/changelog
index 2a0b776..00fd507 100644
--- a/debian/changelog
+++ b/debian/changelog
@@ -1,8 +1,32 @@
-rocalution (5.7.1-3build1) resolute; urgency=medium
-
- * No-change mass rebuild for Ubuntu 26.04 (LP: #2132257)
-
- -- Sebastien Bacher <email address hidden> Thu, 29 Jan 2026 21:42:09 +0100
+rocalution (7.1.0-0ubuntu2) resolute; urgency=medium
+
+ * d/rules: skip failing local_vector_test.extract_coarse_mapping_int test
+
+ -- Bruno Bernardo de Moura <email address hidden> Tue, 03 Feb 2026 16:40:38 -0300
+
+rocalution (7.1.0-0ubuntu1) resolute; urgency=medium
+
+ * New upstream release 7.1.0.
+ * d/p/series: drop outdated and unused patches for 7.1.0
+ - drop d/p/spelling.patch
+ - drop d/p/use-readthedocs-theme.patch
+ - drop d/p/docs-table-of-contents.patch
+ - drop d/p/use-local-mathjax.patch.
+ - drop d/p/fix-host-and-backend-function-documentation.patch
+ * d/control: update standards version to 4.7.2
+ * d/control: update maintainer field for ubuntu archive
+ * d/control: update build-depends for rocm stack
+ * d/control: update SOVERSION (0 -> 1) after API changes
+ * d/rules: fix docs build step with rocm-docs-build
+ * d/rules: fix lintian warnings and improve docs cleanup
+ * d/{librocalution-doc.doc-base,librocalution-doc.docs}:
+ fix doc generation paths
+ * d/shlibs: drop auto generated file
+ * d/{t/control,librocalution0-tests.install}:
+ update test files SOVERSION (0 -> 1)
+ * d/librocalution1.symbols: add symbols file for 7.1.0
+
+ -- Bruno Bernardo de Moura <email address hidden> Wed, 10 Dec 2025 14:27:06 -0300

 rocalution (5.7.1-3) unstable; urgency=medium

review: Needs Fixing

[Andreas Hasenack (ahasenack)]

+ * d/shlibs: drop auto generated file

...
diff --git a/debian/librocalution1.symbols b/debian/librocalution1.symbols
new file mode 100644
index 0000000..2f69eb9
--- /dev/null
+++ b/debian/librocalution1.symbols
@@ -0,0 +1,12060 @@
+librocalution.so.1 librocalution1 #MINVER#
+* Build-Depends-Package: librocalution-dev
+ _ZN10rocalution10BaseMatrixISt7complexIdEE10ItLAnalyseEb@Base 7.1.0
+ _ZN10rocalution10BaseMatrixISt7complexIdEE10ItUAnalyseEb@Base 7.1.0
+ _ZN10rocalution10BaseMatrixISt7complexIdEE10MatMatMultERKS3_S5_@Base 7.1.0
...

These are c++ symbols, and should be demangled with c++filt if you want to go down this route. Note that c++ symbols are notoriously difficult to manage, and that's why the d/shlibs file was there in the first place. Are you sure you want to make these changes?

[Andreas Hasenack (ahasenack)]

d/rues:
@@ -17,16 +17,24 @@ override_dh_auto_configure-arch:

 override_dh_auto_test-arch:
 ifeq (,$(filter nocheck,$(DEB_BUILD_OPTIONS)))
- obj-$(DEB_HOST_GNU_TYPE)/clients/staging/rocalution-test
+ # As described in (https://github.com/ROCm/rocALUTION/issues/322), the test
+ # `local_vector_test.extract_coarse_mapping_int` causes a segmentation fault on
+ # ARM systems. It consistently triggers a SIGSEGV in the extract_coarse_mapping_int
+ # function, even though the same test runs successfully on AMD64 platforms.
+ # Therefore, it is filtered here to enable successful package builds.
+ obj-$(DEB_HOST_GNU_TYPE)/clients/staging/rocalution-test \
+ --gtest_filter=-local_vector_test.extract_coarse_mapping_int

A test is supposedly correctly showing a problem, and you are covering it. Just saying that the test passes on amd64 is not a reason: architectures are indeed different, and just because it passes on one architecture doesn't immediately means it's ok to fail in another.

I see you filed a bug upstream, and there is no response yet. I think this needs a bit more digging, or I didn't quite understand why it's ok to build and release the package on arm64 with a segfault in one of the tests.

review: Needs Information

[Igor Luppi (igorluppi)]

Oh after years someone rebuild v5 of rocalution just this week hehe (https://git.launchpad.net/ubuntu/+source/rocalution/commit/?id=9ae3859003b9b41bff59f3eb8182030193f7ce8a)

It was not a deletion, just a merge conflict because ubuntu/devel diverged after the MP was open, what a bad luck

[Igor Luppi (igorluppi)]

Rebased, now it should be fine :)

[Bruno Bernardo de Moura (bruno-bdmoura)]

@Andreas, thanks for the comments!

Now adressing each of them:

1 - The changelog skip was due to the merge conflict with the rebuild of v5, so it's now already fixed

2- Related to the symbols file: you are absolutely correct, I simply added it due to mistakenly deleting the d/shlibs file after not updating the SONAME version of the package on it.

I've just tested locally and by simply dropping the symbols file and updating the SONAME version on d/shlibs lintian gives no warnings related to symbols. This way, the symbols file will be removed, and the d/shlibs updated accordingly

3- Finally, related to the crashing test on ARM: the AMD team agreed to skip that test due to it not being a mandatory package for the stack itself.

Also, if needed, instead of being a filter on the debian/rules file, the skipping of the test could also be performed as a patch that applies to its file and only runs that specific crashing function in case the architecture is x86 (as is currently performed in hipblaslt, for example). Although this other approach would at least guarantee that the test is passing on x86, the intent of it would be practically the same.

[Andreas Hasenack (ahasenack)]

3) Yes, skipping the test only on arm is less bad for now. But please change the comments, they make it sound like the test is invalid because it passes on amd64. Different architectures behave differently.

[Andreas Hasenack (ahasenack)]

1) ok

2) waiting the the change

[Bruno Bernardo de Moura (bruno-bdmoura)]

@Andreas, back to the topics:

2 - Change done
3 - Transformed it to a patch and description enhanced as dep3 header

[Andreas Hasenack (ahasenack)]

--- rocalution.orig/clients/tests/test_local_vector.cpp
+++ rocalution/clients/tests/test_local_vector.cpp
@@ -435,6 +435,7 @@ GENERATE_P_TEST_CASES(parameterized_loca
                       testing_set_continuous_values,
                       set_continuous_values)

+#if defined(__x86_64__) || defined(__i386__)
 // Test for LocalVector::ExtractCoarseMapping
 TEST_F(local_vector_test, extract_coarse_mapping_int)
 {
@@ -444,6 +445,7 @@ TEST_F(local_vector_test, extract_coarse
     testing_extract_coarse_mapping<int>();
     stop_rocalution();
 }
+#endif

 GENERATE_P_TEST_CASES(parameterized_local_vector_test,
                       testing_move_to_host_async,

I was expecting excluding the test via d/rules, by checking the arch there and calling the test command with or without the gtest_filter argument.

A patch works too, of course, so this part I don't mind that much, keeping in mind the patch might not apply anymore if the surrounding code changes.

But instead of excluding the arch where the test fails, you are listing two arches where it should run (and we don't even build for i386). At the same time we have "Arch any" in d/control, so the above patch will exclude the test in ppc64el, s390x, armhf, arm64, riscv64, thus reducing coverage. Or was the test also not passing there?

[Bruno Bernardo de Moura (bruno-bdmoura)]

@Andreas,

Updated the rules file to skip that test on arm64 instead of applying a patch.

Additionally, the builds can be seen on this ppa: https://launchpad.net/~bruno-bdmoura/+archive/ubuntu/rocm-llvm21-all-archs/+packages

[Andreas Hasenack (ahasenack)]

I see armhf failed to build, but that is already the case with the current package in the archive[1], so no regression.

+# On ARM64, the test local_vector_test.extract_coarse_mapping_int causes
+# segmentation fault, triggering a SIGSEGV in the extract_coarse_mapping
+# function. Since different compilers handle memory allocation differently
+# from one another, the test is being skipped to allow for successfull builds
+# on those ARM plaforms.
+ifeq ($(DEB_HOST_ARCH),arm64)
+ GTEST_FILTER_VAL := *:-local_vector_test.extract_coarse_mapping_int
+endif

Please include the link you had before to the upstream issue about this: https://github.com/ROCm/rocALUTION/issues/322

1. https://launchpad.net/ubuntu/+source/rocalution/5.7.1-3build1

review: Needs Fixing

[Andreas Hasenack (ahasenack)]

This looks wrong:

-override_dh_auto_build-indep:
+execute_after_dh_auto_build-indep: export http_proxy=127.0.0.1:9
+execute_after_dh_auto_build-indep: export https_proxy=127.0.0.1:9
+execute_after_dh_auto_build-indep:
 ifeq (,$(filter nodoc,$(DEB_BUILD_OPTIONS)))
- cd docs/.doxygen; doxygen
- sphinx-build -b html docs html
+ perl -pi -e 's/WARN_AS_ERROR.*=.*YES/WARN_AS_ERROR = NO/' docs/doxygen/Doxyfile
+ rocm-docs-build
 endif

End result in d/rules:

execute_after_dh_auto_build-indep: export http_proxy=127.0.0.1:9
execute_after_dh_auto_build-indep: export https_proxy=127.0.0.1:9
execute_after_dh_auto_build-indep:
ifeq (,$(filter nodoc,$(DEB_BUILD_OPTIONS)))
    perl -pi -e 's/WARN_AS_ERROR.*=.*YES/WARN_AS_ERROR = NO/' docs/doxygen/Doxyfile
    rocm-docs-build
endif

Perhaps take *inspiration* on the rocrand d/rules section. Other rocm packages also had to address this:

# python3-rocm-docs can only build offline with the help of some envvars
override_dh_auto_build-indep:
ifeq (,$(filter nodoc,$(DEB_BUILD_OPTIONS)))
    http_proxy='127.0.0.1:9' \
    https_proxy='127.0.0.1:9' \
    rocm-docs-build -r $(shell dpkg-parsechangelog -S Version | sed 's/[+-].*//')
    rm -rf build/html/_static/fonts
    rm -rf build/html/_static/vendor
endif

Note the "\" continuation marks at the end of each line, and how the proxy variable is being set in the same shell command. Remember that each line is a shell invocation on its own, so if you set a variable in one line, without the "\" continuation, it won't exist in the next line, which is a new shell invocation.

review: Needs Fixing

[Bruno Bernardo de Moura (bruno-bdmoura)]

@Andreas, thanks for the explanation!

The `execute_after_dh_auto_build-indep` rule is now fixed as suggested.

[Andreas Hasenack (ahasenack)]

+1 thanks

review: Approve

[Andreas Hasenack (ahasenack)]

Sponsored:

Uploading rocalution_7.1.0-0ubuntu2.dsc
Uploading rocalution_7.1.0.orig.tar.gz
Uploading rocalution_7.1.0-0ubuntu2.debian.tar.xz
Uploading rocalution_7.1.0-0ubuntu2_source.buildinfo
Uploading rocalution_7.1.0-0ubuntu2_source.changes

This has binary NEW packages, due to the soname change, so it will require an archive admin approval before the binaries land in resolute-proposed.

diff --git a/.azuredevops/rocm-ci.yml b/.azuredevops/rocm-ci.yml
new file mode 100644
index 0000000..b615df9
--- /dev/null
+++ b/.azuredevops/rocm-ci.yml
@@ -0,0 +1,44 @@
+resources:
+  repositories:
+  - repository: pipelines_repo
+    type: github
+    endpoint: ROCm
+    name: ROCm/ROCm
+
+variables:
+- group: common
+- template: /.azuredevops/variables-global.yml@pipelines_repo
+
+trigger:
+  batch: true
+  branches:
+    include:
+    - develop
+    - mainline
+  paths:
+    exclude:
+    - .githooks
+    - .github
+    - .jenkins
+    - docs
+    - '.*.y*ml'
+    - '*.md'
+
+pr:
+  autoCancel: true
+  branches:
+    include:
+    - develop
+    - mainline
+  paths:
+    exclude:
+    - .githooks
+    - .github
+    - .jenkins
+    - docs
+    - '.*.y*ml'
+    - '*.md'
+  drafts: false
+
+jobs:
+  - template: ${{ variables.CI_COMPONENT_PATH }}/rocALUTION.yml@pipelines_repo
diff --git a/.githooks/pre-commit b/.githooks/pre-commit
index 04dabeb..7d23219 100755
--- a/.githooks/pre-commit
+++ b/.githooks/pre-commit
@@ -35,7 +35,7 @@ fi
 for file in $files; do
     if [[ -e $file ]]; then
         /usr/bin/perl -pi -e 'INIT { exit 1 if !-f $ARGV[0] || -B $ARGV[0]; $year = (localtime)[5] + 1900 }
-            s/^([*\/#[:space:]]*)Copyright\s+(?:\(C\)\s*)?(\d+)(?:\s*-\s*\d+)?/qq($1Copyright (c) $2@{[$year != $2 ? "-$year" : ""]})/ie
+            s/^([*\/#[:space:]]*)Copyright\s+(?:\(C\)\s*)?(\d+)(?:\s*-\s*\d+)?/qq($1Copyright (C) $2@{[$year != $2 ? "-$year" : ""]})/ie
             if $. < 10' "$file" && git add -u "$file"
     fi
 done
diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS
old mode 100644
new mode 100755
index 4b7a84c..300b638
--- a/.github/CODEOWNERS
+++ b/.github/CODEOWNERS
@@ -1 +1,6 @@
-* @ntrost57 @YvanMokwinski @jsandham
+* @ntrost57 @YvanMokwinski @jsandham @kliegeois
+# Documentation files
+docs/* @ROCm/rocm-documentation
+*.md @ROCm/rocm-documentation
+*.rst @ROCm/rocm-documentation
+.readthedocs.yaml @ROCm/rocm-documentation
diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md
new file mode 100644
index 0000000..2ab06f0
--- /dev/null
+++ b/.github/CONTRIBUTING.md
@@ -0,0 +1,160 @@
+<head>
+  <meta charset="UTF-8">
+  <meta name="description" content="Contributing to rocALUTION">
+  <meta name="keywords" content="ROCm, contributing, rocALUTION">
+</head>
+
+# Contributing to rocALUTION #
+
+AMD welcomes contributions to rocALUTION from the community. Whether those contributions are bug reports, bug fixes, documentation additions, performance notes, or other improvements, we value collaboration with our users. We can build better solutions together. Please follow these details to help ensure your contributions will be successfully accepted.
+
+Our code contriubtion guidelines closely follow the model of [GitHub pull-requests](https://help.github.com/articles/using-pull-requests/).  This repository follows the [git flow](http://nvie.com/posts/a-successful-git-branching-model/) workflow, which dictates a /master branch where releases are cut, and a /develop branch which serves as an integration branch for new code.
+
+## Issue Discussion ##
+
+Please use the GitHub Issues tab to notify us of issues.
+
+* Use your best judgement for issue creation. If your issue is already listed, upvote the issue and
+  comment or post to provide additional details, such as how you reproduced this issue.
+* If you're not sure if your issue is the same, err on the side of caution and file your issue.
+  You can add a comment to include the issue number (and link) for the similar issue. If we evaluate
+  your issue as being the same as the existing issue, we'll close the duplicate.
+* If your issue doesn't exist, use the issue template to file a new issue.
+  * When filing an issue, be sure to provide as much information as possible, including script output so
+    we can collect information about your configuration. This helps reduce the time required to
+    reproduce your issue.
+  * Check your issue regularly, as we may require additional information to successfully reproduce the
+    issue.
+* You may also open an issue to ask questions to the maintainers about whether a proposed change
+  meets the acceptance criteria, or to discuss an idea pertaining to the library.
+
+## Acceptance Criteria ##
+
+rocALUTION is a sparse linear algebra library with focus on exploring fine-grained parallelism on top of the AMD ROCm runtime and toolchains, targeting modern CPU and GPU platforms. Based on C++ and HIP, it provides a portable, generic and flexible design that allows seamless integration with other scientific software packages.
+
+In rocALUTION we are interested in contributions that:
+* Fix bugs, improve documentation, enhance testing, reduce complexity.
+* Improve the performance of existing routines.
+* Add missing functionality such as new multigrid solvers, iterative solvers, direct solvers, or preconditioners.
+* Extending new or existing functionality to work with MPI or accelerators (such as GPU devices).
+
+We encourage contributors to leverage the GitHub "Issues" tab to discuss possible additions they would like to add.
+
+### Exceptions ###
+
+rocALUTION places a heavy emphasis on being high performance. Because of this, contributions that add new routines (or that modify existing routines) must do so from the perspective that they offer high performance in relation to the hardware they are run on. Furthermore, all routines added to rocalution must have at a minimum a host solution as all routines must have the ability to fall back to a host solution if a GPU accelerator is not avaiable. Because compile times, binary sizes, and general library complexity are important considerations, we reserve the right to make decisions on whether a proposed routine is too niche or specialized to be worth including.
+
+## Code Structure ##
+
+The following is the structure of the rocALUTION library in the GitHub repository. A more detailed description of the directory structure can be found in the rocALUTION [documentation](https://rocm.docs.amd.com/projects/rocALUTION/en/latest/design/orga.html).
+
+The `src/` directory contains the library source code. This is broken up into three sub-directories:
+* `src/base`
+* `src/solvers`
+* `src/utils`
+
+The `src/base` Contains source code related to rocALUTION's vector, matrix, and stencil operator types as well as classes related to parallel management. This directory is further broken up into:
+* `src/base/hip` Contains HIP implementations of vector, matrix, and stencil operators.
+* `src/base/host` Contains host implementations of vector, matrix, and stencil operators.
+
+The `src/solvers` directory contains all the source code related to direct (`src/solvers/direct`), krylov (`src/solvers/krylov`), and multigrid solvers `src/solvers/multigrid`.
+
+The `src/utils` directory contains source code related to logging, memory allocation, math and timing functions.
+
+The `clients/` directory contains the testing and benchmarking code as well as all the samples demonstrating rocALUTION usage.
+
+The `docs/` directory contains all of the documentation files.
+
+## Coding Style ##
+
+In general, follow the style of the surrounding code. C and C++ code is formatted using `clang-format`. Use the clang-format version installed with ROCm (found in the `/opt/rocm/llvm/bin` directory). Please do not use your system's built-in `clang-format`, as this is a different version that may result in incorrect results.
+
+To format a file, use:
+
+```
+/opt/rocm/llvm/bin/clang-format -style=file -i <path-to-source-file>
+```
+
+To format all files, run the following script in rocALUTION directory:
+
+```
+#!/bin/bash
+git ls-files -z *.cc *.cpp *.h *.hpp *.cl *.h.in *.hpp.in *.cpp.in | xargs -0 /opt/rocm/llvm/bin/clang-format  -style=file -i
+```
+
+Also, githooks can be installed to format the code per-commit:
+
+```
+./.githooks/install
+```
+
+## Pull Request Guidelines ##
+
+When you create a pull request, you should target the default branch. Our current default branch is the **develop** branch, which serves as our integration branch.
+
+By submitting a pull request, you acknowlege and agree with the CLA below:
+
+Contribution License Agreement
+1. The code I am contributing is mine, and I have the right to license it.
+2. By submitting a pull request for this project I am granting you a license to distribute said code under the MIT License for the project.
+
+### Deliverables ###
+
+When raising a PR in rocALUTION here are some important things to include:
+
+1. For each new file in the repository, Please include the licensing header
+```
+/* ************************************************************************
+* Copyright (C) 20xx Advanced Micro Devices, Inc. All rights Reserved.
+*
+* Permission is hereby granted, free of charge, to any person obtaining a copy
+* of this software and associated documentation files (the "Software"), to deal
+* in the Software without restriction, including without limitation the rights
+* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+* copies of the Software, and to permit persons to whom the Software is
+* furnished to do so, subject to the following conditions:
+*
+* The above copyright notice and this permission notice shall be included in
+* all copies or substantial portions of the Software.
+*
+* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+* THE SOFTWARE.
+*
+* ************************************************************************ */
+```
+and adjust the date to the current year. When simply modifying a file, the date should automatically be updated pre-commit as long as the githook has been installed (./.githooks/install).
+
+2. When adding a new routine, please make sure you are also adding appropriate testing code. These new unit tests should integrate within the existing [googletest framework](https://github.com/google/googletest/blob/master/googletest/docs/primer.md). This typically involves adding the following files:
+
+* testing_<routine_name>.hpp file in the directory `clients/include/`
+* test_<routine_name>.cpp file in directory `clients/tests/`
+
+See existing tests for guidance when adding your own.
+
+3. When modifiying an existing routine, add appropriate testing to test_<routine_name>.cpp file in directory `clients/tests/`.
+
+4. Tests must have good code coverage.
+
+5. At a minimum, rocALUTION must have a host solution for each direct, iterative, multigrid, or preconditioner. If you add a accelerator solution (say using HIP targetting GPU devices) please also add a fall back host solution.
+
+6. Ensure code builds successfully. This includes making sure that the code can compile, that the code is properly formatted, and that all tests pass.
+
+7. Do not break existing test cases
+
+### Process ###
+
+When a PR is raised targetting the develop branch in rocALUTION, CI will be automatically triggered. This will:
+
+* Test that the PR passes static analysis (i.e ensure clang formatting rules have been followed).
+* Test that the documentation can be properly built
+* Ensure that the PR compiles on different OS and GPU device architecture combinations
+* Ensure that all tests pass on different OS and GPU device architecture combinations
+
+Feel free to ask questions on your PR regarding any CI failures you encounter.
+
+* Reviewers are listed in the CODEOWNERS file
diff --git a/.github/dependabot.yml b/.github/dependabot.yml
index 9cdf2d6..047fac7 100644
--- a/.github/dependabot.yml
+++ b/.github/dependabot.yml
@@ -6,7 +6,14 @@
 version: 2
 updates:
   - package-ecosystem: "pip" # See documentation for possible values
-    directory: "/docs/.sphinx" # Location of package manifests
+    directory: "/docs/sphinx" # Location of package manifests
     open-pull-requests-limit: 10
     schedule:
       interval: "daily"
+    target-branch: "develop"
+    labels:
+      - "documentation"
+      - "dependencies"
+      - "ci:docs-only"
+    reviewers:
+      - "samjwu"
diff --git a/.gitignore b/.gitignore
index 15ec500..50bd106 100644
--- a/.gitignore
+++ b/.gitignore
@@ -38,20 +38,9 @@ tags
 # build-in-source directory
 build
 
-# doc
-docBin
-_build
-
 # Visual Studio stuff
 *.vcxproj.user
 *.suo
 *.sdf
 *.pdb
 *.opensdf
-
-# documentation artifacts
-_build/
-_images/
-_static/
-_templates/
-_toc.yml
diff --git a/.jenkins/codecov.groovy b/.jenkins/codecov.groovy
index 175da9a..76162e8 100644
--- a/.jenkins/codecov.groovy
+++ b/.jenkins/codecov.groovy
@@ -18,7 +18,7 @@ def runCI =
     prj.paths.build_command = './install.sh -cg --codecoverage'
     prj.compiler.compiler_name = 'c++'
     prj.compiler.compiler_path = 'c++'
-    prj.libraryDependencies = ['rocPRIM', 'rocBLAS-internal', 'rocSPARSE-internal', 'rocRAND']
+    prj.libraryDependencies = ['rocPRIM', 'hipBLAS-common', 'hipBLASLt', 'rocBLAS', 'rocSPARSE', 'rocRAND']
     prj.defaults.ccache = false
 
     // Define test architectures, optional rocm version argument is available
@@ -52,12 +52,10 @@ def runCI =
 ci: {
     String urlJobName = auxiliary.getTopJobName(env.BUILD_URL)
 
-    def propertyList = ["compute-rocm-dkms-no-npi":[pipelineTriggers([cron('0 1 * * 6')])],
-                        "compute-rocm-dkms-no-npi-hipclang":[pipelineTriggers([cron('0 1 * * 6')])] ]
+    def propertyList = []
     propertyList = auxiliary.appendPropertyList(propertyList)
 
-    def jobNameList = ["compute-rocm-dkms-no-npi":([ubuntu18:['gfx900']]),
-                       "compute-rocm-dkms-no-npi-hipclang":([ubuntu18:['gfx900']])]
+    def jobNameList = []
     jobNameList = auxiliary.appendJobNameList(jobNameList)
 
     propertyList.each
@@ -81,7 +79,7 @@ ci: {
     {
         properties(auxiliary.addCommonProperties([pipelineTriggers([cron('0 1 * * *')])]))
         stage(urlJobName) {
-            runCI([ubuntu18:['gfx900']], urlJobName)
+            runCI([], urlJobName)
         }
     }
 }
diff --git a/.jenkins/common.groovy b/.jenkins/common.groovy
index 11e8586..1e973d1 100644
--- a/.jenkins/common.groovy
+++ b/.jenkins/common.groovy
@@ -12,7 +12,7 @@ def runCompileCommand(platform, project, boolean sameOrg=false)
     {
         project.libraryDependencies.each
         { libraryName ->
-            getDependenciesCommand += auxiliary.getLibrary(libraryName, platform.jenkinsLabel, null, sameOrg)
+            getDependenciesCommand += auxiliary.getLibrary(libraryName, platform.jenkinsLabel, 'develop', sameOrg)
         }
     }
 
@@ -47,23 +47,29 @@ def runTestCommand (platform, project, gfilter)
                 """
 
     platform.runCommand(this, command)
-    junit "${project.paths.project_build_prefix}/build/release/clients/staging/*.xml"
 }
 
 def runCoverageCommand (platform, project, gfilter, String dirmode = "release")
 {
-    //Temporary workaround due to bug in container
-    String centos7Workaround = platform.jenkinsLabel.contains('centos7') ? 'export LD_LIBRARY_PATH=\$LD_LIBRARY_PATH:/opt/rocm/lib64/' : ''
+    String commitSha
+    String repoUrl
+    (commitSha, repoUrl) = util.getGitHubCommitInformation(project.paths.project_src_prefix)
 
-    def command = """#!/usr/bin/env bash
-                set -x
-                cd ${project.paths.project_build_prefix}/build/${dirmode}
-                export LD_LIBRARY_PATH=/opt/rocm/lib/
-                ${centos7Workaround}
-                GTEST_LISTENER=NO_PASS_LINE_IN_LOG make coverage_cleanup coverage GTEST_FILTER=${gfilter}-*known_bug*
-            """
+    withCredentials([string(credentialsId: "mathlibs-codecov-token-rocalution", variable: 'CODECOV_TOKEN')])
+    {
+        def command = """#!/usr/bin/env bash
+                    set -x
+                    cd ${project.paths.project_build_prefix}/build/${dirmode}
+                    export LD_LIBRARY_PATH=/opt/rocm/lib/
+                    export ROCALUTION_CODE_COVERAGE=1
+                    GTEST_LISTENER=NO_PASS_LINE_IN_LOG make coverage_cleanup coverage GTEST_FILTER=${gfilter}-*known_bug*
+                    curl -Os https://uploader.codecov.io/latest/linux/codecov
+                    chmod +x codecov
+                    ./codecov -v -U \$http_proxy -t ${CODECOV_TOKEN} --file lcoverage/main_coverage.info --name rocALUTION --sha ${commitSha}
+                """
 
-    platform.runCommand(this, command)
+        platform.runCommand(this, command)
+    }
 
     publishHTML([allowMissing: false,
                 alwaysLinkToLastBuild: false,
diff --git a/.jenkins/debug.groovy b/.jenkins/debug.groovy
index 15abd79..fd45e14 100644
--- a/.jenkins/debug.groovy
+++ b/.jenkins/debug.groovy
@@ -18,7 +18,7 @@ def runCI =
     prj.paths.build_command = buildCommand
     prj.compiler.compiler_name = 'c++'
     prj.compiler.compiler_path = 'c++'
-    prj.libraryDependencies = ['rocPRIM', 'rocBLAS-internal', 'rocSPARSE-internal', 'rocRAND']
+    prj.libraryDependencies = ['rocPRIM', 'hipBLAS-common', 'hipBLASLt', 'rocBLAS', 'rocSPARSE', 'rocRAND']
 
     // Define test architectures, optional rocm version argument is available
     def nodes = new dockerNodes(nodeDetails, jobName, prj)
diff --git a/.jenkins/extended.groovy b/.jenkins/extended.groovy
index e40fe5e..d588769 100644
--- a/.jenkins/extended.groovy
+++ b/.jenkins/extended.groovy
@@ -18,7 +18,7 @@ def runCI =
     prj.paths.build_command = buildCommand
     prj.compiler.compiler_name = 'c++'
     prj.compiler.compiler_path = 'c++'
-    prj.libraryDependencies = ['rocPRIM', 'rocBLAS-internal', 'rocSPARSE-internal', 'rocRAND']
+    prj.libraryDependencies = ['rocPRIM', 'hipBLAS-common', 'hipBLASLt', 'rocBLAS', 'rocSPARSE', 'rocRAND']
 
     // Define test architectures, optional rocm version argument is available
     def nodes = new dockerNodes(nodeDetails, jobName, prj)
diff --git a/.jenkins/precheckin.groovy b/.jenkins/precheckin.groovy
index 5c309cb..c198213 100644
--- a/.jenkins/precheckin.groovy
+++ b/.jenkins/precheckin.groovy
@@ -18,7 +18,7 @@ def runCI =
     prj.paths.build_command = buildCommand
     prj.compiler.compiler_name = 'c++'
     prj.compiler.compiler_path = 'c++'
-    prj.libraryDependencies = ['rocPRIM', 'rocBLAS-internal', 'rocSPARSE-internal', 'rocRAND']
+    prj.libraryDependencies = ['rocPRIM', 'hipBLAS-common', 'hipBLASLt', 'rocBLAS', 'rocSPARSE', 'rocRAND']
 
     // Define test architectures, optional rocm version argument is available
     def nodes = new dockerNodes(nodeDetails, jobName, prj)
diff --git a/.jenkins/staticanalysis.groovy b/.jenkins/staticanalysis.groovy
index 8f5c593..4cfcd5e 100644
--- a/.jenkins/staticanalysis.groovy
+++ b/.jenkins/staticanalysis.groovy
@@ -39,8 +39,26 @@ def runCI =
 ci: {
     String urlJobName = auxiliary.getTopJobName(env.BUILD_URL)
 
-    properties(auxiliary.addCommonProperties([pipelineTriggers([cron('0 1 * * 6')])]))
-    stage(urlJobName) {
-        runCI([ubuntu20:['cpu']], urlJobName)
+    def propertyList = ["compute-rocm-dkms-no-npi-hipclang":[pipelineTriggers([cron('0 1 * * 0')])],
+                        "rocm-docker":[]]
+    propertyList = auxiliary.appendPropertyList(propertyList)
+
+    def jobNameList = ["compute-rocm-dkms-no-npi-hipclang":[]]
+    jobNameList = auxiliary.appendJobNameList(jobNameList)
+
+    propertyList.each
+    {
+        jobName, property->
+        if (urlJobName == jobName)
+            properties(auxiliary.addCommonProperties(property))
+    }
+
+    jobNameList.each
+    {
+        jobName, nodeDetails->
+        if (urlJobName == jobName)
+            stage(jobName) {
+                runCI(nodeDetails, jobName)
+            }
     }
 }
diff --git a/.jenkins/staticlibrary.groovy b/.jenkins/staticlibrary.groovy
index da565db..b06385d 100644
--- a/.jenkins/staticlibrary.groovy
+++ b/.jenkins/staticlibrary.groovy
@@ -16,9 +16,9 @@ def runCI =
     def prj = new rocProject('rocALUTION', 'Static Library PreCheckin')
     // customize for project
     prj.paths.build_command = buildCommand
-    prj.compiler.compiler_name = 'hipcc'
-    prj.compiler.compiler_path = '/opt/rocm/bin/hipcc'
-    prj.libraryDependencies = ['rocPRIM', 'rocBLAS-internal', 'rocSPARSE-internal', 'rocRAND']
+    prj.compiler.compiler_name = 'amdclang++'
+    prj.compiler.compiler_path = '/opt/rocm/bin/amdclang++'
+    prj.libraryDependencies = ['rocPRIM', 'hipBLAS-common', 'hipBLASLt', 'rocBLAS', 'rocSPARSE', 'rocRAND']
 
     // Define test architectures, optional rocm version argument is available
     def nodes = new dockerNodes(nodeDetails, jobName, prj)
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index 43a0890..4e3a24d 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -6,9 +6,13 @@ version: 2
 sphinx:
    configuration: docs/conf.py
 
-formats: [htmlzip]
+formats: [htmlzip, epub]
 
 python:
-   version: "3.8"
    install:
-   - requirements: docs/.sphinx/requirements.txt
+   - requirements: docs/sphinx/requirements.txt
+
+build:
+   os: ubuntu-22.04
+   tools:
+      python: "3.10"
diff --git a/CHANGELOG.md b/CHANGELOG.md
index e168a58..fe92d69 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,142 +1,306 @@
-# Change Log for rocALUTION
+# Changelog for rocALUTION
 
-Full documentation for rocALUTION is available at [rocalution.readthedocs.io](https://rocalution.readthedocs.io/en/latest/).
+Full documentation forrocALUTION is available at [https://rocm.docs.amd.com/projects/rocALUTION/en/latest/](https://rocm.docs.amd.com/projects/rocALUTION/en/latest/).
+
+## rocALUTION 4.0.1 for ROCm 7.1.0
 
-## rocALUTION 2.1.11 for ROCm 5.7.0
 ### Added
-- Added support for gfx940, gfx941 and gfx942
-### Improved
-- Fixed OpenMP runtime issue with Windows toolchain
+* Added support for gfx950.
 
-## rocALUTION 2.1.9 for ROCm 5.6.0
-### Improved
-- Fixed synchronization issues in level 1 routines
+### Changed
+* Updated the default build standard to C++17 when compiling rocALUTION from source (previously C++14).
+
+### Optimized
+* Improved and expanded user documentation.
+
+### Resolved issues
+* Fixed a bug in the GPU hashing algorithm that occurred when not compiling with -O2/-O3.
+* Fixed an issue with the SPAI preconditioner when using complex numbers.
+
+## rocALUTION 3.2.3 for ROCm 6.4.1
 
-## rocALUTION 2.1.8 for ROCm 5.5.0
 ### Added
-- Added build support for Navi32
-### Improved
-- Fixed a typo in MPI backend
-- Fixed a bug with the backend when HIP support is disabled
-- Fixed a bug in SAAMG hierarchy building on HIP backend
-- Improved SAAMG hierarchy build performance on HIP backend
+* The `-a` option has been added to the `rmake.py` build script. This option provides a way to select specific architectures when building on Windows.
+
+### Resolved issues
+* Fixed an issue where the `HIP_PATH` environment variable was being ignored when compiling on Windows.
+
+## rocALUTION 3.2.2 for ROCm 6.4.0
+
 ### Changed
-- LocalVector::GetIndexValues(ValueType\*) is deprecated, use LocalVector::GetIndexValues(const LocalVector&, LocalVector\*) instead
-- LocalVector::SetIndexValues(const ValueType\*) is deprecated, use LocalVector::SetIndexValues(const LocalVector&, const LocalVector&) instead
-- LocalMatrix::RSDirectInterpolation(const LocalVector&, const LocalVector&, LocalMatrix\*, LocalMatrix\*) is deprecated, use LocalMatrix::RSDirectInterpolation(const LocalVector&, const LocalVector&, LocalMatrix\*) instead
-- LocalMatrix::RSExtPIInterpolation(const LocalVector&, const LocalVector&, bool, float, LocalMatrix\*, LocalMatrix\*) is deprecated, use LocalMatrix::RSExtPIInterpolation(const LocalVector&, const LocalVector&, bool, LocalMatrix\*) instead
-- LocalMatrix::RugeStueben() is deprecated
-- LocalMatrix::AMGSmoothedAggregation(ValueType, const LocalVector&, const LocalVector&, LocalMatrix\*, LocalMatrix\*, int) is deprecated, use LocalMatrix::AMGAggregation(ValueType, const LocalVector&, const LocalVector&, LocalMatrix\*, int) instead
-- LocalMatrix::AMGAggregation(const LocalVector&, LocalMatrix\*, LocalMatrix\*) is deprecated, use LocalMatrix::AMGAggregation(const LocalVector&, LocalMatrix\*) instead
+* Improved documentation
+
+## rocALUTION 3.2.1 for ROCm 6.3.0
 
-## rocALUTION 2.1.3 for ROCm 5.4.0
 ### Added
-- Added build support for Navi31 and Navi33
-- Added support for non-squared global matrices
-### Improved
-- Fixed a memory leak in MatrixMult on HIP backend
-- Global structures can now be used with a single process
+
+* Support for gfx1200, gfx1201, and gfx1151.
+
 ### Changed
-- Switched GTest death test style to 'threadsafe'
-- GlobalVector::GetGhostSize() is deprecated and will be removed
-- ParallelManager::GetGlobalSize(), ParallelManager::GetLocalSize(), ParallelManager::SetGlobalSize() and ParallelManager::SetLocalSize() are deprecated and will be removed
-- Vector::GetGhostSize() is deprecated and will be removed
-- Multigrid::SetOperatorFormat(unsigned int) is deprecated and will be removed, use Multigrid::SetOperatorFormat(unsigned int, int) instead
-- RugeStuebenAMG::SetCouplingStrength(ValueType) is deprecated and will be removed, use SetStrengthThreshold(float) instead
+
+* Changed the default compiler from `hipcc` to `amdclang` in the installation script and cmake files.
+* Changed the address sanitizer build targets. Now only `gfx908:xnack+`, `gfx90a:xnack+`, `gfx940:xnack+`, `gfx941:xnack+`, and `gfx942:xnack+` are built with `BUILD_ADDRESS_SANITIZER=ON`.
+
+### Resolved issues
+
+* Fix hang in `RS-AMG` for Navi on some specific matrix sparsity patterns.
+* Fix wrong results in `Apply` on multi-GPU setups.
+
+## rocALUTION 3.2.0 for ROCm 6.2.0
+
+### Additions
+* New file I/O based on rocsparse I/O format
+* `GetConvergenceHistory` for ItILU0 preconditioner
+
+### Deprecations
+* `LocalMatrix::ReadFileCSR`
+* `LocalMatrix::WriteFileCSR`
+* `GlobalMatrix::ReadFileCSR`
+* `GlobalMatrix::WriteFileCSR`
+
+## rocALUTION 3.1.1 for ROCm 6.1.0
+
+### Additions
+
+* `TripleMatrixProduct` functionality for `GlobalMatrix`
+* Multi-Node/GPU support for `UA-AMG`, `SA-AMG` and `RS-AMG`
+* Iterative ILU0 preconditioner `ItILU0`
+* Iterative triangular solve, selectable via `SolverDecr` class
+
+### Deprecations
+
+* `LocalMatrix::AMGConnect`
+* `LocalMatrix::AMGAggregate`
+* `LocalMatrix::AMGPMISAggregate`
+* `LocalMatrix::AMGSmoothedAggregation`
+* `LocalMatrix::AMGAggregation`
+* `PairwiseAMG`
+
+### Known Issues
+* `PairwiseAMG` does currently not support matrix sizes that exceed int32 range
+* `PairwiseAMG` might fail building the hierarchy on certain input matrices
+
+## rocALUTION 3.0.3 for ROCm 6.0.0
+
+### Additions
+
+* Support for 64bit integer vectors
+* Inclusive and exclusive sum functionality for vector classes
+* Transpose functionality for `GlobalMatrix` and `LocalMatrix`
+* `TripleMatrixProduct` functionality for `LocalMatrix`
+* `Sort()` function for `LocalVector` class
+* Multiple stream support to the HIP backend
+
+### Optimizations
+
+* `GlobalMatrix::Apply()` now uses multiple streams to better hide communication
+
+### Changes
+
+* Matrix dimensions and number of non-zeros are now stored using 64-bit integers
+* Improved the ILUT preconditioner
+
+### Deprecations
+
+* `LocalVector::GetIndexValues(ValueType*)`
+* `LocalVector::SetIndexValues(const ValueType*)`
+* `LocalMatrix::RSDirectInterpolation(const LocalVector&, const LocalVector&, LocalMatrix*, LocalMatrix*)`
+* `LocalMatrix::RSExtPIInterpolation(const LocalVector&, const LocalVector&, bool, float, LocalMatrix*, LocalMatrix*)`
+* `LocalMatrix::RugeStueben()`
+* `LocalMatrix::AMGSmoothedAggregation(ValueType, const LocalVector&, const LocalVector&, LocalMatrix*, LocalMatrix*, int)`
+* `LocalMatrix::AMGAggregation(const LocalVector&, LocalMatrix*, LocalMatrix*)`
+
+### Fixes
+
+* Unit tests no longer ignore BCSR block dimension
+* Fixed documentation typos
+* Bug in multi-coloring for non-symmetric matrix patterns
+
+## rocALUTION 2.1.11 for ROCm 5.7.0
+
+### Additions
+
+* Support for gfx940, gfx941, and gfx942
+
+### Fixes
+
+* OpenMP runtime issue with Windows toolchain
+
+## rocALUTION 2.1.9 for ROCm 5.6.0
+
+### Fixes
+
+* Synchronization issues in level 1 routines
+
+## rocALUTION 2.1.8 for ROCm 5.5.0
+
+### Additions
+
+* Build support for Navi32
+
+### Fixes
+
+* Typo in MPI backend
+* Bug with the backend when HIP support is disabled
+* Bug in SAAMG hierarchy building on the HIP backend
+* Improved SAAMG hierarchy build performance on the HIP backend
+
+### Deprecations
+
+* `LocalVector::GetIndexValues(ValueType*)`: use
+  `LocalVector::GetIndexValues(const LocalVector&, LocalVector*)` instead
+* `LocalVector::SetIndexValues(const ValueType*)`: use
+  `LocalVector::SetIndexValues(const LocalVector&, const LocalVector&)` instead
+* `LocalMatrix::RSDirectInterpolation(const LocalVector&, const LocalVector&, LocalMatrix*, LocalMatrix*)`:
+  use `LocalMatrix::RSDirectInterpolation(const LocalVector&, const LocalVector&, LocalMatrix*)`
+  instead
+* `LocalMatrix::RSExtPIInterpolation(const LocalVector&, const LocalVector&, bool, float, LocalMatrix*, LocalMatrix*)`:
+  use `LocalMatrix::RSExtPIInterpolation(const LocalVector&, const LocalVector&, bool, LocalMatrix*)`
+  instead
+* `LocalMatrix::RugeStueben()`
+* `LocalMatrix::AMGSmoothedAggregation(ValueType, const LocalVector&, const LocalVector&, LocalMatrix*, LocalMatrix*, int)`:
+  use `LocalMatrix::AMGAggregation(ValueType, const LocalVector&, const LocalVector&, LocalMatrix*, int)`
+  instead
+* `LocalMatrix::AMGAggregation(const LocalVector&, LocalMatrix*, LocalMatrix*)`: use
+  `LocalMatrix::AMGAggregation(const LocalVector&, LocalMatrix*)` instead
+
+## rocALUTION 2.1.3 for ROCm 5.4.0
+
+### Additions
+
+* Build support for Navi31 and Navi33
+* Support for non-squared global matrices
+
+### Fixes
+
+* Memory leak in MatrixMult on HIP backend
+* Global structures can now be used with a single process
+
+### Changes
+
+* Switched GTest death test style to 'threadsafe'
+* Removed the native compiler option that was used during default library compilation
+
+### Deprecations
+
+* `GlobalVector::GetGhostSize()`
+* `ParallelManager::GetGlobalSize(), ParallelManager::GetLocalSize()`, `ParallelManager::SetGlobalSize()`,
+  and `ParallelManager::SetLocalSize()`
+* `Vector::GetGhostSize()`
+* `Multigrid::SetOperatorFormat(unsigned int)`: use `Multigrid::SetOperatorFormat(unsigned int, int)`
+  instead
+* `RugeStuebenAMG::SetCouplingStrength(ValueType)`: use `SetStrengthThreshold(float)` instead
 
 ## rocALUTION 2.1.0 for ROCm 5.3.0
-### Added
-- Benchmarking tool
-- Ext+I Interpolation with sparsify strategies added for RS-AMG
-### Improved
-- ParallelManager
+
+### Additions
+
+* Benchmarking tool
+* Ext+I Interpolation with sparsify strategies added for RS-AMG
+
+### Optimizations
+
+* ParallelManager
 
 ## rocALUTION 2.0.3 for ROCm 5.2.0
-### Added
-- Packages for test and benchmark executables on all supported OSes using CPack.
+
+### Additions
+
+* New packages for test and benchmark executables on all supported operating systems using CPack
 
 ## rocALUTION 2.0.2 for ROCm 5.1.0
-### Added
-- Added out-of-place matrix transpose functionality
-- Added LocalVector<bool>
+
+### Additions
+
+* Added out-of-place matrix transpose functionality
+* Added LocalVector<bool>
 
 ## rocALUTION 2.0.1 for ROCm 5.0.0
-### Changed
-- Removed deprecated GlobalPairwiseAMG class, please use PairwiseAMG instead.
-- Changed to C++ 14 Standard
-### Improved
-- Added sanitizer option
-- Improved documentation
+
+### Changes
+
+* Changed to C++ 14 Standard
+* Added sanitizer option
+* Improved documentation
+
+### Deprecations
+
+* `GlobalPairwiseAMG` class: use `PairwiseAMG` instead
 
 ## rocALUTION 1.13.2 for ROCm 4.5.0
-### Added
-- Address sanitizer build option added
-- Enabling beta support for Windows 10
-### Changed
-- Deprecated GlobalPairwiseAMG, please use PairwiseAMG instead. GlobalPairwiseAMG will be removed in a future major release.
-- Packaging split into a runtime package called rocalution and a development package called rocalution-devel. The development package depends on runtime. The runtime package suggests the development package for all supported OSes except CentOS 7 to aid in the transition. The suggests feature in packaging is introduced as a deprecated feature and will be removed in a future rocm release.
-### Improved
-- (A)MG smoothing and convergence rate improvement
-- Improved IDR Gram-Schmidt process
-- (A)MG solving phase optimization
+
+### Additions
+
+* AddressSanitizer build option
+* Enabled beta support for Windows 10
+
+### Changes
+
+* Packaging has been split into a runtime package (`rocalution`) and a development package
+  (`rocalution-devel`):
+  The development package depends on the runtime package. When installing the runtime package,
+  the package manager will suggest the installation of the development package to aid users
+  transitioning from the previous version's combined package. This suggestion by package manager is
+  for all supported operating systems (except CentOS 7) to aid in the transition. The `suggestion`
+  feature in the runtime package is introduced as a deprecated feature and will be removed in a future
+  ROCm release.
+
+### Deprecations
+
+* `GlobalPairwiseAMG`: use `PairwiseAMG` instead
+
+### Optimizations
+
+* Improved (A)MG smoothing and convergence rate
+* Improved IDR Gram-Schmidt process
+* Optimized (A)MG solving phase
 
 ## rocALUTION 1.12.1 for ROCm 4.3.0
-### Added
-- support for gfx90a target
-- support for gfx1030 target
-### Improved
-- install script
-### Known Issues
-- none
 
-## rocALUTION 1.11.5 for ROCm 4.2.0
-### Added
-- none
-### Known Issues
-- none
+### Additions
 
-## rocALUTION 1.11.5 for ROCm 4.1.0
-### Added
-- none
-### Known Issues
-- none
+* Support for gfx90a target
+* Support for gfx1030 target
+
+### Optimizations
+
+* Install script
 
 ## rocALUTION 1.11.5 for ROCm 4.0.0
-### Added
-- Add changelog
-- Fixing NaN issues
-- update to debian package name
-- bcsr format support.
-- cmake files adjustments.
+
+### Additions
+
+* Changelog
+* Block compressed sparse row (BCSRR) format support
+
+### Changes
+
+* Update to the Debian package name
+* CMake file adjustments
+
+### Fixes
+
+* NaN issues
 
 ## rocALUTION 1.10 for ROCm 3.9
-### Added
-- rocRAND to support GPU sampling of random data.
-### Known Issues
-- none
+
+### Additions
+
+* rocRAND support for GPU sampling of random data
 
 ## rocALUTION 1.9.3 for ROCm 3.8
-### Added
-- csr2dense and dense2csr to HIP backend.
-### Known Issues
-- none
 
-## rocALUTION 1.9.1 for ROCm 3.7
-### Added
-- none
-### Known Issues
-- none
+### Additions
 
-## rocALUTION 1.9.1 for ROCm 3.6
-### Added
-- none
-### Known Issues
-- none
+* `csr2dense` and `dense2csr` to HIP backend
 
 ## rocALUTION 1.9.1 for ROCm 3.5
-### Added
-- static build
-- BCSR matrix format for SpMV
-- Bug fixing in conversion from CSR to HYB format.
-### Known Issues
-- none
+
+### Additions
+
+* Static build
+* BCSR matrix format for SpMV
+
+### Fixes
+
+* Bug in conversion from CSR to HYB format
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 62e94be..8a622f4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,5 +1,5 @@
 # ########################################################################
-# Copyright (C) 2018-2022 Advanced Micro Devices, Inc. All rights Reserved.
+# Copyright (C) 2018-2025 Advanced Micro Devices, Inc. All rights Reserved.
 #
 # Permission is hereby granted, free of charge, to any person obtaining a copy
 # of this software and associated documentation files (the "Software"), to deal
@@ -38,7 +38,8 @@ list(APPEND CMAKE_PREFIX_PATH ${ROCM_PATH}/llvm ${ROCM_PATH})
 list(APPEND CMAKE_MODULE_PATH
      ${CMAKE_CURRENT_SOURCE_DIR}/cmake
      ${ROCM_PATH}/lib/cmake/hip
-     ${ROCM_PATH}/hip/cmake)
+     ${ROCM_PATH}/hip/cmake
+     ${ROCM_PATH}/cmake)
 
 # Set a default build type if none was specified
 if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
@@ -56,9 +57,12 @@ endif()
 project(rocalution LANGUAGES CXX)
 
 # Build flags
-set(CMAKE_CXX_STANDARD 14)
+set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)
+set(CMAKE_HIP_STANDARD 17)
+set(CMAKE_HIP_STANDARD_REQUIRED ON)
+set(CMAKE_HIP_EXTENSIONS OFF)
 
 # Build options
 option(BUILD_SHARED_LIBS "Build rocALUTION as a shared library" ON)
@@ -68,31 +72,67 @@ option(BUILD_CLIENTS_SAMPLES "Build examples" ON)
 option(BUILD_VERBOSE "Output additional build information" OFF)
 option(BUILD_CODE_COVERAGE "Build with code coverage enabled" OFF)
 option(BUILD_ADDRESS_SANITIZER "Build with address sanitizer enabled" OFF)
-option(BUILD_FILE_REORG_BACKWARD_COMPATIBILITY "Build with file/folder reorg with backward compatibility enabled" ON)
+option(BUILD_GLOBALTYPE_64 "Support global number of rows / columns exceeding 32 bits" OFF)
+option(BUILD_LOCALTYPE_64 "Support local number of rows / columns exceeding 32 bits" OFF)
+option(BUILD_PTRTYPE_64 "Support local number of non-zeros exceeding 32 bits" OFF)
+option(BUILD_OPTCPU "Enable all instruction subsets supported by the local machine" OFF)
+option(BUILD_SUPPORT_COMPLEX "Enable complex number support" ON)
 
 # Dependencies
 include(cmake/Dependencies.cmake)
 
-# Availability of rocm_check_target_ids command assures that we can also build
-# for gfx90a target
-if(COMMAND rocm_check_target_ids)
-  set(DEFAULT_AMDGPU_TARGETS "gfx803;gfx900:xnack-;gfx906:xnack-;gfx908:xnack-;gfx1030;gfx90a:xnack-;gfx90a:xnack+;gfx940;gfx941;gfx942;gfx1100;gfx1101;gfx1102")
+if(BUILD_ADDRESS_SANITIZER)
+  set(DEFAULT_AMDGPU_TARGETS
+    gfx908:xnack+
+    gfx90a:xnack+
+    gfx942:xnack+
+  )
 else()
-  set(DEFAULT_AMDGPU_TARGETS "gfx803;gfx900:xnack-;gfx906:xnack-;gfx908:xnack-;gfx1030;gfx1100;gfx1101;gfx1102")
+  set(DEFAULT_AMDGPU_TARGETS
+    gfx803
+    gfx900:xnack-
+    gfx906:xnack-
+    gfx908:xnack-
+    gfx90a:xnack-
+    gfx90a:xnack+
+    gfx942
+    gfx950
+    gfx1030
+    gfx1100
+    gfx1101
+    gfx1102
+    gfx1151
+    gfx1200
+    gfx1201
+  )
 endif()
-set(AMDGPU_TARGETS "${DEFAULT_AMDGPU_TARGETS}" CACHE STRING "List of specific machine types for library to target")
-
-# Find HIP package
-find_package(HIP)
-
-if (NOT HIP_FOUND)
-  message("-- HIP not found. Compiling WITHOUT HIP support.")
+if(AMDGPU_TARGETS AND NOT GPU_TARGETS)
+  message( DEPRECATION "AMDGPU_TARGETS use is deprecated. Use GPU_TARGETS." )
+endif()
+set(AMDGPU_TARGETS "${DEFAULT_AMDGPU_TARGETS}" CACHE STRING "Target default GPUs if AMDGPU_TARGETS is not defined. (Deprecated, prefer GPU_TARGETS)")
+# Don't force, users should be able to override GPU_TARGETS at the command line if desired
+set(GPU_TARGETS "${AMDGPU_TARGETS}" CACHE STRING "GPU architectures to build for")
+
+include(CheckLanguage)
+include(CMakeDependentOption)
+check_language(HIP)
+cmake_dependent_option(USE_HIPCXX "Use CMake HIP language support" OFF CMAKE_HIP_COMPILER OFF)
+if(USE_HIPCXX)
+  enable_language(HIP)
 else()
-  option(SUPPORT_HIP "Compile WITH HIP support." ON)
+  find_package(HIP MODULE) # hip_add_library is only provided by the find module
+  if(NOT HIP_FOUND)
+    message("-- HIP not found. Compiling WITHOUT HIP support.")
+  endif()
 endif()
 
+cmake_dependent_option(SUPPORT_HIP "Compile WITH HIP support" ON "USE_HIPCXX OR HIP_FOUND" OFF)
+
 # HIP related library dependencies
 if(SUPPORT_HIP)
+  if( CMAKE_CXX_COMPILER_ID MATCHES "Clang" )
+    find_package( hip REQUIRED CONFIG PATHS ${HIP_DIR} ${ROCM_PATH} /opt/rocm )
+  endif( )
   find_package(rocblas REQUIRED)
   find_package(rocsparse REQUIRED)
   find_package(rocprim REQUIRED)
@@ -112,9 +152,7 @@ endif()
 
 
 # Setup version
-set(VERSION_STRING "2.1.11")
-rocm_setup_version(VERSION ${VERSION_STRING})
-set(rocalution_SOVERSION 0.1)
+rocm_setup_version(VERSION "4.0.1")
 
 if(BUILD_CLIENTS_SAMPLES OR BUILD_CLIENTS_BENCHMARKS OR BUILD_CLIENTS_TESTS)
   set( BUILD_CLIENTS ON )
@@ -137,7 +175,12 @@ endif()
 
 # Package specific CPACK vars
 if(SUPPORT_HIP)
-  rocm_package_add_dependencies(DEPENDS "hip-rocclr >= 4.0.0" "rocsparse >= 1.12.10" "rocblas >= 2.22.0" "rocrand >= 0.0.1")
+  if(BUILD_ADDRESS_SANITIZER)
+    set(DEPENDS_HIP_RUNTIME "hip-runtime-amd-asan >= 4.5.0" )
+  else()
+    set(DEPENDS_HIP_RUNTIME "hip-runtime-amd >= 4.5.0" )
+  endif()
+  rocm_package_add_dependencies(DEPENDS "${DEPENDS_HIP_RUNTIME}" "rocsparse >= 1.12.10" "rocblas >= 2.22.0" "rocrand >= 2.1.0")
 endif()
 
 set(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_CURRENT_SOURCE_DIR}/LICENSE.md")
@@ -166,7 +209,7 @@ set(ROCALUTION_CONFIG_DIR "\${CPACK_PACKAGING_INSTALL_PREFIX}/${CMAKE_INSTALL_LI
 
 rocm_create_package(
   NAME ${package_name}
-  DESCRIPTION "Radeon Open Compute library for sparse linear systems"
+  DESCRIPTION "ROCm library for sparse linear systems"
   MAINTAINER "rocALUTION Maintainer <[e-mail]>"
   LDCONFIG
   LDCONFIG_DIR ${ROCALUTION_CONFIG_DIR}
@@ -202,18 +245,14 @@ if(BUILD_CODE_COVERAGE)
   add_custom_target(coverage_output
     DEPENDS coverage_analysis
     COMMAND mkdir -p lcoverage
-    COMMAND echo "\\#!/bin/bash" > llvm-gcov.sh
-    COMMAND echo "\\# THIS FILE HAS BEEN GENERATED" >> llvm-gcov.sh
-    COMMAND printf "exec /opt/rocm/llvm/bin/llvm-cov gcov $$\\@" >> llvm-gcov.sh
-    COMMAND chmod +x llvm-gcov.sh
     )
 
   #
   # Generate coverage output.
   #
   add_custom_command(TARGET coverage_output
-    COMMAND lcov --directory . --base-directory . --gcov-tool ${CMAKE_BINARY_DIR}/llvm-gcov.sh --capture -o lcoverage/raw_main_coverage.info
-    COMMAND lcov --remove lcoverage/raw_main_coverage.info "'/opt/*'" "'/usr/*'" -o lcoverage/main_coverage.info
+    COMMAND lcov --directory . --base-directory . --capture -o lcoverage/raw_main_coverage.info
+    COMMAND lcov --remove lcoverage/raw_main_coverage.info "'${CMAKE_SOURCE_DIR}/src/utils/*'" "'${CMAKE_SOURCE_DIR}/src/base/host/host_io.*'" "'${CMAKE_SOURCE_DIR}/clients/*'" "'${CMAKE_SOURCE_DIR}/build/*'" "'/opt/*'" "'/usr/*'" -o lcoverage/main_coverage.info
     COMMAND genhtml lcoverage/main_coverage.info --output-directory lcoverage
     )
 
diff --git a/LICENSE.md b/LICENSE.md
index 7a44e61..e41c5d8 100644
--- a/LICENSE.md
+++ b/LICENSE.md
@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (C) 2018-2023 Advanced Micro Devices, Inc. All rights reserved.
+Copyright (C) 2024 Advanced Micro Devices, Inc. All rights reserved.
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
diff --git a/README.md b/README.md
index b9b6929..d63704f 100644
--- a/README.md
+++ b/README.md
@@ -1,28 +1,52 @@
 # rocALUTION
-rocALUTION is a sparse linear algebra library with focus on exploring fine-grained parallelism on top of AMD's Radeon Open eCosystem Platform [ROCm][] runtime and toolchains, targeting modern CPU and GPU platforms. Based on C++ and [HIP][], it provides a portable, generic and flexible design that allows seamless integration with other scientific software packages.
+
+> [!NOTE]
+> The published documentation is available at [rocALUTION](https://rocm.docs.amd.com/projects/rocALUTION/en/latest/) in an organized, easy-to-read format, with search and a table of contents. The documentation source files reside in the `docs` folder of this repository. As with all ROCm projects, the documentation is open source. For more information on contributing to the documentation, see [Contribute to ROCm documentation](https://rocm.docs.amd.com/en/latest/contribute/contributing.html).
+
+rocALUTION is a sparse linear algebra library that can be used to explore fine-grained parallelism on
+top of the [ROCm](https://github.com/ROCm/ROCm) platform runtime and toolchains.
+Based on C++ and [HIP](https://github.com/ROCm/HIP/), rocALUTION
+provides a portable, generic, and flexible design that allows seamless integration with other scientific
+software packages.
+
+rocALUTION offers various backends for different (parallel) hardware:
+
+* Host
+* [OpenMP](http://www.openmp.org/): Designed for multi-core CPUs
+* [HIP](https://github.com/ROCm/HIP/): Designed for ROCm-compatible devices
+* [MPI](https://www.open-mpi.org/): Designed for multi-node clusters and multi-GPU setups
+
+## Requirements
+
+To use rocALUTION on GPU devices, you must first install the
+[rocBLAS](https://github.com/ROCm/rocBLAS),
+[rocSPARSE](https://github.com/ROCm/rocSPARSE), and
+[rocRAND](https://github.com/ROCm/rocRAND) libraries. You can install these from
+the ROCm repository, the GitHub 'releases' tab, or you can manually compile them.
 
 ## Documentation
-The latest rocALUTION documentation and API description can be found [here][].
 
-### How to build documentation
+Documentation for rocALUTION is available at
+[https://rocm.docs.amd.com/projects/rocALUTION/en/latest/](https://rocm.docs.amd.com/projects/rocALUTION/en/latest/).
 
-Run the following steps to build documentation.
+To build our documentation locally, use the following code:
 
-```
+```bash
 cd docs
 
-pip3 install -r .sphinx/requirements.txt
+pip3 install -r sphinx/requirements.txt
 
 python3 -m sphinx -T -E -b html -d _build/doctrees -D language=en . _build/html
 ```
 
-## Quickstart rocALUTION build
+## Build
 
-#### CMake 3.5 or later
-All compiler specifications are determined automatically. The compilation process can be performed by
-```
+You can compile rocALUTION using CMake 3.5 or later. Note that all compiler specifications are
+determined automatically.
+
+```bash
 # Clone rocALUTION using git
-git clone https://github.com/ROCmSoftwarePlatform/rocALUTION.git
+git clone https://github.com/ROCm/rocALUTION.git
 
 # Go to rocALUTION directory, create and change to build directory
 cd rocALUTION; mkdir build; cd build
@@ -40,69 +64,56 @@ cmake .. -DSUPPORT_HIP=ON -DROCM_PATH=/opt/rocm/
 make
 ```
 
-#### Simple test
-You can test the installation by running a CG solver on a Laplace matrix:
-```
+To test your installation, run a CG solver on a Laplacian matrix:
+
+```bash
 cd rocALUTION; cd build
 wget ftp://math.nist.gov/pub/MatrixMarket2/Harwell-Boeing/laplace/gr_30_30.mtx.gz
 gzip -d gr_30_30.mtx.gz
 ./clients/staging/cg gr_30_30.mtx
 ```
 
-## Overview
-
-#### Backends
-rocALUTION offers various backends for different (parallel) hardware.
-*  Host
-*  [OpenMP][] - designed for multi-core CPUs
-*  [HIP][]    - designed for ROCm compatible devices
-*  [MPI][]    - designed for multi-node clusters and multi-GPU setups
-
-#### Easy to use
-Syntax and structure of the library provide fast learning curves. With the help of the examples, anyone can try out the library - no knowledge in multi-core or GPU programming is required.
-
-#### Requirements
-There are no hardware requirements to install and run rocALUTION. If GPU devices are available, rocALUTION will use them.
-In order to use rocALUTION on GPU devices, you will need to make sure that [rocBLAS][], [rocSPARSE][] and [rocRAND][] libraries are installed on your system. You can install them from ROCm repository, from github releases tab or manually compile them yourself.
-
-#### Iterative solvers
-*  Fixed-Point iteration schemes - Jacobi, (Symmetric) Gauss-Seidel, SOR, SSOR
-*  Krylov subspace methods - CR, CG, BiCGStab, BiCGStab(*l*), GMRES, IDR, QMRCGSTAB, Flexible CG/GMRES
-*  Mixed-precision defect-correction scheme
-*  Chebyshev iteration scheme
-*  Multigrid - geometric and algebraic
-
-#### Preconditioners
-*  Matrix splitting schemes - Jacobi, (multi-colored) (symmetric) Gauss-Seidel, SOR, SSOR
-*  Factorization schemes    - ILU(*0*), ILU(*p*) (based on levels), ILU(*p,q*) (power(*q*)-pattern method), multi-elimination ILU (nested/recursive), ILUT (based on threshold), IC(*0*)
-*  Approximate Inverses - Chebyshev matrix-valued polynomial, SPAI, FSAI, TNS
-*  Diagonal-based preconditioner for Saddle-point problems
-*  Block-type of sub-preconditioners/solvers
-*  (Restricted) Additive Schwarz
-*  Variable type of preconditioners
-
-#### Sparse matrix formats
-*  Compressed Sparse Row (CSR)
-*  Modified Compressed Sparse Row (MCSR)
-*  Dense (DENSE)
-*  Coordinate (COO)
-*  ELL
-*  Diagonal (DIA)
-*  Hybrid ELL+COO (HYB)
-
-#### Generic and robust design
-rocALUTION is based on a generic and robust design, allowing expansion in the direction of new solvers and preconditioners and support for various hardware types. Furthermore, the design of the library allows the use of all solvers as preconditioners in other solvers, for example you can define a CG solver with a multi-elimination preconditioner, where the last-block is preconditioned with another Chebyshev iteration method which is preconditioned with a multi-colored symmetric Gauss-Seidel scheme.
-
-#### Portable code and results
-All code based on rocALUTION is portable and independent of the hardware, it will compile and run on any supported platform. All solvers and preconditioners are based on a single source code implementation, which delivers portable results across all backends (variations are possible due to different rounding modes on the hardware). The only difference which you can see for a hardware change is the performance variation.
-
-
-
-[ROCm]: https://github.com/RadeonOpenCompute/ROCm
-[HIP]: https://github.com/GPUOpen-ProfessionalCompute-Tools/HIP/
-[OpenMP]: http://www.openmp.org/
-[MPI]: https://www.open-mpi.org/
-[rocBLAS]: https://github.com/ROCmSoftwarePlatform/rocBLAS
-[rocSPARSE]: https://github.com/ROCmSoftwarePlatform/rocSPARSE
-[rocRAND]: https://github.com/ROCmSoftwarePlatform/rocRAND
-[here]: https://rocalution.readthedocs.io
+## General information
+
+rocALUTION is based on a generic and robust design that allows expansion in the direction of new
+solvers and preconditioners with support for various hardware types. The library's design allows the
+use of all solvers as preconditioners in other solvers. For example, you can define a CG solver with a
+multi-elimination preconditioner, in which the last-block is preconditioned with another Chebyshev
+iteration method that itself is preconditioned with a multi-colored symmetric Gauss-Seidel scheme.
+
+### Iterative solvers
+
+* Fixed-point iteration schemes: Jacobi, (Symmetric) Gauss-Seidel, SOR, SSOR
+* Krylov subspace methods: CR, CG, BiCGStab, BiCGStab(*l*), GMRES, IDR, QMRCGSTAB,
+  Flexible CG/GMRES
+* Mixed-precision defect correction scheme
+* Chebyshev iteration scheme
+* Multigrid: Geometric and algebraic
+
+### Preconditioners
+
+* Matrix splitting schemes: Jacobi, (multi-colored) (symmetric) Gauss-Seidel, SOR, SSOR
+* Factorization schemes: ILU(*0*), ILU(*p*) (based on levels), ILU(*p,q*) (power(*q*)-pattern method),
+  multi-elimination ILU (nested/recursive), ILUT (based on threshold), IC(*0*)
+* Approximate Inverses: Chebyshev matrix-valued polynomial, SPAI, FSAI, TNS
+* Diagonal-based preconditioner for Saddle-point problems
+* Block-type of sub-preconditioners/solvers
+* Additive Schwarz (restricted)
+* Variable type of preconditioners
+
+### Sparse matrix formats
+
+* Compressed Sparse Row (CSR)
+* Modified Compressed Sparse Row (MCSR)
+* Dense (DENSE)
+* Coordinate (COO)
+* ELL
+* Diagonal (DIA)
+* Hybrid ELL+COO (HYB)
+
+## Portability
+
+All code based on rocALUTION is portable and hardware-independent. It compiles and runs on any
+supported platform. All solvers and preconditioners are based on a single source code implementation
+that delivers portable results across all backends (note that variations are possible due to different
+hardware rounding modes). The only visible difference between hardware is performance variation.
diff --git a/clients/benchmarks/CMakeLists.txt b/clients/benchmarks/CMakeLists.txt
index 781a24c..aa66ca7 100644
--- a/clients/benchmarks/CMakeLists.txt
+++ b/clients/benchmarks/CMakeLists.txt
@@ -1,5 +1,5 @@
 # ########################################################################
-# Copyright (C) 2022 Advanced Micro Devices, Inc. All rights Reserved.
+# Copyright (C) 2022-2023 Advanced Micro Devices, Inc. All rights Reserved.
 #
 # Permission is hereby granted, free of charge, to any person obtaining a copy
 # of this software and associated documentation files (the "Software"), to deal
@@ -29,6 +29,7 @@ set(ROCALUTION_BENCHMARK_SOURCES
   rocalution_bench_solver_results.cpp
   rocalution_enum_coarsening_strategy.cpp
   rocalution_enum_directsolver.cpp
+  rocalution_enum_itilu0_alg.cpp
   rocalution_enum_itsolver.cpp
   rocalution_enum_matrix_init.cpp
   rocalution_enum_preconditioner.cpp
@@ -40,7 +41,7 @@ set(ROCALUTION_BENCHMARK_SOURCES
 add_executable(rocalution-bench ${ROCALUTION_BENCHMARK_SOURCES} ${ROCALUTION_CLIENTS_COMMON})
 
 # Target compile options
-target_compile_options(rocalution-bench PRIVATE -Wno-unused-command-line-argument -Wall)
+target_compile_options(rocalution-bench PRIVATE -Wall)
 
 # Internal common header
 target_include_directories(rocalution-bench PRIVATE $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/../include>)
diff --git a/clients/benchmarks/rocalution_arguments_config.cpp b/clients/benchmarks/rocalution_arguments_config.cpp
index 61cf6eb..692f1ea 100644
--- a/clients/benchmarks/rocalution_arguments_config.cpp
+++ b/clients/benchmarks/rocalution_arguments_config.cpp
@@ -1,6 +1,6 @@
 /*! \file */
 /* ************************************************************************
-* Copyright (C) 2022 Advanced Micro Devices, Inc. All rights Reserved.
+* Copyright (C) 2022-2023 Advanced Micro Devices, Inc. All rights Reserved.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
@@ -82,6 +82,16 @@ void rocalution_arguments_config::set_description(options_description& desc)
             ADD_OPTION(double, e, 0.05, "ilut tolerance");
             break;
         }
+        case rocalution_bench_solver_parameters::itsolve_tol:
+        {
+            ADD_OPTION(double, e, 1e-03, "iterative solve tolerance (see --iterative_solve).");
+            break;
+        }
+        case rocalution_bench_solver_parameters::itilu0_tol:
+        {
+            ADD_OPTION(double, e, 2e-07, "itilu0 tolerance");
+            break;
+        }
         case rocalution_bench_solver_parameters::mcgs_relax:
         {
             ADD_OPTION(double, e, 1.0, "relaxation coefficient");
@@ -155,6 +165,27 @@ void rocalution_arguments_config::set_description(options_description& desc)
             break;
         }
 
+        case rocalution_bench_solver_parameters::itsolve_max_iter:
+        {
+            ADD_OPTION(int,
+                       e,
+                       30,
+                       "iterative solve maximum number of iterations (see "
+                       "--iterative_solve).");
+            break;
+        }
+        case rocalution_bench_solver_parameters::itilu0_max_iter:
+        {
+            ADD_OPTION(int, e, 1000, "itilu0 maximum number of iterations.");
+            break;
+        }
+
+        case rocalution_bench_solver_parameters::itilu0_options:
+        {
+            ADD_OPTION(int, e, 1000, "itilu0 combination of options.");
+            break;
+        }
+
         case rocalution_bench_solver_parameters::mcilu_p:
         {
             ADD_OPTION(int, e, 0, "multicolored ilu parameter p.");
@@ -219,6 +250,11 @@ void rocalution_arguments_config::set_description(options_description& desc)
             ADD_OPTION(std::string, e, "", "coarsening strategy");
             break;
         }
+        case rocalution_bench_solver_parameters::itilu0_alg:
+        {
+            ADD_OPTION(std::string, e, "Default", "ItILU0 algorithm");
+            break;
+        }
         case rocalution_bench_solver_parameters::matrix:
         {
             ADD_OPTION(std::string, e, "", "matrix initialization");
@@ -243,6 +279,11 @@ void rocalution_arguments_config::set_description(options_description& desc)
     {
         switch(e)
         {
+        case rocalution_bench_solver_parameters::iterative_solve:
+        {
+            ADD_OPTION(bool, e, false, "perform triangular iterative solve during solving step.");
+            break;
+        }
         case rocalution_bench_solver_parameters::verbose:
         {
             ADD_OPTION(bool, e, false, "verbose");
@@ -335,6 +376,17 @@ int rocalution_arguments_config::parse(int& argc, char**& argv, options_descript
 
             break;
         }
+        case rocalution_bench_solver_parameters::itilu0_alg:
+        {
+
+            auto itilu0_alg_string = this->Get(rocalution_bench_solver_parameters::itilu0_alg);
+            if(itilu0_alg_string != "")
+            {
+                this->m_enum_itilu0_alg(itilu0_alg_string.c_str());
+            }
+
+            break;
+        }
 
         case rocalution_bench_solver_parameters::matrix:
         {
diff --git a/clients/benchmarks/rocalution_bench_solver_parameters.cpp b/clients/benchmarks/rocalution_bench_solver_parameters.cpp
index d48975d..e6f4268 100644
--- a/clients/benchmarks/rocalution_bench_solver_parameters.cpp
+++ b/clients/benchmarks/rocalution_bench_solver_parameters.cpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2022 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2022-2023 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -35,6 +35,10 @@ rocalution_enum_preconditioner rocalution_bench_solver_parameters::GetEnumPrecon
 {
     return this->m_enum_preconditioner;
 };
+rocalution_enum_itilu0_alg rocalution_bench_solver_parameters::GetEnumItILU0Algorithm() const
+{
+    return this->m_enum_itilu0_alg;
+};
 rocalution_enum_itsolver rocalution_bench_solver_parameters::GetEnumIterativeSolver() const
 {
     return this->m_enum_itsolver;
diff --git a/clients/benchmarks/rocalution_bench_solver_parameters.hpp b/clients/benchmarks/rocalution_bench_solver_parameters.hpp
index 67d5721..debe199 100644
--- a/clients/benchmarks/rocalution_bench_solver_parameters.hpp
+++ b/clients/benchmarks/rocalution_bench_solver_parameters.hpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2022 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2022-2023 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -25,6 +25,7 @@
 
 #include "rocalution_enum_coarsening_strategy.hpp"
 #include "rocalution_enum_directsolver.hpp"
+#include "rocalution_enum_itilu0_alg.hpp"
 #include "rocalution_enum_itsolver.hpp"
 #include "rocalution_enum_matrix_init.hpp"
 #include "rocalution_enum_preconditioner.hpp"
@@ -44,6 +45,11 @@ protected:
     rocalution_enum_matrix_init m_enum_matrix_init{};
 
     //
+    // @brief Which ItILU0 algorithm.
+    //
+    rocalution_enum_itilu0_alg m_enum_itilu0_alg{};
+
+    //
     // @brief Which iterative solver.
     //
     rocalution_enum_itsolver m_enum_itsolver{};
@@ -90,6 +96,10 @@ public:
     //
     rocalution_enum_itsolver GetEnumIterativeSolver() const;
     //
+    // @brief Get which ItILU0 algorithm.
+    //
+    rocalution_enum_itilu0_alg GetEnumItILU0Algorithm() const;
+    //
     // @brief Get which matrix initialization
     //
     rocalution_enum_matrix_init GetEnumMatrixInit() const;
@@ -100,6 +110,7 @@ public:
 
     // clang-format off
 #define PBOOL_TRANSFORM_EACH			\
+  PBOOL_TRANSFORM(iterative_solve) \
   PBOOL_TRANSFORM(verbose)			\
   PBOOL_TRANSFORM(mcilu_use_level)
     // clang-format on
@@ -122,6 +133,9 @@ public:
   PINT_TRANSFORM(krylov_basis)						\
   PINT_TRANSFORM(ndim)							\
   PINT_TRANSFORM(ilut_n)						\
+  PINT_TRANSFORM(itsolve_max_iter)      \
+  PINT_TRANSFORM(itilu0_max_iter)				\
+  PINT_TRANSFORM(itilu0_options)				\
   PINT_TRANSFORM(mcilu_p)						\
   PINT_TRANSFORM(mcilu_q)						\
   PINT_TRANSFORM(max_iter)						\
@@ -153,6 +167,7 @@ public:
   PSTRING_TRANSFORM(coarsening_strategy)				\
   PSTRING_TRANSFORM(direct_solver)					\
   PSTRING_TRANSFORM(iterative_solver)					\
+  PSTRING_TRANSFORM(itilu0_alg)					\
   PSTRING_TRANSFORM(matrix)						\
   PSTRING_TRANSFORM(matrix_filename)					\
   PSTRING_TRANSFORM(preconditioner)					\
@@ -197,6 +212,8 @@ public:
   PDOUBLE_TRANSFORM(div_tol)				\
   PDOUBLE_TRANSFORM(residual_tol)			\
   PDOUBLE_TRANSFORM(ilut_tol)				\
+  PDOUBLE_TRANSFORM(itsolve_tol)    \
+  PDOUBLE_TRANSFORM(itilu0_tol)				\
   PDOUBLE_TRANSFORM(mcgs_relax)				\
   PDOUBLE_TRANSFORM(solver_over_interp)			\
   PDOUBLE_TRANSFORM(solver_coupling_strength) \
diff --git a/clients/benchmarks/rocalution_driver_itsolver.hpp b/clients/benchmarks/rocalution_driver_itsolver.hpp
index 9bf4e36..1dc7ba4 100644
--- a/clients/benchmarks/rocalution_driver_itsolver.hpp
+++ b/clients/benchmarks/rocalution_driver_itsolver.hpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2022 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2022-2023 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -256,7 +256,7 @@ struct rocalution_driver_itsolver_default : rocalution_driver_itsolver_base<ITSO
                 AIChebyshev<rocalution::LocalMatrix<T>, rocalution::LocalVector<T>, T>;
             p->Set(3, lambda_max / 7.0, lambda_max);
             this->m_preconditioner = p;
-            return true;
+            break;
         }
 
         case rocalution_enum_preconditioner::FSAI:
@@ -264,7 +264,7 @@ struct rocalution_driver_itsolver_default : rocalution_driver_itsolver_base<ITSO
             auto* p
                 = new rocalution::FSAI<rocalution::LocalMatrix<T>, rocalution::LocalVector<T>, T>;
             this->m_preconditioner = p;
-            return true;
+            break;
         }
 
         case rocalution_enum_preconditioner::SPAI:
@@ -272,14 +272,14 @@ struct rocalution_driver_itsolver_default : rocalution_driver_itsolver_base<ITSO
             auto* p
                 = new rocalution::SPAI<rocalution::LocalMatrix<T>, rocalution::LocalVector<T>, T>;
             this->m_preconditioner = p;
-            return true;
+            break;
         }
         case rocalution_enum_preconditioner::TNS:
         {
             auto* p
                 = new rocalution::TNS<rocalution::LocalMatrix<T>, rocalution::LocalVector<T>, T>;
             this->m_preconditioner = p;
-            return true;
+            break;
         }
 
         case rocalution_enum_preconditioner::Jacobi:
@@ -290,7 +290,7 @@ struct rocalution_driver_itsolver_default : rocalution_driver_itsolver_base<ITSO
             // no specific parameters
             //
             this->m_preconditioner = p;
-            return true;
+            break;
         }
 
         case rocalution_enum_preconditioner::GS:
@@ -300,7 +300,7 @@ struct rocalution_driver_itsolver_default : rocalution_driver_itsolver_base<ITSO
             // no specific parameters
             //
             this->m_preconditioner = p;
-            return true;
+            break;
         }
 
         case rocalution_enum_preconditioner::SGS:
@@ -311,7 +311,7 @@ struct rocalution_driver_itsolver_default : rocalution_driver_itsolver_base<ITSO
             // no specific parameters
             //
             this->m_preconditioner = p;
-            return true;
+            break;
         }
 
         case rocalution_enum_preconditioner::ILU:
@@ -322,6 +322,53 @@ struct rocalution_driver_itsolver_default : rocalution_driver_itsolver_base<ITSO
             // no specific parameters
             //
             this->m_preconditioner = p;
+            break;
+        }
+        case rocalution_enum_preconditioner::ItILU0:
+        {
+            auto enum_itilu0_alg = parameters.GetEnumItILU0Algorithm();
+            if(enum_itilu0_alg.is_invalid())
+            {
+                rocalution_bench_errmsg << "enum_itilu0_alg is invalid." << std::endl;
+                return false;
+            }
+
+            auto* p
+                = new rocalution::ItILU0<rocalution::LocalMatrix<T>, rocalution::LocalVector<T>, T>;
+            p->SetMaxIter(parameters.Get(params_t::itilu0_max_iter));
+            p->SetTolerance(parameters.Get(params_t::itilu0_tol));
+            p->SetOptions(parameters.Get(params_t::itilu0_options));
+
+            switch(enum_itilu0_alg.value)
+            {
+            case rocalution_enum_itilu0_alg::Default:
+            {
+                p->SetAlgorithm(ItILU0Algorithm::Default);
+                break;
+            }
+            case rocalution_enum_itilu0_alg::AsyncInPlace:
+            {
+                p->SetAlgorithm(ItILU0Algorithm::AsyncInPlace);
+                break;
+            }
+            case rocalution_enum_itilu0_alg::AsyncSplit:
+            {
+                p->SetAlgorithm(ItILU0Algorithm::AsyncSplit);
+                break;
+            }
+            case rocalution_enum_itilu0_alg::SyncSplit:
+            {
+                p->SetAlgorithm(ItILU0Algorithm::SyncSplit);
+                break;
+            }
+            case rocalution_enum_itilu0_alg::SyncSplitFusion:
+            {
+                p->SetAlgorithm(ItILU0Algorithm::SyncSplitFusion);
+                break;
+            }
+            }
+
+            this->m_preconditioner = p;
             return true;
         }
         case rocalution_enum_preconditioner::ILUT:
@@ -331,7 +378,7 @@ struct rocalution_driver_itsolver_default : rocalution_driver_itsolver_base<ITSO
             p->Set(parameters.Get(params_t::ilut_tol), parameters.Get(params_t::ilut_n));
 
             this->m_preconditioner = p;
-            return true;
+            break;
         }
         case rocalution_enum_preconditioner::IC:
         {
@@ -340,7 +387,7 @@ struct rocalution_driver_itsolver_default : rocalution_driver_itsolver_base<ITSO
             // no specific parameters
             //
             this->m_preconditioner = p;
-            return true;
+            break;
         }
         case rocalution_enum_preconditioner::MCGS:
         {
@@ -348,7 +395,7 @@ struct rocalution_driver_itsolver_default : rocalution_driver_itsolver_base<ITSO
                 MultiColoredGS<rocalution::LocalMatrix<T>, rocalution::LocalVector<T>, T>;
             p->SetRelaxation(parameters.Get(params_t::mcgs_relax));
             this->m_preconditioner = p;
-            return true;
+            break;
         }
         case rocalution_enum_preconditioner::MCSGS:
         {
@@ -358,7 +405,7 @@ struct rocalution_driver_itsolver_default : rocalution_driver_itsolver_base<ITSO
             // no specific parameters
             //
             this->m_preconditioner = p;
-            return true;
+            break;
         }
         case rocalution_enum_preconditioner::MCILU:
         {
@@ -370,8 +417,30 @@ struct rocalution_driver_itsolver_default : rocalution_driver_itsolver_base<ITSO
                    parameters.Get(params_t::mcilu_use_level));
 
             this->m_preconditioner = p;
-            return true;
+            break;
+        }
         }
+
+        if(this->m_preconditioner != nullptr)
+        {
+            const auto itsolve = parameters.Get(params_t::iterative_solve);
+
+            rocalution::SolverDescr descr;
+
+            if(itsolve)
+            {
+                descr.SetTriSolverAlg(TriSolverAlg_Iterative);
+                descr.SetIterativeSolverMaxIteration(parameters.Get(params_t::itsolve_max_iter));
+                descr.SetIterativeSolverTolerance(parameters.Get(params_t::itsolve_tol));
+            }
+            else
+            {
+                descr.SetTriSolverAlg(TriSolverAlg_Default);
+            }
+
+            this->m_preconditioner->SetSolverDescriptor(descr);
+
+            return true;
         }
 
         return false;
diff --git a/clients/benchmarks/rocalution_driver_itsolver_uaamg.hpp b/clients/benchmarks/rocalution_driver_itsolver_uaamg.hpp
index 7d83c9c..e513a85 100644
--- a/clients/benchmarks/rocalution_driver_itsolver_uaamg.hpp
+++ b/clients/benchmarks/rocalution_driver_itsolver_uaamg.hpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2022 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2022-2023 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -181,7 +181,7 @@ struct rocalution_driver_itsolver<rocalution_enum_itsolver::uaamg, T>
         preconditioner->SetSolver(cgs);
         preconditioner->SetSmootherPreIter(pre_smooth);
         preconditioner->SetSmootherPostIter(post_smooth);
-        preconditioner->SetOperatorFormat(format);
+        preconditioner->SetOperatorFormat(format, parameters.Get(params_t::blockdim));
         preconditioner->InitMaxIter(1);
         preconditioner->Verbose(0);
 
diff --git a/clients/benchmarks/rocalution_enum_itilu0_alg.cpp b/clients/benchmarks/rocalution_enum_itilu0_alg.cpp
new file mode 100644
index 0000000..c13ba38
--- /dev/null
+++ b/clients/benchmarks/rocalution_enum_itilu0_alg.cpp
@@ -0,0 +1,99 @@
+/*! \file */
+/* ************************************************************************
+* Copyright (C) 2023 Advanced Micro Devices, Inc. All rights Reserved.
+*
+* Permission is hereby granted, free of charge, to any person obtaining a copy
+* of this software and associated documentation files (the "Software"), to deal
+* in the Software without restriction, including without limitation the rights
+* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+* copies of the Software, and to permit persons to whom the Software is
+* furnished to do so, subject to the following conditions:
+*
+* The above copyright notice and this permission notice shall be included in
+* all copies or substantial portions of the Software.
+*
+* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+* THE SOFTWARE.
+*
+* ************************************************************************ */
+#include "rocalution_enum_itilu0_alg.hpp"
+#include <iostream>
+constexpr const char* rocalution_enum_itilu0_alg::names[rocalution_enum_itilu0_alg::size];
+constexpr rocalution_enum_itilu0_alg::value_type rocalution_enum_itilu0_alg::all[];
+
+const char* rocalution_enum_itilu0_alg::to_string() const
+{
+    return rocalution_enum_itilu0_alg::to_string(this->value);
+}
+
+bool rocalution_enum_itilu0_alg::is_invalid() const
+{
+    for(auto v : all)
+    {
+        if(this->value == v)
+        {
+            return false;
+        }
+    }
+    return true;
+}
+
+rocalution_enum_itilu0_alg::rocalution_enum_itilu0_alg(const char* itilu0_alg_name)
+{
+    this->value = (value_type)-1;
+    for(auto v : all)
+    {
+        const char* str = names[v];
+        if(!strcmp(itilu0_alg_name, str))
+        {
+            this->value = v;
+            return;
+        }
+    }
+
+    rocalution_bench_errmsg << "ItILU0 algorithm '" << itilu0_alg_name
+                            << "' is invalid, the list of valid ItILU0 algorithms is" << std::endl;
+    for(auto v : all)
+    {
+        const char* str = names[v];
+        rocalution_bench_errmsg << "      - '" << str << "'" << std::endl;
+    }
+    throw false;
+}
+
+//
+// Default contructor.
+//
+rocalution_enum_itilu0_alg::rocalution_enum_itilu0_alg()
+    : value((value_type)-1){};
+
+//
+//
+//
+rocalution_enum_itilu0_alg& rocalution_enum_itilu0_alg::operator()(const char* itilu0_alg_name)
+{
+    this->value = (value_type)-1;
+    for(auto v : all)
+    {
+        const char* str = names[v];
+        if(!strcmp(itilu0_alg_name, str))
+        {
+            this->value = v;
+            return *this;
+        }
+    }
+
+    rocalution_bench_errmsg << "ItILU0 algorithm '" << itilu0_alg_name
+                            << "' is invalid, the list of valid ItILU0 algorithms is" << std::endl;
+    for(auto v : all)
+    {
+        const char* str = names[v];
+        rocalution_bench_errmsg << "      - '" << str << "'" << std::endl;
+    }
+    throw false;
+}
diff --git a/clients/benchmarks/rocalution_enum_itilu0_alg.hpp b/clients/benchmarks/rocalution_enum_itilu0_alg.hpp
new file mode 100644
index 0000000..eaa6863
--- /dev/null
+++ b/clients/benchmarks/rocalution_enum_itilu0_alg.hpp
@@ -0,0 +1,92 @@
+/*! \file */
+/* ************************************************************************
+* Copyright (C) 2023 Advanced Micro Devices, Inc. All rights Reserved.
+*
+* Permission is hereby granted, free of charge, to any person obtaining a copy
+* of this software and associated documentation files (the "Software"), to deal
+* in the Software without restriction, including without limitation the rights
+* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+* copies of the Software, and to permit persons to whom the Software is
+* furnished to do so, subject to the following conditions:
+*
+* The above copyright notice and this permission notice shall be included in
+* all copies or substantial portions of the Software.
+*
+* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+* THE SOFTWARE.
+*
+* ************************************************************************ */
+#pragma once
+#include "utility.hpp"
+#include <cstring>
+
+//
+// List the enumeration values.
+//
+
+// clang-format off
+#define ROCALUTION_ENUM_ITILU0_ALG_TRANSFORM_EACH		  \
+  ROCALUTION_ENUM_ITILU0_ALG_TRANSFORM(Default)		  \
+  ROCALUTION_ENUM_ITILU0_ALG_TRANSFORM(AsyncInPlace)	  \
+  ROCALUTION_ENUM_ITILU0_ALG_TRANSFORM(AsyncSplit)	  \
+  ROCALUTION_ENUM_ITILU0_ALG_TRANSFORM(SyncSplit)		  \
+  ROCALUTION_ENUM_ITILU0_ALG_TRANSFORM(SyncSplitFusion)
+// clang-format on
+
+struct rocalution_enum_itilu0_alg
+{
+public:
+#define ROCALUTION_ENUM_ITILU0_ALG_TRANSFORM(x_) x_,
+    typedef enum rocalution_enum_itilu0_alg__ : int
+    {
+        ROCALUTION_ENUM_ITILU0_ALG_TRANSFORM_EACH
+    } value_type;
+    static constexpr value_type all[] = {ROCALUTION_ENUM_ITILU0_ALG_TRANSFORM_EACH};
+#undef ROCALUTION_ENUM_ITILU0_ALG_TRANSFORM
+    static constexpr std::size_t size = countof(all);
+    value_type                   value{};
+
+private:
+#define ROCALUTION_ENUM_ITILU0_ALG_TRANSFORM(x_) #x_,
+    static constexpr const char* names[size]{ROCALUTION_ENUM_ITILU0_ALG_TRANSFORM_EACH};
+#undef ROCALUTION_ENUM_ITILU0_ALG_TRANSFORM
+public:
+    operator value_type() const
+    {
+        return this->value;
+    };
+    rocalution_enum_itilu0_alg();
+    rocalution_enum_itilu0_alg& operator()(const char* function);
+    rocalution_enum_itilu0_alg(const char* function);
+    const char*               to_string() const;
+    bool                      is_invalid() const;
+    static inline const char* to_string(rocalution_enum_itilu0_alg::value_type value)
+    {
+        //
+        // switch for checking inconsistency.
+        //
+        switch(value)
+        {
+            ///
+#define ROCALUTION_ENUM_ITILU0_ALG_TRANSFORM(x_) \
+    case x_:                                     \
+    {                                            \
+        if(strcmp(#x_, names[value]))            \
+            return nullptr;                      \
+        break;                                   \
+    }
+
+            ROCALUTION_ENUM_ITILU0_ALG_TRANSFORM_EACH;
+
+#undef ROCALUTION_ENUM_ITILU0_ALG_TRANSFORM
+            ///
+        }
+
+        return names[value];
+    }
+};
diff --git a/clients/benchmarks/rocalution_enum_preconditioner.hpp b/clients/benchmarks/rocalution_enum_preconditioner.hpp
index 2eaacf5..2063476 100644
--- a/clients/benchmarks/rocalution_enum_preconditioner.hpp
+++ b/clients/benchmarks/rocalution_enum_preconditioner.hpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2022 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2022-2023 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -40,6 +40,7 @@ struct rocalution_enum_preconditioner
     ENUM_PRECONDITIONER(GS)                 \
     ENUM_PRECONDITIONER(SGS)                \
     ENUM_PRECONDITIONER(ILU)                \
+    ENUM_PRECONDITIONER(ItILU0)             \
     ENUM_PRECONDITIONER(ILUT)               \
     ENUM_PRECONDITIONER(IC)                 \
     ENUM_PRECONDITIONER(MCGS)               \
diff --git a/clients/include/common.hpp b/clients/include/common.hpp
index 5cdb250..00acd2f 100644
--- a/clients/include/common.hpp
+++ b/clients/include/common.hpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2018-2022 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2018-2023 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -22,11 +22,36 @@
  * ************************************************************************ */
 
 #include <cstring>
+#include <map>
 #include <mpi.h>
 #include <rocalution/rocalution.hpp>
+#include <set>
 
 using namespace rocalution;
 
+static void my_irecv(int* buf, int count, int source, int tag, MPI_Comm comm, MPI_Request* request)
+{
+    MPI_Irecv(buf, count, MPI_INT, source, tag, comm, request);
+}
+
+static void
+    my_irecv(int64_t* buf, int count, int source, int tag, MPI_Comm comm, MPI_Request* request)
+{
+    MPI_Irecv(buf, count, MPI_INT64_T, source, tag, comm, request);
+}
+
+static void
+    my_isend(const int* buf, int count, int dest, int tag, MPI_Comm comm, MPI_Request* request)
+{
+    MPI_Isend(buf, count, MPI_INT, dest, tag, comm, request);
+}
+
+static void
+    my_isend(const int64_t* buf, int count, int dest, int tag, MPI_Comm comm, MPI_Request* request)
+{
+    MPI_Isend(buf, count, MPI_INT64_T, dest, tag, comm, request);
+}
+
 template <typename ValueType>
 void distribute_matrix(const MPI_Comm*          comm,
                        LocalMatrix<ValueType>*  lmat,
@@ -39,11 +64,11 @@ void distribute_matrix(const MPI_Comm*          comm,
     MPI_Comm_rank(*comm, &rank);
     MPI_Comm_size(*comm, &num_procs);
 
-    size_t global_nrow = lmat->GetM();
-    size_t global_ncol = lmat->GetN();
-    size_t global_nnz  = lmat->GetNnz();
+    int64_t global_nrow = lmat->GetM();
+    int64_t global_ncol = lmat->GetN();
+    int64_t global_nnz  = lmat->GetNnz();
 
-    int*       global_row_offset = NULL;
+    PtrType*   global_row_offset = NULL;
     int*       global_col        = NULL;
     ValueType* global_val        = NULL;
 
@@ -74,14 +99,14 @@ void distribute_matrix(const MPI_Comm*          comm,
 
     if(global_nrow % num_procs != 0)
     {
-        for(size_t i = 0; i < global_nrow % num_procs; ++i)
+        for(int i = 0; i < global_nrow % num_procs; ++i)
         {
             ++local_size[i];
         }
     }
 
     // Compute index offsets
-    std::vector<int> index_offset(num_procs + 1);
+    std::vector<PtrType> index_offset(num_procs + 1);
     index_offset[0] = 0;
     for(int i = 0; i < num_procs; ++i)
     {
@@ -89,10 +114,10 @@ void distribute_matrix(const MPI_Comm*          comm,
     }
 
     // Read sub matrix - row_offset
-    int              local_nrow = local_size[rank];
-    std::vector<int> local_row_offset(local_nrow + 1);
+    int                  local_nrow = local_size[rank];
+    std::vector<PtrType> local_row_offset(local_nrow + 1);
 
-    for(int i = index_offset[rank], k = 0; k < local_nrow + 1; ++i, ++k)
+    for(PtrType i = index_offset[rank], k = 0; k < local_nrow + 1; ++i, ++k)
     {
         local_row_offset[k] = global_row_offset[i];
     }
@@ -100,11 +125,11 @@ void distribute_matrix(const MPI_Comm*          comm,
     free_host(&global_row_offset);
 
     // Read sub matrix - col and val
-    int                    local_nnz = local_row_offset[local_nrow] - local_row_offset[0];
+    PtrType                local_nnz = local_row_offset[local_nrow] - local_row_offset[0];
     std::vector<int>       local_col(local_nnz);
     std::vector<ValueType> local_val(local_nnz);
 
-    for(int i = local_row_offset[0], k = 0; k < local_nnz; ++i, ++k)
+    for(PtrType i = local_row_offset[0], k = 0; k < local_nnz; ++i, ++k)
     {
         local_col[k] = global_col[i];
         local_val[k] = global_val[i];
@@ -120,18 +145,18 @@ void distribute_matrix(const MPI_Comm*          comm,
         local_row_offset[i] -= shift;
     }
 
-    int interior_nnz = 0;
-    int ghost_nnz    = 0;
-    int boundary_nnz = 0;
-    int neighbors    = 0;
+    PtrType interior_nnz = 0;
+    PtrType ghost_nnz    = 0;
+    int     boundary_nnz = 0;
+    int     neighbors    = 0;
 
-    std::vector<std::vector<int>>    boundary(num_procs, std::vector<int>());
-    std::vector<bool>                neighbor(num_procs, false);
-    std::vector<std::map<int, bool>> checked(num_procs, std::map<int, bool>());
+    std::vector<std::vector<PtrType>> boundary(num_procs, std::vector<PtrType>());
+    std::vector<bool>                 neighbor(num_procs, false);
+    std::vector<std::map<int, bool>>  checked(num_procs, std::map<int, bool>());
 
     for(int i = 0; i < local_nrow; ++i)
     {
-        for(int j = local_row_offset[i]; j < local_row_offset[i + 1]; ++j)
+        for(PtrType j = local_row_offset[i]; j < local_row_offset[i + 1]; ++j)
         {
 
             // Interior point
@@ -214,7 +239,7 @@ void distribute_matrix(const MPI_Comm*          comm,
         if(neighbor[i] == true)
         {
             // Receive size of boundary from rank i to current rank
-            MPI_Irecv(&(boundary_size[n]), 1, MPI_INT, i, 0, *comm, &mpi_req[n]);
+            my_irecv(&(boundary_size[n]), 1, i, 0, *comm, &mpi_req[n]);
             ++n;
         }
     }
@@ -227,7 +252,7 @@ void distribute_matrix(const MPI_Comm*          comm,
         {
             int size = boundary[i].size();
             // Send size of boundary from current rank to rank i
-            MPI_Isend(&size, 1, MPI_INT, i, 0, *comm, &mpi_req[n]);
+            my_isend(&size, 1, i, 0, *comm, &mpi_req[n]);
             ++n;
         }
     }
@@ -256,7 +281,7 @@ void distribute_matrix(const MPI_Comm*          comm,
     }
 
     // Array to hold boundary for each interface
-    std::vector<std::vector<int>> local_boundary(neighbors);
+    std::vector<std::vector<PtrType>> local_boundary(neighbors);
     for(int i = 0; i < neighbors; ++i)
     {
         local_boundary[i].resize(boundary_size[i]);
@@ -269,8 +294,7 @@ void distribute_matrix(const MPI_Comm*          comm,
         if(neighbor[i] == true)
         {
             // Receive boundary from rank i to current rank
-            MPI_Irecv(
-                local_boundary[n].data(), boundary_size[n], MPI_INT, i, 0, *comm, &mpi_req[n]);
+            my_irecv(local_boundary[n].data(), boundary_size[n], i, 0, *comm, &mpi_req[n]);
             ++n;
         }
     }
@@ -282,7 +306,7 @@ void distribute_matrix(const MPI_Comm*          comm,
         if(boundary[i].size() > 0)
         {
             // Send boundary from current rank to rank i
-            MPI_Isend(&(boundary[i][0]), boundary[i].size(), MPI_INT, i, 0, *comm, &mpi_req[n]);
+            my_isend(&(boundary[i][0]), boundary[i].size(), i, 0, *comm, &mpi_req[n]);
             ++n;
         }
     }
@@ -305,13 +329,13 @@ void distribute_matrix(const MPI_Comm*          comm,
     {
         for(unsigned int j = 0; j < boundary[i].size(); ++j)
         {
-            bnd[k] = boundary[i][j] - index_offset[rank];
+            bnd[k] = static_cast<int>(boundary[i][j] - index_offset[rank]);
             ++k;
         }
     }
 
     // Create boundary index array
-    std::vector<int> boundary_index(nnz_boundary);
+    std::vector<PtrType> boundary_index(nnz_boundary);
 
     k = 0;
     for(int i = 0; i < neighbors; ++i)
@@ -340,11 +364,11 @@ void distribute_matrix(const MPI_Comm*          comm,
     memset(ghost_col, 0, sizeof(int) * ghost_nnz);
     memset(ghost_val, 0, sizeof(ValueType) * ghost_nnz);
 
-    int*       row_offset = new int[local_nrow + 1];
+    PtrType*   row_offset = new PtrType[local_nrow + 1];
     int*       col        = new int[interior_nnz];
     ValueType* val        = new ValueType[interior_nnz];
 
-    memset(row_offset, 0, sizeof(int) * (local_nrow + 1));
+    memset(row_offset, 0, sizeof(PtrType) * (local_nrow + 1));
     memset(col, 0, sizeof(int) * interior_nnz);
     memset(val, 0, sizeof(ValueType) * interior_nnz);
 
@@ -353,7 +377,7 @@ void distribute_matrix(const MPI_Comm*          comm,
     int l         = 0;
     for(int i = 0; i < local_nrow; ++i)
     {
-        for(int j = local_row_offset[i]; j < local_row_offset[i + 1]; ++j)
+        for(PtrType j = local_row_offset[i]; j < local_row_offset[i + 1]; ++j)
         {
 
             // Boundary point -- create ghost part
@@ -403,4 +427,824 @@ void distribute_matrix(const MPI_Comm*          comm,
     gmat->SetParallelManager(*pm);
     gmat->SetLocalDataPtrCSR(&row_offset, &col, &val, "mat", interior_nnz);
     gmat->SetGhostDataPtrCOO(&ghost_row, &ghost_col, &ghost_val, "ghost", ghost_nnz);
+    gmat->Sort();
+}
+
+// This function computes all prime factors of a given number n
+static void compute_prime_factors(int n, std::vector<int>& p)
+{
+    int factor = 2;
+
+    // Factorize
+    while(n > 1)
+    {
+        while(n % factor == 0)
+        {
+            p.push_back(factor);
+            n /= factor;
+        }
+
+        ++factor;
+    }
+}
+
+// This function computes the process distribution for each dimension
+static void compute_2d_process_distribution(int nprocs, int& nprocx, int& nprocy)
+{
+    // Compute prime factors
+    std::vector<int> p;
+    compute_prime_factors(nprocs, p);
+
+    // Compute number of processes in each dimension
+    nprocx = 1;
+    nprocy = 1;
+
+    if(p.size() == 0)
+    {
+        // No entry, this means we have exactly one process
+    }
+    else if(p.size() == 1)
+    {
+        // If we have a single prime number, this is going to be our x dimension
+        nprocx = p[0];
+    }
+    else if(p.size() == 2)
+    {
+        // For two prime numbers, setup x and y
+        nprocx = p[1];
+        nprocy = p[0];
+    }
+    else
+    {
+        // More than two prime numbers
+
+        // #prime numbers
+        int    idx    = 0;
+        size_t nprime = p.size();
+
+        // cubic root
+        double sqroot = std::sqrt(nprocs);
+
+        // Determine x dimension
+        nprocx = p[nprime-- - 1];
+
+        while(nprocx < sqroot && idx < nprime)
+        {
+            nprocx *= p[idx++];
+        }
+
+        // Determine y dimension
+        while(idx < nprime)
+        {
+            nprocy *= p[idx++];
+        }
+    }
+
+    // Number of processes must match
+    assert(nprocx * nprocy == nprocs);
+}
+
+// This function computes the process distribution for each dimension
+static void compute_3d_process_distribution(int nprocs, int& nprocx, int& nprocy, int& nprocz)
+{
+    // Compute prime factors
+    std::vector<int> p;
+    compute_prime_factors(nprocs, p);
+
+    // Compute number of processes in each dimension
+    nprocx = 1;
+    nprocy = 1;
+    nprocz = 1;
+
+    if(p.size() == 0)
+    {
+        // No entry, this means we have exactly one process
+    }
+    else if(p.size() == 1)
+    {
+        // If we have a single prime number, this is going to be our x dimension
+        nprocx = p[0];
+    }
+    else if(p.size() == 2)
+    {
+        // For two prime numbers, setup x and y
+        nprocx = p[1];
+        nprocy = p[0];
+    }
+    else if(p.size() == 3)
+    {
+        // Three prime numbers
+        nprocx = p[2];
+        nprocy = p[1];
+        nprocz = p[0];
+    }
+    else
+    {
+        // More than three prime numbers
+
+        // #prime numbers
+        int    idx    = 0;
+        size_t nprime = p.size();
+
+        // cubic root
+        double qroot = std::cbrt(nprocs);
+
+        // Determine x dimension
+        nprocx = p[nprime-- - 1];
+
+        while(nprocx < qroot && idx < nprime)
+        {
+            nprocx *= p[idx++];
+        }
+
+        // Determine y dimension
+        double sqroot = std::sqrt(nprocs / nprocx);
+
+        nprocy = p[nprime-- - 1];
+
+        while(nprocy < sqroot && idx < nprime)
+        {
+            nprocy *= p[idx++];
+        }
+
+        // Determine z dimension
+        while(idx < nprime)
+        {
+            nprocz *= p[idx++];
+        }
+    }
+
+    // Number of processes must match
+    assert(nprocx * nprocy * nprocz == nprocs);
+}
+
+template <typename ValueType>
+void generate_2d_laplacian(int                      local_dimx,
+                           int                      local_dimy,
+                           const MPI_Comm*          comm,
+                           GlobalMatrix<ValueType>* mat,
+                           ParallelManager*         pm,
+                           int                      rank,
+                           int                      nprocs,
+                           int                      nsten = 9)
+{
+    assert(nsten == 5 || nsten == 9);
+
+    // First, we need to determine process pattern for the unit square
+    int nproc_x;
+    int nproc_y;
+
+    compute_2d_process_distribution(nprocs, nproc_x, nproc_y);
+
+    // Next, determine process index into the unit square
+    int iproc_y = rank / nproc_x;
+    int iproc_x = rank % nproc_x;
+
+    // Global sizes
+    int64_t global_dimx = static_cast<int64_t>(nproc_x) * local_dimx;
+    int64_t global_dimy = static_cast<int64_t>(nproc_y) * local_dimy;
+
+    // Global process entry points
+    int64_t global_iproc_x = iproc_x * local_dimx;
+    int64_t global_iproc_y = iproc_y * local_dimy;
+
+    // Number of rows (global and local)
+    int64_t local_nrow  = local_dimx * local_dimy;
+    int64_t global_nrow = global_dimx * global_dimy;
+
+    // Assemble local CSR matrix row offset pointers
+    PtrType* global_csr_row_ptr = NULL;
+    int64_t* global_csr_col_ind = NULL;
+    int64_t* local2global       = NULL;
+
+    allocate_host(local_nrow + 1, &global_csr_row_ptr);
+    allocate_host(local_nrow * nsten, &global_csr_col_ind);
+    allocate_host(local_nrow, &local2global);
+
+    std::map<int64_t, int> global2local;
+
+    PtrType nnz           = 0;
+    global_csr_row_ptr[0] = 0;
+
+    // Loop over y dimension
+    for(int local_y = 0; local_y < local_dimy; ++local_y)
+    {
+        // Global index into y
+        int64_t global_y = global_iproc_y + local_y;
+
+        // Loop over x dimension
+        for(int local_x = 0; local_x < local_dimx; ++local_x)
+        {
+            // Global index into x
+            int64_t global_x = global_iproc_x + local_x;
+
+            // Local row
+            int local_row = local_y * local_dimx + local_x;
+
+            // Global row
+            int64_t global_row = global_y * global_dimx + global_x;
+
+            // Fill l2g and g2l map
+            local2global[local_row]  = global_row;
+            global2local[global_row] = local_row;
+
+            // 5pt stencil
+            if(nsten == 5)
+            {
+                // Fixed x (leaving out i == j)
+                for(int by = -1; by <= 1; ++by)
+                {
+                    if(global_y + by > -1 && global_y + by < global_dimy && by != 0)
+                    {
+                        // Global column
+                        int64_t global_col = global_row + by * global_dimx;
+
+                        // Fill global CSR column indices
+                        global_csr_col_ind[nnz++] = global_col;
+                    }
+                }
+
+                // Fixed y
+                for(int bx = -1; bx <= 1; ++bx)
+                {
+                    if(global_x + bx > -1 && global_x + bx < global_dimx)
+                    {
+                        // Global column
+                        int64_t global_col = global_row + bx;
+
+                        // Fill global CSR column indices
+                        global_csr_col_ind[nnz++] = global_col;
+                    }
+                }
+            }
+
+            // 9 pt stencil
+            if(nsten == 9)
+            {
+                // Check if current y vertex is on the boundary
+                for(int by = -1; by <= 1; ++by)
+                {
+                    if(global_y + by > -1 && global_y + by < global_dimy)
+                    {
+                        // Check if current x vertex is on the boundary
+                        for(int bx = -1; bx <= 1; ++bx)
+                        {
+                            if(global_x + bx > -1 && global_x + bx < global_dimx)
+                            {
+                                // Global column
+                                int64_t global_col = global_row + by * global_dimx + bx;
+
+                                // Fill global CSR column indices
+                                global_csr_col_ind[nnz++] = global_col;
+                            }
+                        }
+                    }
+                }
+            }
+
+            global_csr_row_ptr[local_row + 1] = nnz;
+        }
+    }
+
+    // Local number of non-zero entries - need to use long long int to make the communication work
+    int64_t local_nnz = global_csr_row_ptr[local_nrow];
+
+    // Total number of non-zeros
+    int64_t global_nnz;
+    MPI_Allreduce(&local_nnz, &global_nnz, 1, MPI_INT64_T, MPI_SUM, *comm);
+
+    // Now, we need to setup the communication pattern
+    std::map<int, std::set<int64_t>> recv_indices;
+    std::map<int, std::set<int64_t>> send_indices;
+
+    // CSR matrix row pointers
+    PtrType* int_csr_row_ptr = NULL;
+    PtrType* gst_csr_row_ptr = NULL;
+
+    allocate_host(local_nrow + 1, &int_csr_row_ptr);
+    allocate_host(local_nrow + 1, &gst_csr_row_ptr);
+
+    int_csr_row_ptr[0] = 0;
+    gst_csr_row_ptr[0] = 0;
+
+    // Determine, which vertices need to be sent / received
+    for(int i = 0; i < local_nrow; ++i)
+    {
+        int_csr_row_ptr[i + 1] = int_csr_row_ptr[i];
+        gst_csr_row_ptr[i + 1] = gst_csr_row_ptr[i];
+
+        int64_t global_row = local2global[i];
+
+        for(PtrType j = global_csr_row_ptr[i]; j < global_csr_row_ptr[i + 1]; ++j)
+        {
+            int64_t global_col = global_csr_col_ind[j];
+
+            // Determine which process owns the vertex
+            int64_t idx_y = global_col / global_dimx;
+            int64_t idx_x = global_col % global_dimx;
+
+            int idx_proc_y = idx_y / local_dimy;
+            int idx_proc_x = idx_x / local_dimx;
+
+            int owner = idx_proc_x + idx_proc_y * nproc_x;
+
+            // If we do not own it, we need to receive it from our neighbor
+            // and also send the current vertex to this neighbor
+            if(owner != rank)
+            {
+                // Store the global column and row id that we have to receive / send from / to a neighbor
+                // We need a set here to eliminate duplicates
+                recv_indices[owner].insert(global_col);
+                send_indices[owner].insert(global_row);
+
+                ++gst_csr_row_ptr[i + 1];
+            }
+            else
+            {
+                ++int_csr_row_ptr[i + 1];
+            }
+        }
+    }
+
+    // Number of processes we communicate with
+    int nrecv = recv_indices.size();
+    int nsend = send_indices.size();
+
+    // Process ids we communicate with
+    std::vector<int> recvs;
+    std::vector<int> sends;
+
+    recvs.reserve(nrecv);
+    sends.reserve(nsend);
+
+    // Index offsets for each neighbor
+    std::vector<int> recv_index_offset;
+    std::vector<int> send_index_offset;
+
+    recv_index_offset.reserve(nrecv + 1);
+    send_index_offset.reserve(nsend + 1);
+
+    recv_index_offset.push_back(0);
+    send_index_offset.push_back(0);
+
+    int                    cnt = 0;
+    std::map<int64_t, int> global2ghost;
+
+    // Go through the recv data
+    for(std::map<int, std::set<int64_t>>::iterator it = recv_indices.begin();
+        it != recv_indices.end();
+        ++it)
+    {
+        recvs.push_back(it->first);
+        recv_index_offset.push_back(it->second.size());
+
+        for(std::set<int64_t>::iterator iit = it->second.begin(); iit != it->second.end(); ++iit)
+        {
+            global2ghost[*iit] = cnt++;
+        }
+    }
+
+    // Go through the send data
+    int boundary_nnz = 0;
+    for(std::map<int, std::set<int64_t>>::iterator it = send_indices.begin();
+        it != send_indices.end();
+        ++it)
+    {
+        sends.push_back(it->first);
+        send_index_offset.push_back(it->second.size());
+        boundary_nnz += it->second.size();
+    }
+
+    // Exclusive sum
+    for(int i = 0; i < nrecv; ++i)
+    {
+        recv_index_offset[i + 1] += recv_index_offset[i];
+    }
+
+    for(int i = 0; i < nsend; ++i)
+    {
+        send_index_offset[i + 1] += send_index_offset[i];
+    }
+
+    // Boundary indices
+    std::vector<int> boundary;
+    boundary.reserve(boundary_nnz);
+
+    for(std::map<int, std::set<int64_t>>::iterator it = send_indices.begin();
+        it != send_indices.end();
+        ++it)
+    {
+        for(std::set<int64_t>::iterator iit = it->second.begin(); iit != it->second.end(); ++iit)
+        {
+            boundary.push_back(global2local[*iit]);
+        }
+    }
+
+    // Initialize manager
+    pm->SetMPICommunicator(comm);
+    pm->SetGlobalNrow(global_nrow);
+    pm->SetGlobalNcol(global_nrow);
+    pm->SetLocalNrow(local_nrow);
+    pm->SetLocalNcol(local_nrow);
+
+    if(nprocs > 1)
+    {
+        pm->SetBoundaryIndex(boundary_nnz, boundary.data());
+        pm->SetReceivers(nrecv, recvs.data(), recv_index_offset.data());
+        pm->SetSenders(nsend, sends.data(), send_index_offset.data());
+    }
+
+    mat->SetParallelManager(*pm);
+
+    // Generate local and ghost matrices
+    local_nnz         = int_csr_row_ptr[local_nrow];
+    int64_t ghost_nnz = gst_csr_row_ptr[local_nrow];
+
+    int*       int_csr_col_ind = NULL;
+    int*       gst_csr_col_ind = NULL;
+    ValueType* int_csr_val     = NULL;
+    ValueType* gst_csr_val     = NULL;
+
+    allocate_host(local_nnz, &int_csr_col_ind);
+    allocate_host(local_nnz, &int_csr_val);
+    allocate_host(ghost_nnz, &gst_csr_col_ind);
+    allocate_host(ghost_nnz, &gst_csr_val);
+
+    // Convert global matrix columns to local columns
+    for(int i = 0; i < local_nrow; ++i)
+    {
+        PtrType local_idx = int_csr_row_ptr[i];
+        PtrType ghost_idx = gst_csr_row_ptr[i];
+
+        int64_t global_row = local2global[i];
+
+        for(PtrType j = global_csr_row_ptr[i]; j < global_csr_row_ptr[i + 1]; ++j)
+        {
+            int64_t global_col = global_csr_col_ind[j];
+
+            // Determine which process owns the vertex
+            int64_t idx_y = global_col / global_dimx;
+            int64_t idx_x = global_col % global_dimx;
+
+            int idx_proc_y = idx_y / local_dimy;
+            int idx_proc_x = idx_x / local_dimx;
+
+            int owner = idx_proc_x + idx_proc_y * nproc_x;
+
+            // If we do not own it, we need to receive it from our neighbor
+            // and also send the current vertex to this neighbor
+            if(owner != rank)
+            {
+                // Store the global column and row id that we have to receive / send from / to a neighbor
+                // We need a set here to eliminate duplicates
+                recv_indices[owner].insert(global_col);
+                send_indices[owner].insert(global_row);
+
+                gst_csr_col_ind[ghost_idx] = global2ghost[global_col];
+                gst_csr_val[ghost_idx]     = -1.0;
+                ++ghost_idx;
+            }
+            else
+            {
+                // This is our part
+                int_csr_col_ind[local_idx] = global2local[global_col];
+                int_csr_val[local_idx]     = (global_col == global_row) ? (nsten - 1.0) : -1.0;
+                ++local_idx;
+            }
+        }
+    }
+
+    free_host(&global_csr_row_ptr);
+    free_host(&global_csr_col_ind);
+    free_host(&local2global);
+
+    mat->SetLocalDataPtrCSR(&int_csr_row_ptr, &int_csr_col_ind, &int_csr_val, "mat", local_nnz);
+    mat->SetGhostDataPtrCSR(&gst_csr_row_ptr, &gst_csr_col_ind, &gst_csr_val, "gst", ghost_nnz);
+    mat->Sort();
+}
+
+template <typename ValueType>
+void generate_3d_laplacian(int                      local_dimx,
+                           int                      local_dimy,
+                           int                      local_dimz,
+                           const MPI_Comm*          comm,
+                           GlobalMatrix<ValueType>* mat,
+                           ParallelManager*         pm,
+                           int                      rank,
+                           int                      nprocs)
+{
+    // First, we need to determine process pattern for the unit cube
+    int nproc_x;
+    int nproc_y;
+    int nproc_z;
+
+    compute_3d_process_distribution(nprocs, nproc_x, nproc_y, nproc_z);
+
+    // Next, determine process index into the unit cube
+    int iproc_z = rank / (nproc_x * nproc_y);
+    int iproc_y = (rank - iproc_z * nproc_x * nproc_y) / nproc_x;
+    int iproc_x = rank % nproc_x;
+
+    // Global sizes
+    int64_t global_dimx = static_cast<int64_t>(nproc_x) * local_dimx;
+    int64_t global_dimy = static_cast<int64_t>(nproc_y) * local_dimy;
+    int64_t global_dimz = static_cast<int64_t>(nproc_z) * local_dimz;
+
+    // Global process entry points
+    int64_t global_iproc_x = iproc_x * local_dimx;
+    int64_t global_iproc_y = iproc_y * local_dimy;
+    int64_t global_iproc_z = iproc_z * local_dimz;
+
+    // Number of rows (global and local)
+    int64_t local_nrow  = local_dimx * local_dimy * local_dimz;
+    int64_t global_nrow = global_dimx * global_dimy * global_dimz;
+
+    // Assemble local CSR matrix row offset pointers
+    std::vector<PtrType> global_csr_row_ptr(local_nrow + 1);
+    std::vector<int64_t> global_csr_col_ind(local_nrow * 27);
+
+    std::vector<int64_t>   local2global(local_nrow);
+    std::map<int64_t, int> global2local;
+
+    PtrType nnz           = 0;
+    global_csr_row_ptr[0] = 0;
+
+    // Loop over z dimension
+    for(int local_z = 0; local_z < local_dimz; ++local_z)
+    {
+        // Global index into z
+        int64_t global_z = global_iproc_z + local_z;
+
+        // Loop over y dimension
+        for(int local_y = 0; local_y < local_dimy; ++local_y)
+        {
+            // Global index into y
+            int64_t global_y = global_iproc_y + local_y;
+
+            // Loop over x dimension
+            for(int local_x = 0; local_x < local_dimx; ++local_x)
+            {
+                // Global index into x
+                int64_t global_x = global_iproc_x + local_x;
+
+                // Local row
+                int local_row = local_z * local_dimx * local_dimy + local_y * local_dimx + local_x;
+
+                // Global row
+                int64_t global_row
+                    = global_z * global_dimx * global_dimy + global_y * global_dimx + global_x;
+
+                // Fill l2g and g2l map
+                local2global[local_row]  = global_row;
+                global2local[global_row] = local_row;
+
+                // Check if current z vertex is on the boundary
+                for(int bz = -1; bz <= 1; ++bz)
+                {
+                    if(global_z + bz > -1 && global_z + bz < global_dimz)
+                    {
+                        // Check if current y vertex is on the boundary
+                        for(int by = -1; by <= 1; ++by)
+                        {
+                            if(global_y + by > -1 && global_y + by < global_dimy)
+                            {
+                                // Check if current x vertex is on the boundary
+                                for(int bx = -1; bx <= 1; ++bx)
+                                {
+                                    if(global_x + bx > -1 && global_x + bx < global_dimx)
+                                    {
+                                        // Global column
+                                        int64_t global_col = global_row
+                                                             + bz * global_dimx * global_dimy
+                                                             + by * global_dimx + bx;
+
+                                        // Fill global CSR column indices
+                                        global_csr_col_ind[nnz++] = global_col;
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+
+                global_csr_row_ptr[local_row + 1] = nnz;
+            }
+        }
+    }
+
+    // Local number of non-zero entries - need to use long long int to make the communication work
+    int64_t local_nnz = global_csr_row_ptr[local_nrow];
+
+    // Total number of non-zeros
+    int64_t global_nnz;
+    MPI_Allreduce(&local_nnz, &global_nnz, 1, MPI_INT64_T, MPI_SUM, *comm);
+
+    // Now, we need to setup the communication pattern
+    std::map<int, std::set<int64_t>> recv_indices;
+    std::map<int, std::set<int64_t>> send_indices;
+
+    // CSR matrix row pointers
+    PtrType* int_csr_row_ptr = NULL;
+    PtrType* gst_csr_row_ptr = NULL;
+
+    allocate_host(local_nrow + 1, &int_csr_row_ptr);
+    allocate_host(local_nrow + 1, &gst_csr_row_ptr);
+
+    int_csr_row_ptr[0] = 0;
+    gst_csr_row_ptr[0] = 0;
+
+    // Determine, which vertices need to be sent / received
+    for(int i = 0; i < local_nrow; ++i)
+    {
+        int_csr_row_ptr[i + 1] = int_csr_row_ptr[i];
+        gst_csr_row_ptr[i + 1] = gst_csr_row_ptr[i];
+
+        int64_t global_row = local2global[i];
+
+        for(PtrType j = global_csr_row_ptr[i]; j < global_csr_row_ptr[i + 1]; ++j)
+        {
+            int64_t global_col = global_csr_col_ind[j];
+
+            // Determine which process owns the vertex
+            int64_t idx_z = global_col / (global_dimx * global_dimy);
+            int64_t idx_y = (global_col - idx_z * global_dimy * global_dimx) / global_dimx;
+            int64_t idx_x = global_col % global_dimx;
+
+            int idx_proc_z = idx_z / local_dimz;
+            int idx_proc_y = idx_y / local_dimy;
+            int idx_proc_x = idx_x / local_dimx;
+
+            int owner = idx_proc_x + idx_proc_y * nproc_x + idx_proc_z * nproc_y * nproc_x;
+
+            // If we do not own it, we need to receive it from our neighbor
+            // and also send the current vertex to this neighbor
+            if(owner != rank)
+            {
+                // Store the global column and row id that we have to receive / send from / to a neighbor
+                // We need a set here to eliminate duplicates
+                recv_indices[owner].insert(global_col);
+                send_indices[owner].insert(global_row);
+
+                ++gst_csr_row_ptr[i + 1];
+            }
+            else
+            {
+                ++int_csr_row_ptr[i + 1];
+            }
+        }
+    }
+
+    // Number of processes we communicate with
+    int nrecv = recv_indices.size();
+    int nsend = send_indices.size();
+
+    // Process ids we communicate with
+    std::vector<int> recvs;
+    std::vector<int> sends;
+
+    recvs.reserve(nrecv);
+    sends.reserve(nsend);
+
+    // Index offsets for each neighbor
+    std::vector<int> recv_index_offset;
+    std::vector<int> send_index_offset;
+
+    recv_index_offset.reserve(nrecv + 1);
+    send_index_offset.reserve(nsend + 1);
+
+    recv_index_offset.push_back(0);
+    send_index_offset.push_back(0);
+
+    int                    cnt = 0;
+    std::map<int64_t, int> global2ghost;
+
+    // Go through the recv data
+    for(std::map<int, std::set<int64_t>>::iterator it = recv_indices.begin();
+        it != recv_indices.end();
+        ++it)
+    {
+        recvs.push_back(it->first);
+        recv_index_offset.push_back(it->second.size());
+
+        for(std::set<int64_t>::iterator iit = it->second.begin(); iit != it->second.end(); ++iit)
+        {
+            global2ghost[*iit] = cnt++;
+        }
+    }
+
+    // Go through the send data
+    int boundary_nnz = 0;
+    for(std::map<int, std::set<int64_t>>::iterator it = send_indices.begin();
+        it != send_indices.end();
+        ++it)
+    {
+        sends.push_back(it->first);
+        send_index_offset.push_back(it->second.size());
+        boundary_nnz += it->second.size();
+    }
+
+    // Exclusive sum
+    for(int i = 0; i < nrecv; ++i)
+    {
+        recv_index_offset[i + 1] += recv_index_offset[i];
+    }
+
+    for(int i = 0; i < nsend; ++i)
+    {
+        send_index_offset[i + 1] += send_index_offset[i];
+    }
+
+    // Boundary indices
+    std::vector<int> boundary;
+    boundary.reserve(boundary_nnz);
+
+    for(std::map<int, std::set<int64_t>>::iterator it = send_indices.begin();
+        it != send_indices.end();
+        ++it)
+    {
+        for(std::set<int64_t>::iterator iit = it->second.begin(); iit != it->second.end(); ++iit)
+        {
+            boundary.push_back(global2local[*iit]);
+        }
+    }
+
+    // Initialize manager
+    pm->SetMPICommunicator(comm);
+    pm->SetGlobalNrow(global_nrow);
+    pm->SetGlobalNcol(global_nrow);
+    pm->SetLocalNrow(local_nrow);
+    pm->SetLocalNcol(local_nrow);
+
+    if(nprocs > 1)
+    {
+        pm->SetBoundaryIndex(boundary_nnz, boundary.data());
+        pm->SetReceivers(nrecv, recvs.data(), recv_index_offset.data());
+        pm->SetSenders(nsend, sends.data(), send_index_offset.data());
+    }
+
+    mat->SetParallelManager(*pm);
+
+    // Generate local and ghost matrices
+    local_nnz         = int_csr_row_ptr[local_nrow];
+    int64_t ghost_nnz = gst_csr_row_ptr[local_nrow];
+
+    int*       int_csr_col_ind = NULL;
+    int*       gst_csr_col_ind = NULL;
+    ValueType* int_csr_val     = NULL;
+    ValueType* gst_csr_val     = NULL;
+
+    allocate_host(local_nnz, &int_csr_col_ind);
+    allocate_host(local_nnz, &int_csr_val);
+    allocate_host(ghost_nnz, &gst_csr_col_ind);
+    allocate_host(ghost_nnz, &gst_csr_val);
+
+    // Convert global matrix columns to local columns
+    for(int i = 0; i < local_nrow; ++i)
+    {
+        PtrType local_idx = int_csr_row_ptr[i];
+        PtrType ghost_idx = gst_csr_row_ptr[i];
+
+        int64_t global_row = local2global[i];
+
+        for(PtrType j = global_csr_row_ptr[i]; j < global_csr_row_ptr[i + 1]; ++j)
+        {
+            int64_t global_col = global_csr_col_ind[j];
+
+            // Determine which process owns the vertex
+            int64_t idx_z = global_col / (global_dimx * global_dimy);
+            int64_t idx_y = (global_col - idx_z * global_dimy * global_dimx) / global_dimx;
+            int64_t idx_x = global_col % global_dimx;
+
+            int idx_proc_z = idx_z / local_dimz;
+            int idx_proc_y = idx_y / local_dimy;
+            int idx_proc_x = idx_x / local_dimx;
+
+            int owner = idx_proc_x + idx_proc_y * nproc_x + idx_proc_z * nproc_y * nproc_x;
+
+            // If we do not own it, we need to receive it from our neighbor
+            // and also send the current vertex to this neighbor
+            if(owner != rank)
+            {
+                // Store the global column and row id that we have to receive / send from / to a neighbor
+                // We need a set here to eliminate duplicates
+                recv_indices[owner].insert(global_col);
+                send_indices[owner].insert(global_row);
+
+                gst_csr_col_ind[ghost_idx] = global2ghost[global_col];
+                gst_csr_val[ghost_idx]     = -1.0;
+                ++ghost_idx;
+            }
+            else
+            {
+                // This is our part
+                int_csr_col_ind[local_idx] = global2local[global_col];
+                int_csr_val[local_idx]     = (global_col == global_row) ? 26.0 : -1.0;
+                ++local_idx;
+            }
+        }
+    }
+
+    mat->SetLocalDataPtrCSR(&int_csr_row_ptr, &int_csr_col_ind, &int_csr_val, "mat", local_nnz);
+    mat->SetGhostDataPtrCSR(&gst_csr_row_ptr, &gst_csr_col_ind, &gst_csr_val, "gst", ghost_nnz);
+    mat->Sort();
 }
diff --git a/clients/include/testing_bicgstab.hpp b/clients/include/testing_bicgstab.hpp
index 924295c..dbee12e 100644
--- a/clients/include/testing_bicgstab.hpp
+++ b/clients/include/testing_bicgstab.hpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2018-2022 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2018-2025 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -44,11 +44,13 @@ static bool check_residual(double res)
 template <typename T>
 bool testing_bicgstab(Arguments argus)
 {
-    int          ndim    = argus.size;
-    std::string  precond = argus.precond;
-    unsigned int format  = argus.format;
+    int          ndim                = argus.size;
+    std::string  precond             = argus.precond;
+    unsigned int format              = argus.format;
+    bool         disable_accelerator = !argus.use_acc;
 
     // Initialize rocALUTION platform
+    disable_accelerator_rocalution(disable_accelerator);
     set_device_rocalution(device);
     init_rocalution();
 
@@ -69,10 +71,13 @@ bool testing_bicgstab(Arguments argus)
     A.SetDataPtrCSR(&csr_ptr, &csr_col, &csr_val, "A", nnz, nrow, nrow);
 
     // Move data to accelerator
-    A.MoveToAccelerator();
-    x.MoveToAccelerator();
-    b.MoveToAccelerator();
-    e.MoveToAccelerator();
+    if(!disable_accelerator)
+    {
+        A.MoveToAccelerator();
+        x.MoveToAccelerator();
+        b.MoveToAccelerator();
+        e.MoveToAccelerator();
+    }
 
     // Allocate x, b and e
     x.Allocate("x", A.GetN());
@@ -124,6 +129,8 @@ bool testing_bicgstab(Arguments argus)
         p = new SGS<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "ILU")
         p = new ILU<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "ItILU0")
+        p = new ItILU0<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "ILUT")
         p = new ILUT<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "IC")
@@ -146,14 +153,20 @@ bool testing_bicgstab(Arguments argus)
         ls.SetPreconditioner(*p);
     }
 
-    ls.Init(1e-8, 0.0, 1e+8, 10000);
+    ls.Init(1e-8, 0.0, 1e+8, 0, 10000);
+    ls.RecordResidualHistory();
+    auto n_iter = ls.GetIterationCount();
     ls.Build();
 
     // Matrix format
-    A.ConvertTo(format, format == BCSR ? 3 : 1);
+    A.ConvertTo(format, format == BCSR ? argus.blockdim : 1);
 
     ls.Solve(b, &x);
 
+    const std::string filename = get_temp_dir() + "test_recorded_history.txt";
+    ls.RecordHistory(filename);
+    std::remove(filename.c_str());
+
     // Verify solution
     x.ScaleAdd(-1.0, e);
     T nrm2 = x.Norm();
@@ -169,6 +182,7 @@ bool testing_bicgstab(Arguments argus)
 
     // Stop rocALUTION platform
     stop_rocalution();
+    disable_accelerator_rocalution(false);
 
     return success;
 }
diff --git a/clients/include/testing_bicgstabl.hpp b/clients/include/testing_bicgstabl.hpp
index 20f7c5d..dc4dfb3 100644
--- a/clients/include/testing_bicgstabl.hpp
+++ b/clients/include/testing_bicgstabl.hpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2018-2022 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2018-2025 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -44,12 +44,14 @@ static bool check_residual(double res)
 template <typename T>
 bool testing_bicgstabl(Arguments argus)
 {
-    int          ndim    = argus.size;
-    std::string  precond = argus.precond;
-    unsigned int format  = argus.format;
-    int          l       = argus.index;
+    int          ndim                = argus.size;
+    std::string  precond             = argus.precond;
+    unsigned int format              = argus.format;
+    int          l                   = argus.index;
+    bool         disable_accelerator = !argus.use_acc;
 
     // Initialize rocALUTION platform
+    disable_accelerator_rocalution(disable_accelerator);
     set_device_rocalution(device);
     init_rocalution();
 
@@ -70,10 +72,13 @@ bool testing_bicgstabl(Arguments argus)
     A.SetDataPtrCSR(&csr_ptr, &csr_col, &csr_val, "A", nnz, nrow, nrow);
 
     // Move data to accelerator
-    A.MoveToAccelerator();
-    x.MoveToAccelerator();
-    b.MoveToAccelerator();
-    e.MoveToAccelerator();
+    if(!disable_accelerator)
+    {
+        A.MoveToAccelerator();
+        x.MoveToAccelerator();
+        b.MoveToAccelerator();
+        e.MoveToAccelerator();
+    }
 
     // Allocate x, b and e
     x.Allocate("x", A.GetN());
@@ -125,6 +130,8 @@ bool testing_bicgstabl(Arguments argus)
         p = new SGS<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "ILU")
         p = new ILU<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "ItILU0")
+        p = new ItILU0<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "ILUT")
         p = new ILUT<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "IC")
@@ -152,7 +159,7 @@ bool testing_bicgstabl(Arguments argus)
     ls.Build();
 
     // Matrix format
-    A.ConvertTo(format, format == BCSR ? 3 : 1);
+    A.ConvertTo(format, format == BCSR ? argus.blockdim : 1);
 
     ls.Solve(b, &x);
 
@@ -171,6 +178,7 @@ bool testing_bicgstabl(Arguments argus)
 
     // Stop rocALUTION platform
     stop_rocalution();
+    disable_accelerator_rocalution(false);
 
     return success;
 }
diff --git a/clients/include/testing_cg.hpp b/clients/include/testing_cg.hpp
index 7c43512..b7404fd 100644
--- a/clients/include/testing_cg.hpp
+++ b/clients/include/testing_cg.hpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2018-2022 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2018-2025 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -44,11 +44,13 @@ static bool check_residual(double res)
 template <typename T>
 bool testing_cg(Arguments argus)
 {
-    int          ndim    = argus.size;
-    std::string  precond = argus.precond;
-    unsigned int format  = argus.format;
+    int          ndim                = argus.size;
+    std::string  precond             = argus.precond;
+    unsigned int format              = argus.format;
+    bool         disable_accelerator = !argus.use_acc;
 
     // Initialize rocALUTION platform
+    disable_accelerator_rocalution(disable_accelerator);
     set_device_rocalution(device);
     init_rocalution();
 
@@ -69,10 +71,13 @@ bool testing_cg(Arguments argus)
     A.SetDataPtrCSR(&csr_ptr, &csr_col, &csr_val, "A", nnz, nrow, nrow);
 
     // Move data to accelerator
-    A.MoveToAccelerator();
-    x.MoveToAccelerator();
-    b.MoveToAccelerator();
-    e.MoveToAccelerator();
+    if(!disable_accelerator)
+    {
+        A.MoveToAccelerator();
+        x.MoveToAccelerator();
+        b.MoveToAccelerator();
+        e.MoveToAccelerator();
+    }
 
     // Allocate x, b and e
     x.Allocate("x", A.GetN());
@@ -124,6 +129,8 @@ bool testing_cg(Arguments argus)
         p = new SGS<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "ILU")
         p = new ILU<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "ItILU0")
+        p = new ItILU0<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "ILUT")
         p = new ILUT<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "IC")
@@ -150,7 +157,7 @@ bool testing_cg(Arguments argus)
     ls.Build();
 
     // Matrix format
-    A.ConvertTo(format, format == BCSR ? 3 : 1);
+    A.ConvertTo(format, format == BCSR ? argus.blockdim : 1);
 
     ls.Solve(b, &x);
 
@@ -169,6 +176,7 @@ bool testing_cg(Arguments argus)
 
     // Stop rocALUTION platform
     stop_rocalution();
+    disable_accelerator_rocalution(false);
 
     return success;
 }
diff --git a/clients/include/testing_chebyshev.hpp b/clients/include/testing_chebyshev.hpp
new file mode 100644
index 0000000..a252f24
--- /dev/null
+++ b/clients/include/testing_chebyshev.hpp
@@ -0,0 +1,218 @@
+/* ************************************************************************
+ * Copyright (C) 2025 Advanced Micro Devices, Inc. All rights Reserved.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ *
+ * ************************************************************************ */
+
+#pragma once
+
+#include "utility.hpp"
+
+#include <rocalution/rocalution.hpp>
+
+template <typename T>
+bool testing_chebyshev(Arguments argus)
+{
+    using namespace rocalution;
+
+    int          ndim                = argus.size;
+    std::string  precond             = argus.precond;
+    unsigned int format              = argus.format;
+    std::string  matrix_type         = argus.matrix_type;
+    bool         rebuildnumeric      = argus.rebuildnumeric;
+    bool         disable_accelerator = !argus.use_acc;
+
+    // Initialize rocALUTION platform
+    disable_accelerator_rocalution(disable_accelerator);
+    set_device_rocalution(device);
+    init_rocalution();
+
+    // rocALUTION structures
+    LocalMatrix<T>  A;
+    LocalVector<T>  x;
+    LocalVector<T>  b;
+    LocalVector<T>  b_old;
+    LocalVector<T>* b_k;
+    LocalVector<T>* b_k1;
+    LocalVector<T>* b_tmp;
+    LocalVector<T>  e;
+    LocalVector<T>  rhs;
+
+    // Generate A
+    int* csr_ptr = NULL;
+    int* csr_col = NULL;
+    T*   csr_val = NULL;
+
+    int nrow = 0;
+    int ncol = 0;
+    if(matrix_type == "Laplacian2D")
+    {
+        nrow = gen_2d_laplacian(ndim, &csr_ptr, &csr_col, &csr_val);
+        ncol = nrow;
+    }
+    else
+    {
+        stop_rocalution();
+        disable_accelerator_rocalution(false);
+        return true;
+    }
+    int nnz = csr_ptr[nrow];
+
+    T* csr_val2 = NULL;
+    if(rebuildnumeric)
+    {
+        csr_val2 = new T[nnz];
+        for(int i = 0; i < nnz; i++)
+        {
+            csr_val2[i] = csr_val[i];
+        }
+    }
+
+    A.SetDataPtrCSR(&csr_ptr, &csr_col, &csr_val, "A", nnz, nrow, nrow);
+
+    // Move data to accelerator
+    if(!disable_accelerator)
+    {
+        A.MoveToAccelerator();
+        x.MoveToAccelerator();
+        rhs.MoveToAccelerator();
+        e.MoveToAccelerator();
+    }
+
+    // Allocate x, b and e
+    x.Allocate("x", A.GetN());
+    rhs.Allocate("b", A.GetM());
+    e.Allocate("e", A.GetN());
+
+    T lambda_min;
+    T lambda_max;
+
+    A.Gershgorin(lambda_min, lambda_max);
+
+    // Chebyshev iteration
+    Chebyshev<LocalMatrix<T>, LocalVector<T>, T> ls;
+
+    // Initialize rhs such that A 1 = rhs
+    e.Ones();
+    A.Apply(e, &rhs);
+
+    // Initial zero guess
+    x.Zeros();
+
+    // Preconditioner
+    Preconditioner<LocalMatrix<T>, LocalVector<T>, T>* p;
+
+    if(precond == "None")
+        p = NULL;
+    else if(precond == "Chebyshev")
+    {
+        // Chebyshev preconditioner
+
+        // Determine min and max eigenvalues
+        T lambda_min;
+        T lambda_max;
+
+        A.Gershgorin(lambda_min, lambda_max);
+
+        AIChebyshev<LocalMatrix<T>, LocalVector<T>, T>* cheb
+            = new AIChebyshev<LocalMatrix<T>, LocalVector<T>, T>;
+        cheb->Set(3, lambda_max / 7.0, lambda_max);
+
+        p = cheb;
+    }
+    else if(precond == "FSAI")
+        p = new FSAI<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "SPAI")
+        p = new SPAI<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "TNS")
+        p = new TNS<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "Jacobi")
+        p = new Jacobi<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "GS")
+        p = new GS<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "SGS")
+        p = new SGS<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "ILU")
+        p = new ILU<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "ItILU0")
+        p = new ItILU0<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "ILUT")
+        p = new ILUT<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "IC")
+        p = new IC<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "MCGS")
+        p = new MultiColoredGS<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "MCSGS")
+        p = new MultiColoredSGS<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "MCILU")
+        p = new MultiColoredILU<LocalMatrix<T>, LocalVector<T>, T>;
+    else
+        return false;
+
+    // Set solver operator
+    ls.SetOperator(A);
+
+    ls.Verbose(0);
+    ls.SetOperator(A);
+
+    // Set preconditioner
+    if(p != NULL)
+    {
+        ls.SetPreconditioner(*p);
+    }
+
+    // Set eigenvalues
+    ls.Set(lambda_min, lambda_max);
+
+    // Build solver
+    ls.Build();
+
+    if(rebuildnumeric)
+    {
+        A.UpdateValuesCSR(csr_val2);
+        delete[] csr_val2;
+
+        A.Apply(e, &rhs);
+
+        ls.ReBuildNumeric();
+        ls.Set(lambda_min, lambda_max);
+    }
+
+    // Solve A x = rhs
+    ls.Solve(rhs, &x);
+
+    // Clear solver
+    ls.Clear();
+    if(p != NULL)
+    {
+        delete p;
+    }
+
+    // Compute error L2 norm
+    e.ScaleAdd(-1.0, x);
+    T error = e.Norm();
+    std::cout << "Chebyshev iteration ||e - x||_2 = " << error << std::endl;
+
+    // Stop rocALUTION platform
+    stop_rocalution();
+    disable_accelerator_rocalution(false);
+
+    return true;
+}
\ No newline at end of file
diff --git a/clients/include/testing_cr.hpp b/clients/include/testing_cr.hpp
index 41658b6..3497337 100644
--- a/clients/include/testing_cr.hpp
+++ b/clients/include/testing_cr.hpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2018-2022 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2018-2025 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -44,11 +44,13 @@ static bool check_residual(double res)
 template <typename T>
 bool testing_cr(Arguments argus)
 {
-    int          ndim    = argus.size;
-    std::string  precond = argus.precond;
-    unsigned int format  = argus.format;
+    int          ndim                = argus.size;
+    std::string  precond             = argus.precond;
+    unsigned int format              = argus.format;
+    bool         disable_accelerator = !argus.use_acc;
 
     // Initialize rocALUTION platform
+    disable_accelerator_rocalution(disable_accelerator);
     set_device_rocalution(device);
     init_rocalution();
 
@@ -69,10 +71,13 @@ bool testing_cr(Arguments argus)
     A.SetDataPtrCSR(&csr_ptr, &csr_col, &csr_val, "A", nnz, nrow, nrow);
 
     // Move data to accelerator
-    A.MoveToAccelerator();
-    x.MoveToAccelerator();
-    b.MoveToAccelerator();
-    e.MoveToAccelerator();
+    if(!disable_accelerator)
+    {
+        A.MoveToAccelerator();
+        x.MoveToAccelerator();
+        b.MoveToAccelerator();
+        e.MoveToAccelerator();
+    }
 
     // Allocate x, b and e
     x.Allocate("x", A.GetN());
@@ -124,6 +129,8 @@ bool testing_cr(Arguments argus)
         p = new SGS<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "ILU")
         p = new ILU<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "ItILU0")
+        p = new ItILU0<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "ILUT")
         p = new ILUT<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "IC")
@@ -150,7 +157,7 @@ bool testing_cr(Arguments argus)
     ls.Build();
 
     // Matrix format
-    A.ConvertTo(format, format == BCSR ? 3 : 1);
+    A.ConvertTo(format, format == BCSR ? argus.blockdim : 1);
 
     ls.Solve(b, &x);
 
@@ -169,6 +176,7 @@ bool testing_cr(Arguments argus)
 
     // Stop rocALUTION platform
     stop_rocalution();
+    disable_accelerator_rocalution(false);
 
     return success;
 }
diff --git a/clients/include/testing_fcg.hpp b/clients/include/testing_fcg.hpp
index 57e37bf..9aade14 100644
--- a/clients/include/testing_fcg.hpp
+++ b/clients/include/testing_fcg.hpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2018-2022 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2018-2025 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -44,11 +44,13 @@ static bool check_residual(double res)
 template <typename T>
 bool testing_fcg(Arguments argus)
 {
-    int          ndim    = argus.size;
-    std::string  precond = argus.precond;
-    unsigned int format  = argus.format;
+    int          ndim                = argus.size;
+    std::string  precond             = argus.precond;
+    unsigned int format              = argus.format;
+    bool         disable_accelerator = !argus.use_acc;
 
     // Initialize rocALUTION platform
+    disable_accelerator_rocalution(disable_accelerator);
     set_device_rocalution(device);
     init_rocalution();
 
@@ -69,10 +71,13 @@ bool testing_fcg(Arguments argus)
     A.SetDataPtrCSR(&csr_ptr, &csr_col, &csr_val, "A", nnz, nrow, nrow);
 
     // Move data to accelerator
-    A.MoveToAccelerator();
-    x.MoveToAccelerator();
-    b.MoveToAccelerator();
-    e.MoveToAccelerator();
+    if(!disable_accelerator)
+    {
+        A.MoveToAccelerator();
+        x.MoveToAccelerator();
+        b.MoveToAccelerator();
+        e.MoveToAccelerator();
+    }
 
     // Allocate x, b and e
     x.Allocate("x", A.GetN());
@@ -124,6 +129,8 @@ bool testing_fcg(Arguments argus)
         p = new SGS<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "ILU")
         p = new ILU<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "ItILU0")
+        p = new ItILU0<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "ILUT")
         p = new ILUT<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "IC")
@@ -150,7 +157,7 @@ bool testing_fcg(Arguments argus)
     ls.Build();
 
     // Matrix format
-    A.ConvertTo(format, format == BCSR ? 3 : 1);
+    A.ConvertTo(format, format == BCSR ? argus.blockdim : 1);
 
     ls.Solve(b, &x);
 
@@ -169,6 +176,7 @@ bool testing_fcg(Arguments argus)
 
     // Stop rocALUTION platform
     stop_rocalution();
+    disable_accelerator_rocalution(false);
 
     return success;
 }
diff --git a/clients/include/testing_fgmres.hpp b/clients/include/testing_fgmres.hpp
index c92f5af..cd39930 100644
--- a/clients/include/testing_fgmres.hpp
+++ b/clients/include/testing_fgmres.hpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2018-2022 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2018-2025 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -34,12 +34,14 @@ using namespace rocalution;
 template <typename T>
 bool testing_fgmres(Arguments argus)
 {
-    int          ndim    = argus.size;
-    int          basis   = argus.index;
-    std::string  precond = argus.precond;
-    unsigned int format  = argus.format;
+    int          ndim                = argus.size;
+    int          basis               = argus.index;
+    std::string  precond             = argus.precond;
+    unsigned int format              = argus.format;
+    bool         disable_accelerator = !argus.use_acc;
 
     // Initialize rocALUTION platform
+    disable_accelerator_rocalution(disable_accelerator);
     set_device_rocalution(device);
     init_rocalution();
 
@@ -60,10 +62,13 @@ bool testing_fgmres(Arguments argus)
     A.SetDataPtrCSR(&csr_ptr, &csr_col, &csr_val, "A", nnz, nrow, nrow);
 
     // Move data to accelerator
-    A.MoveToAccelerator();
-    x.MoveToAccelerator();
-    b.MoveToAccelerator();
-    e.MoveToAccelerator();
+    if(!disable_accelerator)
+    {
+        A.MoveToAccelerator();
+        x.MoveToAccelerator();
+        b.MoveToAccelerator();
+        e.MoveToAccelerator();
+    }
 
     // Allocate x, b and e
     x.Allocate("x", A.GetN());
@@ -115,6 +120,8 @@ bool testing_fgmres(Arguments argus)
         p = new SGS<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "ILU")
         p = new ILU<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "ItILU0")
+        p = new ItILU0<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "ILUT")
         p = new ILUT<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "IC")
@@ -143,7 +150,7 @@ bool testing_fgmres(Arguments argus)
     ls.Build();
 
     // Matrix format
-    A.ConvertTo(format, format == BCSR ? 3 : 1);
+    A.ConvertTo(format, format == BCSR ? argus.blockdim : 1);
 
     ls.Solve(b, &x);
 
@@ -162,6 +169,7 @@ bool testing_fgmres(Arguments argus)
 
     // Stop rocALUTION platform
     stop_rocalution();
+    disable_accelerator_rocalution(false);
 
     return success;
 }
diff --git a/clients/include/testing_global_matrix.hpp b/clients/include/testing_global_matrix.hpp
index de9d7d2..439d789 100644
--- a/clients/include/testing_global_matrix.hpp
+++ b/clients/include/testing_global_matrix.hpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2018-2020 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2018-2023 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -355,17 +355,13 @@ void testing_global_matrix_bad_args(void)
         ParallelManager  pm;
         LocalVector<int> lvint;
         GlobalMatrix<T>* null_mat = nullptr;
-        ParallelManager* null_pm  = nullptr;
         int*             null_int = nullptr;
-        ASSERT_DEATH(mat.CoarsenOperator(
-                         null_mat, &pm, safe_size, safe_size, lvint, safe_size, idata, safe_size),
-                     ".*Assertion.*Ac != (NULL|__null)*");
-        ASSERT_DEATH(mat.CoarsenOperator(
-                         &Ac, null_pm, safe_size, safe_size, lvint, safe_size, idata, safe_size),
-                     ".*Assertion.*pm != (NULL|__null)*");
-        ASSERT_DEATH(mat.CoarsenOperator(
-                         &Ac, &pm, safe_size, safe_size, lvint, safe_size, null_int, safe_size),
-                     ".*Assertion.*rG != (NULL|__null)*");
+        ASSERT_DEATH(
+            mat.CoarsenOperator(null_mat, safe_size, safe_size, lvint, safe_size, idata, safe_size),
+            ".*Assertion.*Ac != (NULL|__null)*");
+        ASSERT_DEATH(
+            mat.CoarsenOperator(&Ac, safe_size, safe_size, lvint, safe_size, null_int, safe_size),
+            ".*Assertion.*rG != (NULL|__null)*");
     }
 
     free_host(&idata);
diff --git a/clients/include/testing_gmres.hpp b/clients/include/testing_gmres.hpp
index 2d875b8..1790cb9 100644
--- a/clients/include/testing_gmres.hpp
+++ b/clients/include/testing_gmres.hpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2018-2022 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2018-2025 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -34,13 +34,15 @@ using namespace rocalution;
 template <typename T>
 bool testing_gmres(Arguments argus, bool expectConvergence = true)
 {
-    int          ndim    = argus.size;
-    int          basis   = argus.index;
-    std::string  matrix  = argus.matrix;
-    std::string  precond = argus.precond;
-    unsigned int format  = argus.format;
+    int          ndim                = argus.size;
+    int          basis               = argus.index;
+    std::string  matrix              = argus.matrix;
+    std::string  precond             = argus.precond;
+    unsigned int format              = argus.format;
+    bool         disable_accelerator = !argus.use_acc;
 
     // Initialize rocALUTION platform
+    disable_accelerator_rocalution(disable_accelerator);
     set_device_rocalution(device);
     init_rocalution();
 
@@ -68,10 +70,13 @@ bool testing_gmres(Arguments argus, bool expectConvergence = true)
     A.SetDataPtrCSR(&csr_ptr, &csr_col, &csr_val, "A", nnz, nrow, nrow);
 
     // Move data to accelerator
-    A.MoveToAccelerator();
-    x.MoveToAccelerator();
-    b.MoveToAccelerator();
-    e.MoveToAccelerator();
+    if(!disable_accelerator)
+    {
+        A.MoveToAccelerator();
+        x.MoveToAccelerator();
+        b.MoveToAccelerator();
+        e.MoveToAccelerator();
+    }
 
     // Allocate x, b and e
     x.Allocate("x", A.GetN());
@@ -123,6 +128,8 @@ bool testing_gmres(Arguments argus, bool expectConvergence = true)
         p = new SGS<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "ILU")
         p = new ILU<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "ItILU0")
+        p = new ItILU0<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "ILUT")
         p = new ILUT<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "IC")
@@ -151,7 +158,7 @@ bool testing_gmres(Arguments argus, bool expectConvergence = true)
     ls.Build();
 
     // Matrix format
-    A.ConvertTo(format, format == BCSR ? 3 : 1);
+    A.ConvertTo(format, format == BCSR ? argus.blockdim : 1);
 
     ls.Solve(b, &x);
 
@@ -170,6 +177,7 @@ bool testing_gmres(Arguments argus, bool expectConvergence = true)
 
     // Stop rocALUTION platform
     stop_rocalution();
+    disable_accelerator_rocalution(false);
 
     return success;
 }
diff --git a/clients/include/testing_idr.hpp b/clients/include/testing_idr.hpp
index 000b8af..35efa3a 100644
--- a/clients/include/testing_idr.hpp
+++ b/clients/include/testing_idr.hpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2018-2022 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2018-2025 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -44,12 +44,14 @@ static bool check_residual(double res)
 template <typename T>
 bool testing_idr(Arguments argus)
 {
-    int          ndim    = argus.size;
-    std::string  precond = argus.precond;
-    unsigned int format  = argus.format;
-    int          l       = argus.index;
+    int          ndim                = argus.size;
+    std::string  precond             = argus.precond;
+    unsigned int format              = argus.format;
+    int          l                   = argus.index;
+    bool         disable_accelerator = !argus.use_acc;
 
     // Initialize rocALUTION platform
+    disable_accelerator_rocalution(disable_accelerator);
     set_device_rocalution(device);
     init_rocalution();
 
@@ -70,10 +72,13 @@ bool testing_idr(Arguments argus)
     A.SetDataPtrCSR(&csr_ptr, &csr_col, &csr_val, "A", nnz, nrow, nrow);
 
     // Move data to accelerator
-    A.MoveToAccelerator();
-    x.MoveToAccelerator();
-    b.MoveToAccelerator();
-    e.MoveToAccelerator();
+    if(!disable_accelerator)
+    {
+        A.MoveToAccelerator();
+        x.MoveToAccelerator();
+        b.MoveToAccelerator();
+        e.MoveToAccelerator();
+    }
 
     // Allocate x, b and e
     x.Allocate("x", A.GetN());
@@ -85,7 +90,7 @@ bool testing_idr(Arguments argus)
     A.Apply(e, &b);
 
     // Random initial guess
-    x.SetRandomUniform(12345ULL, -4.0, 6.0);
+    x.SetRandomUniform(123456ULL, -4.0, 6.0);
 
     // Solver
     IDR<LocalMatrix<T>, LocalVector<T>, T> ls;
@@ -125,6 +130,8 @@ bool testing_idr(Arguments argus)
         p = new SGS<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "ILU")
         p = new ILU<LocalMatrix<T>, LocalVector<T>, T>;
+    else if(precond == "ItILU0")
+        p = new ItILU0<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "ILUT")
         p = new ILUT<LocalMatrix<T>, LocalVector<T>, T>;
     else if(precond == "IC")
@@ -153,7 +160,7 @@ bool testing_idr(Arguments argus)
     ls.Build();
 
     // Matrix format
-    A.ConvertTo(format, format == BCSR ? 3 : 1);
+    A.ConvertTo(format, format == BCSR ? argus.blockdim : 1);
 
     ls.Solve(b, &x);
 
@@ -172,6 +179,7 @@ bool testing_idr(Arguments argus)
 
     // Stop rocALUTION platform
     stop_rocalution();
+    disable_accelerator_rocalution(false);
 
     return success;
 }
diff --git a/clients/include/testing_inversion.hpp b/clients/include/testing_inversion.hpp
index 65cf850..a328af9 100644
--- a/clients/include/testing_inversion.hpp
+++ b/clients/include/testing_inversion.hpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2022 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2022-2025 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -44,9 +44,10 @@ static bool check_residual(double res)
 template <typename T>
 bool testing_inversion(Arguments argus)
 {
-    int          ndim        = argus.size;
-    unsigned int format      = argus.format;
-    std::string  matrix_type = argus.matrix_type;
+    int          ndim             = argus.size;
+    unsigned int format           = argus.format;
+    std::string  matrix_type      = argus.matrix_type;
+    const bool   use_host_and_acc = argus.use_acc;
 
     // Initialize rocALUTION platform
     set_device_rocalution(device);
@@ -83,12 +84,6 @@ bool testing_inversion(Arguments argus)
 
     A.SetDataPtrCSR(&csr_ptr, &csr_col, &csr_val, "A", nnz, nrow, nrow);
 
-    // Move data to accelerator
-    A.MoveToAccelerator();
-    x.MoveToAccelerator();
-    b.MoveToAccelerator();
-    e.MoveToAccelerator();
-
     // Allocate x, b and e
     x.Allocate("x", A.GetN());
     b.Allocate("b", A.GetM());
@@ -110,15 +105,39 @@ bool testing_inversion(Arguments argus)
     dls.Build();
 
     // Matrix format
-    A.ConvertTo(format, format == BCSR ? 3 : 1);
+    A.ConvertTo(format, format == BCSR ? argus.blockdim : 1);
+
+    // Move data to accelerator
+    dls.MoveToAccelerator();
+    A.MoveToAccelerator();
+    x.MoveToAccelerator();
+    b.MoveToAccelerator();
+    e.MoveToAccelerator();
 
     dls.Solve(b, &x);
 
     // Verify solution
     x.ScaleAdd(-1.0, e);
-    T nrm2 = x.Norm();
+    T nrm2_acc = x.Norm();
 
-    bool success = check_residual(nrm2);
+    bool success = check_residual(nrm2_acc);
+
+    if(use_host_and_acc)
+    {
+        dls.MoveToHost();
+        A.MoveToHost();
+        x.MoveToHost();
+        e.MoveToHost();
+        b.MoveToHost();
+
+        dls.Solve(b, &x);
+
+        // Verify solution
+        x.ScaleAdd(-1.0, e);
+        T nrm2_host = x.Norm();
+
+        success = success && check_residual(nrm2_host);
+    }
 
     // Clean up
     dls.Clear();
diff --git a/clients/include/testing_itsolver.hpp b/clients/include/testing_itsolver.hpp
new file mode 100644
index 0000000..7d76972
--- /dev/null
+++ b/clients/include/testing_itsolver.hpp
@@ -0,0 +1,165 @@
+/* ************************************************************************
+ * Copyright (C) 2025 Advanced Micro Devices, Inc. All rights Reserved.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ *
+ * ************************************************************************ */
+
+#pragma once
+
+#include "utility.hpp"
+
+#include <rocalution/rocalution.hpp>
+
+template <typename T>
+bool testing_itsolver(Arguments argus)
+{
+    using namespace rocalution;
+
+    int          ndim                = argus.size;
+    unsigned int format              = argus.format;
+    std::string  matrix_type         = argus.matrix_type;
+    bool         disable_accelerator = !argus.use_acc;
+
+    // Initialize rocALUTION platform
+    disable_accelerator_rocalution(disable_accelerator);
+    set_device_rocalution(device);
+    init_rocalution();
+
+    // rocALUTION structures
+    LocalMatrix<T> A;
+    LocalVector<T> x;
+    LocalVector<T> b;
+    LocalVector<T> e;
+
+    // Generate A
+    int* csr_ptr = NULL;
+    int* csr_col = NULL;
+    T*   csr_val = NULL;
+
+    int nrow = 0;
+    int ncol = 0;
+    if(matrix_type == "Laplacian2D")
+    {
+        nrow = gen_2d_laplacian(ndim, &csr_ptr, &csr_col, &csr_val);
+        ncol = nrow;
+    }
+    else if(matrix_type == "PermutedIdentity")
+    {
+        nrow = gen_permuted_identity(ndim, &csr_ptr, &csr_col, &csr_val);
+        ncol = nrow;
+    }
+    else
+    {
+        stop_rocalution();
+        disable_accelerator_rocalution(false);
+        return true;
+    }
+    int nnz = csr_ptr[nrow];
+
+    A.SetDataPtrCSR(&csr_ptr, &csr_col, &csr_val, "A", nnz, nrow, nrow);
+
+    // Move data to accelerator
+    if(!disable_accelerator)
+    {
+        A.MoveToAccelerator();
+        x.MoveToAccelerator();
+        b.MoveToAccelerator();
+        e.MoveToAccelerator();
+    }
+
+    // Allocate x, b and e
+    x.Allocate("x", A.GetN());
+    b.Allocate("b", A.GetM());
+    e.Allocate("e", A.GetN());
+
+    // Linear Solver
+    FixedPoint<LocalMatrix<T>, LocalVector<T>, T> fp;
+
+    // Preconditioner
+    ItILU0<LocalMatrix<T>, LocalVector<T>, T> p;
+
+    // Set iterative ILU stopping criteria
+    p.SetTolerance(1e-8);
+    p.SetMaxIter(50);
+
+    p.SetAlgorithm(ItILU0Algorithm::SyncSplit);
+
+    // Set up iterative triangular solve
+    SolverDescr descr;
+    descr.SetTriSolverAlg(TriSolverAlg_Iterative);
+    descr.SetIterativeSolverMaxIteration(30);
+    descr.SetIterativeSolverTolerance(1e-8);
+
+    descr.DisableIterativeSolverTolerance();
+    descr.EnableIterativeSolverTolerance();
+    SolverDescr descr_new(descr); // Copy the descriptor
+
+    p.SetSolverDescriptor(descr_new);
+
+    // Initialize b such that A 1 = b
+    e.Ones();
+    A.Apply(e, &b);
+
+    // Initial zero guess
+    x.Zeros();
+
+    // Set solver operator
+    fp.SetOperator(A);
+    // Set solver preconditioner
+    fp.SetPreconditioner(p);
+
+    // Build solver
+    fp.Build();
+
+    // Verbosity output
+    fp.Verbose(1);
+
+    fp.InitMinIter(1);
+    fp.InitMaxIter(1000);
+    fp.InitTol(1e-8, 1e-8, 1e-8);
+
+    // Print matrix info
+    A.Info();
+
+    // Solve A x = b
+    fp.Solve(b, &x);
+
+    int           niter_preconditioner;
+    const double* history = p.GetConvergenceHistory(&niter_preconditioner);
+
+    auto res_final = fp.GetCurrentResidual();
+    //auto res_init      = fp.GetInitialResidual();
+    //auto niter         = fp.GetNumIterations();
+    auto status_solver = fp.GetSolverStatus();
+    auto ind           = fp.GetAmaxResidualIndex();
+    // Clear solver
+    fp.Clear();
+
+    // Compute error L2 norm
+    e.ScaleAdd(-1.0, x);
+    T error = e.Norm();
+    std::cout << "||e - x||_2 = " << error << std::endl;
+
+    // Stop rocALUTION platform
+    stop_rocalution();
+    disable_accelerator_rocalution(false);
+
+    return true;
+}
\ No newline at end of file
diff --git a/clients/include/testing_local_matrix.hpp b/clients/include/testing_local_matrix.hpp
index 3145323..82e21e7 100644
--- a/clients/include/testing_local_matrix.hpp
+++ b/clients/include/testing_local_matrix.hpp
@@ -1,5 +1,5 @@
 /* ************************************************************************
- * Copyright (C) 2018-2022 Advanced Micro Devices, Inc. All rights Reserved.
+ * Copyright (C) 2018-2025 Advanced Micro Devices, Inc. All rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to deal
@@ -41,10 +41,12 @@ void testing_local_matrix_bad_args(void)
     set_device_rocalution(device);
     init_rocalution();
 
-    LocalMatrix<T>   mat1;
-    LocalMatrix<T>   mat2;
-    LocalVector<T>   vec1;
-    LocalVector<int> int1;
+    LocalMatrix<T>       mat1;
+    LocalMatrix<T>       mat2;
+    LocalVector<T>       vec1;
+    LocalVector<bool>    bool1;
+    LocalVector<int>     int1;
+    LocalVector<int64_t> int641;
 
     // null pointers
     int* null_int  = nullptr;
@@ -57,6 +59,10 @@ void testing_local_matrix_bad_args(void)
     allocate_host(safe_size, &vint);
     allocate_host(safe_size, &vdata);
 
+    // Valid matrices
+    LocalMatrix<T> mat3;
+    mat3.AllocateCSR("valid", safe_size, safe_size, safe_size);
+
     // ExtractSubMatrix, ExtractSubMatrices, Extract(Inverse)Diagonal, ExtractL/U
     {
         LocalMatrix<T>*   mat_null  = nullptr;
@@ -114,6 +120,23 @@ void testing_local_matrix_bad_args(void)
         ASSERT_DEATH(mat1.QRSolve(vec1, null_vec), ".*Assertion.*out != (NULL|__null)*");
     }
 
+    // ItLSolve, ItUSolve, ItLLSolve, ItLUSolve
+    {
+        LocalVector<T>* null_vec = nullptr;
+        int             max_iter = 1;
+        double          tol      = 0;
+        ASSERT_DEATH(mat1.ItLSolve(max_iter, tol, true, vec1, null_vec),
+                     ".*Assertion.*out != (NULL|__null)*");
+        ASSERT_DEATH(mat1.ItUSolve(max_iter, tol, true, vec1, null_vec),
+                     ".*Assertion.*out != (NULL|__null)*");
+        ASSERT_DEATH(mat1.ItLLSolve(max_iter, tol, true, vec1, null_vec),
+                     ".*Assertion.*out != (NULL|__null)*");
+        ASSERT_DEATH(mat1.ItLLSolve(max_iter, tol, true, vec1, vec1, null_vec),
+                     ".*Assertion.*out != (NULL|__null)*");
+        ASSERT_DEATH(mat1.ItLUSolve(max_iter, tol, true, vec1, null_vec),
+                     ".*Assertion.*out != (NULL|__null)*");
+    }
+
     // ICFactorize, Householder
     {
         T               val;
@@ -127,12 +150,12 @@ void testing_local_matrix_bad_args(void)
         ASSERT_DEATH(mat1.UpdateValuesCSR(null_data), ".*Assertion.*val != (NULL|__null)*");
         ASSERT_DEATH(mat1.CopyFromCSR(null_int, vint, vdata),
                      ".*Assertion.*row_offsets != (NULL|__null)*");
-        ASSERT_DEATH(mat1.CopyFromCSR(vint, null_int, vdata), ".*Assertion.*col != (NULL|__null)*");
-        ASSERT_DEATH(mat1.CopyFromCSR(vint, vint, null_data), ".*Assertion.*val != (NULL|__null)*");
+        ASSERT_DEATH(mat3.CopyFromCSR(vint, null_int, vdata), ".*Assertion.*col != (NULL|__null)*");
+        ASSERT_DEATH(mat3.CopyFromCSR(vint, vint, null_data), ".*Assertion.*val != (NULL|__null)*");
         ASSERT_DEATH(mat1.CopyToCSR(null_int, vint, vdata),
                      ".*Assertion.*row_offsets != (NULL|__null)*");
-        ASSERT_DEATH(mat1.CopyToCSR(vint, null_int, vdata), ".*Assertion.*col != (NULL|__null)*");
-        ASSERT_DEATH(mat1.CopyToCSR(vint, vint, null_data), ".*Assertion.*val != (NULL|__null)*");
+        ASSERT_DEATH(mat3.CopyToCSR(vint, null_int, vdata), ".*Assertion.*col != (NULL|__null)*");
+        ASSERT_DEATH(mat3.CopyToCSR(vint, vint, null_data), ".*Assertion.*val != (NULL|__null)*");
         ASSERT_DEATH(mat1.CopyFromCOO(null_int, vint, vdata), ".*Assertion.*row != (NULL|__null)*");
         ASSERT_DEATH(mat1.CopyFromCOO(vint, null_int, vdata), ".*Assertion.*col != (NULL|__null)*");
         ASSERT_DEATH(mat1.CopyFromCOO(vint, vint, null_data), ".*Assertion.*val != (NULL|__null)*");
@@ -173,24 +196,27 @@ void testing_local_matrix_bad_args(void)
 
     // AMG
     {
+        int                   val;
+        LocalVector<bool>*    bool_null_vec  = nullptr;
+        LocalVector<int64_t>* int64_null_vec = nullptr;
+        ASSERT_DEATH(mat1.AMGGreedyAggregate(0.1, bool_null_vec, &int641, &int641),
+                     ".*Assertion.*connections != (NULL|__null)*");
+        ASSERT_DEATH(mat1.AMGGreedyAggregate(0.1, &bool1, int64_null_vec, &int641),
+                     ".*Assertion.*aggregates != (NULL|__null)*");
+        ASSERT_DEATH(mat1.AMGGreedyAggregate(0.1, &bool1, &int641, int64_null_vec),
+                     ".*Assertion.*aggregate_root_nodes != (NULL|__null)*");
+
+        LocalMatrix<T>* null_mat = nullptr;
+        ASSERT_DEATH(mat1.AMGSmoothedAggregation(0.1, bool1, int641, int641, null_mat),
+                     ".*Assertion.*prolong != (NULL|__null)*");
+    }
+
+    {
         int               val;
         LocalVector<int>* null_vec = nullptr;
         LocalMatrix<T>*   null_mat = nullptr;
-        ASSERT_DEATH(mat1.AMGConnect(0.1, null_vec), ".*Assertion.*connections != (NULL|__null)*");
-        ASSERT_DEATH(mat1.AMGAggregate(int1, null_vec),
-                     ".*Assertion.*aggregates != (NULL|__null)*");
-        ASSERT_DEATH(mat1.AMGSmoothedAggregation(0.1, int1, int1, null_mat, &mat2),
-                     ".*Assertion.*prolong != (NULL|__null)*");
-        ASSERT_DEATH(mat1.AMGSmoothedAggregation(0.1, int1, int1, &mat2, null_mat),
-                     ".*Assertion.*restrict != (NULL|__null)*");
-        ASSERT_DEATH(mat1.AMGAggregation(int1, null_mat, &mat2),
+        ASSERT_DEATH(mat1.AMGUnsmoothedAggregation(int641, int641, null_mat),
                      ".*Assertion.*prolong != (NULL|__null)*");
-        ASSERT_DEATH(mat1.AMGAggregation(int1, &mat2, null_mat),
-                     ".*Assertion.*restrict != (NULL|__null)*");
-        ASSERT_DEATH(mat1.RugeStueben(0.1, null_mat, &mat2),
-                     ".*Assertion.*prolong != (NULL|__null)*");
-        ASSERT_DEATH(mat1.RugeStueben(0.1, &mat2, null_mat),
-                     ".*Assertion.*restrict != (NULL|__null)*");
         ASSERT_DEATH(mat1.InitialPairwiseAggregation(0.1, val, null_vec, val, &null_int, val, 0),
                      ".*Assertion.*G != (NULL|__null)*");
         ASSERT_DEATH(mat1.InitialPairwiseAggregation(0.1, val, &int1, val, &vint, val, 0),
@@ -208,12 +234,11 @@ void testing_local_matrix_bad_args(void)
                      ".*Assertion.*G != (NULL|__null)*");
         ASSERT_DEATH(mat1.FurtherPairwiseAggregation(mat2, 0.1, val, &int1, val, &null_int, val, 0),
                      ".*Assertion.*rG != (NULL|__null)*");
-        ASSERT_DEATH(mat1.CoarsenOperator(
-                         null_mat, nullptr, safe_size, safe_size, int1, safe_size, vint, safe_size),
-                     ".*Assertion.*Ac != (NULL|__null)*");
         ASSERT_DEATH(
-            mat1.CoarsenOperator(
-                &mat2, nullptr, safe_size, safe_size, int1, safe_size, null_int, safe_size),
+            mat1.CoarsenOperator(null_mat, safe_size, safe_size, int1, safe_size, vint, safe_size),
+            ".*Assertion.*Ac != (NULL|__null)*");
+        ASSERT_DEATH(
+            mat1.CoarsenOperator(&mat2, safe_size, safe_size, int1, safe_size, null_int, safe_size),
             ".*Assertion.*rG != (NULL|__null)*");
     }
 
@@ -478,16 +503,2319 @@ bool testing_local_matrix_allocations(Arguments argus)
     LocalMatrix<T> D;
     D.AllocateDIA("D", nnz, m, n, ndiag);
 
+    LocalMatrix<T> E;
+    E.AllocateMCSR("E", nnz, m, n);
+
     LocalMatrix<T> F;
     F.AllocateELL("F", ell_nnz, m, n, ell_max_row);
 
     LocalMatrix<T> G;
     G.AllocateHYB("G", ell_nnz, coo_nnz, ell_max_row, m, n);
 
+    LocalMatrix<T> H;
+    H.AllocateDENSE("H", m, n);
+
+    // Stop rocALUTION platform
+    stop_rocalution();
+
+    return true;
+}
+
+template <typename T>
+bool testing_local_matrix_zero(Arguments argus)
+{
+    int size     = argus.size;
+    int blockdim = argus.blockdim;
+
+    int m  = size;
+    int n  = size;
+    int mb = (m + blockdim - 1) / blockdim;
+    int nb = (n + blockdim - 1) / blockdim;
+
+    int nnz = 0.05 * m * n;
+    if(nnz == 0)
+    {
+        nnz = m * n;
+    }
+
+    // Initialize rocALUTION
+    set_device_rocalution(device);
+    init_rocalution();
+
+    // Testing Zeros
+    LocalMatrix<T> A;
+    A.AllocateCSR("A", nnz, m, n);
+
+    A.Zeros();
+
+    A.Info();
+
+    // Stop rocALUTION platform
+    stop_rocalution();
+
+    return true;
+}
+
+template <typename T>
+bool testing_local_matrix_set_data_ptr(Arguments argus)
+{
+    int size     = argus.size;
+    int blockdim = argus.blockdim;
+
+    int m  = size;
+    int n  = size;
+    int mb = (m + blockdim - 1) / blockdim;
+    int nb = (n + blockdim - 1) / blockdim;
+
+    int nnz = 0.05 * m * n;
+    if(nnz == 0)
+    {
+        nnz = m * n;
+    }
+
+    int nnzb = 0.01 * mb * nb;
+    if(nnzb == 0)
+    {
+        nnzb = mb * nb;
+    }
+
+    // Initialize rocALUTION
+    set_device_rocalution(device);
+    init_rocalution();
+
+    int ndiag       = 5;
+    int ell_max_row = 6;
+    int ell_nnz     = ell_max_row * m;
+    int coo_nnz     = (nnz - ell_nnz) < 0 ? 0 : nnz - ell_nnz;
+
+    // Testing allocating matrix types
+    {
+        LocalMatrix<T> A;
+        int*           row_offset = NULL;
+        int*           col        = NULL;
+        T*             val        = NULL;
+
+        allocate_host(m + 1, &row_offset);
+        allocate_host(nnz, &col);
+        allocate_host(nnz, &val);
+
+        set_to_zero_host(m + 1, row_offset);
+        set_to_zero_host(nnz, col);
+        set_to_zero_host(nnz, val);
+
+        A.SetDataPtrCSR(&row_offset, &col, &val, "A", nnz, m, n);
+        A.LeaveDataPtrCSR(&row_offset, &col, &val);
+
+        free_host(&row_offset);
+        free_host(&col);
+        free_host(&val);
+    }
+
+    {
+        LocalMatrix<T> B;
+        int*           row_offset = NULL;
+        int*           col        = NULL;
+        T*             val        = NULL;
+
+        allocate_host(mb + 1, &row_offset);
+        allocate_host(nnzb, &col);
+        allocate_host(nnzb, &val);
+
+        set_to_zero_host(mb + 1, row_offset);
+        set_to_zero_host(nnzb, col);
+        set_to_zero_host(nnzb, val);
+
+        B.SetDataPtrBCSR(&row_offset, &col, &val, "C", nnzb, mb, nb, blockdim);
+        B.LeaveDataPtrBCSR(&row_offset, &col, &val, blockdim);
+
+        free_host(&row_offset);
+        free_host(&col);
+        free_host(&val);
+    }
+
+    {
+        LocalMatrix<T> C;
+        int*           row = NULL;
+        int*           col = NULL;
+        T*             val = NULL;
+
+        allocate_host(nnz, &row);
+        allocate_host(nnz, &col);
+        allocate_host(nnz, &val);
+
+        set_to_zero_host(nnz, row);
+        set_to_zero_host(nnz, col);
+        set_to_zero_host(nnz, val);
+
+        C.SetDataPtrCOO(&row, &col, &val, "C", nnz, m, n);
+        C.LeaveDataPtrCOO(&row, &col, &val);
+
+        free_host(&row);
+        free_host(&col);
+        free_host(&val);
+    }
+
+    {
+        LocalMatrix<T> E;
+        int*           row_offset = NULL;
+        int*           col        = NULL;
+        T*             val        = NULL;
+
+        allocate_host(m + 1, &row_offset);
+        allocate_host(nnz, &col);
+        allocate_host(nnz, &val);
+
+        set_to_zero_host(m + 1, row_offset);
+        set_to_zero_host(nnz, col);
+        set_to_zero_host(nnz, val);
+
+        E.SetDataPtrMCSR(&row_offset, &col, &val, "C", nnz, m, n);
+        E.LeaveDataPtrMCSR(&row_offset, &col, &val);
+
+        free_host(&row_offset);
+        free_host(&col);
+        free_host(&val);
+    }
+
+    {
+        LocalMatrix<T> F;
+        int*           col = NULL;
+        T*             val = NULL;
+
+        allocate_host(ell_nnz, &col);
+        allocate_host(ell_nnz, &val);
+
+        set_to_zero_host(ell_nnz, col);
+        set_to_zero_host(ell_nnz, val);
+
+        F.SetDataPtrELL(&col, &val, "C", ell_nnz, m, n, ell_max_row);
+        F.LeaveDataPtrELL(&col, &val, ell_max_row);
+
+        free_host(&col);
+        free_host(&val);
+    }
+
+    {
+        LocalMatrix<T> H;
+        T*             val = NULL;
+
+        allocate_host(m * n, &val);
+
+        set_to_zero_host(m * n, val);
+
+        H.SetDataPtrDENSE(&val, "C", m, n);
+        H.LeaveDataPtrDENSE(&val);
+
+        free_host(&val);
+    }
+
     // Stop rocALUTION platform
     stop_rocalution();
 
     return true;
 }
 
+template <typename T>
+LocalMatrix<T> getTestMatrix()
+{
+    // Create a simple 2x2 CSR matrix
+    LocalMatrix<T> matrix;
+    matrix.AllocateCSR("TestMatrix", 4, 2, 2);
+
+    int row_offsets[3] = {0, 2, 4};
+    int col_indices[4] = {0, 1, 0, 1};
+    T   values[4]      = {1.0, 2.0, 3.0, 4.0};
+    matrix.CopyFromCSR(row_offsets, col_indices, values);
+
+    return matrix;
+}
+
+template <typename T>
+void getTestMatrix(Arguments argus, LocalMatrix<T>& matrix, bool& is_invertible)
+{
+    int         size        = argus.size;
+    int         blockdim    = argus.blockdim;
+    std::string matrix_type = argus.matrix_type;
+
+    // Generate A
+    int* csr_ptr = NULL;
+    int* csr_col = NULL;
+    T*   csr_val = NULL;
+
+    int nrow = 0;
+    int ncol = 0;
+    if(matrix_type == "Laplacian2D")
+    {
+        nrow = gen_2d_laplacian(size, &csr_ptr, &csr_col, &csr_val);
+        ncol = nrow;
+
+        is_invertible = true;
+    }
+    else if(matrix_type == "PermutedIdentity")
+    {
+        nrow = gen_permuted_identity(size, &csr_ptr, &csr_col, &csr_val);
+        ncol = nrow;
+
+        is_invertible = true;
+    }
+    else if(matrix_type == "Random")
+    {
+        nrow = gen_random(100 * size, 50 * size, 6, &csr_ptr, &csr_col, &csr_val);
+        ncol = 50 * size;
+
+        is_invertible = false;
+    }
+    else
+    {
+        is_invertible = true;
+
+        matrix = getTestMatrix<T>();
+        return;
+    }
+
+    int nnz = csr_ptr[nrow];
+
+    matrix.SetDataPtrCSR(&csr_ptr, &csr_col, &csr_val, "TestMatrix", nnz, nrow, ncol);
+}
+
+template <typename T>
+void getTestMatrix(Arguments argus, LocalMatrix<T>& matrix)
+{
+    bool is_invertible;
+    getTestMatrix<T>(argus, matrix, is_invertible);
+}
+
+template <typename T>
+void getMatrixVal(const LocalMatrix<T>& matrix, T* values)
+{
+    // Copy the values from the matrix to the provided array
+    int64_t m   = matrix.GetM();
+    int64_t nnz = matrix.GetNnz();
+
+    int* row_offsets   = new int[m + 1];
+    int* col_indices   = new int[nnz];
+    T*   matrix_values = new T[nnz];
+
+    matrix.CopyToCSR(row_offsets, col_indices, matrix_values);
+    for(int i = 0; i < nnz; ++i)
+    {
+        values[i] = matrix_values[i];
+    }
+
+    delete[] row_offsets;
+    delete[] col_indices;
+    delete[] matrix_values;
+}
+
+template <typename T>
+void getMatrixDiagVal(const LocalMatrix<T>& matrix, T* values)
+{
+    // Copy the values from the matrix to the provided array
+    int64_t m   = matrix.GetM();
+    int64_t nnz = matrix.GetNnz();
+
+    int* row_offsets   = new int[m + 1];
+    int* col_indices   = new int[nnz];
+    T*   matrix_values = new T[nnz];
+
+    matrix.CopyToCSR(row_offsets, col_indices, matrix_values);
+    for(int row = 0; row < m; ++row)
+    {
+        int start = row_offsets[row];
+        int end   = row_offsets[row + 1];
+        for(int i = start; i < end; ++i)
+        {
+            if(col_indices[i] == row) // Diagonal element
+            {
+                values[row] = matrix_values[i];
+                break; // Only one diagonal element per row
+            }
+        }
+    }
+
+    delete[] row_offsets;
+    delete[] col_indices;
+    delete[] matrix_values;
+}
+
+void checkPermutation(const LocalVector<int>& permutation)
+{
+    // Check that permutation is a valid permutation of 0..N-1
+    std::vector<int> seen(permutation.GetSize(), 0);
+    for(int i = 0; i < permutation.GetSize(); ++i)
+    {
+        int idx = permutation[i];
+        EXPECT_GE(idx, 0);
+        EXPECT_LT(idx, permutation.GetSize());
+        seen[idx]++;
+    }
+    for(int i = 0; i < seen.size(); ++i)
+    {
+        EXPECT_EQ(seen[i], 1); // Each index appears exactly once
+    }
+}
+
+template <typename T>
+T getTolerance()
+{
+    // Set tolerance based on the type
+    if(std::is_same<T, float>::value)
+    {
+        return 1e-5f; // Tolerance for float
+    }
+    else
+    {
+        return 1e-10; // Default tolerance for other types
+    }
+}
+
+// Helper to extract dense matrix from LocalMatrix<T>
+template <typename T>
+std::vector<std::vector<T>> extract_dense_matrix(const LocalMatrix<T>& matrix)
+{
+    int                         m   = matrix.GetM();
+    int                         n   = matrix.GetN();
+    int                         nnz = matrix.GetNnz();
+    std::vector<std::vector<T>> dense(m, std::vector<T>(n, static_cast<T>(0)));
+    std::vector<int>            row_offsets(m + 1);
+    std::vector<int>            col_indices(nnz);
+    std::vector<T>              values(nnz);
+
+    matrix.CopyToCSR(row_offsets.data(), col_indices.data(), values.data());
+    for(int row = 0; row < m; ++row)
+    {
+        for(int idx = row_offsets[row]; idx < row_offsets[row + 1]; ++idx)
+        {
+            int col         = col_indices[idx];
+            dense[row][col] = values[idx];
+        }
+    }
+    return dense;
+}
+
+template <typename T>
+void testing_local_allocate()
+{
+    // Test AllocateCSR
+    {
+        LocalMatrix<T> mat;
+        EXPECT_NO_THROW(mat.AllocateCSR("AllocatedMatrix", 4, 2, 2));
+        EXPECT_EQ(mat.GetNnz(), 4);
+        EXPECT_EQ(mat.GetM(), 2);
+        EXPECT_EQ(mat.GetN(), 2);
+    }
+
+    // Test AllocateCOO
+    {
+        LocalMatrix<T> mat;
+        mat.AllocateCOO("AllocatedMatrix", 4, 2, 2);
+        EXPECT_EQ(mat.GetNnz(), 4);
+        EXPECT_EQ(mat.GetM(), 2);
+        EXPECT_EQ(mat.GetN(), 2);
+    }
+
+    // Test AllocateBCSR
+    {
+        LocalMatrix<T> mat;
+        int            nnzb = 4, mb = 2, nb = 2, blockdim = 2;
+        EXPECT_NO_THROW(mat.AllocateBCSR("BCSR", nnzb, mb, nb, blockdim));
+        EXPECT_EQ(mat.GetNnz(), nnzb * blockdim * blockdim);
+        EXPECT_EQ(mat.GetM(), mb * blockdim);
+        EXPECT_EQ(mat.GetN(), nb * blockdim);
+    }
+
+    // Test AllocateDIA
+    {
+        LocalMatrix<T> mat;
+        int            nnz = 6, m = 3, n = 3, ndiag = 2;
+        EXPECT_NO_THROW(mat.AllocateDIA("DIA", nnz, m, n, ndiag));
+        EXPECT_EQ(mat.GetNnz(), nnz);
+        EXPECT_EQ(mat.GetM(), m);
+        EXPECT_EQ(mat.GetN(), n);
+    }
+
+    // Test AllocateMCSR
+    {
+        LocalMatrix<T> mat;
+        int            nnz = 5, m = 3, n = 3;
+        EXPECT_NO_THROW(mat.AllocateMCSR("MCSR", nnz, m, n));
+        EXPECT_EQ(mat.GetNnz(), nnz);
+        EXPECT_EQ(mat.GetM(), m);
+        EXPECT_EQ(mat.GetN(), n);
+    }
+
+    // Test AllocateELL
+    {
+        LocalMatrix<T> mat;
+        int            ell_nnz = 6, m = 3, n = 3, ell_max_row = 2;
+        EXPECT_NO_THROW(mat.AllocateELL("ELL", ell_nnz, m, n, ell_max_row));
+        EXPECT_EQ(mat.GetNnz(), ell_nnz);
+        EXPECT_EQ(mat.GetM(), m);
+        EXPECT_EQ(mat.GetN(), n);
+    }
+
+    // Test AllocateHYB
+    {
+        LocalMatrix<T> mat;
+        int            m = 3, n = 3;
+        int            ell_max_row = 2;
+        int            ell_nnz     = ell_max_row * m; // 2 * 3 = 6
+        int            coo_nnz     = 2;
+
+        EXPECT_NO_THROW(mat.AllocateHYB("HYB", ell_nnz, coo_nnz, ell_max_row, m, n));
+        EXPECT_EQ(mat.GetNnz(), ell_nnz + coo_nnz); // 8
+        EXPECT_EQ(mat.GetM(), m);
+        EXPECT_EQ(mat.GetN(), n);
+    }
+
+    // Test AllocateDENSE
+    {
+        LocalMatrix<T> mat;
+        int            m = 3, n = 3;
+        EXPECT_NO_THROW(mat.AllocateDENSE("DENSE", m, n));
+        EXPECT_EQ(mat.GetNnz(), m * n);
+        EXPECT_EQ(mat.GetM(), m);
+        EXPECT_EQ(mat.GetN(), n);
+    }
+}
+
+template <typename T>
+void testing_check_with_empty_matrix()
+{
+    LocalMatrix<T> empty_matrix;
+    // Check should pass without any issues
+    EXPECT_NO_THROW(empty_matrix.Check());
+    // Info should not throw an error
+    EXPECT_NO_THROW(empty_matrix.Info());
+}
+
+template <typename T>
+void testing_local_copy_from_async()
+{
+    auto           matrix = getTestMatrix<T>();
+    LocalMatrix<T> copy_matrix;
+    // CopyFromAsync should copy the matrix asynchronously (if supported)
+    EXPECT_NO_THROW(copy_matrix.CopyFromAsync(matrix));
+    EXPECT_NO_THROW(matrix.Sync());
+    EXPECT_EQ(copy_matrix.GetM(), matrix.GetM());
+    EXPECT_EQ(copy_matrix.GetN(), matrix.GetN());
+    EXPECT_EQ(copy_matrix.GetNnz(), matrix.GetNnz());
+
+    // Compare dense representations
+    auto dense_orig = extract_dense_matrix(matrix);
+    auto dense_copy = extract_dense_matrix(copy_matrix);
+
+    EXPECT_EQ(dense_orig.size(), dense_copy.size());
+    for(size_t i = 0; i < dense_orig.size(); ++i)
+    {
+        EXPECT_EQ(dense_orig[i].size(), dense_copy[i].size());
+        for(size_t j = 0; j < dense_orig[i].size(); ++j)
+        {
+            EXPECT_EQ(dense_orig[i][j], dense_copy[i][j]);
+        }
+    }
+}
+
+template <typename T>
+void testing_local_update_values_csr()
+{
+    auto matrix = getTestMatrix<T>();
+
+    int64_t nnz = matrix.GetNnz();
+
+    // Use std::vector instead of raw arrays
+    std::vector<T> new_values(nnz);
+    for(int64_t i = 0; i < nnz; ++i)
+    {
+        new_values[i] = static_cast<T>(i + 10); // Fill with some values
+    }
+
+    // UpdateValuesCSR should update the values in the matrix
+    EXPECT_NO_THROW(matrix.UpdateValuesCSR(new_values.data()));
+
+    std::vector<T> check_values(nnz);
+    getMatrixVal(matrix, check_values.data());
+    for(int64_t i = 0; i < nnz; ++i)
+    {
+        EXPECT_EQ(check_values[i], new_values[i]);
+    }
+}
+
+template <typename T>
+void testing_local_move_to_accelerator()
+{
+    auto matrix = getTestMatrix<T>();
+
+    for(int i = 0; i < 10; i++)
+    {
+        EXPECT_NO_THROW(matrix.MoveToAccelerator());
+    }
+    for(int i = 0; i < 10; i++)
+    {
+        EXPECT_NO_THROW(matrix.MoveToHost());
+    }
+    for(int i = 0; i < 10; i++)
+    {
+        EXPECT_NO_THROW(matrix.MoveToHost());
+        EXPECT_NO_THROW(matrix.MoveToAccelerator());
+    }
+    EXPECT_EQ(matrix.Check(), true);
+}
+
+template <typename T>
+void testing_local_move_to_accelerator_async()
+{
+    auto matrix = getTestMatrix<T>();
+    // MoveToAcceleratorAsync should move the matrix asynchronously
+    for(int i = 0; i < 10; i++)
+    {
+        EXPECT_NO_THROW(matrix.MoveToAcceleratorAsync());
+    }
+    for(int i = 0; i < 10; i++)
+    {
+        EXPECT_NO_THROW(matrix.MoveToHostAsync());
+    }
+    for(int i = 0; i < 10; i++)
+    {
+        EXPECT_NO_THROW(matrix.MoveToHostAsync());
+        EXPECT_NO_THROW(matrix.MoveToAcceleratorAsync());
+    }
+    EXPECT_NO_THROW(matrix.Sync());
+    EXPECT_EQ(matrix.Check(), true);
+}
+
+template <typename T>
+void testing_local_move_to_host_async()
+{
+    auto matrix = getTestMatrix<T>();
+    // MoveToHostAsync should move the matrix asynchronously to host
+    EXPECT_NO_THROW(matrix.MoveToHostAsync());
+    EXPECT_NO_THROW(matrix.Sync());
+    EXPECT_EQ(matrix.GetM(), 2);
+    EXPECT_EQ(matrix.GetN(), 2);
+}
+
+template <typename T>
+void testing_local_clear(Arguments argus)
+{
+    LocalMatrix<T> matrix;
+    getTestMatrix<T>(argus, matrix);
+
+    // Clear should remove all data from the matrix
+    // This test checks if the Clear operation is valid
+    // by checking if the number of non-zero entries (nnz) is zero
+    // and the dimensions (m, n) are also zero.
+
+    matrix.Clear();
+    EXPECT_EQ(matrix.GetNnz(), 0);
+    EXPECT_EQ(matrix.GetM(), 0);
+    EXPECT_EQ(matrix.GetN(), 0);
+}
+
+template <typename T>
+void testing_local_zeros(Arguments argus)
+{
+    LocalMatrix<T> matrix;
+    getTestMatrix<T>(argus, matrix);
+
+    // Zeros should set all values in the matrix to zero
+    // This test checks if the Zeros operation is valid
+    // by checking if all values in the matrix are zero after the operation.
+
+    matrix.Zeros();
+    matrix.Check();
+    matrix.Info();
+
+    int64_t nnz          = matrix.GetNnz();
+    T*      check_values = new T[nnz];
+    getMatrixVal(matrix, check_values);
+    for(int64_t i = 0; i < nnz; ++i)
+    {
+        EXPECT_EQ(check_values[i], static_cast<T>(0));
+    }
+    delete[] check_values;
+}
+
+template <typename T>
+void testing_local_copy(Arguments argus)
+{
+    LocalMatrix<T> matrix;
+    getTestMatrix<T>(argus, matrix);
+
+    // CopyFrom should create a copy of the matrix
+    // This test checks if the copy operation is valid
+    // by comparing the original matrix with the copied matrix.
+
+    LocalMatrix<T> copy_matrix;
+    copy_matrix.CopyFrom(matrix);
+    EXPECT_EQ(copy_matrix.GetNnz(), matrix.GetNnz());
+    EXPECT_EQ(copy_matrix.GetM(), matrix.GetM());
+    EXPECT_EQ(copy_matrix.GetN(), matrix.GetN());
+
+    // Compare dense representations
+    auto dense_orig = extract_dense_matrix(matrix);
+    auto dense_copy = extract_dense_matrix(copy_matrix);
+
+    EXPECT_EQ(dense_orig.size(), dense_copy.size());
+    for(size_t i = 0; i < dense_orig.size(); ++i)
+    {
+        EXPECT_EQ(dense_orig[i].size(), dense_copy[i].size());
+        for(size_t j = 0; j < dense_orig[i].size(); ++j)
+        {
+            EXPECT_EQ(dense_orig[i][j], dense_copy[i][j]);
+        }
+    }
+}
+
+template <typename T>
+void testing_local_scale(Arguments argus)
+{
+    LocalMatrix<T> matrix;
+    getTestMatrix<T>(argus, matrix);
+
+    // Scale should multiply all values in the matrix by a scalar
+    // This test checks if the scaling operation is valid
+    // by comparing the scaled values with the expected values.
+    // The expected values are obtained by multiplying the original values
+    // by the scaling factor.
+
+    // Save original dense matrix
+    auto orig_dense = extract_dense_matrix(matrix);
+
+    // Scale the matrix by 2.0
+    matrix.Scale(2.0);
+
+    // Extract new dense matrix
+    auto new_dense = extract_dense_matrix(matrix);
+
+    // Compare each value
+    int m = matrix.GetM();
+    int n = matrix.GetN();
+    for(int i = 0; i < m; ++i)
+        for(int j = 0; j < n; ++j)
+            EXPECT_EQ(new_dense[i][j], orig_dense[i][j] * 2.0);
+}
+
+template <typename T>
+void testing_local_extract_diagonal(Arguments argus)
+{
+    if(argus.matrix_type != "Laplacian2D")
+    {
+        return;
+    }
+
+    LocalMatrix<T> matrix;
+    getTestMatrix<T>(argus, matrix);
+
+    // ExtractDiagonal should extract the diagonal of the matrix
+    // and store it in a LocalVector
+    // This test checks if the diagonal extraction is correct
+    // by comparing the extracted diagonal with the expected values.
+    // The expected values are obtained by iterating through the matrix
+    // and checking the diagonal elements.
+
+    LocalVector<T> diag;
+    matrix.ExtractDiagonal(&diag);
+
+    int64_t m            = matrix.GetM();
+    T*      check_values = new T[m];
+    getMatrixDiagVal(matrix, check_values);
+
+    EXPECT_EQ(diag.GetSize(), m);
+    for(int i = 0; i < m; ++i)
+    {
+        EXPECT_EQ(diag[i], check_values[i]);
+    }
+    delete[] check_values;
+}
+
+template <typename T>
+void testing_local_extract_inverse_diagonal(Arguments argus)
+{
+    if(argus.matrix_type != "Laplacian2D")
+    {
+        return;
+    }
+
+    LocalMatrix<T> matrix;
+    getTestMatrix<T>(argus, matrix);
+
+    // ExtractInverseDiagonal should extract the inverse diagonal of the matrix
+    // and store it in a LocalVector
+    // This test checks if the inverse diagonal extraction is correct
+    // by comparing the extracted inverse diagonal with the expected values.
+    // The expected values are obtained by taking the reciprocal of the diagonal elements.
+
+    LocalVector<T> diag;
+    matrix.ExtractDiagonal(&diag);
+    LocalVector<T> inv_diag;
+    matrix.ExtractInverseDiagonal(&inv_diag);
+    EXPECT_EQ(inv_diag.GetSize(), diag.GetSize());