getrf_level0.hpp

Go to the documentation of this file.

 
//
// Copyright (c) 2021-2023, University of Colorado Denver. All rights reserved.
//
// This file is part of <T>LAPACK.
// <T>LAPACK is free software: you can redistribute it and/or modify it under
// the terms of the BSD 3-Clause license. See the accompanying LICENSE file.
 
#ifndef TLAPACK_GETRF_LV0_HH
#define TLAPACK_GETRF_LV0_HH
 
#include "tlapack/base/utils.hpp"
namespace tlapack {
 
template <TLAPACK_MATRIX matrix_t, TLAPACK_VECTOR piv_t>
int getrf_level0(matrix_t& A, piv_t& piv)
{
    using idx_t = size_type<matrix_t>;
    using T = type_t<matrix_t>;
    using real_t = real_type<T>;
 
    // constants
    const idx_t m = nrows(A);
    const idx_t n = ncols(A);
    const idx_t end = min(m, n);
 
    // check arguments
    tlapack_check((idx_t)size(piv) >= end);
 
    // quick return
    if (m <= 0 || n <= 0) return 0;
 
    for (idx_t j = 0; j < end; j++) {
        // find pivot and swap the row with pivot row
        piv[j] = j;
        for (idx_t i = j + 1; i < m; i++) {
            if (abs1(A(i, j)) > abs1(A(piv[j], j))) piv[j] = i;
        }
 
        // if nonzero pivot does not exist, return
        if (A(piv[j], j) == real_t(0)) {
            return j + 1;
        }
 
        // if the pivot happens to be a piv[j]>j(piv[j] not equal to j), then
        // swap j-th row and piv[j] row of A
        if ((idx_t)piv[j] != j) {
            for (idx_t i = 0; i < n; i++) {
                T tmp = A(j, i);
                A(j, i) = A(piv[j], i);
                A(piv[j], i) = tmp;
            }
        }
 
        // divide below diagonal part of j-th column by the element on the
        // diagonal(A(j,j))
        // This is a level 0 routine, so we do not use tlapack::rscl.
        for (idx_t i = j + 1; i < m; i++) {
            A(i, j) /= A(j, j);
        }
 
        // update the submatrix A(j+1:m-1,j+1:n-1)
        for (idx_t row = j + 1; row < m; row++) {
            for (idx_t col = j + 1; col < n; col++) {
                A(row, col) -= A(row, j) * A(j, col);
            }
        }
    }
 
    return 0;
}
 
}  // namespace tlapack
 
#endif  // TLAPACK_GETRF_LV0_HH