The FEniCS interface allows different ways to access the core
functionality, ranging from very high-level to low-level access. So
far, we have mostly used the high-level call `solve(a == L, u, bc)`

to
solve a variational problem `a == L`

with a certain boundary condition
`bc`

. However, sometimes you may need more fine-grained control of
the solution process. In particular, the call to `solve`

will create
certain objects that are thrown away after the solution has been
computed, and it may be practical or efficient to *reuse* those
objects.

In this section, we will look at an alternative interface to solving
linear variational problems in FEniCS, which may be preferable in
many situations compared to the high-level `solve`

function interface.
This interface uses the two classes `LinearVariationalProblem`

and
`LinearVariationalSolver`

. Using this interface, the equivalent of
`solve(a == L, u, bc)`

looks as follows:

```
u = Function(V)
problem = LinearVariationalProblem(a, L, u, bc)
solver = LinearVariationalSolver(problem)
solver.solve()
```

Many FEniCS objects have an attribute `parameters`

, similar to
the global `parameters`

database,
but local to the object. Here, `solver.parameters`

play that
role. Setting the CG method with ILU preconditioning as the solution
method and specifying solver-specific parameters can be done
like this:

```
solver.parameters.linear_solver = 'gmres'
solver.parameters.preconditioner = 'ilu'
prm = solver.parameters.krylov_solver # short form
prm.absolute_tolerance = 1E-7
prm.relative_tolerance = 1E-4
prm.maximum_iterations = 1000
```

Settings in the global `parameters`

database are
propagated to parameter sets in individual objects, with the
possibility of being overwritten as above. Note that global parameter
values can only affect local parameter values if set before the time
of creation of the local object. Thus, changing the value of the
tolerance in the global parameter database will not affect the
parameters for already created solvers.

As we saw already in the section The Navier--Stokes equations, linear variational
problems can be assembled explicitly in FEniCS into matrices and
vectors using the `assemble`

function. This allows even more
fine-grained control of the solution process compared to using the
high-level `solve`

function or using the classes
`LinearVariationalProblem`

and
`LinearVariationalSolver`

. We will now look more closely into how to
use the `assemble`

function and how to combine this with low-level
calls for solving the assembled linear systems.

Given a variational problem \( a(u,v)=L(v) \), the discrete solution \( u \) is computed by inserting \( u=\sum_{j=1}^N U_j \phi_j \) into \( a(u,v) \) and demanding \( a(u,v)=L(v) \) to be fulfilled for \( N \) test functions \( \hat\phi_1,\ldots,\hat\phi_N \). This implies $$ \begin{equation*} \sum_{j=1}^N a(\phi_j,\hat\phi_i) U_j = L(\hat\phi_i),\quad i=1,\ldots,N, \end{equation*} $$ which is nothing but a linear system, $$ \begin{equation*} AU = b, \end{equation*} $$ where the entries of \( A \) and \( b \) are given by $$ \begin{align*} A_{ij} &= a(\phi_j, \hat{\phi}_i), \\ b_i &= L(\hat\phi_i)\tp \end{align*} $$

The examples so far have specified the left- and right-hand sides of
the variational formulation and then asked FEniCS to assemble the
linear system and solve it. An alternative is to explicitly call
functions for assembling the coefficient matrix \( A \) and the right-hand
side vector \( b \), and then solve the linear system \( AU=b \) for
the vector \( U \). Instead of `solve(a == L, U, b)`

we now write

```
A = assemble(a)
b = assemble(L)
bc.apply(A, b)
u = Function(V)
U = u.vector()
solve(A, U, b)
```

The variables `a`

and `L`

are the same as before; that is, `a`

refers
to the bilinear form involving a `TrialFunction`

object `u`

and a `TestFunction`

object `v`

, and `L`

involves the same `TestFunction`

object `v`

. From `a`

and `L`

, the `assemble`

function can compute
\( A \) and \( b \).

Creating the linear system explicitly in a program can have some advantages in more advanced problem settings. For example, \( A \) may be constant throughout a time-dependent simulation, so we can avoid recalculating \( A \) at every time level and save a significant amount of simulation time.

The matrix \( A \) and vector \( b \) are first assembled without
incorporating essential (Dirichlet) boundary conditions. Thereafter,
the call `bc.apply(A, b)`

performs the necessary modifications of the
linear system such that `u`

is guaranteed to equal the prescribed
boundary values. When we have multiple Dirichlet conditions stored in
a list `bcs`

, we must apply each condition in `bcs`

to the system:

```
for bc in bcs:
bc.apply(A, b)
# Alternative syntax using list comprehension
[bc.apply(A, b) for bc in bcs]
```

Alternatively, we can use the function `assemble_system`

