//------------------------------------------------------------------------------
// gb_norm: compute the norm of a GraphBLAS matrix
//------------------------------------------------------------------------------

// SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2026, All Rights Reserved.
// SPDX-License-Identifier: Apache-2.0

//------------------------------------------------------------------------------

#undef  FREE_ALL
#define FREE_ALL            \
    GrB_Matrix_free (&X) ;  \
    GrB_Vector_free (&t) ;

GrB_Info gb_norm            // compute norm (A,kind)
(
    // output:
    double *s,              // norm of A
    // inputs:
    GrB_Matrix A,
    int64_t norm_kind,      // 0, 1, 2, INT64_MAX, or INT64_MIN
    const int arena,
    char err [ERRLEN]
)
{

    //--------------------------------------------------------------------------
    // get input matrix, select types and operators, and allocate X
    //--------------------------------------------------------------------------

    GrB_Matrix X = NULL ;
    GrB_Vector t = NULL ;
    (*s) = 0 ;

    uint64_t nrows, ncols, nvals ;
    OK (GrB_Matrix_nvals (&nvals, A)) ;
    if (nvals == 0)
    { 
        return (GrB_SUCCESS) ;
    }

    GrB_Type atype, xtype ;
    OK (GrB_Matrix_nrows (&nrows, A)) ;
    OK (GrB_Matrix_ncols (&ncols, A)) ;
    OK (GxB_Matrix_type (&atype, A)) ;

    GrB_UnaryOp absop ;
    GrB_Monoid sumop, maxop, minop ;
    bool is_complex = false ;
    bool is_dense = false ;

    if (atype == GrB_FP32)
    { 
        // if A is FP32, use the FP32 type and operators
        xtype = GrB_FP32 ;
        absop = GrB_ABS_FP32 ;
        sumop = GrB_PLUS_MONOID_FP32 ;
        maxop = GrB_MAX_MONOID_FP32 ;
        minop = GrB_MIN_MONOID_FP32 ;
    }
    else if (atype == GxB_FC32)
    { 
        // if A is FC32, use the FP32/FC32 type and operators
        is_complex = true ;
        xtype = GrB_FP32 ;
        absop = GxB_ABS_FC32 ;
        sumop = GrB_PLUS_MONOID_FP32 ;
        maxop = GrB_MAX_MONOID_FP32 ;
        minop = GrB_MIN_MONOID_FP32 ;
    }
    else if (atype == GxB_FC64)
    { 
        // if A is FC64, use the FP64/FC64 type and operators
        is_complex = true ;
        xtype = GrB_FP64 ;
        absop = GxB_ABS_FC64 ;
        sumop = GrB_PLUS_MONOID_FP64 ;
        maxop = GrB_MAX_MONOID_FP64 ;
        minop = GrB_MIN_MONOID_FP64 ;
    }
    else
    { 
        // otherwise, use FP64 type and operators; this will typecast the 
        // input matrix to FP64 if A is not of that type.
        xtype = GrB_FP64 ;
        absop = GrB_ABS_FP64 ;
        sumop = GrB_PLUS_MONOID_FP64 ;
        maxop = GrB_MAX_MONOID_FP64 ;
        minop = GrB_MIN_MONOID_FP64 ;
    }

    OK (GxB_Matrix_new_arena (&X, xtype, nrows, ncols, arena, arena)) ;

    //--------------------------------------------------------------------------
    // compute the norm
    //--------------------------------------------------------------------------

    if (nrows == 1 || ncols == 1 || norm_kind == 0)
    {

        //----------------------------------------------------------------------
        // vector or Frobenius norm
        //----------------------------------------------------------------------

        switch (norm_kind)
        {

            case 0 :    // Frobenius norm
            case 2 :    // 2-norm

                if (is_complex)
                { 
                    // X = abs (A)
                    OK1 (X, GrB_Matrix_apply (X, NULL, NULL, absop, A, NULL)) ;
                    // X = X.^2
                    if (atype == GxB_FC32)
                    { 
                        OK1 (X, GrB_Matrix_apply_BinaryOp2nd_FP32 (X, NULL,
                            NULL, GxB_POW_FP32, X, (float) 2.0, NULL)) ;
                    }
                    else
                    { 
                        OK1 (X, GrB_Matrix_apply_BinaryOp2nd_FP64 (X, NULL,
                            NULL, GxB_POW_FP64, X, (double) 2.0, NULL)) ;
                    }
                }
                else
                { 
                    // X = A.^2
                    OK1 (X, GrB_Matrix_apply_BinaryOp2nd_FP64 (X, NULL, NULL,
                        GxB_POW_FP64, A, (double) 2.0, NULL)) ;
                }
                // s = sum (X)
                OK (GrB_Matrix_reduce_FP64 (s, NULL, sumop, X, NULL)) ;
                (*s) = sqrt ((*s)) ;
                break ;

            case 1 :    // 1-norm

                // X = abs (A)
                OK1 (X, GrB_Matrix_apply (X, NULL, NULL, absop, A, NULL)) ;
                // s = sum (X)
                OK (GrB_Matrix_reduce_FP64 (s, NULL, sumop, X, NULL)) ;
                break ;

            case INT64_MAX :    // inf-norm

                // X = abs (A)
                OK1 (X, GrB_Matrix_apply (X, NULL, NULL, absop, A, NULL)) ;
                // s = max (X)
                OK (GrB_Matrix_reduce_FP64 (s, NULL, maxop, X, NULL)) ;
                break ;

            case INT64_MIN :    // (-inf)-norm

                OK (gb_is_dense (&is_dense, A, err)) ;
                if (is_dense)
                { 
                    // X = abs (A)
                    OK1 (X, GrB_Matrix_apply (X, NULL, NULL, absop, A, NULL)) ;
                    // s = min (X)
                    OK (GrB_Matrix_reduce_FP64 (s, NULL, minop, X, NULL)) ;
                }
                break ;

            default:

                ERROR ("unknown norm", GrB_INVALID_VALUE) ;
                break ;
        }

    }
    else
    {

        //----------------------------------------------------------------------
        // matrix norm
        //----------------------------------------------------------------------

        switch (norm_kind)
        {

            case 2 :    // 2-norm

                ERROR ("2-norm not available for GrB matrices",
                    GrB_NOT_IMPLEMENTED) ;
                break ;

            case 1 :    // 1-norm:  max sum of columns of abs (A)

                // X = abs (A)
                OK1 (X, GrB_Matrix_apply (X, NULL, NULL, absop, A, NULL)) ;
                // t = zeros (ncols,1)
                OK (GxB_Vector_new_arena (&t, xtype, ncols, arena, arena)) ;
                // t(j) = sum of the ith column, X(:,j)
                OK (GrB_Matrix_reduce_Monoid (t, NULL, NULL, sumop, X,
                    GrB_DESC_T0)) ;
                // s = max (t)
                OK (GrB_Vector_reduce_FP64 (s, NULL, maxop, t, NULL)) ;
                break ;

            case INT64_MAX :    // inf-norm:  max sum of rows of abs (A)

                // X = abs (A)
                OK1 (X, GrB_Matrix_apply (X, NULL, NULL, absop, A, NULL)) ;
                // t = zeros (nrows,1)
                OK (GxB_Vector_new_arena (&t, xtype, nrows, arena, arena)) ;
                // t(i) = sum of the ith row, X(i,:)
                OK (GrB_Matrix_reduce_Monoid (t, NULL, NULL, sumop, X, NULL)) ;
                // s = max (t)
                OK (GrB_Vector_reduce_FP64 (s, NULL, maxop, t, NULL)) ;
                break ;

            case INT64_MIN :

                ERROR ("(-inf)-norm not available for GrB matrices",
                    GrB_NOT_IMPLEMENTED) ;
                break ;

            default :

                ERROR ("unknown norm", GrB_INVALID_VALUE) ;
                break ;
        }
    }

    //--------------------------------------------------------------------------
    // free workspace and return result
    //--------------------------------------------------------------------------

    FREE_ALL ;
    return (GrB_SUCCESS) ;
}

#undef  FREE_ALL
#define FREE_ALL

