Repository: aatmdelissen/pyMOTO (language: python)

Issue:
Convolution filter
I've added a convolution filter that is much faster than the regular matrix-based density filter.
Reasons to add:
- Setup cost is next to none (for H-matrix it can go up to seconds-minutes with a large filter radius and large domain). For a 3D domain of 60^3 with filter radius 5 elements, setup of H-matrix is 12 seconds, while for convolution it is 0.0059.
- Execution cost is low, since convolution makes use of efficient FFT algorithms. For the 3D domain aforementioned, running convolution costs only 0.023s, while H-matrix is almost 10x slower at 0.16s.
- Flexible boundary conditions: padding is required to extend the domain for the convolution to work. This allows different regions to have different behavior (e.g. reflection, constant value, wrap around, ...).


Discussion / hints:

Fix patch:
diff --git a/examples/ex_simple_module.py b/examples/ex_simple_module.py
index fb91db0..65fe8b1 100644
--- a/examples/ex_simple_module.py
+++ b/examples/ex_simple_module.py
@@ -1,65 +1,89 @@
-""" Example: Module template for a generic module
+""" Example: Module templates for generic modules
 It shows possibilities for a generic Module, and also the importance of the sensitivity reset.
 """
 from pymoto import Module, Signal, finite_difference
 
 
-class MyNewModule(Module):
-    """ This is my module that does awesome stuff.
-
-    The (Module) is required for your module to be a module, with correct behaviour.
-    This example is a module with two inputs and two outputs.
-    """
-    def _prepare(self, my_argument, optional_arg='foo'):
-        """
-        This function is called by the initialization of the superclass. It can use the input parameters entered into
-        the Module.create() method.
-        """
-        print('>> Prepare my module')
-        print('\tmy_argument = {}'.format(my_argument))
-        print('\toptional_argument = {}'.format(optional_arg))
-
+class MostSimple(Module):
+    """ Example of a minimal module with two inputs and one output """
     def _response(self, x1, x2):
+        """ Forward-path calculation is put here """
+        print(f'[{type(self).__name__}] Do my response calculation')
+        print(f'\tInputs are {self.sig_in[0].tag} and {self.sig_in[1].tag}')
+        return x1 * x2
+
+    def _sensitivity(self, df_dy):
+        """ Backward-path sensitivity calculation here.
+        In case df_dy is None, the function will automatically not be called
         """
-        This function calculates a response based on a multiple input values, here for example 2. Multiple outputs can
-        easily be added. Also, different response behaviours can be implemented, based on the number of inputs (function
-        overloading). The 'self' object can be used to save state variables.
-
-        :param x1: First input variable
-        :param x2: Second input variable
-        :return: The results of the calculation
-        """
-        print('>> Do my response calculation')
-
-        # Incorrect data
-        if x1 is None or x2 is None:
-            raise RuntimeError("You forgot to set {} and {}".format(self.sig_in[0].tag, self.sig_in[1].tag))
-
+        print(f'[{type(self).__name__}] Do my sensitivity calculation')
+        x1, x2 = [s.state for s in self.sig_in]
+        df_dx1 = df_dy * x2  # Apply chain rule df/dx1 = df/dy * dy/dx1
+        df_dx2 = df_dy * x1
+        return df_dx1, df_dx2  # Return the sensitivities with respect to input signals
+
+
+class WithPrepare(Module):
+    """ This module takes a parameter during initialization
+    Example usage:
+    >>> import pymoto as pym
+    >>> x = pym.Signal('x', state=3.5)
+    >>> y = pym.Signal('y')
+    >>> m = WithPrepare(x, y, 1.2, optional_value='bar')
+    [WithPrepare] Prepare my module
+    value = 1.2
+    optional_value = bar
+    >>> m.response() # doctest: +ELLIPSIS
+    [WithPrepare] Do my response calculation
+    Message = bar, x = 3.5, y=4.2
+    ...
+    >>> y.state
+    4.2
+    """
+    def _prepare(self, value, optional_value='foo'):
+        """ This prepare is called during initialization of the module, and can be used for set-up """
+        print(f'[{type(self).__name__}] Prepare my module')
+        print('value = {}'.format(value))
+        print('optional_value = {}'.format(optional_value))
+        self.value = value
+        self.optional_value = optional_value
+
+    def _response(self, x):
+        print(f'[{type(self).__name__}] Do my response calculation')
+        y = x * self.value
+        print(f"Message = {self.optional_value}, x = {x}, y={y}")
+        return y
+
+    def _sensitivity(self, df_dy):
+        print(f'[{type(self).__name__}] Do my sensitivity calculation')
+        return df_dy * self.value
+
+
+class TwoOutputs(Module):
+    """ This module has two inputs and two outputs """
+    def _response(self, x1, x2):
+        print(f'[{type(self).__name__}]Do my response calculation')
         # Store data, which might be needed for the sensitivity calculation
         self.x1 = x1
         self.x2 = x2
 
         # Calculate two response values
-        v1 = x1 * x2
-        v2 = x1 + x2
+        y1 = x1 * x2
+        y2 = x1 + x2
 
         # Return the results
-        return v1, v2
-
-    def _sensitivity(self, df_dv1, df_dv2):
-        """
-        This function calculate the (backward) sensitivity. It should handle None (zero sensitivity) as incoming adjoint
-        variable. If both are None, the sensitivity will not be called.
+        return y1, y2
 
-        :param df_dv1: The adjoint variable of the first output
-        :param df_dv2: The adjoint variable of the second output
-        :return: The adjoint variables (sensitivities) of the inputs
+    def _sensitivity(self, df_dy1, df_dy2):
+        """ This function calculate the (backward) sensitivity.
+        It should handle None (zero sensitivity) as incoming adjoint variable. If both are None, the sensitivity
+        will not be called.
         """
-        print('>> Do my sensitivity calculation')
+        print(f'[{type(self).__name__}]Do my sensitivity calculation')
 
         # Calculate the gradients with chain-rule
         # First initialize sensitivities with the correct size containing all zeros
-        df_dx1 = self.x1 * 0
+        df_dx1 = self.x1 * 0  # The sensitivity df/dx1 is the same size as x1 (in case of a vector/matrix)
         df_dx2 = self.x2 * 0
 
         # In case the data of x1 and x2 were not stored, it could still be obtained here by directly accessing the state
@@ -71,108 +95,90 @@ def _sensitivity(self, df_dv1, df_dv2):
 
         # If the sensitivity of the output signal is empty, it is None. So we only need to do calculations whenever it
         # is not None. In case both sensitivities of the output signals are None, this function won't be called.
-        if df_dv1 is not None:
-            df_dx1 += df_dv1*self.x2
-            df_dx2 += df_dv1*self.x1
+        if df_dy1 is not None:
+            df_dx1 += df_dy1*self.x2
+            df_dx2 += df_dy1*self.x1
 
-        if df_dv2 is not None:
-            df_dx1 += df_dv2
-            df_dx2 += df_dv2
+        if df_dy2 is not None:
+            df_dx1 += df_dy2
+            df_dx2 += df_dy2
 
         # Return the results
         return df_dx1, df_dx2
 
-    def _reset(self):
-        """
-        This function is used for reset, called whenever a new sensitivity is to be calculated, triggered by calling
-        reset(). This can be used to reset any internal storage required for sensitivity calculation
-        """
-        print('>> Reset my module')
-
 
 if __name__ == "__main__":
     print(__doc__)
     print("_" * 80)
-    print("PART 1: Setup")
-
-    # The print function lists all possible Module types
-    print("Show all possible modules defined, it also lists the new locally defined module")
-    Module.print_children()
+    print("-- Module setup")
 
     # Create signals for the inputs. The argument is the 'tag' of the signal, which is optional.
-    # The signal's 'tag' can be seen as its name, which can be useful for printing or when
+    # The tag of the signal can be seen as its name, which can be useful for printing and debugging
     x1 = Signal("x1")
 
     # Also create a second input signal (as our module has two inputs)
     x2 = Signal("x2")
 
-    # And create two signals for the outputs
-    y1 = Signal("y1")
-    y2 = Signal("y2")
-
-    # The module can be instantiated using the constructor. The first argument is a list of input signals, and the
-    # second argument is a list of output arguments. The other (keyword) arguments are directly passed to the _prepare()
-    # function of the module.
-    print("\nInstantiate directly:")
-    the_mod = MyNewModule([x1, x2], [y1, y2], 'my_arg_data', optional_arg=3.14)
-
-    # It can also be instantiated using Module.create, by referencing its name (case insensitive) as first argument.
-    # All the other arguments remain identical to before.
-    print("\nInstantiate by Module.create:")
-    also_the_mod = Module.create("MyNewModule", [x1, x2], [y1, y2], 'my_arg_data1', optional_arg='bar')
+    # And create a signal for the output
+    y = Signal("y")
 
-    print("_"*80)
-    print("PART 2: Forward analysis")
+    # The module is instantiated using the constructor. The first argument is the input signal (or a list of input
+    # signals), and the second argument is the output signal (or again a list). Different module types are prepared in
+    # this example. Uncomment lines below to see what happens.
+    print("Instantiate Module:")
+    simple_module = MostSimple([x1, x2], y)
+    # simple_module = WithPrepare(x1, y, 3.14, optional_value='bar')  # Module with extra (constant) parameters
+    # simple_module = TwoOutputs([x1, x2], [y, Signal('y2')])  # Module with two outputs
 
+    print("\n-- Forward analysis")
     try:
-        print("\nTrying to call the response without setting initial values results in an error")
-        the_mod.response()
-    except RuntimeError as e:
-        print("ERROR OBTAINED: ", e.__str__())
+        print("Trying to call the response without setting initial values results in an error")
+        simple_module.response()
+    except TypeError as e:
+        print("ERROR OBTAINED: \n TypeError: ", e.__str__())
 
-    # Set the initial values to the state of the input signals using <Signal>.state
-    x1.state = 2.0
+    # To correctly do the forward calculation, a corect state need to be set in x1 and x2
+    x1.state = 2.0  # Values for the for the forward calculation
     x2.state = 3.0
+    print(f"\nState initialized to {x1.tag} = {x1.state}, {x2.tag} = {x2.state}")
 
-    print("\nState initialized to {0} = {1}, {2} = {3}".format(x1.tag, x1.state, x2.tag, x2.state))
+    # Now execute response
     print("Call response")
-    the_mod.response()
-    # The states of the output signals can be accessed using <Signal>.state
-    print("The results: {0} = {1}, {2} = {3}".format(y1.tag, y1.state, y2.tag, y2.state))
+    simple_module.response()
 
-    # Sensitivity analysis
-    print("_"*80)
-    print("PART 3: Backpropagation (sensitivity analysis)")
+    # The state of the output signal can be accessed using <Signal>.state again
+    print(f"The result: {y.tag} = {y.state}")
 
+    print("\n-- Sensitivity analysis by back-propagation")
     # Calculate sensitivities
     print("\nIf no seed is given, no sensitivities will be calculated")
-    the_mod.sensitivity()
-    print("dg/d{} = {}".format(x1.tag, x1.sensitivity))
-    print("dg/d{} = {}".format(x2.tag, x2.sensitivity))
+    simple_module.sensitivity()
+    print(f"dy/d{x1.tag} = {x1.sensitivity}")
+    print(f"dy/d{x2.tag} = {x2.sensitivity}")
 
-    print("\nSeed dg/dy1 = 1.0, so we can calculate dy1/dx1 and dy1/dx2")
-    # An initial 'seed' sensitivity of the last response you're interested in needs to be set. We can do this by setting
+    print("\nSeed dy/dy = 1.0, so we can calculate dy/dx1 and dy/dx2")
+    # An initial 'seed' sensitivity of the response you're interested in needs to be set. We can do this by setting
     # the <Signal>.sensitivity property
-    y1.sensitivity = 1.0
-    the_mod.sensitivity()
-    # The sensitivities of the initial signals can also be accessed by <Signal>.sensitivity
-    print("dy1/d{} = {}".format(x1.tag, x1.sensitivity))
-    print("dy1/d{} = {}".format(x2.tag, x2.sensitivity))
+    y.sensitivity = 1.0
+    simple_module.sensitivity()
+    # The sensitivities of the input signals can now be accessed by <Signal>.sensitivity
+    print(f"dy/d{x1.tag} = {x1.sensitivity}")
+    print(f"dy/d{x2.tag} = {x2.sensitivity}")
 
     print("\nWhen reset is not called after the sensitivity calculation, the results will not be correct.")
-    print("Seed dg/dy1 = 1.0 again (not strictly necessary, since the value already was seeded)")
-    y1.sensitivity = 1.0
-    the_mod.sensitivity()
-    print("Incorrect sensitivity dy1/d{} = {}".format(x1.tag, x1.sensitivity))
-    print("Incorrect sensitivity dy1/d{} = {}".format(x2.tag, x2.sensitivity))
+    print("Seed dy/dy = 1.0 again (not strictly necessary, since the value already was seeded without being cleared)")
+    y.sensitivity = 1.0
+    simple_module.sensitivity()
+    print(f"Incorrect sensitivity dy/d{x1.tag} = {x1.sensitivity}")
+    print(f"Incorrect sensitivity dy/d{x2.tag} = {x2.sensitivity}")
     print("The values are now double of what they're supposed to be, because they're added to what we already had.")
 
-    print("\nRESET! And seed dg/dy2 = 1.0 to calculate the other sensitivities")
-    the_mod.reset()  # !! DON'T FORGET TO RESET, ELSE SENSITIVITIES FROM PREVIOUS RUNS WILL CONTAMINATE YOUR RESULT !!
-    y2.sensitivity = 1.0
-    the_mod.sensitivity()
-    print("dy2/d{} = {}".format(x1.tag, x1.sensitivity))
-    print("dy2/d{} = {}".format(x2.tag, x2.sensitivity))
+    print("\nRESET! And seed dy/dy = 1.0 to calculate the other sensitivities")
+    simple_module.reset()  # !! DON'T FORGET TO RESET, ELSE SENSITIVITIES FROM PREVIOUS RUNS WILL CONTAMINATE YOUR RESULT !!
+    y.sensitivity = 1.0
+    simple_module.sensitivity()
+    print(f"dy/d{x1.tag} = {x1.sensitivity}")
+    print(f"dy/d{x2.tag} = {x2.sensitivity}")
 
     # You can always check your module with finite differencing
-    finite_difference(the_mod)
+    finite_difference(simple_module)
diff --git a/examples/example1.py b/examples/example1.py
index 25499ef..452b17c 100644
--- a/examples/example1.py
+++ b/examples/example1.py
@@ -74,6 +74,7 @@ def _sensitivity(self, df_dy):
     z.state = 3.4
 
     # Create the modules
+    # Depending on how the input and output signals are routed between the modules, different behavior can be created
     ordering = 0
     if ordering == 0:
         # A __
@@ -98,14 +99,14 @@ def _sensitivity(self, df_dy):
     # Create a network of modules
     func = Network(m1, m2, m3)
 
-    print("\nCurrent interconnection:")
+    print("\nCurrent network:")
     print(" -> ".join([type(m).__name__ for m in func.mods]))
 
     # --- FORWARD ANALYSIS ---
     # Perform forward analysis
     func.response()
 
-    print("\nThe response is  g(x={0}, y={1}, z={2}) = {3}".format(x.state, y.state, z.state, g.state))
+    print(f"\nThe response is  g(x={x.state}, y={y.state}, z={z.state}) = {g.state}")
 
     # --- BACKPROPAGATION ---
     # Clear previous sensitivities
@@ -118,9 +119,9 @@ def _sensitivity(self, df_dy):
     func.sensitivity()
 
     print("\nThe sensitivities are:")
-    print("d{0}/d{1} = {2}".format(g.tag, x.tag, x.sensitivity))
-    print("d{0}/d{1} = {2}".format(g.tag, y.tag, y.sensitivity))
-    print("d{0}/d{1} = {2}".format(g.tag, z.tag, z.sensitivity))
+    print(f"d{g.tag}/d{x.tag} = {x.sensitivity}")
+    print(f"d{g.tag}/d{y.tag} = {y.sensitivity}")
+    print(f"d{g.tag}/d{z.tag} = {z.sensitivity}")
 
     # --- Finite difference checks ---
     # On the individual modules
diff --git a/examples/topology_optimization/ex_compliance.py b/examples/topology_optimization/ex_compliance.py
index 9d7a3a7..016de13 100644
--- a/examples/topology_optimization/ex_compliance.py
+++ b/examples/topology_optimization/ex_compliance.py
@@ -3,11 +3,11 @@
 
 import pymoto as pym
 
-nx, ny, nz = 30, 30, 0  # Set nz to zero for the 2D problem
+nx, ny, nz = 120, 40, 0  # Set nz to zero for the 2D problem, nz > 0 runs a 3D problem
 xmin = 1e-9
 filter_radius = 2.0
 volfrac = 0.5
-thermal = False  # True = Static thermal analysis; False = Static mechanical analysis will be done
+thermal = True  # True = Static thermal analysis; False = Static mechanical analysis will be done
 
 if __name__ == "__main__":
     print(__doc__)
@@ -17,12 +17,12 @@
         domain = pym.DomainDefinition(nx, ny)
 
         if thermal:
+            ndof = 1  # Number of dofs per node
             # Get dof numbers at the boundary
             boundary_dofs = domain.get_nodenumber(0, np.arange(ny // 4, (ny + 1) - ny // 4))
 
             # Make a force vector
             force_dofs = domain.get_nodenumber(*np.meshgrid(np.arange(1, nx + 1), np.arange(ny + 1)))
-            ndof = 1  # Number of dofs per node
 
         else:  # Mechanical
             ndof = 2
@@ -42,7 +42,7 @@
             force_dofs = domain.get_nodenumber(*np.meshgrid(np.arange(1, nx+1), np.arange(ny + 1), np.arange(nz + 1))).flatten()
             ndof = 1
 
-        else:
+        else:  # Mechanical
             ndof = 3
             boundary_dofs = np.repeat(boundary_nodes * ndof, ndof, axis=-1) + np.tile(np.arange(ndof), len(boundary_nodes))
             force_dofs = ndof * domain.get_nodenumber(nx, ny // 2, nz // 2) + 2  # Z-direction
@@ -80,9 +80,9 @@
 
     # Linear system solver. The linear solver can be chosen by uncommenting any of the following lines.
     solver = None  # Default (solver is automatically chosen based on matrix properties)
-    # solver = pym.SolverSparsePardiso()  # Requires Intel MKL installed
-    # solver = pym.SolverSparseCholeskyCVXOPT()  # Requires cvxopt installed
-    # solver = pym.SolverSparseCholeskyScikit()  # Requires scikit installed
+    # solver = pym.solvers.SolverSparsePardiso()  # Requires Intel MKL installed
+    # solver = pym.solvers.SolverSparseCholeskyCVXOPT()  # Requires cvxopt installed
+    # solver = pym.solvers.SolverSparseCholeskyScikit()  # Requires scikit installed
     su = func.append(pym.LinSolve([sK, sf], hermitian=True, solver=solver))
 
     # Output the design, deformation, and force field to a Paraview file
@@ -101,7 +101,7 @@
     svol.tag = 'volume'
 
     # Volume constraint
-    sg1 = func.append(pym.MathGeneral(svol, expression='10*(inp0/{} - {})'.format(domain.nel, volfrac)))
+    sg1 = func.append(pym.MathGeneral(svol, expression=f'10*(inp0/{domain.nel} - {volfrac})'))
     sg1.tag = "volume constraint"
 
     # Maybe you want to check the design-sensitivities?
@@ -112,8 +112,8 @@
 
     func.append(pym.PlotIter([sg0, sg1]))  # Plot iteration history
 
-    # Do the optimization with MMA (requires nlopt package)
-    # pym.minimize_mma(func, [sx], [sg0, sg1])  # TODO does not work correctly for the thermal case
+    # Do the optimization with MMA
+    # pym.minimize_mma(func, [sx], [sg0, sg1])
 
     # Do the optimization with OC
     pym.minimize_oc(func, sx, sg0)
diff --git a/examples/topology_optimization/ex_compliance_multigrid.py b/examples/topology_optimization/ex_compliance_multigrid.py
new file mode 100644
index 0000000..d493638
--- /dev/null
+++ b/examples/topology_optimization/ex_compliance_multigrid.py
@@ -0,0 +1,138 @@
+""" Minimal example for a compliance topology optimization with multigrid preconditioned CG as solver """
+import numpy as np
+
+import pymoto as pym
+
+nx, ny, nz = 128, 64, 64  # Set nz to zero for the 2D problem, nz > 0 runs a 3D problem
+xmin = 1e-9
+filter_radius = 2.0
+volfrac = 0.5
+thermal = False  # True = Static thermal analysis; False = Static mechanical analysis will be done
+
+if __name__ == "__main__":
+    print(__doc__)
+
+    if nz == 0:  # 2D analysis
+        # Generate a grid
+        domain = pym.DomainDefinition(nx, ny)
+
+        if thermal:
+            ndof = 1  # Number of dofs per node
+            # Get dof numbers at the boundary
+            boundary_dofs = domain.get_nodenumber(0, np.arange(ny // 4, (ny + 1) - ny // 4))
+
+            # Make a force vector
+            force_dofs = domain.get_nodenumber(*np.meshgrid(np.arange(1, nx + 1), np.arange(ny + 1)))
+
+        else:  # Mechanical
+            ndof = 2
+            # Calculate boundary dof indices
+            boundary_nodes = domain.get_nodenumber(0, np.arange(ny + 1))
+            boundary_dofs = np.repeat(boundary_nodes * ndof, ndof, axis=-1) + np.tile(np.arange(ndof), len(boundary_nodes))
+
+            # Which dofs to put a force on? The 1 is added for a force in y-direction (x-direction would be zero)
+            force_dofs = ndof * domain.get_nodenumber(nx, ny // 2) + 1
+
+    else:
+        domain = pym.DomainDefinition(nx, ny, nz)
+        boundary_nodes = domain.get_nodenumber(*np.meshgrid(0, range(ny + 1), range(nz + 1))).flatten()
+
+        if thermal:
+            boundary_dofs = boundary_nodes
+            force_dofs = domain.get_nodenumber(*np.meshgrid(np.arange(1, nx+1), np.arange(ny + 1), np.arange(nz + 1))).flatten()
+            ndof = 1
+
+        else:  # Mechanical
+            ndof = 3
+            boundary_dofs = np.repeat(boundary_nodes * ndof, ndof, axis=-1) + np.tile(np.arange(ndof), len(boundary_nodes))
+            force_dofs = ndof * domain.get_nodenumber(nx, ny // 2, nz // 2) + 2  # Z-direction
+
+    if domain.nnodes > 1e+6:
+        print("Too many nodes :(")  # Safety, to prevent overloading the memory in your machine
+        exit()
+
+    # Generate a force vector
+    f = np.zeros(domain.nnodes * ndof)
+    f[force_dofs] = 1.0  # Uniform force of 1.0 at all selected dofs
+
+    # Make force and design vector, and fill with initial values
+    sf = pym.Signal('f', state=f)
+    sx = pym.Signal('x', state=np.ones(domain.nel) * volfrac)
+
+    # Start building the modular network
+    func = pym.Network(print_timing=True)
+
+    # Filter
+    sxfilt = func.append(pym.FilterConv(sx, domain=domain, radius=filter_radius))
+    sx_analysis = sxfilt
+
+    # Show the design on the screen as it optimizes (plotting domain in 3D may take significant time)
+    # func.append(pym.PlotDomain(sx_analysis, domain=domain, saveto="out/design", clim=[0, 1]))
+
+    # SIMP material interpolation
+    sSIMP = func.append(pym.MathGeneral(sx_analysis, expression=f"{xmin} + {1.0 - xmin}*inp0^3"))
+
+    # System matrix assembly module
+    if thermal:
+        sK = func.append(pym.AssemblePoisson(sSIMP, domain=domain, bc=boundary_dofs))
+    else:
+        sK = func.append(pym.AssembleStiffness(sSIMP, domain=domain, bc=boundary_dofs))
+
+    # Linear system solver. The linear solver can be chosen by uncommenting any of the following lines.
+    solver = None  # Default (solver is automatically chosen based on matrix properties)
+    # solver = pym.solvers.SolverSparsePardiso()  # Requires Intel MKL installed
+    # solver = pym.solvers.SolverSparseCholeskyCVXOPT()  # Requires cvxopt installed
+    # solver = pym.solvers.SolverSparseCholeskyScikit()  # Requires scikit installed
+
+    ''' Iterative solver: CG with Geometric Multi-grid preconditioning '''
+    # Set up fine level
+    mg1 = pym.solvers.GeometricMultigrid(domain)
+    mgs = [mg1]
+    n_levels = 5  # Number of levels (including fine grid)
+    for i in range(n_levels - 1):  # Setup coarse levels
+        mgs.append(pym.solvers.GeometricMultigrid(mgs[i].sub_domain))
+        mgs[i].inner_level = mgs[i+1]
+
+    # Set up the solver (comment out to use the default factorization, try this to see the difference in time)
+    solver = pym.solvers.CG(preconditioner=mg1, verbosity=1, tol=1e-5)
+
+    ''' From here the rest of the example is identical to ex_compliance.py '''
+    su = func.append(pym.LinSolve([sK, sf], hermitian=True, solver=solver))
+
+    # Output the design, deformation, and force field to a Paraview file
+    func.append(pym.WriteToVTI([sx_analysis, su, sf], domain=domain, saveto='out/dat.vti'))
+
+    # Compliance calculation c = f^T u
+    scompl = func.append(pym.EinSum([su, sf], expression='i,i->'))
+    scompl.tag = 'compliance'
+
+    # MMA needs correct scaling of the objective
+    sg0 = func.append(pym.Scaling(scompl, scaling=100.0))
+    sg0.tag = "objective"
+
+    # Calculate the volume of the domain by adding all design densities together
+    svol = func.append(pym.EinSum(sx_analysis, expression='i->'))
+    svol.tag = 'volume'
+
+    # Volume constraint
+    sg1 = func.append(pym.MathGeneral(svol, expression=f'10*(inp0/{domain.nel} - {volfrac})'))
+    sg1.tag = "volume constraint"
+
+    # Maybe you want to check the design-sensitivities?
+    do_finite_difference = False
+    if do_finite_difference:
+        pym.finite_difference(func, sx, [sg0, sg1], dx=1e-4)
+        exit()
+
+    func.append(pym.PlotIter([sg0, sg1]))  # Plot iteration history
+
+    # Do the optimization with MMA
+    # pym.minimize_mma(func, [sx], [sg0, sg1])
+
+    # Do the optimization with OC
+    pym.minimize_oc(func, sx, sg0)
+
+    # Here you can do some post processing
+    print("The optimization has finished!")
+    print(f"The final compliance value obtained is {scompl.state}")
+    print(f"The maximum {'temperature' if thermal else 'displacement'} is {max(np.absolute(su.state))}")
diff --git a/examples/topology_optimization/ex_compliance_padding_filter.py b/examples/topology_optimization/ex_compliance_padding_filter.py
new file mode 100644
index 0000000..8370912
--- /dev/null
+++ b/examples/topology_optimization/ex_compliance_padding_filter.py
@@ -0,0 +1,163 @@
+""" Compliance topology optimization with a padded-domain filter, which avoids the design 'sticking' to the boundary """
+import numpy as np
+
+import pymoto as pym
+
+nx, ny, nz = 120, 40, 0  # Set nz to zero for the 2D problem
+xmin = 1e-9
+filter_radius = 5.0
+volfrac = 0.5
+
+
+class Continuation(pym.Module):
+    """ Module that generates a continuated value """
+    def _prepare(self, start=0.0, stop=1.0, nsteps=80, stepstart=10):
+        self.startval = start
+        self.endval = stop
+        self.dval = (stop - start) / nsteps
+        self.nstart = stepstart
+        self.iter = -1
+        self.val = self.startval
+
+    def _response(self):
+        if (self.val < self.endval) and (self.iter > self.nstart):
+            self.val += self.dval
+
+        self.val = np.clip(self.val, min(self.startval, self.endval), max(self.startval, self.endval))
+        print(self.sig_out[0].tag, ' = ', self.val)
+        self.iter += 1
+        return self.val
+
+    def _sensitivity(self, *args):
+        pass
+
+
+if __name__ == "__main__":
+    print(__doc__)
+
+    if nz == 0:  # 2D analysis
+        ndof = 2
+
+        # Generate a grid
+        domain = pym.DomainDefinition(nx, ny)
+
+        # Get node numbers for the boundary condition
+        boundary_nodes = domain.get_nodenumber(0, np.arange(ny // 3, (ny+1) - ny//3))
+
+        # Calculate boundary dof indices
+        boundary_dofs = np.repeat(boundary_nodes * ndof, ndof, axis=-1) + np.tile(np.arange(ndof), len(boundary_nodes))
+
+        # Which dofs to put a force on? The 1 is added for a force in y-direction (x-direction would be zero)
+        force_dofs = ndof * domain.get_nodenumber(nx, ny // 2) + 1
+
+    else:
+        ndof = 3
+        domain = pym.DomainDefinition(nx, ny, nz)
+        boundary_nodes = domain.get_nodenumber(*np.meshgrid(0, np.arange(ny // 3, (ny+1) - ny//3), np.arange(nz // 3, (nz+1) - nz//3))).flatten()
+        boundary_dofs = np.repeat(boundary_nodes * ndof, ndof, axis=-1) + np.tile(np.arange(ndof), len(boundary_nodes))
+        force_dofs = ndof * domain.get_nodenumber(nx, ny // 2, nz // 2) + 2  # Z-direction
+
+    if domain.nnodes > 1e+6:
+        print("Too many nodes :(")  # Safety, to prevent overloading the memory in your machine
+        exit()
+
+    # Generate a force vector
+    f = np.zeros(domain.nnodes * ndof)
+    f[force_dofs] = 1.0  # Uniform force of 1.0 at all selected dofs
+
+    # Make force and design vector, and fill with initial values
+    sf = pym.Signal('f', state=f)
+    sx = pym.Signal('x', state=np.ones(domain.nel) * volfrac)
+
+    # Start building the modular network
+    func = pym.Network(print_timing=False)
+
+    do_padding = True  # Switch padding on or off
+    if do_padding:
+        ''' Filter with padding to avoid the design 'sticking' to the boundaries
+        Set densities just outside the boundaries to 0, except for the xmax. The xmax boundary is kept symmetric, which 
+        is the standard behavior of a density filter. This is done to allow full solid to be formed here in case the 
+        mechanical load is at the boundary.
+        '''
+        sxfilt = func.append(pym.FilterConv(sx, domain=domain, radius=filter_radius,
+                                            xmin_bc=0, xmax_bc='symmetric',
+                                            ymin_bc=0, ymax_bc=0,
+                                            zmin_bc=0, zmax_bc=0))
+        m_filt = func[-1]
+        ''' Furthermore, a density of 1 is set just outside the boundary condition at xmin for correct filtering '''
+        xrange = np.arange(m_filt.pad_sizes[0])
+        yrange = m_filt.pad_sizes[1] + np.arange(domain.nely//3, domain.nely - domain.nely//3)
+        zrange = np.array([0]) if domain.dim < 3 else m_filt.pad_sizes[2] + np.arange(domain.nelz//3, domain.nelz - domain.nelz//3)
+        ex, ey, ez = np.meshgrid(xrange, yrange, zrange)
+        m_filt.override_padded_values((ex, ey, ez), 1.0)
+        ''' Everything from here is the same as in ex_compliance.py '''
+    else:
+        # Use regular density filter (without padding) to see the difference
+        sxfilt = func.append(pym.DensityFilter(sx, domain=domain, radius=filter_radius))
+
+    # Heaviside projections
+    etaDi, etaNo, etaEr = 0.3, 0.5, 0.7
+    sBeta = pym.Signal("beta")
+    func.append(Continuation([], sBeta, start=1.0, stop=20.0, stepstart=10))
+
+    sxNom = pym.Signal("xnominal")
+    sxEr = pym.Signal("xeroded")
+    sxDi = pym.Signal("xdilated")
+
+    heaviside = "(tanh(inp1 * {0}) + tanh(inp1 * (inp0 - {0}))) / (tanh(inp1 * {0}) + tanh(inp1 * (1 - {0})))"
+    func.append(pym.MathGeneral([sxfilt, sBeta], sxNom, expression=heaviside.format(etaNo)))
+    func.append(pym.MathGeneral([sxfilt, sBeta], sxEr, expression=heaviside.format(etaEr)))
+    func.append(pym.MathGeneral([sxfilt, sBeta], sxDi, expression=heaviside.format(etaDi)))
+
+    sx_analysis = sxEr
+
+    # Show the design on the screen as it optimizes
+    func.append(pym.PlotDomain(sx_analysis, domain=domain, saveto="out/design", clim=[0, 1]))
+
+    # SIMP material interpolation
+    sSIMP = func.append(pym.MathGeneral(sx_analysis, expression=f"{xmin} + {1.0 - xmin}*inp0^3"))
+
+    # System matrix assembly module
+    sK = func.append(pym.AssembleStiffness(sSIMP, domain=domain, bc=boundary_dofs))
+
+    # Linear system solver. The linear solver can be chosen by uncommenting any of the following lines.
+    solver = None  # Default (solver is automatically chosen based on matrix properties)
+    # solver = pym.solvers.SolverSparsePardiso()  # Requires Intel MKL installed
+    # solver = pym.solvers.SolverSparseCholeskyCVXOPT()  # Requires cvxopt installed
+    # solver = pym.solvers.SolverSparseCholeskyScikit()  # Requires scikit installed
+    su = func.append(pym.LinSolve([sK, sf], hermitian=True, solver=solver))
+
+    # Compliance calculation c = f^T u
+    scompl = func.append(pym.EinSum([su, sf], expression='i,i->'))
+    scompl.tag = 'compliance'
+
+    # MMA needs correct scaling of the objective
+    sg0 = func.append(pym.Scaling(scompl, scaling=100.0))
+    sg0.tag = "objective"
+
+    # Calculate the volume of the domain by adding all design densities together
+    svol = func.append(pym.EinSum(sxNom, expression='i->'))
+    svol.tag = 'volume'
+
+    # Volume constraint
+    sg1 = func.append(pym.MathGeneral(svol, expression=f'10*(inp0/{domain.nel} - {volfrac})'))
+    sg1.tag = "volume constraint"
+
+    # Maybe you want to check the design-sensitivities?
+    do_finite_difference = False
+    if do_finite_difference:
+        pym.finite_difference(func, sx, [sg0, sg1], dx=1e-4)
+        exit()
+
+    func.append(pym.PlotIter([sg0, sg1]))  # Plot iteration history
+
+    # Do the optimization with MMA
+    pym.minimize_mma(func, [sx], [sg0, sg1])
+
+    # Do the optimization with OC
+    # pym.minimize_oc(func, sx, sg0)
+
+    # Here you can do some post processing
+    print("The optimization has finished!")
+    print(f"The final compliance value obtained is {scompl.state}")
+    print(f"The maximum displacement is {max(np.absolute(su.state))}")
diff --git a/examples/topology_optimization/ex_volume_stress.py b/examples/topology_optimization/ex_volume_stress.py
index af1a604..ca60eec 100644
--- a/examples/topology_optimization/ex_volume_stress.py
+++ b/examples/topology_optimization/ex_volume_stress.py
@@ -28,54 +28,6 @@
 max_displacement = 20.0
 
 
-class Stress(pym.Module):
-    def _prepare(self, E=1, nu=0.3, plane='strain', domain=pym.DomainDefinition, *args, **kwargs):
-        siz = domain.element_size
-        self.domain = domain
-
-        # Constitutive model
-        self.D = siz[2] * get_D(E, nu, plane.lower())
-
-        # Numerical integration
-        self.B = np.zeros((3, 8), dtype=float)
-        w = np.prod(siz[:domain.dim] / 2)
-        for n in domain.node_numbering:
-            pos = n * (siz / 2) / np.sqrt(3)  # Sampling point
-            dN_dx = domain.eval_shape_fun_der(pos)
-            self.B += w * get_B(dN_dx)
-
-        self.dofconn = domain.get_dofconnectivity(2)
-
-    def _response(self, u):
-        self.elemental_strain = self.B.dot(u[self.dofconn].transpose())
-        self.elemental_strain[2, :] *= 2  # voigt notation
-        return self.D.dot(self.elemental_strain).transpose()
-
-    def _sensitivity(self, dfdv):
-        dgdsstrainmat = np.einsum('jk,kl->jl', dfdv, self.D)
-        dgdsstrainmat[:, 2] *= 2
-        dgdue = np.einsum('ij,jl->il', dgdsstrainmat, self.B)
-
-        y = np.zeros(self.domain.nnodes * 2)
-        for i in range(0, self.domain.nel):
-            y[self.dofconn[i, :]] += dgdue[i, :]
-        return y
-
-
-class VonMises(pym.Module):
-    def _prepare(self, *args, **kwargs):
-        # Vandermonde matrix
-        self.V = np.array([[1, -0.5, 0], [-0.5, 1, 0], [0, 0, 3]])
-
-    def _response(self, x):
-        self.x = x
-        self.y = (x.dot(self.V) * x).sum(1)
-        return np.sqrt(self.y)
-
-    def _sensitivity(self, dfdv):
-        return dfdv[:, np.newaxis] * (self.y ** (-0.5))[:, np.newaxis] * self.x.dot(self.V)
-
-
 class ConstraintAggregation(pym.Module):
     """
     Unified aggregation and relaxation.
