! 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) !)