# To unbundle, sh this file
echo f2c.h 1>&2
cat >f2c.h <<'End of f2c.h'
/* f2c.h  --  Standard Fortran to C header file */

/**  barf  [ba:rf]  2.  "He suggested using FORTRAN, and everybody barfed."

	- From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */

#ifndef F2C_INCLUDE
#define F2C_INCLUDE

typedef long int integer;
typedef char *address;
typedef short int shortint;
typedef float real;
typedef double doublereal;
typedef struct { real r, i; } complex;
typedef struct { doublereal r, i; } doublecomplex;
typedef long int logical;
typedef short int shortlogical;
typedef char logical1;
typedef char integer1;

#define TRUE_ (1)
#define FALSE_ (0)

/* Extern is for use with -E */
#ifndef Extern
#define Extern extern
#endif

/* I/O stuff */

#ifdef f2c_i2
/* for -i2 */
typedef short flag;
typedef short ftnlen;
typedef short ftnint;
#else
typedef long flag;
typedef long ftnlen;
typedef long ftnint;
#endif

/*external read, write*/
typedef struct
{	flag cierr;
	ftnint ciunit;
	flag ciend;
	char *cifmt;
	ftnint cirec;
} cilist;

/*internal read, write*/
typedef struct
{	flag icierr;
	char *iciunit;
	flag iciend;
	char *icifmt;
	ftnint icirlen;
	ftnint icirnum;
} icilist;

/*open*/
typedef struct
{	flag oerr;
	ftnint ounit;
	char *ofnm;
	ftnlen ofnmlen;
	char *osta;
	char *oacc;
	char *ofm;
	ftnint orl;
	char *oblnk;
} olist;

/*close*/
typedef struct
{	flag cerr;
	ftnint cunit;
	char *csta;
} cllist;

/*rewind, backspace, endfile*/
typedef struct
{	flag aerr;
	ftnint aunit;
} alist;

/* inquire */
typedef struct
{	flag inerr;
	ftnint inunit;
	char *infile;
	ftnlen infilen;
	ftnint	*inex;	/*parameters in standard's order*/
	ftnint	*inopen;
	ftnint	*innum;
	ftnint	*innamed;
	char	*inname;
	ftnlen	innamlen;
	char	*inacc;
	ftnlen	inacclen;
	char	*inseq;
	ftnlen	inseqlen;
	char 	*indir;
	ftnlen	indirlen;
	char	*infmt;
	ftnlen	infmtlen;
	char	*inform;
	ftnint	informlen;
	char	*inunf;
	ftnlen	inunflen;
	ftnint	*inrecl;
	ftnint	*innrec;
	char	*inblank;
	ftnlen	inblanklen;
} inlist;

#define VOID void

union Multitype {	/* for multiple entry points */
	shortint h;
	integer i;
	real r;
	doublereal d;
	complex c;
	doublecomplex z;
	};

typedef union Multitype Multitype;

typedef long Long;	/* No longer used; formerly in Namelist */

struct Vardesc {	/* for Namelist */
	char *name;
	char *addr;
	ftnlen *dims;
	int  type;
	};
typedef struct Vardesc Vardesc;

struct Namelist {
	char *name;
	Vardesc **vars;
	int nvars;
	};
typedef struct Namelist Namelist;

#define abs(x) ((x) >= 0 ? (x) : -(x))
#define dabs(x) (doublereal)abs(x)
#define min(a,b) ((a) <= (b) ? (a) : (b))
#define max(a,b) ((a) >= (b) ? (a) : (b))
#define dmin(a,b) (doublereal)min(a,b)
#define dmax(a,b) (doublereal)max(a,b)

/* procedure parameter types for -A and -C++ */

#define F2C_proc_par_types 1
#ifdef __cplusplus
typedef int /* Unknown procedure type */ (*U_fp)(...);
typedef shortint (*J_fp)(...);
typedef integer (*I_fp)(...);
typedef real (*R_fp)(...);
typedef doublereal (*D_fp)(...), (*E_fp)(...);
typedef /* Complex */ VOID (*C_fp)(...);
typedef /* Double Complex */ VOID (*Z_fp)(...);
typedef logical (*L_fp)(...);
typedef shortlogical (*K_fp)(...);
typedef /* Character */ VOID (*H_fp)(...);
typedef /* Subroutine */ int (*S_fp)(...);
#else
typedef int /* Unknown procedure type */ (*U_fp)();
typedef shortint (*J_fp)();
typedef integer (*I_fp)();
typedef real (*R_fp)();
typedef doublereal (*D_fp)(), (*E_fp)();
typedef /* Complex */ VOID (*C_fp)();
typedef /* Double Complex */ VOID (*Z_fp)();
typedef logical (*L_fp)();
typedef shortlogical (*K_fp)();
typedef /* Character */ VOID (*H_fp)();
typedef /* Subroutine */ int (*S_fp)();
#endif
/* E_fp is for real functions when -R is not specified */
typedef VOID C_f;	/* complex function */
typedef VOID H_f;	/* character function */
typedef VOID Z_f;	/* double complex function */
typedef doublereal E_f;	/* real function with -R not specified */

/* undef any lower-case symbols that your C compiler predefines, e.g.: */

#ifndef Skip_f2c_Undefs
#undef cray
#undef gcos
#undef mc68010
#undef mc68020
#undef mips
#undef pdp11
#undef sgi
#undef sparc
#undef sun
#undef sun2
#undef sun3
#undef sun4
#undef u370
#undef u3b
#undef u3b2
#undef u3b5
#undef unix
#undef vax
#endif
#endif
End of f2c.h
echo matrix.h 1>&2
cat >matrix.h <<'End of matrix.h'
/*************************************************************************
**************************************************************************
**************************************************************************
******                                                                ****
******  Matrix. A C++ matrix class                                    ****
******                                                                ****
******  Copyright (c) 1993  Universities Research Association, Inc.   ****
******                All Rights Reserved                             ****
******                                                                ****
******  Author:    James Holt                                         ****
******                                                                ****
******   Taken and modified from the book C++ Inside and out          ****
******   by Bruce Eckel                                               ****
******                                                                ****
******                                                                ****
******                                                                ****
******   Eigen Vector/Value routine taken from eispack                ****
******                                                                ****
******             Fermilab                                           ****
******             P.O.Box 500                                        ****
******             Mail Stop 345                                      ****
******             Batavia, IL   60510                                ****
******                                                                ****
******             Phone: (708) 840 8385                              ****
******             Email: holt@calvin.fnal.gov                        ****
******                                                                ****
******  Release    Version 1.0 : October  4, 1993                     ****
******                                                                ****
**************************************************************************
**************************************************************************
******                                                                ****
******  This material resulted from work developed under a            ****
******  Government Contract and is subject to the following           ****
******  license: The Government retains a paid-up, nonexclusive,      ****
******  irrevocable worldwide license to reproduce, prepare           ****
******  derivative works, perform publicly and display publicly       ****
******  by or for the Government, including the right to distribute   ****
******  to other Government contractors.  Neither the United          ****
******  States nor the United States Department of Energy, nor        ****
******  any of their employees, makes any warranty, express or        ****
******  implied, or assumes any legal liability or responsibility     ****
******  for the accuracy, completeness or usefulness of any           ****
******  information, apparatus, product, or process disclosed, or     ****
******  represents that its use would not infringe privately owned    ****
******  rights.                                                       ****
******                                                                ****
******  These files are made avaliable for use by specific            ****
******  individuals or a specific group.  They are not to be          ****
******  reproduced or redistributed.                                  ****
******                                                                ****
******  These files are provided "as is" with no warranties of any    ****
******  kind including the warranties of design,                      ****
******  merchantibility and fitness for a particular purpose,         ****
******  or arising from a course of dealing, usage or trade           ****
******  practice.                                                     ****
******                                                                ****
******  These files are provided with no support or obligation on     ****
******  the part of U.R.A. or Fermilab to assist in their use,        ****
******  correction, modification or enhancement.                      ****
******                                                                ****
******  Neither U.R.A. nor Fermilab shall be liable with              ****
******  respect to the infringement of copyrights, trade              ****
******  secrets or any patents by these files or any part             ****
******  thereof.                                                      ****
******                                                                ****
******  In no event will U.R.A., Fermilab, or the author(s) be        ****
******  liable for any lost revenue or profits or other               ****
******  special, indirect and consequential damages, even if          ****
******  they have been advised of the possibility of such             ****
******  damages.                                                      ****
******                                                                ****
**************************************************************************
*************************************************************************/

#ifndef MATRIX_H_
#define MATRIX_H_


#include <stdlib.h>
#include <iostream.h>
#include <fstream.h>
#include <iomanip.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <complex.h>


class complexmatrix;

struct doublematrep {
  double **m;       // pointer to the doublematrix
  int r,c;          // number of rows and columns
  int n;            // reference count
};


class doublematrix {
 private:
  doublematrep *p;  
  void error(char * msgl, char * msg2 = " ")const ;  // private function

// functions used by inverse() and determinant {

  void switch_columns(int col1, int col2);
  void switch_rows(int row1, int row2);
  void copy_column(doublematrix& m, int from_col, int to_col);
  doublematrix scale(); // Scale a matrix (used in L-U decomposition
  void deepcopy(doublematrix& from, doublematrix& to);  // make an image
//
// Returns the L-U decomposition of a matrix
//
  doublematrix lu_decompose (int* indx,  int& d ) ;
//
// Uses L-U decomposition for matrix inverse.
//
  void lu_back_subst (int* indx, doublematrix& b);

// Used by matrix functions which know they aren't
// exceeding the boundaries
//
  double & mval(int row, int col){
     return (p->m[row][col]);
   }


friend class complexmatrix;

public:

doublematrix(int mrows = 1, int columns = 1, double initval = 0);
doublematrix(int mrows, int columns, double* initvalues);
doublematrix(char * flag, int dimension);// create an ident matrix
					 // or a symplectic matrix
doublematrix(char * matfile);           // read from a "standard" doublematrix
doublematrix(doublematrix& X);          // copy constructor
~doublematrix ( ) ;                     // destructor

int rows() const {  return p->r; };    // rows in matrix
int cols() const {  return p->c; };    // cols in matrix

// Write a "standard" matrix file:
void write_standard (char * filename, char * msg = " ");

				// operators

double& operator()(int i,int j);
double& operator()(int i);
doublematrix operator= (doublematrix& rval); // matrix assignement
doublematrix operator-(); // unary minus
friend doublematrix operator+ (const doublematrix &,const doublematrix &);
friend doublematrix operator+ (const doublematrix &,const double); 
friend doublematrix operator+ (const double ,const doublematrix &); 
friend doublematrix operator-(const doublematrix &,const doublematrix&); 
friend doublematrix operator-(const doublematrix &,const double); 
friend doublematrix operator-(const double,const doublematrix &); 
friend doublematrix operator*(const doublematrix &,const doublematrix&); 
friend complexmatrix operator*(const doublematrix &,const complexmatrix&); 
friend complexmatrix operator*(const complexmatrix &,const doublematrix&); 
friend doublematrix operator*(const doublematrix &,const double);
friend doublematrix operator*(const double,const doublematrix &);
//
// *** the following can be used to read or write an element.
//
double & val(int row, int col); // element selection;

doublematrix transpose(); // transpose a square matrix
double determinant();
doublematrix inverse();
complexmatrix& eigenValues();
complexmatrix& eigenVectors();
double mmin();  // find minimum element in the matrix
double mmax();  // find maximum element in the matrix

double mean(); //  average all the elements of the matrix
double variance(); // statistical variance of all elements
void print(char * msg = " "); // print matrix with a message

};

ostream& operator<<(ostream&, doublematrix&);

struct complexmatrep {
  complex **m;       // pointer to the complexmatrix
  int r,c;          // number of rows and columns
  int n;            // reference count
};


class complexmatrix {

  complexmatrep *p;  
  void error(char * msgl, char * msg2 = " ") const ;  // private function
// functions used by inverse() and determinant 
void switch_columns(int col1, int col2);
void switch_rows(int row1, int row2);
void copy_column(complexmatrix& m, int from_col, int to_col);
complexmatrix scale(); // Scale a matrix (used in L-U decomposition
void deepcopy(complexmatrix& from, complexmatrix& to);  // make an image
//
// Returns the L-U decomposition of a matrix
//
complexmatrix lu_decompose (int* indx,  int& d ) ;
//
// Uses L-U decomposition for matrix inverse
//
void lu_back_subst (int* indx, complexmatrix& b);
//
// used by matrix functions which know they aren't
// exceeding the boundaries
//
complex & mval(int row, int col){ return (p->m[row][col]);}

friend class doublematrix;
public:

complexmatrix(doublematrix& x);
complexmatrix(int mrows , int columns , double* initvalues);
complexmatrix(int mrows = 1, int columns = 1, complex initval = complex_zero);
complexmatrix(int mrows, int columns, complex* initvalues);
complexmatrix(char * flag, int dimension);  // create an ident matrix
complexmatrix(char * matfile);          // read from a "standard" complexmatrix
complexmatrix(complexmatrix& X);        // copy constructor
~complexmatrix ( ) ;                    // destructor

// Write a "standard" matrix file:
void write_standard (char * filename, char * msg = " ");

int rows() const {  return p->r; };    // rows in matrix
int cols() const {  return p->c; };    // cols in matrix

				       // operators

complex& operator()(int i,int j);
complex& operator()(int i);
complexmatrix operator= (complexmatrix& rval); // matrix assignement
complexmatrix operator-(); // unary minus

friend complexmatrix operator+ (const complexmatrix &,const complexmatrix &);
friend complexmatrix operator+ (const complexmatrix &,const complex); 
friend complexmatrix operator+ (const complex ,const complexmatrix &); 
friend complexmatrix operator-(const complexmatrix &,const complexmatrix&); 
friend complexmatrix operator-(const complexmatrix &,const complex); 
friend complexmatrix operator-(const complex,const complexmatrix &); 
friend complexmatrix operator*(const complexmatrix &,const complexmatrix&); 
friend complexmatrix operator*(const doublematrix &,const complexmatrix&); 
friend complexmatrix operator*(const complexmatrix &,const doublematrix&); 
friend complexmatrix operator*(const complexmatrix &,const complex);
friend complexmatrix operator*(const complex,const complexmatrix &);

//
// *** can be used to read or write an element.
//
complex & val(int row, int col); // element selection;

complexmatrix transpose(); // transpose a square matrix
complex determinant();
complexmatrix inverse();
complex mmin();  // find minimum element in the matrix
complex mmax();  // find maximum element in the matrix
complexmatrix& eigenValues();
complexmatrix& eigenVectors();

void print(char * msg = " "); // print matrix with a message

};

ostream& operator<<(ostream&, complexmatrix&);

#endif // MATRIX_H_
End of matrix.h
echo cg.c 1>&2
cat >cg.c <<'End of cg.c'
/* cg.f -- translated by f2c (version of 3 February 1990  3:36:42).
   from the eispack library
*/

#include "f2c.h"

/* Subroutine */ int cg_(nm, n, ar, ai, wr, wi, matz, zr, zi, fv1, fv2, fv3, 
	ierr)
integer *nm, *n;
doublereal *ar, *ai, *wr, *wi;
integer *matz;
doublereal *zr, *zi, *fv1, *fv2, *fv3;
integer *ierr;
{
    /* System generated locals */
    integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset, 
	    zi_dim1, zi_offset;

    /* Local variables */
    extern /* Subroutine */ int cbal_(), corth_(), comqr_(), cbabk2_(), 
	    comqr2_();
    static integer is1, is2;

    /* Parameter adjustments */
    ar_dim1 = *nm;
    ar_offset = ar_dim1 + 1;
    ar -= ar_offset;
    ai_dim1 = *nm;
    ai_offset = ai_dim1 + 1;
    ai -= ai_offset;
    --wr;
    --wi;
    zr_dim1 = *nm;
    zr_offset = zr_dim1 + 1;
    zr -= zr_offset;
    zi_dim1 = *nm;
    zi_offset = zi_dim1 + 1;
    zi -= zi_offset;
    --fv1;
    --fv2;
    --fv3;

    /* Function Body */


/*     this subroutine calls the recommended sequence of */
/*     subroutines from the eigensystem subroutine package (eispack) */
/*     to find the eigenvalues and eigenvectors (if desired) */
/*     of a complex general matrix. */

/*     on input */

/*        nm  must be set to the row dimension of the two-dimensional */
/*        array parameters as declared in the calling program */
/*        dimension statement. */

/*        n  is the order of the matrix  a=(ar,ai). */

/*        ar  and  ai  contain the real and imaginary parts, */
/*        respectively, of the complex general matrix. */

/*        matz  is an integer variable set equal to zero if */
/*        only eigenvalues are desired.  otherwise it is set to */
/*        any non-zero integer for both eigenvalues and eigenvectors. */

/*     on output */

/*        wr  and  wi  contain the real and imaginary parts, */
/*        respectively, of the eigenvalues. */

/*        zr  and  zi  contain the real and imaginary parts, */
/*        respectively, of the eigenvectors if matz is not zero. */

/*        ierr  is an integer output variable set equal to an error */
/*           completion code described in the documentation for comqr */
/*           and comqr2.  the normal completion code is zero. */

/*        fv1, fv2, and  fv3  are temporary storage arrays. */

/*     questions and comments should be directed to burton s. garbow, */
/*     mathematics and computer science div, argonne national laboratory 
*/

/*     this version dated august 1983. */

/*     ------------------------------------------------------------------ 
*/

    if (*n <= *nm) {
	goto L10;
    }
    *ierr = *n * 10;
    goto L50;

L10:
    cbal_(nm, n, &ar[ar_offset], &ai[ai_offset], &is1, &is2, &fv1[1]);
    corth_(nm, n, &is1, &is2, &ar[ar_offset], &ai[ai_offset], &fv2[1], &fv3[1]
	    );
    if (*matz != 0) {
	goto L20;
    }
/*     .......... find eigenvalues only .......... */
    comqr_(nm, n, &is1, &is2, &ar[ar_offset], &ai[ai_offset], &wr[1], &wi[1], 
	    ierr);
    goto L50;
/*     .......... find both eigenvalues and eigenvectors .......... */
L20:
    comqr2_(nm, n, &is1, &is2, &fv2[1], &fv3[1], &ar[ar_offset], &ai[
	    ai_offset], &wr[1], &wi[1], &zr[zr_offset], &zi[zi_offset], ierr);

    if (*ierr != 0) {
	goto L50;
    }
    cbabk2_(nm, n, &is1, &is2, &fv1[1], n, &zr[zr_offset], &zi[zi_offset]);
L50:
    return 0;
} /* cg_ */

/* Subroutine */ int cbabk2_(nm, n, low, igh, scale, m, zr, zi)
integer *nm, *n, *low, *igh;
doublereal *scale;
integer *m;
doublereal *zr, *zi;
{
    /* System generated locals */
    integer zr_dim1, zr_offset, zi_dim1, zi_offset, i_1, i_2;

    /* Local variables */
    static integer i, j, k;
    static doublereal s;
    static integer ii;

    /* Parameter adjustments */
    --scale;
    zr_dim1 = *nm;
    zr_offset = zr_dim1 + 1;
    zr -= zr_offset;
    zi_dim1 = *nm;
    zi_offset = zi_dim1 + 1;
    zi -= zi_offset;

    /* Function Body */


/*     this subroutine is a translation of the algol procedure */
/*     cbabk2, which is a complex version of balbak, */
/*     num. math. 13, 293-304(1969) by parlett and reinsch. */
/*     handbook for auto. comp., vol.ii-linear algebra, 315-326(1971). */

/*     this subroutine forms the eigenvectors of a complex general */
/*     matrix by back transforming those of the corresponding */
/*     balanced matrix determined by  cbal. */

/*     on input */

/*        nm must be set to the row dimension of two-dimensional */
/*          array parameters as declared in the calling program */
/*          dimension statement. */

/*        n is the order of the matrix. */

/*        low and igh are integers determined by  cbal. */

/*        scale contains information determining the permutations */
/*          and scaling factors used by  cbal. */

/*        m is the number of eigenvectors to be back transformed. */

/*        zr and zi contain the real and imaginary parts, */
/*          respectively, of the eigenvectors to be */
/*          back transformed in their first m columns. */

/*     on output */

/*        zr and zi contain the real and imaginary parts, */
/*          respectively, of the transformed eigenvectors */
/*          in their first m columns. */

/*     questions and comments should be directed to burton s. garbow, */
/*     mathematics and computer science div, argonne national laboratory 
*/

/*     this version dated august 1983. */

/*     ------------------------------------------------------------------ 
*/

    if (*m == 0) {
	goto L200;
    }
    if (*igh == *low) {
	goto L120;
    }

    i_1 = *igh;
    for (i = *low; i <= i_1; ++i) {
	s = scale[i];
/*     .......... left hand eigenvectors are back transformed */
/*                if the foregoing statement is replaced by */
/*                s=1.0d0/scale(i). .......... */
	i_2 = *m;
	for (j = 1; j <= i_2; ++j) {
	    zr[i + j * zr_dim1] *= s;
	    zi[i + j * zi_dim1] *= s;
/* L100: */
	}

/* L110: */
    }
/*     .......... for i=low-1 step -1 until 1, */
/*                igh+1 step 1 until n do -- .......... */
L120:
    i_1 = *n;
    for (ii = 1; ii <= i_1; ++ii) {
	i = ii;
	if (i >= *low && i <= *igh) {
	    goto L140;
	}
	if (i < *low) {
	    i = *low - ii;
	}
	k = (integer) scale[i];
	if (k == i) {
	    goto L140;
	}

	i_2 = *m;
	for (j = 1; j <= i_2; ++j) {
	    s = zr[i + j * zr_dim1];
	    zr[i + j * zr_dim1] = zr[k + j * zr_dim1];
	    zr[k + j * zr_dim1] = s;
	    s = zi[i + j * zi_dim1];
	    zi[i + j * zi_dim1] = zi[k + j * zi_dim1];
	    zi[k + j * zi_dim1] = s;
/* L130: */
	}

L140:
    ;}

L200:
    return 0;
} /* cbabk2_ */

/* Subroutine */ int comqr2_(nm, n, low, igh, ortr, orti, hr, hi, wr, wi, zr, 
	zi, ierr)
