! Interface to some Lapack functions
! ==================================
!
! packs the dgees and zgees Lapack subroutines into a single, overloaded
! wrapper CPL function eigv(A) providing an ARRAY of eigenvalues
! of a general real or complex matrix as its result.
! This example is a copy of the lapack-eigenvalues library provided in the
! main package. It may serve as a template to include more Lapack functions.

USE complex

#link "-llapack -lblas"
! directive to the linker channeled through the C compiler.
! On Debian the liblapack-dev package or one equivalent must have been installed

INLINE ARRAY(A.LO..A.HI) OF COMPLEX FUNCTION eigv(COMPLEX A(*,*))
  ! This function will produce a constant array if assigned to an undeclared id
  INTEGER SDIM,info
  LA=LENGTH(A)
  IF LA#LENGTH(A(LO)) THEN ERROR "zgees: matrix is not square"
  IF STRIDEOF[A(LO)]#1 THEN ERROR "zgees: matrix rows are not contiguous"
  ! A CPL matrix is allowed to be not contiguous, but a FORTRAN matrix must be.
  LWORK=3*LA
  COMPLEX WORK(1..LWORK)
  REAL RWORK(1..LA)
  FORTRANCALL zgees("N","N",NULL,LA,A,STRIDEOF(A),SDIM,RESULT,NULL,LA,WORK,LWORK,RWORK,NULL,info)
  ! interface to FORTRAN call, applied without any prototype checking.
  ! However, a prototype must be included if this file is to be USEd by icpl;
  ! see lapack-eigenvalues.cpl for an example.
  IF info#0 THEN ERROR "zgees: " info
END eigv   ! encapsulating zgees inside eigv turns on CPL prototype checking.

INLINE ARRAY(A.LO..A.HI) OF COMPLEX FUNCTION eigv(REAL A(*,*))
  ! same function name as above: correct instance of overloaded function
  ! will be selected based on argument type
  INTEGER SDIM,info
  LA=LENGTH(A)
  IF LA#LENGTH(A(LO)) THEN ERROR "dgees: matrix is not square"
  IF STRIDEOF[A(LO)]#1 THEN ERROR "dgees: matrix rows are not contiguous" 
  ! A CPL matrix is allowed to be not contiguous, but a FORTRAN matrix must be.
  LWORK=3*LA
  REAL WORK(1..LWORK)
  REAL WR(1..LA),WI(1..LA)
  FORTRANCALL dgees("N","N",NULL,LA,A,STRIDEOF(A),SDIM,WR,WI,NULL,LA,WORK,LWORK,NULL,info)
  ! interface to FORTRAN call, applied without any prototype checking.
  IF info#0 THEN ERROR "dgees: " info
  RESULT(*+LO-1).REAL=WR; RESULT(*+LO-1).IMAG=WI
END eigv   ! encapsulating dgees inside eigv turns on CPL prototype checking.

! Example: uncomment the text below, compile this file and execute.
! Alternately type similar commands in interactive cpl, but USE lapack-eigenvalues instead.
!(
  A=[(1.,2.),(3.,4.)]
  WRITE eigv(A)
!)