Ungeblockte LU-Zerlegung

Bei der ungeblockten LU-Zerlegung mit partieller Pivotisierung soll die \(m \times n\) Matrix \(A\) mit der Zerlegung kompakt überschrieben werden. Die Zeilenvertauschungen sollen in einem Pivot-Vektor festgehalten werden. Nur wenn die Matrix \(U\) bei der Faktorisierung exakt singulär ist, soll das Verfahren abbrechen. In diesem Fall soll die Funktion den kleinsten Index \(i\) mit \(U(i,i)\) zurückliefern. Andernfalls wird durch einen Rückgabewert \(-1\) signalisiert, dass die Zerlegung erfolgreich durchgeführt werden konnte.

Aufgabe

Implementiert die ungeblockte LU-Zerlegung hpc::ulmlapack::getf2 nach dem in der Vorlesung vorgestellten Algorithmus. Dabei seien die Matrizen aber bei Null beginnend indiziert:

\(mn \leftarrow \min\{m, n\}\)
\(\text{For}\; j=0,\dots,mn\)
- \(i:\;|x_i| = \max_{k = j,\dots,m-1} |x_k|\)
- \(\text{If}\; i \neq j\)
  - \(A_{i,k} \leftrightarrow A_{j,k} \quad (0\leq k < n)\)
- \(p_j = i\)
- \(\text{If}\; A_{j,j} = 0\)
  - \(\text{return:}\; j\)
- \(A_{k,j} \leftarrow \frac{A_{k,j}}{A_{j,j}} \quad (j+1 \leq k < m)\)
- \(A_{k, l} \leftarrow A_{k, l} - A_{k,j} A_{j,l} \quad (j+1 \leq k,l < m)\)
\(\text{return:}\; -1\)

Natürlich sollen dabei möglichst die auf der vorigen Seite entwickelten Bausteine verwendet werden. Als Signatur legen wir fest:

#ifndef HPC_ULMLAPACK_GETF2_H
#define HPC_ULMLAPACK_GETF2_H 1

namespace hpc { namespace ulmlapack {

template <typename Index, typename TA, typename TP>
Index
getf2(Index m, Index n, TA *A, Index incRowA, Index incColA, TP *p, Index incP)
{
    /* ... */
}

} } // namespace ulmblas, hpc

#endif // HPC_ULMLAPACK_GETF2_H

Zum Test kann folgendes Programm benutzt werden:

session17/hpc/tests/lu.cc

#include <cassert>
#include <random>
#include <hpc/aux/walltime.h>
#include <hpc/matvec/copy.h>
#include <hpc/matvec/gematrix.h>
#include <hpc/matvec/mm.h>
#include <hpc/matvec/print.h>
#include <hpc/ulmlapack/getf2.h>

using namespace hpc::matvec;

template <typename T, typename I>
double
checkLU(const GeMatrixView<T, I> &A,
        GeMatrixView<T, I> &A_,
        const GeMatrixView<I, I> &p)
{
    assert(A.numRows==A_.numRows);
    assert(A.numCols==A_.numCols);
    assert(A.numCols>=A.numRows);
    assert(A.numRows==p.numRows);

    typedef typename GeMatrix<T, I>::Index  Index;

    Index m = A.numRows;
    Index n = A.numRows;

    // Copy L form A_ and set A_ to U
    GeMatrix<T, I> L(m, m);
    for (Index j=0; j<m; ++j) {
        for (Index i=0; i<m; ++i) {
            if (i==j) {
                L(i,i) = 1;
            } else if (i>j) {
                L(i,j)  = A_(i,j);
                A_(i,j) = 0;
            } else if (j>i) {
                L(i,j)  = 0;
            }
        }
    }

    // Compute LU = L*U
    GeMatrix<T, I> LU(m, n);
    for (Index j=0; j<n; ++j) {
        for (Index i=0; i<m; ++i) {
            LU(i,j) = 0;
        }
    }
    hpc::matvec::mm(T(1), L, A_, T(0), LU);

    // Apply P
    for (Index i=m; i>0; --i) {
        if (i-1!=p(i-1,0)) {
            hpc::ulmblas::swap(n,
                               &LU(i-1,0), LU.incCol,
                               &LU(p(i-1,0),0), LU.incCol);
        }
    }

    double diff = 0;
    for (Index j=0; j<n; ++j) {
        for (Index i=0; i<m; ++i) {
            diff += std::abs(A(i,j)-LU(i,j));
        }
    }
    return diff;
}

//
//  Random initializer for general matrices: real and complex valued
//
template <typename Index, typename T>
void
randomInit(Index m, Index n, T *A, Index incRowA, Index incColA)
{
    std::random_device                  random;
    std::default_random_engine          mt(random());
    std::uniform_real_distribution<T>   uniform(-100,100);

    for (Index i=0; i<m; ++i) {
        for (Index j=0; j<n; ++j) {
            A[i*incRowA+j*incColA] = uniform(mt);
        }
    }
}

int
main()
{
    const std::ptrdiff_t   m_min = 10;
    const std::ptrdiff_t   n_min = 10;

    const std::ptrdiff_t   m_max = 6000;
    const std::ptrdiff_t   n_max = 8000;

    const std::ptrdiff_t   m_inc = 10;
    const std::ptrdiff_t   n_inc = 10;

    GeMatrix<double, std::ptrdiff_t>      A(m_max, n_max);
    GeMatrix<double, std::ptrdiff_t>      A_org(m_max, n_max);
    GeMatrix<std::ptrdiff_t, std::ptrdiff_t>        p(m_max, 1);
    hpc::aux::WallTime<double>  timer;

    for (std::ptrdiff_t m=m_min, n=n_min; m<=m_max && n<=n_max; m+=m_inc, n+=n_inc) {

        randomInit(m, n, A.data, A.incRow, A.incCol);
        hpc::ulmblas::gecopy(m, n,
                             A.data, A.incRow, A.incCol,
                             A_org.data, A_org.incRow, A_org.incCol);
        //print(A_org, "A_org");

        timer.tic();
        std::ptrdiff_t info = hpc::ulmlapack::getf2(m, n,
                                                    A.data, A.incRow, A.incCol,
                                                    p.data, p.incRow);
        double t = timer.toc();
        //print(A, "A");
        //print(p, "p");

        auto A_ = A_org(0, 0, m, n);
        auto LU = A(0, 0, m, n);
        auto p_ = p(0, 0, m, 1);

        double diff = checkLU(A_, LU, p_);
        //print(A, "LU");
        std::printf("%4ld %4ld %4ld %16.3e %16.5lf\n", m, n, info, diff, t);
    }
}