integer *nm, *n, *low, *igh;
doublereal *ortr, *orti, *hr, *hi, *wr, *wi, *zr, *zi;
integer *ierr;
{
    /* System generated locals */
    integer hr_dim1, hr_offset, hi_dim1, hi_offset, zr_dim1, zr_offset, 
	    zi_dim1, zi_offset, i_1, i_2, i_3;
    doublereal d_1, d_2, d_3, d_4;

    /* Local variables */
    static integer iend;
    extern /* Subroutine */ int cdiv_();
    static doublereal norm;
    static integer i, j, k, l, m, ii, en, jj, ll, nn;
    static doublereal si, ti, xi, yi, sr, tr, xr, yr;
    extern doublereal pythag_();
    extern /* Subroutine */ int csroot_();
    static integer ip1, lp1, itn, its;
    static doublereal zzi, zzr;
    static integer enm1;
    static doublereal tst1, tst2;

    /* Parameter adjustments */
    --ortr;
    --orti;
    hr_dim1 = *nm;
    hr_offset = hr_dim1 + 1;
    hr -= hr_offset;
    hi_dim1 = *nm;
    hi_offset = hi_dim1 + 1;
    hi -= hi_offset;
    --wr;
    --wi;
    zr_dim1 = *nm;
    zr_offset = zr_dim1 + 1;
    zr -= zr_offset;
    zi_dim1 = *nm;
    zi_offset = zi_dim1 + 1;
    zi -= zi_offset;

    /* Function Body */
/*  MESHED overflow control WITH vectors of isolated roots (10/19/89 BSG) 
*/
/*  MESHED overflow control WITH triangular multiply (10/30/89 BSG) */


/*     this subroutine is a translation of a unitary analogue of the */
/*     algol procedure  comlr2, num. math. 16, 181-204(1970) by peters */
/*     and wilkinson. */
/*     handbook for auto. comp., vol.ii-linear algebra, 372-395(1971). */
/*     the unitary analogue substitutes the qr algorithm of francis */
/*     (comp. jour. 4, 332-345(1962)) for the lr algorithm. */

/*     this subroutine finds the eigenvalues and eigenvectors */
/*     of a complex upper hessenberg matrix by the qr */
/*     method.  the eigenvectors of a complex general matrix */
/*     can also be found if  corth  has been used to reduce */
/*     this general matrix to hessenberg form. */

/*     on input */

/*        nm must be set to the row dimension of two-dimensional */
/*          array parameters as declared in the calling program */
/*          dimension statement. */

/*        n is the order of the matrix. */

/*        low and igh are integers determined by the balancing */
/*          subroutine  cbal.  if  cbal  has not been used, */
/*          set low=1, igh=n. */

/*        ortr and orti contain information about the unitary trans- */
/*          formations used in the reduction by  corth, if performed. */
/*          only elements low through igh are used.  if the eigenvectors 
*/
/*          of the hessenberg matrix are desired, set ortr(j) and */
/*          orti(j) to 0.0d0 for these elements. */

/*        hr and hi contain the real and imaginary parts, */
/*          respectively, of the complex upper hessenberg matrix. */
/*          their lower triangles below the subdiagonal contain further */

/*          information about the transformations which were used in the 
*/
/*          reduction by  corth, if performed.  if the eigenvectors of */
/*          the hessenberg matrix are desired, these elements may be */
/*          arbitrary. */

/*     on output */

/*        ortr, orti, and the upper hessenberg portions of hr and hi */
/*          have been destroyed. */

/*        wr and wi contain the real and imaginary parts, */
/*          respectively, of the eigenvalues.  if an error */
/*          exit is made, the eigenvalues should be correct */
/*          for indices ierr+1,...,n. */

/*        zr and zi contain the real and imaginary parts, */
/*          respectively, of the eigenvectors.  the eigenvectors */
/*          are unnormalized.  if an error exit is made, none of */
/*          the eigenvectors has been found. */

/*        ierr is set to */
/*          zero       for normal return, */
/*          j          if the limit of 30*n iterations is exhausted */
/*                     while the j-th eigenvalue is being sought. */

/*     calls cdiv for complex division. */
/*     calls csroot for complex square root. */
/*     calls pythag for  dsqrt(a*a + b*b) . */

/*     questions and comments should be directed to burton s. garbow, */
/*     mathematics and computer science div, argonne national laboratory 
*/

/*     this version dated october 1989. */

/*     ------------------------------------------------------------------ 
*/

    *ierr = 0;
/*     .......... initialize eigenvector matrix .......... */
    i_1 = *n;
    for (j = 1; j <= i_1; ++j) {

	i_2 = *n;
	for (i = 1; i <= i_2; ++i) {
	    zr[i + j * zr_dim1] = 0.;
	    zi[i + j * zi_dim1] = 0.;
/* L100: */
	}
	zr[j + j * zr_dim1] = 1.;
/* L101: */
    }
/*     .......... form the matrix of accumulated transformations */
/*                from the information left by corth .......... */
    iend = *igh - *low - 1;
    if (iend < 0) {
	goto L180;
    } else if (iend == 0) {
	goto L150;
    } else {
	goto L105;
    }
/*     .......... for i=igh-1 step -1 until low+1 do -- .......... */
L105:
    i_1 = iend;
    for (ii = 1; ii <= i_1; ++ii) {
	i = *igh - ii;
	if (ortr[i] == 0. && orti[i] == 0.) {
	    goto L140;
	}
	if (hr[i + (i - 1) * hr_dim1] == 0. && hi[i + (i - 1) * hi_dim1] == 
		0.) {
	    goto L140;
	}
/*     .......... norm below is negative of h formed in corth 
.......... */
	norm = hr[i + (i - 1) * hr_dim1] * ortr[i] + hi[i + (i - 1) * hi_dim1]
		 * orti[i];
	ip1 = i + 1;

	i_2 = *igh;
	for (k = ip1; k <= i_2; ++k) {
	    ortr[k] = hr[k + (i - 1) * hr_dim1];
	    orti[k] = hi[k + (i - 1) * hi_dim1];
/* L110: */
	}

	i_2 = *igh;
	for (j = i; j <= i_2; ++j) {
	    sr = 0.;
	    si = 0.;

	    i_3 = *igh;
	    for (k = i; k <= i_3; ++k) {
		sr = sr + ortr[k] * zr[k + j * zr_dim1] + orti[k] * zi[k + j *
			 zi_dim1];
		si = si + ortr[k] * zi[k + j * zi_dim1] - orti[k] * zr[k + j *
			 zr_dim1];
/* L115: */
	    }

	    sr /= norm;
	    si /= norm;

	    i_3 = *igh;
	    for (k = i; k <= i_3; ++k) {
		zr[k + j * zr_dim1] = zr[k + j * zr_dim1] + sr * ortr[k] - si 
			* orti[k];
		zi[k + j * zi_dim1] = zi[k + j * zi_dim1] + sr * orti[k] + si 
			* ortr[k];
/* L120: */
	    }

/* L130: */
	}

L140:
    ;}
/*     .......... create real subdiagonal elements .......... */
L150:
    l = *low + 1;

    i_1 = *igh;
    for (i = l; i <= i_1; ++i) {
/* Computing MAX */
	i_2 = i + 1;
	ll = min(*igh,i_2);
	if (hi[i + (i - 1) * hi_dim1] == 0.) {
	    goto L170;
	}
	norm = pythag_(&hr[i + (i - 1) * hr_dim1], &hi[i + (i - 1) * hi_dim1])
		;
	yr = hr[i + (i - 1) * hr_dim1] / norm;
	yi = hi[i + (i - 1) * hi_dim1] / norm;
	hr[i + (i - 1) * hr_dim1] = norm;
	hi[i + (i - 1) * hi_dim1] = 0.;

	i_2 = *n;
	for (j = i; j <= i_2; ++j) {
	    si = yr * hi[i + j * hi_dim1] - yi * hr[i + j * hr_dim1];
	    hr[i + j * hr_dim1] = yr * hr[i + j * hr_dim1] + yi * hi[i + j * 
		    hi_dim1];
	    hi[i + j * hi_dim1] = si;
/* L155: */
	}

	i_2 = ll;
	for (j = 1; j <= i_2; ++j) {
	    si = yr * hi[j + i * hi_dim1] + yi * hr[j + i * hr_dim1];
	    hr[j + i * hr_dim1] = yr * hr[j + i * hr_dim1] - yi * hi[j + i * 
		    hi_dim1];
	    hi[j + i * hi_dim1] = si;
/* L160: */
	}

	i_2 = *igh;
	for (j = *low; j <= i_2; ++j) {
	    si = yr * zi[j + i * zi_dim1] + yi * zr[j + i * zr_dim1];
	    zr[j + i * zr_dim1] = yr * zr[j + i * zr_dim1] - yi * zi[j + i * 
		    zi_dim1];
	    zi[j + i * zi_dim1] = si;
/* L165: */
	}

L170:
    ;}
/*     .......... store roots isolated by cbal .......... */
L180:
    i_1 = *n;
    for (i = 1; i <= i_1; ++i) {
	if (i >= *low && i <= *igh) {
	    goto L200;
	}
	wr[i] = hr[i + i * hr_dim1];
	wi[i] = hi[i + i * hi_dim1];
L200:
    ;}

    en = *igh;
    tr = 0.;
    ti = 0.;
    itn = *n * 30;
/*     .......... search for next eigenvalue .......... */
L220:
    if (en < *low) {
	goto L680;
    }
    its = 0;
    enm1 = en - 1;
/*     .......... look for single small sub-diagonal element */
/*                for l=en step -1 until low do -- .......... */
L240:
    i_1 = en;
    for (ll = *low; ll <= i_1; ++ll) {
	l = en + *low - ll;
	if (l == *low) {
	    goto L300;
	}
	tst1 = (d_1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[l 
		- 1 + (l - 1) * hi_dim1], abs(d_2)) + (d_3 = hr[l + l * 
		hr_dim1], abs(d_3)) + (d_4 = hi[l + l * hi_dim1], abs(d_4));
	tst2 = tst1 + (d_1 = hr[l + (l - 1) * hr_dim1], abs(d_1));
	if (tst2 == tst1) {
	    goto L300;
	}
/* L260: */
    }
/*     .......... form shift .......... */
L300:
    if (l == en) {
	goto L660;
    }
    if (itn == 0) {
	goto L1000;
    }
    if (its == 10 || its == 20) {
	goto L320;
    }
    sr = hr[en + en * hr_dim1];
    si = hi[en + en * hi_dim1];
    xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1];
    xi = hi[enm1 + en * hi_dim1] * hr[en + enm1 * hr_dim1];
    if (xr == 0. && xi == 0.) {
	goto L340;
    }
    yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
    yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
/* Computing 2nd power */
    d_2 = yr;
/* Computing 2nd power */
    d_3 = yi;
    d_1 = d_2 * d_2 - d_3 * d_3 + xr;
    d_4 = yr * 2. * yi + xi;
    csroot_(&d_1, &d_4, &zzr, &zzi);
    if (yr * zzr + yi * zzi >= 0.) {
	goto L310;
    }
    zzr = -zzr;
    zzi = -zzi;
L310:
    d_1 = yr + zzr;
    d_2 = yi + zzi;
    cdiv_(&xr, &xi, &d_1, &d_2, &xr, &xi);
    sr -= xr;
    si -= xi;
    goto L340;
/*     .......... form exceptional shift .......... */
L320:
    sr = (d_1 = hr[en + enm1 * hr_dim1], abs(d_1)) + (d_2 = hr[enm1 + (en - 2)
	     * hr_dim1], abs(d_2));
    si = 0.;

L340:
    i_1 = en;
    for (i = *low; i <= i_1; ++i) {
	hr[i + i * hr_dim1] -= sr;
	hi[i + i * hi_dim1] -= si;
/* L360: */
    }

    tr += sr;
    ti += si;
    ++its;
    --itn;
/*     .......... reduce to triangle (rows) .......... */
    lp1 = l + 1;

    i_1 = en;
    for (i = lp1; i <= i_1; ++i) {
	sr = hr[i + (i - 1) * hr_dim1];
	hr[i + (i - 1) * hr_dim1] = 0.;
	d_1 = pythag_(&hr[i - 1 + (i - 1) * hr_dim1], &hi[i - 1 + (i - 1) * 
		hi_dim1]);
	norm = pythag_(&d_1, &sr);
	xr = hr[i - 1 + (i - 1) * hr_dim1] / norm;
	wr[i - 1] = xr;
	xi = hi[i - 1 + (i - 1) * hi_dim1] / norm;
	wi[i - 1] = xi;
	hr[i - 1 + (i - 1) * hr_dim1] = norm;
	hi[i - 1 + (i - 1) * hi_dim1] = 0.;
	hi[i + (i - 1) * hi_dim1] = sr / norm;

	i_2 = *n;
	for (j = i; j <= i_2; ++j) {
	    yr = hr[i - 1 + j * hr_dim1];
	    yi = hi[i - 1 + j * hi_dim1];
	    zzr = hr[i + j * hr_dim1];
	    zzi = hi[i + j * hi_dim1];
	    hr[i - 1 + j * hr_dim1] = xr * yr + xi * yi + hi[i + (i - 1) * 
		    hi_dim1] * zzr;
	    hi[i - 1 + j * hi_dim1] = xr * yi - xi * yr + hi[i + (i - 1) * 
		    hi_dim1] * zzi;
	    hr[i + j * hr_dim1] = xr * zzr - xi * zzi - hi[i + (i - 1) * 
		    hi_dim1] * yr;
	    hi[i + j * hi_dim1] = xr * zzi + xi * zzr - hi[i + (i - 1) * 
		    hi_dim1] * yi;
/* L490: */
	}

/* L500: */
    }

    si = hi[en + en * hi_dim1];
    if (si == 0.) {
	goto L540;
    }
    norm = pythag_(&hr[en + en * hr_dim1], &si);
    sr = hr[en + en * hr_dim1] / norm;
    si /= norm;
    hr[en + en * hr_dim1] = norm;
    hi[en + en * hi_dim1] = 0.;
    if (en == *n) {
	goto L540;
    }
    ip1 = en + 1;

    i_1 = *n;
    for (j = ip1; j <= i_1; ++j) {
	yr = hr[en + j * hr_dim1];
	yi = hi[en + j * hi_dim1];
	hr[en + j * hr_dim1] = sr * yr + si * yi;
	hi[en + j * hi_dim1] = sr * yi - si * yr;
/* L520: */
    }
/*     .......... inverse operation (columns) .......... */
L540:
    i_1 = en;
    for (j = lp1; j <= i_1; ++j) {
	xr = wr[j - 1];
	xi = wi[j - 1];

	i_2 = j;
	for (i = 1; i <= i_2; ++i) {
	    yr = hr[i + (j - 1) * hr_dim1];
	    yi = 0.;
	    zzr = hr[i + j * hr_dim1];
	    zzi = hi[i + j * hi_dim1];
	    if (i == j) {
		goto L560;
	    }
	    yi = hi[i + (j - 1) * hi_dim1];
	    hi[i + (j - 1) * hi_dim1] = xr * yi + xi * yr + hi[j + (j - 1) * 
		    hi_dim1] * zzi;
L560:
	    hr[i + (j - 1) * hr_dim1] = xr * yr - xi * yi + hi[j + (j - 1) * 
		    hi_dim1] * zzr;
	    hr[i + j * hr_dim1] = xr * zzr + xi * zzi - hi[j + (j - 1) * 
		    hi_dim1] * yr;
	    hi[i + j * hi_dim1] = xr * zzi - xi * zzr - hi[j + (j - 1) * 
		    hi_dim1] * yi;
/* L580: */
	}

	i_2 = *igh;
	for (i = *low; i <= i_2; ++i) {
	    yr = zr[i + (j - 1) * zr_dim1];
	    yi = zi[i + (j - 1) * zi_dim1];
	    zzr = zr[i + j * zr_dim1];
	    zzi = zi[i + j * zi_dim1];
	    zr[i + (j - 1) * zr_dim1] = xr * yr - xi * yi + hi[j + (j - 1) * 
		    hi_dim1] * zzr;
	    zi[i + (j - 1) * zi_dim1] = xr * yi + xi * yr + hi[j + (j - 1) * 
		    hi_dim1] * zzi;
	    zr[i + j * zr_dim1] = xr * zzr + xi * zzi - hi[j + (j - 1) * 
		    hi_dim1] * yr;
	    zi[i + j * zi_dim1] = xr * zzi - xi * zzr - hi[j + (j - 1) * 
		    hi_dim1] * yi;
/* L590: */
	}

/* L600: */
    }

    if (si == 0.) {
	goto L240;
    }

    i_1 = en;
    for (i = 1; i <= i_1; ++i) {
	yr = hr[i + en * hr_dim1];
	yi = hi[i + en * hi_dim1];
	hr[i + en * hr_dim1] = sr * yr - si * yi;
	hi[i + en * hi_dim1] = sr * yi + si * yr;
/* L630: */
    }

    i_1 = *igh;
    for (i = *low; i <= i_1; ++i) {
	yr = zr[i + en * zr_dim1];
	yi = zi[i + en * zi_dim1];
	zr[i + en * zr_dim1] = sr * yr - si * yi;
	zi[i + en * zi_dim1] = sr * yi + si * yr;
/* L640: */
    }

    goto L240;
/*     .......... a root found .......... */
L660:
    hr[en + en * hr_dim1] += tr;
    wr[en] = hr[en + en * hr_dim1];
    hi[en + en * hi_dim1] += ti;
    wi[en] = hi[en + en * hi_dim1];
    en = enm1;
    goto L220;
/*     .......... all roots found.  backsubstitute to find */
/*                vectors of upper triangular form .......... */
L680:
    norm = 0.;

    i_1 = *n;
    for (i = 1; i <= i_1; ++i) {

	i_2 = *n;
	for (j = i; j <= i_2; ++j) {
	    tr = (d_1 = hr[i + j * hr_dim1], abs(d_1)) + (d_2 = hi[i + j * 
		    hi_dim1], abs(d_2));
	    if (tr > norm) {
		norm = tr;
	    }
/* L720: */
	}
    }

    if (*n == 1 || norm == 0.) {
	goto L1001;
    }
/*     .......... for en=n step -1 until 2 do -- .......... */
    i_2 = *n;
    for (nn = 2; nn <= i_2; ++nn) {
	en = *n + 2 - nn;
	xr = wr[en];
	xi = wi[en];
	hr[en + en * hr_dim1] = 1.;
	hi[en + en * hi_dim1] = 0.;
	enm1 = en - 1;
/*     .......... for i=en-1 step -1 until 1 do -- .......... */
	i_1 = enm1;
	for (ii = 1; ii <= i_1; ++ii) {
	    i = en - ii;
	    zzr = 0.;
	    zzi = 0.;
	    ip1 = i + 1;

	    i_3 = en;
	    for (j = ip1; j <= i_3; ++j) {
		zzr = zzr + hr[i + j * hr_dim1] * hr[j + en * hr_dim1] - hi[i 
			+ j * hi_dim1] * hi[j + en * hi_dim1];
		zzi = zzi + hr[i + j * hr_dim1] * hi[j + en * hi_dim1] + hi[i 
			+ j * hi_dim1] * hr[j + en * hr_dim1];
/* L740: */
	    }

	    yr = xr - wr[i];
	    yi = xi - wi[i];
	    if (yr != 0. || yi != 0.) {
		goto L765;
	    }
	    tst1 = norm;
	    yr = tst1;
L760:
	    yr *= .01;
	    tst2 = norm + yr;
	    if (tst2 > tst1) {
		goto L760;
	    }
L765:
	    cdiv_(&zzr, &zzi, &yr, &yi, &hr[i + en * hr_dim1], &hi[i + en * 
		    hi_dim1]);
/*     .......... overflow control .......... */
	    tr = (d_1 = hr[i + en * hr_dim1], abs(d_1)) + (d_2 = hi[i + en * 
		    hi_dim1], abs(d_2));
	    if (tr == 0.) {
		goto L780;
	    }
	    tst1 = tr;
	    tst2 = tst1 + 1. / tst1;
	    if (tst2 > tst1) {
		goto L780;
	    }
	    i_3 = en;
	    for (j = i; j <= i_3; ++j) {
		hr[j + en * hr_dim1] /= tr;
		hi[j + en * hi_dim1] /= tr;
/* L770: */
	    }

L780:
	;}

/* L800: */
    }
/*     .......... end backsubstitution .......... */
/*     .......... vectors of isolated roots .......... */
    i_2 = *n;
    for (i = 1; i <= i_2; ++i) {
	if (i >= *low && i <= *igh) {
	    goto L840;
	}

	i_1 = *n;
	for (j = i; j <= i_1; ++j) {
	    zr[i + j * zr_dim1] = hr[i + j * hr_dim1];
	    zi[i + j * zi_dim1] = hi[i + j * hi_dim1];
/* L820: */
	}

L840:
    ;}
/*     .......... multiply by transformation matrix to give */
/*                vectors of original full matrix. */
/*                for j=n step -1 until low do -- .......... */
    i_2 = *n;
    for (jj = *low; jj <= i_2; ++jj) {
	j = *n + *low - jj;
	m = min(j,*igh);

	i_1 = *igh;
	for (i = *low; i <= i_1; ++i) {
	    zzr = 0.;
	    zzi = 0.;

	    i_3 = m;
	    for (k = *low; k <= i_3; ++k) {
		zzr = zzr + zr[i + k * zr_dim1] * hr[k + j * hr_dim1] - zi[i 
			+ k * zi_dim1] * hi[k + j * hi_dim1];
		zzi = zzi + zr[i + k * zr_dim1] * hi[k + j * hi_dim1] + zi[i 
			+ k * zi_dim1] * hr[k + j * hr_dim1];
/* L860: */
	    }

	    zr[i + j * zr_dim1] = zzr;
	    zi[i + j * zi_dim1] = zzi;
/* L880: */
	}
    }

    goto L1001;
/*     .......... set error -- all eigenvalues have not */
/*                converged after 30*n iterations .......... */
L1000:
    *ierr = en;
L1001:
    return 0;
} /* comqr2_ */