, which takes
the boundary conditions into account when assembling the linear
system. This method preserves the symmetry of the linear system for a
symmetric bilinear form. Even if the matrix `A`

that comes out
of the call to `assemble`

is symmetric (for a symmetric bilinear form
`a`

), the call to `bc.apply`

will break the symmetry. Preserving the
symmetry of a variational problem is important when using particular
linear solvers designed for symmetric systems, such as the conjugate
gradient method.

Once the linear system has been assembled, we need to compute the
solution \( U=A^{-1}b \) and store the solution \( U \) in the vector
`U = u.vector()`

. In the same way as linear variational problems can be
programmed using different interfaces in FEniCS—the high-level
`solve`

function, the class `LinearVariationalSolver`

, and the
low-level `assemble`

function—linear systems can also be programmed
using different interfaces in FEniCS. The high-level interface to
solving a linear system in FEniCS is also named `solve`

:

```
solve(A, U, b)
```

By default, `solve(A, U, b)`

uses sparse LU decomposition to compute
the solution. Specification of an iterative solver and preconditioner
can be made through two optional arguments:

```
solve(A, U, b, 'cg', 'ilu')
```

Appropriate names of solvers and preconditioners are found in the section List of linear solver methods and preconditioners.

This high-level interface is useful for many applications, but
sometimes more fine-grained control is needed. One can then create one
or more `KrylovSolver`

objects that are then used to solve linear
systems. Each different solver object can have its own set of
parameters and selection of iterative method and preconditioner. Here
is an example:

```
solver = KrylovSolver('cg', 'ilu')
prm = solver.parameters
prm.absolute_tolerance = 1E-7
prm.relative_tolerance = 1E-4
prm.maximum_iterations = 1000
u = Function(V)
U = u.vector()
solver.solve(A, U, b)
```

The function `solver_linalg`

in the program file
`ft10_poisson_extended.py`
implements such a solver.

The choice of start vector for the iterations in a linear solver is
often important. By default, the values of `u`

and thus the vector ```
U
= u.vector()
```

will be initialized to zero. If we instead wanted to
initialize `U`

with random numbers in the interval \( [-100,100] \) this
can be done as follows:

```
n = u.vector().array().size
U = u.vector()
U[:] = numpy.random.uniform(-100, 100, n)
solver.parameters.nonzero_initial_guess = True
solver.solve(A, U, b)
```

Note that we must both turn off the default behavior of setting the start
vector ("initial guess") to zero, and also set the values of the
vector `U`

to nonzero values. This is useful if we happen to
know a good initial guess for the solution.

Using a nonzero initial guess can be particularly important for
time-dependent problems or when solving a linear system as part of a
nonlinear iteration, since then the previous solution vector `U`

will
often be a good initial guess for the solution in the next time step
or iteration. In this case, the values in the vector `U`

will
naturally be initialized with the previous solution vector (if we just
used it to solve a linear system), so the only extra step necessary is
to set the parameter `nonzero_initial_guess`

to `True`

.

When calling `A = assemble(a)`

and `b = assemble(L)`

, the object `A`

will be of type `Matrix`

, while `b`

and `u.vector()`

are of type
`Vector`

. To examine the values, we may convert the matrix and vector
data to `numpy`

arrays by calling the `array`

method as shown
before. For example, if you wonder how essential boundary conditions are
incorporated into linear systems, you can print out `A`

and `b`

before and after the `bc.apply(A, b)`

call:

```
A = assemble(a)
b = assemble(L)
if mesh.num_cells() < 16: # print for small meshes only
print(A.array())
print(b.array())
bc.apply(A, b)
if mesh.num_cells() < 16:
print(A.array())
print(b.array())
```

With access to the elements in `A`

through a `numpy`

array, we can easily
perform computations on this matrix, such as computing the eigenvalues
(using the `eig`

function in `numpy.linalg`

). We can alternatively dump
`A.array()`

and `b.array()`

to file in MATLAB format and invoke
MATLAB or Octave to analyze the linear system.
Dumping the arrays to MATLAB format is done by

```
import scipy.io
scipy.io.savemat('Ab.mat', {'A': A.array(), 'b': b.array()})
```

Writing `load Ab.mat`

in MATLAB or Octave will then make
the array variables `A`

and `b`

available for computations.

Matrix processing in Python or MATLAB/Octave is only feasible for
small PDE problems since the `numpy`

arrays or matrices in MATLAB file
format are dense matrices. FEniCS also has an interface to the
eigensolver package SLEPc, which is the preferred tool for computing
the eigenvalues of large, sparse matrices of the type encountered in
PDE problems (see `demo/documented/eigenvalue/python/`

in the
FEniCS/DOLFIN source code tree for a demo).