/*
LAPACK interface to solve Generalized Eigenvalue problem

Copyright (C) 2005 Stefan van der Walt

This is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version.

Octave is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with Octave; see the file COPYING.  If not, write to the Free
Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/

#include <octave/oct.h>
#include <f77-fcn.h>

extern "C"
{
   F77_RET_T
   F77_FUNC (zggev, ZGGEV) (F77_CONST_CHAR_ARG_DECL,
			    F77_CONST_CHAR_ARG_DECL,
			    const int& N,
			    Complex* A, const int& LDA,
			    Complex* B, const int& LDB,
			    Complex* ALPHA, Complex* BETA,
			    Complex* VL, const int& LDVL,
			    Complex* VR, const int& LDVR,
			    Complex* WORK,
			    const int& LWORK, double* RWORK, int& INFO);
}

DEFUN_DLD(gev, args, nargout, "Generalised Eigenvalue Test\n\
Usage: gev(A,B)")
{
    octave_value_list retval;

    if (args.length() != 2) {
	print_usage ("gev");
	return retval;
    }

    ComplexMatrix A = args(0).complex_matrix_value();
    ComplexMatrix B = args(1).complex_matrix_value();
    if (error_state) {
	error ("gev: invalid matrix arguments");
	return retval;
    }

    int n = A.rows();
    if ((n != A.columns()) || (n != B.columns()) || (n != B.rows())) {
	error ("gev: requires two square matrices as input");
	return retval;
    }

    Complex* a_tmp = A.fortran_vec ();
    Complex* b_tmp = B.fortran_vec ();

    Array<Complex> alpha (n);
    Complex* p_alpha = alpha.fortran_vec ();

    Array<Complex> beta (n);
    Complex* p_beta = beta.fortran_vec ();

    Array<double> rwork (8*n);
    double* p_rwork = rwork.fortran_vec ();

    int idummy = 1;
    Complex* dummy = 0;

    Complex work_length;
    
    int lwork = -1; // workspace query
    int info = 0;

    F77_FUNC (zggev, ZGGEV) (F77_CONST_CHAR_ARG2 ("N", 1),
     			     F77_CONST_CHAR_ARG2 ("N", 1),
 			     n,
			     a_tmp, n,
			     b_tmp, n,
 			     p_alpha, p_beta,
 			     dummy, idummy,
 			     dummy, idummy,
 			     &work_length, lwork, p_rwork,
 			     info);

    if (! f77_exception_encountered && info == 0) {
	lwork = (int)work_length.real();

	Array<Complex> work (lwork);
	Complex* p_work = work.fortran_vec ();

	F77_FUNC (zggev, ZGGEV) (F77_CONST_CHAR_ARG2 ("N", 1),
				 F77_CONST_CHAR_ARG2 ("N", 1),
				 n,
				 a_tmp, n,
				 b_tmp, n,
				 p_alpha, p_beta,
				 dummy, idummy,
				 dummy, idummy,
				 p_work, lwork, p_rwork,
				 info);

	if (f77_exception_encountered || info < 0) {
	    error("gev: unrecoverable error in zggev");
	    return retval;
	}

	else if (info >= 1 && info <= n) {
	    error("gev: QZ iteration failed in zggev");
	}

	else if (info == n+1) {
	    error("gev: iteration failed in dhgeqz called from zggev");
	}

	else if (info == n+2) {
	    error("gev: error return from dtgevc called from zggev");
	}

	ComplexColumnVector lambda (n);

	for (int i = 0; i < n; i++) {
	    lambda(i) = alpha(i) / beta(i);
	}

	retval.append (lambda);
	
	
    } else {
	error("gev: zggev workspace query failed");
	return retval;
    }

    return retval;
}

/*
%!test
%! a=[2,i,4;3,3*i,2;4,6,7];
%! b=[4,4,4;3,6*i,2;4,i,7];
%! gev(a,b)
*/