@@ -101,7 +53,7 @@ def _response(self, x):
         self.x = x
         self.y = self.x + 1
         self.z = self.y ** self.P
-        z = ((1 / len(self.x)) * np.sum(self.z)) ** (1 / self.P)  # P-mean aggregation function
+        z = (np.sum(self.z) / self.n) ** (1 / self.P)  # P-mean aggregation function
         return z - 1
 
     def _sensitivity(self, dfdc):
@@ -139,14 +91,18 @@ def _sensitivity(self, dfdc):
     s_force = pym.Signal('f', state=f)
     s_displacement = fn.append(pym.LinSolve([s_K, s_force]))
 
-    # Calculate stress
-    s_stress = fn.append(Stress([s_displacement], domain=domain))
-    s_stress_vm = fn.append(VonMises([s_stress]))
-    s_stress_constraints = fn.append(pym.Scaling([s_stress_vm], maxval=maximum_vm_stress, scaling=1.0))
+    # Calculate stress components
+    s_stress = fn.append(pym.Stress(s_displacement, domain=domain))
 
-    s_stress_constraints_scaled = fn.append(
-        pym.EinSum([s_filtered_variables, s_stress_constraints], expression='i,i->i'))
+    # Calculate Von-Mises stress
+    V = np.array([[1, -0.5, 0], [-0.5, 1, 0], [0, 0, 3]])  # Vandermonde matrix
+    s_stress_vm2 = fn.append(pym.EinSum([s_stress, pym.Signal(state=V), s_stress], expression='ij,ik,kj->j'))
+    s_stress_vm = fn.append(pym.MathGeneral(s_stress_vm2, expression='sqrt(inp0)'))
 
+    # Stress constraint
+    s_stress_constraints = fn.append(pym.Scaling(s_stress_vm, maxval=maximum_vm_stress, scaling=1.0))
+    s_stress_constraints_scaled = fn.append(pym.EinSum([s_filtered_variables, s_stress_constraints],
+                                                       expression='i,i->i'))
     s_stress_constraint = fn.append(ConstraintAggregation([s_stress_constraints_scaled], P=10))
     s_stress_constraint.tag = "Stress constraint"
 
diff --git a/pymoto/__init__.py b/pymoto/__init__.py
index 823cded..c0af7d0 100644
--- a/pymoto/__init__.py
+++ b/pymoto/__init__.py
@@ -5,21 +5,23 @@
 # Imports from common
 from .common.dyadcarrier import DyadCarrier
 from .common.mma import MMA
-from .common.solvers import matrix_is_complex, matrix_is_diagonal, matrix_is_symmetric, matrix_is_hermitian, LinearSolver, LDAWrapper
-from .common.solvers_dense import SolverDiagonal, SolverDenseQR, SolverDenseLU, SolverDenseCholesky, SolverDenseLDL
-from .common.solvers_sparse import SolverSparsePardiso, SolverSparseLU, SolverSparseCholeskyScikit, SolverSparseCholeskyCVXOPT
+
+# Import solvers
+from . import solvers
 
 # Modular inports
 from .core_objects import Signal, Module, Network, make_signals
 
 # Import modules
 from .modules.assembly import AssembleGeneral, AssembleStiffness, AssembleMass, AssemblePoisson
+from .modules.assembly import ElementOperation, Strain, Stress
 from .modules.autodiff import AutoMod
 from .modules.complex import MakeComplex, RealPart, ImagPart, ComplexNorm
 from .modules.filter import FilterConv, Filter, DensityFilter, OverhangFilter
 from .modules.generic import MathGeneral, EinSum, ConcatSignal
-from .modules.io import PlotDomain, PlotGraph, PlotIter, WriteToVTI
+from .modules.io import FigModule, PlotDomain, PlotGraph, PlotIter, WriteToVTI
 from .modules.linalg import Inverse, LinSolve, EigenSolve, SystemOfEquations, StaticCondensation
+from .modules.aggregation import AggScaling, AggActiveSet, Aggregation, PNorm, SoftMinMax, KSFunction
 from .modules.scaling import Scaling
 
 # Further helper routines
@@ -28,21 +30,25 @@
 __all__ = [
     'Signal', 'Module', 'Network', 'make_signals',
     'finite_difference', 'minimize_oc', 'minimize_mma',
+
     # Common
     'MMA',
     'DyadCarrier',
     'DomainDefinition',
-    'matrix_is_complex', 'matrix_is_diagonal', 'matrix_is_symmetric', 'matrix_is_hermitian',
-    'LinearSolver', 'LDAWrapper',
-    'SolverDiagonal', 'SolverDenseQR', 'SolverDenseLU', 'SolverDenseCholesky', 'SolverDenseLDL',
-    'SolverSparsePardiso', 'SolverSparseLU', 'SolverSparseCholeskyScikit', 'SolverSparseCholeskyCVXOPT',
+    'solvers',
+
+    # Helpers
+    "AggScaling", "AggActiveSet",
+
     # Modules
     "MathGeneral", "EinSum", "ConcatSignal",
     "Inverse", "LinSolve", "EigenSolve", "SystemOfEquations", "StaticCondensation",
     "AssembleGeneral", "AssembleStiffness", "AssembleMass", "AssemblePoisson",
+    "ElementOperation", "Strain", "Stress",
     "FilterConv", "Filter", "DensityFilter", "OverhangFilter",
-    "PlotDomain", "PlotGraph", "PlotIter", "WriteToVTI",
+    "FigModule", "PlotDomain", "PlotGraph", "PlotIter", "WriteToVTI",
     "MakeComplex", "RealPart", "ImagPart", "ComplexNorm",
     "AutoMod",
-    "Scaling"
+    "Aggregation", "PNorm", "SoftMinMax", "KSFunction",
+    "Scaling",
 ]
diff --git a/pymoto/common/domain.py b/pymoto/common/domain.py
index facd32b..95e2883 100644
--- a/pymoto/common/domain.py
+++ b/pymoto/common/domain.py
@@ -5,6 +5,31 @@
 import warnings
 from typing import Union
 import numpy as np
+from matplotlib.patches import PathPatch
+from matplotlib.path import Path
+
+
+def plot_deformed_element(ax, x, y, **kwargs):
+    codes, verts = zip(*[
+        (Path.MOVETO, [x[0], y[0]]),
+        (Path.LINETO, [x[1], y[1]]),
+        (Path.LINETO, [x[3], y[3]]),
+        (Path.LINETO, [x[2], y[2]]),
+        (Path.CLOSEPOLY, [x[0], y[0]])])
+    path = Path(verts, codes)
+    patch = PathPatch(path, **kwargs)
+    ax.add_artist(patch)
+    return patch
+
+
+def get_path(x, y):
+    codes, verts = zip(*[
+        (Path.MOVETO, [x[0], y[0]]),
+        (Path.LINETO, [x[1], y[1]]),
+        (Path.LINETO, [x[3], y[3]]),
+        (Path.LINETO, [x[2], y[2]]),
+        (Path.CLOSEPOLY, [x[0], y[0]])])
+    return Path(verts, codes)
 
 
 class DomainDefinition:
