potrf2.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_POTRF2_HH
#define TLAPACK_POTRF2_HH
 
#include "tlapack/base/utils.hpp"
#include "tlapack/blas/herk.hpp"
#include "tlapack/blas/trsm.hpp"
 
namespace tlapack {
 
template <TLAPACK_UPLO uplo_t, TLAPACK_SMATRIX matrix_t>
int potrf2(uplo_t uplo, matrix_t& A, const EcOpts& opts = {})
{
    using T = type_t<matrix_t>;
    using real_t = real_type<T>;
    using idx_t = size_type<matrix_t>;
    using range = pair<idx_t, idx_t>;
 
    // Constants
    const real_t one(1);
    const real_t zero(0);
    const idx_t n = nrows(A);
 
    // check arguments
    tlapack_check_false(uplo != Uplo::Lower && uplo != Uplo::Upper);
    tlapack_check_false(nrows(A) != ncols(A));
 
    // Quick return
    if (n <= 0) return 0;
 
    // Stop recursion
    else if (n == 1) {
        const real_t a00 = real(A(0, 0));
        if (a00 > zero) {
            A(0, 0) = sqrt(a00);
            return 0;
        }
        else {
            tlapack_error_if(
                opts.ec.internal, 1,
                "The leading minor of order 1 is not positive definite,"
                " and the factorization could not be completed.");
            return 1;
        }
    }
 
    // Recursive code
    else {
        const idx_t n1 = n / 2;
 
        // Define A11 and A22
        auto A11 = slice(A, range{0, n1}, range{0, n1});
        auto A22 = slice(A, range{n1, n}, range{n1, n});
 
        // Factor A11
        int info = potrf2(uplo, A11, NO_ERROR_CHECK);
        if (info != 0) {
            tlapack_error_if(
                opts.ec.internal, info,
                "The leading minor of the reported order is not positive "
                "definite,"
                " and the factorization could not be completed.");
            return info;
        }
 
        if (uplo == Uplo::Upper) {
            // Update and scale A12
            auto A12 = slice(A, range{0, n1}, range{n1, n});
            trsm(LEFT_SIDE, Uplo::Upper, CONJ_TRANS, NON_UNIT_DIAG, one, A11,
                 A12);
 
            // Update A22
            herk(UPPER_TRIANGLE, CONJ_TRANS, -one, A12, one, A22);
        }
        else {
            // Update and scale A21
            auto A21 = slice(A, range{n1, n}, range{0, n1});
            trsm(RIGHT_SIDE, Uplo::Lower, CONJ_TRANS, NON_UNIT_DIAG, one, A11,
                 A21);
 
            // Update A22
            herk(LOWER_TRIANGLE, NO_TRANS, -one, A21, one, A22);
        }
 
        // Factor A22
        info = potrf2(uplo, A22, NO_ERROR_CHECK);
        if (info == 0)
            return 0;
        else {
            tlapack_error_if(
                opts.ec.internal, info + n1,
                "The leading minor of the reported order is not positive "
                "definite,"
                " and the factorization could not be completed.");
            return info + n1;
        }
    }
}
 
}  // namespace tlapack
 
#endif  // TLAPACK_POTRF2_HH