/* Subroutine */ int comqr_(nm, n, low, igh, hr, hi, wr, wi, ierr)
integer *nm, *n, *low, *igh;
doublereal *hr, *hi, *wr, *wi;
integer *ierr;
{
    /* System generated locals */
    integer hr_dim1, hr_offset, hi_dim1, hi_offset, i_1, i_2;
    doublereal d_1, d_2, d_3, d_4;

    /* Local variables */
    extern /* Subroutine */ int cdiv_();
    static doublereal norm;
    static integer i, j, l, en, ll;
    static doublereal si, ti, xi, yi, sr, tr, xr, yr;
    extern doublereal pythag_();
    extern /* Subroutine */ int csroot_();
    static integer lp1, itn, its;
    static doublereal zzi, zzr;
    static integer enm1;
    static doublereal tst1, tst2;

    /* Parameter adjustments */
    hr_dim1 = *nm;
    hr_offset = hr_dim1 + 1;
    hr -= hr_offset;
    hi_dim1 = *nm;
    hi_offset = hi_dim1 + 1;
    hi -= hi_offset;
    --wr;
    --wi;

    /* Function Body */


/*     this subroutine is a translation of a unitary analogue of the */
/*     algol procedure  comlr, num. math. 12, 369-376(1968) by martin */
/*     and wilkinson. */
/*     handbook for auto. comp., vol.ii-linear algebra, 396-403(1971). */
/*     the unitary analogue substitutes the qr algorithm of francis */
/*     (comp. jour. 4, 332-345(1962)) for the lr algorithm. */

/*     this subroutine finds the eigenvalues of a complex */
/*     upper hessenberg matrix by the qr method. */

/*     on input */

/*        nm must be set to the row dimension of two-dimensional */
/*          array parameters as declared in the calling program */
/*          dimension statement. */

/*        n is the order of the matrix. */

/*        low and igh are integers determined by the balancing */
/*          subroutine  cbal.  if  cbal  has not been used, */
/*          set low=1, igh=n. */

/*        hr and hi contain the real and imaginary parts, */
/*          respectively, of the complex upper hessenberg matrix. */
/*          their lower triangles below the subdiagonal contain */
/*          information about the unitary transformations used in */
/*          the reduction by  corth, if performed. */

/*     on output */

/*        the upper hessenberg portions of hr and hi have been */
/*          destroyed.  therefore, they must be saved before */
/*          calling  comqr  if subsequent calculation of */
/*          eigenvectors is to be performed. */

/*        wr and wi contain the real and imaginary parts, */
/*          respectively, of the eigenvalues.  if an error */
/*          exit is made, the eigenvalues should be correct */
/*          for indices ierr+1,...,n. */

/*        ierr is set to */
/*          zero       for normal return, */
/*          j          if the limit of 30*n iterations is exhausted */
/*                     while the j-th eigenvalue is being sought. */

/*     calls cdiv for complex division. */
/*     calls csroot for complex square root. */
/*     calls pythag for  dsqrt(a*a + b*b) . */

/*     questions and comments should be directed to burton s. garbow, */
/*     mathematics and computer science div, argonne national laboratory 
*/

/*     this version dated august 1983. */

/*     ------------------------------------------------------------------ 
*/

    *ierr = 0;
    if (*low == *igh) {
	goto L180;
    }
/*     .......... create real subdiagonal elements .......... */
    l = *low + 1;

    i_1 = *igh;
    for (i = l; i <= i_1; ++i) {
/* Computing MAX */
	i_2 = i + 1;
	ll = min(*igh,i_2);
	if (hi[i + (i - 1) * hi_dim1] == 0.) {
	    goto L170;
	}
	norm = pythag_(&hr[i + (i - 1) * hr_dim1], &hi[i + (i - 1) * hi_dim1])
		;
	yr = hr[i + (i - 1) * hr_dim1] / norm;
	yi = hi[i + (i - 1) * hi_dim1] / norm;
	hr[i + (i - 1) * hr_dim1] = norm;
	hi[i + (i - 1) * hi_dim1] = 0.;

	i_2 = *igh;
	for (j = i; j <= i_2; ++j) {
	    si = yr * hi[i + j * hi_dim1] - yi * hr[i + j * hr_dim1];
	    hr[i + j * hr_dim1] = yr * hr[i + j * hr_dim1] + yi * hi[i + j * 
		    hi_dim1];
	    hi[i + j * hi_dim1] = si;
/* L155: */
	}

	i_2 = ll;
	for (j = *low; j <= i_2; ++j) {
	    si = yr * hi[j + i * hi_dim1] + yi * hr[j + i * hr_dim1];
	    hr[j + i * hr_dim1] = yr * hr[j + i * hr_dim1] - yi * hi[j + i * 
		    hi_dim1];
	    hi[j + i * hi_dim1] = si;
/* L160: */
	}

L170:
    ;}
/*     .......... store roots isolated by cbal .......... */
L180:
    i_1 = *n;
    for (i = 1; i <= i_1; ++i) {
	if (i >= *low && i <= *igh) {
	    goto L200;
	}
	wr[i] = hr[i + i * hr_dim1];
	wi[i] = hi[i + i * hi_dim1];
L200:
    ;}

    en = *igh;
    tr = 0.;
    ti = 0.;
    itn = *n * 30;
/*     .......... search for next eigenvalue .......... */
L220:
    if (en < *low) {
	goto L1001;
    }
    its = 0;
    enm1 = en - 1;
/*     .......... look for single small sub-diagonal element */
/*                for l=en step -1 until low d0 -- .......... */
L240:
    i_1 = en;
    for (ll = *low; ll <= i_1; ++ll) {
	l = en + *low - ll;
	if (l == *low) {
	    goto L300;
	}
	tst1 = (d_1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[l 
		- 1 + (l - 1) * hi_dim1], abs(d_2)) + (d_3 = hr[l + l * 
		hr_dim1], abs(d_3)) + (d_4 = hi[l + l * hi_dim1], abs(d_4));
	tst2 = tst1 + (d_1 = hr[l + (l - 1) * hr_dim1], abs(d_1));
	if (tst2 == tst1) {
	    goto L300;
	}
/* L260: */
    }
/*     .......... form shift .......... */
L300:
    if (l == en) {
	goto L660;
    }
    if (itn == 0) {
	goto L1000;
    }
    if (its == 10 || its == 20) {
	goto L320;
    }
    sr = hr[en + en * hr_dim1];
    si = hi[en + en * hi_dim1];
    xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1];
    xi = hi[enm1 + en * hi_dim1] * hr[en + enm1 * hr_dim1];
    if (xr == 0. && xi == 0.) {
	goto L340;
    }
    yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
    yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
/* Computing 2nd power */
    d_2 = yr;
/* Computing 2nd power */
    d_3 = yi;
    d_1 = d_2 * d_2 - d_3 * d_3 + xr;
    d_4 = yr * 2. * yi + xi;
    csroot_(&d_1, &d_4, &zzr, &zzi);
    if (yr * zzr + yi * zzi >= 0.) {
	goto L310;
    }
    zzr = -zzr;
    zzi = -zzi;
L310:
    d_1 = yr + zzr;
    d_2 = yi + zzi;
    cdiv_(&xr, &xi, &d_1, &d_2, &xr, &xi);
    sr -= xr;
    si -= xi;
    goto L340;
/*     .......... form exceptional shift .......... */
L320:
    sr = (d_1 = hr[en + enm1 * hr_dim1], abs(d_1)) + (d_2 = hr[enm1 + (en - 2)
	     * hr_dim1], abs(d_2));
    si = 0.;

L340:
    i_1 = en;
    for (i = *low; i <= i_1; ++i) {
	hr[i + i * hr_dim1] -= sr;
	hi[i + i * hi_dim1] -= si;
/* L360: */
    }

    tr += sr;
    ti += si;
    ++its;
    --itn;
/*     .......... reduce to triangle (rows) .......... */
    lp1 = l + 1;

    i_1 = en;
    for (i = lp1; i <= i_1; ++i) {
	sr = hr[i + (i - 1) * hr_dim1];
	hr[i + (i - 1) * hr_dim1] = 0.;
	d_1 = pythag_(&hr[i - 1 + (i - 1) * hr_dim1], &hi[i - 1 + (i - 1) * 
		hi_dim1]);
	norm = pythag_(&d_1, &sr);
	xr = hr[i - 1 + (i - 1) * hr_dim1] / norm;
	wr[i - 1] = xr;
	xi = hi[i - 1 + (i - 1) * hi_dim1] / norm;
	wi[i - 1] = xi;
	hr[i - 1 + (i - 1) * hr_dim1] = norm;
	hi[i - 1 + (i - 1) * hi_dim1] = 0.;
	hi[i + (i - 1) * hi_dim1] = sr / norm;

	i_2 = en;
	for (j = i; j <= i_2; ++j) {
	    yr = hr[i - 1 + j * hr_dim1];
	    yi = hi[i - 1 + j * hi_dim1];
	    zzr = hr[i + j * hr_dim1];
	    zzi = hi[i + j * hi_dim1];
	    hr[i - 1 + j * hr_dim1] = xr * yr + xi * yi + hi[i + (i - 1) * 
		    hi_dim1] * zzr;
	    hi[i - 1 + j * hi_dim1] = xr * yi - xi * yr + hi[i + (i - 1) * 
		    hi_dim1] * zzi;
	    hr[i + j * hr_dim1] = xr * zzr - xi * zzi - hi[i + (i - 1) * 
		    hi_dim1] * yr;
	    hi[i + j * hi_dim1] = xr * zzi + xi * zzr - hi[i + (i - 1) * 
		    hi_dim1] * yi;
/* L490: */
	}

/* L500: */
    }

    si = hi[en + en * hi_dim1];
    if (si == 0.) {
	goto L540;
    }
    norm = pythag_(&hr[en + en * hr_dim1], &si);
    sr = hr[en + en * hr_dim1] / norm;
    si /= norm;
    hr[en + en * hr_dim1] = norm;
    hi[en + en * hi_dim1] = 0.;
/*     .......... inverse operation (columns) .......... */
L540:
    i_1 = en;
    for (j = lp1; j <= i_1; ++j) {
	xr = wr[j - 1];
	xi = wi[j - 1];

	i_2 = j;
	for (i = l; i <= i_2; ++i) {
	    yr = hr[i + (j - 1) * hr_dim1];
	    yi = 0.;
	    zzr = hr[i + j * hr_dim1];
	    zzi = hi[i + j * hi_dim1];
	    if (i == j) {
		goto L560;
	    }
	    yi = hi[i + (j - 1) * hi_dim1];
	    hi[i + (j - 1) * hi_dim1] = xr * yi + xi * yr + hi[j + (j - 1) * 
		    hi_dim1] * zzi;
L560:
	    hr[i + (j - 1) * hr_dim1] = xr * yr - xi * yi + hi[j + (j - 1) * 
		    hi_dim1] * zzr;
	    hr[i + j * hr_dim1] = xr * zzr + xi * zzi - hi[j + (j - 1) * 
		    hi_dim1] * yr;
	    hi[i + j * hi_dim1] = xr * zzi - xi * zzr - hi[j + (j - 1) * 
		    hi_dim1] * yi;
/* L580: */
	}

/* L600: */
    }

    if (si == 0.) {
	goto L240;
    }

    i_1 = en;
    for (i = l; i <= i_1; ++i) {
	yr = hr[i + en * hr_dim1];
	yi = hi[i + en * hi_dim1];
	hr[i + en * hr_dim1] = sr * yr - si * yi;
	hi[i + en * hi_dim1] = sr * yi + si * yr;
/* L630: */
    }

    goto L240;
/*     .......... a root found .......... */
L660:
    wr[en] = hr[en + en * hr_dim1] + tr;
    wi[en] = hi[en + en * hi_dim1] + ti;
    en = enm1;
    goto L220;
/*     .......... set error -- all eigenvalues have not */
/*                converged after 30*n iterations .......... */
L1000:
    *ierr = en;
L1001:
    return 0;
} /* comqr_ */


/* Subroutine */ int csroot_(xr, xi, yr, yi)
doublereal *xr, *xi, *yr, *yi;
{
    /* Builtin functions */
    double sqrt();

    /* Local variables */
    static doublereal s, ti, tr;
    extern doublereal pythag_();


/*     (yr,yi) = complex dsqrt(xr,xi) */
/*     branch chosen so that yr .ge. 0.0 and sign(yi) .eq. sign(xi) */

    tr = *xr;
    ti = *xi;
    s = sqrt((pythag_(&tr, &ti) + abs(tr)) * .5);
    if (tr >= 0.) {
	*yr = s;
    }
    if (ti < 0.) {
	s = -s;
    }
    if (tr <= 0.) {
	*yi = s;
    }
    if (tr < 0.) {
	*yr = ti / *yi * .5;
    }
    if (tr > 0.) {
	*yi = ti / *yr * .5;
    }
    return 0;
} /* csroot_ */

/* Subroutine */ int corth_(nm, n, low, igh, ar, ai, ortr, orti)
integer *nm, *n, *low, *igh;
doublereal *ar, *ai, *ortr, *orti;
{
    /* System generated locals */
    integer ar_dim1, ar_offset, ai_dim1, ai_offset, i_1, i_2, i_3;
    doublereal d_1, d_2;

    /* Builtin functions */
    double sqrt();

    /* Local variables */
    static doublereal f, g, h;
    static integer i, j, m;
    static doublereal scale;
    static integer la;
    static doublereal fi;
    static integer ii, jj;
    static doublereal fr;
    static integer mp;
    extern doublereal pythag_();
    static integer kp1;

    /* Parameter adjustments */
    ar_dim1 = *nm;
    ar_offset = ar_dim1 + 1;
    ar -= ar_offset;
    ai_dim1 = *nm;
    ai_offset = ai_dim1 + 1;
    ai -= ai_offset;
    --ortr;
    --orti;

    /* Function Body */


/*     this subroutine is a translation of a complex analogue of */
/*     the algol procedure orthes, num. math. 12, 349-368(1968) */
/*     by martin and wilkinson. */
/*     handbook for auto. comp., vol.ii-linear algebra, 339-358(1971). */

/*     given a complex general matrix, this subroutine */
/*     reduces a submatrix situated in rows and columns */
/*     low through igh to upper hessenberg form by */
/*     unitary similarity transformations. */

/*     on input */

/*        nm must be set to the row dimension of two-dimensional */
/*          array parameters as declared in the calling program */
/*          dimension statement. */

/*        n is the order of the matrix. */

/*        low and igh are integers determined by the balancing */
/*          subroutine  cbal.  if  cbal  has not been used, */
/*          set low=1, igh=n. */

/*        ar and ai contain the real and imaginary parts, */
/*          respectively, of the complex input matrix. */

/*     on output */

/*        ar and ai contain the real and imaginary parts, */
/*          respectively, of the hessenberg matrix.  information */
/*          about the unitary transformations used in the reduction */
/*          is stored in the remaining triangles under the */
/*          hessenberg matrix. */

/*        ortr and orti contain further information about the */
/*          transformations.  only elements low through igh are used. */

/*     calls pythag for  dsqrt(a*a + b*b) . */

/*     questions and comments should be directed to burton s. garbow, */
/*     mathematics and computer science div, argonne national laboratory 
*/

/*     this version dated august 1983. */

/*     ------------------------------------------------------------------ 
*/

    la = *igh - 1;
    kp1 = *low + 1;
    if (la < kp1) {
	goto L200;
    }

    i_1 = la;
    for (m = kp1; m <= i_1; ++m) {
	h = 0.;
	ortr[m] = 0.;
	orti[m] = 0.;
	scale = 0.;
/*     .......... scale column (algol tol then not needed) .......... 
*/
	i_2 = *igh;
	for (i = m; i <= i_2; ++i) {
/* L90: */
	    scale = scale + (d_1 = ar[i + (m - 1) * ar_dim1], abs(d_1)) + (
		    d_2 = ai[i + (m - 1) * ai_dim1], abs(d_2));
	}

	if (scale == 0.) {
	    goto L180;
	}
	mp = m + *igh;
/*     .......... for i=igh step -1 until m do -- .......... */
	i_2 = *igh;
	for (ii = m; ii <= i_2; ++ii) {
	    i = mp - ii;
	    ortr[i] = ar[i + (m - 1) * ar_dim1] / scale;
	    orti[i] = ai[i + (m - 1) * ai_dim1] / scale;
	    h = h + ortr[i] * ortr[i] + orti[i] * orti[i];
/* L100: */
	}

	g = sqrt(h);
	f = pythag_(&ortr[m], &orti[m]);
	if (f == 0.) {
	    goto L103;
	}
	h += f * g;
	g /= f;
	ortr[m] = (g + 1.) * ortr[m];
	orti[m] = (g + 1.) * orti[m];
	goto L105;

L103:
	ortr[m] = g;
	ar[m + (m - 1) * ar_dim1] = scale;
/*     .......... form (i-(u*ut)/h) * a .......... */
L105:
	i_2 = *n;
	for (j = m; j <= i_2; ++j) {
	    fr = 0.;
	    fi = 0.;
/*     .......... for i=igh step -1 until m do -- .......... */
	    i_3 = *igh;
	    for (ii = m; ii <= i_3; ++ii) {
		i = mp - ii;
		fr = fr + ortr[i] * ar[i + j * ar_dim1] + orti[i] * ai[i + j *
			 ai_dim1];
		fi = fi + ortr[i] * ai[i + j * ai_dim1] - orti[i] * ar[i + j *
			 ar_dim1];
/* L110: */
	    }

	    fr /= h;
	    fi /= h;

	    i_3 = *igh;
	    for (i = m; i <= i_3; ++i) {
		ar[i + j * ar_dim1] = ar[i + j * ar_dim1] - fr * ortr[i] + fi 
			* orti[i];
		ai[i + j * ai_dim1] = ai[i + j * ai_dim1] - fr * orti[i] - fi 
			* ortr[i];
/* L120: */
	    }

/* L130: */
	}
/*     .......... form (i-(u*ut)/h)*a*(i-(u*ut)/h) .......... */
	i_2 = *igh;
	for (i = 1; i <= i_2; ++i) {
	    fr = 0.;
	    fi = 0.;
/*     .......... for j=igh step -1 until m do -- .......... */
	    i_3 = *igh;
	    for (jj = m; jj <= i_3; ++jj) {
		j = mp - jj;
		fr = fr + ortr[j] * ar[i + j * ar_dim1] - orti[j] * ai[i + j *
			 ai_dim1];
		fi = fi + ortr[j] * ai[i + j * ai_dim1] + orti[j] * ar[i + j *
			 ar_dim1];
/* L140: */
	    }

	    fr /= h;
	    fi /= h;

	    i_3 = *igh;
	    for (j = m; j <= i_3; ++j) {
		ar[i + j * ar_dim1] = ar[i + j * ar_dim1] - fr * ortr[j] - fi 
			* orti[j];
		ai[i + j * ai_dim1] = ai[i + j * ai_dim1] + fr * orti[j] - fi 
			* ortr[j];
/* L150: */
	    }

/* L160: */
	}

	ortr[m] = scale * ortr[m];
	orti[m] = scale * orti[m];
	ar[m + (m - 1) * ar_dim1] = -g * ar[m + (m - 1) * ar_dim1];
	ai[m + (m - 1) * ai_dim1] = -g * ai[m + (m - 1) * ai_dim1];
L180:
    ;}

L200:
    return 0;
} /* corth_ */

doublereal pythag_(a, b)
doublereal *a, *b;
{
    /* System generated locals */
    doublereal ret_val, d_1, d_2, d_3;

    /* Local variables */
    static doublereal p, r, s, t, u;


/*     finds dsqrt(a**2+b**2) without overflow or destructive underflow */


/* Computing MAX */
    d_1 = abs(*a), d_2 = abs(*b);
    p = max(d_2,d_1);
    if (p == 0.) {
	goto L20;
    }
/* Computing MAX */
    d_2 = abs(*a), d_3 = abs(*b);
/* Computing 2nd power */
    d_1 = min(d_3,d_2) / p;
    r = d_1 * d_1;
L10:
    t = r + 4.;
    if (t == 4.) {
	goto L20;
    }
    s = r / t;
    u = s * 2. + 1.;
    p = u * p;
/* Computing 2nd power */
    d_1 = s / u;
    r = d_1 * d_1 * r;
    goto L10;
L20:
    ret_val = p;
    return ret_val;
} /* pythag_ */

/* Subroutine */ int cbal_(nm, n, ar, ai, low, igh, scale)
integer *nm, *n;
doublereal *ar, *ai;
integer *low, *igh;
doublereal *scale;
{
    /* System generated locals */
    integer ar_dim1, ar_offset, ai_dim1, ai_offset, i_1, i_2;
    doublereal d_1, d_2;

    /* Local variables */
    static integer iexc;
    static doublereal c, f, g;
    static integer i, j, k, l, m;
    static doublereal r, s, radix, b2;
    static integer jj;
    static logical noconv;

    /* Parameter adjustments */
    ar_dim1 = *nm;
    ar_offset = ar_dim1 + 1;
    ar -= ar_offset;
    ai_dim1 = *nm;
    ai_offset = ai_dim1 + 1;
    ai -= ai_offset;
    --scale;

    /* Function Body */


/*     this subroutine is a translation of the algol procedure */
/*     cbalance, which is a complex version of balance, */
/*     num. math. 13, 293-304(1969) by parlett and reinsch. */
/*     handbook for auto. comp., vol.ii-linear algebra, 315-326(1971). */

/*     this subroutine balances a complex matrix and isolates */
/*     eigenvalues whenever possible. */

/*     on input */

/*        nm must be set to the row dimension of two-dimensional */
/*          array parameters as declared in the calling program */
/*          dimension statement. */

/*        n is the order of the matrix. */

/*        ar and ai contain the real and imaginary parts, */
/*          respectively, of the complex matrix to be balanced. */

/*     on output */

/*        ar and ai contain the real and imaginary parts, */
/*          respectively, of the balanced matrix. */

/*        low and igh are two integers such that ar(i,j) and ai(i,j) */
/*          are equal to zero if */
/*           (1) i is greater than j and */
/*           (2) j=1,...,low-1 or i=igh+1,...,n. */

/*        scale contains information determining the */
/*           permutations and scaling factors used. */

/*     suppose that the principal submatrix in rows low through igh */
/*     has been balanced, that p(j) denotes the index interchanged */
/*     with j during the permutation step, and that the elements */
/*     of the diagonal matrix used are denoted by d(i,j).  then */
/*        scale(j) = p(j),    for j = 1,...,low-1 */
/*                 = d(j,j)       j = low,...,igh */
/*                 = p(j)         j = igh+1,...,n. */
/*     the order in which the interchanges are made is n to igh+1, */
/*     then 1 to low-1. */

/*     note that 1 is returned for igh if igh is zero formally. */

/*     the algol procedure exc contained in cbalance appears in */
/*     cbal  in line.  (note that the algol roles of identifiers */
/*     k,l have been reversed.) */

/*     arithmetic is real throughout. */

/*     questions and comments should be directed to burton s. garbow, */
/*     mathematics and computer science div, argonne national laboratory 
*/

/*     this version dated august 1983. */

/*     ------------------------------------------------------------------ 
*/

    radix = 16.;

    b2 = radix * radix;
    k = 1;
    l = *n;
    goto L100;
/*     .......... in-line procedure for row and */
/*                column exchange .......... */
L20:
    scale[m] = (doublereal) j;
    if (j == m) {
	goto L50;
    }

    i_1 = l;
    for (i = 1; i <= i_1; ++i) {
	f = ar[i + j * ar_dim1];
	ar[i + j * ar_dim1] = ar[i + m * ar_dim1];
	ar[i + m * ar_dim1] = f;
	f = ai[i + j * ai_dim1];
	ai[i + j * ai_dim1] = ai[i + m * ai_dim1];
	ai[i + m * ai_dim1] = f;
/* L30: */
    }

    i_1 = *n;
    for (i = k; i <= i_1; ++i) {
	f = ar[j + i * ar_dim1];
	ar[j + i * ar_dim1] = ar[m + i * ar_dim1];
	ar[m + i * ar_dim1] = f;
	f = ai[j + i * ai_dim1];
	ai[j + i * ai_dim1] = ai[m + i * ai_dim1];
	ai[m + i * ai_dim1] = f;
/* L40: */
    }

L50:
    switch (iexc) {
	case 1:  goto L80;
	case 2:  goto L130;
    }
/*     .......... search for rows isolating an eigenvalue */
/*                and push them down .......... */
L80:
    if (l == 1) {
	goto L280;
    }
    --l;
/*     .......... for j=l step -1 until 1 do -- .......... */
L100:
    i_1 = l;
    for (jj = 1; jj <= i_1; ++jj) {
	j = l + 1 - jj;

	i_2 = l;
	for (i = 1; i <= i_2; ++i) {
	    if (i == j) {
		goto L110;
	    }
	    if (ar[j + i * ar_dim1] != 0. || ai[j + i * ai_dim1] != 0.) {
		goto L120;
	    }
L110:
	;}

	m = l;
	iexc = 1;
	goto L20;
L120:
    ;}

    goto L140;
/*     .......... search for columns isolating an eigenvalue */
/*                and push them left .......... */
L130:
    ++k;

L140:
    i_1 = l;
    for (j = k; j <= i_1; ++j) {

	i_2 = l;
	for (i = k; i <= i_2; ++i) {
	    if (i == j) {
		goto L150;
	    }
	    if (ar[i + j * ar_dim1] != 0. || ai[i + j * ai_dim1] != 0.) {
		goto L170;
	    }
L150:
	;}

	m = k;
	iexc = 2;
	goto L20;
L170:
    ;}
/*     .......... now balance the submatrix in rows k to l .......... */
    i_1 = l;
    for (i = k; i <= i_1; ++i) {
/* L180: */
	scale[i] = 1.;
    }
/*     .......... iterative loop for norm reduction .......... */
L190:
    noconv = FALSE_;

    i_1 = l;
    for (i = k; i <= i_1; ++i) {
	c = 0.;
	r = 0.;

	i_2 = l;
	for (j = k; j <= i_2; ++j) {
	    if (j == i) {
		goto L200;
	    }
	    c = c + (d_1 = ar[j + i * ar_dim1], abs(d_1)) + (d_2 = ai[j + i * 
		    ai_dim1], abs(d_2));
	    r = r + (d_1 = ar[i + j * ar_dim1], abs(d_1)) + (d_2 = ai[i + j * 
		    ai_dim1], abs(d_2));
L200:
	;}
/*     .......... guard against zero c or r due to underflow 
.......... */
	if (c == 0. || r == 0.) {
	    goto L270;
	}
	g = r / radix;
	f = 1.;
	s = c + r;
L210:
	if (c >= g) {
	    goto L220;
	}
	f *= radix;
	c *= b2;
	goto L210;
L220:
	g = r * radix;
L230:
	if (c < g) {
	    goto L240;
	}
	f /= radix;
	c /= b2;
	goto L230;
/*     .......... now balance .......... */
L240:
	if ((c + r) / f >= s * .95) {
	    goto L270;
	}
	g = 1. / f;
	scale[i] *= f;
	noconv = TRUE_;

	i_2 = *n;
	for (j = k; j <= i_2; ++j) {
	    ar[i + j * ar_dim1] *= g;
	    ai[i + j * ai_dim1] *= g;
/* L250: */
	}

	i_2 = l;
	for (j = 1; j <= i_2; ++j) {
	    ar[j + i * ar_dim1] *= f;
	    ai[j + i * ai_dim1] *= f;
/* L260: */
	}

L270:
    ;}

    if (noconv) {
	goto L190;
    }

L280:
    *low = k;
    *igh = l;
    return 0;
} /* cbal_ */

