next up previous
Next: Conclusions Up: Applications Previous: Assembling the matrix

Parallel Jacobi preconditioning

When preconditioning a matrix, for example in a Conjugate Gradients solver, one may choose the Jacobi scheme for preconditioning. The preconditioned vector $\tilde v$ is computed from the vector v using the following relationship:

\begin{displaymath}\tilde v_i = \frac 1{a_{ii}} v_i,
\end{displaymath}

where aii are the diagonal elements of the matrix which we are presently inverting. As is obvious, the result of preconditioning one element of v is entirely independent of all other elements, so this operation is trivially parallelizable. In practice, this is done by splitting the interval [0,n)into equal parts $[n_i,n_{i+1}), i=0,\dots,p-1$, where n is the size of the matrix, and p is the number of processors. Obviously, n0=0, np=n, and ni<ni+1.

Just like for splitting a range of iterators using the function split_range used above, there is a function

    vector<pair<unsigned int, unsigned int> >
    split_interval (const unsigned int &begin, const unsigned int &end,
                    const unsigned int n_intervals);
that splits the interval [begin,end) into n_intervals equal parts. This function will be used to assign each processor its share of elements vi.

Furthermore, we will use some functionality provided by the MultithreadInfo class in deal.II. Upon start-up of the library, the static variable multithread_info.n_cpus is set to the number of processors in the computer the program is presently running on. multithread_info is a global variable of type MultithreadInfo available in all parts of the library. Furthermore, there is a variable multithread_info.n_default_threads, which by default is set to n_cpus, but which can be changed by the user; it denotes the default number of threads which the library shall use whenever multi-threading is implemented for some operation. We will use this variable to decide how many threads shall be used to precondition the vector.

The implementation of the preconditioning function then looks like this:

                        // define an abbreviatory data type for an interval
    typedef pair<unsigned int, unsigned int> Interval;

    void Preconditioner::precondition_jacobi (const Matrix &m,
                                              const Vector &v,
                                              Vector       &v_tilde) {
                            // define an abbreviation to the number
                            // of threads which we will use
      const unsigned int n_threads = multithread_info.n_default_threads;
                            // first split the interval into equal pieces
      vector<Interval> intervals = Threads::split_interval (0, m.rows(),
                                                            n_threads);
 
                            // then define a thread manager
      ACE_Thread_Manager thread_manager;
                            // and finally start all the threads:
      for (unsigned int i=0; i<n_threads; ++i)
        Threads::spawn (thread_manager,
                   Threads::encapsulate (&Preconditioner::threaded_jacobi)
                          .collect_args (this, m, v, v_tilde, intervals[i]));
  
                            // wait for all the threads to finish
      thread_manager.wait ();
    };


    void Preconditioner::threaded_jacobi (const Matrix   &m,
                                          const Vector   &v,
                                          Vector         &v_tilde,
                                          const Interval &interval) {
                           // apply the preconditioner in the given interval
      for (unsigned int i=interval.first; i<interval.second; ++i)
        v_tilde(i) = v(i) / m(i,i);
    };

It is noted, however, that more practical preconditioners are usually not easily parallelized. However, matrix-vector and vector-vector operations can often be reduced to independent parts and can then be implemented using multiple threads.


next up previous
Next: Conclusions Up: Applications Previous: Assembling the matrix
Wolfgang Bangerth
2000-04-20