@@ -100,6 +125,13 @@ def __init__(self, nelx: int, nely: int, nelz: int = 0, unitx: float = 1.0, unit
         self.conn = np.zeros((self.nel, self.elemnodes), dtype=int)
         self.conn[el, :] = self.get_elemconnectivity(elx, ely, elz)
 
+        # Helper for element slicing
+        eli, elj, elk = np.meshgrid(np.arange(self.nelx), np.arange(self.nely), np.arange(self.nelz), indexing='ij')
+        self.elements = self.get_elemnumber(eli, elj, elk)
+
+        ndi, ndj, ndk = np.meshgrid(np.arange(self.nelx+1), np.arange(self.nely+1), np.arange(self.nelz+1), indexing='ij')
+        self.nodes = self.get_nodenumber(ndi, ndj, ndk)
+
     def get_elemnumber(self, eli: Union[int, np.ndarray], elj: Union[int, np.ndarray], elk: Union[int, np.ndarray] = 0):
         """ Gets the element number(s) for element(s) with given Cartesian indices (i, j, k)
 
@@ -126,7 +158,7 @@ def get_nodenumber(self, nodi: Union[int, np.ndarray], nodj: Union[int, np.ndarr
         """
         return (nodk * (self.nely + 1) + nodj) * (self.nelx + 1) + nodi
 
-    def get_node_indices(self, nod_idx: Union[int, np.ndarray]):
+    def get_node_indices(self, nod_idx: Union[int, np.ndarray] = None):
         """ Gets the Cartesian index (i, j, k) for given node number(s)
 
         Args:
@@ -135,16 +167,18 @@ def get_node_indices(self, nod_idx: Union[int, np.ndarray]):
         Returns:
             i, j, k for requested node(s); k is only returned in 3D
         """
+        if nod_idx is None:
+            nod_idx = np.arange(self.nnodes)
         nodi = nod_idx % (self.nelx + 1)
         nodj = (nod_idx // (self.nelx + 1)) % (self.nely + 1)
         if self.dim == 2:
-            return nodi, nodj
+            return np.stack([nodi, nodj], axis=0)
         nodk = nod_idx // ((self.nelx + 1)*(self.nely + 1))
-        return nodi, nodj, nodk
+        return np.stack([nodi, nodj, nodk], axis=0)
 
-    def get_node_position(self, nod_idx: Union[int, np.ndarray]):
+    def get_node_position(self, nod_idx: Union[int, np.ndarray] = None):
         ijk = self.get_node_indices(nod_idx)
-        return [idx * self.element_size[ii] for ii, idx in enumerate(ijk)]
+        return (self.element_size[:self.dim] * ijk.T).T
 
     def get_elemconnectivity(self, i: Union[int, np.ndarray], j: Union[int, np.ndarray], k: Union[int, np.ndarray] = 0):
         """ Get the connectivity for element identified with Cartesian indices (i, j, k)
@@ -230,6 +264,27 @@ def eval_shape_fun_der(self, pos: np.ndarray):
             dN_dx[i, :] *= np.array([n[i] for n in self.node_numbering])  # Flip +/- signs according to node position
         return dN_dx
 
+    def plot(self, ax, deformation=None, scaling=None):
+        patches = []
+        for e in range(self.nel):
+            n = self.conn[e]
+            x, y = self.get_node_position(n)
+            u, v = deformation[n * 2], deformation[n * 2 + 1]
+            color = (1 - scaling[e], 1 - scaling[e], 1 - scaling[e]) if scaling is not None else 'grey'
+            patch = plot_deformed_element(ax, x + u, v + y, linewidth=0.1, color=color)
+            patches.append(patch)
+        return patches
+
+    def update_plot(self, patches, deformation=None, scaling=None):
+        for e in range(self.nel):
+            patch = patches[e]
+            n = self.conn[e]
+            x, y = self.get_node_position(n)
+            u, v = deformation[n * 2], deformation[n * 2 + 1]
+            color = (1 - scaling[e], 1 - scaling[e], 1 - scaling[e]) if scaling is not None else 'grey'
+            patch.set_color(color)
+            patch.set_path(self.get_path(x + u, y + v))
+
     # flake8: noqa: C901
     def write_to_vti(self, vectors: dict, filename="out.vti", scale=1.0, origin=(0.0, 0.0, 0.0)):
         """ Write all given vectors to a Paraview (VTI) file
diff --git a/pymoto/common/dyadcarrier.py b/pymoto/common/dyadcarrier.py
index 776d1ae..53ac7c1 100644
--- a/pymoto/common/dyadcarrier.py
+++ b/pymoto/common/dyadcarrier.py
@@ -1,8 +1,8 @@
-from typing import Union, Iterable
+from typing import Union, Iterable, List
 import warnings
 import numpy as np
 from numpy.typing import NDArray
-from scipy.sparse import spmatrix
+from scipy.sparse import spmatrix, coo_matrix
 from ..utils import _parse_to_list
 try:  # Import fast optimized einsum
     from opt_einsum import contract as einsum
@@ -31,7 +31,7 @@ def isnullslice(x):
 
 
 class DyadCarrier(object):
-    """ Efficient storage for dyadic or rank-N matrix
+    r""" Efficient storage for dyadic or rank-N matrix
 
     Stores only the vectors instead of creating a full rank-N matrix
     :math:`\mathbf{A} = \sum_k^N \mathbf{u}_k\otimes\mathbf{v}_k`
@@ -67,7 +67,7 @@ def size(self):
             return self.ulen * self.vlen
 
     def add_dyad(self, u: Iterable, v: Iterable = None, fac: float = None):
-        """ Adds a list of vectors to the dyad carrier
+        r""" Adds a list of vectors to the dyad carrier
 
         Checks for conforming sizes of `u` and `v`. The data inside the vectors are copied.
 
@@ -373,6 +373,31 @@ def contract(self, mat: Union[NDArray, spmatrix] = None, rows: NDArray[int] = No
 
         return val
 
+    def contract_multi(self, mats: List[spmatrix], dtype=None):
+        """ Faster version of contraction for a list of sparse matrices """
+        if dtype is None:
+            dtype = np.result_type(self.dtype, mats[0].dtype)
+        val = np.zeros(len(mats), dtype=dtype)
+
+        if len(self.u) == 0 or len(self.v) == 0:
+            return val
+        U = np.array(self.u).T
+        V = np.array(self.v).T
+
+        for i, m in enumerate(mats):
+            if m is None:
+                vali = 0.0
+            else:
+                try:
+                    if not isinstance(m, coo_matrix):
+                        warnings.warn("Inefficiency: Matrix must be converted to coo_matrix for contraction")
+                    mat_coo = m.tocoo()
+                    vali = np.einsum('ij,i,ij->', U[mat_coo.row, :], mat_coo.data, V[mat_coo.col, :])
+                except AttributeError:
+                    vali = self.contract(m)
+            val[i] = vali
+        return val
+
     def todense(self):
         """ Returns a full (dense) matrix from the DyadCarrier matrix """
         warning_size = 100e+6  # Bytes
@@ -387,6 +412,10 @@ def todense(self):
 
         return val
 
+    def toarray(self):
+        """ Convert to array, same as todense(). To be consistent with scipy.sparse """
+        return self.todense()
+
     def iscomplex(self):
         """ Check if the DyadCarrier is of complex type """
         return np.iscomplexobj(np.array([], dtype=self.dtype))
diff --git a/pymoto/common/mma.py b/pymoto/common/mma.py
index aa98687..38b4c6c 100644
--- a/pymoto/common/mma.py
+++ b/pymoto/common/mma.py
@@ -219,11 +219,16 @@ class MMA:
         xmin: Minimum design variable (can be a vector)
         xmax: Maximum design variable (can be a vector)
         fn_callback: A function that is called just before calling the response() in each iteration
-        verbosity: 0 - No prints, 1 - Only convergence message, 2 - Convergence and iteration info, 3 - Extended info
+        verbosity: Level of information to print
+          0 - No prints
+          1 - Only convergence message
+          2 - Convergence and iteration info (default)
+          3 - Additional info on variables
+          4 - Additional info on sensitivity information
 
     """
 
-    def __init__(self, function, variables, responses, tolx=1e-4, tolf=0.0, move=0.1, maxit=100, xmin=0.0, xmax=1.0, fn_callback=None, verbosity=0, **kwargs):
+    def __init__(self, function, variables, responses, tolx=1e-4, tolf=0.0, move=0.1, maxit=100, xmin=0.0, xmax=1.0, fn_callback=None, verbosity=2, **kwargs):
         self.funbl = function
         self.verbosity = verbosity
 
@@ -301,8 +306,7 @@ def response(self):
                 self.xmin[self.cumlens[i]:self.cumlens[i+1]] = xminvals[i]
 
         if len(self.xmin) != self.n:
-            raise RuntimeError(
-                "Length of the xmin vector not correct ({} != {})".format(len(self.xmin), self.n))
+            raise RuntimeError(f"Length of the xmin vector ({len(self.xmin)}) should be equal to # design variables ({self.n})")
 
         if not hasattr(self.xmax, '__len__'):
             self.xmax = self.xmax * np.ones_like(xval)
@@ -313,17 +317,19 @@ def response(self):
                 self.xmax[self.cumlens[i]:self.cumlens[i + 1]] = xmaxvals[i]
 
         if len(self.xmax) != self.n:
-            raise RuntimeError("Length of the xmax vector not correct ({} != {})".format(len(self.xmax), self.n))
+            raise RuntimeError(f"Length of the xmax vector ({len(self.xmax)}) should be equal to # design variables ({self.n})")
 
-        # Set movelimit in case of multiple
         if hasattr(self.move, '__len__'):
-            if len(self.move) == len(self.variables):
-                movevals = self.move
+            # Set movelimit in case of multiple are given
+            move_input = np.asarray(self.move).copy()
+            if move_input.size == len(self.variables):
                 self.move = np.zeros_like(xval)
-                for i in range(len(movevals)):
-                    self.move[self.cumlens[i]:self.cumlens[i + 1]] = movevals[i]
+                for i in range(move_input.size):
+                    self.move[self.cumlens[i]:self.cumlens[i + 1]] = move_input[i]
             elif len(self.move) != self.n:
-                raise RuntimeError("Length of the move vector not correct ({} != {})".format(len(self.move), self.n))
+                raise RuntimeError(f"Length of the move vector ({len(self.move)}) should be equal to number of "
+                                   f"design variable signals ({len(self.variables)}) or "
+                                   f"total number of design variables ({self.n}).")
 
         fcur = 0.0
         while self.iter < self.maxIt:
@@ -333,12 +339,9 @@ def response(self):
             # Set the new states
             for i, s in enumerate(self.variables):
                 if self.cumlens[i+1]-self.cumlens[i] == 1:
-                    try:
-                        s.state[:] = xval[self.cumlens[i]]
-                    except TypeError:
-                        s.state = xval[self.cumlens[i]]
+                    s.state = xval[self.cumlens[i]]
                 else:
-                    s.state[:] = xval[self.cumlens[i]:self.cumlens[i+1]]
+                    s.state = xval[self.cumlens[i]:self.cumlens[i+1]]
 
             if self.fn_callback is not None:
                 self.fn_callback()
@@ -346,19 +349,11 @@ def response(self):
             # Calculate response
             self.funbl.response()
 
-            # Update the states
-            for i, s in enumerate(self.variables):
-                if self.cumlens[i+1]-self.cumlens[i] == 1:
-                    try:
-                        xval[self.cumlens[i]] = s.state[:]
-                    except (TypeError, IndexError):
-                        xval[self.cumlens[i]] = s.state
-                else:
-                    xval[self.cumlens[i]:self.cumlens[i+1]] = s.state[:]
-
             # Save response
             f = ()
             for s in self.responses:
+                if not np.isscalar(s.state):
+                    raise TypeError("State of responses must be scalar.")
                 f += (s.state, )
 
             # Check function change convergence criterion
@@ -388,27 +383,48 @@ def response(self):
                 # Reset sensitivities for the next response
                 self.funbl.reset()
 
-            # Display info on variables
             if self.verbosity >= 3:
+                # Display info on variables
+                show_sensitivities = self.verbosity >= 4
+                msg = ""
                 for i, s in enumerate(self.variables):
-                    isscal = self.cumlens[i + 1] - self.cumlens[i] == 1
-                    msg = "{0:>10s} = ".format(s.tag)
-                    if isscal:
-                        try:
-                            msg += "         {0: .3e}         ".format(s.state)
-                        except TypeError:
-                            msg += "         {0: .3e}         ".format(s.state[0])
+                    if show_sensitivities:
+                        msg += "{0:>10s} = ".format(s.tag[:10])
                     else:
-                        msg += "[{0: .3e} ... {1: .3e}] ".format(min(s.state), max(s.state))
-                    for j, s_out in enumerate(self.responses):
-                        msg += "| {0:>10s}/{1:10s} = ".format("d"+s_out.tag, "d"+s.tag)
-                        if isscal:
-                            msg += "         {0: .3e}         ".format(df[j][self.cumlens[i]])
+                        msg += f"{s.tag} = "
+
+                    # Display value range
+                    fmt = '% .2e'
+                    minval, maxval = np.min(s.state), np.max(s.state)
+                    mintag, maxtag = fmt % minval, fmt % maxval
+                    if mintag == maxtag:
+                        if show_sensitivities:
+                            msg += f"       {mintag}      "
                         else:
-                            msg += "[{0: .3e} ... {1: .3e}] ".format(min(df[j][self.cumlens[i]:self.cumlens[i+1]]), max(df[j][self.cumlens[i]:self.cumlens[i+1]]))
-                    print(msg)
+                            msg += f" {mintag}"
+                    else:
+                        sep = '…' if len(s.state) > 2 else ','
+                        msg += f"[{mintag}{sep}{maxtag}]"
+                        if show_sensitivities:
+                            msg += " "
+
+                    if show_sensitivities:
+                        # Display info on sensivity values
+                        for j, s_out in enumerate(self.responses):
+                            msg += "| {0:s}/{1:11s} = ".format("d" + s_out.tag, "d" + s.tag[:10])
+                            minval = np.min(df[j][self.cumlens[i]:self.cumlens[i+1]])
+                            maxval = np.max(df[j][self.cumlens[i]:self.cumlens[i+1]])
+                            mintag, maxtag = fmt % minval, fmt % maxval
+                            if mintag == maxtag:
+                                msg += f"       {mintag}      "
+                            else:
+                                sep = '…' if self.cumlens[i + 1] - self.cumlens[i] > 2 else ','
+                                msg += f"[{mintag}{sep}{maxtag}] "
+                        msg += '\n'
+                    elif i != len(self.variables)-1:
+                        msg += ', '
+                print(msg)
 
-            self.iter += 1
             xnew, change = self.mmasub(xval.copy(), np.hstack(f), np.vstack(df))
 
             # Stopping criteria on step size
@@ -419,6 +435,7 @@ def response(self):
                 break
 
             xval = xnew
+            self.iter += 1
 
     def mmasub(self, xval, g, dg):
         if self.dx is None:
@@ -535,20 +552,33 @@ def mmasub(self, xval, g, dg):
         change = np.average(abs(xval - xmma))
 
         if self.verbosity >= 2:
+            # Display iteration status message
             msgs = ["g{0:d}({1:s}): {2:+.4e}".format(i, s.tag, g[i]) for i, s in enumerate(self.responses)]
-            print("It. {0: 4d}, {1}".format(self.iter, ", ".join(msgs)))
+            max_infeasibility = max(g[1:])
+            is_feasible = max_infeasibility <= 0
+
+            feasibility_tag = 'f' if is_feasible else ' '
+            print("It. {0: 4d}, [{1:1s}] {2}".format(self.iter, feasibility_tag, ", ".join(msgs)))
 
-        if self.verbosity >=3:
-            # Print changes
-            printstr = "Changes: "
+        if self.verbosity >= 3:
+            # Report design feasibility
+            iconst_max = np.argmax(g[1:])
+            print(f"  | {np.sum(g[1:]>0)} / {len(g)-1} violated constraints, "
+                  f"max. violation ({self.responses[iconst_max+1].tag}) = {'%.2g'%g[iconst_max+1]}")
+
+            # Print design changes
+            change_msgs = []
             for i, s in enumerate(self.variables):
-                isscal = self.cumlens[i + 1] - self.cumlens[i] == 1
-                if isscal:
-                    chg = abs(xval[self.cumlens[i]] - xmma[self.cumlens[i]])
+                minchg = np.min(abs(xval[self.cumlens[i]:self.cumlens[i + 1]] - xmma[self.cumlens[i]:self.cumlens[i + 1]]))
+                maxchg = np.max(abs(xval[self.cumlens[i]:self.cumlens[i + 1]] - xmma[self.cumlens[i]:self.cumlens[i + 1]]))
+                fmt = '%.2g'
+                mintag, maxtag = fmt % minchg, fmt % maxchg
+
+                if mintag == maxtag:
+                    change_msgs.append(f"Δ({s.tag}) = {mintag}")
                 else:
-                    chg = np.average(abs(xval[self.cumlens[i]:self.cumlens[i + 1]] - xmma[self.cumlens[i]:self.cumlens[i + 1]]))
+                    change_msgs.append(f"Δ({s.tag}) = {mintag}…{maxtag}")
 
-                printstr += "{0:s} = {1:.3e}   ".format("Δ_"+s.tag, chg)
-            print(printstr)
+            print(f"  | Changes: {', '.join(change_msgs)}")
 
-        return xmma, change
\ No newline at end of file
+        return xmma, change
diff --git a/pymoto/common/solvers.py b/pymoto/common/solvers.py
deleted file mode 100644
index 4bb5a64..0000000
--- a/pymoto/common/solvers.py
+++ /dev/null
@@ -1,236 +0,0 @@
-import numpy as np
-import scipy.sparse as sps
-try:
-    import cvxopt
-    _has_cvxopt = True
-except ImportError:
-    _has_cvxopt = False
-
-
-def is_cvxopt_spmatrix(A):
-    """ Checks if the argument is a cvxopt sparse matrix """
-    return isinstance(A, cvxopt.spmatrix) if _has_cvxopt else False
-
-
-def matrix_is_complex(A):
-    """ Checks if the matrix is complex """
-    if is_cvxopt_spmatrix(A):
-        return A.typecode == 'z'
-    else:
-        return np.iscomplexobj(A)
-
-
-def matrix_is_diagonal(A):
-    """ Checks if the matrix is diagonal"""
-    if sps.issparse(A):
-        if isinstance(A, sps.dia_matrix):
-            return len(A.offsets) == 1 and A.offsets[0] == 0
-        else:
-            return np.allclose((A - sps.spdiags(A.diagonal(), 0, *A.shape)).data, 0.0)
-    elif is_cvxopt_spmatrix(A):
-        return max(abs(A.I - A.J)) == 0
-    else:
-        return np.allclose(A, np.diag(np.diag(A)))
-
-
-def matrix_is_symmetric(A):
-    """ Checks whether a matrix is numerically symmetric """
-    if sps.issparse(A):
-        return np.allclose((A-A.T).data, 0)
-    elif is_cvxopt_spmatrix(A):
-        return np.isclose(max(abs(A-A.T)), 0.0)
-    else:
-        return np.allclose(A, A.T)
-
-
-def matrix_is_hermitian(A):
-    """ Checks whether a matrix is numerically Hermitian """
-    if matrix_is_complex(A):
-        if sps.issparse(A):
-            return np.allclose((A-A.T.conj()).data, 0)
-        elif is_cvxopt_spmatrix(A):
-            return np.isclose(max(abs(A-A.ctrans())), 0.0)
-        else:
-            return np.allclose(A, A.T.conj())
-    else:
-        return matrix_is_symmetric(A)
-
-
-class LinearSolver:
-    """ Base class of all linear solvers
-
-    Keyword Args:
-        A (matrix): Optionally provide a matrix, which is used in :method:`update` right away.
-
-    Attributes:
-        defined (bool): Flag if the solver is able to run, e.g. false if some dependent library is not available
-    """
-
-    defined = True
-    _err_msg = ""
-
-    def __init__(self, A=None):
-        if A is not None:
-            self.update(A)
-
-    def update(self, A):
-        """ Updates with a new matrix of the same structure
-
-        Args:
-            A (matrix): The new matrix of size ``(N, N)``
-
-        Returns:
-            self
-        """
-        raise NotImplementedError(f"Solver not implemented {self._err_msg}")
-
-    def solve(self, rhs):
-        r""" Solves the linear system of equations :math:`\mathbf{A} \mathbf{x} = \mathbf{b}`
-
-        Args:
-            rhs: Right hand side :math:`\mathbf{b}` of shape ``(N)`` or ``(N, K)`` for multiple right-hand-sides
-
-        Returns:
-            Solution vector :math:`\mathbf{x}` of same shape as :math:`\mathbf{b}`
-        """
-        raise NotImplementedError(f"Solver not implemented {self._err_msg}")
-
-    def adjoint(self, rhs):
-        r""" Solves the adjoint linear system of equations
-
-        The system of equations is :math:`\mathbf{A}^\text{H} \mathbf{x} = \mathbf{b}` (conjugate transpose) in case of
-        complex matrix or :math:`\mathbf{A}^\text{T} \mathbf{x} = \mathbf{b}` for a real-valued matrix.
-
-        Args:
-            rhs: Right hand side :math:`\mathbf{b}` of shape ``(N)`` or ``(N, K)`` for multiple right-hand-sides
-
-        Returns:
-            Solution vector :math:`\mathbf{x}` of same shape as :math:`\mathbf{b}`
-        """
-        raise NotImplementedError(f"Solver not implemented {self._err_msg}")
-
-    @staticmethod
-    def residual(A, x, b):
-        r""" Calculates the (relative) residual of the linear system of equations
-
-        The residual is calculated as
-        :math:`r = \frac{\left| \mathbf{A} \mathbf{x} - \mathbf{b} \right|}{\left| \mathbf{b} \right|}`
-
-        Args:
-            A: The matrix
-            x: Solution vector
-            b: Right-hand side
-
-        Returns:
-            Residual value
-        """
-        return np.linalg.norm(A@x - b) / np.linalg.norm(b)
-
-
-class LDAWrapper(LinearSolver):
-    r""" Linear dependency aware solver (LDAS)
-
-    This solver uses previous solutions of the system :math:`\mathbf{A} \mathbf{x} = \mathbf{b}` to reduce computational
-    effort. In case the solution :math:`\mathbf{x}` is linearly dependent on the previous solutions, the solution
-    will be nearly free of cost.
-
-    Args:
-        solver: The internal solver to be used
-        tol (optional): Residual tolerance above which the internal solver is used to add a new solution vector.
-        A (optional): The matrix :math:`\mathbf{A}`
-
-    References:
-
-    Koppen, S., van der Kolk, M., van den Boom, S., & Langelaar, M. (2022).
-    Efficient computation of states and sensitivities for compound structural optimisation problems using a Linear Dependency Aware Solver (LDAS).
-    Structural and Multidisciplinary Optimization, 65(9), 273.
-    DOI: 10.1007/s00158-022-03378-8
-    """
-    def __init__(self, solver: LinearSolver, tol=1e-8, A=None, hermitian=False, symmetric=False):
-        self.solver = solver
-        self.tol = tol
-        self.x_stored = []
-        self.b_stored = []
-        self.xadj_stored = []
-        self.badj_stored = []
-        self.A = None
-        self._did_solve = False  # For debugging purposes
-        self._last_rtol = 0.
-        self.hermitian = hermitian
-        self.symmetric = symmetric
-        super().__init__(A)
-
-    def update(self, A):
-        """ Clear the internal stored solution vectors and update the internal ``solver`` """
-        self.A = A
-        self.x_stored.clear()
-        self.b_stored.clear()
-        self.xadj_stored.clear()
-        self.badj_stored.clear()
-        self.solver.update(A)
-
-    def _do_solve_1rhs(self, A, rhs, x_data, b_data, solve_fn):
-        rhs_loc = rhs.copy()
-        sol = 0
-
-        # Check linear dependencies in the rhs using modified Gram-Schmidt
-        for (x, b) in zip(x_data, b_data):
-            alpha = rhs_loc.conj() @ b / (b.conj() @ b)
-            rhs_loc -= alpha * b
-            sol += alpha * x
-
-        # Check tolerance
-        self._last_rtol = 1.0 if len(x_data) == 0 else self.residual(A, sol, rhs)
-
-        if self._last_rtol > self.tol:
-            # Calculate a new solution
-            xnew = solve_fn(rhs_loc)
-            x_data.append(xnew)
-            b_data.append(rhs_loc)
-            sol += xnew
-            self._did_solve = True
-        else:
-            self._did_solve = False
-
-        return sol
-
-    def _solve_1x(self, b):
-        return self._do_solve_1rhs(self.A, b, self.x_stored, self.b_stored, self.solver.solve)
-
-    def _adjoint_1x(self, b):
-        return self._do_solve_1rhs(self.A.conj().T, b, self.xadj_stored, self.badj_stored, self.solver.adjoint)
-
-    def solve(self, rhs):
-        r""" Solves the linear system of equations :math:`\mathbf{A} \mathbf{x} = \mathbf{b}` by performing a modified
-        Gram-Schmidt over the previously calculated solutions :math:`\mathbf{U}` and corresponding right-hand-sides
-        :math:`\mathbf{F}`. This is used to construct an approximate solution
-        :math:`\tilde{\mathbf{x}} = \sum_k \alpha_k \mathbf{u}_k` in the subspace of :math:`\mathbf{U}`.
-        If the residual of :math:`\mathbf{A} \tilde{\mathbf{x}} = \mathbf{b}` is above the tolerance, a new solution
-        :math:`\mathbf{u}_{k+1}` will be added to the database such that
-        :math:`\mathbf{x} = \tilde{\mathbf{x}}+\mathbf{u}_{k+1}` is the solution to the system
-        :math:`\mathbf{A} \mathbf{x} = \mathbf{b}`.
-
-        The right-hand-side :math:`\mathbf{b}` can be of size ``(N)`` or ``(N, K)``, where ``N`` is the size of matrix
-        :math:`\mathbf{A}` and ``K`` is the number of right-hand sides.
-        """
-        if rhs.ndim == 1:
-            return self._solve_1x(rhs)
-        else:  # Multiple rhs
-            sol = []
-            for i in range(rhs.shape[-1]):
-                sol.append(self._solve_1x(rhs[..., i]))
-            return np.stack(sol, axis=-1)
-
-    def adjoint(self, rhs):
-        if self.hermitian:
-            return self.solve(rhs)
-        elif self.symmetric:
-            return self.solve(rhs.conj()).conj()
-        else:
-            if rhs.ndim == 1:
-                return self._adjoint_1x(rhs)
-            else:  # Multiple rhs
-                sol = []
-                for i in range(rhs.shape[-1]):
-                    sol.append(self._adjoint_1x(rhs[..., i]))
-                return np.stack(sol, axis=-1)
diff --git a/pymoto/core_objects.py b/pymoto/core_objects.py
index dc07632..2755e48 100644
--- a/pymoto/core_objects.py
+++ b/pymoto/core_objects.py
@@ -1,9 +1,11 @@
-from typing import Union, List, Any
+import sys
 import warnings
 import inspect
 import time
-from .utils import _parse_to_list, _concatenate_to_array, _split_from_array
+import copy
+from typing import Union, List, Any
 from abc import ABC, abstractmethod
+from .utils import _parse_to_list, _concatenate_to_array, _split_from_array
 
 
 # Local helper functions
@@ -11,7 +13,7 @@ def err_fmt(*args):
     """ Format error strings for locating Modules and Signals"""
     err_str = ""
     for a in args:
-        err_str += f"\n\t[ {a} ]"
+        err_str += f"\n\t| {a}"
     return err_str
 
 
@@ -76,16 +78,20 @@ class Signal:
     >> Signal(tag='x2')
 
     """
-    def __init__(self, tag: str = "", state: Any = None, sensitivity: Any = None):
+    def __init__(self, tag: str = "", state: Any = None, sensitivity: Any = None, min: Any = None, max: Any = None):
         """
-
-        :param tag: The name of the signal (string)
-        :param state: The initialized state (optional)
-        :param sensitivity: The initialized sensitivity (optional)
+        Keyword Args:
+            tag: The name of the signal
+            state: The initialized state
+            sensitivity: The initialized sensitivity
+            min: Minimum allowed value
+            max: Maximum allowed value
         """
         self.tag = tag
         self.state = state
         self.sensitivity = sensitivity
+        self.min = min
+        self.max = max
         self.keep_alloc = sensitivity is not None
 
         # Save error string to location where it is initialized
@@ -95,11 +101,12 @@ def _err_str(self):
         return err_fmt(f"Signal \'{self.tag}\', initialized in {self._init_loc}")
 
     def add_sensitivity(self, ds: Any):
+        """ Add a new term to internal sensitivity """
         try:
             if ds is None:
                 return
             if self.sensitivity is None:
-                self.sensitivity = ds
+                self.sensitivity = copy.deepcopy(ds)
             else:
                 self.sensitivity += ds
             return self
@@ -116,8 +123,12 @@ def add_sensitivity(self, ds: Any):
     def reset(self, keep_alloc: bool = None):
         """ Reset the sensitivities to zero or None
         This must be called to clear internal memory of subsequent sensitivity calculations.
-        :param keep_alloc: Keep the sensitivity allocation intact?
-        :return: self
+
+        Args:
+            keep_alloc: Keep the sensitivity allocation intact?
+
+        Returns:
+            self
         """
         if self.sensitivity is None:
             return self
@@ -138,11 +149,34 @@ def reset(self, keep_alloc: bool = None):
 
     def __getitem__(self, item):
         """ Obtain a sliced signal, for using its partial contents.
-        :param item: Slice indices
-        :return: Sliced signal (SignalSlice)
+
+        Args:
+            item: Slice indices
+
+        Returns:
+            Sliced signal (SignalSlice)
         """
         return SignalSlice(self, item)
 
+    def __str__(self):
+        state_msg = f"state {self.state}" if self.state is not None else "empty state"
+        state_msg = state_msg.split('\n')
+        if len(state_msg) > 1:
+            state_msg = state_msg[0] + ' ... ' + state_msg[-1]
+        else:
+            state_msg = state_msg[0]
+        return f"Signal \"{self.tag}\" with {state_msg}"
+
+    def __repr__(self):
+        state_msg = f"state {self.state}" if self.state is not None else "empty state"
+        state_msg = state_msg.split('\n')
+        if len(state_msg) > 1:
+            state_msg = state_msg[0] + ' ... ' + state_msg[-1]
+        else:
+            state_msg = state_msg[0]
+        sens_msg = 'empty sensitivity' if self.sensitivity is None else 'non-empty sensitivity'
+        return f"Signal \"{self.tag}\" with {state_msg} and {sens_msg} at {hex(id(self))}"
+
 
 class SignalSlice(Signal):
     """ Slice operator for a Signal
@@ -169,7 +203,8 @@ def state(self):
             return None if self.orig_signal.state is None else self.orig_signal.state[self.slice]
         except Exception as e:
             # Possibilities: Unslicable object (TypeError) or Wrong dimensions or out of range (IndexError)
-            raise type(e)("SignalSlice.state (getter)" + self._err_str()) from e
+            raise type(e)(str(e) + "\n\t| Above error was raised in SignalSlice.state (getter). Signal details:" +
+                      self._err_str()).with_traceback(sys.exc_info()[2])
 
     @state.setter
     def state(self, new_state):
@@ -177,7 +212,8 @@ def state(self, new_state):
             self.orig_signal.state[self.slice] = new_state
         except Exception as e:
             # Possibilities: Unslicable object (TypeError) or Wrong dimensions or out of range (IndexError)
-            raise type(e)("SignalSlice.state (setter)" + self._err_str()) from e
+            raise type(e)(str(e) + "\n\t| Above error was raised in SignalSlice.state (setter). Signal details:" +
+                          self._err_str()).with_traceback(sys.exc_info()[2])
 
     @property
     def sensitivity(self):
@@ -185,7 +221,8 @@ def sensitivity(self):
             return None if self.orig_signal.sensitivity is None else self.orig_signal.sensitivity[self.slice]
         except Exception as e:
             # Possibilities: Unslicable object (TypeError) or Wrong dimensions or out of range (IndexError)
-            raise type(e)("SignalSlice.sensitivity (getter)" + self._err_str()) from e
+            raise type(e)(str(e) + "\n\t| Above error was raised in SignalSlice.sensitivity (getter). Signal details:" +
+                          self._err_str()).with_traceback(sys.exc_info()[2])
 
     @sensitivity.setter
     def sensitivity(self, new_sens):
@@ -207,7 +244,8 @@ def sensitivity(self, new_sens):
             self.orig_signal.sensitivity[self.slice] = new_sens
         except Exception as e:
             # Possibilities: Unslicable object (TypeError) or Wrong dimensions or out of range (IndexError)
-            raise type(e)("SignalSlice.sensitivity (setter)" + self._err_str()) from e
+            raise type(e)(str(e) + "\n\t| Above error was raised in SignalSlice.state (setter). Signal details:" +
+                          self._err_str()).with_traceback(sys.exc_info()[2])
 
     def reset(self, keep_alloc: bool = None):
         """ Reset the sensitivities to zero or None
@@ -231,7 +269,7 @@ def make_signals(*args):
     return ret
 
 
-def _check_valid_signal(sig: Any):
+def _is_valid_signal(sig: Any):
     """ Checks if the argument is a valid Signal object
     :param sig: The object to check
     :return: True if it is a valid Signal
@@ -240,10 +278,10 @@ def _check_valid_signal(sig: Any):
         return True
     if all([hasattr(sig, f) for f in ["state", "sensitivity", "add_sensitivity", "reset"]]):
         return True
-    raise TypeError(f"Given argument with type \'{type(sig).__name__}\' is not a valid Signal")
+    return False
 
 
-def _check_valid_module(mod: Any):
+def _is_valid_module(mod: Any):
     """ Checks if the argument is a valid Module object
     :param mod: The object to check
     :return: True if it is a valid Module
@@ -252,7 +290,7 @@ def _check_valid_module(mod: Any):
         return True
     if hasattr(mod, "response") and hasattr(mod, "sensitivity") and hasattr(mod, "reset"):
         return True
-    raise TypeError(f"Given argument with type \'{type(mod).__name__}\' is not a valid Module")
+    return False
 
 
 def _check_function_signature(fn, signals):
@@ -370,60 +408,48 @@ class Module(ABC, RegisteredClass):
     >> Module(sig_in=[inputs], sig_out=[outputs]
     """
 
-    def _err_str(self, init: bool = True, add_signal: bool = True, fn=None):
+    def _err_str(self, module_signature: bool = True, init: bool = True, fn=None):
         str_list = []
-        if init:
-            str_list.append(f"Module \'{type(self).__name__}\', initialized in {self._init_loc}")
-        if add_signal:
+
+        if module_signature:
             inp_str = "Inputs: " + ", ".join([s.tag if hasattr(s, 'tag') else 'N/A' for s in self.sig_in]) if len(self.sig_in) > 0 else "No inputs"
             out_str = "Outputs: " + ", ".join([s.tag if hasattr(s, 'tag') else 'N/A' for s in self.sig_out]) if len(self.sig_out) > 0 else "No outputs"
-            str_list.append(inp_str + " --> " + out_str)
+            str_list.append(f"Module \'{type(self).__name__}\'( " + inp_str + " ) --> " + out_str)
+        if init:
+            str_list.append(f"Used in {self._init_loc}")
         if fn is not None:
             name = f"{fn.__self__.__class__.__name__}.{fn.__name__}{inspect.signature(fn)}"
             lineno = inspect.getsourcelines(fn)[1]
             filename = inspect.getfile(fn)
-            str_list.append(f"Implemented in File \"{filename}\", line {lineno}, in {name}")
+            str_list.append(f"Implementation in File \"{filename}\", line {lineno}, in {name}")
         return err_fmt(*str_list)
 
     # flake8: noqa: C901
     def __init__(self, sig_in: Union[Signal, List[Signal]] = None, sig_out: Union[Signal, List[Signal]] = None,
                  *args, **kwargs):
-        # TODO: Reduce complexity of this init
         self._init_loc = get_init_str()
 
         self.sig_in = _parse_to_list(sig_in)
         self.sig_out = _parse_to_list(sig_out)
         for i, s in enumerate(self.sig_in):
-            try:
-                _check_valid_signal(s)
-            except Exception as e:
-                earg0 = e.args[0] if len(e.args) > 0 else ''
-                earg1 = e.args[1:] if len(e.args) > 1 else ()
-                raise type(e)(f"Invalid input signal #{i+1} - " + str(earg0) + self._err_str(), *earg1) from None
+            if not _is_valid_signal(s):
+                tag = f" (\'{s.tag}\')" if hasattr(s, 'tag') else ''
+                raise TypeError(f"Input {i}{tag} is not a valid signal, type=\'{type(s).__name__}\'.")
 
         for i, s in enumerate(self.sig_out):
-            try:
-                _check_valid_signal(s)
-            except Exception as e:
-                earg0 = e.args[0] if len(e.args) > 0 else ''
-                earg1 = e.args[1:] if len(e.args) > 1 else ()
-                raise type(e)(f"Invalid output signal #{i+1} - " + str(earg0) + self._err_str(), *earg1) from None
+            if not _is_valid_signal(s):
+                tag = f" (\'{s.tag}\')" if hasattr(s, 'tag') else ''
+                raise TypeError(f"Output {i}{tag} is not a valid signal, type=\'{type(s).__name__}\'.")
 
-        try:
-            # Call preparation of submodule with remaining arguments
-            self._prepare(*args, **kwargs)
-        except Exception as e:
-            earg0 = e.args[0] if len(e.args) > 0 else ''
-            earg1 = e.args[1:] if len(e.args) > 1 else ()
-            raise type(e)("_prepare() - " + str(earg0) + self._err_str(fn=self._prepare), *earg1) from e
+        # Call preparation of submodule with remaining arguments
+        self._prepare(*args, **kwargs)
 
         try:
             # Check if the signals match _response() signature
             _check_function_signature(self._response, self.sig_in)
         except Exception as e:
-            earg0 = e.args[0] if len(e.args) > 0 else ''
-            earg1 = e.args[1:] if len(e.args) > 1 else ()
-            raise type(e)(str(earg0) + self._err_str(fn=self._response), *earg1) from None
+            raise type(e)(str(e) + "\n\t| Module details:" +
+                          self._err_str(fn=self._response)).with_traceback(sys.exc_info()[2])
 
         try:
             # If no output signals are given, but are required, try to initialize them here
@@ -441,9 +467,8 @@ def __init__(self, sig_in: Union[Signal, List[Signal]] = None, sig_out: Union[Si
             # Check if signals match _sensitivity() signature
             _check_function_signature(self._sensitivity, self.sig_out)
         except Exception as e:
-            earg0 = e.args[0] if len(e.args) > 0 else ''
-            earg1 = e.args[1:] if len(e.args) > 1 else ()
-            raise type(e)(str(earg0) + self._err_str(fn=self._sensitivity), *earg1) from None
+            raise type(e)(str(e) + "\n\t| Module details:" +
+                          self._err_str(fn=self._sensitivity)).with_traceback(sys.exc_info()[2])
 
     def response(self):
         """ Calculate the response from sig_in and output this to sig_out """
@@ -461,9 +486,9 @@ def response(self):
                 self.sig_out[i].state = val
             return self
         except Exception as e:
-            earg0 = e.args[0] if len(e.args) > 0 else ''
-            earg1 = e.args[1:] if len(e.args) > 1 else ()
-            raise type(e)("response() - " + str(earg0) + self._err_str(fn=self._response), *earg1) from e
+            # https://stackoverflow.com/questions/6062576/adding-information-to-an-exception
+            raise type(e)(str(e) + "\n\t| Above error was raised when calling response(). Module details:" +
+                          self._err_str(fn=self._response)).with_traceback(sys.exc_info()[2])
 
     def __call__(self):
         return self.response()
@@ -494,9 +519,8 @@ def sensitivity(self):
 
             return self
         except Exception as e:
-            earg0 = e.args[0] if len(e.args) > 0 else ''
-            earg1 = e.args[1:] if len(e.args) > 1 else ()
-            raise type(e)("sensitivity() - " + str(earg0) + self._err_str(fn=self._sensitivity), *earg1) from e
+            raise type(e)(str(e) + "\n\t| Above error was raised when calling sensitivity(). Module details:" +
+                          self._err_str(fn=self._sensitivity)).with_traceback(sys.exc_info()[2])
 
     def reset(self):
         """ Reset the state of the sensitivities (they are set to zero or to None) """
@@ -506,9 +530,8 @@ def reset(self):
             self._reset()
             return self
         except Exception as e:
-            earg0 = e.args[0] if len(e.args) > 0 else ''
-            earg1 = e.args[1:] if len(e.args) > 1 else ()
-            raise type(e)("reset() - " + str(earg0) + self._err_str(fn=self._reset), *earg1) from e
+            raise type(e)(str(e) + "\n\t| Above error was raised when calling reset(). Module details:" +
+                          self._err_str(fn=self._response)).with_traceback(sys.exc_info()[2])
 
     # METHODS TO BE DEFINED BY USER
     def _prepare(self, *args, **kwargs):
@@ -542,35 +565,33 @@ class Network(Module):
     """
     def __init__(self, *args, print_timing=False):
         self._init_loc = get_init_str()
-        try:
-            # Obtain the internal blocks
-            self.mods = _parse_to_list(*args)
-
-            # Check if the blocks are initialized, else create them
-            for i, b in enumerate(self.mods):
-                if isinstance(b, dict):
-                    exclude_keys = ['type']
-                    b_ex = {k: b[k] for k in set(list(b.keys())) - set(exclude_keys)}
-                    self.mods[i] = Module.create(b['type'], **b_ex)
-
-            # Check validity of modules
-            [_check_valid_module(m) for m in self.mods]
-
-            # Gather all the input and output signals of the internal blocks
-            all_in = set()
-            all_out = set()
-            [all_in.update(b.sig_in) for b in self.mods]
-            [all_out.update(b.sig_out) for b in self.mods]
-            in_unique = all_in - all_out
-
-            # Initialize the parent module, with correct inputs and outputs
-            super().__init__(list(in_unique), list(all_out))
-
-            self.print_timing = print_timing
-        except Exception as e:
-            earg0 = e.args[0] if len(e.args) > 0 else ''
-            earg1 = e.args[1:] if len(e.args) > 1 else ()
-            raise type(e)(str(earg0) + self._err_str(add_signal=False), *earg1) from None
+
+        # Obtain the internal blocks
+        self.mods = _parse_to_list(*args)
+
+        # Check if the blocks are initialized, else create them
+        for i, b in enumerate(self.mods):
+            if isinstance(b, dict):
+                exclude_keys = ['type']
+                b_ex = {k: b[k] for k in set(list(b.keys())) - set(exclude_keys)}
+                self.mods[i] = Module.create(b['type'], **b_ex)
+
+        # Check validity of modules
+        for m in self.mods:
+            if not _is_valid_module(m):
+                raise TypeError(f"Argument is not a valid Module, type=\'{type(mod).__name__}\'.")
+
+        # Gather all the input and output signals of the internal blocks
+        all_in = set()
+        all_out = set()
+        [all_in.update(b.sig_in) for b in self.mods]
+        [all_out.update(b.sig_out) for b in self.mods]
+        in_unique = all_in - all_out
+
+        # Initialize the parent module, with correct inputs and outputs
+        super().__init__(list(in_unique), list(all_out))
+
+        self.print_timing = print_timing
 
     def timefn(self, fn):
         start_t = time.time()
@@ -611,12 +632,9 @@ def append(self, *newmods):
         modlist = _parse_to_list(*newmods)
 
         # Check if the blocks are initialized, else create them
-        for i, b in enumerate(modlist):
-            try:  # Check validity of modules
-                _check_valid_module(b)
-            except Exception as e:
-                raise type(e)("append() - Trying to append invalid module " + str(e.args[0])
-                              + self._err_str(add_signal=False), *e.args[1:]) from None
+        for i, m in enumerate(modlist):
+            if not _is_valid_module(m):
+                raise TypeError(f"Argument #{i} is not a valid module, type=\'{type(mod).__name__}\'.")
 
         # Obtain the internal blocks
         self.mods.extend(modlist)
@@ -624,25 +642,11 @@ def append(self, *newmods):
         # Gather all the input and output signals of the internal blocks
         all_in = set()
         all_out = set()
-        [all_in.update(b.sig_in) for b in self.mods]
-        [all_out.update(b.sig_out) for b in self.mods]
+        [all_in.update(m.sig_in) for m in self.mods]
+        [all_out.update(m.sig_out) for m in self.mods]
         in_unique = all_in - all_out
 
         self.sig_in = _parse_to_list(in_unique)
-        try:
-            [_check_valid_signal(s) for s in self.sig_in]
-        except Exception as e:
-            earg0 = e.args[0] if len(e.args) > 0 else ''
-            earg1 = e.args[1:] if len(e.args) > 1 else ()
-            raise type(e)("append() - Invalid input signals " + str(earg0)
-                          + self._err_str(add_signal=False), *earg1) from None
         self.sig_out = _parse_to_list(all_out)
-        try:
-            [_check_valid_signal(s) for s in self.sig_out]
-        except Exception as e:
-            earg0 = e.args[0] if len(e.args) > 0 else ''
-            earg1 = e.args[1:] if len(e.args) > 1 else ()
-            raise type(e)("append() - Invalid output signals " + str(earg0)
-                          + self._err_str(add_signal=False), *earg1) from None
 
         return modlist[-1].sig_out[0] if len(modlist[-1].sig_out) == 1 else modlist[-1].sig_out  # Returns the output signal
diff --git a/pymoto/modules/aggregation.py b/pymoto/modules/aggregation.py
new file mode 100644
index 0000000..e7fb367
--- /dev/null
+++ b/pymoto/modules/aggregation.py
@@ -0,0 +1,207 @@
+import warnings
+import abc
+import numpy as np
+import scipy.special as spsp
+from pymoto import Module
+
+
+class AggActiveSet:
+    """ Determine active set by discarding lower or upper fraction of a set of values
+
+    Args:
+       lower_rel: Fraction of values closest to minimum to discard (based on value)
+       upper_rel: Fraction of values closest to maximum to discard (based on value)
+       lower_amt: Fraction of lowest values to discard (based on sorting)
+       upper_amt: Fraction of highest values to discard (based on sorting)
+    """
+    def __init__(self, lower_rel=0.0, upper_rel=1.0, lower_amt=0.0, upper_amt=1.0):
+        assert upper_rel > lower_rel, "Upper must be larger than lower to keep values in the set"
+        assert upper_amt > lower_amt, "Upper must be larger than lower to keep values in the set"
+        self.lower_rel, self.upper_rel = lower_rel, upper_rel
+        self.lower_amt, self.upper_amt = lower_amt, upper_amt
+
+    def __call__(self, x):
+        """ Generate an active set for given array """
+        xmin, xmax = np.min(x), np.max(x)
+        if (xmax - xmin) == 0:  # All values are the same, so no active set can be taken
+            return Ellipsis
+
+        sel = np.ones_like(x, dtype=bool)
+
+        # Select based on value
+        xrel = (x - xmin) / (xmax - xmin)  # Normalize between 0 and 1
+        if self.lower_rel > 0:
+            sel = np.logical_and(sel, xrel >= self.lower_rel)
+        if self.upper_rel < 1:
+            sel = np.logical_and(sel, xrel <= self.upper_rel)
+
+        # Remove lowest and highest N values
+        i_sort = np.argsort(x)
+        if self.lower_amt > 0:
+            n_lower_amt = int(x.size * self.lower_amt)
+            sel[i_sort[:n_lower_amt]] = False
+
+        if self.upper_amt < 1:
+            n_upper_amt = int(x.size * (1 - self.upper_amt))
+            sel[i_sort[-n_upper_amt:]] = False
+
+        return sel
+
+
+class AggScaling:
+    """ Scaling strategy to absolute minimum or maximum
+
+    Args:
+        which: Scale to `min` or `max`
+    """
+    def __init__(self, which: str, damping=0.0):
+        self.damping = damping
+        if which.lower() == 'min':
+            self.f = np.min
+        elif which.lower() == 'max':
+            self.f = np.max
+        else:
+            raise ValueError("Argument `which` can only be 'min' or 'max'")
+        self.sf = None
+
+    def __call__(self, x, fx_approx):
+        """ Determine scaling factor
+
+        Args:
+            x: Set of values
+            fx_approx: Approximated minimum / maximum
+
+        Returns:
+            Scaling factor
+        """
+        trueval = self.f(x)
+        scale = trueval / fx_approx
+        if self.sf is None:
+            self.sf = scale
+        else:
+            self.sf = self.damping * self.sf + (1 - self.damping) * scale
+        return self.sf
+
+
+class Aggregation(Module):
+    """ Generic Aggregation module (cannot be used directly, but can only be used as superclass)
+
+    Keyword Args:
+        scaling: Scaling strategy to improve approximation
+        active_set: Active set strategy to improve approximation
+    """
+    def _prepare(self, scaling: AggScaling = None, active_set: AggActiveSet = None):
+        # This prepare function MUST be called in the _prepare function of sub-classes
+        self.scaling = scaling
+        self.active_set = active_set
+        self.sf = 1.0
+
+    @abc.abstractmethod
+    def aggregation_function(self, x):
+        """ Calculates f(x) """
+        raise NotImplementedError()
+
+    @abc.abstractmethod
+    def aggregation_derivative(self, x):
+        """" Calculates df(x) / dx """
+        raise NotImplementedError()
+
+    def _response(self, x):
+        # Determine active set
+        if self.active_set is not None:
+            self.select = self.active_set(x)
+        else:
+            self.select = Ellipsis
+
+        # Get aggregated value
+        xagg = self.aggregation_function(x[self.select])
+
+        # Scale
+        if self.scaling is not None:
+            self.sf = self.scaling(x[self.select], xagg)
+        return self.sf * xagg
+
+    def _sensitivity(self, dfdy):
+        x = self.sig_in[0].state
+        dydx = self.aggregation_derivative(x[self.select])
+        dx = np.zeros_like(x)
+        dx[self.select] += self.sf * dfdy * dydx
+        return dx
+
+
+class PNorm(Aggregation):
+    r""" P-norm aggregration
+
+    :math:`S_p(x_1, x_2, \dotsc, x_n) = \left( \sum_i (|x_i|^p) \right)^{1/p}
+
+    Only valid for positive :math:`x_i` when approximating the minimum or maximum
+
+    Args:
+        p: Power of the p-norm. Approximate maximum for `p>0` and minimum for `p<0`
+        scaling: Scaling strategy to improve approximation
+        active_set: Active set strategy to improve approximation
+    """
+    def _prepare(self, p=2, scaling: AggScaling = None, active_set: AggActiveSet = None):
+        self.p = p
+        self.y = None
+        super()._prepare(scaling, active_set)
+
+    def aggregation_function(self, x):
+        if np.min(x) < 0:
+            warnings.warn("PNorm is only valid for positive x")
+
+        # Get p-norm
+        return np.sum(np.abs(x) ** self.p) ** (1/self.p)
+
+    def aggregation_derivative(self, x):
+        pval = np.sum(np.abs(x) ** self.p) ** (1 / self.p - 1)
+        return pval * np.sign(x) * np.abs(x)**(self.p - 1)
+
+
+class SoftMinMax(Aggregation):
+    r""" Soft maximum/minimum function
+
+    :math:`S_a(x_1, x_2, \dotsc, x_n) = \frac{\sum_i (x_i \exp(a x_i))}{\sum_i (\exp(a x_i))}`
+
+    When using as maximum, it underestimates the maximum
+    It is exact however when :math:`x_1=x_2=\dotsc=x_n`
+
+    Args:
+        alpha: Scaling factor of the soft function. Approximate maximum for `alpha>0` and minimum for `alpha<0`
+        scaling: Scaling strategy to improve approximation
+        active_set: Active set strategy to improve approximation
+    """
+    def _prepare(self, alpha=1.0, scaling: AggScaling = None, active_set: AggActiveSet = None):
+        self.alpha = alpha
+        self.y = None
+        super()._prepare(scaling, active_set)
+
+    def aggregation_function(self, x):
+        self.y = np.sum(x * spsp.softmax(self.alpha * x))
+        return self.y
+
+    def aggregation_derivative(self, x):
+        return spsp.softmax(self.alpha * x) * (1 + self.alpha * (x - self.y))
+
+
+class KSFunction(Aggregation):
+    r""" Kreisselmeier and Steinhauser function from 1979
+
+    :math:`S_\rho(x_1, x_2, \dotsc, x_n) = \frac{1}{\rho} \ln \left( \sum_i \exp(\rho x_i) \right)`
+
+    Args:
+        rho: Scaling factor of the KS function. Approximate maximum for `alpha>0` and minimum for `alpha<0`
+        scaling: Scaling strategy to improve approximation
+        active_set: Active set strategy to improve approximation
+    """
+    def _prepare(self, rho=1.0, scaling: AggScaling = None, active_set: AggActiveSet = None):
+        self.rho = rho
+        self.y = None
+        super()._prepare(scaling, active_set)
+
+    def aggregation_function(self, x):
+        return 1/self.rho * np.log(np.sum(np.exp(self.rho * x)))
+
+    def aggregation_derivative(self, x):
+        erx = np.exp(self.rho * x)
+        return erx / np.sum(erx)
diff --git a/pymoto/modules/assembly.py b/pymoto/modules/assembly.py
index af05006..bf7d936 100644
--- a/pymoto/modules/assembly.py
+++ b/pymoto/modules/assembly.py
@@ -49,8 +49,8 @@ def _prepare(self, domain: DomainDefinition, element_matrix: np.ndarray, bc=None
         self.dofconn = domain.get_dofconnectivity(self.ndof)
 
         # Row and column indices for the matrix
-        self.rows = np.kron(self.dofconn, np.ones((domain.elemnodes*self.ndof, 1), dtype=int)).flatten()
-        self.cols = np.kron(self.dofconn, np.ones((1, domain.elemnodes*self.ndof), dtype=int)).flatten()
+        self.rows = np.kron(self.dofconn, np.ones((1, domain.elemnodes*self.ndof), dtype=int)).flatten()
+        self.cols = np.kron(self.dofconn, np.ones((domain.elemnodes * self.ndof, 1), dtype=int)).flatten()
         self.matrix_type = matrix_type
 
         # Boundary conditions
@@ -68,6 +68,8 @@ def _prepare(self, domain: DomainDefinition, element_matrix: np.ndarray, bc=None
         self.add_constant = add_constant
 
     def _response(self, xscale: np.ndarray):
+        nel = self.dofconn.shape[0]
+        assert xscale.size == nel, f"Input vector wrong size ({xscale.size}), must be of size #nel ({nel})"
         scaled_el = ((self.elmat.flatten()[np.newaxis]).T * xscale).flatten(order='F')
 
         # Set boundary conditions
@@ -305,3 +307,116 @@ def _prepare(self, domain: DomainDefinition, *args, material_property: float = 1
             self.poisson_element += w * self.material_property * Bn.T @ Bn  # Add contribution
 
         super()._prepare(domain, self.poisson_element, *args, **kwargs)
+
+
+class ElementOperation(Module):
+    r""" Generic module for element-wise operations based on nodal information
+
+    :math:`y_e = \mathbf{B} \mathbf{u}_e`
+
+    Input Signal:
+        - ``u``: Nodal vector of size ``(n_dof_per_element * #nodes)``
+
+    Output Signal:
+        - ``y``: Elemental output data of size ``(..., #elements)``
+
+    Args:
+        domain: The domain defining element and nodal connectivity
+        element_matrix: The element operator matrix :math:`\mathbf{B}` of size ``(..., n_dof_per_element)``
+    """
+    def _prepare(self, domain: DomainDefinition, element_matrix: np.ndarray):
+        if element_matrix.shape[-1] % domain.elemnodes != 0:
+            raise IndexError("Size of last dimension of element operator matrix is not compatible with mesh. "
+                             "Must be dividable by the number of nodes.")
+
+        ndof = element_matrix.shape[-1] // domain.elemnodes
+
+        self.element_matrix = element_matrix
+        self.dofconn = domain.get_dofconnectivity(ndof)
+        self.usiz = ndof * domain.nnodes
+
+    def _response(self, u):
+        assert u.size == self.usiz
+        return einsum('...k, lk -> ...l', self.element_matrix, u[self.dofconn], optimize=True)
+
+    def _sensitivity(self, dy):
+        du_el = einsum('...k, ...l -> lk', self.element_matrix, dy, optimize=True)
+        du = np.zeros_like(self.sig_in[0].state)
+        np.add.at(du, self.dofconn, du_el)
+        return du
+
+
+class Strain(ElementOperation):
+    r""" Evaluate average mechanical strains in solid elements based on deformation
+
+    The strains are returned in Voigt notation.
+    :math:`\mathbf{\epsilon}_e = \mathbf{B} \mathbf{u}_e`
+
+    Each integration point in the element has different strain values. Here, the average is returned.
+
+    The returned strain is either
+    :math:`\mathbf{\epsilon} = \begin{bmatrix}\epsilon_{xx} & \epsilon_{yy} & \epsilon_{xy} \end{bmatrix}`
+    in case ``voigt = False`` or
+    :math:`\mathbf{\epsilon} = \begin{bmatrix}\epsilon_{xx} & \epsilon_{yy} & \gamma_{xy} \end{bmatrix}`
+    in case ``voigt = True``, for which :math:`\gamma_{xy}=2\epsilon_{xy}`.
+
+    Input Signal:
+       - ``u``: Nodal vector of size ``(#dof per element * #nodes)``
+
+    Output Signal:
+       - ``e``: Strain matrix of size ``(#strains per element, #elements)``
+
+    Args:
+        domain: The domain defining element and nodal connectivity
+
+    Keyword Args:
+        voigt: Use Voigt strain notation (2x off-diagonal strain contribution)
+    """
+    def _prepare(self, domain: DomainDefinition, voigt: bool = True):
+        # Numerical integration
+        B = None
+        siz = domain.element_size
+        w = 1/domain.elemnodes  # Average strain at the integration points
+        for n in domain.node_numbering:
+            pos = n * (siz / 2) / np.sqrt(3)  # Sampling point
+            dN_dx = domain.eval_shape_fun_der(pos)
+            B_add = w*get_B(dN_dx)  # Add contribution
+            if B is None:
+                B = B_add
+            else:
+                B += B_add
+
+        if voigt:
+            idx_shear = np.count_nonzero(B, axis=1) == 2*domain.elemnodes  # Shear is combination of two displacements
+            B[idx_shear, :] *= 2  # Voigt notation with engineering strain
+
+        super()._prepare(domain, B)
+
+
+class Stress(Strain):
+    """ Calculate the average stresses per element
+
+    Input Signal:
+       - ``u``: Nodal vector of size ``(#dof per element * #nodes)``
+
+    Output Signal:
+       - ``s``: Stress matrix of size ``(#stresses per element, #elements)``
+
+    Args:
+        domain: The domain defining element and nodal connectivity
+
+    Keyword Args:
+        e_modulus: Young's modulus
+        poisson_ratio: Poisson ratio
+        plane: Plane 'strain' or 'stress'
+    """
+    def _prepare(self, domain: DomainDefinition,
+                 e_modulus: float = 1.0, poisson_ratio: float = 0.3, plane: str = 'strain'):
+
+        super()._prepare(domain, voigt=True)
+
+        # Get material relation
+        D = get_D(e_modulus, poisson_ratio, '3d' if domain.dim == 3 else plane.lower())
+        if domain.dim == 2:
+            D *= domain.element_size[2]
+        self.element_matrix = D @ self.element_matrix
diff --git a/pymoto/modules/complex.py b/pymoto/modules/complex.py
index 4224daf..abfbfef 100644
--- a/pymoto/modules/complex.py
+++ b/pymoto/modules/complex.py
@@ -73,7 +73,7 @@ def _sensitivity(self, dz):
 
 
 class RealPart(Module):
-    """ Takes the real part of a complex value :math:`x = \\text{Re}(z)`
+    r""" Takes the real part of a complex value :math:`x = \text{Re}(z)`
 
     Input Signal:
         - ``z``: Complex value
@@ -89,7 +89,7 @@ def _sensitivity(self, dx):
 
 
 class ImagPart(Module):
-    """ Takes the imaginary part of a complex value :math:`y = \\text{Im}(z)`
+    r""" Takes the imaginary part of a complex value :math:`y = \text{Im}(z)`
 
     Input Signal:
         - ``z``: Complex value
@@ -105,7 +105,7 @@ def _sensitivity(self, dy):
 
 
 class ComplexNorm(Module):
-    """ Takes the complex norm :math:`A = z z^*`
+    r""" Takes the complex norm :math:`A = \sqrt(z z^*)`
 
     Input Signal:
         - ``z``: Complex value
diff --git a/pymoto/modules/filter.py b/pymoto/modules/filter.py
index b3a7a44..fc54f03 100644
--- a/pymoto/modules/filter.py
+++ b/pymoto/modules/filter.py
@@ -1,46 +1,247 @@
-from pymoto import Module
-from .assembly import DomainDefinition
+from pymoto import Module, DomainDefinition
 import numpy as np
 from scipy.sparse import coo_matrix
-from scipy.ndimage import convolve
+from scipy.signal import convolve, correlate
+from numbers import Number
 
 
 class FilterConv(Module):
-    def _prepare(self, domain: DomainDefinition, weights: np.ndarray, mode: str = 'reflect'):
+    def _prepare(self, domain: DomainDefinition, radius: float = None, relative_units: bool = True, weights: np.ndarray = None,
+                 xmin_bc='symmetric', xmax_bc='symmetric',
+                 ymin_bc='symmetric', ymax_bc='symmetric',
+                 zmin_bc='symmetric', zmax_bc='symmetric'):
         self.domain = domain
-        self.weights = weights
-        self.mode = mode
-
-    def set_filter_radius(self, radius: float, relative_units: bool = False):
+        self.weights = None
+        if (weights is None and radius is None) or (weights is not None and radius is not None):
+            raise ValueError("Only one of arguments 'filter_radius' or 'weights' must be provided.")
+        elif weights is not None:
+            self.weights = weights.copy()
+            while self.weights.ndim < 3:
+                self.weights = np.expand_dims(self.weights, axis=-1)
+            for i in range(self.weights.ndim):
+                assert self.weights.shape[i] % 2 == 1, "Size of weights must be uneven"
+        elif radius is not None:
+            self.set_filter_radius(radius, relative_units)
+
+        # Process padding
+        self.overrides = []
+        self.pad_sizes = [v//2 for v in self.weights.shape]
+
+        domain_sizes = [self.domain.nelx, self.domain.nely, self.domain.nelz]
+        el_x, el_y, el_z = np.meshgrid(*[np.arange(max(1, s)) for s in domain_sizes], indexing='ij')
+        self.el3d_orig = self.domain.get_elemnumber(el_x, el_y, el_z)
+
+        padx = self._process_padding(self.el3d_orig, xmin_bc, xmax_bc, 0, self.pad_sizes[0])
+        pady = self._process_padding(padx, ymin_bc, ymax_bc, 1, self.pad_sizes[1])
+        self.el3d_pad = self._process_padding(pady, zmin_bc, zmax_bc, 2, self.pad_sizes[2])
+        # # Process padding of the domain
+        # if left == 'wrap' and right == 'wrap':
+        #     pad1a = np.pad(self.el3d_orig, [(self.pad_sizes[0], self.pad_sizes[0]), (0, 0), (0, 0)], mode='wrap')
+        # elif left == 'wrap':
+        #     pad1a = np.pad(self.el3d_orig, [(self.pad_sizes[0], 0), (0, 0), (0, 0)], mode='wrap')
+        # elif right == 'wrap':
+        #     pad1a = np.pad(self.el3d_orig, [(0, self.pad_sizes[0]), (0, 0), (0, 0)], mode='wrap')
+        # else:
+        #     pad1a = self.el3d_orig
+        #
+        # if right == 'edge':
+        #     pad1b = np.pad(pad1a, [(0, self.pad_sizes[0]), (0, 0), (0, 0)], mode='edge')
+        # elif right == 'symmetric':
+        #     pad1b = np.pad(pad1a, [(0, self.pad_sizes[0]), (0, 0), (0, 0)], mode='symmetric')
+        # elif isinstance(right, Number):
+        #     pad1b = np.pad(pad1a, [(0, self.pad_sizes[0]), (0, 0), (0, 0)], mode='constant', constant_values=0)
+        #     xrange = self.pad_sizes[0] + domain.nelx + np.arange(max(1, self.pad_sizes[0]))
+        #     yrange = np.arange(max(1, domain.nely + 2*self.pad_sizes[1]))
+        #     zrange = np.arange(max(1, domain.nelz + 2*self.pad_sizes[2]))
+        #     el_x, el_y, el_z = np.meshgrid(xrange, yrange, zrange, indexing='ij')
+        #     self.override_padded_values((el_x, el_y, el_z), right)
+        # else:
+        #     pad1b = pad1a
+        #
+        # if left == 'edge':
+        #     pad1 = np.pad(pad1b, [(self.pad_sizes[0], 0), (0, 0), (0, 0)], mode='edge')
+        # elif left == 'symmetric':
+        #     pad1 = np.pad(pad1b, [(self.pad_sizes[0], 0), (0, 0), (0, 0)], mode='symmetric')
+        # elif isinstance(left, Number):
+        #     pad1 = np.pad(pad1b, [(self.pad_sizes[0], 0), (0, 0), (0, 0)], mode='constant', constant_values=0)
+        #     xrange = np.arange(max(1, self.pad_sizes[0]))
+        #     yrange = np.arange(max(1, domain.nely + 2 * self.pad_sizes[1]))
+        #     zrange = np.arange(max(1, domain.nelz + 2 * self.pad_sizes[2]))
+        #     el_x, el_y, el_z = np.meshgrid(xrange, yrange, zrange, indexing='ij')
+        #     self.override_padded_values((el_x, el_y, el_z), left)
+        # else:
+        #     pad1 = pad1b
+        #
+        # # Y-direction
+        # if bottom == 'wrap' and top == 'wrap':
+        #     pad2a = np.pad(pad1, [(0, 0), (self.pad_sizes[1], self.pad_sizes[1]), (0, 0)], mode='wrap')
+        # elif bottom == 'wrap':
+        #     pad2a = np.pad(pad1, [(0, 0), (self.pad_sizes[1], 0), (0, 0)], mode='wrap')
+        # elif top == 'wrap':
+        #     pad2a = np.pad(pad1, [(0, 0), (0, self.pad_sizes[1]), (0, 0)], mode='wrap')
+        # else:
+        #     pad2a = pad1
+        #
+        # if top == 'edge':
+        #     pad2b = np.pad(pad2a, [(0, 0), (0, self.pad_sizes[1]), (0, 0)], mode='edge')
+        # elif top == 'symmetric':
+        #     pad2b = np.pad(pad2a, [(0, 0), (0, self.pad_sizes[1]), (0, 0)], mode='symmetric')
+        # elif isinstance(top, Number):
+        #     pad2b = np.pad(pad2a, [(0, 0), (0, self.pad_sizes[1]), (0, 0)], mode='constant', constant_values=0)
+        #     xrange = np.arange(max(1, domain.nelx + 2*self.pad_sizes[0]))
+        #     yrange = self.pad_sizes[1] + domain.nely + np.arange(max(1, self.pad_sizes[1]))
+        #     zrange = np.arange(max(1, domain.nelz + 2*self.pad_sizes[2]))
+        #     el_x, el_y, el_z = np.meshgrid(xrange, yrange, zrange, indexing='ij')
+        #     self.override_padded_values((el_x, el_y, el_z), top)
+        # else:
+        #     pad2b = pad2a
+        #
+        # if bottom == 'edge':
+        #     pad2 = np.pad(pad2b, [(0, 0), (self.pad_sizes[1], 0), (0, 0)], mode='edge')
+        # elif bottom == 'symmetric':
+        #     pad2 = np.pad(pad2b, [(0, 0), (self.pad_sizes[1], 0), (0, 0)], mode='symmetric')
+        # elif isinstance(bottom, Number):
+        #     pad2 = np.pad(pad2b, [(0, 0), (self.pad_sizes[1], 0), (0, 0)], mode='constant', constant_values=0)
+        #     xrange = np.arange(max(1, domain.nelx + 2 * self.pad_sizes[0]))
+        #     yrange = np.arange(max(1, self.pad_sizes[1]))
+        #     zrange = np.arange(max(1, domain.nelz + 2 * self.pad_sizes[2]))
+        #     el_x, el_y, el_z = np.meshgrid(xrange, yrange, zrange, indexing='ij')
+        #     self.override_padded_values((el_x, el_y, el_z), bottom)
+        # else:
+        #     pad2 = pad2b
+
+
+        # if isinstance(left, Number):
+            # mode = 'constant'
+            # mode = 'edge'
+            # mode = 'symmetric'
+            # mode = 'wrap'
+
+
+        # Do wrap first
+
+        # np.pad(..., **bottom)
+        # self.el3d_pad = pad2  #np.pad(self.el3d_orig, [(d, d) for d in pad_sizes], mode='reflect')
+
+    def _process_padding(self, indices, type_edge0, type_edge1, direction: int, pad_size: int):
+        # First process wrapped padding
+        wrap_size = (0, 0)
+        do_wrap = False
+        if type_edge0 == 'wrap':
+            do_wrap = True
+            wrap_size = (pad_size, wrap_size[1])
+        if type_edge1 == 'wrap':
+            do_wrap = True
+            wrap_size = (wrap_size[0], pad_size)
+
+        if do_wrap:
+            pad_width = [(0, 0) for _ in range(indices.ndim)]
+            pad_width[direction] = wrap_size
+            pad1a = np.pad(indices, pad_width, mode='wrap')
+        else:
+            pad1a = indices
+
+        domain_sizes = [self.domain.nelx, self.domain.nely, self.domain.nelz]
+        padded_sizes = [self.domain.nelx + 2 * self.pad_sizes[0],
+                        self.domain.nely + 2 * self.pad_sizes[1],
+                        self.domain.nelz + 2 * self.pad_sizes[2]]
+
+        # Process edge 1
+        pad_width = [(0, 0) for _ in range(indices.ndim)]
+        pad_width[direction] = (0, pad_size)
+        if type_edge1 == 'edge':
+            pad1b = np.pad(pad1a, pad_width, mode='edge')
+        elif type_edge1 == 'symmetric':
+            pad1b = np.pad(pad1a, pad_width, mode='symmetric')
+        elif isinstance(type_edge1, Number):  # Constant
+            value = type_edge1
+            pad1b = np.pad(pad1a, pad_width, mode='constant', constant_values=0)
+
+            n_range = [np.arange(max(1, s)) for s in padded_sizes]
+            n_range[direction] = pad_size + domain_sizes[direction] + np.arange(pad_size)
+
+            el_nos = np.meshgrid(*n_range, indexing='ij')
+            self.override_padded_values(tuple(el_nos), value)
+        else:
+            pad1b = pad1a
+
+        # Process edge 0
+        pad_width = [(0, 0) for _ in range(indices.ndim)]
+        pad_width[direction] = (pad_size, 0)
+        if type_edge0 == 'edge':
+            pad1 = np.pad(pad1b, pad_width, mode='edge')
+        elif type_edge0 == 'symmetric':
+            pad1 = np.pad(pad1b, pad_width, mode='symmetric')
+        elif isinstance(type_edge0, Number):
+            value = type_edge0
+            pad1 = np.pad(pad1b, pad_width, mode='constant', constant_values=0)
+
+            n_range = [np.arange(max(1, s)) for s in padded_sizes]
+            n_range[direction] = np.arange(pad_size)
+
+            el_nos = np.meshgrid(*n_range, indexing='ij')
+            self.override_padded_values(tuple(el_nos), value)
+        else:
+            pad1 = pad1b
+        return pad1
+
+    @property
+    def padded_domain(self):
+        domain_sizes = [self.domain.nelx, self.domain.nely, self.domain.nelz]
+        nx, ny, nz = [n + 2*p for n, p in zip(domain_sizes, self.pad_sizes)]
+        lx, ly, lz = self.domain.element_size
+        return DomainDefinition(nx, ny, nz, unitx=lx, unity=ly, unitz=lz)
+
+    def override_padded_values(self, index, value):
+        if all([np.asarray(i).size == 0 for i in index]):
+            # Don't add empty sets
+            return
+        self.overrides.append((index, value))
+
+    def override_values(self, index, value):
+        # Change index to extended domain
+        xrange = self.pad_sizes[0] + np.arange(max(1, self.domain.nelx))
+        yrange = self.pad_sizes[1] + np.arange(max(1, self.domain.nely))
+        zrange = self.pad_sizes[2] + np.arange(max(1, self.domain.nelz))
+        el_x, el_y, el_z = np.meshgrid(xrange, yrange, zrange, indexing='ij')
+        self.overrides.append(((el_x[index], el_y[index], el_z[index]), value))
+
+    def get_padded_vector(self, x):
+        xpad = x[self.el3d_pad]
+        for index, value in self.overrides:
+            xpad[index] = value
+        return xpad
+
+    def set_filter_radius(self, radius: float, relative_units: bool = True):
         if relative_units:
             dx, dy, dz = 1.0, 1.0, 1.0
         else:
             dx, dy, dz = self.domain.element_size
-        delemx, delemy, delemz = int((radius-1e-10*dx)/dx), int((radius-1e-10*dy)/dy), int((radius-1e-10*dz)/dz)
-        if self.domain.dim < 3:
-            delemz = 0
+        nx, ny, nz = self.domain.nelx, self.domain.nely, self.domain.nelz
+        delemx = min(nx, int((radius-1e-10*dx)/dx))
+        delemy = min(ny, int((radius-1e-10*dy)/dy))
+        delemz = min(nz, int((radius-1e-10*dz)/dz))
         xrange = np.arange(-delemx, delemx+1)*dx
         yrange = np.arange(-delemy, delemy+1)*dy
         zrange = np.arange(-delemz, delemz+1)*dz
-        coords_x, coords_y, coords_z = np.meshgrid(xrange, yrange, zrange)
+        coords_x, coords_y, coords_z = np.meshgrid(xrange, yrange, zrange, indexing='ij')
         self.weights = np.maximum(0.0, radius - np.sqrt(coords_x*coords_x + coords_y*coords_y + coords_z*coords_z))
         self.weights /= np.sum(self.weights)  # Volume preserving
-        if self.domain.dim < 3:
-            self.weights = self.weights[:, :, 0]
 
     def _response(self, x):
-        domain_sizes = [self.domain.nelx, self.domain.nely]
-        if self.domain.dim >= 3:
-            domain_sizes.append(self.domain.nely)
-        xbox = x.reshape(*domain_sizes, order='F').T  # TODO 3d?
-        return convolve(xbox, self.weights, mode=self.mode).flatten()
+        xpad = self.get_padded_vector(x)
+        y3d = convolve(xpad, self.weights, mode='valid')
+        y = np.zeros_like(x)
+        np.add.at(y, self.el3d_orig, y3d)
+        return y
 
     def _sensitivity(self, dfdv):
-        domain_sizes = [self.domain.nelx, self.domain.nely]
-        if self.domain.dim == 3:
-            domain_sizes.append(self.domain.nely)
-        ybox = dfdv.reshape(*domain_sizes, order='F').T  # TODO 3d
-        return convolve(ybox, self.weights, mode=self.mode).flatten()
+        dx3d = correlate(dfdv[self.el3d_orig], self.weights, mode='full')
+        for index, _ in self.overrides:
+            dx3d[index] = 0
+        dx = np.zeros_like(self.sig_in[0].state)
+        np.add.at(dx, self.el3d_pad, dx3d)
+        return dx
 
 
 class Filter(Module):
@@ -167,6 +368,8 @@ def _calculate_h(domain: DomainDefinition, radius=2.0):
         nwindy = yupp - ylow + 1
         nwindz = zupp - zlow + 1
         nwind = nwindx * nwindy * nwindz
+        if np.sum(nwind) < 0:
+            raise OverflowError("Filter size too large for this mesh size")
 
         # Total number of window elements
         ncum = np.cumsum(nwind)
diff --git a/pymoto/modules/generic.py b/pymoto/modules/generic.py
index 5cb0426..1fe4428 100644
--- a/pymoto/modules/generic.py
+++ b/pymoto/modules/generic.py
@@ -105,10 +105,21 @@ def _sensitivity(self, df_dy):
             if np.isrealobj(dg_dx[i]) and np.iscomplexobj(dg_dx_add):
                 dg_dx_add = np.real(dg_dx_add)
 
-            # Add the contribution
+            # Add the contribution according to broadcasting rules of NumPy
+            # https://numpy.org/doc/stable/user/basics.broadcasting.html
             if (not hasattr(dg_dx[i], '__len__')) or (hasattr(dg_dx[i], 'ndim') and dg_dx[i].ndim == 0):
                 # Scalar type or 0-dimensional array
                 dg_dx[i] += np.sum(dg_dx_add)
+            elif dg_dx[i].shape != dg_dx_add.shape:
+                # Reverse broadcast https://stackoverflow.com/questions/76002989/numpy-is-there-a-reverse-broadcast
+                n_leading_dims = dg_dx_add.ndim - dg_dx[i].ndim
+                broadcasted_dims = tuple(range(n_leading_dims))
+                for ii in range(dg_dx_add.ndim - n_leading_dims):
+                    if dg_dx[i].shape[ii] == 1 and dg_dx_add.shape[ii+n_leading_dims] != 1:
+                        broadcasted_dims = (*broadcasted_dims, n_leading_dims+ii)
+
+                dg_dx_add1 = np.add.reduce(dg_dx_add, axis=broadcasted_dims, keepdims=True)  # Sum broadcasted axis
+                dg_dx[i] += np.squeeze(dg_dx_add1, axis=tuple(range(n_leading_dims)))  # Squeeze out singleton axis
             else:
                 dg_dx[i] += dg_dx_add
 
diff --git a/pymoto/modules/io.py b/pymoto/modules/io.py
index 4608d49..760275b 100644
--- a/pymoto/modules/io.py
+++ b/pymoto/modules/io.py
@@ -14,7 +14,7 @@
 from .assembly import DomainDefinition
 
 
-class _FigModule(Module):
+class FigModule(Module):
     """ Abstract base class for any module which produces a figure
 
     Keyword Args:
@@ -55,7 +55,7 @@ def __del__(self):
         plt.close(self.fig)
 
 
-class PlotDomain(_FigModule):
+class PlotDomain(FigModule):
     """ Plots the densities of a domain (2D or 3D)
 
     Input Signal:
@@ -96,10 +96,16 @@ def _plot_2d(self, x):
             self.im.set_data(data)
         else:
             ax = self.fig.add_subplot(111)
-            self.im = ax.imshow(data, origin='lower', cmap=self.cmap)
+            Lx = self.domain.nelx * self.domain.unitx
+            Ly = self.domain.nely * self.domain.unity
+            self.im = ax.imshow(data, cmap=self.cmap, origin='lower', extent=(0.0, Lx, 0.0, Ly))
             self.cbar = self.fig.colorbar(self.im, orientation='horizontal')
             ax.set(xlabel='x', ylabel='y')
-        clim = [np.min(data), np.max(data)] if self.clim is None else self.clim
+        vmin, vmax = np.min(data), np.max(data)
+        if vmin < 0:
+            vabs = max(abs(vmin), abs(vmax))
+            vmin, vmax = -vabs, vabs
+        clim = [vmin, vmax] if self.clim is None else self.clim
         self.im.set_clim(vmin=clim[0], vmax=clim[1])
 
     def _plot_3d(self, x):
@@ -138,7 +144,7 @@ def _plot_3d(self, x):
         self.fac = ax.voxels(sel, facecolors=colors, linewidth=0.5, edgecolors='k')
 
 
-class PlotGraph(_FigModule):
+class PlotGraph(FigModule):
     """ Plot an X-Y graph
 
     Input Signals:
@@ -175,14 +181,15 @@ def _response(self, x, *ys):
         for i, y in enumerate(ys):
             self.line[i].set_xdata(x)
             self.line[i].set_ydata(y)
-            ymin, ymax = min(ymin, min(y)), max(ymax, max(y))
-        self.ax.set_xlim([min(x), max(x)])
-        self.ax.set_ylim([ymin, ymax])
+            ymin, ymax = min(ymin, np.min(y)), max(ymax, np.max(y))
+        self.ax.set_xlim([np.min(x), np.max(x)])
+        dy = ymax - ymin
+        self.ax.set_ylim([ymin-0.05*dy, ymax+0.05*dy])
 
         self._update_fig()
 
 
-class PlotIter(_FigModule):
+class PlotIter(FigModule):
     """ Plot iteration history of one or more variables
 
     Input Signals:
@@ -193,10 +200,12 @@ class PlotIter(_FigModule):
         overwrite (bool): Overwrite saved image every time the figure is updated, else prefix ``_0000`` is added to the
           filename (default = ``False``)
         show (bool): Show the figure on the screen
+        ylim: Provide y-axis limits for the plot
     """
-    def _prepare(self):
+    def _prepare(self, ylim=None):
         self.minlim = 1e+200
         self.maxlim = -1e+200
+        self.ylim = ylim
 
     def _response(self, *args):
         if not hasattr(self, 'ax'):
@@ -227,7 +236,9 @@ def _response(self, *args):
         dy = max((self.maxlim - self.minlim)*0.05, sys.float_info.min)
 
         self.ax.set_xlim([-0.5, self.iter+0.5])
-        if np.isfinite(self.minlim) and np.isfinite(self.maxlim):
+        if self.ylim is not None:
+            self.ax.set_ylim(self.ylim)
+        elif np.isfinite(self.minlim) and np.isfinite(self.maxlim):
             self.ax.set_ylim([self.minlim - dy, self.maxlim + dy])
 
         self._update_fig()
diff --git a/pymoto/modules/linalg.py b/pymoto/modules/linalg.py
index 968dd2f..56a149e 100644
--- a/pymoto/modules/linalg.py
+++ b/pymoto/modules/linalg.py
@@ -7,10 +7,9 @@
 import scipy.sparse as sps
 import scipy.sparse.linalg as spsla
 
-from pymoto import Signal, Module, DyadCarrier, LDAWrapper
-from pymoto import SolverDenseLU, SolverDenseLDL, SolverDenseCholesky, SolverDiagonal, SolverDenseQR
-from pymoto import SolverSparseLU, SolverSparseCholeskyCVXOPT, SolverSparsePardiso, SolverSparseCholeskyScikit
-from pymoto import matrix_is_symmetric, matrix_is_hermitian, matrix_is_diagonal
+from pymoto import Signal, Module, DyadCarrier
+from pymoto.solvers import auto_determine_solver
+from pymoto.solvers import matrix_is_hermitian, LDAWrapper
 
 
 class StaticCondensation(Module):
@@ -172,7 +171,7 @@ def _sensitivity(self, dgdx, dgdb):
             adjoint_load += self.Afp * dgdb[self.p, ...]
 
         lam = np.zeros_like(self.x)
-        lamf = -1.0 * self.module_LinSolve.solver.adjoint(adjoint_load)
+        lamf = -1.0 * self.module_LinSolve.solver.solve(adjoint_load, trans='T')
         lam[self.f, ...] = lamf
 
         if dgdb is not None:
@@ -219,105 +218,6 @@ def _sensitivity(self, dB):
         return dA if np.iscomplexobj(A) else np.real(dA)
 
 
-# flake8: noqa: C901
-def auto_determine_solver(A, isdiagonal=None, islowertriangular=None, isuppertriangular=None,
-                          ishermitian=None, issymmetric=None, ispositivedefinite=None):
-    """
-    Uses parts of Matlab's scheme https://nl.mathworks.com/help/matlab/ref/mldivide.html
-    :param A: The matrix
-    :param isdiagonal: Manual override for diagonal matrix
-    :param islowertriangular: Override for lower triangular matrix
-    :param isuppertriangular: Override for upper triangular matrix
-    :param ishermitian: Override for hermitian matrix (prevents check)
-    :param issymmetric: Override for symmetric matrix (prevents check). Is the same as hermitian for a real matrix
-    :param ispositivedefinite: Manual override for positive definiteness
-    :return: LinearSolver which should be 'best' for the matrix
-    """
-    issparse = sps.issparse(A)  # Check if the matrix is sparse
-    issquare = A.shape[0] == A.shape[1]  # Check if the matrix is square
-
-    if not issquare:
-        if issparse:
-            sps.SparseEfficiencyWarning("Only a dense version of QR solver is available")  # TODO
-        return SolverDenseQR()
-
-    # l_bw, u_bw = spla.bandwidth(A) # TODO Get bandwidth (implemented in scipy version > 1.8.0)
-
-    if isdiagonal is None:  # Check if matrix is diagonal
-        # TODO: This could be improved to check other sparse matrix types as well
-        isdiagonal = matrix_is_diagonal(A)
-    if isdiagonal:
-        return SolverDiagonal()
-
-    # Check if the matrix is triangular
-    # TODO Currently only for dense matrices
-    if islowertriangular is None:  # Check if matrix is lower triangular
-        islowertriangular = False if issparse else np.allclose(A, np.tril(A))
-    if islowertriangular:
-        warnings.WarningMessage("Lower triangular solver not implemented",
-                                UserWarning, getframeinfo(currentframe()).filename, getframeinfo(currentframe()).lineno)
-
-    if isuppertriangular is None:  # Check if matrix is upper triangular
-        isuppertriangular = False if issparse else np.allclose(A, np.triu(A))
-    if isuppertriangular:
-        warnings.WarningMessage("Upper triangular solver not implemented",
-                                UserWarning, getframeinfo(currentframe()).filename, getframeinfo(currentframe()).lineno)
-
-    ispermutedtriangular = False
-    if ispermutedtriangular:
-        warnings.WarningMessage("Permuted triangular solver not implemented",
-                                UserWarning, getframeinfo(currentframe()).filename, getframeinfo(currentframe()).lineno)
-
-    # Check if the matrix is complex-valued
-    iscomplex = np.iscomplexobj(A)
-    if iscomplex:
-        # Detect if the matrix is hermitian and/or symmetric
-        if ishermitian is None:
-            ishermitian = matrix_is_hermitian(A)
-        if issymmetric is None:
-            issymmetric = matrix_is_symmetric(A)
-    else:
-        if ishermitian is None and issymmetric is None:
-            # Detect if the matrix is symmetric
-            issymmetric = matrix_is_symmetric(A)
-            ishermitian = issymmetric
-        elif ishermitian is not None and issymmetric is not None:
-            assert ishermitian == issymmetric, "For real-valued matrices, symmetry and hermitian must be equal"
-        elif ishermitian is None:
-            ishermitian = issymmetric
-        elif issymmetric is None:
-            issymmetric = ishermitian
-
-    if issparse:
-        # Prefer Intel Pardiso solver as it can solve any matrix TODO: Check for complex matrix
-        if SolverSparsePardiso.defined and not iscomplex:
-            # TODO check for positive definiteness?  np.all(A.diagonal() > 0) or np.all(A.diagonal() < 0)
-            return SolverSparsePardiso(symmetric=issymmetric, hermitian=ishermitian, positive_definite=ispositivedefinite)
-
-        if ishermitian:
-            # Check if diagonal is all positive or all negative -> Cholesky
-            if np.all(A.diagonal() > 0) or np.all(A.diagonal() < 0):  # TODO what about the complex case?
-                if SolverSparseCholeskyScikit.defined:
-                    return SolverSparseCholeskyScikit()
-                if SolverSparseCholeskyCVXOPT.defined:
-                    return SolverSparseCholeskyCVXOPT()
-
-        return SolverSparseLU()  # Default to LU, which should be possible for any non-singular square matrix
-
-    else:  # Dense branch
-        if ishermitian:
-            # Check if diagonal is all positive or all negative
-            if np.all(A.diagonal() > 0) or np.all(A.diagonal() < 0):
-                return SolverDenseCholesky()
-            else:
-                return SolverDenseLDL(hermitian=ishermitian)
-        elif issymmetric:
-            return SolverDenseLDL(hermitian=ishermitian)
-        else:
-            # TODO: Detect if the matrix is Hessenberg
-            return SolverDenseLU()
-
-
 class LinSolve(Module):
     r""" Solves linear system of equations :math:`\mathbf{A}\mathbf{x}=\mathbf{b}`
 
@@ -347,6 +247,7 @@ def _prepare(self, dep_tol=1e-5, hermitian=None, symmetric=None, solver=None):
         self.ishermitian = hermitian
         self.issymmetric = symmetric
         self.solver = solver
+        self.u = None  # Solution storage
 
     def _response(self, mat, rhs):
         # Do some detections on the matrix type
@@ -365,19 +266,23 @@ def _response(self, mat, rhs):
         if self.solver is None:
             self.solver = auto_determine_solver(mat, ishermitian=self.ishermitian)
         if not isinstance(self.solver, LDAWrapper) and self.use_lda_solver:
-            self.solver = LDAWrapper(self.solver, hermitian=self.ishermitian, symmetric=self.issymmetric)
+            lda_kwargs = dict(hermitian=self.ishermitian, symmetric=self.issymmetric)
+            if hasattr(self.solver, 'tol'):
+                lda_kwargs['tol'] = self.solver.tol * 2
+            self.solver = LDAWrapper(self.solver, **lda_kwargs)
 
         # Update solver with new matrix
         self.solver.update(mat)
 
         # Solution
-        self.u = self.solver.solve(rhs)
+        self.u = self.solver.solve(rhs, x0=self.u)
 
         return self.u
 
     def _sensitivity(self, dfdv):
         mat, rhs = [s.state for s in self.sig_in]
-        lam = self.solver.adjoint(dfdv.conj()).conj()
+        # lam = self.solver.solve(dfdv.conj(), trans='H').conj()
+        lam = self.solver.solve(dfdv, trans='T')
 
         if self.issparse:
             if self.u.ndim > 1:
@@ -440,6 +345,7 @@ def _prepare(self, sorting_func=lambda W, Q: np.argsort(W), hermitian=None, nmod
         self.sigma = sigma
         self.mode = mode
         self.Ainv = None
+        self.do_solve = False
 
     def _response(self, A, *args):
         B = args[0] if len(args) > 0 else None
@@ -500,9 +406,14 @@ def _sparse_eigs(self, A, B=None):
         # Use shift-and-invert, so make inverse operator
         if self.Ainv is None:
             self.Ainv = auto_determine_solver(mat_shifted, ishermitian=self.is_hermitian)
-        self.Ainv.update(mat_shifted)
-
-        AinvOp = spsla.LinearOperator(mat_shifted.shape, matvec=self.Ainv.solve, rmatvec=self.Ainv.adjoint)
+            self.do_solve = True
+        if self.sigma != 0:
+            self.do_solve = True
+        if self.do_solve:
+            self.Ainv.update(mat_shifted)
+
+        AinvOp = spsla.LinearOperator(mat_shifted.shape, matvec=self.Ainv.solve,
+                                      rmatvec=lambda b: self.Ainv.solve(b, trans='H'))
 
         if self.is_hermitian:
             return spsla.eigsh(A, M=B, k=self.nmodes, OPinv=AinvOp, sigma=self.sigma, mode=self.mode)
diff --git a/pymoto/modules/scaling.py b/pymoto/modules/scaling.py
index 1a48255..74d644b 100644
--- a/pymoto/modules/scaling.py
+++ b/pymoto/modules/scaling.py
@@ -1,12 +1,12 @@
 from pymoto import Module
-
+import numpy as np
 
 class Scaling(Module):
     r""" Scales (scalar) input for different response functions in optimization (objective / constraints).
     This is useful, for instance, for MMA where the objective must be scaled in a certain way for good convergence.
 
-    Objective scaling (`minval` and `maxval` are both undefined):
-    :math:`y^{(i)} = s \frac{x^{(i)}}{x^{(0)}}`
+    Objective scaling using absolute value or vector norm (`minval` and `maxval` are both undefined):
+    :math:`y^{(i)} = s \frac{x^{(i)}}{||x^{(0)}||}`
 
     For the constraints, the negative null-form convention is used, which means the constraint is :math:`y(x) \leq 0`.
 
@@ -39,7 +39,7 @@ def _prepare(self, scaling: float = 100.0, minval: float = None, maxval: float =
 
     def _response(self, x):
         if not hasattr(self, 'sf'):
-            self.sf = self.scaling/x
+            self.sf = self.scaling/np.linalg.norm(x)
         if self.minval is not None:
             g = 1 - x/self.minval
         elif self.maxval is not None:
diff --git a/pymoto/routines.py b/pymoto/routines.py
index a48a684..4ed553a 100644
--- a/pymoto/routines.py
+++ b/pymoto/routines.py
@@ -1,8 +1,10 @@
+import warnings
 import numpy as np
 from .utils import _parse_to_list, _concatenate_to_array
 from .core_objects import Signal, SignalSlice, Module, Network
 from .common.mma import MMA
 from typing import List, Iterable, Union, Callable
+from scipy.sparse import issparse
 
 
 def _has_signal_overlap(sig1: List[Signal], sig2: List[Signal]):
@@ -171,7 +173,10 @@ def finite_difference(blk: Module, fromsig: Union[Signal, Iterable[Signal]] = No
                 fp = Sout.state
 
                 # Finite difference sensitivity
-                df = (fp - f0[Iout])/(dx*sf)
+                if issparse(fp):
+                    df = (fp.toarray() - f0[Iout].toarray()) / (dx * sf)
+                else:
+                    df = (fp - f0[Iout])/(dx*sf)
 
                 dgdx_fd = np.real(np.sum(df*df_an[Iout]))
 
@@ -223,7 +228,10 @@ def finite_difference(blk: Module, fromsig: Union[Signal, Iterable[Signal]] = No
                     fp = Sout.state
 
                     # Finite difference sensitivity
-                    df = (fp - f0[Iout])/(dx*1j*sf)
+                    if issparse(fp):
+                        df = (fp.toarray() - f0[Iout].toarray()) / (dx * 1j * sf)
+                    else:
+                        df = (fp - f0[Iout])/(dx*1j*sf)
                     dgdx_fd = np.imag(np.sum(df*df_an[Iout]))
 
                     if dx_an[Iout][Iin] is not None:
@@ -335,8 +343,9 @@ def minimize_oc(function, variables, objective: Signal,
         function.sensitivity()
         dfdx, _ = _concatenate_to_array(obtain_sensitivities(variables))
         maxdfdx = max(dfdx)
-        if maxdfdx > 0:
-            raise RuntimeError(f"OC only works for negative sensitivities: max(dfdx) = {maxdfdx}")
+        if maxdfdx > 1e-15:
+            warnings.warn(f"OC only works for negative sensitivities: max(dfdx) = {maxdfdx}. Clipping positive values.")
+        dfdx = np.minimum(dfdx, 0)
 
         # Do OC update
         l1, l2 = l1init, l2init
diff --git a/pymoto/solvers/__init__.py b/pymoto/solvers/__init__.py
new file mode 100644
index 0000000..636a423
--- /dev/null
+++ b/pymoto/solvers/__init__.py
@@ -0,0 +1,14 @@
+from .solvers import LinearSolver, LDAWrapper
+from .matrix_checks import matrix_is_complex, matrix_is_diagonal, matrix_is_symmetric, matrix_is_hermitian
+from .dense import SolverDiagonal, SolverDenseQR, SolverDenseLU, SolverDenseCholesky, SolverDenseLDL
+from .sparse import SolverSparsePardiso, SolverSparseLU, SolverSparseCholeskyScikit, SolverSparseCholeskyCVXOPT
+from .iterative import Preconditioner, CG, DampedJacobi, SOR, ILU, GeometricMultigrid
+from .auto_determine import auto_determine_solver
+
+__all__ = ['matrix_is_complex', 'matrix_is_diagonal', 'matrix_is_symmetric', 'matrix_is_hermitian',
+           'LinearSolver', 'LDAWrapper',
+           'SolverDiagonal', 'SolverDenseQR', 'SolverDenseLU', 'SolverDenseCholesky', 'SolverDenseLDL',
+           'SolverSparsePardiso', 'SolverSparseLU', 'SolverSparseCholeskyScikit', 'SolverSparseCholeskyCVXOPT',
+           'Preconditioner', 'CG', 'DampedJacobi', 'SOR', 'ILU', 'GeometricMultigrid',
+           'auto_determine_solver',
+           ]
\ No newline at end of file
diff --git a/pymoto/solvers/auto_determine.py b/pymoto/solvers/auto_determine.py
new file mode 100644
index 0000000..7ad9709
--- /dev/null
+++ b/pymoto/solvers/auto_determine.py
@@ -0,0 +1,108 @@
+import warnings
+import scipy.sparse as sps
+from inspect import currentframe, getframeinfo
+
+from .dense import *
+from .sparse import *
+from .matrix_checks import *
+
+
+# flake8: noqa: C901
+def auto_determine_solver(A, isdiagonal=None, islowertriangular=None, isuppertriangular=None,
+                          ishermitian=None, issymmetric=None, ispositivedefinite=None):
+    """
+    Uses parts of Matlab's scheme https://nl.mathworks.com/help/matlab/ref/mldivide.html
+    :param A: The matrix
+    :param isdiagonal: Manual override for diagonal matrix
+    :param islowertriangular: Override for lower triangular matrix
+    :param isuppertriangular: Override for upper triangular matrix
+    :param ishermitian: Override for hermitian matrix (prevents check)
+    :param issymmetric: Override for symmetric matrix (prevents check). Is the same as hermitian for a real matrix
+    :param ispositivedefinite: Manual override for positive definiteness
+    :return: LinearSolver which should be 'best' for the matrix
+    """
+    issparse = sps.issparse(A)  # Check if the matrix is sparse
+    issquare = A.shape[0] == A.shape[1]  # Check if the matrix is square
+
+    if not issquare:
+        if issparse:
+            sps.SparseEfficiencyWarning("Only a dense version of QR solver is available")  # TODO
+        return SolverDenseQR()
+
+    # l_bw, u_bw = spla.bandwidth(A) # TODO Get bandwidth (implemented in scipy version > 1.8.0)
+
+    if isdiagonal is None:  # Check if matrix is diagonal
+        # TODO: This could be improved to check other sparse matrix types as well
+        isdiagonal = matrix_is_diagonal(A)
+    if isdiagonal:
+        return SolverDiagonal()
+
+    # Check if the matrix is triangular
+    # TODO Currently only for dense matrices
+    if islowertriangular is None:  # Check if matrix is lower triangular
+        islowertriangular = False if issparse else np.allclose(A, np.tril(A))
+    if islowertriangular:
+        warnings.WarningMessage("Lower triangular solver not implemented",
+                                UserWarning, getframeinfo(currentframe()).filename, getframeinfo(currentframe()).lineno)
+
+    if isuppertriangular is None:  # Check if matrix is upper triangular
+        isuppertriangular = False if issparse else np.allclose(A, np.triu(A))
+    if isuppertriangular:
+        warnings.WarningMessage("Upper triangular solver not implemented",
+                                UserWarning, getframeinfo(currentframe()).filename, getframeinfo(currentframe()).lineno)
+
+    ispermutedtriangular = False
+    if ispermutedtriangular:
+        warnings.WarningMessage("Permuted triangular solver not implemented",
+                                UserWarning, getframeinfo(currentframe()).filename, getframeinfo(currentframe()).lineno)
+
+    # Check if the matrix is complex-valued
+    iscomplex = np.iscomplexobj(A)
+    if iscomplex:
+        # Detect if the matrix is hermitian and/or symmetric
+        if ishermitian is None:
+            ishermitian = matrix_is_hermitian(A)
+        if issymmetric is None:
+            issymmetric = matrix_is_symmetric(A)
+    else:
+        if ishermitian is None and issymmetric is None:
+            # Detect if the matrix is symmetric
+            issymmetric = matrix_is_symmetric(A)
+            ishermitian = issymmetric
+        elif ishermitian is not None and issymmetric is not None:
+            assert ishermitian == issymmetric, "For real-valued matrices, symmetry and hermitian must be equal"
+        elif ishermitian is None:
+            ishermitian = issymmetric
+        elif issymmetric is None:
+            issymmetric = ishermitian
+
+    if issparse:
+        # Prefer Intel Pardiso solver as it can solve any matrix TODO: Check for complex matrix
+        if SolverSparsePardiso.defined and not iscomplex:
+            # TODO check for positive definiteness?  np.all(A.diagonal() > 0) or np.all(A.diagonal() < 0)
+            return SolverSparsePardiso(symmetric=issymmetric, hermitian=ishermitian, positive_definite=ispositivedefinite)
+
+        if ishermitian:
+            # Check if diagonal is all positive or all negative -> Cholesky
+            if ispositivedefinite is None:
+                ispositivedefinite = np.all(A.diagonal() > 0) or np.all(A.diagonal() < 0)
+            if ispositivedefinite:  # TODO what about the complex case?
+                if SolverSparseCholeskyScikit.defined:
+                    return SolverSparseCholeskyScikit()
+                if SolverSparseCholeskyCVXOPT.defined:
+                    return SolverSparseCholeskyCVXOPT()
+
+        return SolverSparseLU()  # Default to LU, which should be possible for any non-singular square matrix
+
+    else:  # Dense branch
+        if ishermitian:
+            # Check if diagonal is all positive or all negative
+            if np.all(A.diagonal() > 0) or np.all(A.diagonal() < 0):
+                return SolverDenseCholesky()
+            else:
+                return SolverDenseLDL(hermitian=ishermitian)
+        elif issymmetric:
+            return SolverDenseLDL(hermitian=ishermitian)
+        else:
+            # TODO: Detect if the matrix is Hessenberg
+            return SolverDenseLU()
diff --git a/pymoto/common/solvers_dense.py b/pymoto/solvers/dense.py
similarity index 54%
rename from pymoto/common/solvers_dense.py
rename to pymoto/solvers/dense.py
index 3fd8c35..201c989 100644
--- a/pymoto/common/solvers_dense.py
+++ b/pymoto/solvers/dense.py
@@ -1,7 +1,8 @@
 import warnings
 import numpy as np
 import scipy.linalg as spla  # Dense matrix solvers
-from .solvers import matrix_is_hermitian, matrix_is_diagonal, LinearSolver
+from .matrix_checks import matrix_is_hermitian, matrix_is_diagonal
+from .solvers import LinearSolver
 
 
 class SolverDiagonal(LinearSolver):
@@ -11,24 +12,17 @@ def update(self, A):
         self.diag = A.diagonal()
         return self
 
-    def solve(self, rhs):
+    def solve(self, rhs, x0=None, trans='N'):
         r""" Solve using the diagonal only, by :math:`x_i = b_i / A_{ii}`
 
         The right-hand-side :math:`\mathbf{b}` can be of size ``(N)`` or ``(N, K)``, where ``N`` is the size of matrix
         :math:`\mathbf{A}` and ``K`` is the number of right-hand sides.
         """
+        d = self.diag.conj() if trans == 'H' else self.diag
         if rhs.ndim == 1:
-            return rhs / self.diag
+            return rhs / d
         else:
-            return rhs / self.diag[..., None]
-
-    def adjoint(self, rhs):
-        r""" Solve using the diagonal only, by :math:`x_i = b_i / A_{ii}^*`
-
-        The right-hand-side :math:`\mathbf{b}` can be of size ``(N)`` or ``(N, K)``, where ``N`` is the size of matrix
-        :math:`\mathbf{A}` and ``K`` is the number of right-hand sides.
-        """
-        return self.solve(rhs.conj()).conj()
+            return rhs / d[..., None]
 
 
 # Dense QR solver
@@ -42,23 +36,31 @@ def update(self, A):
         self.q, self.r = spla.qr(A)
         return self
 
-    def solve(self, rhs):
-        r""" Solves the linear system of equations :math:`\mathbf{A} \mathbf{x} = \mathbf{b}` by backward substitution of
-        :math:`\mathbf{x} = \mathbf{R}^{-1}\mathbf{Q}^\text{H}\mathbf{b}`.
-
-        The right-hand-side :math:`\mathbf{b}` can be of size ``(N)`` or ``(N, K)``, where ``N`` is the size of matrix
-        :math:`\mathbf{A}` and ``K`` is the number of right-hand sides.
-        """
-        return spla.solve_triangular(self.r, self.q.T.conj()@rhs)
+    def solve(self, rhs, x0=None, trans='N'):
+        r""" Solves the linear system of equations using the QR factorization.
 
-    def adjoint(self, rhs):
-        r""" Solves the linear system of equations :math:`\mathbf{A}^\text{H}\mathbf{x} = \mathbf{b}` by
-        forward substitution of :math:`\mathbf{x} = \mathbf{Q}\mathbf{R}^{-H}\mathbf{b}`.
+        ======= ================= =====================
+        `trans`     Equation      Solution of :math:`x`
+        ------- ----------------- ---------------------
+          `N`   :math:`A x = b`   :math:`R^{-1} Q^H b`
+          `T`   :math:`A^T x = b` :math:`Q^* R^{-T} b`
+          `H`   :math:`A^H x = b` :math:`Q R^{-H} b`
+        ======= ================= =====================
 
         The right-hand-side :math:`\mathbf{b}` can be of size ``(N)`` or ``(N, K)``, where ``N`` is the size of matrix
         :math:`\mathbf{A}` and ``K`` is the number of right-hand sides.
         """
-        return self.q@spla.solve_triangular(self.r, rhs, trans='C')
+        if trans == 'N':
+            # A = Q R -> inv(A) = inv(R) inv(Q) = inv(R) Q^H
+            return spla.solve_triangular(self.r, self.q.T.conj() @ rhs)
+        elif trans == 'T':
+            # A^T = R^T Q^T -> inv(A^T) = inv(Q^T) inv(R^T) = conj(Q) inv(R^T)
+            return self.q.conj() @ spla.solve_triangular(self.r, rhs, trans='T')
+        elif trans == 'H':
+            # A^H = R^H Q^H -> inv(A^H) = inv(Q^H) inv(R^H) = Q inv(R^H)
+            return self.q @ spla.solve_triangular(self.r, rhs, trans='C')
+        else:
+            raise TypeError("Only N, T, and H transposition is possible")
 
 
 # Dense LU solver
@@ -71,24 +73,34 @@ def update(self, A):
         self.p, self.l, self.u = spla.lu(A)
         return self
 
-    def solve(self, rhs):
-        r""" Solves the linear system of equations :math:`\mathbf{A} \mathbf{x} = \mathbf{b}` by forward and backward
-        substitution of :math:`\mathbf{x} = \mathbf{U}^{-1}\mathbf{L}^{-1}\mathbf{b}`.
+    def solve(self, rhs, x0=None, trans='N'):
+        r""" Solves the linear system of equations using the LU factorization.
 
-        The right-hand-side :math:`\mathbf{b}` can be of size ``(N)`` or ``(N, K)``, where ``N`` is the size of matrix
-        :math:`\mathbf{A}` and ``K`` is the number of right-hand sides.
-        """
-        return spla.solve_triangular(self.u, spla.solve_triangular(self.l, self.p.T@rhs, lower=True))
+        :math:`\mathbf{A} \mathbf{x} = \mathbf{b}` by forward and backward
+        substitution of :math:`\mathbf{x} = \mathbf{U}^{-1}\mathbf{L}^{-1}\mathbf{b}`.
 
-    def adjoint(self, rhs):
-        r""" Solves the linear system of equations :math:`\mathbf{A}^\text{H}\mathbf{x} = \mathbf{b}` by forward and
-        backward substitution of :math:`\mathbf{x} = \mathbf{L}^{-\text{H}}\mathbf{U}^{-\text{H}}\mathbf{b}`.
+        ======= ================= =========================
+        `trans`     Equation        Solution of :math:`x`
+        ------- ----------------- -------------------------
+          `N`   :math:`A x = b`   :math:`x = U^{-1} L^{-1}`
+          `T`   :math:`A^T x = b` :math:`x = L^{-1} U^{-1}`
+          `H`   :math:`A^H x = b` :math:`x = L^{-*} U^{-*}`
+        ======= ================= =========================
 
         The right-hand-side :math:`\mathbf{b}` can be of size ``(N)`` or ``(N, K)``, where ``N`` is the size of matrix
         :math:`\mathbf{A}` and ``K`` is the number of right-hand sides.
         """
-        return self.p@spla.solve_triangular(self.l, spla.solve_triangular(self.u, rhs, trans='C'),
-                                            lower=True, trans='C')  # TODO permutation
+        if trans == 'N':
+            # A = P L U -> x = U^-1 L^-1 P^T b
+            return spla.solve_triangular(self.u, spla.solve_triangular(self.l, self.p.T@rhs, lower=True))
+        elif trans == 'T':
+            return self.p @ spla.solve_triangular(self.l, spla.solve_triangular(self.u, rhs, trans='T'),
+                                                  lower=True, trans='T')
+        elif trans == 'H':
+            return self.p @ spla.solve_triangular(self.l, spla.solve_triangular(self.u, rhs, trans='C'),
+                                                  lower=True, trans='C')
+        else:
+            raise TypeError("Only N, T, and H transposition is possible")
 
 
 # Dense Cholesky solver
@@ -106,7 +118,7 @@ def update(self, A):
         upper triangular matrix.
         """
         try:
-            self.u = spla.cholesky(A)
+            self.U = spla.cholesky(A)
             self.success = True
         except np.linalg.LinAlgError as err:
             warnings.warn(f"{type(self).__name__}: {err} -- using {type(self.backup_solver).__name__} instead")
@@ -114,7 +126,7 @@ def update(self, A):
             self.success = False
         return self
 
-    def solve(self, rhs):
+    def solve(self, rhs, x0=None, trans='N'):
         r""" Solves the linear system of equations :math:`\mathbf{A} \mathbf{x} = \mathbf{b}` by forward and backward
         substitution of :math:`\mathbf{x} = \mathbf{U}^{-1}\mathbf{U}^{-\text{H}}\mathbf{b}`.
 
@@ -124,21 +136,16 @@ def solve(self, rhs):
         # TODO When Cholesky factorization A = U^T U is used, symmetric complex matrices can also be solved, but this is
         #  not implemented in scipy
         if self.success:
-            return spla.solve_triangular(self.u, spla.solve_triangular(self.u, rhs, trans='C'))
-        else:
-            return self.backup_solver.solve(rhs)
-
-    def adjoint(self, rhs):
-        r""" A Hermitian matrix is self-adjoint (:math:`\mathbf{A}=\mathbf{A}^\text{H}`), so this is equal to the
-        regular solution.
-
-        The right-hand-side :math:`\mathbf{b}` can be of size ``(N)`` or ``(N, K)``, where ``N`` is the size of matrix
-        :math:`\mathbf{A}` and ``K`` is the number of right-hand sides.
-        """
-        if self.success:
-            return self.solve(rhs)
+            if trans == 'N' or trans == 'H':
+                # A = U^H U -> A^-1 = U^-1 U^-H
+                return spla.solve_triangular(self.U, spla.solve_triangular(self.U, rhs, trans='C'))
+            elif trans == 'T':
+                # A^T = U^T conj(U) -> A^-T = conj(U^-1) U^-T
+                return spla.solve_triangular(self.U, spla.solve_triangular(self.U, rhs, trans='T').conj()).conj()
+            else:
+                raise TypeError("Only N, T, and H transposition is possible")
         else:
-            return self.backup_solver.adjoint(rhs)
+            return self.backup_solver.solve(rhs, trans=trans)
 
 
 # Dense LDL solver
@@ -171,23 +178,13 @@ def update(self, A):
         self.lp = self.l[self.p, :]
         return self
 
-    def solve(self, rhs):
+    def solve(self, rhs, x0=None, trans='N'):
         r""" Solves the linear system of equations :math:`\mathbf{A} \mathbf{x} = \mathbf{b}` by forward and backward
         substitution of :math:`\mathbf{x} = \mathbf{L}^{-\text{H}}\mathbf{D}^{-1}\mathbf{L}^{-1}\mathbf{b}` in the
         Hermitian case or as :math:`\mathbf{x} = \mathbf{L}^{-\text{T}}\mathbf{D}^{-1}\mathbf{L}^{-1}\mathbf{b}` in the
         symmetric case.
 
-        The right-hand-side :math:`\mathbf{b}` can be of size ``(N)`` or ``(N, K)``, where ``N`` is the size of matrix
-        :math:`\mathbf{A}` and ``K`` is the number of right-hand sides.
-        """
-        u1 = spla.solve_triangular(self.lp, rhs[self.p], lower=True, unit_diagonal=True)
-        u2 = self.dinv(u1)
-        u = np.zeros_like(rhs, dtype=u2.dtype)
-        u[self.p] = spla.solve_triangular(self.lp, u2, trans='C' if self.hermitian else 'T', lower=True, unit_diagonal=True)
-        return u
-
-    def adjoint(self, rhs):
-        r""" Solves the linear system of equations :math:`\mathbf{A}^\text{H} \mathbf{x} = \mathbf{b}` by forward and
+        The adjoint system of equations :math:`\mathbf{A}^\text{H} \mathbf{x} = \mathbf{b}` is solved by forward and
         backward substitution of
         :math:`\mathbf{x} = \mathbf{L}^{-\text{H}}\mathbf{D}^{-\text{H}}\mathbf{L}^{-1}\mathbf{b}` in the  Hermitian
         case or as :math:`\mathbf{x} = \mathbf{L}^{-\text{H}}\mathbf{D}^{-\text{H}}\mathbf{L}^{-*}\mathbf{b}`
@@ -196,11 +193,34 @@ def adjoint(self, rhs):
         The right-hand-side :math:`\mathbf{b}` can be of size ``(N)`` or ``(N, K)``, where ``N`` is the size of matrix
         :math:`\mathbf{A}` and ``K`` is the number of right-hand sides.
         """
-        if not self.hermitian:
-            u1 = spla.solve_triangular(self.lp, rhs[self.p].conj(), lower=True, unit_diagonal=True).conj()
-        else:
+        if trans == 'N':
+            # Hermitian matrix A: A = L D L^H -> inv(A) = inv(L^H) inv(D) inv(L)
+            # Symmetric matrix A: A = L D L^T -> inv(A) = inv(L^T) inv(D) inv(L)
             u1 = spla.solve_triangular(self.lp, rhs[self.p], lower=True, unit_diagonal=True)
-        u2 = self.dinvH(u1)
-        u = np.zeros_like(rhs, dtype=u2.dtype)
-        u[self.p] = spla.solve_triangular(self.lp, u2, trans='C', lower=True, unit_diagonal=True)
+            u2 = self.dinv(u1)
+            u = np.zeros_like(rhs, dtype=u2.dtype)
+            u[self.p] = spla.solve_triangular(self.lp, u2, trans='C' if self.hermitian else 'T', lower=True, unit_diagonal=True)
+        elif trans == 'T':
+            # Hermitian matrix A^T: A = conj(L) D^T L^T -> inv(A) = inv(L^T) inv(D^T) inv(L^*)
+            # Symmetric matrix A^T: A = L D^T L^T -> inv(A) = inv(L^T) inv(D^T) inv(L)
+            if self.hermitian:
+                u1 = spla.solve_triangular(self.lp, rhs[self.p].conj(), lower=True, unit_diagonal=True).conj()
+            else:
+                u1 = spla.solve_triangular(self.lp, rhs[self.p], lower=True, unit_diagonal=True)
+
+            u2 = self.dinvH(u1.conj()).conj()
+            u = np.zeros_like(rhs, dtype=u2.dtype)
+            u[self.p] = spla.solve_triangular(self.lp, u2, trans='T', lower=True, unit_diagonal=True)
+        elif trans == 'H':
+            # Hermitian matrix A: inv(A^H) = inv(L^H) inv(D^H) inv(L)
+            # Symmetric matrix A: inv(A^H) = inv(L^H) inv(D^H) inv(L^*)
+            if not self.hermitian:
+                u1 = spla.solve_triangular(self.lp, rhs[self.p].conj(), lower=True, unit_diagonal=True).conj()
+            else:
+                u1 = spla.solve_triangular(self.lp, rhs[self.p], lower=True, unit_diagonal=True)
+            u2 = self.dinvH(u1)
+            u = np.zeros_like(rhs, dtype=u2.dtype)
+            u[self.p] = spla.solve_triangular(self.lp, u2, trans='C', lower=True, unit_diagonal=True)
+        else:
+            raise TypeError("Only N, T, and H transposition is possible")
         return u
diff --git a/pymoto/solvers/iterative.py b/pymoto/solvers/iterative.py
new file mode 100644
index 0000000..952a3ba
--- /dev/null
+++ b/pymoto/solvers/iterative.py
@@ -0,0 +1,361 @@
+import warnings
+import time
+import numpy as np
+import scipy.sparse as sps
+from scipy.sparse.linalg import splu, spilu
+from .solvers import LinearSolver
+from .auto_determine import auto_determine_solver
+from pymoto import DomainDefinition
+
+
+class Preconditioner(LinearSolver):
+    """ Abstract base class for preconditioners to inexact solvers """
+    def update(self, A):
+        pass
+
+    def solve(self, rhs, x0=None, trans='N'):
+        return rhs.copy()
+
+
+class DampedJacobi(Preconditioner):
+    r""" Damped Jacobi preconditioner
+    :math:`M = \frac{1}{\omega} D`
+    Args:
+        A (optional): The matrix
+        w (optional): Weight factor :math:`0 < \omega \leq 1`
+    """
+    def __init__(self, A=None, w=1.0):
+        assert 0 < w <= 1, 'w must be between 0 and 1'
+        self.w = w
+        self.D = None
+        super().__init__(A)
+
+    def update(self, A):
+        self.D = A.diagonal()
+
+    def solve(self, rhs, x0=None, trans='N'):
+        if trans == 'N' or trans == 'T':
+            return self.w * (rhs.T/self.D).T
+        elif trans == 'H':
+            return self.w * (rhs.T/self.D.conj()).T
+        else:
+            raise TypeError("Only N, T, or H transposition is possible")
+
+
+class SOR(Preconditioner):
+    r""" Successive over-relaxation preconditioner
+    The matrix :math:`A = L + D + U` is split into a lower triangular, diagonal, and upper triangular part.
+    :math:`M = \left(\frac{D}{\omega} + L\right) \frac{\omega D^{-1}}{2-\omega} \left(\frac{D}{\omega} + U\right)`
+
+    Args:
+        A (optional): The matrix
+        w (optional): Weight factor :math:`0 < \omega < 2`
+    """
+    def __init__(self, A=None, w=1.0):
+        assert 0 < w < 2, 'w must be between 0 and 2'
+        self.w = w
+        self.L = None
+        self.U = None
+        self.Dw = None
+        super().__init__(A)
+
+    def update(self, A):
+        diag = A.diagonal()
+        diagw = sps.diags(diag)/self.w
+        self.L = splu(sps.tril(A, k=-1) + diagw)  # Lower triangular part including diagonal
+        self.U = splu(sps.triu(A, k=1) + diagw)
+
+        self.Dw = diag * (2 - self.w) / self.w
+
+    def solve(self, rhs, x0=None, trans='N'):
+        if trans == 'N':
+            # M = (D/w + L) wD^-1 / (2-w) (D/w + U)
+            # from scipy.sparse.linalg import spsolve_triangular
+            # u1 = spsolve_triangular(self.L, rhs, lower=True, overwrite_A=False)  # Solve triangular is still very slow :(
+            u1 = self.L.solve(rhs)
+            u1 *= self.Dw[:, None]
+            # u2 = spsolve_triangular(self.U, u1, lower=False, overwrite_A=False, overwrite_b=True)
+            u2 = self.U.solve(u1)
+            return u2
+        elif trans == 'T':
+            u1 = self.U.solve(rhs, trans='T')
+            u1 *= self.Dw[:, None]
+            u2 = self.L.solve(u1, trans='T')
+            return u2
+        elif trans == 'H':
+            u1 = self.U.solve(rhs, trans='H')
+            u1 *= self.Dw[:, None].conj()
+            u2 = self.L.solve(u1, trans='H')
+            return u2
+        else:
+            raise TypeError("Only N, T, or H transposition is possible")
+
+
+class ILU(Preconditioner):
+    """ Incomplete LU factorization
+
+    Args:
+        A (optional): The matrix
+        **kwargs (optional): Keyword arguments passed to `scipy.sparse.linalg.spilu`
+    """
+    def __init__(self, A=None, **kwargs):
+        self.kwargs = kwargs
+        self.ilu = None
+        super().__init__(A)
+
+    def update(self, A):
+        self.ilu = spilu(A, **self.kwargs)
+
+    def solve(self, rhs, x0=None, trans='N'):
+        return self.ilu.solve(rhs, trans=trans)
+
+
+class GeometricMultigrid(Preconditioner):
+    """ Geometric multigrid preconditioner
+
+    Args:
+        domain: The `DomainDefinition` used for the geometry
+        A (optional): The matrix
+        inner_level (optional): Inner solver for the coarse grid, for instance, a direct solver or another MG level.
+            The default is an automatically determined direct solver.
+        smoother (optional): Smoother to use to smooth the residual and solution before and after coarse level.
+            The default is `DampedJacobi(w=0.5)`.
+        smooth_steps (optional): Number of smoothing steps to execute
+    """
+    _available_cycles = ['v', 'w']
+
+    def __init__(self, domain: DomainDefinition, A=None, cycle='V', inner_level=None, smoother=None, smooth_steps=5):
+        assert domain.nelx % 2 == 0 and domain.nely % 2 == 0 and domain.nelz % 2 == 0, \
+            f"Domain sizes {domain.nelx, domain.nely, domain.nelz} must be divisible by 2"
+        self.domain = domain
+        self.A = A
+        assert cycle.lower() in self._available_cycles, f"Cycle ({cycle}) is not available. Options are {self._available_cycles}"
+        self.cycle = cycle
+        self.inner_level = None if inner_level is None else inner_level
+        self.smoother = DampedJacobi(w=0.5) if smoother is None else None
+        self.smooth_steps = smooth_steps
+        self.R = None
+        self.sub_domain = DomainDefinition(domain.nelx // 2, domain.nely // 2, domain.nelz // 2,
+                                           domain.unitx * 2, domain.unity * 2, domain.unitz * 2)
+
+        super().__init__(A)
+
+    def update(self, A):
+        if self.R is None:
+            self.setup_interpolation(A)
+        self.A = A
+        self.smoother.update(A)
+        Ac = self.R.T @ A @ self.R
+        if self.inner_level is None:
+            self.inner_level = auto_determine_solver(Ac)
+        self.inner_level.update(Ac)
+
+    def setup_interpolation(self, A):
+        assert A.shape[0] % self.domain.nnodes == 0
+        ndof = int(A.shape[0] / self.domain.nnodes)  # Number of dofs per node
+
+        w = np.ones((3, 3, 3))*0.125
+        w[1, :, :] = 0.25
+        w[:, 1, :] = 0.25
+        w[:, :, 1] = 0.25
+        w[1, 1, :] = 0.5
+        w[1, :, 1] = 0.5
+        w[:, 1, 1] = 0.5
+        w[1, 1, 1] = 1.0
+
+        rows = []
+        cols = []
+        vals = []
+        for i in [-1, 0, 1]:
+            imin, imax = max(-i, 0), min(self.sub_domain.nelx + 1 - i, self.sub_domain.nelx + 1)
+            ix = np.arange(imin, imax)
+            for j in [-1, 0, 1]:
+                jmin, jmax = max(-j, 0), min(self.sub_domain.nely + 1 - j, self.sub_domain.nely + 1)
+                iy = np.arange(jmin, jmax)
+                for k in ([-1, 0, 1] if self.domain.dim == 3 else [0]):
+                    # Coarse node cartesian indices
+                    kmin, kmax = max(-k, 0), min(self.sub_domain.nelz + 1 - k, self.sub_domain.nelz + 1)
+                    iz = np.arange(kmin, kmax)
+                    # Coarse node numbers
+                    nod_c = self.sub_domain.get_nodenumber(*np.meshgrid(ix, iy, iz, indexing='ij')).flatten()
+                    # Fine node numbers with offset
+                    ixc, iyc, izc = ix * 2 + i, iy * 2 + j, iz * 2 + k
+                    nod_f = self.domain.get_nodenumber(*np.meshgrid(ixc, iyc, izc, indexing='ij')).flatten()
+                    for d in range(ndof):
+                        rows.append(nod_f * ndof + d)
+                        cols.append(nod_c * ndof + d)
+                        vals.append(np.ones_like(rows[-1], dtype=w.dtype) * w[1+i, 1+j, 1+k])
+
+        rows = np.concatenate(rows)
+        cols = np.concatenate(cols)
+        vals = np.concatenate(vals)
+        nfine = ndof * self.domain.nnodes
+        ncoarse = ndof * self.sub_domain.nnodes
+        self.R = sps.coo_matrix((vals, (rows, cols)), shape=(nfine, ncoarse))
+        self.R = type(A)(self.R)  # Convert to correct matrix type
+
+    def solve(self, rhs, x0=None, trans='N'):
+        if trans == 'N':
+            A = self.A
+        elif trans == 'T':
+            A = self.A.T
+        elif trans == 'H':
+            A = self.A.conj().T
+        else:
+            raise TypeError("Only N, T, or H transposition is possible")
+
+        # Pre-smoothing
+        if x0 is None:
+            u_f = self.smoother.solve(rhs, trans=trans)
+        else:
+            r = rhs - self.A @ x0
+            u_f = x0 + self.smoother.solve(r, trans=trans)
+        for i in range(self.smooth_steps-1):
+            r = rhs - self.A @ u_f
+            u_f += self.smoother.solve(r, trans=trans)
+
+        r = rhs - A @ u_f
+        # Restrict residual to coarse level
+        r_c = self.R.T @ r
+
+        # Solve at coarse level
+        u_c = self.inner_level.solve(r_c)
+
+        # Interpolate and correct
+        u_f += self.R @ u_c
+
+        # Post-smoothing
+        for i in range(self.smooth_steps):
+            r = rhs - self.A @ u_f
+            u_f += self.smoother.solve(r, trans=trans)
+        return u_f
+
+
+def orth(u, normalize=True, zero_rtol=1e-15):
+    """ Create orthogonal basis from a set of vectors
+
+    Args:
+        u: Set of vectors of size (#dof, #vectors)
+        normalize: Also normalize the basis vectors
+        zero_rtol: Relative tolerance for detection of zero vectors (in case of a rank-deficient basis)
+
+    Returns:
+        v: Orthogonal basis vectors (#dof, #non-zero-vectors)
+    """
+    if u.ndim == 1:
+        return u
+    elif u.ndim > 2:
+        raise TypeError("Only valid for 1D or 2D matrix")
+
+    def dot(a, b):  # Define inner product
+        return a @ b.conj()
+
+    orth_vecs = []
+    for i in range(u.shape[-1]):
+        vi = np.copy(u[..., i])
+        beta_i = dot(vi, vi)
+        for vj in orth_vecs:
+            alpha_ij = dot(vi, vj)
+            alpha_jj = 1.0 if normalize else dot(vj, vj)
+            vi -= vj * alpha_ij / alpha_jj
+        beta_i_new = dot(vi, vi)
+        if beta_i_new / beta_i < zero_rtol:  # Detect zero vector
+            continue
+        if normalize:
+            vi /= np.sqrt(beta_i_new)
+        orth_vecs.append(vi)
+    return np.stack(orth_vecs, axis=-1)
+
+
+class CG(LinearSolver):
+    """ Preconditioned conjugate gradient method
+    Works for positive-definite self-adjoint matrices (:math:`A=A^H`)
+
+    References:
+        Ji & Li (2017), A breakdown-free BCG method. DOI 10.1007/s10543-016-0631-z
+          https://www.cs.odu.edu/~yaohang/portfolio/BIT2017.pdf
+        Shewchuck (1994), Introduction to CG method without the agonzing pain.
+          https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf
+
+    Args:
+        A: The matrix
+        preconditioner: Preconditioner to use
+        tol: Convergence tolerance
+        maxit: Maximum number of iterations
+        restart: Restart every Nth iteration
+        verbosity: Log level
+    """
+    def __init__(self, A=None, preconditioner=Preconditioner(), tol=1e-7, maxit=10000, restart=50, verbosity=0):
+        self.preconditioner = preconditioner
+        self.A = A
+        self.tol = tol
+        self.maxit = maxit
+        self.restart = restart
+        self.verbosity = verbosity
+        super().__init__(A)
+
+    def update(self, A):
+        tstart = time.perf_counter()
+        self.A = A
+        self.preconditioner.update(A)
+        if self.verbosity >= 1:
+            print(f"Preconditioner set up in {np.round(time.perf_counter() - tstart,3)}s")
+
+    def solve(self, rhs, x0=None, trans='N'):
+        if trans == 'N':
+            A = self.A
+        elif trans == 'T':
+            A = self.A.T
+        elif trans == 'H':
+            A = self.A.conj().T
+        else:
+            raise TypeError("Only N, T, or H transposition is possible")
+
+        tstart = time.perf_counter()
+        if rhs.ndim == 1:
+            b = rhs.reshape((rhs.size, 1))
+        else:
+            b = rhs
+        x = np.zeros_like(rhs, dtype=np.result_type(rhs, A)) if x0 is None else x0.copy()
+        if x.ndim == 1:
+            x = x.reshape((x.size, 1))
+
+        r = b - A@x
+        z = self.preconditioner.solve(r, trans=trans)
+        p = orth(z, normalize=True)
+        if self.verbosity >= 2:
+            print(f"Initial residual = {np.linalg.norm(r, axis=0) / np.linalg.norm(b, axis=0)}")
+
+        for i in range(self.maxit):
+            q = A @ p
+            pq = p.conj().T @ q
+            pq_inv = np.linalg.inv(pq)
+            alpha = pq_inv @ (p.conj().T @ r)
+
+            x += p @ alpha
+            if i % 50 == 0:  # Explicit restart
+                r = b - A@x
+            else:
+                r -= q @ alpha
+
+            if self.verbosity >= 2:
+                print(f"i = {i}, residuals = {np.linalg.norm(r, axis=0) / np.linalg.norm(b, axis=0)}")
+
+            tval = np.linalg.norm(r)/np.linalg.norm(b)
+            if tval <= self.tol:
+                break
+
+            z = self.preconditioner.solve(r, trans=trans)
+
+            beta = -pq_inv @ (q.conj().T @ z)
+            p = orth(z + p@beta, normalize=False)
+
+        if tval > self.tol:
+            warnings.warn(f'Maximum iterations ({self.maxit}) reached, with final residual {tval}')
+        elif self.verbosity >= 1:
+            print(f"Converged in {i} iterations and {np.round(time.perf_counter() - tstart, 3)}s, with final residuals {np.linalg.norm(r, axis=0) / np.linalg.norm(b, axis=0)}")
+
+        if rhs.ndim == 1:
+            return x.flatten()
+        else:
+            return x
diff --git a/pymoto/solvers/matrix_checks.py b/pymoto/solvers/matrix_checks.py
new file mode 100644
index 0000000..a94aead
--- /dev/null
+++ b/pymoto/solvers/matrix_checks.py
@@ -0,0 +1,56 @@
+import numpy as np
+import scipy.sparse as sps
+try:
+    import cvxopt
+    _has_cvxopt = True
+except ImportError:
+    _has_cvxopt = False
+
+
+def is_cvxopt_spmatrix(A):
+    """ Checks if the argument is a cvxopt sparse matrix """
+    return isinstance(A, cvxopt.spmatrix) if _has_cvxopt else False
+
+
+def matrix_is_complex(A):
+    """ Checks if the matrix is complex """
+    if is_cvxopt_spmatrix(A):
+        return A.typecode == 'z'
+    else:
+        return np.iscomplexobj(A)
+
+
+def matrix_is_diagonal(A):
+    """ Checks if the matrix is diagonal"""
+    if sps.issparse(A):
+        if isinstance(A, sps.dia_matrix):
+            return len(A.offsets) == 1 and A.offsets[0] == 0
+        else:
+            return np.allclose((A - sps.spdiags(A.diagonal(), 0, *A.shape)).data, 0.0)
+    elif is_cvxopt_spmatrix(A):
+        return max(abs(A.I - A.J)) == 0
+    else:
+        return np.allclose(A, np.diag(np.diag(A)))
+
+
+def matrix_is_symmetric(A):
+    """ Checks whether a matrix is numerically symmetric """
+    if sps.issparse(A):
+        return np.allclose((A-A.T).data, 0)
+    elif is_cvxopt_spmatrix(A):
+        return np.isclose(max(abs(A-A.T)), 0.0)
+    else:
+        return np.allclose(A, A.T)
+
+
+def matrix_is_hermitian(A):
+    """ Checks whether a matrix is numerically Hermitian """
+    if matrix_is_complex(A):
+        if sps.issparse(A):
+            return np.allclose((A-A.T.conj()).data, 0)
+        elif is_cvxopt_spmatrix(A):
+            return np.isclose(max(abs(A-A.ctrans())), 0.0)
+        else:
+            return np.allclose(A, A.T.conj())
+    else:
+        return matrix_is_symmetric(A)
diff --git a/pymoto/solvers/solvers.py b/pymoto/solvers/solvers.py
new file mode 100644
index 0000000..a5b6b64
--- /dev/null
+++ b/pymoto/solvers/solvers.py
@@ -0,0 +1,253 @@
+import warnings
+import numpy as np
+from .matrix_checks import matrix_is_hermitian, matrix_is_symmetric, matrix_is_complex
+
+
+class LinearSolver:
+    """ Base class of all linear solvers
+
+    Keyword Args:
+        A (matrix): Optionally provide a matrix, which is used in :method:`update` right away.
+
+    Attributes:
+        defined (bool): Flag if the solver is able to run, e.g. false if some dependent library is not available
+    """
+
+    defined = True
+    _err_msg = ""
+
+    def __init__(self, A=None):
+        if A is not None:
+            self.update(A)
+
+    def update(self, A):
+        """ Updates with a new matrix of the same structure
+
+        Args:
+            A (matrix): The new matrix of size ``(N, N)``
+
+        Returns:
+            self
+        """
+        raise NotImplementedError(f"Solver not implemented {self._err_msg}")
+
+    def solve(self, rhs, x0=None, trans='N'):
+        r""" Solves the linear system of equations :math:`\mathbf{A} \mathbf{x} = \mathbf{b}`
+
+        Args:
+            rhs: Right hand side :math:`\mathbf{b}` of shape ``(N)`` or ``(N, K)`` for multiple right-hand-sides
+            x0 (optional): Initial guess for the solution
+            trans (optional): Option to transpose matrix
+                'N':   A   @ x == rhs   (default)   Normal matrix
+                'T':   A^T @ x == rhs               Transposed matrix
+                'H':   A^H @ x == rhs               Hermitian transposed matrix (conjugate transposed)
+
+        Returns:
+            Solution vector :math:`\mathbf{x}` of same shape as :math:`\mathbf{b}`
+        """
+        raise NotImplementedError(f"Solver not implemented {self._err_msg}")
+
+    @staticmethod
+    def residual(A, x, b, trans='N'):
+        r""" Calculates the (relative) residual of the linear system of equations
+
+        The residual is calculated as
+        :math:`r = \frac{\left| \mathbf{A} \mathbf{x} - \mathbf{b} \right|}{\left| \mathbf{b} \right|}`
+
+        Args:
+            A: The matrix
+            x: Solution vector
+            b: Right-hand side
+            trans (optional): Matrix tranformation (`N` is normal, `T` is transposed, `H` is hermitian transposed)
+
+        Returns:
+            Residual value
+        """
+        assert x.shape == b.shape
+        if trans == 'N':
+            mat = A
+        elif trans == 'T':
+            mat = A.T
+        elif trans == 'H':
+            mat = A.conj().T
+        else:
+            raise TypeError("Only N, T, or H transposition is possible")
+        return np.linalg.norm(mat@x - b, axis=0) / np.linalg.norm(b, axis=0)
+
+
+class LDAWrapper(LinearSolver):
+    r""" Linear dependency aware solver (LDAS)
+
+    This solver uses previous solutions of the system :math:`\mathbf{A} \mathbf{x} = \mathbf{b}` to reduce computational
+    effort. In case the solution :math:`\mathbf{x}` is linearly dependent on the previous solutions, the solution
+    will be nearly free of cost.
+
+    Args:
+        solver: The internal solver to be used
+        tol (optional): Residual tolerance above which the internal solver is used to add a new solution vector.
+        A (optional): The matrix :math:`\mathbf{A}`
+        symmetric (optional): Flag to indicate a symmetric matrix :math:`A=A^T`
+        hermitian (optional): Flag to indicate a Hermitian matrix :math:`A=A^H`
+
+    References:
+    Koppen, S., van der Kolk, M., van den Boom, S., & Langelaar, M. (2022).
+    Efficient computation of states and sensitivities for compound structural optimisation problems using a Linear Dependency Aware Solver (LDAS).
+    Structural and Multidisciplinary Optimization, 65(9), 273.
+    DOI: 10.1007/s00158-022-03378-8
+    """
+    def __init__(self, solver: LinearSolver, tol=1e-7, A=None, symmetric=None, hermitian=None):
+        self.solver = solver
+        self.tol = tol
+        # Storage for solution vectors (solutions of A x = b)
+        self.x_stored = []
+        self.b_stored = []
+        # Storage for adjoint solution vectors (solutions of A^H x = b)
+        self.xadj_stored = []
+        self.badj_stored = []
+        self.A = None
+        self._did_solve = False  # For debugging purposes
+        self._last_rtol = 0.
+        self.hermitian = hermitian
+        self.symmetric = symmetric
+        self.complex = None
+        super().__init__(A)
+
+    def update(self, A):
+        """ Clear the internal stored solution vectors and update the internal ``solver`` """
+        if self.symmetric is None:
+            self.symmetric = matrix_is_symmetric(A)
+
+        if self.hermitian is None:
+            if not matrix_is_complex(A):
+                self.hermitian = self.symmetric
+            self.hermitian = matrix_is_hermitian(A)
+
+        self.A = A
+        self.x_stored.clear()
+        self.b_stored.clear()
+        self.xadj_stored.clear()
+        self.badj_stored.clear()
+        self.solver.update(A)
+
+    def _do_solve_1rhs(self, A, rhs, x_data, b_data, solve_fn, x0=None):
+        dtype = np.result_type(A, rhs)
+        if rhs.ndim == 1:
+            rhs_loc = np.zeros((rhs.size, 1), dtype=dtype)
+            rhs_loc[:, 0] = rhs
+        else:
+            rhs_loc = np.zeros_like(rhs, dtype=dtype)
+            rhs_loc[:] = rhs
+        sol = np.zeros_like(rhs_loc, dtype=dtype)
+
+        # Check linear dependencies in the rhs using modified Gram-Schmidt
+        for (x, b) in zip(x_data, b_data):
+            assert x.ndim == b.ndim == 1
+            assert x.size == b.size == rhs_loc.shape[0]
+            alpha = rhs_loc.T @ b.conj() / (b.conj() @ b)
+
+            rem_rhs = alpha * b[:, None]
+            add_sol = alpha * x[:, None]
+
+            if np.iscomplexobj(rem_rhs) and not np.iscomplexobj(rhs_loc):
+                if np.linalg.norm(np.imag(rem_rhs)) < 1e-10*np.linalg.norm(np.real(rem_rhs)):
+                    rem_rhs = np.real(rem_rhs)
+                else:
+                    warnings.warn('LDAS: Complex vector cannot be subtracted from real rhs')
+                    continue
+
+            if np.iscomplexobj(add_sol) and not np.iscomplexobj(sol):
+                if np.linalg.norm(np.imag(add_sol)) < 1e-10*np.linalg.norm(np.real(add_sol)):
+                    add_sol = np.real(add_sol)
+                else:
+                    warnings.warn('LDAS: Complex vector cannot be added to real solution')
+                    continue
+
+            rhs_loc -= rem_rhs
+            sol += add_sol
+
+        # Check tolerance
+        self._last_rtol = np.ones(rhs_loc.shape[-1]) if len(x_data) == 0 else self.residual(A, sol, rhs if rhs.ndim > 1 else rhs.reshape(-1, 1))
+        self._did_solve = self._last_rtol > self.tol
+        if np.any(self._did_solve):
+            if x0 is not None:
+                if x0.ndim == 1:
+                    x0_loc = x0.reshape(-1, 1).copy()
+                else:
+                    x0_loc = x0[..., self._did_solve].copy()
+                for x in x_data:
+                    beta = x0_loc.T @ x.conj() / (x.conj() @ x)
+                    x0_loc -= beta * x
+            else:
+                x0_loc = None
+
+            # Calculate a new solution
+            xnew = solve_fn(rhs_loc[..., self._did_solve], x0_loc)
+            self.residual(A, xnew, rhs_loc[..., self._did_solve])
+            sol[..., self._did_solve] += xnew
+
+            # Add to database
+            for i in range(xnew.shape[-1]):
+                # Remove all previous components that are already in the database (orthogonalize)
+                xadd = xnew[..., i]
+                badd = A @ xadd
+                for x, b in zip(x_data, b_data):
+                    beta = badd @ b.conj() / (b.conj() @ b)
+                    badd -= beta * b
+                    xadd -= beta * x
+                bnrm = np.linalg.norm(badd)
+                if not np.isfinite(bnrm) or bnrm == 0:
+                    continue
+                badd /= bnrm
+                xadd /= bnrm
+                x_data.append(xadd)
+                b_data.append(badd)
+
+        if rhs.ndim == 1:
+            return sol.flatten()
+        else:
+            return sol
+
+    def solve(self, rhs, x0=None, trans='N'):
+        r""" Solves the linear system of equations :math:`\mathbf{A} \mathbf{x} = \mathbf{b}` by performing a modified
+        Gram-Schmidt over the previously calculated solutions :math:`\mathbf{U}` and corresponding right-hand-sides
+        :math:`\mathbf{F}`. This is used to construct an approximate solution
+        :math:`\tilde{\mathbf{x}} = \sum_k \alpha_k \mathbf{u}_k` in the subspace of :math:`\mathbf{U}`.
+        If the residual of :math:`\mathbf{A} \tilde{\mathbf{x}} = \mathbf{b}` is above the tolerance, a new solution
+        :math:`\mathbf{u}_{k+1}` will be added to the database such that
+        :math:`\mathbf{x} = \tilde{\mathbf{x}}+\mathbf{u}_{k+1}` is the solution to the system
+        :math:`\mathbf{A} \mathbf{x} = \mathbf{b}`.
+
+        The right-hand-side :math:`\mathbf{b}` can be of size ``(N)`` or ``(N, K)``, where ``N`` is the size of matrix
+        :math:`\mathbf{A}` and ``K`` is the number of right-hand sides.
+        """
+        ''' Required    Symm. matrix    Herm. matrix    Any matrix (uses adjoint storage)
+                        A^T = A         A^H = A
+            A x = b
+            A^T x = b   A x = b         A x^* = b^*     A^H x^* = b^*
+            A^H x = b   A x^* = b^*     A x = b         A^H x = b
+            
+            For symmetric or Hermitian matrices, only the normal storage is required. For any other matrix, the `T` and 
+            `H` mode will require the adjoint storage space.
+        '''
+        if trans not in ['N', 'T', 'H']:
+            raise TypeError("Only N, T, or H transposition is possible")
+
+        # Use adjoint storage?
+        adjoint_mode = trans != 'N' and not (self.symmetric or self.hermitian)
+        # Use conjugation?
+        conj_mode = self.symmetric and trans == 'H' or not self.symmetric and trans == 'T'
+
+        storage = 'H' if adjoint_mode else 'N'
+        rhs = rhs.conj() if conj_mode else rhs
+        if storage == 'N':  # Use the normal storage
+            ret = self._do_solve_1rhs(self.A, rhs,
+                                      self.x_stored, self.b_stored,
+                                      lambda b, x_init: self.solver.solve(b, trans='N', x0=x_init),
+                                      x0=x0)
+        elif storage == 'H':  # Use adjoint storage
+            ret = self._do_solve_1rhs(self.A.conj().T, rhs,
+                                      self.xadj_stored, self.badj_stored,
+                                      lambda b, x_init: self.solver.solve(b, trans='H', x0=x_init),
+                                      x0=x0)
+
+        return ret.conj() if conj_mode else ret
diff --git a/pymoto/common/solvers_sparse.py b/pymoto/solvers/sparse.py
similarity index 86%
rename from pymoto/common/solvers_sparse.py
rename to pymoto/solvers/sparse.py
index 0fc56a1..e4c48a1 100644
--- a/pymoto/common/solvers_sparse.py
+++ b/pymoto/solvers/sparse.py
@@ -2,7 +2,8 @@
 import numpy as np
 import scipy.sparse as sps
 from scipy.sparse import SparseEfficiencyWarning
-from .solvers import matrix_is_hermitian, matrix_is_complex, matrix_is_symmetric, LinearSolver
+from .matrix_checks import matrix_is_hermitian, matrix_is_complex, matrix_is_symmetric
+from .solvers import LinearSolver
 
 # ------------------------------------ Pardiso Solver -----------------------------------
 try:
@@ -126,26 +127,31 @@ def update(self, A):
 
         self._pardiso_solver.factorize(A)
 
-    def solve(self, b):
+    def solve(self, b, x0=None, trans='N'):
         """ solve Ax=b for x
 
         Args:
             A (scipy.sparse.csr.csr_matrix): sparse square CSR matrix , CSC matrix also possible
             b (numpy.ndarray): right-hand side(s), b.shape[0] needs to be the same as A.shape[0]
+            x0 (unused)
+            trans (optional): Indicate which system to solve (Normal, Transposed, or Hermitian transposed)
 
         Returns:
             Solution of the system of linear equations, same shape as input b
         """
-        return self._pardiso_solver.solve(self.A, b)
-
-    def adjoint(self, b):
-        # Cannot use _pardiso_solver.solve because it changes flag 12 internally
-        iparm_prev = self._pardiso_solver.get_iparm(12)
-        self._pardiso_solver.set_iparm(12, int(not iparm_prev))  # Adjoint solver (transpose)
-        b = self._pardiso_solver._check_b(self.A, b)
-        x = self._pardiso_solver._call_pardiso(self.A, b)
-        self._pardiso_solver.set_iparm(12, iparm_prev)  # Revert back to normal solver
-        return x
+        if trans == 'N':
+            return self._pardiso_solver.solve(self.A, b)
+        elif trans == 'T' or trans == 'H':
+            # T and H are the same, because only real matrix is supported
+            # Cannot use _pardiso_solver.solve because it changes flag 12 internally
+            iparm_prev = self._pardiso_solver.get_iparm(12)
+            self._pardiso_solver.set_iparm(12, int(not iparm_prev))  # Adjoint solver (transpose)
+            b = self._pardiso_solver._check_b(self.A, b)
+            x = self._pardiso_solver._call_pardiso(self.A, b)
+            self._pardiso_solver.set_iparm(12, iparm_prev)  # Revert back to normal solver
+            return x
+        else:
+            raise TypeError("Only N, T, or H transposition is possible")
 
     def _print_iparm(self):
         """ Print all iparm settings to console """
@@ -219,17 +225,16 @@ def update(self, A):
         self.inv = splu(A)
         return self
 
-    def solve(self, rhs):
+    def solve(self, rhs, x0=None, trans='N'):
         r""" Solves the linear system of equations :math:`\mathbf{A} \mathbf{x} = \mathbf{b}` by forward and backward
         substitution of :math:`\mathbf{x} = \mathbf{U}^{-1}\mathbf{L}^{-1}\mathbf{b}`.
-        """
-        return self.inv.solve(rhs)
 
-    def adjoint(self, rhs):
-        r""" Solves the linear system of equations :math:`\mathbf{A}^\text{H}\mathbf{x} = \mathbf{b}` by forward and
-        backward substitution of :math:`\mathbf{x} = \mathbf{L}^{-\text{H}}\mathbf{U}^{-\text{H}}\mathbf{b}`.
+        Adjoint system solves the linear system of equations :math:`\mathbf{A}^\text{H}\mathbf{x} = \mathbf{b}` by
+        forward and backward substitution of :math:`\mathbf{x} = \mathbf{L}^{-\text{H}}\mathbf{U}^{-\text{H}}\mathbf{b}`
         """
-        return self.inv.solve(rhs, trans=('H' if self.iscomplex else 'T'))
+        if trans not in ['N', 'T', 'H']:
+            raise TypeError("Only N, T, or H transposition is possible")
+        return self.inv.solve(rhs, trans=trans)
 
 
 # ------------------------------------ Cholesky Solver scikit-sparse -----------------------------------
@@ -281,7 +286,7 @@ def update(self, A):
 
         return self
 
-    def solve(self, rhs):
+    def solve(self, rhs, x0=None, trans='N'):
         r""" Solves the linear system of equations :math:`\mathbf{A} \mathbf{x} = \mathbf{b}` by forward and backward
         substitution of :math:`\mathbf{x} = \mathbf{L}^{-\text{H}}\mathbf{L}^{-1}\mathbf{b}` in case of an
         Hermitian matrix.
@@ -290,10 +295,12 @@ def solve(self, rhs):
         :math:`\mathbf{A}` and ``K`` is the number of right-hand sides.
 
         """
-        return self.inv(rhs)
-
-    def adjoint(self, rhs):
-        return self.solve(rhs)
+        if trans not in ['N', 'T', 'H']:
+            raise TypeError("Only N, T, or H transposition is possible")
+        if trans == 'T':
+            return self.inv(rhs.conj()).conj()
+        else:
+            return self.inv(rhs)
 
 
 # ------------------------------------ Cholesky Solver cvxopt -----------------------------------
@@ -344,15 +351,20 @@ def update(self, A):
 
         return self
 
-    def solve(self, rhs):
+    def solve(self, rhs, x0=None, trans='N'):
         r""" Solves the linear system of equations :math:`\mathbf{A} \mathbf{x} = \mathbf{b}` by forward and backward
         substitution of :math:`\mathbf{x} = \mathbf{L}^{-\text{H}}\mathbf{L}^{-1}\mathbf{b}`. """
+        if trans not in ['N', 'T', 'H']:
+            raise TypeError("Only N, T, or H transposition is possible")
         if rhs.dtype != self._dtype:
             warnings.warn(f"{type(self).__name__}: Type warning: rhs value type ({rhs.dtype}) is converted to {self._dtype}")
-        B = cvxopt.matrix(rhs.astype(self._dtype))
+        if trans == 'T':
+            B = cvxopt.matrix(rhs.conj().astype(self._dtype))
+        else:
+            B = cvxopt.matrix(rhs.astype(self._dtype))
         nrhs = 1 if rhs.ndim == 1 else rhs.shape[1]
+
         cvxopt.cholmod.solve(self.inv, B, nrhs=nrhs)
-        return np.array(B).flatten() if nrhs == 1 else np.array(B)
 
-    def adjoint(self, rhs):
-        return self.solve(rhs)
+        x = np.array(B).flatten() if rhs.ndim == 1 else np.array(B)
+        return x.conj() if trans == 'T' else x


Test patch:
diff --git a/.github/workflows/commit_test.yaml b/.github/workflows/commit_test.yaml
index 2c17292..1d528f2 100644
--- a/.github/workflows/commit_test.yaml
+++ b/.github/workflows/commit_test.yaml
@@ -7,7 +7,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ["3.7", "3.8", "3.9", "3.10", "3.11"]
+        python-version: ["3.7", "3.8", "3.9", "3.10", "3.11", "3.12"]
 
     steps:
       - uses: actions/checkout@v3
@@ -26,7 +26,7 @@ jobs:
         run: |
           python -c "import sys; print(f'Python {sys.version}')"
           python -c "import pip; print(f'pip {pip.__version__}')"
-          python -c "import setuptools; print(f'setuptools {setuptools.__version__}')"
+          #python -c "import setuptools; print(f'setuptools {setuptools.__version__}')"
           python -m pip freeze
       - name: Lint with flake8
         run: |
diff --git a/tests/test_aggregration.py b/tests/test_aggregration.py
new file mode 100644
index 0000000..b70ae07
--- /dev/null
+++ b/tests/test_aggregration.py
@@ -0,0 +1,205 @@
+import unittest
+import pymoto as pym
+import numpy as np
+import numpy.testing as npt
+
+
+class TestActiveSet(unittest.TestCase):
+    def test_inactive(self):
+        """ No options given means all values are passed """
+        np.random.seed(0)
+        x = np.random.rand(5)
+        actset = pym.AggActiveSet()
+        self.assertEqual(x[actset(x)].size, 5)
+
+    def test_all_equal(self):
+        """ In case all values are equal they all must be passed """
+        x = np.ones(5)
+        actset = pym.AggActiveSet(0.1, 0.9, 0.1, 0.9)
+        self.assertEqual(x[actset(x)].size, 5)
+
+    def test_remove_below_value(self):
+        x = np.arange(101)
+        actset = pym.AggActiveSet(lower_rel=0.1)
+        self.assertEqual(min(x[actset(x)]), 10)
+        self.assertEqual(max(x[actset(x)]), 100)
+
+    def test_remove_above_value(self):
+        x = np.arange(101)
+        actset = pym.AggActiveSet(upper_rel=0.9)
+        self.assertEqual(min(x[actset(x)]), 0)
+        self.assertEqual(max(x[actset(x)]), 90)
+
+    def test_remove_above_and_below_value(self):
+        x = np.arange(101)
+        actset = pym.AggActiveSet(lower_rel=0.1, upper_rel=0.9)
+        self.assertEqual(min(x[actset(x)]), 10)
+        self.assertEqual(max(x[actset(x)]), 90)
+
+    def test_remove_lowest_amount(self):
+        x = np.arange(101)
+        actset = pym.AggActiveSet(lower_amt=0.15)
+        self.assertEqual(min(x[actset(x)]), 15)
+        self.assertEqual(max(x[actset(x)]), 100)
+
+    def test_remove_highest_amount(self):
+        x = np.arange(101)
+        actset = pym.AggActiveSet(upper_amt=0.95)
+        self.assertEqual(min(x[actset(x)]), 0)
+        self.assertEqual(max(x[actset(x)]), 95)
+
+    def test_remove_lowest_and_highest_amount(self):
+        x = np.arange(101)
+        actset = pym.AggActiveSet(lower_amt=0.15, upper_amt=0.95)
+        self.assertEqual(min(x[actset(x)]), 15)
+        self.assertEqual(max(x[actset(x)]), 95)
+
+    def test_remove_both_value_and_amount(self):
+        x = np.arange(101)
+        actset = pym.AggActiveSet(lower_rel=0.1, upper_rel=0.9, lower_amt=0.15, upper_amt=0.95)
+        self.assertEqual(min(x[actset(x)]), 15)
+        self.assertEqual(max(x[actset(x)]), 90)
+
+
+class TestScaling(unittest.TestCase):
+    def test_min_no_damping(self):
+        np.random.seed(0)
+        x = np.random.rand(5)
+        sc = pym.AggScaling('min')
+        sf = sc(x, 1.0)
+        self.assertEqual(min(x), sf)
+        x = np.random.rand(5)
+        sf = sc(x, 2.0)
+        self.assertEqual(min(x), 2.0*sf)
+
+    def test_max_no_damping(self):
+        np.random.seed(0)
+        x = np.random.rand(5)
+        sc = pym.AggScaling('max')
+        sf = sc(x, 1.0)
+        self.assertEqual(max(x), sf)
+        x = np.random.rand(5)
+        sf = sc(x, 2.0)
+        self.assertEqual(max(x), 2.0 * sf)
+
+    def test_min_full_damping(self):
+        np.random.seed(0)
+        x = np.random.rand(5)
+        sc = pym.AggScaling('min', damping=1.0)
+        sf0 = sc(x, 1.0)
+        self.assertEqual(min(x), sf0)
+        x = np.random.rand(5)
+        sf = sc(x, 2.0)
+        self.assertEqual(sf, sf0)
+
+    def test_max_full_damping(self):
+        np.random.seed(0)
+        x = np.random.rand(5)
+        sc = pym.AggScaling('max', damping=1.0)
+        sf0 = sc(x, 1.0)
+        self.assertEqual(max(x), sf0)
+        x = np.random.rand(5)
+        sf = sc(x, 2.0)
+        self.assertEqual(sf, sf0)
+
+    def test_min_half_damping(self):
+        np.random.seed(0)
+        x = np.random.rand(5)
+        sc = pym.AggScaling('min', damping=0.5)
+        sf0 = sc(x, 1.0)
+        self.assertEqual(min(x), sf0)
+        x = np.random.rand(5)
+
+        for i in range(10):  # Converges toward minimum
+            sf = sc(x, 1.0)
+            self.assertLess(sf - min(x), sf0 - min(x))
+            sf0 = sf
+
+    def test_max_half_damping(self):
+        np.random.seed(0)
+        x = np.random.rand(5)
+        sc = pym.AggScaling('max', damping=0.5)
+        sf0 = sc(x, 1.0)
+        self.assertEqual(max(x), sf0)
+        x = np.random.rand(5)
+
+        for i in range(10):  # Converges toward minimum
+            sf = sc(x, 1.0)
+            self.assertLess(max(x) - sf, max(x) - sf0)
+            sf0 = sf
+
+
+class TestPNorm(unittest.TestCase):
+    @staticmethod
+    def fd_testfn(x0, dx, df_an, df_fd):
+        npt.assert_allclose(df_an, df_fd, rtol=1e-7, atol=1e-5)
+
+    def test_PNorm_max_fd(self):
+        np.random.seed(0)
+        sx = pym.Signal('x', np.random.rand(100))
+        m = pym.PNorm(sx, p=2)
+
+        pym.finite_difference(m, test_fn=self.fd_testfn)
+
+    def test_PNorm_max1_fd(self):
+        np.random.seed(0)
+        sx = pym.Signal('x', np.random.rand(100))
+        m = pym.PNorm(sx, p=20)
+
+        pym.finite_difference(m, test_fn=self.fd_testfn)
+
+    def test_PNorm_min_fd(self):
+        np.random.seed(0)
+        sx = pym.Signal('x', np.random.rand(100))
+        m = pym.PNorm(sx, p=-4)
+
+        pym.finite_difference(m, test_fn=self.fd_testfn)
+
+
+class TestSoftMinMax(unittest.TestCase):
+    @staticmethod
+    def fd_testfn(x0, dx, df_an, df_fd):
+        npt.assert_allclose(df_an, df_fd, rtol=1e-7, atol=1e-5)
+
+    def test_soft_max_fd(self):
+        np.random.seed(0)
+        sx = pym.Signal('x', np.random.rand(100))
+        m = pym.SoftMinMax(sx, alpha=2)
+
+        pym.finite_difference(m, test_fn=self.fd_testfn)
+
+    def test_soft_min_fd(self):
+        np.random.seed(0)
+        sx = pym.Signal('x', np.random.rand(100))
+        m = pym.SoftMinMax(sx, alpha=-20)
+
+        pym.finite_difference(m, test_fn=self.fd_testfn)
+
+
+class TestKS(unittest.TestCase):
+    @staticmethod
+    def fd_testfn(x0, dx, df_an, df_fd):
+        npt.assert_allclose(df_an, df_fd, rtol=1e-7, atol=1e-5)
+
+    def test_KS_max_fd(self):
+        np.random.seed(0)
+        sx = pym.Signal('x', np.random.rand(100))
+        m = pym.KSFunction(sx, rho=2)
+        m.response()
+        self.assertGreaterEqual(m.sig_out[0].state, max(sx.state))
+
+        pym.finite_difference(m, test_fn=self.fd_testfn)
+
+        # Converge to actual maximum
+        y_2 = m.sig_out[0].state
+        m.rho = 3
+        m.response()
+        y_3 = m.sig_out[0].state
+        self.assertLessEqual(y_3 - max(sx.state), y_2 - max(sx.state))
+
+    def test_KS_min_fd(self):
+        np.random.seed(0)
+        sx = pym.Signal('x', np.random.rand(100))
+        m = pym.KSFunction(sx, rho=-20)
+
+        pym.finite_difference(m, test_fn=self.fd_testfn)
\ No newline at end of file
diff --git a/tests/test_assembly.py b/tests/test_assembly.py
index e17e429..4371985 100644
--- a/tests/test_assembly.py
+++ b/tests/test_assembly.py
@@ -5,6 +5,20 @@
 
 
 class TestAssembleStiffness(unittest.TestCase):
+    def test_rows_columns(self):
+        """ Check if element rows and columns are implemented correctly """
+        domain = pym.DomainDefinition(1, 1)
+        elmat = np.arange(4*4).reshape((4, 4))
+
+        s_x = pym.Signal('x', state=np.ones(domain.nel))
+
+        # Assemble stiffness matrix
+        m = pym.AssembleGeneral(s_x, domain=domain, element_matrix=elmat)
+        m.response()
+        A = m.sig_out[0].state
+
+        npt.assert_allclose(A.toarray(), elmat)
+
     def test_FEA_pure_tensile_2d_one_element(self):
         Lx, Ly, Lz = 0.1, 0.2, 0.3
         domain = pym.DomainDefinition(1, 1, unitx=Lx, unity=Ly, unitz=Lz)
diff --git a/tests/test_complex.py b/tests/test_complex.py
index a85bebd..bd8f5cc 100644
--- a/tests/test_complex.py
+++ b/tests/test_complex.py
@@ -1,7 +1,7 @@
 import unittest
 import pymoto as pym
 import numpy as np
-np.random.seed(0)
+import numpy.testing as npt
 
 
 class ComplexVecDot(pym.Module):
@@ -14,7 +14,12 @@ def _sensitivity(self, dy):
 
 
 class TestComplex(unittest.TestCase):
+    @staticmethod
+    def fd_testfn(x0, dx, df_an, df_fd):
+        npt.assert_allclose(df_an, df_fd, rtol=1e-7, atol=1e-5)
+
     def test_real_to_real_via_complex(self):
+        np.random.seed(0)
         N = 4
         s_ur, s_ui = pym.Signal('ur', np.random.rand(N)), pym.Signal('ui', np.random.rand(N))
         s_vr, s_vi = pym.Signal('vr', np.random.rand(N)), pym.Signal('vi', np.random.rand(N))
@@ -24,9 +29,11 @@ def test_real_to_real_via_complex(self):
         m1 = ComplexVecDot([ma.sig_out[0], mb.sig_out[0]], pym.Signal('uv'))
         m2 = pym.ComplexNorm(m1.sig_out, pym.Signal("|uv|"))
 
-        pym.finite_difference(pym.Network(ma, mb, m1, m2), [s_ur, s_ui, s_vr, s_vi], m2.sig_out, dx=1e-5)
+        fn = pym.Network(ma, mb, m1, m2)
+        pym.finite_difference(fn, [s_ur, s_ui, s_vr, s_vi], m2.sig_out, dx=1e-5, test_fn=self.fd_testfn)
 
     def test_real_to_complex(self):
+        np.random.seed(0)
         N = 4
         s_ur, s_ui = pym.Signal('ur', np.random.rand(N)), pym.Signal('ui', np.random.rand(N) * 0)
         s_vr, s_vi = pym.Signal('vr', np.random.rand(N)), pym.Signal('vi', np.random.rand(N) * 0)
@@ -36,36 +43,40 @@ def test_real_to_complex(self):
         m1 = ComplexVecDot([ma.sig_out[0], mb.sig_out[0]], pym.Signal('u.v'))
 
         pym.finite_difference(pym.Network(ma, mb, m1), [s_ur, s_ui, s_vr, s_vi], m1.sig_out, dx=1e-5, random=False,
-                              keep_zero_structure=False)
+                              keep_zero_structure=False, test_fn=self.fd_testfn)
 
     def test_real_to_complex1(self):
+        np.random.seed(0)
         N = 4
         s_ur, s_ui = pym.Signal('ur', np.random.rand(N)), pym.Signal('ui', np.random.rand(N))
         ma = pym.MakeComplex([s_ur, s_ui], pym.Signal("uc"))
-        pym.finite_difference(ma, dx=1e-5, random=False)
+        pym.finite_difference(ma, dx=1e-5, random=False, test_fn=self.fd_testfn)
 
     def test_complex_to_real(self):
         s_z = pym.Signal('z', np.array([0.3 + 0.4*1j]))
         m1 = pym.RealPart(s_z, pym.Signal('ur'))
         m2 = pym.ImagPart(s_z, pym.Signal('ui'))
-        pym.finite_difference(pym.Network(m1, m2), s_z, [*m1.sig_out, *m2.sig_out], dx=1e-5, random=False)
+        pym.finite_difference(pym.Network(m1, m2), s_z, [*m1.sig_out, *m2.sig_out], dx=1e-5, random=False,
+                              test_fn=self.fd_testfn)
 
     def test_complex_to_real1(self):
         s_z = pym.Signal('z', np.array([0.3 + 0.4*1j]))
         m = pym.ComplexNorm(s_z, pym.Signal('|z|'))
-        pym.finite_difference(m, dx=1e-5, random=False)
+        pym.finite_difference(m, dx=1e-5, random=False, test_fn=self.fd_testfn)
 
     def test_complex_to_real2(self):
+        np.random.seed(0)
         N = 4
         u = np.random.rand(N) + 1j * np.random.rand(N)
         v = np.random.rand(N) + 1j * np.random.rand(N)
         m1 = ComplexVecDot([pym.Signal('u', u), pym.Signal('v', v)], pym.Signal('u'))
         m2 = pym.ComplexNorm(m1.sig_out, pym.Signal("ur"))
-        pym.finite_difference(pym.Network(m1, m2), m1.sig_in, m2.sig_out, dx=1e-5)
+        pym.finite_difference(pym.Network(m1, m2), m1.sig_in, m2.sig_out, dx=1e-5, test_fn=self.fd_testfn)
 
     def test_complex_to_complex(self):
+        np.random.seed(0)
         N = 4
         u = np.random.rand(N) + 1j * np.random.rand(N)
         v = np.random.rand(N) + 1j * np.random.rand(N)
         m1 = ComplexVecDot([pym.Signal('u', u), pym.Signal('v', v)], pym.Signal('u.v'))
-        pym.finite_difference(m1, dx=1e-5)
+        pym.finite_difference(m1, dx=1e-5, test_fn=self.fd_testfn)
diff --git a/tests/test_core.py b/tests/test_core.py
index 27f4881..7017e01 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -553,6 +553,24 @@ def _sensitivity(self, dc):
         self.assertRaises(ValueError, m1.sensitivity)
         # m1.sensitivity()
 
+    def test_identical_sensitivity(self):
+        """ Check to see if identically variables returned in _sensitivity are handled correctly """
+        class Add(pym.Module):
+            def _response(self, A, B):
+                return A + B
+            def _sensitivity(self, dC):
+                return dC, dC
+
+        s_A = pym.Signal('A', np.array([1, 2, 3]))
+        s_B = pym.Signal('B', np.array([4, 5, 6]))
+        m = Add([s_A, s_B])
+        m.response()
+        s_C = m.sig_out[0]
+        s_C.sensitivity = np.array([1, 1, 1])
+        m.sensitivity()
+        self.assertFalse(np.may_share_memory(s_A.sensitivity, s_B.sensitivity))
+        self.assertFalse(np.may_share_memory(s_A.sensitivity, s_C.sensitivity))
+
 
 class TestNetwork(unittest.TestCase):
     def test_correct_network(self):
diff --git a/tests/test_domain.py b/tests/test_domain.py
index 2cefa9a..e14dfe7 100644
--- a/tests/test_domain.py
+++ b/tests/test_domain.py
@@ -4,6 +4,10 @@
 import pymoto as pym
 
 
+def fd_testfn(x0, dx, df_an, df_fd):
+    npt.assert_allclose(df_an, df_fd, rtol=1e-7, atol=1e-5)
+
+
 class TestDomainDefinition(unittest.TestCase):
     def test_node_numbering_2D(self):
         Nx, Ny = 100, 142
@@ -62,6 +66,28 @@ def test_node_numbering_array_3D(self):
         npt.assert_allclose(j_pos, j_nod * 0.2)
         npt.assert_allclose(k_pos, k_nod * 0.3)
 
+    def test_shape_fn_2D(self):
+        unitx, unity = 0.8, 0.3
+        domain = pym.DomainDefinition(1, 1, unitx=unitx, unity=unity)
+
+        for i, n in enumerate(domain.node_numbering):
+            N_chk = np.zeros(domain.elemnodes)
+            N_chk[i] = 1.0
+            pos = np.array([n[0]*unitx/2, n[1]*unity/2])
+            npt.assert_allclose(domain.eval_shape_fun(pos), N_chk)
+    def test_shape_fn_derivatives_2D(self):
+        unitx, unity = 0.8, 0.3
+        domain = pym.DomainDefinition(1, 1, unitx=unitx, unity=unity)
+        pos = np.array([0.2, 0.1])
+
+        class ShapeFn(pym.Module):
+            def _response(self, pos):
+                return domain.eval_shape_fun(pos)
+            def _sensitivity(self, dN):
+                return domain.eval_shape_fun_der(pos) @ dN
+
+        pym.finite_difference(ShapeFn(pym.Signal('pos', state=pos)), test_fn=fd_testfn)
+
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/test_dyadcarrier.py b/tests/test_dyadcarrier.py
index a9ea9f0..c985c0c 100644
--- a/tests/test_dyadcarrier.py
+++ b/tests/test_dyadcarrier.py
@@ -442,6 +442,31 @@ def test_contract_sparse(self):
 
         self.assertAlmostEqual(a.contract(S), np.sum(u1 * diag * v1 + u2 * diag * v2), delta=1e-10)
 
+    def test_contract_sparse_performance(self):
+        # Test contraction with a sparse matrix
+        np.random.seed(0)
+        n = 10000
+
+        a = pym.DyadCarrier([np.random.rand(n) for _ in range(30)], [np.random.rand(n) for _ in range(30)])
+        diag = np.random.rand(1, n)
+        S = spsp.random(n, n, density=200/(n*n)).tocsc()  # Size must be given for scipy==1.7
+        # for i in range(100):
+        #     a.contract(S)
+
+        # Generate random matrices
+        mats = []
+        for _ in range(500):
+            N = np.random.randint(10, 200)
+            rows = np.random.randint(0, n, N)
+            cols = np.random.randint(0, n, N)
+            data = np.random.rand(N)
+            mats.append(spsp.coo_matrix((data, (rows, cols)), shape=(n, n)).tocoo())
+
+        for _ in range(10):
+            a.contract_multi(mats)
+
+        # self.assertAlmostEqual(a.contract(S), np.sum(u1 * diag * v1 + u2 * diag * v2), delta=1e-10)
+
     def test_contract_sparse_slice(self):
         # Test contraction with a sliced sparse matrix
         n = 10
diff --git a/tests/test_element_operations.py b/tests/test_element_operations.py
new file mode 100644
index 0000000..b0e5613
--- /dev/null
+++ b/tests/test_element_operations.py
@@ -0,0 +1,185 @@
+import unittest
+import numpy as np
+import pymoto as pym
+import numpy.testing as npt
+
+np.random.seed(0)
+
+def fd_testfn(x0, dx, df_an, df_fd):
+    npt.assert_allclose(df_an, df_fd, rtol=1e-7, atol=1e-5)
+
+class TestElementOperations(unittest.TestCase):
+    def test_0D_output(self):
+        domains = dict()
+        domains.update(one_element_2D=pym.DomainDefinition(1, 1, unitx=0.7, unity=0.8))
+        domains.update(multi_element_2D=pym.DomainDefinition(2, 2, unitx=0.7, unity=0.8))
+        domains.update(one_element_3D=pym.DomainDefinition(1, 1, 1, unitx=0.7, unity=0.8, unitz=0.9))
+        domains.update(multi_element_3D=pym.DomainDefinition(2, 2, 2, unitx=0.7, unity=0.8, unitz=0.9))
+        for k, domain in domains.items():
+            with self.subTest(k):
+                u = 0.1*np.arange(domain.nnodes * domain.dim)
+                s_u = pym.Signal('u', state=u)
+
+                dofconn = domain.get_dofconnectivity(domain.dim)
+
+                # Average x displacement
+                avg_ux = np.zeros((domain.dim*domain.elemnodes))
+                avg_ux[::2] = 1/domain.elemnodes
+                m_avgx = pym.ElementOperation(s_u, domain=domain, element_matrix=avg_ux)
+                m_avgx.response()
+
+                npt.assert_allclose(m_avgx.sig_out[0].state.shape, (domain.nel))
+                ux_chk = np.sum(u[dofconn[:, ::2]], axis=-1)/domain.elemnodes
+                npt.assert_allclose(m_avgx.sig_out[0].state, ux_chk)
+
+                pym.finite_difference(m_avgx, test_fn=fd_testfn)
+
+                # Average y displacement
+                avg_uy = np.zeros_like(avg_ux)
+                avg_uy[1::2] = 1/domain.elemnodes
+                m_avgy = pym.ElementOperation(s_u, domain=domain, element_matrix=avg_uy)
+                m_avgy.response()
+
+                npt.assert_allclose(m_avgy.sig_out[0].state.shape, (domain.nel))
+                uy_chk = np.sum(u[dofconn[:, 1::2]], axis=-1)/domain.elemnodes
+                npt.assert_allclose(m_avgy.sig_out[0].state, uy_chk)
+
+                pym.finite_difference(m_avgy, test_fn=fd_testfn)
+
+    def test_1D_output(self):
+        domains = dict()
+        domains.update(one_element_2D=pym.DomainDefinition(1, 1, unitx=0.7, unity=0.8))
+        domains.update(multi_element_2D=pym.DomainDefinition(2, 2, unitx=0.7, unity=0.8))
+        domains.update(one_element_3D=pym.DomainDefinition(1, 1, 1, unitx=0.7, unity=0.8, unitz=0.9))
+        domains.update(multi_element_3D=pym.DomainDefinition(2, 2, 2, unitx=0.7, unity=0.8, unitz=0.9))
+        for k, domain in domains.items():
+            with self.subTest(k):
+                u = 0.1*np.arange(domain.nnodes * domain.dim)
+                s_u = pym.Signal('u', state=u)
+
+                dofconn = domain.get_dofconnectivity(domain.dim)
+
+                # Average (x, y) displacement
+                avg_u = np.zeros((2, domain.dim*domain.elemnodes))
+                avg_u[0, 0::2] = 1 / domain.elemnodes
+                avg_u[1, 1::2] = 1 / domain.elemnodes
+                m_avg = pym.ElementOperation(s_u, domain=domain, element_matrix=avg_u)
+                m_avg.response()
+
+                npt.assert_allclose(m_avg.sig_out[0].state.shape, (2, domain.nel))
+                ux_chk = np.sum(u[dofconn[:, ::2]], axis=-1)/domain.elemnodes
+                uy_chk = np.sum(u[dofconn[:, 1::2]], axis=-1) / domain.elemnodes
+                npt.assert_allclose(m_avg.sig_out[0].state[0, :], ux_chk)
+                npt.assert_allclose(m_avg.sig_out[0].state[1, :], uy_chk)
+
+                pym.finite_difference(m_avg, test_fn=fd_testfn)
+
+    def test_2D_output(self):
+        domains = dict()
+        domains.update(one_element_2D=pym.DomainDefinition(1, 1, unitx=0.7, unity=0.8))
+        domains.update(multi_element_2D=pym.DomainDefinition(2, 2, unitx=0.7, unity=0.8))
+        domains.update(one_element_3D=pym.DomainDefinition(1, 1, 1, unitx=0.7, unity=0.8, unitz=0.9))
+        domains.update(multi_element_3D=pym.DomainDefinition(2, 2, 2, unitx=0.7, unity=0.8, unitz=0.9))
+        for k, domain in domains.items():
+            with self.subTest(k):
+                u = 0.1*np.arange(domain.nnodes * domain.dim)
+                s_u = pym.Signal('u', state=u)
+
+                dofconn = domain.get_dofconnectivity(domain.dim)
+
+                # Average (x, y) displacement
+                avg_u = np.zeros((2, 2, domain.dim*domain.elemnodes))
+                avg_u[0, 0, 0::2] = 1 / domain.elemnodes
+                avg_u[1, 1, 1::2] = 1 / domain.elemnodes
+                m_avg = pym.ElementOperation(s_u, domain=domain, element_matrix=avg_u)
+                m_avg.response()
+
+                npt.assert_allclose(m_avg.sig_out[0].state.shape, (2, 2, domain.nel))
+                ux_chk = np.sum(u[dofconn[:, ::2]], axis=-1)/domain.elemnodes
+                uy_chk = np.sum(u[dofconn[:, 1::2]], axis=-1) / domain.elemnodes
+                npt.assert_allclose(m_avg.sig_out[0].state[0, 0, :], ux_chk)
+                npt.assert_allclose(m_avg.sig_out[0].state[0, 1, :], 0)
+                npt.assert_allclose(m_avg.sig_out[0].state[1, 0, :], 0)
+                npt.assert_allclose(m_avg.sig_out[0].state[1, 1, :], uy_chk)
+
+                pym.finite_difference(m_avg, test_fn=fd_testfn)
+
+    def test_3D_output(self):
+        domains = dict()
+        domains.update(one_element_2D=pym.DomainDefinition(1, 1, unitx=0.7, unity=0.8))
+        domains.update(multi_element_2D=pym.DomainDefinition(2, 2, unitx=0.7, unity=0.8))
+        domains.update(one_element_3D=pym.DomainDefinition(1, 1, 1, unitx=0.7, unity=0.8, unitz=0.9))
+        domains.update(multi_element_3D=pym.DomainDefinition(2, 2, 2, unitx=0.7, unity=0.8, unitz=0.9))
+        for k, domain in domains.items():
+            with self.subTest(k):
+                u = 0.1 * np.arange(domain.nnodes * domain.dim)
+                s_u = pym.Signal('u', state=u)
+
+                dofconn = domain.get_dofconnectivity(domain.dim)
+
+                # Displacement of each node per element
+                avg_u = np.zeros((2, 2, 2, domain.dim * domain.elemnodes))
+                for i, n in enumerate(domain.node_numbering):
+                    if i==4:
+                        break
+                    avg_u[(n[0]+1)//2, (n[1]+1)//2, 0, i*domain.dim+0] = 1.2
+                    avg_u[(n[0]+1)//2, (n[1]+1)//2, 1, i*domain.dim+1] = 1.1
+                m_avg = pym.ElementOperation(s_u, domain=domain, element_matrix=avg_u)
+                m_avg.response()
+
+                npt.assert_allclose(m_avg.sig_out[0].state.shape, (2, 2, 2, domain.nel))
+                for i, n in enumerate(domain.node_numbering):
+                    if i==4:
+                        break
+                    npt.assert_allclose(m_avg.sig_out[0].state[(n[0]+1)//2, (n[1]+1)//2, 0, :], u[dofconn[:, i * domain.dim + 0]] * 1.2)
+                    npt.assert_allclose(m_avg.sig_out[0].state[(n[0]+1)//2, (n[1]+1)//2, 1, :], u[dofconn[:, i * domain.dim + 1]] * 1.1)
+
+                pym.finite_difference(m_avg, test_fn=fd_testfn)
+
+
+class TestStressStrain(unittest.TestCase):
+    def test_strain_xx(self):
+        domain = pym.DomainDefinition(1, 1, unitx=0.6, unity=0.7)
+
+        u = np.zeros(domain.elemnodes * domain.dim)
+        n_right = domain.get_nodenumber(1, np.arange(2))
+
+        idx_x = n_right * domain.dim
+        disp = 0.1
+        u[idx_x] = disp
+
+        m_strain = pym.Strain(pym.Signal(state=u), domain=domain)
+        m_strain.response()
+
+        exx_chk = disp/domain.unitx
+        npt.assert_allclose(m_strain.sig_out[0].state[:, 0], np.array([exx_chk, 0, 0]), atol=1e-16)
+
+        E = 67.0
+        m_stress = pym.Stress(pym.Signal(state=u), domain=domain, e_modulus=E, poisson_ratio=0.0)
+        m_stress.response()
+        sxx_chk = E * exx_chk
+        npt.assert_allclose(m_stress.sig_out[0].state[:, 0], np.array([sxx_chk, 0, 0]), atol=1e-16)
+
+    def test_pure_shear(self):
+        domain = pym.DomainDefinition(1, 1, unitx=0.6, unity=0.7)
+
+        u = np.zeros(domain.elemnodes * domain.dim)
+        n_top = domain.get_nodenumber(np.arange(2), 1)
+
+        idx_x = n_top * domain.dim
+        disp = 0.1
+        u[idx_x] = disp
+
+        m_strain = pym.Strain(pym.Signal(state=u), domain=domain)
+        m_strain.response()
+
+        gam_xy_chk = disp/domain.unity
+        npt.assert_allclose(m_strain.sig_out[0].state[:, 0], np.array([0, 0, 2*gam_xy_chk]), atol=1e-16)
+
+        E, nu = 67.0, 0.3
+        m_stress = pym.Stress(pym.Signal(state=u), domain=domain, e_modulus=E, poisson_ratio=nu)
+        m_stress.response()
+
+        G = E / (2*(1+nu))
+        sxy_chk = 2 * G * gam_xy_chk
+        npt.assert_allclose(m_stress.sig_out[0].state[:, 0], np.array([0, 0, sxy_chk]), atol=1e-16)
diff --git a/tests/test_filter.py b/tests/test_filter.py
new file mode 100644
index 0000000..7d7bd6f
--- /dev/null
+++ b/tests/test_filter.py
@@ -0,0 +1,231 @@
+import unittest
+import numpy as np
+import pymoto as pym
+import numpy.testing as npt
+import matplotlib.pyplot as plt
+import time
+
+
+def fd_testfn(x0, dx, df_an, df_fd, rtol=1e-5, atol=1e-5):
+    npt.assert_allclose(df_an, df_fd, rtol=rtol, atol=atol)
+
+
+class TestConvolutionFilter(unittest.TestCase):
+    def test_2D_dot(self):
+        """ Test one element in the middle without boundary effects """
+        np.random.seed(0)
+        domain = pym.DomainDefinition(10, 12)
+
+        ix, iy = 5, 6
+
+        x = np.zeros(domain.nel)
+        x[domain.get_elemnumber(ix, iy)] = 1.0
+        sx = pym.Signal('x', state=x)
+
+        m = pym.FilterConv(sx, domain=domain, radius=2, relative_units=True)
+        m.response()
+
+        y = m.sig_out[0].state
+        w = m.weights
+
+        npt.assert_allclose(y[domain.get_elemnumber(ix, iy)], w[1,1,0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix+1, iy)], w[0,1,0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix, iy+1)], w[1,0,0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix-1, iy)], w[2,1,0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix, iy-1)], w[1,2,0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix+1, iy+1)], w[0, 0, 0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix-1, iy+1)], w[2, 0, 0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix-1, iy-1)], w[2, 2, 0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix+1, iy-1)], w[0, 2, 0])
+
+    def test_2D_edge_xmin_symmetric(self):
+        """ Test one element at the edge to test xmin boundary effect """
+        np.random.seed(0)
+        domain = pym.DomainDefinition(10, 12)
+
+        ix, iy = 0, 6
+
+        x = np.zeros(domain.nel)
+        x[domain.get_elemnumber(ix, iy)] = 1.0
+        sx = pym.Signal('x', state=x)
+
+        m = pym.FilterConv(sx, domain=domain, radius=2, relative_units=True, xmin_bc='symmetric')
+        m.response()
+
+        y = m.sig_out[0].state
+        w = m.weights
+
+        npt.assert_allclose(y[domain.get_elemnumber(ix, iy - 1)], w[1, 2, 0] + w[2, 2, 0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix, iy)],     w[1, 1, 0] + w[2, 1, 0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix, iy + 1)], w[1, 0, 0] + w[2, 0, 0])
+
+        npt.assert_allclose(y[domain.get_elemnumber(ix + 1, iy - 1)], w[0, 2, 0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix + 1, iy)],     w[0, 1, 0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix + 1, iy + 1)], w[0, 0, 0])
+
+    def test_2D_edge_xmin_constval(self):
+        """ Test one element at the edge to test xmin boundary effect """
+        np.random.seed(0)
+        domain = pym.DomainDefinition(10, 12)
+
+        ix, iy = 0, 6
+
+        x = np.zeros(domain.nel)
+        x[domain.get_elemnumber(ix, iy)] = 1.0
+        sx = pym.Signal('x', state=x)
+
+        m = pym.FilterConv(sx, domain=domain, radius=2, relative_units=True, xmin_bc=0.0)
+        m.response()
+
+        y = m.sig_out[0].state
+        w = m.weights
+
+        npt.assert_allclose(y[domain.get_elemnumber(ix, iy - 1)], w[1, 2, 0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix, iy)],     w[1, 1, 0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix, iy + 1)], w[1, 0, 0])
+
+        npt.assert_allclose(y[domain.get_elemnumber(ix + 1, iy - 1)], w[0, 2, 0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix + 1, iy)],     w[0, 1, 0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix + 1, iy + 1)], w[0, 0, 0])
+
+    def test_2D_edge_xmax_symmetric(self):
+        """ Test one element at the edge to test xmin boundary effect """
+        np.random.seed(0)
+        domain = pym.DomainDefinition(10, 12)
+
+        ix, iy = domain.nelx-1, 6
+
+        x = np.zeros(domain.nel)
+        x[domain.get_elemnumber(ix, iy)] = 1.0
+        sx = pym.Signal('x', state=x)
+
+        m = pym.FilterConv(sx, domain=domain, radius=2, relative_units=True, xmin_bc='symmetric')
+        m.response()
+
+        y = m.sig_out[0].state
+        w = m.weights
+
+        npt.assert_allclose(y[domain.get_elemnumber(ix, iy - 1)], w[1, 2, 0] + w[0, 2, 0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix, iy)],     w[1, 1, 0] + w[0, 1, 0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix, iy + 1)], w[1, 0, 0] + w[0, 0, 0])
+
+        npt.assert_allclose(y[domain.get_elemnumber(ix - 1, iy - 1)], w[2, 2, 0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix - 1, iy)],     w[2, 1, 0])
+        npt.assert_allclose(y[domain.get_elemnumber(ix - 1, iy + 1)], w[2, 0, 0])
+
+    def test_2D_fd_symmetric_kernel(self):
+        np.random.seed(0)
+        domain = pym.DomainDefinition(10, 12, unitx=0.5, unity=1.0)
+
+        sx = pym.Signal('x', state=np.random.rand(domain.nel))
+
+        m = pym.FilterConv(sx, domain=domain, radius=5.3, relative_units=False)
+
+        pym.finite_difference(m, test_fn=fd_testfn)
+
+    def test_2D_fd_asymmetric_kernel(self):
+        np.random.seed(0)
+        domain = pym.DomainDefinition(10, 21)
+
+        x = np.random.rand(domain.nel)
+
+        sx = pym.Signal('x', state=x)
+        nx, ny = 3, 1
+        weights = 1+np.arange(nx*ny).reshape((nx, ny, 1))
+        m = pym.FilterConv(sx, domain=domain, weights=weights)
+
+        pym.finite_difference(m, test_fn=fd_testfn)
+
+    def test_3D_dot(self):
+        np.random.seed(0)
+        domain = pym.DomainDefinition(10, 11, 12)
+
+        ix, iy, iz = 5, 6, 7
+
+        x = np.zeros(domain.nel)
+        x[domain.get_elemnumber(ix, iy, iz)] = 1.0
+        sx = pym.Signal('x', state=x)
+
+        m = pym.FilterConv(sx, domain=domain, radius=2, relative_units=True)
+        m.response()
+
+        y = m.sig_out[0].state
+        w = m.weights
+
+        for i in range(3):
+            for j in range(3):
+                for k in range(3):
+                    npt.assert_allclose(y[domain.get_elemnumber(ix-1+i, iy-1+j, iz-1+k)], w[i,j,k])
+
+    def test_3D_dot_symm_at_zmin(self):
+        np.random.seed(0)
+        domain = pym.DomainDefinition(10, 11, 12)
+
+        ix, iy, iz = 5, 6, 0
+
+        x = np.zeros(domain.nel)
+        x[domain.get_elemnumber(ix, iy, iz)] = 1.0
+        sx = pym.Signal('x', state=x)
+
+        m = pym.FilterConv(sx, domain=domain, radius=2, relative_units=True)
+        m.response()
+
+        y = m.sig_out[0].state
+        w = m.weights
+        ysel = y[domain.elements[ix + np.arange(-1, 2), :, :][:, iy + np.arange(-1, 2), :][:, :, :2]]
+
+        # Layer affected by symmetry
+        npt.assert_allclose(ysel[0, 0, 0], w[0, 0, 1] + w[0, 0, 0])
+        npt.assert_allclose(ysel[1, 0, 0], w[1, 0, 1] + w[1, 0, 0])
+        npt.assert_allclose(ysel[2, 0, 0], w[2, 0, 1] + w[2, 0, 0])
+
+        npt.assert_allclose(ysel[0, 1, 0], w[0, 1, 1] + w[0, 1, 0])
+        npt.assert_allclose(ysel[1, 1, 0], w[1, 1, 1] + w[1, 1, 0])
+        npt.assert_allclose(ysel[2, 1, 0], w[2, 1, 1] + w[2, 1, 0])
+
+        npt.assert_allclose(ysel[0, 2, 0], w[0, 2, 1] + w[0, 2, 0])
+        npt.assert_allclose(ysel[1, 2, 0], w[1, 2, 1] + w[1, 2, 0])
+        npt.assert_allclose(ysel[2, 2, 0], w[2, 2, 1] + w[2, 2, 0])
+
+        # Layer unaffected by symmetry
+        npt.assert_allclose(ysel[0, 0, 1], w[0, 0, 2])
+        npt.assert_allclose(ysel[1, 0, 1], w[1, 0, 2])
+        npt.assert_allclose(ysel[2, 0, 1], w[2, 0, 2])
+
+        npt.assert_allclose(ysel[0, 1, 1], w[0, 1, 2])
+        npt.assert_allclose(ysel[1, 1, 1], w[1, 1, 2])
+        npt.assert_allclose(ysel[2, 1, 1], w[2, 1, 2])
+
+        npt.assert_allclose(ysel[0, 2, 1], w[0, 2, 2])
+        npt.assert_allclose(ysel[1, 2, 1], w[1, 2, 2])
+        npt.assert_allclose(ysel[2, 2, 1], w[2, 2, 2])
+
+    def test_3D_symmetric(self):
+        np.random.seed(0)
+        domain = pym.DomainDefinition(10, 12, 13, unitx=0.5, unity=1.0, unitz=1.2)
+
+        sx = pym.Signal('x', state=np.random.rand(domain.nel))
+
+        m = pym.FilterConv(sx, domain=domain, radius=5.3, relative_units=False)
+
+        pym.finite_difference(m, test_fn=fd_testfn)
+
+    def test_3D_symmetric1(self):
+        np.random.seed(0)
+        domain = pym.DomainDefinition(6, 6, 6, unitx=0.5, unity=1.0, unitz=1.2)
+
+        sx = pym.Signal('x', state=np.random.rand(domain.nel))
+        start = time.time()
+        m1 = pym.FilterConv(sx, domain=domain, radius=5, relative_units=True)
+        print(f"Convolution setup = {time.time() - start} s")
+        start = time.time()
+        m1.response()
+        print(f"Convolution elapsed = {time.time() - start} s")
+
+        start = time.time()
+        m2 = pym.DensityFilter(sx, domain=domain, radius=5)
+        print(f"H-matrix setup = {time.time() - start} s")
+        start = time.time()
+        m2.response()
+        print(f"H-matrix elapsed = {time.time() - start} s")
+
diff --git a/tests/test_module_mathgeneral.py b/tests/test_module_mathgeneral.py
index 9efb6f2..765da46 100644
--- a/tests/test_module_mathgeneral.py
+++ b/tests/test_module_mathgeneral.py
@@ -1,5 +1,6 @@
 import unittest
 import numpy as np
+import numpy.testing as npt
 import pymoto as pym
 
 
@@ -183,6 +184,53 @@ def testVec_npScalar_complex_real_input(self):
         def tfn(x0, dx, df_an, df_fd): self.assertTrue(np.allclose(df_an, df_fd, rtol=1e-7, atol=1e-5))
         pym.finite_difference(mod, test_fn=tfn)
 
+    def test_numpy_arrays_with_broadcasting(self):
+        """ Check if broadcasting works for sensitivity calculation """
+        np.random.seed(0)
+        sv1 = pym.Signal("vec", np.random.rand(15))
+        sv2 = pym.Signal("vec", np.random.rand(2, 15))
+        sv3 = pym.Signal("vec", np.random.rand(2, 2, 15))
+        s_scalar = pym.Signal("scalar", 3.5)
+
+        sRes = pym.Signal("result")
+
+        mod = pym.MathGeneral([sv1, sv2, sv3, s_scalar], sRes, expression="inp0*inp1*inp2*inp3")
+        mod.response()
+
+        # Check value of response
+        self.assertEqual(sRes.state.shape, (2, 2, 15))
+        npt.assert_allclose(sRes.state, sv1.state * sv2.state * sv3.state * s_scalar.state)
+
+        # Check sensitivities
+        sRes.sensitivity = np.random.rand(*sRes.state.shape)
+        mod.sensitivity()
+
+        def tfn(x0, dx, df_an, df_fd): self.assertTrue(np.allclose(df_an, df_fd, rtol=1e-7, atol=1e-5))
+        pym.finite_difference(mod, test_fn=tfn)
+
+    def test_numpy_arrays_with_broadcasting1(self):
+        """ Broadcast with singleton axes """
+        np.random.seed(0)
+        sv1 = pym.Signal("vec", np.random.rand(36, 1, 15))
+        sv2 = pym.Signal("vec", np.random.rand(2, 1, 1, 15))
+        sv3 = pym.Signal("vec", np.random.rand(2, 1, 2, 15))
+        s_scalar = pym.Signal("scalar", 3.5)
+
+        sRes = pym.Signal("result")
+
+        mod = pym.MathGeneral([sv1, sv2, sv3, s_scalar], sRes, expression="inp0*inp1*inp2*inp3")
+        mod.response()
+
+        # Check value of response
+        self.assertEqual(sRes.state.shape, (2, 36, 2, 15))
+        npt.assert_allclose(sRes.state, sv1.state * sv2.state * sv3.state * s_scalar.state)
+
+        # Check sensitivities
+        sRes.sensitivity = np.random.rand(*sRes.state.shape)
+        mod.sensitivity()
+
+        def tfn(x0, dx, df_an, df_fd): self.assertTrue(np.allclose(df_an, df_fd, rtol=1e-7, atol=1e-5))
+        pym.finite_difference(mod, test_fn=tfn)
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/test_solvers_dense.py b/tests/test_solvers_dense.py
index de928aa..8b9f60f 100644
--- a/tests/test_solvers_dense.py
+++ b/tests/test_solvers_dense.py
@@ -1,7 +1,8 @@
 import unittest
 import numpy as np
+import numpy.testing as npt
 import pymoto as pym
-from pymoto.modules.linalg import auto_determine_solver
+from pymoto.solvers import auto_determine_solver
 import sys
 import inspect
 np.random.seed(0)
@@ -87,7 +88,7 @@ def test_all_matrices(self):
 
 
 class TestIsComplex(TestGenericUtility):
-    fn = staticmethod(pym.matrix_is_complex)
+    fn = staticmethod(pym.solvers.matrix_is_complex)
     data = [
         (mat_real_diagonal, False),
         (mat_real_symm, False),
@@ -103,7 +104,7 @@ class TestIsComplex(TestGenericUtility):
 
 
 class TestIsDiagonal(TestGenericUtility):
-    fn = staticmethod(pym.matrix_is_diagonal)
+    fn = staticmethod(pym.solvers.matrix_is_diagonal)
     data = [
         (mat_real_diagonal, True),
         (mat_real_symm, False),
@@ -119,7 +120,7 @@ class TestIsDiagonal(TestGenericUtility):
 
 
 class TestIsSymmetric(TestGenericUtility):
-    fn = staticmethod(pym.matrix_is_symmetric)
+    fn = staticmethod(pym.solvers.matrix_is_symmetric)
     data = [
         (mat_real_diagonal, True),
         (mat_real_symm, True),
@@ -135,7 +136,7 @@ class TestIsSymmetric(TestGenericUtility):
 
 
 class TestIsHermitian(TestGenericUtility):
-    fn = staticmethod(pym.matrix_is_hermitian)
+    fn = staticmethod(pym.solvers.matrix_is_hermitian)
     data = [
         (mat_real_diagonal, True),
         (mat_real_symm, True),
@@ -154,16 +155,16 @@ class TestAutoSolver(TestGenericUtility):
     """ Check if the auto_determine returns correct solver types """
     fn = staticmethod(lambda A: type(auto_determine_solver(A)))
     data = [
-        (mat_real_diagonal, pym.SolverDiagonal),
-        (mat_real_symm, pym.SolverDenseCholesky),
-        (mat_real_symm_pos_def, pym.SolverDenseCholesky),
-        (mat_real_asymm, pym.SolverDenseLU),
-        (mat_complex_diagonal, pym.SolverDiagonal),
-        (mat_complex_herm, pym.SolverDenseCholesky),
-        (mat_complex_herm_pos_def, pym.SolverDenseCholesky),
-        (mat_complex_symm, pym.SolverDenseLDL),
-        (mat_complex_symm_pos_def, pym.SolverDenseLDL),
-        (mat_complex_asymm, pym.SolverDenseLU),
+        (mat_real_diagonal, pym.solvers.SolverDiagonal),
+        (mat_real_symm, pym.solvers.SolverDenseCholesky),
+        (mat_real_symm_pos_def, pym.solvers.SolverDenseCholesky),
+        (mat_real_asymm, pym.solvers.SolverDenseLU),
+        (mat_complex_diagonal, pym.solvers.SolverDiagonal),
+        (mat_complex_herm, pym.solvers.SolverDenseCholesky),
+        (mat_complex_herm_pos_def, pym.solvers.SolverDenseCholesky),
+        (mat_complex_symm, pym.solvers.SolverDenseLDL),
+        (mat_complex_symm_pos_def, pym.solvers.SolverDenseLDL),
+        (mat_complex_asymm, pym.solvers.SolverDenseLU),
     ]
 
 
@@ -195,60 +196,109 @@ def get_tags(self):
         start = src.find('matrices')
         start += src[start:].find('=') + 1
         finish = start + src[start:].find(']')
-        return [t.strip(' ').replace('\n', '').replace(' ', '') for t in src[start:finish].strip(' []').split(',')]
+        tags_orig = src[start:finish].strip(' []').split(',')
+        tags_proc = [t.strip(' ').replace('\n', '').replace(' ', '') for t in tags_orig]
+        # Remove commented or empty
+        tags_final = []
+        for t in tags_proc:
+            if len(t) == 0:
+                continue
+            if t[0] == '#':
+                continue
+            tags_final.append(t)
+        return tags_final
 
     def run_solver(self, solver, A, b):
+        atol = 1e-12
         # Get reference solutions
-        Aadj = A.conj().T
-        xref = np.linalg.solve(A, b)
-        xrefadj = np.linalg.solve(Aadj, b)
-
-        # Do the solution and adjoint using provided solver
+        A_N = A
+        A_T = A.T
+        A_H = A.conj().T
+        xref_N = np.linalg.solve(A_N, b)
+        xref_T = np.linalg.solve(A_T, b)
+        xref_H = np.linalg.solve(A_H, b)
+
+        # Do the solutions using provided solver
         solver.update(A)
-        x = solver.solve(b)
-        xadj = solver.adjoint(b)
+        x_N = solver.solve(b, trans='N')
+        x_T = solver.solve(b, trans='T')
+        x_H = solver.solve(b, trans='H')
+
+        # Check solution N
+        npt.assert_allclose(x_N, xref_N)
+        npt.assert_allclose(A_N @ x_N - b, 0.0, atol=atol)
 
-        # Check solution
-        self.assertTrue(np.allclose(x, xref))
-        self.assertTrue(np.allclose(A @ x - b, 0.0))
+        # Check transpose solution T
+        npt.assert_allclose(x_T, xref_T)
+        npt.assert_allclose(A_T @ x_T - b, 0.0, atol=atol)
 
-        # Check adjoint solution
-        self.assertTrue(np.allclose(xadj, xrefadj))
-        self.assertTrue(np.allclose(Aadj @ xadj - b, 0.0))
+        # Check adjoint solution H
+        npt.assert_allclose(x_H, xref_H)
+        npt.assert_allclose(A_H @ x_H - b, 0.0, atol=atol)
 
         # Run with LDAWrapper
-        is_herm = pym.matrix_is_hermitian(A)
-        is_symm = pym.matrix_is_symmetric(A)
-        LDAsolver = pym.LDAWrapper(solver, hermitian=is_herm, symmetric=is_symm)
+        is_herm = pym.solvers.matrix_is_hermitian(A)
+        is_symm = pym.solvers.matrix_is_symmetric(A)
+        is_complex = np.iscomplexobj(b)
+        LDAsolver = pym.solvers.LDAWrapper(solver)
         LDAsolver.update(A)
 
         # First solve (do the solution)
-        x = LDAsolver.solve(b)
-        self.assertTrue(np.allclose(x, xref))
-        self.assertTrue(np.allclose(A @ x - b, 0.0))
-        self.assertTrue(LDAsolver._did_solve)  # the solution must have been done here
+        x_N = LDAsolver.solve(b, trans='N')
+        npt.assert_allclose(x_N, xref_N)
+        npt.assert_allclose(A_N @ x_N - b, 0.0, atol=atol)
+        self.assertTrue(np.all(LDAsolver._did_solve))  # the solution must have been done here
 
         # Second solve
-        x = LDAsolver.solve(b)
-        self.assertTrue(np.allclose(x, xref))
-        self.assertTrue(np.allclose(A @ x - b, 0.0))
-        self.assertFalse(LDAsolver._did_solve)  # the solution is already done
+        x_N = LDAsolver.solve(2.5*b, trans='N')
+        npt.assert_allclose(x_N, 2.5*xref_N)
+        npt.assert_allclose(A_N @ x_N - 2.5*b, 0.0, atol=atol)
+        self.assertFalse(np.any(LDAsolver._did_solve))  # the solution is already done
+
+        # Transpose solve
+        x_T = LDAsolver.solve(3.8*b, trans='T')
+        npt.assert_allclose(x_T, 3.8*xref_T)
+        npt.assert_allclose(A_T @ x_T - 3.8*b, 0.0, atol=atol)
+        if (is_herm and not is_complex) or is_symm:
+            self.assertFalse(np.any(LDAsolver._did_solve))
+        else:
+            self.assertTrue(np.all(LDAsolver._did_solve))
 
-        # Adjoint solve
-        xadj = LDAsolver.adjoint(b)
-        self.assertTrue(np.allclose(xadj, xrefadj))
-        self.assertTrue(np.allclose(Aadj @ xadj - b, 0.0))
+        # Second transpose solve
+        x_T = LDAsolver.solve(4.1*b, trans='T')
+        npt.assert_allclose(x_T, 4.1*xref_T)
+        npt.assert_allclose(A_T @ x_T - 4.1*b, 0.0, atol=atol)
+        self.assertFalse(np.any(LDAsolver._did_solve))
 
-        if is_herm or (is_symm and np.isrealobj(b)):
-            self.assertFalse(LDAsolver._did_solve)
+        # Adjoint solve
+        x_H = LDAsolver.solve(5.3*b, trans='H')
+        npt.assert_allclose(x_H, 5.3*xref_H)
+        npt.assert_allclose(A_H @ x_H - 5.3*b, 0.0, atol=atol)
+        if (is_herm and not is_symm) or (is_symm and not is_complex) or not is_complex:
+            self.assertFalse(np.any(LDAsolver._did_solve))
         else:
-            self.assertTrue(LDAsolver._did_solve)
+            self.assertTrue(np.all(LDAsolver._did_solve))
+
+        # Second adjoint solve
+        x_H = LDAsolver.solve(6.3 * b, trans='H')
+        npt.assert_allclose(x_H, 6.3 * xref_H)
+        npt.assert_allclose(A_H @ x_H - 6.3 * b, 0.0, atol=atol)
+        self.assertFalse(np.any(LDAsolver._did_solve))
+
+        # Normal solve with complex scaling
+        LDAsolver = pym.solvers.LDAWrapper(solver)
+        LDAsolver.update(A)
+        x_N = LDAsolver.solve((6.3 + 1j*3.1) * b, trans='N')
+        npt.assert_allclose(x_N, (6.3 + 1j*3.1) * xref_N)
+        npt.assert_allclose(A_N @ x_N - (6.3 + 1j*3.1) * b, 0.0, atol=atol)
+        self.assertTrue(np.all(LDAsolver._did_solve))
+
+        # Normal solve with real scaling again
+        x_N = LDAsolver.solve(0.9 * b, trans='N')
+        npt.assert_allclose(x_N, 0.9 * xref_N)
+        npt.assert_allclose(A_N @ x_N - 0.9 * b, 0.0, atol=atol)
+        self.assertFalse(np.any(LDAsolver._did_solve))
 
-        # Again adjoint
-        xadj = LDAsolver.adjoint(b)
-        self.assertTrue(np.allclose(xadj, xrefadj))
-        self.assertTrue(np.allclose(Aadj @ xadj - b, 0.0))
-        self.assertFalse(LDAsolver._did_solve)  # the solution is already done
 
     def test_all_matrices(self):
         """ Run the tests on all given matrices """
@@ -257,12 +307,12 @@ def test_all_matrices(self):
             """ Test for different type of right-hand-sides """
             sys.stdout.write(f"Test \"{self.solver.__name__}\" for matrix \"{t}\"\n")
             N = A.shape[0]
-            with self.subTest(msg=f"{t}.real-rhs"):
-                b = np.random.rand(N)
-                self.run_solver(self.solver(), A, b)
-            with self.subTest(msg=f"{t}.complex-rhs"):
-                b = np.random.rand(N) + 1j * np.random.rand(N)
-                self.run_solver(self.solver(), A, b)
+            # with self.subTest(msg=f"{t}.real-rhs"):
+            #     b = np.random.rand(N)
+            #     self.run_solver(self.solver(), A, b)
+            # with self.subTest(msg=f"{t}.complex-rhs"):
+            #     b = np.random.rand(N) + 1j * np.random.rand(N)
+            #     self.run_solver(self.solver(), A, b)
             with self.subTest(msg=f"{t}.multi-real-rhs"):
                 b = np.random.rand(N, 3)
                 self.run_solver(self.solver(), A, b)
@@ -272,7 +322,7 @@ def test_all_matrices(self):
 
 
 class TestDenseDiagonal(GenericTestDenseSolvers):
-    solver = pym.SolverDiagonal
+    solver = pym.solvers.SolverDiagonal
     matrices = [
         mat_real_diagonal,
         mat_complex_diagonal,
@@ -280,7 +330,7 @@ class TestDenseDiagonal(GenericTestDenseSolvers):
 
 
 class TestDenseQR(GenericTestDenseSolvers):
-    solver = pym.SolverDenseQR
+    solver = pym.solvers.SolverDenseQR
     matrices = [
         mat_real_diagonal,
         mat_real_symm,
@@ -296,7 +346,7 @@ class TestDenseQR(GenericTestDenseSolvers):
 
 
 class TestDenseLU(GenericTestDenseSolvers):
-    solver = pym.SolverDenseLU
+    solver = pym.solvers.SolverDenseLU
     matrices = [
         mat_real_diagonal,
         mat_real_symm,
@@ -313,18 +363,18 @@ class TestDenseLU(GenericTestDenseSolvers):
 
 class TestDenseCholesky(GenericTestDenseSolvers):
     # The indefinite matrices use LDL instead as a backup
-    solver = pym.SolverDenseCholesky
+    solver = pym.solvers.SolverDenseCholesky
     matrices = [
         mat_real_diagonal,
-        mat_real_symm,
+        mat_real_symm,  # Uses backup solver
         mat_real_symm_pos_def,
-        mat_complex_herm,
+        mat_complex_herm,  # Uses backup solver
         mat_complex_herm_pos_def,
     ]
 
 
 class TestDenseLDL(GenericTestDenseSolvers):
-    solver = pym.SolverDenseLDL
+    solver = pym.solvers.SolverDenseLDL
     matrices = [
         mat_real_diagonal,
         mat_real_symm,
@@ -377,14 +427,17 @@ def get_tags(self):
         start = src.find('matrices')
         start += src[start:].find('=') + 1
         finish = start + src[start:].find(']')
-
-        def parse_str(t: str):
-            rmv = ['\n', ' ', '(', ',', 'True', 'False']
-            for c in rmv:
-                t = t.replace(c, '')
-            return t
-
-        return [parse_str(t) for t in src[start:finish].strip(' []').split(')')]
+        tags_orig = src[start:finish].strip(' []').split(',')
+        tags_proc = [t.strip(' ').replace('\n', '').replace(' ', '') for t in tags_orig]
+        # Remove commented or empty
+        tags_final = []
+        for t in tags_proc:
+            if len(t) == 0:
+                continue
+            if t[0] == '#':
+                continue
+            tags_final.append(t)
+        return tags_final
 
     def run_solver(self, A, b):
         sA = pym.Signal("A", A)
@@ -392,8 +445,8 @@ def run_solver(self, A, b):
 
         sx = pym.Signal("x")
         fn = pym.Network()
-        symmetry = pym.matrix_is_symmetric(A)
-        hermitian = pym.matrix_is_hermitian(A)
+        symmetry = pym.solvers.matrix_is_symmetric(A)
+        hermitian = pym.solvers.matrix_is_hermitian(A)
 
         if symmetry or hermitian:
             sAsys = pym.Signal("Asym", A)
diff --git a/tests/test_solvers_iterative.py b/tests/test_solvers_iterative.py
new file mode 100644
index 0000000..4441e2f
--- /dev/null
+++ b/tests/test_solvers_iterative.py
@@ -0,0 +1,96 @@
+import inspect
+import pathlib  # For importing files
+import sys
+import unittest
+
+import numpy as np
+import numpy.testing as npt
+
+import pymoto as pym
+
+
+class TestMultigrid(unittest.TestCase):
+    def test_interpolation_2D(self):
+        domain = pym.DomainDefinition(10, 10)
+        mg1 = pym.solvers.GeometricMultigrid(domain)
+        bc_nodes = domain.nodes[:, 0].flatten()
+        bc = np.concatenate([bc_nodes*2, bc_nodes*2+1])
+
+        sx = pym.Signal('x', np.ones(domain.nel))
+        m = pym.AssembleStiffness(sx, domain=domain, bc=bc)
+        m.response()
+        K = m.sig_out[0].state
+
+        mg1.update(K)
+
+        # Test restriction fine -> coarse
+        uf = np.ones(domain.nnodes * 2)
+        uc = mg1.R.T @ uf
+        uc = uc.reshape(int(domain.nelx/2)+1, int(domain.nely/2)+1, 2)
+        npt.assert_allclose(uc[1:-1, :, :][:, 1:-1, :], 4)
+        npt.assert_allclose(uc[0, 1:-1, :], 3)
+        npt.assert_allclose(uc[-1, 1:-1, :], 3)
+        npt.assert_allclose(uc[1:-1, 0, :], 3)
+        npt.assert_allclose(uc[1:-1, -1, :], 3)
+
+        npt.assert_allclose(uc[0, 0, :], 2.25)
+        npt.assert_allclose(uc[0, -1, :], 2.25)
+        npt.assert_allclose(uc[-1, -1, :], 2.25)
+        npt.assert_allclose(uc[-1, 0, :], 2.25)
+
+        # Test interpolation coarse -> fine
+        uc = np.ones(int(domain.nelx / 2 + 1) * int(domain.nely / 2 + 1) * 2)
+        uf = mg1.R @ uc
+        npt.assert_allclose(uf, 1.0)
+
+    def test_interpolation_3D(self):
+        domain = pym.DomainDefinition(10, 10, 10)
+        mg1 = pym.solvers.GeometricMultigrid(domain)
+        bc_nodes = domain.nodes[:, :, 0].flatten()
+        bc = np.concatenate([bc_nodes * 3, bc_nodes * 3 + 1, bc_nodes * 3 + 2])
+
+        sx = pym.Signal('x', np.ones(domain.nel))
+        m = pym.AssembleStiffness(sx, domain=domain, bc=bc)
+        m.response()
+        K = m.sig_out[0].state
+
+        mg1.update(K)
+
+        # Test restriction fine -> coarse
+        uf = np.ones(domain.nnodes * 3)
+        uc = mg1.R.T @ uf
+        uc = uc.reshape(int(domain.nelx/2)+1, int(domain.nely/2)+1, int(domain.nelz/2)+1, 3)
+        npt.assert_allclose(uc[1:-1, :, :, :][:, 1:-1, :, :][:, :, 1:-1, :], 8)
+        npt.assert_allclose(uc[0, 1:-1, :, :][:, 1:-1, :], 6)
+        npt.assert_allclose(uc[-1, 1:-1, :, :][:, 1:-1, :], 6)
+        npt.assert_allclose(uc[1:-1, 0, :, :][:, 1:-1, :], 6)
+        npt.assert_allclose(uc[1:-1, -1, :, :][:, 1:-1, :], 6)
+        npt.assert_allclose(uc[1:-1, :, 0, :][:, 1:-1, :], 6)
+        npt.assert_allclose(uc[1:-1, :, -1, :][:, 1:-1, :], 6)
+
+        npt.assert_allclose(uc[0, 0, 1:-1, :], 4.5)
+        npt.assert_allclose(uc[0, -1, 1:-1, :], 4.5)
+        npt.assert_allclose(uc[-1, -1, 1:-1, :], 4.5)
+        npt.assert_allclose(uc[-1, 0, 1:-1, :], 4.5)
+        npt.assert_allclose(uc[0, 1:-1, 0, :], 4.5)
+        npt.assert_allclose(uc[0, 1:-1, -1, :], 4.5)
+        npt.assert_allclose(uc[-1, 1:-1, -1, :], 4.5)
+        npt.assert_allclose(uc[-1, 1:-1, 0, :], 4.5)
+        npt.assert_allclose(uc[1:-1, 0, 0, :], 4.5)
+        npt.assert_allclose(uc[1:-1, 0, -1, :], 4.5)
+        npt.assert_allclose(uc[1:-1, -1, -1, :], 4.5)
+        npt.assert_allclose(uc[1:-1, -1, 0, :], 4.5)
+
+        npt.assert_allclose(uc[0, 0, 0, :], 3.375)
+        npt.assert_allclose(uc[-1, 0, 0, :], 3.375)
+        npt.assert_allclose(uc[0, -1, 0, :], 3.375)
+        npt.assert_allclose(uc[-1, -1, 0, :], 3.375)
+        npt.assert_allclose(uc[0, 0, -1, :], 3.375)
+        npt.assert_allclose(uc[-1, 0, -1, :], 3.375)
+        npt.assert_allclose(uc[0, -1, -1, :], 3.375)
+        npt.assert_allclose(uc[-1, -1, -1, :], 3.375)
+
+        # Test interpolation coarse -> fine
+        uc = np.ones(int(domain.nelx / 2 + 1) * int(domain.nely / 2 + 1) * int(domain.nelz / 2 + 1) * 3)
+        uf = mg1.R @ uc
+        npt.assert_allclose(uf, 1.0)
diff --git a/tests/test_solvers_sparse.py b/tests/test_solvers_sparse.py
index 6055988..15a6f4d 100644
--- a/tests/test_solvers_sparse.py
+++ b/tests/test_solvers_sparse.py
@@ -5,6 +5,7 @@
 from math import isclose
 
 import numpy as np
+import numpy.testing as npt
 import scipy as sp
 import scipy.sparse as spsp
 import scipy.sparse.linalg as spspla
@@ -49,7 +50,7 @@ def generate_symm_indef():
     sigma = -0.1  # To prevent exact singularity
     As = A - sigma*spsp.eye(A.shape[0])
     # Calculate eigenvalues
-    # solver = pym.SolverSparseLU().update(As)
+    # solver = pym.solvers.SolverSparseLU().update(As)
     # opinv = spspla.LinearOperator(A.shape, matvec=solver.solve, rmatvec=solver.adjoint)
     # w, _ = spspla.eigsh(A, sigma=sigma, k=100, OPinv=opinv, ncv=200)]
     return As.tocoo()
@@ -59,8 +60,8 @@ def generate_dynamic(complex=False, indef=False):
     """ Generate a dynamic stiffness matrix, either complex or not, and indefinite or not"""
     K = load_matrix("bcsstk26.mtx.gz")
     M = load_matrix("bcsstm26.mtx.gz")
-    solver = pym.SolverSparseLU().update(K.tocsc())
-    opinv = spspla.LinearOperator(K.shape, matvec=solver.solve, rmatvec=solver.adjoint)
+    solver = pym.solvers.SolverSparseLU().update(K.tocsc())
+    opinv = spspla.LinearOperator(K.shape, matvec=solver.solve, rmatvec=lambda b: solver.adjoint(b, trans='H'))
     w, _ = spspla.eigsh(K, M=M, sigma=0, k=100, OPinv=opinv, ncv=200)
     eigfreqs = np.sqrt(w)
     if indef:
@@ -77,8 +78,8 @@ def generate_dynamic(complex=False, indef=False):
 def generate_saddlepoint():
     K = load_matrix("bcsstk26.mtx.gz")
     M = load_matrix("bcsstm26.mtx.gz")
-    solver = pym.SolverSparseLU().update(K.tocsc())
-    opinv = spspla.LinearOperator(K.shape, matvec=solver.solve, rmatvec=solver.adjoint)
+    solver = pym.solvers.SolverSparseLU().update(K.tocsc())
+    opinv = spspla.LinearOperator(K.shape, matvec=solver.solve, rmatvec=lambda b: solver.adjoint(b, trans='H'))
     w, V = spspla.eigsh(K, M=M, sigma=0, k=100, OPinv=opinv, ncv=200)
     wi, vi = w[0], V[:, [0]]
     Mv = M@vi
@@ -172,7 +173,7 @@ def test_all_matrices(self):
 
 
 class TestIsComplex(TestGenericUtility):
-    fn = staticmethod(pym.matrix_is_complex)
+    fn = staticmethod(pym.solvers.matrix_is_complex)
     data = [
         (mat_real_diagonal, False),
         (mat_real_spdiag, False),
@@ -192,7 +193,7 @@ class TestIsComplex(TestGenericUtility):
 
 
 class TestIsDiagonal(TestGenericUtility):
-    fn = staticmethod(pym.matrix_is_diagonal)
+    fn = staticmethod(pym.solvers.matrix_is_diagonal)
     data = [
         (mat_real_diagonal, True),
         (mat_real_spdiag, True),
@@ -212,7 +213,7 @@ class TestIsDiagonal(TestGenericUtility):
 
 
 class TestIsSymmetric(TestGenericUtility):
-    fn = staticmethod(pym.matrix_is_symmetric)
+    fn = staticmethod(pym.solvers.matrix_is_symmetric)
     data = [
         (mat_real_diagonal, True),
         (mat_real_spdiag, True),
@@ -232,7 +233,7 @@ class TestIsSymmetric(TestGenericUtility):
 
 
 class TestIsHermitian(TestGenericUtility):
-    fn = staticmethod(pym.matrix_is_hermitian)
+    fn = staticmethod(pym.solvers.matrix_is_hermitian)
     data = [
         (mat_real_diagonal, True),
         (mat_real_spdiag, True),
@@ -263,6 +264,7 @@ class GenericTestSolvers(unittest.TestCase):
     """ Generic test runner for any solver and list of matrices """
     solver = None  # The solver to use
     matrices = []  # The list of matrices the solver should be able to handle
+    kwargs = dict()
 
     @classmethod
     def setUpClass(cls) -> None:
@@ -278,30 +280,40 @@ def get_tags(self):
         start = src.find('matrices')
         start += src[start:].find('=') + 1
         finish = start + src[start:].find(']')
-        return [t.strip(' ').replace('\n', '').replace(' ', '') for t in src[start:finish].strip(' []').split(',')]
+        tags_orig = src[start:finish].strip(' []').split(',')
+        tags_proc = [t.strip(' ').replace('\n', '').replace(' ', '') for t in tags_orig]
+        # Remove commented or empty
+        tags_final = []
+        for t in tags_proc:
+            if len(t) == 0:
+                continue
+            if t[0] == '#':
+                continue
+            tags_final.append(t)
+        return tags_final
 
     def run_solver(self, solver, A, b):
         """ Run the actual test for given solver, matrix and right-hand-side """
         # Reference solution
-        x_ref = pym.SolverSparseLU(A).solve(b)
-        xadj_ref = pym.SolverSparseLU(A.conj().T).solve(b)
-
-        # Calculate the solution and the adjoint solution
+        mats = dict()
+        mats['N'] = A
+        mats['T'] = A.T
+        mats['H'] = A.conj().T
         solver.update(A)
-        x = solver.solve(b)
-        xadj = solver.adjoint(b)
 
-        # Check residual and solution
-        residual = A @ x - b
-        r = np.linalg.norm(residual) / np.linalg.norm(b)
-        self.assertTrue(r < 1e-4, msg=f"Residual is {r}")
-        self.assertTrue(np.allclose(x, x_ref))
+        for k, Ai in mats.items():
+            print(f'Run mode = {k}')
+            xref = pym.solvers.SolverSparseLU(Ai).solve(b)
 
-        # Check adjoint residual and solution
-        residual_adj = (A.conj().T) @ xadj - b
-        r_adj = np.linalg.norm(residual_adj) / np.linalg.norm(b)
-        self.assertTrue(r_adj < 1e-4, msg=f"Adjoint residual is {r_adj} (normal residual = {r})")
-        self.assertTrue(np.allclose(xadj, xadj_ref))
+            # Calculate the solution and the adjoint solution
+            x = solver.solve(b, trans=k)
+
+            # Check residual and solution
+            r = np.linalg.norm(Ai @ x - b) / np.linalg.norm(b)
+            self.assertTrue(r < 1e-4, msg=f"Residual is {r}")
+            npt.assert_allclose(np.real(x), np.real(xref), atol=1e-10, rtol=1e-5)
+            if np.iscomplexobj(x) or np.iscomplexobj(xref):
+                npt.assert_allclose(np.imag(x), np.imag(xref), atol=1e-8, rtol=1e-5)
 
     def run_sparse_tests(self, solver, Acoo, rhs):
         """ Run the sparse tests for different matrix types """
@@ -318,24 +330,28 @@ def test_all_matrices(self):
             N = A.shape[0]
             with self.subTest(msg=f"{t}.real-rhs"):
                 b = np.random.rand(N)
-                self.run_sparse_tests(self.solver(), A, b)
+                self.run_sparse_tests(self.solver(**self.kwargs), A, b)
             with self.subTest(msg=f"{t}.multi-real-rhs"):
                 b = np.random.rand(N, 3)
-                self.run_sparse_tests(self.solver(), A, b)
+                self.run_sparse_tests(self.solver(**self.kwargs), A, b)
+            with self.subTest(msg=f"{t}.multi-real-rhs-lineardependent"):
+                b = np.random.rand(N, 3)
+                b[:, 1] = b[:, 2]
+                self.run_sparse_tests(self.solver(**self.kwargs), A, b)
 
-            if not pym.matrix_is_complex(A):
+            if not pym.solvers.matrix_is_complex(A):
                 continue
 
             with self.subTest(msg=f"{t}.complex-rhs"):
                 b = np.random.rand(N) + 1j * np.random.rand(N)
-                self.run_sparse_tests(self.solver(), A, b)
+                self.run_sparse_tests(self.solver(**self.kwargs), A, b)
             with self.subTest(msg=f"{t}.multi-complex-rhs"):
                 b = np.random.rand(N, 3) + 1j * np.random.rand(N, 3)
-                self.run_sparse_tests(self.solver(), A, b)
+                self.run_sparse_tests(self.solver(**self.kwargs), A, b)
 
 
 class TestSparseDiagonal(GenericTestSolvers):
-    solver = pym.SolverDiagonal
+    solver = pym.solvers.SolverDiagonal
     matrices = [
         mat_real_diagonal,
         mat_real_spdiag,
@@ -345,7 +361,7 @@ class TestSparseDiagonal(GenericTestSolvers):
 
 
 class TestSparseLU(GenericTestSolvers):
-    solver = pym.SolverSparseLU
+    solver = pym.solvers.SolverSparseLU
     matrices = [
         mat_real_diagonal,
         mat_real_spdiag,
@@ -367,7 +383,7 @@ class TestSparseLU(GenericTestSolvers):
 
 class TestCholeskyScikit(GenericTestSolvers):
     # I don't know why, but scikit is able to solve indefinite matrix as well. Maybe they do some LDL inside?
-    solver = pym.SolverSparseCholeskyScikit
+    solver = pym.solvers.SolverSparseCholeskyScikit
     matrices = [
         mat_real_diagonal,
         mat_real_spdiag,
@@ -380,7 +396,7 @@ class TestCholeskyScikit(GenericTestSolvers):
 
 
 class TestCholeskyCVXOPT(GenericTestSolvers):
-    solver = pym.SolverSparseCholeskyCVXOPT
+    solver = pym.solvers.SolverSparseCholeskyCVXOPT
     matrices = [
         mat_real_diagonal,
         mat_real_spdiag,
@@ -392,7 +408,7 @@ class TestCholeskyCVXOPT(GenericTestSolvers):
 
 
 class TestPardiso(GenericTestSolvers):
-    solver = pym.SolverSparsePardiso
+    solver = pym.solvers.SolverSparsePardiso
     matrices = [
         mat_real_diagonal,
         mat_real_spdiag,
@@ -410,5 +426,29 @@ class TestPardiso(GenericTestSolvers):
     ]
 
 
+class TestSparseCG(GenericTestSolvers):
+    solver = pym.solvers.CG
+    # kwargs = dict(preconditioner=pym.solvers.ILU(), verbosity=1, tol=1e-9)
+    kwargs = dict(preconditioner=pym.solvers.SOR(w=1.0), verbosity=1, tol=1e-9)
+    # kwargs = dict(preconditioner=pym.solvers.DampedJacobi(), verbosity=1, tol=1e-9)
+    matrices = [
+        mat_real_diagonal,
+        mat_real_spdiag,
+        mat_real_symm_pos_def,
+        mat_real_symm_pos_def_dynamic,
+        # mat_real_symm_indef,
+        # mat_real_symm_indef_dynamic,
+        # mat_real_asymm,
+        # mat_real_symm_saddle,
+        mat_complex_diagonal,
+        mat_complex_spdiag,
+        # mat_complex_symm_pos_def_dynamic, # Works but needs many iterations
+        # mat_complex_symm_indef_dynamic,
+        mat_complex_hermitian_1,
+        # mat_complex_hermitian_indef,
+        mat_complex_hermitian_pos_def,
+    ]
+
+
 if __name__ == '__main__':
     unittest.main()