End of cg.c
echo rg.c 1>&2
cat >rg.c <<'End of rg.c'
/* rg.f -- translated by f2c (version of 3 February 1990  3:36:42).
   from the eispack library
*/

#include "f2c.h"

/* Table of constant values */

static doublereal c_b60 = 0.;

/* Subroutine */ int rg_(nm, n, a, wr, wi, matz, z, iv1, fv1, ierr)
integer *nm, *n;
doublereal *a, *wr, *wi;
integer *matz;
doublereal *z;
integer *iv1;
doublereal *fv1;
integer *ierr;
{
    /* System generated locals */
    integer a_dim1, a_offset, z_dim1, z_offset;

    /* Local variables */
    extern /* Subroutine */ int balbak_(), balanc_(), elmhes_(), eltran_();
    static integer is1, is2;
    extern /* Subroutine */ int hqr_(), hqr2_();

    /* Parameter adjustments */
    a_dim1 = *nm;
    a_offset = a_dim1 + 1;
    a -= a_offset;
    --wr;
    --wi;
    z_dim1 = *nm;
    z_offset = z_dim1 + 1;
    z -= z_offset;
    --iv1;
    --fv1;

    /* Function Body */


/*     this subroutine calls the recommended sequence of */
/*     subroutines from the eigensystem subroutine package (eispack) */
/*     to find the eigenvalues and eigenvectors (if desired) */
/*     of a real general matrix. */

/*     on input */

/*        nm  must be set to the row dimension of the two-dimensional */
/*        array parameters as declared in the calling program */
/*        dimension statement. */

/*        n  is the order of the matrix  a. */

/*        a  contains the real general matrix. */

/*        matz  is an integer variable set equal to zero if */
/*        only eigenvalues are desired.  otherwise it is set to */
/*        any non-zero integer for both eigenvalues and eigenvectors. */

/*     on output */

/*        wr  and  wi  contain the real and imaginary parts, */
/*        respectively, of the eigenvalues.  complex conjugate */
/*        pairs of eigenvalues appear consecutively with the */
/*        eigenvalue having the positive imaginary part first. */

/*        z  contains the real and imaginary parts of the eigenvectors */
/*        if matz is not zero.  if the j-th eigenvalue is real, the */
/*        j-th column of  z  contains its eigenvector.  if the j-th */
/*        eigenvalue is complex with positive imaginary part, the */
/*        j-th and (j+1)-th columns of  z  contain the real and */
/*        imaginary parts of its eigenvector.  the conjugate of this */
/*        vector is the eigenvector for the conjugate eigenvalue. */

/*        ierr  is an integer output variable set equal to an error */
/*           completion code described in the documentation for hqr */
/*           and hqr2.  the normal completion code is zero. */

/*        iv1  and  fv1  are temporary storage arrays. */

/*     questions and comments should be directed to burton s. garbow, */
/*     mathematics and computer science div, argonne national laboratory 
*/

/*     this version dated august 1983. */

/*     ------------------------------------------------------------------ 
*/

    if (*n <= *nm) {
	goto L10;
    }
    *ierr = *n * 10;
    goto L50;

L10:
    balanc_(nm, n, &a[a_offset], &is1, &is2, &fv1[1]);
    elmhes_(nm, n, &is1, &is2, &a[a_offset], &iv1[1]);
    if (*matz != 0) {
	goto L20;
    }
/*     .......... find eigenvalues only .......... */
    hqr_(nm, n, &is1, &is2, &a[a_offset], &wr[1], &wi[1], ierr);
    goto L50;
/*     .......... find both eigenvalues and eigenvectors .......... */
L20:
    eltran_(nm, n, &is1, &is2, &a[a_offset], &iv1[1], &z[z_offset]);
    hqr2_(nm, n, &is1, &is2, &a[a_offset], &wr[1], &wi[1], &z[z_offset], ierr)
	    ;
    if (*ierr != 0) {
	goto L50;
    }
    balbak_(nm, n, &is1, &is2, &fv1[1], n, &z[z_offset]);
L50:
    return 0;
} /* rg_ */

/* Subroutine */ int balbak_(nm, n, low, igh, scale, m, z)
integer *nm, *n, *low, *igh;
doublereal *scale;
integer *m;
doublereal *z;
{
    /* System generated locals */
    integer z_dim1, z_offset, i_1, i_2;

    /* Local variables */
    static integer i, j, k;
    static doublereal s;
    static integer ii;

    /* Parameter adjustments */
    --scale;
    z_dim1 = *nm;
    z_offset = z_dim1 + 1;
    z -= z_offset;

    /* Function Body */


/*     this subroutine is a translation of the algol procedure balbak, */
/*     num. math. 13, 293-304(1969) by parlett and reinsch. */
/*     handbook for auto. comp., vol.ii-linear algebra, 315-326(1971). */

/*     this subroutine forms the eigenvectors of a real general */
/*     matrix by back transforming those of the corresponding */
/*     balanced matrix determined by  balanc. */

/*     on input */

/*        nm must be set to the row dimension of two-dimensional */
/*          array parameters as declared in the calling program */
/*          dimension statement. */

/*        n is the order of the matrix. */

/*        low and igh are integers determined by  balanc. */

/*        scale contains information determining the permutations */
/*          and scaling factors used by  balanc. */

/*        m is the number of columns of z to be back transformed. */

/*        z contains the real and imaginary parts of the eigen- */
/*          vectors to be back transformed in its first m columns. */

/*     on output */

/*        z contains the real and imaginary parts of the */
/*          transformed eigenvectors in its first m columns. */

/*     questions and comments should be directed to burton s. garbow, */
/*     mathematics and computer science div, argonne national laboratory 
*/

/*     this version dated august 1983. */

/*     ------------------------------------------------------------------ 
*/

    if (*m == 0) {
	goto L200;
    }
    if (*igh == *low) {
	goto L120;
    }

    i_1 = *igh;
    for (i = *low; i <= i_1; ++i) {
	s = scale[i];
/*     .......... left hand eigenvectors are back transformed */
/*                if the foregoing statement is replaced by */
/*                s=1.0d0/scale(i). .......... */
	i_2 = *m;
	for (j = 1; j <= i_2; ++j) {
/* L100: */
	    z[i + j * z_dim1] *= s;
	}

/* L110: */
    }
/*     ......... for i=low-1 step -1 until 1, */
/*               igh+1 step 1 until n do -- .......... */
L120:
    i_1 = *n;
    for (ii = 1; ii <= i_1; ++ii) {
	i = ii;
	if (i >= *low && i <= *igh) {
	    goto L140;
	}
	if (i < *low) {
	    i = *low - ii;
	}
	k = (integer) scale[i];
	if (k == i) {
	    goto L140;
	}

	i_2 = *m;
	for (j = 1; j <= i_2; ++j) {
	    s = z[i + j * z_dim1];
	    z[i + j * z_dim1] = z[k + j * z_dim1];
	    z[k + j * z_dim1] = s;
/* L130: */
	}

L140:
    ;}

L200:
    return 0;
} /* balbak_ */

/* Subroutine */ int hqr2_(nm, n, low, igh, h, wr, wi, z, ierr)
integer *nm, *n, *low, *igh;
doublereal *h, *wr, *wi, *z;
integer *ierr;
{
    /* System generated locals */
    integer h_dim1, h_offset, z_dim1, z_offset, i_1, i_2, i_3;
    doublereal d_1, d_2, d_3, d_4;

    /* Builtin functions */
    double sqrt(), d_sign();

    /* Local variables */
    extern /* Subroutine */ int cdiv_();
    static doublereal norm;
    static integer i, j, k, l, m;
    static doublereal p, q, r, s, t, w, x, y;
    static integer na, ii, en, jj;
    static doublereal ra, sa;
    static integer ll, mm, nn;
    static doublereal vi, vr, zz;
    static logical notlas;
    static integer mp2, itn, its, enm2;
    static doublereal tst1, tst2;

    /* Parameter adjustments */
    h_dim1 = *nm;
    h_offset = h_dim1 + 1;
    h -= h_offset;
    --wr;
    --wi;
    z_dim1 = *nm;
    z_offset = z_dim1 + 1;
    z -= z_offset;

    /* Function Body */


/*     this subroutine is a translation of the algol procedure hqr2, */
/*     num. math. 16, 181-204(1970) by peters and wilkinson. */
/*     handbook for auto. comp., vol.ii-linear algebra, 372-395(1971). */

/*     this subroutine finds the eigenvalues and eigenvectors */
/*     of a real upper hessenberg matrix by the qr method.  the */
/*     eigenvectors of a real general matrix can also be found */
/*     if  elmhes  and  eltran  or  orthes  and  ortran  have */
/*     been used to reduce this general matrix to hessenberg form */
/*     and to accumulate the similarity transformations. */

/*     on input */

/*        nm must be set to the row dimension of two-dimensional */
/*          array parameters as declared in the calling program */
/*          dimension statement. */

/*        n is the order of the matrix. */

/*        low and igh are integers determined by the balancing */
/*          subroutine  balanc.  if  balanc  has not been used, */
/*          set low=1, igh=n. */

/*        h contains the upper hessenberg matrix. */

/*        z contains the transformation matrix produced by  eltran */
/*          after the reduction by  elmhes, or by  ortran  after the */
/*          reduction by  orthes, if performed.  if the eigenvectors */
/*          of the hessenberg matrix are desired, z must contain the */
/*          identity matrix. */

/*     on output */

/*        h has been destroyed. */

/*        wr and wi contain the real and imaginary parts, */
/*          respectively, of the eigenvalues.  the eigenvalues */
/*          are unordered except that complex conjugate pairs */
/*          of values appear consecutively with the eigenvalue */
/*          having the positive imaginary part first.  if an */
/*          error exit is made, the eigenvalues should be correct */
/*          for indices ierr+1,...,n. */

/*        z contains the real and imaginary parts of the eigenvectors. */
/*          if the i-th eigenvalue is real, the i-th column of z */
/*          contains its eigenvector.  if the i-th eigenvalue is complex 
*/
/*          with positive imaginary part, the i-th and (i+1)-th */
/*          columns of z contain the real and imaginary parts of its */
/*          eigenvector.  the eigenvectors are unnormalized.  if an */
/*          error exit is made, none of the eigenvectors has been found. 
*/

/*        ierr is set to */
/*          zero       for normal return, */
/*          j          if the limit of 30*n iterations is exhausted */
/*                     while the j-th eigenvalue is being sought. */

/*     calls cdiv for complex division. */

/*     questions and comments should be directed to burton s. garbow, */
/*     mathematics and computer science div, argonne national laboratory 
*/

/*     this version dated august 1983. */

/*     ------------------------------------------------------------------ 
*/

    *ierr = 0;
    norm = 0.;
    k = 1;
/*     .......... store roots isolated by balanc */
/*                and compute matrix norm .......... */
    i_1 = *n;
    for (i = 1; i <= i_1; ++i) {

	i_2 = *n;
	for (j = k; j <= i_2; ++j) {
/* L40: */
	    norm += (d_1 = h[i + j * h_dim1], abs(d_1));
	}

	k = i;
	if (i >= *low && i <= *igh) {
	    goto L50;
	}
	wr[i] = h[i + i * h_dim1];
	wi[i] = 0.;
L50:
    ;}

    en = *igh;
    t = 0.;
    itn = *n * 30;
/*     .......... search for next eigenvalues .......... */
L60:
    if (en < *low) {
	goto L340;
    }
    its = 0;
    na = en - 1;
    enm2 = na - 1;
/*     .......... look for single small sub-diagonal element */
/*                for l=en step -1 until low do -- .......... */
L70:
    i_1 = en;
    for (ll = *low; ll <= i_1; ++ll) {
	l = en + *low - ll;
	if (l == *low) {
	    goto L100;
	}
	s = (d_1 = h[l - 1 + (l - 1) * h_dim1], abs(d_1)) + (d_2 = h[l + l * 
		h_dim1], abs(d_2));
	if (s == 0.) {
	    s = norm;
	}
	tst1 = s;
	tst2 = tst1 + (d_1 = h[l + (l - 1) * h_dim1], abs(d_1));
	if (tst2 == tst1) {
	    goto L100;
	}
/* L80: */
    }
/*     .......... form shift .......... */
L100:
    x = h[en + en * h_dim1];
    if (l == en) {
	goto L270;
    }
    y = h[na + na * h_dim1];
    w = h[en + na * h_dim1] * h[na + en * h_dim1];
    if (l == na) {
	goto L280;
    }
    if (itn == 0) {
	goto L1000;
    }
    if (its != 10 && its != 20) {
	goto L130;
    }
/*     .......... form exceptional shift .......... */
    t += x;

    i_1 = en;
    for (i = *low; i <= i_1; ++i) {
/* L120: */
	h[i + i * h_dim1] -= x;
    }

    s = (d_1 = h[en + na * h_dim1], abs(d_1)) + (d_2 = h[na + enm2 * h_dim1], 
	    abs(d_2));
    x = s * .75;
    y = x;
    w = s * -.4375 * s;
L130:
    ++its;
    --itn;
/*     .......... look for two consecutive small */
/*                sub-diagonal elements. */
/*                for m=en-2 step -1 until l do -- .......... */
    i_1 = enm2;
    for (mm = l; mm <= i_1; ++mm) {
	m = enm2 + l - mm;
	zz = h[m + m * h_dim1];
	r = x - zz;
	s = y - zz;
	p = (r * s - w) / h[m + 1 + m * h_dim1] + h[m + (m + 1) * h_dim1];
	q = h[m + 1 + (m + 1) * h_dim1] - zz - r - s;
	r = h[m + 2 + (m + 1) * h_dim1];
	s = abs(p) + abs(q) + abs(r);
	p /= s;
	q /= s;
	r /= s;
	if (m == l) {
	    goto L150;
	}
	tst1 = abs(p) * ((d_1 = h[m - 1 + (m - 1) * h_dim1], abs(d_1)) + abs(
		zz) + (d_2 = h[m + 1 + (m + 1) * h_dim1], abs(d_2)));
	tst2 = tst1 + (d_1 = h[m + (m - 1) * h_dim1], abs(d_1)) * (abs(q) + 
		abs(r));
	if (tst2 == tst1) {
	    goto L150;
	}
/* L140: */
    }

L150:
    mp2 = m + 2;

    i_1 = en;
    for (i = mp2; i <= i_1; ++i) {
	h[i + (i - 2) * h_dim1] = 0.;
	if (i == mp2) {
	    goto L160;
	}
	h[i + (i - 3) * h_dim1] = 0.;
L160:
    ;}
/*     .......... double qr step involving rows l to en and */
/*                columns m to en .......... */
    i_1 = na;
    for (k = m; k <= i_1; ++k) {
	notlas = k != na;
	if (k == m) {
	    goto L170;
	}
	p = h[k + (k - 1) * h_dim1];
	q = h[k + 1 + (k - 1) * h_dim1];
	r = 0.;
	if (notlas) {
	    r = h[k + 2 + (k - 1) * h_dim1];
	}
	x = abs(p) + abs(q) + abs(r);
	if (x == 0.) {
	    goto L260;
	}
	p /= x;
	q /= x;
	r /= x;
L170:
	d_1 = sqrt(p * p + q * q + r * r);
	s = d_sign(&d_1, &p);
	if (k == m) {
	    goto L180;
	}
	h[k + (k - 1) * h_dim1] = -s * x;
	goto L190;
L180:
	if (l != m) {
	    h[k + (k - 1) * h_dim1] = -h[k + (k - 1) * h_dim1];
	}
L190:
	p += s;
	x = p / s;
	y = q / s;
	zz = r / s;
	q /= p;
	r /= p;
	if (notlas) {
	    goto L225;
	}
/*     .......... row modification .......... */
	i_2 = *n;
	for (j = k; j <= i_2; ++j) {
	    p = h[k + j * h_dim1] + q * h[k + 1 + j * h_dim1];
	    h[k + j * h_dim1] -= p * x;
	    h[k + 1 + j * h_dim1] -= p * y;
/* L200: */
	}

/* Computing MAX */
	i_2 = en, i_3 = k + 3;
	j = min(i_3,i_2);
/*     .......... column modification .......... */
	i_2 = j;
	for (i = 1; i <= i_2; ++i) {
	    p = x * h[i + k * h_dim1] + y * h[i + (k + 1) * h_dim1];
	    h[i + k * h_dim1] -= p;
	    h[i + (k + 1) * h_dim1] -= p * q;
/* L210: */
	}
/*     .......... accumulate transformations .......... */
	i_2 = *igh;
	for (i = *low; i <= i_2; ++i) {
	    p = x * z[i + k * z_dim1] + y * z[i + (k + 1) * z_dim1];
	    z[i + k * z_dim1] -= p;
	    z[i + (k + 1) * z_dim1] -= p * q;
/* L220: */
	}
	goto L255;
L225:
/*     .......... row modification .......... */
	i_2 = *n;
	for (j = k; j <= i_2; ++j) {
	    p = h[k + j * h_dim1] + q * h[k + 1 + j * h_dim1] + r * h[k + 2 + 
		    j * h_dim1];
	    h[k + j * h_dim1] -= p * x;
	    h[k + 1 + j * h_dim1] -= p * y;
	    h[k + 2 + j * h_dim1] -= p * zz;
/* L230: */
	}

/* Computing MAX */
	i_2 = en, i_3 = k + 3;
	j = min(i_3,i_2);
/*     .......... column modification .......... */
	i_2 = j;
	for (i = 1; i <= i_2; ++i) {
	    p = x * h[i + k * h_dim1] + y * h[i + (k + 1) * h_dim1] + zz * h[
		    i + (k + 2) * h_dim1];
	    h[i + k * h_dim1] -= p;
	    h[i + (k + 1) * h_dim1] -= p * q;
	    h[i + (k + 2) * h_dim1] -= p * r;
/* L240: */
	}
/*     .......... accumulate transformations .......... */
	i_2 = *igh;
	for (i = *low; i <= i_2; ++i) {
	    p = x * z[i + k * z_dim1] + y * z[i + (k + 1) * z_dim1] + zz * z[
		    i + (k + 2) * z_dim1];
	    z[i + k * z_dim1] -= p;
	    z[i + (k + 1) * z_dim1] -= p * q;
	    z[i + (k + 2) * z_dim1] -= p * r;
/* L250: */
	}
L255:

L260:
    ;}

    goto L70;
/*     .......... one root found .......... */
L270:
    h[en + en * h_dim1] = x + t;
    wr[en] = h[en + en * h_dim1];
    wi[en] = 0.;
    en = na;
    goto L60;
/*     .......... two roots found .......... */
L280:
    p = (y - x) / 2.;
    q = p * p + w;
    zz = sqrt((abs(q)));
    h[en + en * h_dim1] = x + t;
    x = h[en + en * h_dim1];
    h[na + na * h_dim1] = y + t;
    if (q < 0.) {
	goto L320;
    }
/*     .......... real pair .......... */
    zz = p + d_sign(&zz, &p);
    wr[na] = x + zz;
    wr[en] = wr[na];
    if (zz != 0.) {
	wr[en] = x - w / zz;
    }
    wi[na] = 0.;
    wi[en] = 0.;
    x = h[en + na * h_dim1];
    s = abs(x) + abs(zz);
    p = x / s;
    q = zz / s;
    r = sqrt(p * p + q * q);
    p /= r;
    q /= r;
/*     .......... row modification .......... */
    i_1 = *n;
    for (j = na; j <= i_1; ++j) {
	zz = h[na + j * h_dim1];
	h[na + j * h_dim1] = q * zz + p * h[en + j * h_dim1];
	h[en + j * h_dim1] = q * h[en + j * h_dim1] - p * zz;
/* L290: */
    }
/*     .......... column modification .......... */
    i_1 = en;
    for (i = 1; i <= i_1; ++i) {
	zz = h[i + na * h_dim1];
	h[i + na * h_dim1] = q * zz + p * h[i + en * h_dim1];
	h[i + en * h_dim1] = q * h[i + en * h_dim1] - p * zz;
/* L300: */
    }
/*     .......... accumulate transformations .......... */
    i_1 = *igh;
    for (i = *low; i <= i_1; ++i) {
	zz = z[i + na * z_dim1];
	z[i + na * z_dim1] = q * zz + p * z[i + en * z_dim1];
	z[i + en * z_dim1] = q * z[i + en * z_dim1] - p * zz;
/* L310: */
    }

    goto L330;
/*     .......... complex pair .......... */
L320:
    wr[na] = x + p;
    wr[en] = x + p;
    wi[na] = zz;
    wi[en] = -zz;
L330:
    en = enm2;
    goto L60;
/*     .......... all roots found.  backsubstitute to find */
/*                vectors of upper triangular form .......... */
L340:
    if (norm == 0.) {
	goto L1001;
    }
/*     .......... for en=n step -1 until 1 do -- .......... */
    i_1 = *n;
    for (nn = 1; nn <= i_1; ++nn) {
	en = *n + 1 - nn;
	p = wr[en];
	q = wi[en];
	na = en - 1;
	if (q < 0.) {
	    goto L710;
	} else if (q == 0) {
	    goto L600;
	} else {
	    goto L800;
	}
/*     .......... real vector .......... */
L600:
	m = en;
	h[en + en * h_dim1] = 1.;
	if (na == 0) {
	    goto L800;
	}
/*     .......... for i=en-1 step -1 until 1 do -- .......... */
	i_2 = na;
	for (ii = 1; ii <= i_2; ++ii) {
	    i = en - ii;
	    w = h[i + i * h_dim1] - p;
	    r = 0.;

	    i_3 = en;
	    for (j = m; j <= i_3; ++j) {
/* L610: */
		r += h[i + j * h_dim1] * h[j + en * h_dim1];
	    }

	    if (wi[i] >= 0.) {
		goto L630;
	    }
	    zz = w;
	    s = r;
	    goto L700;
L630:
	    m = i;
	    if (wi[i] != 0.) {
		goto L640;
	    }
	    t = w;
	    if (t != 0.) {
		goto L635;
	    }
	    tst1 = norm;
	    t = tst1;
L632:
	    t *= .01;
	    tst2 = norm + t;
	    if (tst2 > tst1) {
		goto L632;
	    }
L635:
	    h[i + en * h_dim1] = -r / t;
	    goto L680;
/*     .......... solve real equations .......... */
L640:
	    x = h[i + (i + 1) * h_dim1];
	    y = h[i + 1 + i * h_dim1];
	    q = (wr[i] - p) * (wr[i] - p) + wi[i] * wi[i];
	    t = (x * s - zz * r) / q;
	    h[i + en * h_dim1] = t;
	    if (abs(x) <= abs(zz)) {
		goto L650;
	    }
	    h[i + 1 + en * h_dim1] = (-r - w * t) / x;
	    goto L680;
L650:
	    h[i + 1 + en * h_dim1] = (-s - y * t) / zz;

/*     .......... overflow control .......... */
L680:
	    t = (d_1 = h[i + en * h_dim1], abs(d_1));
	    if (t == 0.) {
		goto L700;
	    }
	    tst1 = t;
	    tst2 = tst1 + 1. / tst1;
	    if (tst2 > tst1) {
		goto L700;
	    }
	    i_3 = en;
	    for (j = i; j <= i_3; ++j) {
		h[j + en * h_dim1] /= t;
/* L690: */
	    }

L700:
	;}
/*     .......... end real vector .......... */
	goto L800;
/*     .......... complex vector .......... */
L710:
	m = na;
/*     .......... last vector component chosen imaginary so that */
/*                eigenvector matrix is triangular .......... */
	if ((d_1 = h[en + na * h_dim1], abs(d_1)) <= (d_2 = h[na + en * 
		h_dim1], abs(d_2))) {
	    goto L720;
	}
	h[na + na * h_dim1] = q / h[en + na * h_dim1];
	h[na + en * h_dim1] = -(h[en + en * h_dim1] - p) / h[en + na * h_dim1]
		;
	goto L730;
L720:
	d_1 = -h[na + en * h_dim1];
	d_2 = h[na + na * h_dim1] - p;
	cdiv_(&c_b60, &d_1, &d_2, &q, &h[na + na * h_dim1], &h[na + en * 
		h_dim1]);
L730:
	h[en + na * h_dim1] = 0.;
	h[en + en * h_dim1] = 1.;
	enm2 = na - 1;
	if (enm2 == 0) {
	    goto L800;
	}
/*     .......... for i=en-2 step -1 until 1 do -- .......... */
	i_2 = enm2;
	for (ii = 1; ii <= i_2; ++ii) {
	    i = na - ii;
	    w = h[i + i * h_dim1] - p;
	    ra = 0.;
	    sa = 0.;

	    i_3 = en;
	    for (j = m; j <= i_3; ++j) {
		ra += h[i + j * h_dim1] * h[j + na * h_dim1];
		sa += h[i + j * h_dim1] * h[j + en * h_dim1];
/* L760: */
	    }

	    if (wi[i] >= 0.) {
		goto L770;
	    }
	    zz = w;
	    r = ra;
	    s = sa;
	    goto L795;
L770:
	    m = i;
	    if (wi[i] != 0.) {
		goto L780;
	    }
	    d_1 = -ra;
	    d_2 = -sa;
	    cdiv_(&d_1, &d_2, &w, &q, &h[i + na * h_dim1], &h[i + en * h_dim1]
		    );
	    goto L790;
/*     .......... solve complex equations .......... */
L780:
	    x = h[i + (i + 1) * h_dim1];
	    y = h[i + 1 + i * h_dim1];
	    vr = (wr[i] - p) * (wr[i] - p) + wi[i] * wi[i] - q * q;
	    vi = (wr[i] - p) * 2. * q;
	    if (vr != 0. || vi != 0.) {
		goto L784;
	    }
	    tst1 = norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(zz));
	    vr = tst1;
L783:
	    vr *= .01;
	    tst2 = tst1 + vr;
	    if (tst2 > tst1) {
		goto L783;
	    }
L784:
	    d_1 = x * r - zz * ra + q * sa;
	    d_2 = x * s - zz * sa - q * ra;
	    cdiv_(&d_1, &d_2, &vr, &vi, &h[i + na * h_dim1], &h[i + en * 
		    h_dim1]);
	    if (abs(x) <= abs(zz) + abs(q)) {
		goto L785;
	    }
	    h[i + 1 + na * h_dim1] = (-ra - w * h[i + na * h_dim1] + q * h[i 
		    + en * h_dim1]) / x;
	    h[i + 1 + en * h_dim1] = (-sa - w * h[i + en * h_dim1] - q * h[i 
		    + na * h_dim1]) / x;
	    goto L790;
L785:
	    d_1 = -r - y * h[i + na * h_dim1];
	    d_2 = -s - y * h[i + en * h_dim1];
	    cdiv_(&d_1, &d_2, &zz, &q, &h[i + 1 + na * h_dim1], &h[i + 1 + en 
		    * h_dim1]);

/*     .......... overflow control .......... */
L790:
/* Computing MAX */
	    d_3 = (d_1 = h[i + na * h_dim1], abs(d_1)), d_4 = (d_2 = h[i + en 
		    * h_dim1], abs(d_2));
	    t = max(d_4,d_3);
	    if (t == 0.) {
		goto L795;
	    }
	    tst1 = t;
	    tst2 = tst1 + 1. / tst1;
	    if (tst2 > tst1) {
		goto L795;
	    }
	    i_3 = en;
	    for (j = i; j <= i_3; ++j) {
		h[j + na * h_dim1] /= t;
		h[j + en * h_dim1] /= t;
/* L792: */
	    }

L795:
	;}
/*     .......... end complex vector .......... */
L800:
    ;}
