Skip to content

Commit

Permalink
Added reference implementation of Hermitian EVD.
Browse files Browse the repository at this point in the history
Details:
- Defined a new level-4 operation, hevd, which performs a Hermitian
  eigenvalue decomposition (EVD) via the QR algorithm (i.e., Givens
  rotations). The operation has two companion functions:
  - rhevd: a so-called "reverse" Hermitian EVD in which the eigenvector
    and eigenvalue output from an hEVD is used to recompute the original
    Hermitian matrix.
  - hevpinv: an operation which explicitly computes the pseudo-inverse
    of a Hermitian matrix via hevd, inverttv, and rhevd, where inverttv
    inverts each element of a vector, unless the absolute value of the
    element is less than a given threshold, in which case the value is
    set to 0.
  The implementation of hevd is, for now, based on code that was
  converted via f2c. (At some point in the future, the implementation
  will likely be replaced (or at least augmented) by a native
  implementation based on the so-called Restructured QR algorithm.)
  This f2c'ed code lives in frame/4/f2c and is based on netlib LAPACK
  3.9.1.
- Added a testsuite module to test the new hevd operation. It is tested
  by default for all four input.operations.* files.
- Added scripts, called via a Makefile, that gradually process the
  output of f2c into something that can compile and link in the BLIS
  build system (without reliance on libf2c). These scripts (and
  Makefile) live in frame/4/other/staging and the original Fortran code
  is in frame/4/other/staging/netlib.
- Implemented fscanm and fscanv operations, which will read in a matrix
  or vector, respectively, from a C file stream, as well as their
  stdin-assumping counterparts, scanm and scanv. For now, these
  operations live in frame/base/scan.
- Implemented inverttv, which inverts each element of a vector unless
  the absolute value of the element in question is less than a given
  threshold, in which case the value is set to 0 (instead of being
  inverted).
- Defined dlamch_() and slamch_() in frame/compat/bla_lamch.c, which
  are implemented as wrappers to their pre-existing bli_?lamch()
  counterparts.
- Changed various BLAS code in frame/compat to call bla_lsame_()
  instead of lsame_(). Also relocated frame/compat/f2c/bla_lsame.c to
  'util' subdir.
- Defined f2c_[cz]dot[cu]() functions, which serve as f2c-compatible
  wrappers to the corresponding BLAS functions. (By f2c-compatible, we
  mean that they "return" their result in a hidden extra parameter that
  is emitted by f2c.)
- Added various missing operations to frame/compate/f2c/util, such as:
  bla_r_real(), bla_d_real(), bla_i_len(), bla_nint(), bla_pow_??(),
  bla_s_cmp(), and bla_s_copy(). (In addition, bla_a_max() and
  bla_a_min() were implemented as macros.) These functions serve to
  provide siloed libf2c functionality in a BLIS-specific namespace.
- Modified existing bla_r_abs(), bla_d_abs(), and bla_r_imag(), to use
  call-by-value conventions. Also return double in all cases, even if
  input is a float (bla_real).
- Defined PASTEF2C macros that works like PASTEF77 except that it also
  prepends f2c_ to the final token.
- Whitespace changes.
  • Loading branch information
fgvanzee committed Jul 14, 2023
1 parent 1232831 commit 59e6151
Show file tree
Hide file tree
Showing 338 changed files with 68,981 additions and 325 deletions.
6 changes: 6 additions & 0 deletions frame/4/bli_l4.h
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,8 @@

#include "bli_l4_check.h"

#include "f2c_lapack.h"

// Define function types.
#include "bli_l4_ft_opt.h"
#include "bli_l4_ft.h"
Expand All @@ -55,3 +57,7 @@
#include "bli_trinv.h"
#include "bli_ttmm.h"
#include "bli_hpdinv.h"

#include "bli_hevd.h"
#include "bli_hevpinv.h"

153 changes: 153 additions & 0 deletions frame/4/bli_l4_check.c
Original file line number Diff line number Diff line change
Expand Up @@ -186,4 +186,157 @@ void bli_hpdinv_check
bli_check_error_code( e_val );
}

