/*
LAPACK interface to solve generalised eigenvalue problems

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
			    F77_CHAR_ARG_LEN_DECL
			    F77_CHAR_ARG_LEN_DECL
			    );
}

DEFUN_DLD(gev, args, nargout,
"-*- texinfo -*-\n\
@deftypefn {Loadable Function} address@hidden =} gev (@var{A}, @var{B}))\n\
Solve the generalized eigenvalue problem.\n\
\n\
Usage: E = gev(A,B)\n\
Usage: [V, D] = gev(A,B)\n\
@end deftypefn\n\
\n\
@seealso{eig, qz}")
{
    octave_value_list retval;

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

    bool calc_evec = (nargout > 1);

    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 ldvr = calc_evec ? n : 1;
    ComplexMatrix vr (ldvr, ldvr);
    Complex* p_vr = vr.fortran_vec ();

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

    int idummy = 1;
    Complex* dummy = 0;

    F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
     			     F77_CONST_CHAR_ARG2 (calc_evec ? "V" : "N", 1),
 			     n,
			     a_tmp, n,
			     b_tmp, n,
 			     p_alpha, p_beta,
 			     dummy, idummy,
 			     p_vr, n,
 			     &work_length, lwork, p_rwork,
 			     info
			     F77_CHAR_ARG_LEN (1)
			     F77_CHAR_ARG_LEN (1)
			     ));

    if (error_state || f77_exception_encountered || info != 0) {
	error("gev: zggev workspace query failed");
	return retval;
    }

    lwork = (int)work_length.real();
    
    Array<Complex> work (lwork);
    Complex* p_work = work.fortran_vec ();
    
    F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
			     F77_CONST_CHAR_ARG2 (calc_evec ? "V" : "N", 1),
			     n,
			     a_tmp, n,
			     b_tmp, n,
			     p_alpha, p_beta,
			     dummy, idummy,
			     p_vr, n,
			     p_work, lwork, p_rwork,
			     info
			     F77_CHAR_ARG_LEN (1)
			     F77_CHAR_ARG_LEN (1)
			     ));
    
    if (error_state || 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);
    }

    if (calc_evec) {
	retval.append (vr);
	retval.append (ComplexDiagMatrix (lambda) );
    } else {
	retval.append (lambda); 
    }

   
    return retval;
}

/*
%!shared A, B
%!
%!test
%! A = [2,i,4; 3,3*i,2; 4,6,7];
%! B = [4,4,4; 3,6*i,2; 4,i,7];
%! lambda = gev(A,B);
%! assert(lambda, [-0.56498 - 0.16150i;
%!                  0.90696 + 0.29594i;
%!                  1.00000 - 0.00000i], 1e-4);
%!
%!test
%! [E,V] = gev(A,B);
%! assert(diag(V), [-0.56498 - 0.16150i;
%!                   0.90696 + 0.29594i;
%!                   1.00000 - 0.00000i], 1e-4);
%! assert(E, [-0.71544 + 0.28456i,  0.27236 + 0.72764i,  0.00000 + 0.00000i;
%!            -0.20817 - 0.25941i, -0.00880 - 0.20488i,  0.00000 + 0.00000i;
%!             0.52058 - 0.02514i, -0.68547 - 0.13278i, -0.15558 - 0.84442i], 1e-4);
%!
*/

/*
%!demo
%!
%! # Solve the generalised eigenvalue problem
%!
%! A = [2,i,4; 3,3*i,2; 4,6,7]
%! B = [4,4,4; 3,6*i,2; 4,i,7]
%! [E, V] = gev(A,B)
*/