/*     .......... end back substitution. */
/*                vectors of isolated roots .......... */
    i_1 = *n;
    for (i = 1; i <= i_1; ++i) {
	if (i >= *low && i <= *igh) {
	    goto L840;
	}

	i_2 = *n;
	for (j = i; j <= i_2; ++j) {
/* L820: */
	    z[i + j * z_dim1] = h[i + j * h_dim1];
	}

L840:
    ;}
/*     .......... multiply by transformation matrix to give */
/*                vectors of original full matrix. */
/*                for j=n step -1 until low do -- .......... */
    i_1 = *n;
    for (jj = *low; jj <= i_1; ++jj) {
	j = *n + *low - jj;
	m = min(j,*igh);

	i_2 = *igh;
	for (i = *low; i <= i_2; ++i) {
	    zz = 0.;

	    i_3 = m;
	    for (k = *low; k <= i_3; ++k) {
/* L860: */
		zz += z[i + k * z_dim1] * h[k + j * h_dim1];
	    }

	    z[i + j * z_dim1] = zz;
/* L880: */
	}
    }

    goto L1001;
/*     .......... set error -- all eigenvalues have not */
/*                converged after 30*n iterations .......... */
L1000:
    *ierr = en;
L1001:
    return 0;
} /* hqr2_ */

/* Subroutine */ int cdiv_(ar, ai, br, bi, cr, ci)
doublereal *ar, *ai, *br, *bi, *cr, *ci;
{
    /* System generated locals */
    doublereal d_1, d_2;

    /* Local variables */
    static doublereal s, ais, bis, ars, brs;


/*     complex division, (cr,ci) = (ar,ai)/(br,bi) */

    s = abs(*br) + abs(*bi);
    ars = *ar / s;
    ais = *ai / s;
    brs = *br / s;
    bis = *bi / s;
/* Computing 2nd power */
    d_1 = brs;
/* Computing 2nd power */
    d_2 = bis;
    s = d_1 * d_1 + d_2 * d_2;
    *cr = (ars * brs + ais * bis) / s;
    *ci = (ais * brs - ars * bis) / s;
    return 0;
} /* cdiv_ */

/* Subroutine */ int eltran_(nm, n, low, igh, a, int_, z)
integer *nm, *n, *low, *igh;
doublereal *a;
integer *int_;
doublereal *z;
{
    /* System generated locals */
    integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2;

    /* Local variables */
    static integer i, j, kl, mm, mp, mp1;

    /* Parameter adjustments */
    a_dim1 = *nm;
    a_offset = a_dim1 + 1;
    a -= a_offset;
    --int_;
    z_dim1 = *nm;
    z_offset = z_dim1 + 1;
    z -= z_offset;

    /* Function Body */


/*     this subroutine is a translation of the algol procedure elmtrans, 
*/
/*     num. math. 16, 181-204(1970) by peters and wilkinson. */
/*     handbook for auto. comp., vol.ii-linear algebra, 372-395(1971). */

/*     this subroutine accumulates the stabilized elementary */
/*     similarity transformations used in the reduction of a */
/*     real general matrix to upper hessenberg form by  elmhes. */

/*     on input */

/*        nm must be set to the row dimension of two-dimensional */
/*          array parameters as declared in the calling program */
/*          dimension statement. */

/*        n is the order of the matrix. */

/*        low and igh are integers determined by the balancing */
/*          subroutine  balanc.  if  balanc  has not been used, */
/*          set low=1, igh=n. */

/*        a contains the multipliers which were used in the */
/*          reduction by  elmhes  in its lower triangle */
/*          below the subdiagonal. */

/*        int contains information on the rows and columns */
/*          interchanged in the reduction by  elmhes. */
/*          only elements low through igh are used. */

/*     on output */

/*        z contains the transformation matrix produced in the */
/*          reduction by  elmhes. */

/*     questions and comments should be directed to burton s. garbow, */
/*     mathematics and computer science div, argonne national laboratory 
*/

/*     this version dated august 1983. */

/*     ------------------------------------------------------------------ 
*/

/*     .......... initialize z to identity matrix .......... */
    i_1 = *n;
    for (j = 1; j <= i_1; ++j) {

	i_2 = *n;
	for (i = 1; i <= i_2; ++i) {
/* L60: */
	    z[i + j * z_dim1] = 0.;
	}

	z[j + j * z_dim1] = 1.;
/* L80: */
    }

    kl = *igh - *low - 1;
    if (kl < 1) {
	goto L200;
    }
/*     .......... for mp=igh-1 step -1 until low+1 do -- .......... */
    i_1 = kl;
    for (mm = 1; mm <= i_1; ++mm) {
	mp = *igh - mm;
	mp1 = mp + 1;

	i_2 = *igh;
	for (i = mp1; i <= i_2; ++i) {
/* L100: */
	    z[i + mp * z_dim1] = a[i + (mp - 1) * a_dim1];
	}

	i = int_[mp];
	if (i == mp) {
	    goto L140;
	}

	i_2 = *igh;
	for (j = mp; j <= i_2; ++j) {
	    z[mp + j * z_dim1] = z[i + j * z_dim1];
	    z[i + j * z_dim1] = 0.;
/* L130: */
	}

	z[i + mp * z_dim1] = 1.;
L140:
    ;}

L200:
    return 0;
} /* eltran_ */

/* Subroutine */ int hqr_(nm, n, low, igh, h, wr, wi, ierr)
integer *nm, *n, *low, *igh;
doublereal *h, *wr, *wi;
integer *ierr;
{
    /* System generated locals */
    integer h_dim1, h_offset, i_1, i_2, i_3;
    doublereal d_1, d_2;

    /* Builtin functions */
    double sqrt(), d_sign();

    /* Local variables */
    static doublereal norm;
    static integer i, j, k, l, m;
    static doublereal p, q, r, s, t, w, x, y;
    static integer na, en, ll, mm;
    static doublereal zz;
    static logical notlas;
    static integer mp2, itn, its, enm2;
    static doublereal tst1, tst2;

    /* Parameter adjustments */
    h_dim1 = *nm;
    h_offset = h_dim1 + 1;
    h -= h_offset;
    --wr;
    --wi;

    /* Function Body */
/*  RESTORED CORRECT INDICES OF LOOPS (200,210,230,240). (9/29/89 BSG) */


/*     this subroutine is a translation of the algol procedure hqr, */
/*     num. math. 14, 219-231(1970) by martin, peters, and wilkinson. */
/*     handbook for auto. comp., vol.ii-linear algebra, 359-371(1971). */

/*     this subroutine finds the eigenvalues of a real */
/*     upper hessenberg matrix by the qr method. */

/*     on input */

/*        nm must be set to the row dimension of two-dimensional */
/*          array parameters as declared in the calling program */
/*          dimension statement. */

/*        n is the order of the matrix. */

/*        low and igh are integers determined by the balancing */
/*          subroutine  balanc.  if  balanc  has not been used, */
/*          set low=1, igh=n. */

/*        h contains the upper hessenberg matrix.  information about */
/*          the transformations used in the reduction to hessenberg */
/*          form by  elmhes  or  orthes, if performed, is stored */
/*          in the remaining triangle under the hessenberg matrix. */

/*     on output */

/*        h has been destroyed.  therefore, it must be saved */
/*          before calling  hqr  if subsequent calculation and */
/*          back transformation of eigenvectors is to be performed. */

/*        wr and wi contain the real and imaginary parts, */
/*          respectively, of the eigenvalues.  the eigenvalues */
/*          are unordered except that complex conjugate pairs */
/*          of values appear consecutively with the eigenvalue */
/*          having the positive imaginary part first.  if an */
/*          error exit is made, the eigenvalues should be correct */
/*          for indices ierr+1,...,n. */

/*        ierr is set to */
/*          zero       for normal return, */
/*          j          if the limit of 30*n iterations is exhausted */
/*                     while the j-th eigenvalue is being sought. */

/*     questions and comments should be directed to burton s. garbow, */
/*     mathematics and computer science div, argonne national laboratory 
*/

/*     this version dated september 1989. */

/*     ------------------------------------------------------------------ 
*/

    *ierr = 0;
    norm = 0.;
    k = 1;
/*     .......... store roots isolated by balanc */
/*                and compute matrix norm .......... */
    i_1 = *n;
    for (i = 1; i <= i_1; ++i) {

	i_2 = *n;
	for (j = k; j <= i_2; ++j) {
/* L40: */
	    norm += (d_1 = h[i + j * h_dim1], abs(d_1));
	}

	k = i;
	if (i >= *low && i <= *igh) {
	    goto L50;
	}
	wr[i] = h[i + i * h_dim1];
	wi[i] = 0.;
L50:
    ;}

    en = *igh;
    t = 0.;
    itn = *n * 30;
/*     .......... search for next eigenvalues .......... */
L60:
    if (en < *low) {
	goto L1001;
    }
    its = 0;
    na = en - 1;
    enm2 = na - 1;
/*     .......... look for single small sub-diagonal element */
/*                for l=en step -1 until low do -- .......... */
L70:
    i_1 = en;
    for (ll = *low; ll <= i_1; ++ll) {
	l = en + *low - ll;
	if (l == *low) {
	    goto L100;
	}
	s = (d_1 = h[l - 1 + (l - 1) * h_dim1], abs(d_1)) + (d_2 = h[l + l * 
		h_dim1], abs(d_2));
	if (s == 0.) {
	    s = norm;
	}
	tst1 = s;
	tst2 = tst1 + (d_1 = h[l + (l - 1) * h_dim1], abs(d_1));
	if (tst2 == tst1) {
	    goto L100;
	}
/* L80: */
    }
/*     .......... form shift .......... */
L100:
    x = h[en + en * h_dim1];
    if (l == en) {
	goto L270;
    }
    y = h[na + na * h_dim1];
    w = h[en + na * h_dim1] * h[na + en * h_dim1];
    if (l == na) {
	goto L280;
    }
    if (itn == 0) {
	goto L1000;
    }
    if (its != 10 && its != 20) {
	goto L130;
    }
/*     .......... form exceptional shift .......... */
    t += x;

    i_1 = en;
    for (i = *low; i <= i_1; ++i) {
/* L120: */
	h[i + i * h_dim1] -= x;
    }

    s = (d_1 = h[en + na * h_dim1], abs(d_1)) + (d_2 = h[na + enm2 * h_dim1], 
	    abs(d_2));
    x = s * .75;
    y = x;
    w = s * -.4375 * s;
L130:
    ++its;
    --itn;
/*     .......... look for two consecutive small */
/*                sub-diagonal elements. */
/*                for m=en-2 step -1 until l do -- .......... */
    i_1 = enm2;
    for (mm = l; mm <= i_1; ++mm) {
	m = enm2 + l - mm;
	zz = h[m + m * h_dim1];
	r = x - zz;
	s = y - zz;
	p = (r * s - w) / h[m + 1 + m * h_dim1] + h[m + (m + 1) * h_dim1];
	q = h[m + 1 + (m + 1) * h_dim1] - zz - r - s;
	r = h[m + 2 + (m + 1) * h_dim1];
	s = abs(p) + abs(q) + abs(r);
	p /= s;
	q /= s;
	r /= s;
	if (m == l) {
	    goto L150;
	}
	tst1 = abs(p) * ((d_1 = h[m - 1 + (m - 1) * h_dim1], abs(d_1)) + abs(
		zz) + (d_2 = h[m + 1 + (m + 1) * h_dim1], abs(d_2)));
	tst2 = tst1 + (d_1 = h[m + (m - 1) * h_dim1], abs(d_1)) * (abs(q) + 
		abs(r));
	if (tst2 == tst1) {
	    goto L150;
	}
/* L140: */
    }

L150:
    mp2 = m + 2;

    i_1 = en;
    for (i = mp2; i <= i_1; ++i) {
	h[i + (i - 2) * h_dim1] = 0.;
	if (i == mp2) {
	    goto L160;
	}
	h[i + (i - 3) * h_dim1] = 0.;
L160:
    ;}
/*     .......... double qr step involving rows l to en and */
/*                columns m to en .......... */
    i_1 = na;
    for (k = m; k <= i_1; ++k) {
	notlas = k != na;
	if (k == m) {
	    goto L170;
	}
	p = h[k + (k - 1) * h_dim1];
	q = h[k + 1 + (k - 1) * h_dim1];
	r = 0.;
	if (notlas) {
	    r = h[k + 2 + (k - 1) * h_dim1];
	}
	x = abs(p) + abs(q) + abs(r);
	if (x == 0.) {
	    goto L260;
	}
	p /= x;
	q /= x;
	r /= x;
L170:
	d_1 = sqrt(p * p + q * q + r * r);
	s = d_sign(&d_1, &p);
	if (k == m) {
	    goto L180;
	}
	h[k + (k - 1) * h_dim1] = -s * x;
	goto L190;
L180:
	if (l != m) {
	    h[k + (k - 1) * h_dim1] = -h[k + (k - 1) * h_dim1];
	}
L190:
	p += s;
	x = p / s;
	y = q / s;
	zz = r / s;
	q /= p;
	r /= p;
	if (notlas) {
	    goto L225;
	}
/*     .......... row modification .......... */
	i_2 = en;
	for (j = k; j <= i_2; ++j) {
	    p = h[k + j * h_dim1] + q * h[k + 1 + j * h_dim1];
	    h[k + j * h_dim1] -= p * x;
	    h[k + 1 + j * h_dim1] -= p * y;
/* L200: */
	}

/* Computing MAX */
	i_2 = en, i_3 = k + 3;
	j = min(i_3,i_2);
/*     .......... column modification .......... */
	i_2 = j;
	for (i = l; i <= i_2; ++i) {
	    p = x * h[i + k * h_dim1] + y * h[i + (k + 1) * h_dim1];
	    h[i + k * h_dim1] -= p;
	    h[i + (k + 1) * h_dim1] -= p * q;
/* L210: */
	}
	goto L255;
L225:
/*     .......... row modification .......... */
	i_2 = en;
	for (j = k; j <= i_2; ++j) {
	    p = h[k + j * h_dim1] + q * h[k + 1 + j * h_dim1] + r * h[k + 2 + 
		    j * h_dim1];
	    h[k + j * h_dim1] -= p * x;
	    h[k + 1 + j * h_dim1] -= p * y;
	    h[k + 2 + j * h_dim1] -= p * zz;
/* L230: */
	}

/* Computing MAX */
	i_2 = en, i_3 = k + 3;
	j = min(i_3,i_2);
/*     .......... column modification .......... */
	i_2 = j;
	for (i = l; i <= i_2; ++i) {
	    p = x * h[i + k * h_dim1] + y * h[i + (k + 1) * h_dim1] + zz * h[
		    i + (k + 2) * h_dim1];
	    h[i + k * h_dim1] -= p;
	    h[i + (k + 1) * h_dim1] -= p * q;
	    h[i + (k + 2) * h_dim1] -= p * r;
/* L240: */
	}
L255:

L260:
    ;}

    goto L70;
/*     .......... one root found .......... */
L270:
    wr[en] = x + t;
    wi[en] = 0.;
    en = na;
    goto L60;
/*     .......... two roots found .......... */
L280:
    p = (y - x) / 2.;
    q = p * p + w;
    zz = sqrt((abs(q)));
    x += t;
    if (q < 0.) {
	goto L320;
    }
/*     .......... real pair .......... */
    zz = p + d_sign(&zz, &p);
    wr[na] = x + zz;
    wr[en] = wr[na];
    if (zz != 0.) {
	wr[en] = x - w / zz;
    }
    wi[na] = 0.;
    wi[en] = 0.;
    goto L330;
/*     .......... complex pair .......... */
L320:
    wr[na] = x + p;
    wr[en] = x + p;
    wi[na] = zz;
    wi[en] = -zz;
L330:
    en = enm2;
    goto L60;
/*     .......... set error -- all eigenvalues have not */
/*                converged after 30*n iterations .......... */
L1000:
    *ierr = en;
L1001:
    return 0;
} /* hqr_ */