void bli_hevd_check
(
const obj_t* a,
const obj_t* v,
const obj_t* e,
const cntx_t* cntx
)
{
err_t e_val;

// Check object datatypes.

e_val = bli_check_floating_object( a );
bli_check_error_code( e_val );

e_val = bli_check_floating_object( v );
bli_check_error_code( e_val );

e_val = bli_check_floating_object( e );
bli_check_error_code( e_val );

e_val = bli_check_real_object( e );
bli_check_error_code( e_val );

e_val = bli_check_consistent_object_datatypes( a, v );
bli_check_error_code( e_val );

e_val = bli_check_object_real_proj_of( v, e );
bli_check_error_code( e_val );

// Check object dimensions.

e_val = bli_check_matrix_object( a );
bli_check_error_code( e_val );

e_val = bli_check_matrix_object( v );
bli_check_error_code( e_val );

e_val = bli_check_vector_object( e );
bli_check_error_code( e_val );

// Check matrix squareness.

e_val = bli_check_square_object( a );
bli_check_error_code( e_val );

e_val = bli_check_square_object( v );
bli_check_error_code( e_val );

// Check object dimensions.

e_val = bli_check_conformal_dims( a, v );
bli_check_error_code( e_val );

e_val = bli_check_vector_dim_equals( e, bli_obj_length( v ) );
bli_check_error_code( e_val );

// Check object structure.

bool is_herm = bli_obj_is_hermitian( a );
bool is_real_symm = bli_obj_is_symmetric( a ) && bli_obj_is_real( a );

if ( !is_herm && !is_real_symm )
{
e_val = BLIS_EXPECTED_HERMITIAN_OBJECT;
bli_check_error_code( e_val );
}

// Check object buffers (for non-NULLness).

e_val = bli_check_object_buffer( a );
bli_check_error_code( e_val );

e_val = bli_check_object_buffer( v );
bli_check_error_code( e_val );

e_val = bli_check_object_buffer( e );
bli_check_error_code( e_val );
}

void bli_rhevd_check
(
const obj_t* v,
const obj_t* e,
const obj_t* a,
const cntx_t* cntx
)
{
bli_hevd_check( a, v, e, cntx );
}

void bli_hevpinv_check
(
double thresh,
const obj_t* a,
const obj_t* p,
const cntx_t* cntx
)
{
err_t e_val;

// Check object datatypes.

e_val = bli_check_floating_object( a );
bli_check_error_code( e_val );

e_val = bli_check_floating_object( p );
bli_check_error_code( e_val );

e_val = bli_check_consistent_object_datatypes( a, p );
bli_check_error_code( e_val );

// Check object dimensions.

e_val = bli_check_matrix_object( a );
bli_check_error_code( e_val );

e_val = bli_check_matrix_object( p );
bli_check_error_code( e_val );

// Check matrix squareness.

e_val = bli_check_square_object( a );
bli_check_error_code( e_val );

e_val = bli_check_square_object( p );
bli_check_error_code( e_val );

// Check object dimensions.

e_val = bli_check_conformal_dims( a, p );
bli_check_error_code( e_val );

// Check object structure.

bool is_herm = bli_obj_is_hermitian( a );
bool is_real_symm = bli_obj_is_symmetric( a ) && bli_obj_is_real( a );

if ( !is_herm && !is_real_symm )
{
e_val = BLIS_EXPECTED_HERMITIAN_OBJECT;
bli_check_error_code( e_val );
}

// Check object buffers (for non-NULLness).

e_val = bli_check_object_buffer( a );
bli_check_error_code( e_val );

e_val = bli_check_object_buffer( p );
bli_check_error_code( e_val );
}

#endif
24 changes: 24 additions & 0 deletions frame/4/bli_l4_check.h
Original file line number Diff line number Diff line change
Expand Up @@ -64,4 +64,28 @@ void bli_hpdinv_check
const cntx_t* cntx
);

void bli_hevd_check
(
const obj_t* a,
const obj_t* v,
const obj_t* e,
const cntx_t* cntx
);

void bli_rhevd_check
(
const obj_t* v,
const obj_t* e,
const obj_t* a,
const cntx_t* cntx
);

void bli_hevpinv_check
(
double thresh,
const obj_t* a,
const obj_t* p,
const cntx_t* cntx
);

