zggev.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/d3/d47/zggev_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 alpha_ptr_real, alpha_ptr_imag
      mwpointer beta_ptr_real, beta_ptr_imag
      mwpointer lvecs_ptr_real, lvecs_ptr_imag
      mwpointer rvecs_ptr_real, rvecs_ptr_imag 

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

C     Arguments for computational routine:      
      complex*16, allocatable :: H_eff(:,:), B_eff(:,:)
      complex*16, allocatable :: ALPHA(:), BETA(:)
      complex*16, allocatable :: l_eigvec(:,:)
      complex*16, allocatable :: r_eigvec(:,:)     

      complex*16 Work1(1)
      complex*16, allocatable :: Work(:)
      DOUBLE PRECISION, allocatable :: RWork(:)
      
      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. 4) 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. 1) .or.           
     +        (mxiscomplex(prhs(2)) .ne. 1)) then
         call mexerrmsgidandtxt ('MATLAB:convec:NonComplex',
     +                           'Inputs must be complex.')
      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(alpha(brows))
      allocate(beta(brows))
      allocate(l_eigvec(brows,bcols))
      allocate(r_eigvec(brows,bcols))
      allocate(rwork(8*hrows))
C
C     Check the allocated arrays
      if (.not. ( allocated(h_eff) .and.
     + allocated(h_eff) .and.
     + allocated(b_eff) .and.
     + allocated(alpha) .and.
     + allocated(beta) .and.
     + allocated(l_eigvec) .and.
     + allocated(r_eigvec) .and.
     + allocated(rwork) )) then
         call mexerrmsgidandtxt ('MATLAB:zggev:BadAllocation',
     +       'The allocation of variables was unsuccesfull.')
      endif
C
C     Copy input data to the Fortran arrays
      call mxcopyptrtocomplex16(mxgetpr(prhs(1)),
     +                          mxgetpi(prhs(1)),h_eff,sizeh)
     
      call mxcopyptrtocomplex16(mxgetpr(prhs(2)),
     +                          mxgetpi(prhs(2)),b_eff,sizeb)     

      lwork = -1
C     Query for the workspace
      call zggev( 'V', 'V', hcols, h_eff, hrows,
     + b_eff, brows, 
     + alpha, beta, l_eigvec, hrows, r_eigvec, hrows, 
     + work1, lwork, rwork, 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 zggev( 'V', 'V', hcols, h_eff, hrows,
     + b_eff, hrows, 
     + alpha, beta, l_eigvec, hrows, r_eigvec, hrows, 
     + work, lwork, rwork, 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,1) ! ALPHA
      plhs(2) = mxcreatedoublematrix(hrows,cols,1) ! BETA
      plhs(3) = mxcreatedoublematrix(hrows,hcols,1) ! l_eigvec
      plhs(4) = mxcreatedoublematrix(hrows,hcols,1) ! r_eigvec
      
C     Load the data into to MATLAB outputs.      
      alpha_ptr_real = mxgetpr(plhs(1))
      alpha_ptr_imag = mxgetpi(plhs(1)) 
      call mxcopycomplex16toptr(alpha,alpha_ptr_real,alpha_ptr_imag,
     +                                                       hrows) 
C      
      beta_ptr_real = mxgetpr(plhs(2))
      beta_ptr_imag = mxgetpi(plhs(2))  
      call mxcopycomplex16toptr(beta,beta_ptr_real,beta_ptr_imag,
     +                                                     hrows)  
C      
      lvecs_ptr_real = mxgetpr(plhs(3))
      lvecs_ptr_imag = mxgetpi(plhs(3))      
      call mxcopycomplex16toptr(l_eigvec,lvecs_ptr_real,
     +                                     lvecs_ptr_imag,sizeb) 
C      
      rvecs_ptr_real = mxgetpr(plhs(4))
      rvecs_ptr_imag = mxgetpi(plhs(4))
      call mxcopycomplex16toptr(r_eigvec,rvecs_ptr_real,
     +                                      rvecs_ptr_imag,sizeb)
  
     

C     Free the allocated matrices
      deallocate(h_eff)
      deallocate(b_eff)
      deallocate(alpha)
      deallocate(beta)
      deallocate(l_eigvec)
      deallocate(r_eigvec)
      deallocate(work)
      deallocate(rwork)


      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