/* Subroutine */ int elmhes_(nm, n, low, igh, a, int_)
integer *nm, *n, *low, *igh;
doublereal *a;
integer *int_;
{
    /* System generated locals */
    integer a_dim1, a_offset, i_1, i_2, i_3;
    doublereal d_1;

    /* Local variables */
    static integer i, j, m;
    static doublereal x, y;
    static integer la, mm1, kp1, mp1;

    /* Parameter adjustments */
    a_dim1 = *nm;
    a_offset = a_dim1 + 1;
    a -= a_offset;
    --int_;

    /* Function Body */


/*     this subroutine is a translation of the algol procedure elmhes, */
/*     num. math. 12, 349-368(1968) by martin and wilkinson. */
/*     handbook for auto. comp., vol.ii-linear algebra, 339-358(1971). */

/*     given a real general matrix, this subroutine */
/*     reduces a submatrix situated in rows and columns */
/*     low through igh to upper hessenberg form by */
/*     stabilized elementary similarity transformations. */

/*     on input */

/*        nm must be set to the row dimension of two-dimensional */
/*          array parameters as declared in the calling program */
/*          dimension statement. */

/*        n is the order of the matrix. */

/*        low and igh are integers determined by the balancing */
/*          subroutine  balanc.  if  balanc  has not been used, */
/*          set low=1, igh=n. */

/*        a contains the input matrix. */

/*     on output */

/*        a contains the hessenberg matrix.  the multipliers */
/*          which were used in the reduction are stored in the */
/*          remaining triangle under the hessenberg matrix. */

/*        int contains information on the rows and columns */
/*          interchanged in the reduction. */
/*          only elements low through igh are used. */

/*     questions and comments should be directed to burton s. garbow, */
/*     mathematics and computer science div, argonne national laboratory 
*/

/*     this version dated august 1983. */

/*     ------------------------------------------------------------------ 
*/

    la = *igh - 1;
    kp1 = *low + 1;
    if (la < kp1) {
	goto L200;
    }

    i_1 = la;
    for (m = kp1; m <= i_1; ++m) {
	mm1 = m - 1;
	x = 0.;
	i = m;

	i_2 = *igh;
	for (j = m; j <= i_2; ++j) {
	    if ((d_1 = a[j + mm1 * a_dim1], abs(d_1)) <= abs(x)) {
		goto L100;
	    }
	    x = a[j + mm1 * a_dim1];
	    i = j;
L100:
	;}

	int_[m] = i;
	if (i == m) {
	    goto L130;
	}
/*     .......... interchange rows and columns of a .......... */
	i_2 = *n;
	for (j = mm1; j <= i_2; ++j) {
	    y = a[i + j * a_dim1];
	    a[i + j * a_dim1] = a[m + j * a_dim1];
	    a[m + j * a_dim1] = y;
/* L110: */
	}

	i_2 = *igh;
	for (j = 1; j <= i_2; ++j) {
	    y = a[j + i * a_dim1];
	    a[j + i * a_dim1] = a[j + m * a_dim1];
	    a[j + m * a_dim1] = y;
/* L120: */
	}
/*     .......... end interchange .......... */
L130:
	if (x == 0.) {
	    goto L180;
	}
	mp1 = m + 1;

	i_2 = *igh;
	for (i = mp1; i <= i_2; ++i) {
	    y = a[i + mm1 * a_dim1];
	    if (y == 0.) {
		goto L160;
	    }
	    y /= x;
	    a[i + mm1 * a_dim1] = y;

	    i_3 = *n;
	    for (j = m; j <= i_3; ++j) {
/* L140: */
		a[i + j * a_dim1] -= y * a[m + j * a_dim1];
	    }

	    i_3 = *igh;
	    for (j = 1; j <= i_3; ++j) {
/* L150: */
		a[j + m * a_dim1] += y * a[j + i * a_dim1];
	    }

L160:
	;}

L180:
    ;}

L200:
    return 0;
} /* elmhes_ */

/* Subroutine */ int balanc_(nm, n, a, low, igh, scale)
integer *nm, *n;
doublereal *a;
integer *low, *igh;
doublereal *scale;
{
    /* System generated locals */
    integer a_dim1, a_offset, i_1, i_2;
    doublereal d_1;

    /* Local variables */
    static integer iexc;
    static doublereal c, f, g;
    static integer i, j, k, l, m;
    static doublereal r, s, radix, b2;
    static integer jj;
    static logical noconv;

    /* Parameter adjustments */
    a_dim1 = *nm;
    a_offset = a_dim1 + 1;
    a -= a_offset;
    --scale;

    /* Function Body */


/*     this subroutine is a translation of the algol procedure balance, */

/*     num. math. 13, 293-304(1969) by parlett and reinsch. */
/*     handbook for auto. comp., vol.ii-linear algebra, 315-326(1971). */

/*     this subroutine balances a real matrix and isolates */
/*     eigenvalues whenever possible. */

/*     on input */

/*        nm must be set to the row dimension of two-dimensional */
/*          array parameters as declared in the calling program */
/*          dimension statement. */

/*        n is the order of the matrix. */

/*        a contains the input matrix to be balanced. */

/*     on output */

/*        a contains the balanced matrix. */

/*        low and igh are two integers such that a(i,j) */
/*          is equal to zero if */
/*           (1) i is greater than j and */
/*           (2) j=1,...,low-1 or i=igh+1,...,n. */

/*        scale contains information determining the */
/*           permutations and scaling factors used. */

/*     suppose that the principal submatrix in rows low through igh */
/*     has been balanced, that p(j) denotes the index interchanged */
/*     with j during the permutation step, and that the elements */
/*     of the diagonal matrix used are denoted by d(i,j).  then */
/*        scale(j) = p(j),    for j = 1,...,low-1 */
/*                 = d(j,j),      j = low,...,igh */
/*                 = p(j)         j = igh+1,...,n. */
/*     the order in which the interchanges are made is n to igh+1, */
/*     then 1 to low-1. */

/*     note that 1 is returned for igh if igh is zero formally. */

/*     the algol procedure exc contained in balance appears in */
/*     balanc  in line.  (note that the algol roles of identifiers */
/*     k,l have been reversed.) */

/*     questions and comments should be directed to burton s. garbow, */
/*     mathematics and computer science div, argonne national laboratory 
*/

/*     this version dated august 1983. */

/*     ------------------------------------------------------------------ 
*/

    radix = 16.;

    b2 = radix * radix;
    k = 1;
    l = *n;
    goto L100;
/*     .......... in-line procedure for row and */
/*                column exchange .......... */
L20:
    scale[m] = (doublereal) j;
    if (j == m) {
	goto L50;
    }

    i_1 = l;
    for (i = 1; i <= i_1; ++i) {
	f = a[i + j * a_dim1];
	a[i + j * a_dim1] = a[i + m * a_dim1];
	a[i + m * a_dim1] = f;
/* L30: */
    }

    i_1 = *n;
    for (i = k; i <= i_1; ++i) {
	f = a[j + i * a_dim1];
	a[j + i * a_dim1] = a[m + i * a_dim1];
	a[m + i * a_dim1] = f;
/* L40: */
    }

L50:
    switch (iexc) {
	case 1:  goto L80;
	case 2:  goto L130;
    }
/*     .......... search for rows isolating an eigenvalue */
/*                and push them down .......... */
L80:
    if (l == 1) {
	goto L280;
    }
    --l;
/*     .......... for j=l step -1 until 1 do -- .......... */
L100:
    i_1 = l;
    for (jj = 1; jj <= i_1; ++jj) {
	j = l + 1 - jj;

	i_2 = l;
	for (i = 1; i <= i_2; ++i) {
	    if (i == j) {
		goto L110;
	    }
	    if (a[j + i * a_dim1] != 0.) {
		goto L120;
	    }
L110:
	;}

	m = l;
	iexc = 1;
	goto L20;
L120:
    ;}

    goto L140;
/*     .......... search for columns isolating an eigenvalue */
/*                and push them left .......... */
L130:
    ++k;

L140:
    i_1 = l;
    for (j = k; j <= i_1; ++j) {

	i_2 = l;
	for (i = k; i <= i_2; ++i) {
	    if (i == j) {
		goto L150;
	    }
	    if (a[i + j * a_dim1] != 0.) {
		goto L170;
	    }
L150:
	;}

	m = k;
	iexc = 2;
	goto L20;
L170:
    ;}
/*     .......... now balance the submatrix in rows k to l .......... */
    i_1 = l;
    for (i = k; i <= i_1; ++i) {
/* L180: */
	scale[i] = 1.;
    }
/*     .......... iterative loop for norm reduction .......... */
L190:
    noconv = FALSE_;

    i_1 = l;
    for (i = k; i <= i_1; ++i) {
	c = 0.;
	r = 0.;

	i_2 = l;
	for (j = k; j <= i_2; ++j) {
	    if (j == i) {
		goto L200;
	    }
	    c += (d_1 = a[j + i * a_dim1], abs(d_1));
	    r += (d_1 = a[i + j * a_dim1], abs(d_1));
L200:
	;}
/*     .......... guard against zero c or r due to underflow 
.......... */
	if (c == 0. || r == 0.) {
	    goto L270;
	}
	g = r / radix;
	f = 1.;
	s = c + r;
L210:
	if (c >= g) {
	    goto L220;
	}
	f *= radix;
	c *= b2;
	goto L210;
L220:
	g = r * radix;
L230:
	if (c < g) {
	    goto L240;
	}
	f /= radix;
	c /= b2;
	goto L230;
/*     .......... now balance .......... */
L240:
	if ((c + r) / f >= s * .95) {
	    goto L270;
	}
	g = 1. / f;
	scale[i] *= f;
	noconv = TRUE_;

	i_2 = *n;
	for (j = k; j <= i_2; ++j) {
/* L250: */
	    a[i + j * a_dim1] *= g;
	}

	i_2 = l;
	for (j = 1; j <= i_2; ++j) {
/* L260: */
	    a[j + i * a_dim1] *= f;
	}

L270:
    ;}

    if (noconv) {
	goto L190;
    }

L280:
    *low = k;
    *igh = l;
    return 0;
} /* balanc_ */

End of rg.c
echo matrix.cc 1>&2
cat >matrix.cc <<'End of matrix.cc'
/*************************************************************************
**************************************************************************
**************************************************************************
******                                                                ****
******  Matrix. A C++ matrix class                                    ****
******                                                                ****
******  Copyright (c) 1993  Universities Research Association, Inc.   ****
******                All Rights Reserved                             ****
******                                                                ****
******  Author:    James Holt                                         ****
******                                                                ****
******   Taken and modified from the book C++ Inside and out          ****
******   by Bruce Eckel                                               ****
******                                                                ****
******                                                                ****
******                                                                ****
******   Eigen Vector/Value routine taken from eispack                ****
******                                                                ****
******             Fermilab                                           ****
******             P.O.Box 500                                        ****
******             Mail Stop 345                                      ****
******             Batavia, IL   60510                                ****
******                                                                ****
******             Phone: (708) 840 8385                              ****
******             Email: holt@calvin.fnal.gov                        ****
******                                                                ****
******  Release    Version 1.0 : October  4, 1993                     ****
**************************************************************************
**************************************************************************
******                                                                ****
******  This material resulted from work developed under a            ****
******  Government Contract and is subject to the following           ****
******  license: The Government retains a paid-up, nonexclusive,      ****
******  irrevocable worldwide license to reproduce, prepare           ****
******  derivative works, perform publicly and display publicly       ****
******  by or for the Government, including the right to distribute   ****
******  to other Government contractors.  Neither the United          ****
******  States nor the United States Department of Energy, nor        ****
******  any of their employees, makes any warranty, express or        ****
******  implied, or assumes any legal liability or responsibility     ****
******  for the accuracy, completeness or usefulness of any           ****
******  information, apparatus, product, or process disclosed, or     ****
******  represents that its use would not infringe privately owned    ****
******  rights.                                                       ****
******                                                                ****
******  These files are made avaliable for use by specific            ****
******  individuals or a specific group.  They are not to be          ****
******  reproduced or redistributed.                                  ****
******                                                                ****
******  These files are provided "as is" with no warranties of any    ****
******  kind including the warranties of design,                      ****
******  merchantibility and fitness for a particular purpose,         ****
******  or arising from a course of dealing, usage or trade           ****
******  practice.                                                     ****
******                                                                ****
******  These files are provided with no support or obligation on     ****
******  the part of U.R.A. or Fermilab to assist in their use,        ****
******  correction, modification or enhancement.                      ****
******                                                                ****
******  Neither U.R.A. nor Fermilab shall be liable with              ****
******  respect to the infringement of copyrights, trade              ****
******  secrets or any patents by these files or any part             ****
******  thereof.                                                      ****
******                                                                ****
******  In no event will U.R.A., Fermilab, or the author(s) be        ****
******  liable for any lost revenue or profits or other               ****
******  special, indirect and consequential damages, even if          ****
******  they have been advised of the possibility of such             ****
******  damages.                                                      ****
******                                                                ****
**************************************************************************
*************************************************************************/

#include "matrix.h"

extern "C" {strstr(char*,char*);} // this may need to be commented out
				  // on some systems.
extern "C" { void rg_(int*, int*, double*, double*,double*, int*, double*,
		 int*, double*, int*); }

extern "C" { void cg_(int*, int*, double*, double*,double*,double*, int*, 
		      double*,double*,double*, double*, double*,int*); }

const double tiny = double(1e-20);
const double limit = double(1e-14);

#define M_SMALL  1.0e-14

#define WID setw(8)
#define PREC setprecision(4)

void doublematrix::error(char * msg1, char * msg2) const {
  cerr << "doublematrix error: " << msg1 << " " << msg2  << endl;
  exit(1);
}

doublematrix::doublematrix(int mrows, int columns, double initval) {
  // create the structure:
  p = new doublematrep;
  p->r = mrows;
  p->c = columns;
  // allocate memory for the actual matrix:
  p->m = new double *[mrows];
  for (int x = 0; x < mrows; x++)
    p->m[x] = new double[columns];
  p->n = 1;  // so far, there's one reference to this data
  for (int i=0; i< mrows; i++) {
    for (int j = 0; j < columns; j++)
      mval(i,j) = initval;
  }
}

doublematrix::doublematrix(int mrows, int columns, double * initvalues){
  // create the structure:
  p    = new doublematrep;
  p->r = mrows;
  p->c = columns;
  // allocate memory for the actual matrix:
  p->m = new double *[mrows];
  for (int x = 0; x < mrows; x++)
    p->m[x] = new double[columns];
  p->n = 1;  // so far, there's one reference to this data:
  int c = 0;
  for (int i=0; i< mrows; i++) {
    for (int j = 0; j < columns; j++)
      mval(i,j) =  initvalues[c++];
  }
}
// create an identity matrix:
doublematrix::doublematrix(char * flag, int dimension) {
  int i,j;
  float tmp_float = dimension/2.0;
  
  if(((dimension/2) != tmp_float) && (flag[0] == 'J'))
    error("dimension must be even","for J matrix");
  if (flag[0]  == 'I') {
    p = new doublematrep;
    p->r = dimension;
    p->c = dimension;
    p->m = new double *[dimension];
    for (int x = 0; x < dimension; x++)
      p->m[x] = new double[dimension];
    p->n = 1;
    for (i = 0; i< dimension; i++) {
      for (j = 0; j < dimension; j++)
	mval(i,j) = (i == j ? 1.0 : 0.0);
    }
  } else if (flag[0] == 'J') {
    p = new doublematrep;
    p->r = dimension;
    p->c = dimension;
    p->m = new double *[dimension];
    for (int x = 0; x < dimension; x++)
      p->m[x] = new double[dimension];
    p->n = 1;
    for ( i = 0; i< dimension; i++) {
      for ( j = 0; j < dimension; j++)
	mval(i,j) = 0.0;
    }
    for (i = dimension/2; i< dimension; i++) {
      mval(i-dimension/2,i) = 1;
      mval(i,i-dimension/2) = -1;
    }
  } else
    error("Unknown flag: ",flag);
}

// error message when trying to read a "standard"
// matrix file:

static char nonstandard[] =
" is a 'non-standard' file. A 'standard' matrix file must \n"
"start with the dimensions of the matrix, i.e.:\n"
"\t rows=12 columns=14\n or abbreviated:\n\t r=12 c=14 \n"
"Notice rows appear before columns, and chars are lower case\n"
"comments follow '#' signs to end of line, data follows :::\n";

// read from "standard" matrix file:
doublematrix::doublematrix(char * initfile) {
  ifstream from(initfile);
  if(!from)
    error("cannot open doublematrix initializer file", initfile);
  const bsize = 120;
  char buf[bsize], *cp, *cp2;
  int rfound = 0, cfound = 0, colonsfound = 0;
  p = new doublematrep;
  // Parse file initialization header 
  while(from.getline(buf, bsize)) { // for each header line
    // Remove comments with ANSI C library function "strprk":
    if( ( cp = strpbrk(buf,"#"))  != 0 ) // look for comments
      *cp = '\0';  // terminate string at comment
    if( ( cp = strpbrk(buf,"r") )  != 0 )
      if ( ( cp2 = strpbrk(cp,  "="))  !=0)
	if ( ( cp = strpbrk(cp2,  "0123456789"))  != 0 ) {
	 p->r = atoi(cp);
	 rfound++;  // flag to say rows were found
      }
    if( ( cp = strpbrk(buf,"c") )  != 0)
      if ( ( cp2 = strpbrk(cp,  "="))  != 0)
	if ( ( cp = strpbrk(cp2,  "0123456789"))  != 0 ) {
	 p->c = atoi(cp);
	 cfound++;  // flag to say cols were found
	}
    if ( strstr(buf,":::")  != 0 )  {
      colonsfound++;
      break; // ... out of "while" loop
    }
  }

  if ( !rfound || !cfound || !colonsfound )  {
    cerr << initfile << nonstandard;
    exit(1);
  }

  p->m = new double *[p->r];
  for (int x = 0; x < p->r; x++)
  p->m[x] = new double[p->c];
  p->n = 1;  // so far, there's one reference to this data
  for (int row = 0; row < p->r; row++)  {
    for(int col = 0; col < p->c; col++)  {
      const nbsz = 20; char nb[nbsz];
      from >> nb; // scan for space-delimited string
      mval(row,col) = atof(nb); // convert int to a double
      if(from.bad())
	error ( "problem with doublematrix initializer file", initfile);
    }
  }
}

doublematrix::doublematrix(doublematrix& x) {
  p = new doublematrep;
  p->r = x.p->r;
  p->c = x.p->c;
  // allocate memory for the actual matrix:
  p->m = new double *[x.p->r];
  for (int i = 0; i < x.p->r; i++)
    p->m[i] = new double[x.p->c];
  p->n = 1;  // so far, there's one reference to this data
  deepcopy(x,*this);
}


doublematrix::~doublematrix()  {
  if (--p->n == 0) { // if reference count goes to 0
    for (int x = 0; x < rows(); x++)
      delete []p->m[x];
    delete []p->m; // delete data
    delete p;
  }
}

doublematrix doublematrix::operator=(doublematrix& x) {

  if(this == &x) return *this;

  p = new doublematrep;
  p->r = x.p->r;
  p->c = x.p->c;
  // allocate memory for the actual matrix:
  p->m = new double *[x.p->r];
  for (int i = 0; i < x.p->r; i++)
    p->m[i] = new double[x.p->c];
  p->n = 1;  // so far, there's one reference to this data
  deepcopy(x,*this);

  return *this;
}

doublematrix doublematrix::operator-() {
  doublematrix unaryminus(rows(),cols());
  for (int i=0; i< rows(); i++) {
    for (int j = 0; j < cols(); j++)
      unaryminus.mval(i,j) = -mval(i,j);
  }
  return unaryminus;
}

void doublematrix::write_standard(char * filename, char * msg) {
  ofstream to(filename);
  if(!to)
    error("cannot open or create doublematrix output file", filename);
  to << "# " << filename
     << ": doublematrix file written in \"standard\" format" <<endl;
//  time_t clock;
//  time(&clock);
//  to <<  "# " << asctime(localtime(&clock));
  to <<  "# " << msg << endl;
  to << "rows= " << rows() << " columns= " << cols() << endl;
  to <<  ":::" << endl;
  for (int row = 0; row < rows(); row++) {
    for(int col = 0; col < cols(); col++)
      to << mval(row,col) << "  ";
    to << endl;
  }
}

doublematrix operator+(const doublematrix& x, const doublematrix& y)  {
  if(( x.rows()  != y.rows())  ||  ( x.cols()  != y.cols()))
    x.error("must have equal dimensions for addition!");
  doublematrix* sum = new doublematrix(x.rows(),x.cols());
  for (int i=0; i< x.rows(); i++) {
    for (int j = 0; j < x.cols(); j++) {
      sum->p->m[i][j] = x.p->m[i][j] + y.p->m[i][j];
      if( fabs(sum->p->m[i][j]) < M_SMALL*fabs(y.p->m[i][j])  ) 
	sum->p->m[i][j] = 0.0;
    }
  }
  return *sum; // see note for operator*()
}

doublematrix operator+(const doublematrix& x, const double arg)  {
  doublematrix* sum = new doublematrix(x.rows(),x.cols());
  for (int i=0; i< x.rows(); i++) {
    for (int j = 0; j < x.cols(); j++) {
      sum->p->m[i][j] = x.p->m[i][j] + arg;
      if( fabs(sum->p->m[i][j]) < M_SMALL*fabs(arg)  ) 
	sum->p->m[i][j] = 0.0;
    }
  }
  return *sum; // see note for operator*()
}

doublematrix operator+( const double arg,const doublematrix& x)  {
  doublematrix* sum = new doublematrix(x.rows(),x.cols());
  for (int i=0; i< x.rows(); i++) {
    for (int j = 0; j < x.cols(); j++) {
      sum->p->m[i][j] = x.p->m[i][j] + arg;
      if( fabs(sum->p->m[i][j]) < M_SMALL*fabs(arg)  ) 
	sum->p->m[i][j] = 0.0;
    }
  }
  return *sum; // see note for operator*()
}

doublematrix operator-(const doublematrix& x,const doublematrix& y)  {
  if(( x.rows()  != y.rows())  || ( x.cols()  != y.cols()))
    x.error("must have equal dimensions for subtraction!");
  doublematrix* sum = new doublematrix(x.rows(),x.cols());
  for (int i=0; i< x.rows(); i++) {
    for (int j = 0; j < x.cols(); j++) {
      sum->p->m[i][j] = x.p->m[i] [j] - y.p->m[i] [j];
      if( fabs(sum->p->m[i][j]) < M_SMALL*fabs(y.p->m[i][j])  ) 
	sum->p->m[i][j] = 0.0;
    }
  }
  return *sum; // see note for operator*()
}

doublematrix operator-(const doublematrix& x, const double arg)  {
  doublematrix* sum = new doublematrix(x.rows(),x.cols());
  for (int i=0; i < x.rows(); i++) {
    for (int j = 0; j < x.cols(); j++) {
      sum->p->m[i][j] = x.p->m[i][j] - arg;
      if( fabs(sum->p->m[i][j]) < M_SMALL*fabs(arg)  ) 
	sum->p->m[i][j] = 0.0;
    }
  }
  return *sum; // see note for operator*()
}

doublematrix operator-( const double arg, const doublematrix& x)  {
  doublematrix* sum = new doublematrix(x.rows(),x.cols());
  for (int i=0; i < x.rows(); i++) {
    for (int j = 0; j < x.cols(); j++) {
      sum->p->m[i][j] = x.p->m[i][j] - arg;
      if( fabs(sum->p->m[i][j]) < M_SMALL*fabs(arg)  ) 
	sum->p->m[i][j] = 0.0;
    }
  }
  return *sum; // see note for operator*()
}


doublematrix operator*(const doublematrix& x, const doublematrix& y)  {
  double tmp;
  if( x.cols()  != y.rows())
    x.error("# rows of second mat must equal "
	  "# cols of first for multiply!");
  doublematrix* result =  new doublematrix(x.rows(),y.cols());
  for(int row = 0; row < x.rows(); row++) {
    for(int col = 0; col < y.cols(); col++){
      double sum = 0;
      for(int i = 0; i < x.cols(); i++) {
	tmp = x.p->m[row][i] * y.p->m[i][col];
	sum += tmp;
	if(fabs(sum) < M_SMALL*fabs(tmp))
	  sum = 0.0;
      }
      result->p->m[row][col] = sum;
    }
  }
  return *result;
}