#endif
25 changes: 25 additions & 0 deletions frame/4/bli_l4_fpa.c
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,31 @@

#ifdef BLIS_ENABLE_LEVEL4

//
// Define function pointer query interfaces for typed APIs.
//

#undef GENFRONT
#define GENFRONT( opname ) \
\
GENARRAY_FPA( PASTECH2(opname,BLIS_TAPI_EX_SUF,_vft), \
PASTECH(opname,BLIS_TAPI_EX_SUF) ); \
\
PASTECH2(opname,BLIS_TAPI_EX_SUF,_vft) \
PASTEMAC2(opname,BLIS_TAPI_EX_SUF,_qfp)( num_t dt ) \
{ \
return PASTECH2(opname,BLIS_TAPI_EX_SUF,_fpa)[ dt ]; \
}


//GENFRONT( chol )
//GENFRONT( trinv )
//GENFRONT( ttmm )
GENFRONT( hevd )
GENFRONT( rhevd )
GENFRONT( hevpinv )


//
// Define function pointer query interfaces for implementations.
//
Expand Down
18 changes: 18 additions & 0 deletions frame/4/bli_l4_fpa.h
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,24 @@
#ifndef BLIS_L4_FPA_H
#define BLIS_L4_FPA_H

//
// Prototype function pointer query interfaces for typed APIs.
//

#undef GENPROT
#define GENPROT( opname ) \
\
PASTECH2(opname,BLIS_TAPI_EX_SUF,_vft) \
PASTEMAC2(opname,BLIS_TAPI_EX_SUF,_qfp)( num_t dt );

//GENPROT( chol )
//GENPROT( trinv )
//GENPROT( ttmm )
GENPROT( hevd )
GENPROT( rhevd )
GENPROT( hevpinv )


//
// Prototype function pointer query interfaces for implementations.
//
Expand Down
60 changes: 60 additions & 0 deletions frame/4/bli_l4_ft.h
Original file line number Diff line number Diff line change
Expand Up @@ -75,4 +75,64 @@ typedef err_t (*PASTECH2(ch,opname,tsuf)) \

INSERT_GENTDEF( trinv )

// hevd

#undef GENTDEFR
#define GENTDEFR( ctype, ctype_r, ch, chr, opname, tsuf ) \
\
typedef err_t (*PASTECH3(ch,opname,BLIS_TAPI_EX_SUF,tsuf)) \
( \
bool comp_evecs, \
uplo_t uploa, \
dim_t m, \
const ctype* a, inc_t rs_a, inc_t cs_a, \
ctype* v, inc_t rs_v, inc_t cs_v, \
ctype_r* e, inc_t ince, \
ctype* work, \
dim_t lwork, \
ctype_r* rwork, \
const cntx_t* cntx, \
const rntm_t* rntm \
);

INSERT_GENTDEFR( hevd )

// rhevd

#undef GENTDEFR
#define GENTDEFR( ctype, ctype_r, ch, chr, opname, tsuf ) \
\
typedef err_t (*PASTECH3(ch,opname,BLIS_TAPI_EX_SUF,tsuf)) \
( \
uplo_t uploa, \
dim_t m, \
const ctype* v, inc_t rs_v, inc_t cs_v, \
ctype_r* e, inc_t ince, \
ctype* a, inc_t rs_a, inc_t cs_a, \
const cntx_t* cntx, \
const rntm_t* rntm \
);

INSERT_GENTDEFR( rhevd )

// hevpinv

#undef GENTDEFR
#define GENTDEFR( ctype, ctype_r, ch, chr, opname, tsuf ) \
\
typedef err_t (*PASTECH3(ch,opname,BLIS_TAPI_EX_SUF,tsuf)) \
( \
double thresh, \
uplo_t uploa, \
dim_t m, \
const ctype* a, inc_t rs_a, inc_t cs_a, \
ctype* p, inc_t rs_p, inc_t cs_p, \
const cntx_t* cntx, \
const rntm_t* rntm \
);

INSERT_GENTDEFR( hevpinv )


#endif

Loading

0 comments on commit 59e6151

Please sign in to comment.