core_clrdbg.c

Go to the documentation of this file.

/**
 *
 * @file core_clrdbg.c
 *
 * PaStiX low-rank kernel debug routines that may be call within gdb.
 *
 * @copyright 2016-2024 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
 *                      Univ. Bordeaux. All rights reserved.
 *
 * @version 6.4.0
 * @author Gregoire Pichon
 * @author Mathieu Faverge
 * @author Pierre Ramet
 * @date 2024-07-05
 * @generated from /builds/2mk6rsew/0/solverstack/pastix/kernels/core_zlrdbg.c, normal z -> c, Tue Feb 25 14:34:50 2025
 *
 **/
#include "common.h"
#include <cblas.h>
#include <lapacke.h>
#include "pastix_clrcores.h"
#include "c_nan_check.h"
 
/**
 *******************************************************************************
 *
 * @brief Print the svd values of the given matrix.
 *
 *******************************************************************************
 *
 * @param[in] M
 *          The number of rows of the matrix A.
 *
 * @param[in] N
 *          The number of columns of the matrix A.
 *
 * @param[in] A
 *          The matrix A to study of size lda-by-N
 *
 * @param[in] lda
 *          The leading dimension of the matrix A. lda = max( 1, M )
 *
 *******************************************************************************/
void
core_clrdbg_printsvd( pastix_int_t              M,
                      pastix_int_t              N,
                      const pastix_complex32_t *A,
                      pastix_int_t              lda )
{
    pastix_int_t i, ret;
    pastix_int_t minMN = pastix_imin( M, N );
    size_t lrwork =  2 * minMN;
    size_t lzwork =  M * N;
    pastix_complex32_t *W;
    float *s, *superb;
 
    W = malloc( lzwork * sizeof(pastix_complex32_t) + lrwork * sizeof(float) );
    s = (float*)(W + M*N);
    superb = s + minMN;
 
    ret = LAPACKE_clacpy_work( LAPACK_COL_MAJOR, 'A', M, N, A, lda, W, M );
    assert( ret == 0 );
    ret = LAPACKE_cgesvd(LAPACK_COL_MAJOR, 'N', 'N', M, N, W, M, s, NULL, 1, NULL, 1, superb );
    assert( ret == 0 );
 
    for(i=0; i<minMN; i++) {
        fprintf( stderr, "%e ", s[i] );
    }
    fprintf(stderr, "\n");
 
    (void)ret;
    free(W);
}
 
/**
 *******************************************************************************
 *
 * @brief Check the orthogonality of the matrix A
 *
 *******************************************************************************
 *
 * @param[in] M
 *          The number of rows of the matrix A.
 *
 * @param[in] N
 *          The number of columns of the matrix A.
 *
 * @param[in] A
 *          The matrix A to study of size lda-by-N
 *
 * @param[in] lda
 *          The leading dimension of the matrix A. lda = max( 1, M )
 *
 *******************************************************************************
 *
 * @retval 0 if the matrix A is orthogonal
 * @retval 1 if the matrix A is not orthogonal
 *
 *******************************************************************************/
int
core_clrdbg_check_orthogonality( pastix_int_t              M,
                                 pastix_int_t              N,
                                 const pastix_complex32_t *A,
                                 pastix_int_t              lda )
{
    pastix_complex32_t *Id;
    float alpha, beta;
    float normQ, res;
    pastix_int_t info_ortho, ret;
    pastix_int_t minMN = pastix_imin(M, N);
    pastix_int_t maxMN = pastix_imax(M, N);
    float eps = LAPACKE_slamch_work('e');
 
    alpha = 1.0;
    beta  = -1.0;
 
    /* Build the identity matrix */
    Id = malloc( minMN * minMN * sizeof(pastix_complex32_t) );
    ret = LAPACKE_claset_work( LAPACK_COL_MAJOR, 'A', minMN, minMN,
                              0., 1., Id, minMN );
    assert( ret == 0 );
 
    if (M > N) {
        /* Perform Id - Q'Q */
        cblas_cherk(CblasColMajor, CblasUpper, CblasConjTrans, N, M, alpha, A, lda, beta, Id, minMN);
    }
    else {
        /* Perform Id - QQ' */
        cblas_cherk(CblasColMajor, CblasUpper, CblasNoTrans,   M, N, alpha, A, lda, beta, Id, minMN);
    }
 
    normQ = LAPACKE_clanhe_work( LAPACK_COL_MAJOR, 'M', 'U', minMN, Id, minMN, NULL );
    res = normQ / (maxMN * eps);
 
    if ( isnan(res) || isinf(res) || (res > 60.0) ) {
        fprintf(stderr, "Check Orthogonality: || I - Q*Q' || = %e, ||Id-Q'*Q||_oo / (N*eps) = %e : \n",
                normQ, res );
        info_ortho = 1;
    }
    else {
        info_ortho = 0;
    }
 
    free(Id);
    (void)ret;
    return info_ortho;
}
 
/**
 *******************************************************************************
 *
 * @brief Check the orthogonality of the matrix A relatively to the matrix B
 *
 * Check that A^t B = 0
 *
 *******************************************************************************
 *
 * @param[in] M
 *          The number of rows of the matrix A.
 *
 * @param[in] NA
 *          The number of columns of the matrix A.
 *
 * @param[in] NB
 *          The number of columns of the matrix B.
 *
 * @param[in] A
 *          The matrix A to study of size lda-by-NA
 *
 * @param[in] lda
 *          The leading dimension of the matrix A. lda = max( 1, M )
 *
 * @param[in] B
 *          The matrix B to study of size ldb-by-NB
 *
 * @param[in] ldb
 *          The leading dimension of the matrix B. ldb = max( 1, M )
 *
 *******************************************************************************
 *
 * @retval 0 if the matrices A and B are orthogonal
 * @retval 1 if the matrices A anb B are not orthogonal
 *
 *******************************************************************************/
int
core_clrdbg_check_orthogonality_AB( pastix_int_t              M,
                                    pastix_int_t              NA,
                                    pastix_int_t              NB,
                                    const pastix_complex32_t *A,
                                    pastix_int_t              lda,
                                    const pastix_complex32_t *B,
                                    pastix_int_t              ldb )
{
    pastix_complex32_t *Zero;
    float norm, res;
    pastix_int_t info_ortho, ret;
    float eps = LAPACKE_slamch_work('e');
    pastix_complex32_t cone = 1.0;
    pastix_complex32_t czero = 0.0;
 
    /* Build the null matrix */
    Zero = malloc( NA * NB * sizeof(pastix_complex32_t) );
    ret = LAPACKE_claset_work( LAPACK_COL_MAJOR, 'A', NA, NB,
                               0., 0., Zero, NA );
    assert( ret == 0 );
 
    cblas_cgemm(CblasColMajor, CblasConjTrans, CblasNoTrans,
                NA, NB, M,
                CBLAS_SADDR(cone),  A, lda,
                                    B, ldb,
                CBLAS_SADDR(czero), Zero, NA);
 
    norm = LAPACKE_clange_work( LAPACK_COL_MAJOR, 'M', NA, NB, Zero, NA, NULL );
    res = norm / (M * eps);
 
    if ( isnan(res) || isinf(res) || (res > 60.0) ) {
        fprintf(stderr, "Check Orthogonality: || A' B || = %e, || A' B ||_oo / (M*eps) = %e : \n",
                norm, res );
        info_ortho = 1;
    }
    else {
        info_ortho = 0;
    }
 
    free(Zero);
    (void)ret;
    return info_ortho;
}