doublematrix operator*(const doublematrix& x,const double arg) {
  doublematrix* result = new doublematrix(x.rows(),x.cols());
  for (int i=0; i< x.rows(); i++) {
    for (int j = 0; j < x.cols(); j++)
      result->p->m[i][j] = x.p->m[i][j] * arg;
  }
  return *result;
}

doublematrix operator*(const double arg,const doublematrix& x) {
  doublematrix* result = new doublematrix(x.rows(),x.cols());
  for (int i=0; i< x.rows(); i++) {
    for (int j = 0; j < x.cols(); j++)
      result->p->m[i][j] = x.p->m[i][j] * arg;
  }
  return *result;
}

double & doublematrix::val(int row, int col)  {
  //  if (row >= 0 && row < rows() && col >= 0 && col << cols())
  if (row < 0 || row >= rows()  || col < 0 || col >= cols())
    error("index out of range");
  return (mval(row,col));
}

doublematrix doublematrix::transpose() {
  if(rows()  != cols())
    error("doublematrix must be square to transpose!\n");
  doublematrix trans(rows(),cols());
  for (int row = 0; row < rows(); row++)  {
    for(int col = 0; col < cols(); col++)
      trans.mval(col,row) = mval(row,col);
  }
  return trans;
}

double doublematrix::mmin()  {
  double temp;
  if(rows() <= 0 || cols() <= 0)
    error("bad doublematrix size for min()");
  double minimum = mval(0,0);
  for (int row = 0; row < rows(); row++)  {
    for(int col = 0; col < cols(); col++)
      if ((temp = mval(row,col)) < minimum)
	minimum = temp;
  }
  return minimum;
}

double doublematrix::mmax()  {
  double temp;
  if(rows() <= 0 || cols() <= 0)
    error("bad doublematrix size for max()");
  double maximum = mval(0,0);
  for (int row = 0; row < rows(); row++)  {
    for(int col = 0; col < cols(); col++){
      if ((temp = mval(row,col)) > maximum)
	maximum = temp;
    }
  }
  return maximum;
}


double doublematrix::mean()  {
  int row, col;
  double sum = 0;
  for (row = 0; row < rows(); row++)  {
    for(col = 0; col < cols(); col++){
      sum += fabs(mval(row,col));
      if(fabs(sum) < M_SMALL*fabs(mval(row,col)))
	sum = 0.0;
    }
  }
  return sum/(row * col);
}

double doublematrix::variance()  {
  int row, col;
  double s_squared = 0;
  double mn = mean();
  for (row = 0; row < rows(); row++)  {
    for(col = 0; col < cols(); col++) {
      double temp = mval(row,col) - mn;
      temp *= temp;
      s_squared += temp;
    }
  }
  s_squared /= row * col -1; // number of element
  return s_squared;
}

double doublematrix::determinant()  {
  if(rows()  != cols())
    error("doublematrix must be square for determinant()");
  int* indx = new int[cols()]; // create the "index vecto
  doublematrix B(cols());    // see pp 38. in Numerical Recipes
  int d;
  // perform the decomposition once:
  doublematrix decomp = lu_decompose(indx,d);
  double determinant = d;
  for(int i=0; i < cols() ; i++)
    determinant *= decomp.mval(i,i);
  delete []indx;
  return determinant;
}

doublematrix doublematrix::inverse() {
  if(rows()  != cols())
    error("doublematrix must be square for inverse()");
  doublematrix Y("I",rows()); // create an identity doublematrix
  int* indx = new int[cols()];  // create the "index vector"
  doublematrix B(cols());     // see Press & Flannery
  int d;
  // perform the decomposition once:
  doublematrix decomp = lu_decompose(indx,d);
  for(int col = 0; col < cols(); col++){
    B.copy_column(Y,col,0);
    decomp.lu_back_subst(indx,B);
    Y.copy_column(B,0,col);
  }
  delete []indx;
  return Y;
}

complexmatrix& doublematrix::eigenValues() {
  int nm = rows();
  int n  = cols();
  int matz = 0;
  double* wr = new double[cols()];
  double* wi = new double[cols()];
  double z1;
  int* iv1 = new int[cols()];
  double* fv1 = new double[cols()];
  double* b = new double[rows()*cols()];
  int ierr;
  complexmatrix* tmp = new complexmatrix(1,cols());
  int i,j;
  int k = 0;
  for( i = 0; i < rows(); i++) {
    for(j = 0; j< cols(); j++) {
      b[k] = mval(j,i);		// the rg_ routine uses the transpose
      k++;
    }
  }
  rg_(&nm,&n,b,wr,wi,&matz,&z1,iv1,fv1,&ierr);
  if(ierr != 0) {
    cerr << "doublematrix: Error in eigenvalue routine. error = ";
    cerr << ierr<< endl;
    delete []wr;
    delete []wi;
    delete []iv1;
    delete []fv1;
    delete []b;
    return *tmp;
  }

  int* fill = new int[n];
  for(i = 0; i<n; i++)
    fill[i] = 0;
  i = 0;
  k = 0;
  j = 0;
  int tmpInt;
  while( i< nm) {
    if(wi[i] == 0){
      if(fill[k])
	while((k<n) && (fill[k] == 1))
	  k++;
      tmp->mval(0,k) = complex(wr[i],0.0);
      fill[k] = 1;
      k++;
      i++;
    } else {
      tmp->mval(0,k) = complex(wr[i],wi[i]);
      if(k<n/2) {
	tmpInt = k+n/2;
	tmp->mval(0,tmpInt) = complex(wr[i],-wi[i]);
      } else {
	tmp->mval(0,k+1) = complex(wr[i],-wi[i]);
      }
      fill[k] = 1;
      if(k<n/2)
	fill[k+n/2] = 1;
      else
	fill[k+1] = 1;
      k++;
      i += 2;
    }
  }

  delete []fill;
  delete []wr;
  delete []wi;
  delete []iv1;
  delete []fv1;
  delete []b;

  return *tmp;

}

complexmatrix& doublematrix::eigenVectors() {
  int nm = rows();
  int n  = cols();
  int matz = 1;
  double* wr = new double[cols()];
  double* wi = new double[cols()];
  int* iv1 = new int[cols()];
  double* fv1 = new double[cols()];
  double* b = new double[rows()*cols()];
  double* c = new double[rows()*cols()];
  int ierr,i,j,k;
  complexmatrix* tmp = new complexmatrix(rows(),cols());

  k = 0;
  for(i = 0; i < rows(); i++) {
    for(j = 0; j< cols(); j++) {
      b[k] = mval(j,i);		// the rg_ routine uses the transpose
      k++;
    }
  }
  rg_(&nm,&n,b,wr,wi,&matz,c,iv1,fv1,&ierr);

  if(ierr != 0) {
    cerr << "doublematrix: Error in eigenvalue routine. error = ";
    cerr << ierr<< end;
    delete []wr;
    delete []wi;
    delete []iv1;
    delete []fv1;
    delete []b;
    delete []c;
    return *tmp;
  }
  doublematrix z1(nm,n);
  k = 0;
  for(i = 0; i < rows(); i++) {
    for(j = 0; j< cols(); j++) {
      z1(i,j) = c[k];
      k++;
    }
  }

  int* fill = new int[n];
  for(i = 0; i<n; i++)
    fill[i] = 0;
  i = 0;
  k = 0;
  j = 0;
  int tmpInt;
  while( i< nm) {
    if(wi[i] == 0){
      if(fill[k])
	while((k<n) && (fill[k] == 1))
	  k++;
      for(j=0; j<n; j++) {
	tmp->mval(j,k) = complex(z1(i,j),0.0);
      }
      fill[k] = 1;
      k++;
      i++;
    } else {
      for(j=0; j<n; j++) {
	tmp->mval(j,k) = complex(z1(i,j),z1(i+1,j));
	if(k<n/2) {
	  tmpInt = k+n/2;
	  tmp->mval(j,tmpInt) = complex(z1(i,j),-z1(i+1,j));
	} else {
	  tmp->mval(j,k+1) = complex(z1(i,j),-z1(i+1,j));
	}
      }
      fill[k] = 1;
      if(k<n/2)
	fill[k+n/2] = 1;
      else
	fill[k+1] = 1;
      k++;
      i += 2;
    }
  }

  delete []fill;
  delete []wr;
  delete []wi;
  delete []iv1;
  delete []fv1;
  delete []b;
  delete []c;

  return *tmp;

}

// operators

double& doublematrix::operator()(int i, int j) {
  return p->m[i][j];
}

double& doublematrix::operator()(int i) {
  if( p->r == 1 ) return p->m[0][i];
  if( p->c == 1 ) return p->m[i][0];
  cerr << "\n*** ERROR ***                                   \n" << endl;
  cerr << "*** ERROR *** double& doublematrix::operator()(int i) \n" << endl;
  cerr << "*** ERROR *** Can only be used with single row  \n" << endl;
  cerr << "*** ERROR *** or column matrices.               \n" << endl;
  cerr << "*** ERROR ***                                   \n" << endl;
  exit(1);
}

ostream& operator<<(ostream& os, doublematrix& x)
{
  int i,j;
  
  for(i=0; i< x.rows(); i++) {
    os << "( ";
    for(j=0; j< x.cols(); j++) {
      os << x.val(i,j) << ", ";
    }
    os << " )\n";
  }
  return os << endl;
}

void doublematrix::print(char *msg) {
  if (*msg) cout << msg << ":" << endl;
  for (int row=0; row< rows(); row++) {
    for (int col = 0; col < cols(); col++) {
      double result = mval(row, col);
      if (fabs(result) < limit) result=0;
      cout <<  result << " ";
    }
    cout << endl;
  }
}


/********************************************************
The private support functions for determinant & inverse:
*********************************************************/

// copy the from_col of mm to the to_col of "this"
void doublematrix::copy_column(doublematrix& mm, int from_col, int to_col){
  if(rows()  != mm.rows())
    error("number of rows must be equal for copy_column()");
  for(int row=0; row < rows(); row++)
    mval(row,to_col) = mm.mval(row,from_col);
}

void doublematrix::switch_columns(int col1, int col2) {
  doublematrix temp(rows());
  for(int row = 0; row < rows(); row++)
    // temporarily store col 1:
    temp.mval(row,0) = mval(row,col1);
  for(row = 0; row < rows(); row++)
    mval(row,col1) = mval(row,col2); // move col2 to col1
  for(row = 0; row < rows(); row++)
    mval(row,col2) = temp.mval(row,0); // move temp to col2
}

void doublematrix::switch_rows(int row1, int row2) {
  doublematrix temp(cols());
  for(int col = 0; col < cols(); col++)
    // temporarily store row 1:
    temp.mval(col,0) = mval(row1,col);
  for(col = 0; col < cols(); col++)
    mval(row1,col) = mval(row2,col); // move row2 to row1
  for(col = 0; col < cols(); col++)
    mval(row2,col) = temp.mval(col,0); // move temp to row2
}

// make an image of a matrix (used in L-U decomposition)
void doublematrix::deepcopy(doublematrix& from, doublematrix& to) {
  if(from.rows()  != to.rows()  || from.cols()  != to.cols())
    error("matrices must be equal dimensions for deepcopy()");
  for(int row = 0; row < from.rows(); row++) {
    for(int col = 0; col < from.cols(); col++)
      to.mval(row,col) = from.mval(row,col);
  }
}


// scale a matrix (used in L-U decomposition)
doublematrix doublematrix::scale()  {
  double temp;
  if (rows() <= 0 || cols() <= 0)
    error(" bad doublematrix for scale()");
  if (rows() != cols())
    error(" doublematrix must be square for scale()");   
  doublematrix scale_vector(cols());
  
  for (int row = 0; row < rows(); row++){
    double maximum = 0.0;
    for(int col = 0; col < cols(); col++) {
      temp = mval(row,col);
      if (fabs(temp) > maximum)
	maximum = fabs(temp); 
    }
    if(maximum == 0.0)
      error("singular doublematrix in scale()");
    // save scaling
    scale_vector.mval(row,0) = 1.0/maximum; 
  }
  return scale_vector;
}

doublematrix doublematrix::lu_decompose(int* indx, int& d ) {
/*
 Returns the L-U decomposition of a matrix. indx is an ouput
 vector which records the row permutation affected by the
 partial pivoting, d is output as +-1 depending on whether the
 number of row interchanges was even or odd, respectively.
 This routine is used in combination with lu_back_subst to
 solve linear equations or invert a matrix.
*/
  if(rows()  != cols())
    error("Matrix must be square to L-U decompose!\n");
  d = 1; // parity check
  int j,i,k,col_max; // counters
  double dum; // from the book - I don't know signj
  double sum;
  double maximum;
  doublematrix lu_decomp(rows(),cols());
  double tmp;

  // make a direct copy of the original matrix:
  deepcopy(*this,lu_decomp);
  doublematrix scale_vector = lu_decomp.scale(); // scale the matrix
  
  // The loop over columns of Crout's method:
  for(j = 0; j < rows(); j++) {
    if (j > 0) {
      // eqn 2.3.12 except for j=i:
      for (i = 0; i <= j-1; i++) {
        sum = lu_decomp.mval(i,j);
        if(i > 0) {
	  for(k = 0; k <= i-1; k++) {
	    tmp = lu_decomp.mval(i,k)*lu_decomp.mval(k,j);
	    sum -= tmp;
	    if(fabs(sum) < M_SMALL*fabs(tmp))
	      sum = 0.0;
	  }
          lu_decomp.mval(i,j) = sum;
        }
      }
    }
    
// Initialize the search for the largest pivot element:
    maximum = 0.0;
    // i=j of eq 2.3.12 & i=j+I..N of 2.3.13:
    for(i=j; i <= cols()-1; i++) {
      sum = lu_decomp.mval(i,j);
      
      if(j > 0) {
	for(k=0; k <= j-1; k++) {
	  tmp =  lu_decomp.mval(i,k) * lu_decomp.mval(k,j);
	  sum -= tmp;
	  if(fabs(sum) < M_SMALL*fabs(tmp))
	    sum = 0.0;
	}
	lu_decomp.mval(i,j) = sum;
      }
      // figure of merit for pivot:
      dum = scale_vector.mval(i,0) * fabs(sum);      
      if (dum >= maximum){
	// is it better than the best so far ?
	col_max = i;
	maximum = dum;
      }
    }
    // Do we need to interchange rows?

    if(j != col_max)  {
      lu_decomp.switch_rows(col_max,j); // Yes, do so ...
      d *= -1;  // ... and change the parity of d
      // also interchange the scale factor:
      dum = scale_vector.mval(col_max,0);
      scale_vector.mval(col_max,0) = scale_vector.mval(j,0);
      scale_vector.mval(j,0) = dum;
    }

    indx[j] = col_max;
    // Now, finally, divide by the pivot element:
    if(j != rows() -1){
      // If the pivot element is .zero the matrix is
      // singular (at least to the precision of the
      // algorithm). For some applications on singular
      // matrices, it is desirable to substitute tiny for 0
      if(lu_decomp.mval(j,j) == 0)
	lu_decomp.mval(j,j) = tiny;
      dum = 1/lu_decomp.mval(j,j);
      for(i=j+1; i <= cols()-1; i++)
	lu_decomp.mval(i,j) *= dum;
    }

  }
  if(lu_decomp.mval(rows()-1,cols()-1) == 0)
    lu_decomp.mval(rows()-1,cols()-1) = tiny;

  return lu_decomp;
}


void doublematrix::lu_back_subst(int* indx, doublematrix& b)  {

// Solves the set of N linear equations A*X = B. Here "this"
// is the LU-decomposition of the matrix A, determined by the
// routine lu_decompose(). Indx is input as the permutation
// vector returned  by lu_decompose(). B is input as the
// right-hand side vector B, and returns with the solution
// vector X. This routine takes into  account the possibility
// that B will begin with many zero elements, so it is efficient
// for use in matrix inversion. See pp 36-37 in
// Press & Flannery.
  
  if(rows()  != cols())
    error ("non-square lu_decomp matrix in lu_back_subst()");
  if(rows()  != b.rows()) 
    error("wrong size B vector passed to lu_back_subst()");
//  if(rows()  != indx.rows())
//    error("wrong size indx vector passed to lu_back_subst()");
  int i,j,ip;
  int ii = -1;
  double sum;
  double tmp;

  for(i=0;i < rows(); i++){
    ip = indx[i];
    sum = b.mval(ip,0);
    b.mval(ip,0) = b.mval(i,0);
    if (ii >= 0) {
      for(j = ii; j <= i-1; j++) {
	tmp = mval(i,j) * b.mval(j,0);
	sum -= tmp;
	if(fabs(sum) < M_SMALL*fabs(tmp))
	  sum = 0.0;
      }
    } else if(sum != 0)
      ii = i;
    b.mval(i,0) = sum;
  }
  for(i = rows() -1; i >= 0; i--){
    sum = b.mval(i,0);
    if (i < cols() -1) {
      for (j = i + 1; j <= rows()-1; j++) {
	tmp = mval(i,j) * b.mval(j,0);
	sum -= tmp;
	if(fabs(sum) < M_SMALL*fabs(tmp))
	  sum = 0.0;
      }
    }
    // store a component of the soln vector X:
    b.mval(i,0) = sum/mval(i,i);
  }
}

// Now do complexmatrix methods.

void complexmatrix::error(char * msg1, char * msg2) const{
  cerr << "complexmatrix error: " << msg1 << " " << msg2  << endl;
  exit(1);
}

complexmatrix::complexmatrix(int mrows, int columns, complex initval) {
  // create the structure:
  p = new complexmatrep;
  p->r = mrows;
  p->c = columns;
  // allocate memory for the actual matrix:
  p->m = new complex *[mrows];
  for (int x = 0; x < mrows; x++)
    p->m[x] = new complex[columns];
  p->n = 1;  // so far, there's one reference to this data
  for (int i=0; i< mrows; i++) {
    for (int j = 0; j < columns; j++)
      mval(i,j) = initval;
  }
}

complexmatrix::complexmatrix(int mrows, int columns, complex * initvalues){
  // create the structure:
  p    = new complexmatrep;
  p->r = mrows;
  p->c = columns;
  // allocate memory for the actual matrix:
  p->m = new complex *[mrows];
  for (int x = 0; x < mrows; x++)
    p->m[x] = new complex[columns];
  p->n = 1;  // so far, there's one reference to this data:
  int c = 0;
  for (int i=0; i< mrows; i++) {
    for (int j = 0; j < columns; j++)
      mval(i,j) =  initvalues[c++];
  }
}

// create an identity or symplectic matrix:
complexmatrix::complexmatrix(char * flag, int dimension) {
  int i,j;
  float tmp_float = dimension/2.0;
  
  if(((dimension/2) != tmp_float) && (flag[0] == 'J'))
    error("dimension must be even","for J matrix");
  if (flag[0]  == 'I') {
    p = new complexmatrep;
    p->r = dimension;
    p->c = dimension;
    p->m = new complex *[dimension];
    for (int x = 0; x < dimension; x++)
      p->m[x] = new complex[dimension];
    p->n = 1;
    for (i = 0; i< dimension; i++) {
      for (j = 0; j < dimension; j++)
	mval(i,j) = (i == j ? complex(1.0,0.0) : complex_zero);
    }
  } else if (flag[0] == 'J') {
    p = new complexmatrep;
    p->r = dimension;
    p->c = dimension;
    p->m = new complex *[dimension];
    for (int x = 0; x < dimension; x++)
      p->m[x] = new complex[dimension];
    p->n = 1;
    for ( i = 0; i< dimension; i++) {
      for ( j = 0; j < dimension; j++)
	mval(i,j) = complex_zero;
    }
    for (i = dimension/2; i< dimension; i++) {
      mval(i-dimension/2,i) = 1;
      mval(i,i-dimension/2) = -1;
    }
  } else
    error("Unknown flag: ",flag);
}

complexmatrix::complexmatrix(int mrows, int columns, double * initvalues){
  // create the structure:
  p    = new complexmatrep;
  p->r = mrows;
  p->c = columns;
  // allocate memory for the actual matrix:
  p->m = new complex *[mrows];
  for (int x = 0; x < mrows; x++)
    p->m[x] = new complex[columns];
  p->n = 1;  // so far, there is one reference to this data:
  int c = 0;
  for (int i=0; i< mrows; i++) {
    for (int j = 0; j < columns; j++)
      mval(i,j) =  complex(initvalues[c++],0.0);
  }
}

complexmatrix::complexmatrix(doublematrix& x) {
  p = new complexmatrep;
  p->r = x.p->r;
  p->c = x.p->c;
  // allocate memory for the actual matrix:
  p->m = new complex *[x.p->r];
  for (int i = 0; i < x.p->r; i++)
    p->m[i] = new complex[x.p->c];
  p->n = 1;  // so far, there is one reference to this data
  for(int row = 0; row < x.rows(); row++) {
    for(int col = 0; col < x.cols(); col++)
      mval(row,col) = complex(x.mval(row,col),0.0);
  }
}

// read from "standard" matrix file:
complexmatrix::complexmatrix(char * initfile) {
  ifstream from(initfile);
  if(!from)
    error("cannot open complexmatrix initializer file", initfile);
  const bsize = 120;
  char buf[bsize], *cp, *cp2;
  int rfound = 0, cfound = 0, colonsfound = 0;
  p = new complexmatrep;
  /* Parse file initialization header  */
  while(from.getline(buf, bsize)) { // for each header line
    // Remove comments with ANSI C library function "strprk":
    if( ( cp = strpbrk(buf,"#"))  != 0 ) // look for comments
      *cp = '\0';  // terminate string at comment
    if( ( cp = strpbrk(buf,"r") )  != 0 )
      if ( ( cp2 = strpbrk(cp,  "="))  !=0)
	if ( ( cp = strpbrk(cp2,  "0123456789"))  != 0 ) {
	  p->r = atoi(cp);
	  rfound++;  // flag to say rows were found
	}
    if( ( cp = strpbrk(buf,"c") )  != 0)
      if ( ( cp2 = strpbrk(cp,  "="))  != 0)
	if ( ( cp = strpbrk(cp2,  "0123456789"))  != 0 ) {
	  p->c = atoi(cp);
	  cfound++;  // flag to say cols were found
	}
    if ( strstr(buf,":::")  != 0 )  {
      colonsfound++;
      break; // ... out of "while" loop
    }
  }
  if ( !rfound || !cfound || !colonsfound )  {
    cerr << initfile << nonstandard;
    exit(1);
  }
  p->m = new complex *[p->r];
  for (int x = 0; x < p->r; x++)
    p->m[x] = new complex[p->c];
  p->n = 1;  // so far, there's one reference to this data
  for (int row = 0; row < p->r; row++)  {
    for(int col = 0; col < p->c; col++)  {
      const nbsz = 20; char nb[nbsz];
      from >> nb; // scan for space-delimited string
      mval(row,col) = atof(nb); // convert int to a double
      if(from.bad())
	error ( "problem with complexmatrix initializer file", initfile);
    }
  }
}

