! 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 icpl, but USE lapack-eigenvalues instead.
!(
A=[(1.,2.),(3.,4.)]
WRITE eigv(A)
!)