dggev.F

Go to the documentation of this file.

#include "fintrf.h"
C======================================================================
C     Gateway routine to call the zggev.f function from LAPACK package.
C     (See: http://www.netlib.org/lapack/explore-html/d9/d52/dggev_8f.html for details)
C     Copyright (C) 2009-2015 Peter Rakyta, Ph.D.
C
C     This program is free software: you can redistribute it and/or modify
C     it under the terms of the GNU General Public License as published by
C     the Free Software Foundation, either version 3 of the License, or
C     (at your option) any later version.
C 
C     This program is distributed in the hope that it will be useful,
C     but WITHOUT ANY WARRANTY; without even the implied warranty of
C     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C     GNU General Public License for more details.
C
C     You should have received a copy of the GNU General Public License
C     along with this program.  If not, see http://www.gnu.org/licenses/.
C     
C======================================================================
C     Gateway routine
      subroutine mexfunction(nlhs, plhs, nrhs, prhs)

C     Declarations
      implicit none

C     mexFunction arguments:
      mwpointer plhs(*), prhs(*)
      integer nlhs, nrhs


C     Function declarations:
      mwpointer mxgetpr
      mwpointer mxgetpi
      mwpointer mxcreatedoublematrix
      mwpointer mxgetm, mxgetn
      integer*4 mxIsNumeric
      integer*4 mxIsComplex

      

C     Pointers to input/output mxArrays:
      mwpointer alphar_ptr
      mwpointer alphai_ptr
      mwpointer beta_ptr
      mwpointer lvecs_ptr
      mwpointer rvecs_ptr

C     Array information:
      mwsize hrows, hcols, brows, bcols, cols
      mwsize sizeh, sizeb

C     Arguments for computational routine:      
      real*8, allocatable :: h_eff(:,:), b_eff(:,:)
      real*8, allocatable :: alphar(:), alphai(:)
      real*8, allocatable :: beta(:)
      real*8, allocatable :: l_eigvec(:,:)
      real*8, allocatable :: r_eigvec(:,:)     

      real*8 work1(1)
      real*8, allocatable :: work(:)
      
      mwsignedindex lwork, info


C-----------------------------------------------------------------------
C     Check for proper number of arguments. 
      if(nrhs .ne. 2) then
         call mexerrmsgidandtxt ('MATLAB:timestwo:nInput',
     +                           'Two inputs required.')
      elseif(nlhs .gt. 5) then
         call mexerrmsgidandtxt ('MATLAB:timestwo:nOutput',
     +                           'Too many output arguments.')
      endif

C     Validate inputs
      
C     Check that the input is a number.
      if ((mxisnumeric(prhs(1)) .eq. 0)  .or.           
     +        (mxisnumeric(prhs(2)) .ne. 1)) then
         call mexerrmsgidandtxt ('MATLAB:timestwo:NonNumeric',
     +                           'Input must be a number.')
      endif
      
C     Check if inputs are complex.
      if ((mxiscomplex(prhs(1)) .ne. 0) .or.           
     +        (mxiscomplex(prhs(2)) .ne. 0)) then
         call mexerrmsgidandtxt ('MATLAB:convec:NonComplex',
     +                           'Inputs must be real.')
      endif    

C     Get the size of the input array. (MATLAB)
      hrows = mxgetm(prhs(1))
      hcols = mxgetn(prhs(1))
      brows = mxgetm(prhs(2))
      bcols = mxgetn(prhs(2))
      sizeh = hrows*hcols
      sizeb = brows*bcols
      
      
      if ((hrows .ne. brows) .or.           
     +        (hcols .ne. bcols)) then
         call mexerrmsgidandtxt ('MATLAB:zggev:WrongSize',
     +       'The dimensions of the matrices must equal.')
      endif            
                    
     
C     Allocating memory for the matrices
      allocate(h_eff(hrows,hcols))
      allocate(b_eff(brows,bcols))
      allocate(alphar(brows))
      allocate(alphai(brows))
      allocate(beta(brows))
      allocate(l_eigvec(brows,bcols))
      allocate(r_eigvec(brows,bcols))
C
C     Check the allocated arrays
      if (.not. ( allocated(h_eff) .and.
     + allocated(h_eff) .and.
     + allocated(b_eff) .and.
     + allocated(alphar) .and.
     + allocated(alphai) .and.
     + allocated(beta) .and.
     + allocated(l_eigvec) .and.
     + allocated(r_eigvec) )) then
         call mexerrmsgidandtxt ('MATLAB:zggev:BadAllocation',
     +       'The allocation of variables was unsuccesfull.')
      endif
C
C     Copy input data to the Fortran arrays
      call mxcopyptrtoreal8(mxgetpr(prhs(1)),
     +                          h_eff,sizeh)
     
      call mxcopyptrtoreal8(mxgetpr(prhs(2)),
     +                          b_eff,sizeb)     

      lwork = -1
C     Query for the workspace
      call dggev( 'V', 'V', hcols, h_eff, hrows,
     + b_eff, brows, 
     + alphar, alphai, beta, l_eigvec, hrows, r_eigvec, hrows, 
     + work1, lwork, info )
C   
C     Test the resulted INFO flag
      call checkinfo( info )
C     
      lwork = work1(1)
      allocate( work(lwork))
      if (.not. ( allocated(work) )) then
         call mexerrmsgidandtxt ('MATLAB:zggev:BadAllocation',
     +       'The allocation of variable Work was unsuccesfull.')
      endif
C
C     Call the computational subroutine.     
      call dggev( 'V', 'V', hcols, h_eff, hrows,
     + b_eff, hrows, 
     + alphar, alphai, beta, l_eigvec, hrows, r_eigvec, hrows, 
     + work, lwork, info )
C   
C     Test the resulted INFO flag
      call checkinfo( info )
     
C     Create matrices for the return arguments.
      cols = 1
      plhs(1) = mxcreatedoublematrix(hrows,cols,0) ! ALPHAR
      plhs(2) = mxcreatedoublematrix(hrows,cols,0) ! ALPHAI
      plhs(3) = mxcreatedoublematrix(hrows,cols,0) ! BETA
      plhs(4) = mxcreatedoublematrix(hrows,hcols,1) ! l_eigvec
      plhs(5) = mxcreatedoublematrix(hrows,hcols,1) ! r_eigvec
      
C     Load the data into to MATLAB outputs.      
      alphar_ptr = mxgetpr(plhs(1))
      call mxcopyreal8toptr(alphar,alphar_ptr, hrows)
C
      alphai_ptr = mxgetpr(plhs(2))
      call mxcopyreal8toptr(alphai,alphai_ptr, hrows) 
C      
      beta_ptr = mxgetpr(plhs(3)) 
      call mxcopyreal8toptr(beta,beta_ptr, hrows)  
C      
      lvecs_ptr = mxgetpr(plhs(4))     
      call mxcopyreal8toptr(l_eigvec, lvecs_ptr,sizeb) 
C      
      rvecs_ptr = mxgetpr(plhs(5))
      call mxcopyreal8toptr(r_eigvec, rvecs_ptr,sizeb)
  
     

C     Free the allocated matrices
      deallocate(h_eff)
      deallocate(b_eff)
      deallocate(alphar)
      deallocate(alphai)
      deallocate(beta)
      deallocate(l_eigvec)
      deallocate(r_eigvec)
      deallocate(work)


      return
      end



C=========================================================
C      Test the output INFO flag of the LAPACK function
      subroutine checkinfo( INFO )

      mwsignedindex info
      character(len=100) :: msg
      
      if (info<0) then
        msg = 'Illegal argumemnt in the LAPACK function'
        call mexerrmsgidandtxt ('MATLAB:zggev:IllegalLAPACKArgument',
     +       msg)
      elseif (info>0) then
        msg = 'Calculations failed in the LAPACK function'
        call mexerrmsgidandtxt ('MATLAB:zggev:FailedLAPACK',
     +       msg)
      endif

      return
      end subroutine