complexmatrix::complexmatrix(complexmatrix& x) {
  p = new complexmatrep;
  p->r = x.p->r;
  p->c = x.p->c;
  // allocate memory for the actual matrix:
  p->m = new complex *[x.p->r];
  for (int i = 0; i < x.p->r; i++)
    p->m[i] = new complex[x.p->c];
  p->n = 1;  // so far, there's one reference to this data
  deepcopy(x,*this);
}

complexmatrix::~complexmatrix()  {
  if (--p->n == 0) { // if reference count goes to 0
    for (int x = 0; x < rows(); x++)
      delete []p->m[x];
    delete []p->m; // delete data
    delete p;
  }
}

complexmatrix complexmatrix::operator=(complexmatrix& x) {

  if(this == &x) return *this;

  p = new complexmatrep;
  p->r = x.p->r;
  p->c = x.p->c;
  // allocate memory for the actual matrix:
  p->m = new complex *[x.p->r];
  for (int i = 0; i < x.p->r; i++)
    p->m[i] = new complex[x.p->c];
  p->n = 1;  // so far, there's one reference to this data
  deepcopy(x,*this);

  return *this;
}

complexmatrix complexmatrix::operator-() {
  complexmatrix unaryminus(rows(),cols());
  for (int i=0; i< rows(); i++) {
    for (int j = 0; j < cols(); j++)
      unaryminus.mval(i,j) = -mval(i,j);
  }
  return unaryminus;
}

void complexmatrix::write_standard(char * filename, char * msg) {
  ofstream to(filename);
  if(!to)
    error("cannot open or create complexmatrix output file", filename);
  to << "# " << filename
     << ": complexmatrix file written in \"standard\" format" <<endl;
//  time_t clock;
//  time(&clock);
//  to <<  "# " << asctime(localtime(&clock));
  to <<  "# " << msg << endl;
  to << "rows= " << rows() << " columns= " << cols() << endl;
  to <<  ":::" << endl;
  for (int row = 0; row < rows(); row++) {
    for(int col = 0; col < cols(); col++)
      to << mval(row,col) << "  ";
    to << endl;
  }
}

complexmatrix operator+(const complexmatrix& x, const complexmatrix& y)  {
  if(( x.rows()  != y.rows())  ||  ( x.cols()  != y.cols()))
    x.error("must have equal dimensions for addition!");
  complexmatrix* sum = new complexmatrix(x.rows(),x.cols());
  for (int i=0; i< x.rows(); i++) {
    for (int j = 0; j < x.cols(); j++) {
      sum->p->m[i][j] = x.p->m[i][j] + y.p->m[i][j];
      if( abs(sum->p->m[i][j]) < M_SMALL*abs(y.p->m[i][j])  ) 
	sum->p->m[i][j] = 0.0;
    }
  }
  return *sum; // see note for operator*()
}

complexmatrix operator+(const complexmatrix& x, const complex arg)  {
  complexmatrix* sum = new complexmatrix(x.rows(),x.cols());
  for (int i=0; i< x.rows(); i++) {
    for (int j = 0; j < x.cols(); j++) {
      sum->p->m[i][j] = x.p->m[i][j] + arg;
      if( abs(sum->p->m[i][j]) < M_SMALL*abs(arg)  ) 
	sum->p->m[i][j] = 0.0;
    }
  }
  return *sum; // see note for operator*()
}

complexmatrix operator+( const complex arg,const complexmatrix& x)  {
  complexmatrix* sum = new complexmatrix(x.rows(),x.cols());
  for (int i=0; i< x.rows(); i++) {
    for (int j = 0; j < x.cols(); j++) {
      sum->p->m[i][j] = x.p->m[i][j] + arg;
      if( abs(sum->p->m[i][j]) < M_SMALL*abs(arg)  ) 
	sum->p->m[i][j] = 0.0;
    }
  }
  return *sum; // see note for operator*()
}

complexmatrix operator-(const complexmatrix& x,const complexmatrix& y)  {
  if(( x.rows()  != y.rows())  || ( x.cols()  != y.cols()))
    x.error("must have equal dimensions for subtraction!");
  complexmatrix* sum = new complexmatrix(x.rows(),x.cols());
  for (int i=0; i< x.rows(); i++) {
    for (int j = 0; j < x.cols(); j++) {
      sum->p->m[i][j] = x.p->m[i] [j] - y.p->m[i] [j];
      if( abs(sum->p->m[i][j]) < M_SMALL*abs(y.p->m[i][j])  ) 
	sum->p->m[i][j] = 0.0;
    }
  }
  return *sum; // see note for operator*()
}

complexmatrix operator-(const complexmatrix& x, const complex arg)  {
  complexmatrix* sum = new complexmatrix(x.rows(),x.cols());
  for (int i=0; i < x.rows(); i++) {
    for (int j = 0; j < x.cols(); j++) {
      sum->p->m[i][j] = x.p->m[i][j] - arg;
      if( abs(sum->p->m[i][j]) < M_SMALL*abs(arg)  ) 
	sum->p->m[i][j] = 0.0;
    }
  }
  return *sum; // see note for operator*()
}

complexmatrix operator-( const complex arg, const complexmatrix& x)  {
  complexmatrix* sum = new complexmatrix(x.rows(),x.cols());
  for (int i=0; i < x.rows(); i++) {
    for (int j = 0; j < x.cols(); j++) {
      sum->p->m[i][j] = x.p->m[i][j] - arg;
      if( abs(sum->p->m[i][j]) < M_SMALL*abs(arg)  ) 
	sum->p->m[i][j] = 0.0;
    }
  }
  return *sum; // see note for operator*()
}


complexmatrix operator*(const complexmatrix& x, const complexmatrix& y)  {
  complex tmp;
  if( x.cols()  != y.rows())
    x.error("# rows of second mat must equal "
	  "# cols of first for multiply!");
  complexmatrix* result =  new complexmatrix(x.rows(),x.cols());
  for(int row = 0; row < x.rows(); row++) {
    for(int col = 0; col < x.cols(); col++){
      complex sum = 0;
      for(int i = 0; i < x.cols(); i++) {
	tmp = x.p->m[row][i] * y.p->m[i][col];
	sum += tmp;
	if(abs(sum) < M_SMALL*abs(tmp))
	  sum = 0.0;
      }
      result->p->m[row][col] = sum;
    }
  }
  return *result;
}

complexmatrix operator*(const complexmatrix& x,const complex arg) {
  complexmatrix* result = new complexmatrix(x.rows(),x.cols());
  for (int i=0; i< x.rows(); i++) {
    for (int j = 0; j < x.cols(); j++)
      result->p->m[i][j] = x.p->m[i][j] * arg;
  }
  return *result;
}

complexmatrix operator*(const complex arg,const complexmatrix& x) {
  complexmatrix* result = new complexmatrix(x.rows(),x.cols());
  for (int i=0; i< x.rows(); i++) {
    for (int j = 0; j < x.cols(); j++)
      result->p->m[i][j] = x.p->m[i][j] * arg;
  }
  return *result;
}

complexmatrix operator*(const complexmatrix& x, const doublematrix& y)  {
  complex tmp;
  if( x.cols()  != y.rows())
    x.error("# rows of second mat must equal "
	  "# cols of first for multiply!");
  complexmatrix* result =  new complexmatrix(x.rows(),y.cols());
  for(int row = 0; row < x.rows(); row++) {
    for(int col = 0; col < y.cols(); col++){
      complex sum = 0;
      for(int i = 0; i < x.cols(); i++) {
	tmp = x.p->m[row][i] * complex(y.p->m[i][col],0.0);
	sum += tmp;
	if(abs(sum) < M_SMALL*abs(tmp))
	  sum = 0.0;
      }
      result->p->m[row][col] = sum;
    }
  }
  return *result;
}

complexmatrix operator*( const doublematrix& x,const complexmatrix& y)  {
  complex tmp;
  if( x.cols()  != y.rows())
    x.error("# rows of second mat must equal "
	  "# cols of first for multiply!");
  complexmatrix* result =  new complexmatrix(x.rows(),y.cols());
  for(int row = 0; row < x.rows(); row++) {
    for(int col = 0; col < y.cols(); col++){
      complex sum = 0;
      for(int i = 0; i < x.cols(); i++) {
	tmp = complex(x.p->m[row][i],0.0) * y.p->m[i][col];
	sum += tmp;
	if(abs(sum) < M_SMALL*abs(tmp))
	  sum = 0.0;
      }
      result->p->m[row][col] = sum;
    }
  }
  return *result;
}

complex & complexmatrix::val(int row, int col)  {
  //  if (row >= 0 && row < rows() && col >= 0 && col << cols())
  if (row < 0 || row >= rows()  || col < 0 || col >= cols())
    error("index out of range");
  return (mval(row,col));
}

complexmatrix complexmatrix::transpose() {
  if(rows()  != cols())
    error("complexmatrix must be square to transpose!\n");
  complexmatrix trans(rows(),cols());
  for (int row = 0; row < rows(); row++)  {
    for(int col = 0; col < cols(); col++)
      trans.mval(col,row) = mval(row,col);
  }
  return trans;
}

complex complexmatrix::mmin()  {
  complex temp;
  if(rows() <= 0 || cols() <= 0)
    error("bad complexmatrix size for min()");
  complex minimum = mval(0,0);
  for (int row = 0; row < rows(); row++)  {
    for(int col = 0; col < cols(); col++)
      if (abs((temp = mval(row,col))) < abs(minimum))
	minimum = temp;
  }
  return minimum;
}

complex complexmatrix::mmax()  {
  complex temp;
  if(rows() <= 0 || cols() <= 0)
    error("bad complexmatrix size for max()");
  complex maximum = mval(0,0);
  for (int row = 0; row < rows(); row++)  {
    for(int col = 0; col < cols(); col++){
      if (abs((temp = mval(row,col))) > abs(maximum))
	maximum = temp;
    }
  }
  return maximum;
}

complexmatrix& complexmatrix::eigenValues() {
  int nm = rows();
  int n  = cols();
  int matz = 0;
  double* wr = new double[cols()];
  double* wi = new double[cols()];
  double zr,zi;
  double* fv1 = new double[cols()];
  double* fv2 = new double[cols()];
  double* fv3 = new double[cols()];
  double* br = new double[rows()*cols()];
  double* bi = new double[rows()*cols()];
  int ierr;
  complexmatrix* tmp = new complexmatrix(1,cols());

  int k = 0;
  for(int i = 0; i < rows(); i++) {
    for(int j = 0; j< cols(); j++) {
      br[k] = real(mval(j,i));  	// the cg_ routine uses the transpose
      bi[k] = imag(mval(j,i));  	// the cg_ routine uses the transpose
      k++;
    }
  }

  cg_(&nm,&n,br,bi,wr,wi,&matz,&zr,&zi,fv1,fv2,fv3,&ierr);
  if(ierr != 0) {
    cerr << "complexmatrix: Error in eigenvalue routine. error = ";
    cerr << ierr<< end;
    delete []wr;
    delete []wi;
    delete []fv1;
    delete []fv2;
    delete []fv3;
    delete []br;
    delete []bi;
    return *tmp;
  }
  for(i = 0; i< cols(); i++)
    tmp->mval(0,i) = complex(wr[i],wi[i]);

  delete []wr;
  delete []wi;
  delete []fv1;
  delete []fv2;
  delete []fv3;
  delete []br;
  delete []bi;

  return *tmp;
}

complexmatrix& complexmatrix::eigenVectors() {
  int nm = rows();
  int n  = cols();
  int matz = 1;
  double* wr = new double[cols()];
  double* wi = new double[cols()];
  double* fv1 = new double[cols()];
  double* fv2 = new double[cols()];
  double* fv3 = new double[cols()];
  double* br = new double[rows()*cols()];
  double* bi = new double[rows()*cols()];
  double* cr = new double[rows()*cols()];
  double* ci = new double[rows()*cols()];
  int ierr,i,j,k;
  complexmatrix* tmp = new complexmatrix(rows(),cols());

  k = 0;
  for(i = 0; i < rows(); i++) {
    for(j = 0; j< cols(); j++) {
      br[k] = real(mval(j,i));	// the cg_ routine uses the transpose
      bi[k] = imag(mval(j,i));		// the cg_ routine uses the transpose
      k++;
    }
  }
  cg_(&nm,&n,br,bi,wr,wi,&matz,cr,ci,fv1,fv2,fv3,&ierr);
  if(ierr != 0) {
    cerr << "complexmatrix: Error in eigenvector routine. error = ";
    cerr << ierr<< end;
    delete []wr;
    delete []wi;
    delete []fv1;
    delete []fv2;
    delete []fv3;
    delete []br;
    delete []bi;
    delete []cr;
    delete []ci;
    return *tmp;
  }
  complexmatrix z1(nm,n);
  k = 0;
  for(i = 0; i < rows(); i++) {
    for(j = 0; j< cols(); j++) {
      z1(i,j) = complex(cr[k],ci[k]);
      k++;
    }
  }
  *tmp = z1;
  delete []wr;
  delete []wi;
  delete []fv1;
  delete []fv2;
  delete []fv3;
  delete []br;
  delete []bi;
  delete []cr;
  delete []ci;
  
  return *tmp;
}


complex complexmatrix::determinant()  {
  if(rows()  != cols())
    error("complexmatrix must be square for determinant()");
  int* indx = new int[cols()]; // create the "index vecto
  complexmatrix B(cols());    // see pp 38. in Numerical Recipes
  int d;
  // perform the decomposition once:
  complexmatrix decomp = lu_decompose(indx,d);
  complex determinant = d;
  for(int i=0; i < cols() ; i++)
    determinant *= decomp.mval(i,i);
  delete []indx;
  return determinant;
}

complexmatrix complexmatrix::inverse() {
  if(rows()  != cols())
    error("complexmatrix must be square for inverse()");
  complexmatrix Y("I",rows()); // create an identity complexmatrix
  int* indx = new int[cols()];  // create the "index vector"
  complexmatrix B(cols());     // see Press & Flannery
  int d;
  // perform the decomposition once:
  complexmatrix decomp = lu_decompose(indx,d);
  for(int col = 0; col < cols(); col++){
    B.copy_column(Y,col,0);
    decomp.lu_back_subst(indx,B);
    Y.copy_column(B,0,col);
  }
  delete []indx;
  return Y;
}

complex& complexmatrix::operator()(int i, int j) {
  return p->m[i][j];
}

complex& complexmatrix::operator()(int i) {
  if( p->r == 1 ) return p->m[0][i];
  if( p->c == 1 ) return p->m[i][0];
  cerr << "\n*** ERROR ***                                   \n" << endl;
  cerr << "*** ERROR *** complex& complexmatrix::operator()(int i) \n" << endl;
  cerr << "*** ERROR *** Can only be used with single row  \n" << endl;
  cerr << "*** ERROR *** or column matrices.               \n" << endl;
  cerr << "*** ERROR ***                                   \n" << endl;
  exit(1);
}

ostream& operator<<(ostream& os, complexmatrix& x)
{
  int i,j;
  
  for(i=0; i< x.rows(); i++) {
    os << "( ";
    for(j=0; j< x.cols(); j++) {
      os << x.val(i,j) << ", ";
    }
    os << " )\n";
  }
  return os << endl;
}

void complexmatrix::print(char *msg) {
  if (*msg) cout << msg << ":" << endl;
  for (int row=0; row< rows(); row++) {
    for (int col = 0; col < cols(); col++) {
      complex result = mval(row, col);
      if (abs(result) < limit) result=0;
      cout << result << " ";
    }
    cout << endl;
  }
}


/********************************************************
The private support functions for determinant & inverse:
*********************************************************/

// copy the from_col of mm to the to_col of "this"
void complexmatrix::copy_column(complexmatrix& mm, int from_col, int to_col){
  if(rows()  != mm.rows())
    error("number of rows must be equal for copy_column()");
  for(int row=0; row < rows(); row++)
    mval(row,to_col) = mm.mval(row,from_col);
}

void complexmatrix::switch_columns(int col1, int col2) {
  complexmatrix temp(rows());
  for(int row = 0; row < rows(); row++)
    // temporarily store col 1:
    temp.mval(row,0) = mval(row,col1);
  for(row = 0; row < rows(); row++)
    mval(row,col1) = mval(row,col2); // move col2 to col1
  for(row = 0; row < rows(); row++)
    mval(row,col2) = temp.mval(row,0); // move temp to col2
}

void complexmatrix::switch_rows(int row1, int row2) {
  complexmatrix temp(cols());
  for(int col = 0; col < cols(); col++)
    // temporarily store row 1:
    temp.mval(col,0) = mval(row1,col);
  for(col = 0; col < cols(); col++)
    mval(row1,col) = mval(row2,col); // move row2 to row1
  for(col = 0; col < cols(); col++)
    mval(row2,col) = temp.mval(col,0); // move temp to row2
}

// make an image of a matrix (used in L-U decomposition)
void complexmatrix::deepcopy(complexmatrix& from, complexmatrix& to) {
  if(from.rows()  != to.rows()  || from.cols()  != to.cols())
    error("matrices must be equal dimensions for deepcopy()");
  for(int row = 0; row < from.rows(); row++) {
    for(int col = 0; col < from.cols(); col++)
      to.mval(row,col) = from.mval(row,col);
  }
}

// scale a matrix (used in L-U decomposition)
complexmatrix complexmatrix::scale()  {
  complex temp;
  if (rows() <= 0 || cols() <= 0)
    error(" bad complexmatrix for scale()");
  if (rows() != cols())
    error(" complexmatrix must be square for scale()");   
  complexmatrix scale_vector(cols());
  
  for (int row = 0; row < rows(); row++){
    double maximum = 0.0;
    for(int col = 0; col < cols(); col++) {
      temp = mval(row,col);
      if (abs(temp)  > maximum)
	maximum = abs(temp); 
    }
    if(maximum == complex_zero)
      error("singular complexmatrix in scale()");
    // save scaling
    scale_vector.mval(row,0) = complex(1.0,0.0)/maximum; 
  }
  return scale_vector;
}

complexmatrix complexmatrix::lu_decompose(int* indx, int& d ) {
/*
 Returns the L-U decomposition of a matrix. indx is an ouput
 vector which records the row permutation affected by the
 partial pivoting, d is output as +-1 depending on whether the
 number of row interchanges was even or odd, respectively.
 This routine is used in combination with lu_back_subst to
 solve linear equations or invert a matrix.
*/
  if(rows()  != cols())
    error("Matrix must be square to L-U decompose!\n");
  d = 1; // parity check
  int j,i,k,col_max; // counters
  complex dum; // from the book - I don't know signj
  complex sum;
  complex maximum;
  complexmatrix lu_decomp(rows(),cols());
  complex tmp;

  // make a direct copy of the original matrix:
  deepcopy(*this,lu_decomp);
  complexmatrix scale_vector = lu_decomp.scale(); // scale the matrix

  // The loop over columns of Crout's method:
  for(j = 0; j < rows(); j++) {
    if (j > 0) {
      // eqn 2.3.12 except for j=i:
      for (i = 0; i <= j-1; i++) {
        sum = lu_decomp.mval(i,j);
        if(i > 0) {
	  for(k = 0; k <= i-1; k++) {
	    tmp = lu_decomp.mval(i,k)*lu_decomp.mval(k,j);
	    sum -= tmp;
	    if(abs(sum) < M_SMALL*abs(tmp))
	      sum = complex_zero;
	  }
          lu_decomp.mval(i,j) = sum;
        }
      }
    }

// Initialize the search for the largest pivot element:
    maximum = 0.0;
    // i=j of eq 2.3.12 & i=j+I..N of 2.3.13:
    for(i=j; i <= cols()-1; i++) {
      sum = lu_decomp.mval(i,j);
      if(j > 0) {
	for(k=0; k <= j-1; k++) {
	  tmp =  lu_decomp.mval(i,k) * lu_decomp.mval(k,j);
	  sum -= tmp;
	  if(abs(sum) < M_SMALL*abs(tmp))
	    sum = complex_zero;
	  }
	lu_decomp.mval(i,j) = sum;
      }
      // figure of merit for pivot:
      dum = scale_vector.mval(i,0) * sum;
      if (abs(dum) >= abs(maximum)){
	// is it better than the best so far ?
	col_max = i;
	maximum = dum;
      }
    }
    // Do we need to interchange rows?

    if(j != col_max)  {
      lu_decomp.switch_rows(col_max,j); // Yes, do so ...
      d *= -1;  // ... and change the parity of d
      // also interchange the scale factor:
      dum = scale_vector.mval(col_max,0);
      scale_vector.mval(col_max,0) = scale_vector.mval(j,0);
      scale_vector.mval(j,0) = dum;
    }
    indx[j] = col_max;
    // Now, finally, divide by the pivot element:
    if(j != rows() -1){
      // If the pivot element is .zero the matrix is
      // singular (at least to the precision of the
      // algorithm). For some applications on singular
      // matrices, it is desirable to substitute tiny for 0
      if(lu_decomp.mval(j,j) == complex_zero)
	lu_decomp.mval(j,j) = complex(tiny,0.0);
      dum = complex(1.0,0.0)/lu_decomp.mval(j,j);
      for(i=j+1; i <= cols()-1; i++)
	lu_decomp.mval(i,j) *= dum;
    }
  }
  if(lu_decomp.mval(rows()-1,cols()-1) == complex_zero)
    lu_decomp.mval(rows()-1,cols()-1) = complex(tiny,0.0);

  return lu_decomp;
}

void complexmatrix::lu_back_subst(int* indx, complexmatrix& b)  {
// Solves the set of N linear equations A*X = B. Here "this"
// is the LU-decomposition of the matrix A, determined by the
// routine lu_decompose(). Indx is input as the permutation
// vector returned  by lu_decompose(). B is input as the
// right-hand side vector B, and returns with the solution
// vector X. This routine takes into  account the possibility
// that B will begin with many zero elements, so it is efficient
// for use in matrix inversion. See pp 36-37 in
// Press & Flannery.

  if(rows()  != cols())
    error ("non-square lu_decomp matrix in lu_back_subst()");
  if(rows()  != b.rows()) 
    error("wrong size B vector passed to lu_back_subst()");
//  if(rows()  != indx.rows())
//    error("wrong size indx vector passed to lu_back_subst()");
  int i,j,ip;
  int ii = -1;
  complex sum;
  complex tmp;
  for(i=0;i < rows(); i++){
    ip = indx[i];
    sum = b.mval(ip,0);
    b.mval(ip,0) = b.mval(i,0);
    if (ii >= 0) {
      for(j = ii; j <= i-1; j++) {
	tmp = mval(i,j) * b.mval(j,0);
	sum -= tmp;
	if(abs(sum) < M_SMALL*abs(tmp))
	  sum = complex_zero;
      }
    } else if(sum != 0)
      ii = i;
    b.mval(i,0) = sum;
  }
  for(i = rows() -1; i >= 0; i--){
    sum = b.mval(i,0);
    if (i < cols() -1)
      for (j = i + 1; j <= rows()-1; j++) {
	tmp = mval(i,j) * b.mval(j,0);
	sum -= tmp;
	if(abs(sum) < M_SMALL*abs(tmp))
	  sum = complex_zero;
      }
    // store a component of the soln vector X:
    b.mval(i,0) = sum/mval(i,i);
  }
}




End of matrix.cc
