From bug-octave-request at bevo dot che dot wisc dot edu  Tue Jan 20 15:41:01 1998
Subject: Re: c2d.m/expm
From: "A. Scottedward Hodel" <A dot S dot Hodel at Eng dot Auburn dot EDU>
To: schoenh at topas dot flo dot sh dot bosch dot de
Cc: bug-octave at bevo dot che dot wisc dot edu
Date: Tue, 20 Jan 1998 15:40:50 -0600

----------
X-Sun-Data-Type: text
X-Sun-Data-Description: text
X-Sun-Data-Name: text
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 17

Thanks for the c2d/expm example.  Problem and solution are as follows:

Problem: the trace normalization step in Ward's 1977 algorithm for
         computing the matrix exponential causes problems when computing
         the matrix exponential for stiff systems.
Solution: only apply trace normalization when it will not result in
         a destabilizing shift  A := A - sI (s = trace(A)/rows(a)).

This works for your example.

John will release updated code in octave 2.0.10.  For now, here's
updated files from octave2.0.9; these should be placed in
octave-2.0.9/liboctave.  I've include my RCS directory for these two
files in case you need to build patches.

A S Hodel Dept. of Elect Eng, 200 Broun Hall Auburn Univ, AL 36849 (334)8441854 
FAX:-1809       http://www.eng.auburn.edu/users/scotte    scotte at eng dot auburn dot edu 
----------
X-Sun-Data-Type: default
X-Sun-Data-Description: default
X-Sun-Data-Name: CMatrix.cc
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 3851

// Matrix manipulations.
/*

Copyright (C) 1996 John W. Eaton

This file is part of Octave.

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

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

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

*/

#if defined (__GNUG__)
#pragma implementation
#endif

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <cfloat>

#include <iostream.h>

// XXX FIXME XXX
#ifdef HAVE_SYS_TYPES_H
#include <sys/types.h>
#endif

#include "CmplxAEPBAL.h"
#include "CmplxDET.h"
#include "CmplxSCHUR.h"
#include "CmplxSVD.h"
#include "f77-fcn.h"
#include "lo-error.h"
#include "lo-ieee.h"
#include "lo-mappers.h"
#include "lo-utils.h"
#include "mx-base.h"
#include "mx-inlines.cc"
#include "oct-cmplx.h"

// Fortran functions we call.

extern "C"
{
  int F77_FCN (zgemm, ZGEMM) (const char*, const char*, const int&,
			      const int&, const int&, const Complex&,
			      const Complex*, const int&,
			      const Complex*, const int&,
			      const Complex&, Complex*, const int&, 
			      long, long);

  int F77_FCN (zgeco, ZGECO) (Complex*, const int&, const int&, int*,
			      double&, Complex*);

  int F77_FCN (zgedi, ZGEDI) (Complex*, const int&, const int&, int*,
			      Complex*, Complex*, const int&);

  int F77_FCN (zgesl, ZGESL) (Complex*, const int&, const int&, int*,
			      Complex*, const int&);

  int F77_FCN (zgelss, ZGELSS) (const int&, const int&, const int&,
				Complex*, const int&, Complex*,
				const int&, double*, double&, int&,
				Complex*, const int&, double*, int&);

  // Note that the original complex fft routines were not written for
  // double complex arguments.  They have been modified by adding an
  // implicit double precision (a-h,o-z) statement at the beginning of
  // each subroutine.

  int F77_FCN (cffti, CFFTI) (const int&, Complex*);

  int F77_FCN (cfftf, CFFTF) (const int&, Complex*, Complex*);

  int F77_FCN (cfftb, CFFTB) (const int&, Complex*, Complex*);

  int F77_FCN (zlartg, ZLARTG) (const Complex&, const Complex&,
				double&, Complex&, Complex&);

  int F77_FCN (ztrsyl, ZTRSYL) (const char*, const char*, const int&,
				const int&, const int&,
				const Complex*, const int&,
				const Complex*, const int&, 
				const Complex*, const int&, double&,
				int&, long, long);

  double F77_FCN (zlange, ZLANGE) (const char*, const int&,
				   const int&, const Complex*,
				   const int&, double*); 
}

static const Complex Complex_NaN_result (octave_NaN, octave_NaN);

// Complex Matrix class

ComplexMatrix::ComplexMatrix (const Matrix& a)
  : MArray2<Complex> (a.rows (), a.cols ())
{
  for (int j = 0; j < cols (); j++)
    for (int i = 0; i < rows (); i++)
      elem (i, j) = a.elem (i, j);
}

ComplexMatrix::ComplexMatrix (const RowVector& rv)
  : MArray2<Complex> (1, rv.length (), 0.0)
{
  for (int i = 0; i < rv.length (); i++)
    elem (0, i) = rv.elem (i);
}

ComplexMatrix::ComplexMatrix (const ColumnVector& cv)
  : MArray2<Complex> (cv.length (), 1, 0.0)
{
  for (int i = 0; i < cv.length (); i++)
    elem (i, 0) = cv.elem (i);
}

ComplexMatrix::ComplexMatrix (const DiagMatrix& a)
  : MArray2<Complex> (a.rows (), a.cols (), 0.0)
{
  for (int i = 0; i < a.length (); i++)
    elem (i, i) = a.elem (i, i);
}

ComplexMatrix::ComplexMatrix (const ComplexRowVector& rv)
  : MArray2<Complex> (1, rv.length (), 0.0)
{
  for (int i = 0; i < rv.length (); i++)
    elem (0, i) = rv.elem (i);
}

ComplexMatrix::ComplexMatrix (const ComplexColumnVector& cv)
  : MArray2<Complex> (cv.length (), 1, 0.0)
{
  for (int i = 0; i < cv.length (); i++)
    elem (i, 0) = cv.elem (i);
}

ComplexMatrix::ComplexMatrix (const ComplexDiagMatrix& a)
  : MArray2<Complex> (a.rows (), a.cols (), 0.0)
{
  for (int i = 0; i < a.length (); i++)
    elem (i, i) = a.elem (i, i);
}

// XXX FIXME XXX -- could we use a templated mixed-type copy function
// here?

ComplexMatrix::ComplexMatrix (const charMatrix& a)
{
  for (int i = 0; i < a.cols (); i++)
    for (int j = 0; j < a.rows (); j++)
      elem (i, j) = a.elem (i, j);
}

bool
ComplexMatrix::operator == (const ComplexMatrix& a) const
{
  if (rows () != a.rows () || cols () != a.cols ())
    return false;

  return equal (data (), a.data (), length ());
}

bool
ComplexMatrix::operator != (const ComplexMatrix& a) const
{
  return !(*this == a);
}

// destructive insert/delete/reorder operations

ComplexMatrix&
ComplexMatrix::insert (const Matrix& a, int r, int c)
{
  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
    {
      (*current_liboctave_error_handler) ("range error for insert");
      return *this;
    }

  for (int j = 0; j < a_nc; j++)
    for (int i = 0; i < a_nr; i++)
      elem (r+i, c+j) = a.elem (i, j);

  return *this;
}

ComplexMatrix&
ComplexMatrix::insert (const RowVector& a, int r, int c)
{
  int a_len = a.length ();
  if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
    {
      (*current_liboctave_error_handler) ("range error for insert");
      return *this;
    }

  for (int i = 0; i < a_len; i++)
    elem (r, c+i) = a.elem (i);

  return *this;
}

ComplexMatrix&
ComplexMatrix::insert (const ColumnVector& a, int r, int c)
{
  int a_len = a.length ();
  if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
    {
      (*current_liboctave_error_handler) ("range error for insert");
      return *this;
    }

  for (int i = 0; i < a_len; i++)
    elem (r+i, c) = a.elem (i);

  return *this;
}

ComplexMatrix&
ComplexMatrix::insert (const DiagMatrix& a, int r, int c)
{
  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
    {
      (*current_liboctave_error_handler) ("range error for insert");
      return *this;
    }

  fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);

  for (int i = 0; i < a.length (); i++)
    elem (r+i, c+i) = a.elem (i, i);

  return *this;
}

ComplexMatrix&
ComplexMatrix::insert (const ComplexMatrix& a, int r, int c)
{
  Array2<Complex>::insert (a, r, c);
  return *this;
}

ComplexMatrix&
ComplexMatrix::insert (const ComplexRowVector& a, int r, int c)
{
  int a_len = a.length ();
  if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
    {
      (*current_liboctave_error_handler) ("range error for insert");
      return *this;
    }

  for (int i = 0; i < a_len; i++)
    elem (r, c+i) = a.elem (i);

  return *this;
}

ComplexMatrix&
ComplexMatrix::insert (const ComplexColumnVector& a, int r, int c)
{
  int a_len = a.length ();
  if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
    {
      (*current_liboctave_error_handler) ("range error for insert");
      return *this;
    }

  for (int i = 0; i < a_len; i++)
    elem (r+i, c) = a.elem (i);

  return *this;
}

ComplexMatrix&
ComplexMatrix::insert (const ComplexDiagMatrix& a, int r, int c)
{
  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
    {
      (*current_liboctave_error_handler) ("range error for insert");
      return *this;
    }

  fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);

  for (int i = 0; i < a.length (); i++)
    elem (r+i, c+i) = a.elem (i, i);

  return *this;
}

ComplexMatrix&
ComplexMatrix::fill (double val)
{
  int nr = rows ();
  int nc = cols ();
  if (nr > 0 && nc > 0)
    for (int j = 0; j < nc; j++)
      for (int i = 0; i < nr; i++)
	elem (i, j) = val;

  return *this;
}

ComplexMatrix&
ComplexMatrix::fill (const Complex& val)
{
  int nr = rows ();
  int nc = cols ();
  if (nr > 0 && nc > 0)
    for (int j = 0; j < nc; j++)
      for (int i = 0; i < nr; i++)
	elem (i, j) = val;

  return *this;
}

ComplexMatrix&
ComplexMatrix::fill (double val, int r1, int c1, int r2, int c2)
{
  int nr = rows ();
  int nc = cols ();
  if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
      || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
    {
      (*current_liboctave_error_handler) ("range error for fill");
      return *this;
    }

  if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }

  for (int j = c1; j <= c2; j++)
    for (int i = r1; i <= r2; i++)
      elem (i, j) = val;

  return *this;
}

ComplexMatrix&
ComplexMatrix::fill (const Complex& val, int r1, int c1, int r2, int c2)
{
  int nr = rows ();
  int nc = cols ();
  if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
      || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
    {
      (*current_liboctave_error_handler) ("range error for fill");
      return *this;
    }

  if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }

  for (int j = c1; j <= c2; j++)
    for (int i = r1; i <= r2; i++)
      elem (i, j) = val;

  return *this;
}

ComplexMatrix
ComplexMatrix::append (const Matrix& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nr != a.rows ())
    {
      (*current_liboctave_error_handler) ("row dimension mismatch for append");
      return *this;
    }

  int nc_insert = nc;
  ComplexMatrix retval (nr, nc + a.cols ());
  retval.insert (*this, 0, 0);
  retval.insert (a, 0, nc_insert);
  return retval;
}

ComplexMatrix
ComplexMatrix::append (const RowVector& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nr != 1)
    {
      (*current_liboctave_error_handler) ("row dimension mismatch for append");
      return *this;
    }

  int nc_insert = nc;
  ComplexMatrix retval (nr, nc + a.length ());
  retval.insert (*this, 0, 0);
  retval.insert (a, 0, nc_insert);
  return retval;
}

ComplexMatrix
ComplexMatrix::append (const ColumnVector& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nr != a.length ())
    {
      (*current_liboctave_error_handler) ("row dimension mismatch for append");
      return *this;
    }

  int nc_insert = nc;
  ComplexMatrix retval (nr, nc + 1);
  retval.insert (*this, 0, 0);
  retval.insert (a, 0, nc_insert);
  return retval;
}

ComplexMatrix
ComplexMatrix::append (const DiagMatrix& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nr != a.rows ())
    {
      (*current_liboctave_error_handler) ("row dimension mismatch for append");
      return *this;
    }

  int nc_insert = nc;
  ComplexMatrix retval (nr, nc + a.cols ());
  retval.insert (*this, 0, 0);
  retval.insert (a, 0, nc_insert);
  return retval;
}

ComplexMatrix
ComplexMatrix::append (const ComplexMatrix& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nr != a.rows ())
    {
      (*current_liboctave_error_handler) ("row dimension mismatch for append");
      return *this;
    }

  int nc_insert = nc;
  ComplexMatrix retval (nr, nc + a.cols ());
  retval.insert (*this, 0, 0);
  retval.insert (a, 0, nc_insert);
  return retval;
}

ComplexMatrix
ComplexMatrix::append (const ComplexRowVector& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nr != 1)
    {
      (*current_liboctave_error_handler) ("row dimension mismatch for append");
      return *this;
    }

  int nc_insert = nc;
  ComplexMatrix retval (nr, nc + a.length ());
  retval.insert (*this, 0, 0);
  retval.insert (a, 0, nc_insert);
  return retval;
}

ComplexMatrix
ComplexMatrix::append (const ComplexColumnVector& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nr != a.length ())
    {
      (*current_liboctave_error_handler) ("row dimension mismatch for append");
      return *this;
    }

  int nc_insert = nc;
  ComplexMatrix retval (nr, nc + 1);
  retval.insert (*this, 0, 0);
  retval.insert (a, 0, nc_insert);
  return retval;
}

ComplexMatrix
ComplexMatrix::append (const ComplexDiagMatrix& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nr != a.rows ())
    {
      (*current_liboctave_error_handler) ("row dimension mismatch for append");
      return *this;
    }

  int nc_insert = nc;
  ComplexMatrix retval (nr, nc + a.cols ());
  retval.insert (*this, 0, 0);
  retval.insert (a, 0, nc_insert);
  return retval;
}

ComplexMatrix
ComplexMatrix::stack (const Matrix& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nc != a.cols ())
    {
      (*current_liboctave_error_handler)
	("column dimension mismatch for stack");
      return *this;
    }

  int nr_insert = nr;
  ComplexMatrix retval (nr + a.rows (), nc);
  retval.insert (*this, 0, 0);
  retval.insert (a, nr_insert, 0);
  return retval;
}

ComplexMatrix
ComplexMatrix::stack (const RowVector& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nc != a.length ())
    {
      (*current_liboctave_error_handler)
	("column dimension mismatch for stack");
      return *this;
    }

  int nr_insert = nr;
  ComplexMatrix retval (nr + 1, nc);
  retval.insert (*this, 0, 0);
  retval.insert (a, nr_insert, 0);
  return retval;
}

ComplexMatrix
ComplexMatrix::stack (const ColumnVector& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nc != 1)
    {
      (*current_liboctave_error_handler)
	("column dimension mismatch for stack");
      return *this;
    }

  int nr_insert = nr;
  ComplexMatrix retval (nr + a.length (), nc);
  retval.insert (*this, 0, 0);
  retval.insert (a, nr_insert, 0);
  return retval;
}

ComplexMatrix
ComplexMatrix::stack (const DiagMatrix& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nc != a.cols ())
    {
      (*current_liboctave_error_handler)
	("column dimension mismatch for stack");
      return *this;
    }

  int nr_insert = nr;
  ComplexMatrix retval (nr + a.rows (), nc);
  retval.insert (*this, 0, 0);
  retval.insert (a, nr_insert, 0);
  return retval;
}

ComplexMatrix
ComplexMatrix::stack (const ComplexMatrix& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nc != a.cols ())
    {
      (*current_liboctave_error_handler)
	("column dimension mismatch for stack");
      return *this;
    }

  int nr_insert = nr;
  ComplexMatrix retval (nr + a.rows (), nc);
  retval.insert (*this, 0, 0);
  retval.insert (a, nr_insert, 0);
  return retval;
}

ComplexMatrix
ComplexMatrix::stack (const ComplexRowVector& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nc != a.length ())
    {
      (*current_liboctave_error_handler)
	("column dimension mismatch for stack");
      return *this;
    }

  int nr_insert = nr;
  ComplexMatrix retval (nr + 1, nc);
  retval.insert (*this, 0, 0);
  retval.insert (a, nr_insert, 0);
  return retval;
}

ComplexMatrix
ComplexMatrix::stack (const ComplexColumnVector& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nc != 1)
    {
      (*current_liboctave_error_handler)
	("column dimension mismatch for stack");
      return *this;
    }

  int nr_insert = nr;
  ComplexMatrix retval (nr + a.length (), nc);
  retval.insert (*this, 0, 0);
  retval.insert (a, nr_insert, 0);
  return retval;
}

ComplexMatrix
ComplexMatrix::stack (const ComplexDiagMatrix& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nc != a.cols ())
    {
      (*current_liboctave_error_handler)
	("column dimension mismatch for stack");
      return *this;
    }

  int nr_insert = nr;
  ComplexMatrix retval (nr + a.rows (), nc);
  retval.insert (*this, 0, 0);
  retval.insert (a, nr_insert, 0);
  return retval;
}

ComplexMatrix
ComplexMatrix::hermitian (void) const
{
  int nr = rows ();
  int nc = cols ();
  ComplexMatrix result;
  if (length () > 0)
    {
      result.resize (nc, nr);
      for (int j = 0; j < nc; j++)
	for (int i = 0; i < nr; i++)
	  result.elem (j, i) = conj (elem (i, j));
    }
  return result;
}

ComplexMatrix
ComplexMatrix::transpose (void) const
{
  int nr = rows ();
  int nc = cols ();
  ComplexMatrix result (nc, nr);
  if (length () > 0)
    {
      for (int j = 0; j < nc; j++)
	for (int i = 0; i < nr; i++)
	  result.elem (j, i) = elem (i, j);
    }
  return result;
}

ComplexMatrix
conj (const ComplexMatrix& a)
{
  int a_len = a.length ();
  ComplexMatrix retval;
  if (a_len > 0)
    retval = ComplexMatrix (conj_dup (a.data (), a_len), a.rows (),
			    a.cols ());
  return retval;
}

// resize is the destructive equivalent for this one

ComplexMatrix
ComplexMatrix::extract (int r1, int c1, int r2, int c2) const
{
  if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }

  int new_r = r2 - r1 + 1;
  int new_c = c2 - c1 + 1;

  ComplexMatrix result (new_r, new_c);

  for (int j = 0; j < new_c; j++)
    for (int i = 0; i < new_r; i++)
      result.elem (i, j) = elem (r1+i, c1+j);

  return result;
}

// extract row or column i.

ComplexRowVector
ComplexMatrix::row (int i) const
{
  int nc = cols ();
  if (i < 0 || i >= rows ())
    {
      (*current_liboctave_error_handler) ("invalid row selection");
      return ComplexRowVector ();
    }

  ComplexRowVector retval (nc);
  for (int j = 0; j < cols (); j++)
    retval.elem (j) = elem (i, j);

  return retval;
}

ComplexRowVector
ComplexMatrix::row (char *s) const
{
  if (! s)
    {
      (*current_liboctave_error_handler) ("invalid row selection");
      return ComplexRowVector ();
    }

  char c = *s;
  if (c == 'f' || c == 'F')
    return row (0);
  else if (c == 'l' || c == 'L')
    return row (rows () - 1);
  else
    {
      (*current_liboctave_error_handler) ("invalid row selection");
      return ComplexRowVector ();
    }
}

ComplexColumnVector
ComplexMatrix::column (int i) const
{
  int nr = rows ();
  if (i < 0 || i >= cols ())
    {
      (*current_liboctave_error_handler) ("invalid column selection");
      return ComplexColumnVector ();
    }

  ComplexColumnVector retval (nr);
  for (int j = 0; j < nr; j++)
    retval.elem (j) = elem (j, i);

  return retval;
}

ComplexColumnVector
ComplexMatrix::column (char *s) const
{
  if (! s)
    {
      (*current_liboctave_error_handler) ("invalid column selection");
      return ComplexColumnVector ();
    }

  char c = *s;
  if (c == 'f' || c == 'F')
    return column (0);
  else if (c == 'l' || c == 'L')
    return column (cols () - 1);
  else
    {
      (*current_liboctave_error_handler) ("invalid column selection");
      return ComplexColumnVector ();
    }
}

ComplexMatrix
ComplexMatrix::inverse (void) const
{
  int info;
  double rcond;
  return inverse (info, rcond);
}

ComplexMatrix
ComplexMatrix::inverse (int& info) const
{
  double rcond;
  return inverse (info, rcond);
}

ComplexMatrix
ComplexMatrix::inverse (int& info, double& rcond, int force) const
{
  ComplexMatrix retval;

  int nr = rows ();
  int nc = cols ();

  if (nr != nc)
    (*current_liboctave_error_handler) ("inverse requires square matrix");
  else
    {
      info = 0;

      Array<int> ipvt (nr);
      int *pipvt = ipvt.fortran_vec ();

      Array<Complex> z (nr);
      Complex *pz = z.fortran_vec ();

      retval = *this;
      Complex *tmp_data = retval.fortran_vec ();

      F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nc, pipvt, rcond, pz));

      if (f77_exception_encountered)
	(*current_liboctave_error_handler) ("unrecoverable error in zgeco");
      else
	{
	  volatile double rcond_plus_one = rcond + 1.0;

	  if (rcond_plus_one == 1.0)
	    info = -1;

	  if (info == -1 && ! force)
	    retval = *this;  // Restore contents.
	  else
	    {
	      Complex *dummy = 0;

	      F77_XFCN (zgedi, ZGEDI, (tmp_data, nr, nc, pipvt, dummy,
				       pz, 1));

	      if (f77_exception_encountered)
		(*current_liboctave_error_handler)
		  ("unrecoverable error in zgedi");
	    }
	}
    }

  return retval;
}

ComplexMatrix
ComplexMatrix::pseudo_inverse (double tol)
{
  ComplexMatrix retval;

  ComplexSVD result (*this);

  DiagMatrix S = result.singular_values ();
  ComplexMatrix U = result.left_singular_matrix ();
  ComplexMatrix V = result.right_singular_matrix ();

  ColumnVector sigma = S.diag ();

  int r = sigma.length () - 1;
  int nr = rows ();
  int nc = cols ();

  if (tol <= 0.0)
    {
      if (nr > nc)
	tol = nr * sigma.elem (0) * DBL_EPSILON;
      else
	tol = nc * sigma.elem (0) * DBL_EPSILON;
    }

  while (r >= 0 && sigma.elem (r) < tol)
    r--;

  if (r < 0)
    retval = ComplexMatrix (nc, nr, 0.0);
  else
    {
      ComplexMatrix Ur = U.extract (0, 0, nr-1, r);
      DiagMatrix D = DiagMatrix (sigma.extract (0, r)) . inverse ();
      ComplexMatrix Vr = V.extract (0, 0, nc-1, r);
      retval = Vr * D * Ur.hermitian ();
    }

  return retval;
}

ComplexMatrix
ComplexMatrix::fourier (void) const
{
  ComplexMatrix retval;

  int nr = rows ();
  int nc = cols ();

  int npts, nsamples;

  if (nr == 1 || nc == 1)
    {
      npts = nr > nc ? nr : nc;
      nsamples = 1;
    }
  else
    {
      npts = nr;
      nsamples = nc;
    }

  int nn = 4*npts+15;

  Array<Complex> wsave (nn);
  Complex *pwsave = wsave.fortran_vec ();

  retval = *this;
  Complex *tmp_data = retval.fortran_vec ();

  F77_FCN (cffti, CFFTI) (npts, pwsave);

  for (int j = 0; j < nsamples; j++)
    F77_FCN (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);

  return retval;
}

ComplexMatrix
ComplexMatrix::ifourier (void) const
{
  ComplexMatrix retval;

  int nr = rows ();
  int nc = cols ();

  int npts, nsamples;

  if (nr == 1 || nc == 1)
    {
      npts = nr > nc ? nr : nc;
      nsamples = 1;
    }
  else
    {
      npts = nr;
      nsamples = nc;
    }

  int nn = 4*npts+15;

  Array<Complex> wsave (nn);
  Complex *pwsave = wsave.fortran_vec ();

  retval = *this;
  Complex *tmp_data = retval.fortran_vec ();

  F77_FCN (cffti, CFFTI) (npts, pwsave);

  for (int j = 0; j < nsamples; j++)
    F77_FCN (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);

  for (int j = 0; j < npts*nsamples; j++)
    tmp_data[j] = tmp_data[j] / (double) npts;

  return retval;
}

ComplexMatrix
ComplexMatrix::fourier2d (void) const
{
  ComplexMatrix retval;

  int nr = rows ();
  int nc = cols ();

  int npts, nsamples;

  if (nr == 1 || nc == 1)
    {
      npts = nr > nc ? nr : nc;
      nsamples = 1;
    }
  else
    {
      npts = nr;
      nsamples = nc;
    }

  int nn = 4*npts+15;

  Array<Complex> wsave (nn);
  Complex *pwsave = wsave.fortran_vec ();

  retval = *this;
  Complex *tmp_data = retval.fortran_vec ();

  F77_FCN (cffti, CFFTI) (npts, pwsave);

  for (int j = 0; j < nsamples; j++)
    F77_FCN (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);

  npts = nc;
  nsamples = nr;
  nn = 4*npts+15;

  wsave.resize (nn);
  pwsave = wsave.fortran_vec ();

  Array<Complex> row (npts);
  Complex *prow = row.fortran_vec ();

  F77_FCN (cffti, CFFTI) (npts, pwsave);

  for (int j = 0; j < nsamples; j++)
    {
      for (int i = 0; i < npts; i++)
	prow[i] = tmp_data[i*nr + j];

      F77_FCN (cfftf, CFFTF) (npts, prow, pwsave);

      for (int i = 0; i < npts; i++)
	tmp_data[i*nr + j] = prow[i];
    }

  return retval;
}

ComplexMatrix
ComplexMatrix::ifourier2d (void) const
{
  ComplexMatrix retval;

  int nr = rows ();
  int nc = cols ();

  int npts, nsamples;

  if (nr == 1 || nc == 1)
    {
      npts = nr > nc ? nr : nc;
      nsamples = 1;
    }
  else
    {
      npts = nr;
      nsamples = nc;
    }

  int nn = 4*npts+15;

  Array<Complex> wsave (nn);
  Complex *pwsave = wsave.fortran_vec ();

  retval = *this;
  Complex *tmp_data = retval.fortran_vec ();

  F77_FCN (cffti, CFFTI) (npts, pwsave);

  for (int j = 0; j < nsamples; j++)
    F77_FCN (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);

  for (int j = 0; j < npts*nsamples; j++)
    tmp_data[j] = tmp_data[j] / (double) npts;

  npts = nc;
  nsamples = nr;
  nn = 4*npts+15;

  wsave.resize (nn);
  pwsave = wsave.fortran_vec ();

  Array<Complex> row (npts);
  Complex *prow = row.fortran_vec ();

  F77_FCN (cffti, CFFTI) (npts, pwsave);

  for (int j = 0; j < nsamples; j++)
    {
      for (int i = 0; i < npts; i++)
	prow[i] = tmp_data[i*nr + j];

      F77_FCN (cfftb, CFFTB) (npts, prow, pwsave);

      for (int i = 0; i < npts; i++)
	tmp_data[i*nr + j] = prow[i] / (double) npts;
    }

  return retval;
}

ComplexDET
ComplexMatrix::determinant (void) const
{
  int info;
  double rcond;
  return determinant (info, rcond);
}

ComplexDET
ComplexMatrix::determinant (int& info) const
{
  double rcond;
  return determinant (info, rcond);
}

ComplexDET
ComplexMatrix::determinant (int& info, double& rcond) const
{
  ComplexDET retval;

  int nr = rows ();
  int nc = cols ();

  if (nr == 0 || nc == 0)
    {
      Complex d[2];
      d[0] = 1.0;
      d[1] = 0.0;
      retval = ComplexDET (d);
    }
  else
    {
      info = 0;

      Array<int> ipvt (nr);
      int *pipvt = ipvt.fortran_vec ();

      Array<Complex> z (nr);
      Complex *pz = z.fortran_vec ();

      ComplexMatrix atmp = *this;
      Complex *tmp_data = atmp.fortran_vec ();

      F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz));

      if (f77_exception_encountered)
	(*current_liboctave_error_handler) ("unrecoverable error in zgeco");
      else
	{
	  volatile double rcond_plus_one = rcond + 1.0;

	  if (rcond_plus_one == 1.0)
	    {
	      info = -1;
	      retval = ComplexDET ();
	    }
	  else
	    {
	      Complex d[2];

	      F77_XFCN (zgedi, ZGEDI, (tmp_data, nr, nr, pipvt, d, pz, 10));

	      if (f77_exception_encountered)
		(*current_liboctave_error_handler)
		  ("unrecoverable error in dgedi");
	      else
		retval = ComplexDET (d);
	    }
	}
    }

  return retval;
}

ComplexMatrix
ComplexMatrix::solve (const Matrix& b) const
{
  int info;
  double rcond;
  return solve (b, info, rcond);
}

ComplexMatrix
ComplexMatrix::solve (const Matrix& b, int& info) const
{
  double rcond;
  return solve (b, info, rcond);
}

ComplexMatrix
ComplexMatrix::solve (const Matrix& b, int& info, double& rcond) const
{
  ComplexMatrix tmp (b);
  return solve (tmp, info, rcond);
}

ComplexMatrix
ComplexMatrix::solve (const ComplexMatrix& b) const
{
  int info;
  double rcond;
  return solve (b, info, rcond);
}

ComplexMatrix
ComplexMatrix::solve (const ComplexMatrix& b, int& info) const
{
  double rcond;
  return solve (b, info, rcond);
}
ComplexMatrix
ComplexMatrix::solve (const ComplexMatrix& b, int& info, double& rcond) const
{
  ComplexMatrix retval;

  int nr = rows ();
  int nc = cols ();

  if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
    (*current_liboctave_error_handler)
      ("matrix dimension mismatch in solution of linear equations");
  else
    {
      info = 0;

      Array<int> ipvt (nr);
      int *pipvt = ipvt.fortran_vec ();

      Array<Complex> z (nr);
      Complex *pz = z.fortran_vec ();

      ComplexMatrix atmp = *this;
      Complex *tmp_data = atmp.fortran_vec ();

      F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz));

      if (f77_exception_encountered)
	(*current_liboctave_error_handler) ("unrecoverable error in zgeco");
      else
	{
	  volatile double rcond_plus_one = rcond + 1.0;

	  if (rcond_plus_one == 1.0)
	    {
	      info = -2;
	    }
	  else
	    {
	      retval = b;
	      Complex *result = retval.fortran_vec ();

	      int b_nc = b.cols ();

	      for (volatile int j = 0; j < b_nc; j++)
		{
		  F77_XFCN (zgesl, ZGESL, (tmp_data, nr, nr, pipvt,
					   &result[nr*j], 0));

		  if (f77_exception_encountered)
		    {
		      (*current_liboctave_error_handler)
			("unrecoverable error in dgesl");

		      break;
		    }
		}
	    }
	}
    }

  return retval;
}

ComplexColumnVector
ComplexMatrix::solve (const ComplexColumnVector& b) const
{
  int info;
  double rcond;
  return solve (b, info, rcond);
}

ComplexColumnVector
ComplexMatrix::solve (const ComplexColumnVector& b, int& info) const
{
  double rcond;
  return solve (b, info, rcond);
}

ComplexColumnVector
ComplexMatrix::solve (const ComplexColumnVector& b, int& info,
		      double& rcond) const
{
  ComplexColumnVector retval;

  int nr = rows ();
  int nc = cols ();

  if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
    (*current_liboctave_error_handler)
      ("matrix dimension mismatch in solution of linear equations");
  else
    {
      info = 0;

      Array<int> ipvt (nr);
      int *pipvt = ipvt.fortran_vec ();

      Array<Complex> z (nr);
      Complex *pz = z.fortran_vec ();

      ComplexMatrix atmp = *this;
      Complex *tmp_data = atmp.fortran_vec ();

      F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz));

      if (f77_exception_encountered)
	(*current_liboctave_error_handler)
	  ("unrecoverable error in dgeco");
      else
	{
	  volatile double rcond_plus_one = rcond + 1.0;

	  if (rcond_plus_one == 1.0)
	    {
	      info = -2;
	    }
	  else
	    {
	      retval = b;
	      Complex *result = retval.fortran_vec ();

	      F77_XFCN (zgesl, ZGESL, (tmp_data, nr, nr, pipvt, result, 0));

	      if (f77_exception_encountered)
		(*current_liboctave_error_handler)
		  ("unrecoverable error in dgesl");
	    }
	}
    }

  return retval;
}

ComplexMatrix
ComplexMatrix::lssolve (const ComplexMatrix& b) const
{
  int info;
  int rank;
  return lssolve (b, info, rank);
}

ComplexMatrix
ComplexMatrix::lssolve (const ComplexMatrix& b, int& info) const
{
  int rank;
  return lssolve (b, info, rank);
}

ComplexMatrix
ComplexMatrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const
{
  ComplexMatrix retval;

  int nrhs = b.cols ();

  int m = rows ();
  int n = cols ();

  if (m == 0 || n == 0 || m != b.rows ())
    (*current_liboctave_error_handler)
      ("matrix dimension mismatch solution of linear equations");
  else
    {
      ComplexMatrix atmp = *this;
      Complex *tmp_data = atmp.fortran_vec ();

      int nrr = m > n ? m : n;
      ComplexMatrix result (nrr, nrhs);

      for (int j = 0; j < nrhs; j++)
	for (int i = 0; i < m; i++)
	  result.elem (i, j) = b.elem (i, j);

      Complex *presult = result.fortran_vec ();

      int len_s = m < n ? m : n;
      Array<double> s (len_s);
      double *ps = s.fortran_vec ();

      double rcond = -1.0;

      int lwork;
      if (m < n)
	lwork = 2*m + (nrhs > n ? nrhs : n);
      else
	lwork = 2*n + (nrhs > m ? nrhs : m);

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

      int lrwork = (5 * (m < n ? m : n)) - 4;
      lrwork = lrwork > 1 ? lrwork : 1;
      Array<double> rwork (lrwork);
      double *prwork = rwork.fortran_vec ();

      F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult,
				 nrr, ps, rcond, rank, pwork, lwork,
				 prwork, info));

      if (f77_exception_encountered)
	(*current_liboctave_error_handler) ("unrecoverable error in zgelss");
      else
	{
	  retval.resize (n, nrhs);
	  for (int j = 0; j < nrhs; j++)
	    for (int i = 0; i < n; i++)
	      retval.elem (i, j) = result.elem (i, j);
	}
    }

  return retval;
}

ComplexColumnVector
ComplexMatrix::lssolve (const ComplexColumnVector& b) const
{
  int info;
  int rank;
  return lssolve (b, info, rank);
}

ComplexColumnVector
ComplexMatrix::lssolve (const ComplexColumnVector& b, int& info) const
{
  int rank;
  return lssolve (b, info, rank);
}

ComplexColumnVector
ComplexMatrix::lssolve (const ComplexColumnVector& b, int& info,
			int& rank) const
{
  ComplexColumnVector retval;

  int nrhs = 1;

  int m = rows ();
  int n = cols ();

  if (m == 0 || n == 0 || m != b.length ())
    (*current_liboctave_error_handler)
      ("matrix dimension mismatch solution of least squares problem");
  else
    {
      ComplexMatrix atmp = *this;
      Complex *tmp_data = atmp.fortran_vec ();

      int nrr = m > n ? m : n;
      ComplexColumnVector result (nrr);

      for (int i = 0; i < m; i++)
	result.elem (i) = b.elem (i);

      Complex *presult = result.fortran_vec ();

      int len_s = m < n ? m : n;
      Array<double> s (len_s);
      double *ps = s.fortran_vec ();

      double rcond = -1.0;

      int lwork;
      if (m < n)
	lwork = 2*m + (nrhs > n ? nrhs : n);
      else
	lwork = 2*n + (nrhs > m ? nrhs : m);

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

      int lrwork = (5 * (m < n ? m : n)) - 4;
      lrwork = lrwork > 1 ? lrwork : 1;
      Array<double> rwork (lrwork);
      double *prwork = rwork.fortran_vec ();

      F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult,
				 nrr, ps, rcond, rank, pwork, lwork,
				 prwork, info));

      if (f77_exception_encountered)
	(*current_liboctave_error_handler) ("unrecoverable error in zgelss");
      else
	{
	  retval.resize (n);
	  for (int i = 0; i < n; i++)
	    retval.elem (i) = result.elem (i);
	}
    }

  return retval;
}

// Constants for matrix exponential calculation.

static double padec [] =
{
  5.0000000000000000e-1,
  1.1666666666666667e-1,
  1.6666666666666667e-2,
  1.6025641025641026e-3,
  1.0683760683760684e-4,
  4.8562548562548563e-6,
  1.3875013875013875e-7,
  1.9270852604185938e-9,
};

ComplexMatrix
ComplexMatrix::expm (void) const
{
  ComplexMatrix retval;

  ComplexMatrix m = *this;

  int nc = columns ();

  // trace shift value
  Complex trshift = 0.0;

  // Preconditioning step 1: trace normalization.

  for (int i = 0; i < nc; i++)
    trshift += m.elem (i, i);

  trshift /= nc;

  if(trshift.real() < 0)
    trshift = 0.0 + trshift.imag(); 
    // don't want stable eigenvalues to become unstable

  for (int i = 0; i < nc; i++)
    m.elem (i, i) -= trshift;

  // Preconditioning step 2: eigenvalue balancing.

  ComplexAEPBALANCE mbal (m, "B");
  m = mbal.balanced_matrix ();
  ComplexMatrix d = mbal.balancing_matrix ();

  // Preconditioning step 3: scaling.

  ColumnVector work (nc);
  double inf_norm
    = F77_FCN (zlange, ZLANGE) ("I", nc, nc, m.fortran_vec (), nc,
				work.fortran_vec ());

  int sqpow = (int) (inf_norm > 0.0
		     ? (1.0 + log (inf_norm) / log (2.0))
		     : 0.0);

  // Check whether we need to square at all.

  if (sqpow < 0)
    sqpow = 0;

  if (sqpow > 0)
    {
      double scale_factor = 1.0;
      for (int i = 0; i < sqpow; i++)
	scale_factor *= 2.0;

      m = m / scale_factor;
    }

  // npp, dpp: pade' approx polynomial matrices.

  ComplexMatrix npp (nc, nc, 0.0);
  ComplexMatrix dpp = npp;

  // Now powers a^8 ... a^1.

  int minus_one_j = -1;
  for (int j = 7; j >= 0; j--)
    {
      npp = m * npp + m * padec[j];
      dpp = m * dpp + m * (minus_one_j * padec[j]);
      minus_one_j *= -1;
    }

  // Zero power.

  dpp = -dpp;
  for (int j = 0; j < nc; j++)
    {
      npp.elem (j, j) += 1.0;
      dpp.elem (j, j) += 1.0;
    }

  // Compute pade approximation = inverse (dpp) * npp.

  retval = dpp.solve (npp);
	
  // Reverse preconditioning step 3: repeated squaring.

  while (sqpow)
    {
      retval = retval * retval;
      sqpow--;
    }

  // Reverse preconditioning step 2: inverse balancing.
  // XXX FIXME XXX -- should probably do this with Lapack calls
  // instead of a complete matrix inversion.

  retval = retval.transpose ();
  d = d.transpose ();
  retval = retval * d;
  retval = d.solve (retval);
  retval = retval.transpose ();

  // Reverse preconditioning step 1: fix trace normalization.

  return retval * exp (trshift);
}

// column vector by row vector -> matrix operations

ComplexMatrix
operator * (const ColumnVector& v, const ComplexRowVector& a)
{
  ComplexColumnVector tmp (v);
  return tmp * a;
}

ComplexMatrix
operator * (const ComplexColumnVector& a, const RowVector& b)
{
  ComplexRowVector tmp (b);
  return a * tmp;
}

ComplexMatrix
operator * (const ComplexColumnVector& v, const ComplexRowVector& a)
{
  ComplexMatrix retval;

  int len = v.length ();
  int a_len = a.length ();

  if (len != a_len)
    gripe_nonconformant ("operator *", len, 1, 1, a_len);
  else
    {
      if (len != 0)
	{
	  retval.resize (len, a_len);
	  Complex *c = retval.fortran_vec ();

	  F77_XFCN (zgemm, ZGEMM, ("N", "N", len, a_len, 1, 1.0,
				   v.data (), len, a.data (), 1, 0.0,
				   c, len, 1L, 1L)); 

	  if (f77_exception_encountered)
	    (*current_liboctave_error_handler)
	      ("unrecoverable error in zgemm");
	}
    }

  return retval;
}

// diagonal matrix by scalar -> matrix operations

ComplexMatrix
operator + (const DiagMatrix& a, const Complex& s)
{
  ComplexMatrix tmp (a.rows (), a.cols (), s);
  return a + tmp;
}

ComplexMatrix
operator - (const DiagMatrix& a, const Complex& s)
{
  ComplexMatrix tmp (a.rows (), a.cols (), -s);
  return a + tmp;
}

ComplexMatrix
operator + (const ComplexDiagMatrix& a, double s)
{
  ComplexMatrix tmp (a.rows (), a.cols (), s);
  return a + tmp;
}

ComplexMatrix
operator - (const ComplexDiagMatrix& a, double s)
{
  ComplexMatrix tmp (a.rows (), a.cols (), -s);
  return a + tmp;
}

ComplexMatrix
operator + (const ComplexDiagMatrix& a, const Complex& s)
{
  ComplexMatrix tmp (a.rows (), a.cols (), s);
  return a + tmp;
}

ComplexMatrix
operator - (const ComplexDiagMatrix& a, const Complex& s)
{
  ComplexMatrix tmp (a.rows (), a.cols (), -s);
  return a + tmp;
}

// scalar by diagonal matrix -> matrix operations

ComplexMatrix
operator + (const Complex& s, const DiagMatrix& a)
{
  ComplexMatrix tmp (a.rows (), a.cols (), s);
  return tmp + a;
}

ComplexMatrix
operator - (const Complex& s, const DiagMatrix& a)
{
  ComplexMatrix tmp (a.rows (), a.cols (), s);
  return tmp - a;
}

ComplexMatrix
operator + (double s, const ComplexDiagMatrix& a)
{
  ComplexMatrix tmp (a.rows (), a.cols (), s);
  return tmp + a;
}

ComplexMatrix
operator - (double s, const ComplexDiagMatrix& a)
{
  ComplexMatrix tmp (a.rows (), a.cols (), s);
  return tmp - a;
}

ComplexMatrix
operator + (const Complex& s, const ComplexDiagMatrix& a)
{
  ComplexMatrix tmp (a.rows (), a.cols (), s);
  return tmp + a;
}

ComplexMatrix
operator - (const Complex& s, const ComplexDiagMatrix& a)
{
  ComplexMatrix tmp (a.rows (), a.cols (), s);
  return tmp - a;
}

// matrix by diagonal matrix -> matrix operations

ComplexMatrix&
ComplexMatrix::operator += (const DiagMatrix& a)
{
  int nr = rows ();
  int nc = cols ();

  int a_nr = rows ();
  int a_nc = cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
      return *this;
    }

  for (int i = 0; i < a.length (); i++)
    elem (i, i) += a.elem (i, i);

  return *this;
}

ComplexMatrix&
ComplexMatrix::operator -= (const DiagMatrix& a)
{
  int nr = rows ();
  int nc = cols ();

  int a_nr = rows ();
  int a_nc = cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
      return *this;
    }

  for (int i = 0; i < a.length (); i++)
    elem (i, i) -= a.elem (i, i);

  return *this;
}

ComplexMatrix&
ComplexMatrix::operator += (const ComplexDiagMatrix& a)
{
  int nr = rows ();
  int nc = cols ();

  int a_nr = rows ();
  int a_nc = cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
      return *this;
    }

  for (int i = 0; i < a.length (); i++)
    elem (i, i) += a.elem (i, i);

  return *this;
}

ComplexMatrix&
ComplexMatrix::operator -= (const ComplexDiagMatrix& a)
{
  int nr = rows ();
  int nc = cols ();

  int a_nr = rows ();
  int a_nc = cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
      return *this;
    }

  for (int i = 0; i < a.length (); i++)
    elem (i, i) -= a.elem (i, i);

  return *this;
}

ComplexMatrix
operator + (const Matrix& m, const ComplexDiagMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  ComplexMatrix result (m);
  for (int i = 0; i < a.length (); i++)
    result.elem (i, i) += a.elem (i, i);

  return result;
}

ComplexMatrix
operator - (const Matrix& m, const ComplexDiagMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  ComplexMatrix result (m);
  for (int i = 0; i < a.length (); i++)
    result.elem (i, i) -= a.elem (i, i);

  return result;
}

ComplexMatrix
operator * (const Matrix& m, const ComplexDiagMatrix& a)
{
  ComplexMatrix retval;

  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nc != a_nr)
    gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  else
    {
      if (nr == 0 || nc == 0 || a_nc == 0)
	retval.resize (nr, a_nc, 0.0);
      else
	{
	  retval.resize (nr, a_nc);
	  Complex *c = retval.fortran_vec ();

	  Complex *ctmp = 0;

	  for (int j = 0; j < a.length (); j++)
	    {
	      int idx = j * nr;
	      ctmp = c + idx;
	      if (a.elem (j, j) == 1.0)
		{
		  for (int i = 0; i < nr; i++)
		    ctmp[i] = m.elem (i, j);
		}
	      else if (a.elem (j, j) == 0.0)
		{
		  for (int i = 0; i < nr; i++)
		    ctmp[i] = 0.0;
		}
	      else
		{
		  for (int i = 0; i < nr; i++)
		    ctmp[i] = a.elem (j, j) * m.elem (i, j);
		}
	    }

	  if (a_nr < a_nc)
	    {
	      for (int i = nr * nc; i < nr * a_nc; i++)
		ctmp[i] = 0.0;
	    }
	}
    }

  return retval;
}

// diagonal matrix by matrix -> matrix operations

ComplexMatrix
operator + (const DiagMatrix& m, const ComplexMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  ComplexMatrix result (a);
  for (int i = 0; i < m.length (); i++)
    result.elem (i, i) += m.elem (i, i);

  return result;
}

ComplexMatrix
operator - (const DiagMatrix& m, const ComplexMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  ComplexMatrix result (-a);
  for (int i = 0; i < m.length (); i++)
    result.elem (i, i) += m.elem (i, i);

  return result;
}

ComplexMatrix
operator * (const DiagMatrix& m, const ComplexMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nc != a_nr)
    {
      gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0 || a_nc == 0)
    return ComplexMatrix (nr, nc, 0.0);

  ComplexMatrix c (nr, a_nc);

  for (int i = 0; i < m.length (); i++)
    {
      if (m.elem (i, i) == 1.0)
	{
	  for (int j = 0; j < a_nc; j++)
	    c.elem (i, j) = a.elem (i, j);
	}
      else if (m.elem (i, i) == 0.0)
	{
	  for (int j = 0; j < a_nc; j++)
	    c.elem (i, j) = 0.0;
	}
      else
	{
	  for (int j = 0; j < a_nc; j++)
	    c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
	}
    }

  if (nr > nc)
    {
      for (int j = 0; j < a_nc; j++)
	for (int i = a_nr; i < nr; i++)
	  c.elem (i, j) = 0.0;
    }

  return c;
}

ComplexMatrix
operator + (const ComplexDiagMatrix& m, const Matrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  ComplexMatrix result (a);
  for (int i = 0; i < m.length (); i++)
    result.elem (i, i) += m.elem (i, i);

  return result;
}

ComplexMatrix
operator - (const ComplexDiagMatrix& m, const Matrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  ComplexMatrix result (-a);
  for (int i = 0; i < m.length (); i++)
    result.elem (i, i) += m.elem (i, i);

  return result;
}

ComplexMatrix
operator * (const ComplexDiagMatrix& m, const Matrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nc != a_nr)
    {
      gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0 || a_nc == 0)
    return ComplexMatrix (nr, a_nc, 0.0);

  ComplexMatrix c (nr, a_nc);

  for (int i = 0; i < m.length (); i++)
    {
      if (m.elem (i, i) == 1.0)
	{
	  for (int j = 0; j < a_nc; j++)
	    c.elem (i, j) = a.elem (i, j);
	}
      else if (m.elem (i, i) == 0.0)
	{
	  for (int j = 0; j < a_nc; j++)
	    c.elem (i, j) = 0.0;
	}
      else
	{
	  for (int j = 0; j < a_nc; j++)
	    c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
	}
    }

  if (nr > nc)
    {
      for (int j = 0; j < a_nc; j++)
	for (int i = a_nr; i < nr; i++)
	  c.elem (i, j) = 0.0;
    }

  return c;
}

ComplexMatrix
operator + (const ComplexDiagMatrix& m, const ComplexMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  ComplexMatrix result (a);
  for (int i = 0; i < m.length (); i++)
    result.elem (i, i) += m.elem (i, i);

  return result;
}

ComplexMatrix
operator - (const ComplexDiagMatrix& m, const ComplexMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  ComplexMatrix result (-a);
  for (int i = 0; i < m.length (); i++)
    result.elem (i, i) += m.elem (i, i);

  return result;
}

ComplexMatrix
operator * (const ComplexDiagMatrix& m, const ComplexMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nc != a_nr)
    {
      gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0 || a_nc == 0)
    return ComplexMatrix (nr, a_nc, 0.0);

  ComplexMatrix c (nr, a_nc);

  for (int i = 0; i < m.length (); i++)
    {
      if (m.elem (i, i) == 1.0)
	{
	  for (int j = 0; j < a_nc; j++)
	    c.elem (i, j) = a.elem (i, j);
	}
      else if (m.elem (i, i) == 0.0)
	{
	  for (int j = 0; j < a_nc; j++)
	    c.elem (i, j) = 0.0;
	}
      else
	{
	  for (int j = 0; j < a_nc; j++)
	    c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
	}
    }

  if (nr > nc)
    {
      for (int j = 0; j < a_nc; j++)
	for (int i = a_nr; i < nr; i++)
	  c.elem (i, j) = 0.0;
    }

  return c;
}

// matrix by matrix -> matrix operations

ComplexMatrix&
ComplexMatrix::operator += (const Matrix& a)
{
  int nr = rows ();
  int nc = cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
      return *this;
    }

  if (nr == 0 || nc == 0)
    return *this;

  Complex *d = fortran_vec (); // Ensures only one reference to my privates!

  add2 (d, a.data (), length ());
  return *this;
}

ComplexMatrix&
ComplexMatrix::operator -= (const Matrix& a)
{
  int nr = rows ();
  int nc = cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
      return *this;
    }

  if (nr == 0 || nc == 0)
    return *this;

  Complex *d = fortran_vec (); // Ensures only one reference to my privates!

  subtract2 (d, a.data (), length ());
  return *this;
}

ComplexMatrix&
ComplexMatrix::operator += (const ComplexMatrix& a)
{
  int nr = rows ();
  int nc = cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
      return *this;
    }

  if (nr == 0 || nc == 0)
    return *this;

  Complex *d = fortran_vec (); // Ensures only one reference to my privates!

  add2 (d, a.data (), length ());
  return *this;
}

ComplexMatrix&
ComplexMatrix::operator -= (const ComplexMatrix& a)
{
  int nr = rows ();
  int nc = cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
      return *this;
    }

  if (nr == 0 || nc == 0)
    return *this;

  Complex *d = fortran_vec (); // Ensures only one reference to my privates!

  subtract2 (d, a.data (), length ());
  return *this;
}

// unary operations

Matrix
ComplexMatrix::operator ! (void) const
{
  return Matrix (not (data (), length ()), rows (), cols ());
}

// matrix by scalar -> matrix operations

ComplexMatrix
operator + (const Matrix& a, const Complex& s)
{
  return ComplexMatrix (add (a.data (), a.length (), s),
			a.rows (), a.cols ());
}

ComplexMatrix
operator - (const Matrix& a, const Complex& s)
{
  return ComplexMatrix (subtract (a.data (), a.length (), s),
			a.rows (), a.cols ());
}

ComplexMatrix
operator * (const Matrix& a, const Complex& s)
{
  return ComplexMatrix (multiply (a.data (), a.length (), s),
			a.rows (), a.cols ());
}

ComplexMatrix
operator / (const Matrix& a, const Complex& s)
{
  return ComplexMatrix (divide (a.data (), a.length (), s),
			a.rows (), a.cols ());
}

ComplexMatrix
operator + (const ComplexMatrix& a, double s)
{
  return ComplexMatrix (add (a.data (), a.length (), s),
			a.rows (), a.cols ());
}

ComplexMatrix
operator - (const ComplexMatrix& a, double s)
{
  return ComplexMatrix (subtract (a.data (), a.length (), s),
			a.rows (), a.cols ());
}

ComplexMatrix
operator * (const ComplexMatrix& a, double s)
{
  return ComplexMatrix (multiply (a.data (), a.length (), s),
			a.rows (), a.cols ());
}

ComplexMatrix
operator / (const ComplexMatrix& a, double s)
{
  return ComplexMatrix (divide (a.data (), a.length (), s),
			a.rows (), a.cols ());
}

// scalar by matrix -> matrix operations

ComplexMatrix
operator + (double s, const ComplexMatrix& a)
{
  return ComplexMatrix (add (a.data (), a.length (), s), a.rows (),
			a.cols ());
}

ComplexMatrix
operator - (double s, const ComplexMatrix& a)
{
  return ComplexMatrix (subtract (s, a.data (), a.length ()),
			a.rows (), a.cols ());
}

ComplexMatrix
operator * (double s, const ComplexMatrix& a)
{
  return ComplexMatrix (multiply (a.data (), a.length (), s),
			a.rows (), a.cols ());
}

ComplexMatrix
operator / (double s, const ComplexMatrix& a)
{
  return ComplexMatrix (divide (s, a.data (), a.length ()),
			a.rows (), a.cols ());
}

ComplexMatrix
operator + (const Complex& s, const Matrix& a)
{
  return ComplexMatrix (add (s, a.data (), a.length ()),
			a.rows (), a.cols ());
}

ComplexMatrix
operator - (const Complex& s, const Matrix& a)
{
  return ComplexMatrix (subtract (s, a.data (), a.length ()),
			a.rows (), a.cols ());
}

ComplexMatrix
operator * (const Complex& s, const Matrix& a)
{
  return ComplexMatrix (multiply (a.data (), a.length (), s),
			a.rows (), a.cols ());
}

ComplexMatrix
operator / (const Complex& s, const Matrix& a)
{
  return ComplexMatrix (divide (s, a.data (), a.length ()),
			a.rows (), a.cols ());
}

// matrix by diagonal matrix -> matrix operations

ComplexMatrix
operator + (const ComplexMatrix& m, const DiagMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  ComplexMatrix result (m);
  for (int i = 0; i < a.length (); i++)
    result.elem (i, i) += a.elem (i, i);

  return result;
}

ComplexMatrix
operator - (const ComplexMatrix& m, const DiagMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  ComplexMatrix result (m);
  for (int i = 0; i < a.length (); i++)
    result.elem (i, i) -= a.elem (i, i);

  return result;
}

ComplexMatrix
operator * (const ComplexMatrix& m, const DiagMatrix& a)
{
  ComplexMatrix retval;

  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nc != a_nr)
    gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  else
    {
      if (nr == 0 || nc == 0 || a_nc == 0)
	retval.resize (nr, nc, 0.0);
      else
	{
	  retval.resize (nr, a_nc);
	  Complex *c = retval.fortran_vec ();
	  Complex *ctmp = 0;

	  for (int j = 0; j < a.length (); j++)
	    {
	      int idx = j * nr;
	      ctmp = c + idx;
	      if (a.elem (j, j) == 1.0)
		{
		  for (int i = 0; i < nr; i++)
		    ctmp[i] = m.elem (i, j);
		}
	      else if (a.elem (j, j) == 0.0)
		{
		  for (int i = 0; i < nr; i++)
		    ctmp[i] = 0.0;
		}
	      else
		{
		  for (int i = 0; i < nr; i++)
		    ctmp[i] = a.elem (j, j) * m.elem (i, j);
		}
	    }

	  if (a.rows () < a_nc)
	    {
	      for (int i = nr * nc; i < nr * a_nc; i++)
		ctmp[i] = 0.0;
	    }
	}
    }

  return retval;
}

ComplexMatrix
operator + (const ComplexMatrix& m, const ComplexDiagMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  ComplexMatrix result (m);
  for (int i = 0; i < a.length (); i++)
    result.elem (i, i) += a.elem (i, i);

  return result;
}

ComplexMatrix
operator - (const ComplexMatrix& m, const ComplexDiagMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  ComplexMatrix result (m);
  for (int i = 0; i < a.length (); i++)
    result.elem (i, i) -= a.elem (i, i);

  return result;
}

ComplexMatrix
operator * (const ComplexMatrix& m, const ComplexDiagMatrix& a)
{
  ComplexMatrix retval;

  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nc != a_nr)
    gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  else
    {
      if (nr == 0 || nc == 0 || a_nc == 0)
	retval.resize (nr, nc, 0.0);
      else
	{
	  retval.resize (nr, nc);
	  Complex *c = retval.fortran_vec ();
	  Complex *ctmp = 0;

	  for (int j = 0; j < a.length (); j++)
	    {
	      int idx = j * nr;
	      ctmp = c + idx;
	      if (a.elem (j, j) == 1.0)
		{
		  for (int i = 0; i < nr; i++)
		    ctmp[i] = m.elem (i, j);
		}
	      else if (a.elem (j, j) == 0.0)
		{
		  for (int i = 0; i < nr; i++)
		    ctmp[i] = 0.0;
		}
	      else
		{
		  for (int i = 0; i < nr; i++)
		    ctmp[i] = a.elem (j, j) * m.elem (i, j);
		}
	    }

	  if (a.rows () < a_nc)
	    {
	      for (int i = nr * nc; i < nr * a_nc; i++)
		ctmp[i] = 0.0;
	    }
	}
    }

  return retval;
}

// matrix by matrix -> matrix operations

ComplexMatrix
operator + (const ComplexMatrix& m, const Matrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  return ComplexMatrix (add (m.data (), a.data (), m.length ()), nr, nc);
}

ComplexMatrix
operator - (const ComplexMatrix& m, const Matrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  return ComplexMatrix (subtract (m.data (), a.data (), m.length ()), nr, nc);
}

ComplexMatrix
operator + (const Matrix& m, const ComplexMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  return ComplexMatrix (add (m.data (), a.data (), m.length ()), nr, nc);
}

ComplexMatrix
operator - (const Matrix& m, const ComplexMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  return ComplexMatrix (subtract (m.data (), a.data (), m.length ()), nr, nc);
}

ComplexMatrix
operator * (const ComplexMatrix& m, const Matrix& a)
{
  ComplexMatrix tmp (a);
  return m * tmp;
}

ComplexMatrix
operator * (const Matrix& m, const ComplexMatrix& a)
{
  ComplexMatrix tmp (m);
  return tmp * a;
}

ComplexMatrix
operator * (const ComplexMatrix& m, const ComplexMatrix& a)
{
  ComplexMatrix retval;

  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nc != a_nr)
    gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  else
    {
      if (nr == 0 || nc == 0 || a_nc == 0)
	retval.resize (nr, nc, 0.0);
      else
	{
	  int ld  = nr;
	  int lda = a.rows ();

	  retval.resize (nr, a_nc);
	  Complex *c = retval.fortran_vec ();

	  F77_XFCN (zgemm, ZGEMM, ("N", "N", nr, a_nc, nc, 1.0,
				   m.data (), ld, a.data (), lda, 0.0,
				   c, nr, 1L, 1L));

	  if (f77_exception_encountered)
	    (*current_liboctave_error_handler)
	      ("unrecoverable error in zgemm");
	}
    }

  return retval;
}

ComplexMatrix
product (const ComplexMatrix& m, const Matrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("product", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  return ComplexMatrix (multiply (m.data (), a.data (), m.length ()), nr, nc);
}

ComplexMatrix
quotient (const ComplexMatrix& m, const Matrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("quotient", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  return ComplexMatrix (divide (m.data (), a.data (), m.length ()), nr, nc);
}

ComplexMatrix
product (const Matrix& m, const ComplexMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("product", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  return ComplexMatrix (multiply (m.data (), a.data (), m.length ()), nr, nc);
}

ComplexMatrix
quotient (const Matrix& m, const ComplexMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("quotient", nr, nc, a_nr, a_nc);
      return ComplexMatrix ();
    }

  if (nr == 0 || nc == 0)
    return ComplexMatrix (nr, nc);

  return ComplexMatrix (divide (m.data (), a.data (), m.length ()), nr, nc);
}

// other operations

ComplexMatrix
ComplexMatrix::map (c_c_Mapper f) const
{
  ComplexMatrix b (*this);
  return b.apply (f);
}

Matrix
ComplexMatrix::map (d_c_Mapper f) const
{
  const Complex *d = data ();

  Matrix retval (rows (), columns ());

  double *r = retval.fortran_vec ();

  for (int i = 0; i < length (); i++)
    r[i] = f (d[i]);

  return retval;
}

ComplexMatrix&
ComplexMatrix::apply (c_c_Mapper f)
{
  Complex *d = fortran_vec (); // Ensures only one reference to my privates!

  for (int i = 0; i < length (); i++)
    d[i] = f (d[i]);

  return *this;
}

bool
ComplexMatrix::any_element_is_inf_or_nan (void) const
{
  int nr = rows ();
  int nc = cols ();

  for (int j = 0; j < nc; j++)
    for (int i = 0; i < nr; i++)
      {
	Complex val = elem (i, j);
	if (xisinf (val) || xisnan (val))
	  return true;
      }

  return false;
}

// Return true if no elements have imaginary components.

bool
ComplexMatrix::all_elements_are_real (void) const
{
  int nr = rows ();
  int nc = cols ();

  for (int j = 0; j < nc; j++)
    for (int i = 0; i < nr; i++)
      if (imag (elem (i, j)) != 0.0)
	return false;

  return true;
}

// Return nonzero if any element of CM has a non-integer real or
// imaginary part.  Also extract the largest and smallest (real or
// imaginary) values and return them in MAX_VAL and MIN_VAL. 

bool
ComplexMatrix::all_integers (double& max_val, double& min_val) const
{
  int nr = rows ();
  int nc = cols ();

  if (nr > 0 && nc > 0)
    {
      Complex val = elem (0, 0);

      double r_val = real (val);
      double i_val = imag (val);

      max_val = r_val;
      min_val = r_val;

      if (i_val > max_val)
	max_val = i_val;

      if (i_val < max_val)
	min_val = i_val;
    }
  else
    return false;

  for (int j = 0; j < nc; j++)
    for (int i = 0; i < nr; i++)
      {
	Complex val = elem (i, j);

	double r_val = real (val);
	double i_val = imag (val);

	if (r_val > max_val)
	  max_val = r_val;

	if (i_val > max_val)
	  max_val = i_val;

	if (r_val < min_val)
	  min_val = r_val;

	if (i_val < min_val)
	  min_val = i_val;

	if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val)
	  return false;
      }

  return true;
}

bool
ComplexMatrix::too_large_for_float (void) const
{
  int nr = rows ();
  int nc = cols ();

  for (int j = 0; j < nc; j++)
    for (int i = 0; i < nr; i++)
      {
	Complex val = elem (i, j);

	double r_val = real (val);
	double i_val = imag (val);

	if (r_val > FLT_MAX
	    || i_val > FLT_MAX
	    || r_val < FLT_MIN
	    || i_val < FLT_MIN)
	  return true;
      }

  return false;
}

Matrix
ComplexMatrix::all (void) const
{
  int nr = rows ();
  int nc = cols ();
  Matrix retval;
  if (nr > 0 && nc > 0)
    {
      if (nr == 1)
	{
	  retval.resize (1, 1);
	  retval.elem (0, 0) = 1.0;
	  for (int j = 0; j < nc; j++)
	    {
	      if (elem (0, j) == 0.0)
		{
		  retval.elem (0, 0) = 0.0;
		  break;
		}
	    }
	}
      else if (nc == 1)
	{
	  retval.resize (1, 1);
	  retval.elem (0, 0) = 1.0;
	  for (int i = 0; i < nr; i++)
	    {
	      if (elem (i, 0) == 0.0)
		{
		  retval.elem (0, 0) = 0.0;
		  break;
		}
	    }
	}
      else
	{
	  retval.resize (1, nc);
	  for (int j = 0; j < nc; j++)
	    {
	      retval.elem (0, j) = 1.0;
	      for (int i = 0; i < nr; i++)
		{
		  if (elem (i, j) == 0.0)
		    {
		      retval.elem (0, j) = 0.0;
		      break;
		    }
		}
	    }
	}
    }
  return retval;
}

Matrix
ComplexMatrix::any (void) const
{
  int nr = rows ();
  int nc = cols ();
  Matrix retval;
  if (nr > 0 && nc > 0)
    {
      if (nr == 1)
	{
	  retval.resize (1, 1);
	  retval.elem (0, 0) = 0.0;
	  for (int j = 0; j < nc; j++)
	    {
	      if (elem (0, j) != 0.0)
		{
		  retval.elem (0, 0) = 1.0;
		  break;
		}
	    }
	}
      else if (nc == 1)
	{
	  retval.resize (1, 1);
	  retval.elem (0, 0) = 0.0;
	  for (int i = 0; i < nr; i++)
	    {
	      if (elem (i, 0) != 0.0)
		{
		  retval.elem (0, 0) = 1.0;
		  break;
		}
	    }
	}
      else
	{
	  retval.resize (1, nc);
	  for (int j = 0; j < nc; j++)
	    {
	      retval.elem (0, j) = 0.0;
	      for (int i = 0; i < nr; i++)
		{
		  if (elem (i, j) != 0.0)
		    {
		      retval.elem (0, j) = 1.0;
		      break;
		    }
		}
	    }
	}
    }
  return retval;
}

ComplexMatrix
ComplexMatrix::cumprod (void) const
{
  int nr = rows ();
  int nc = cols ();
  ComplexMatrix retval;
  if (nr > 0 && nc > 0)
    {
      if (nr == 1)
	{
	  retval.resize (1, nc);
	  Complex prod = elem (0, 0);
	  for (int j = 0; j < nc; j++)
	    {
	      retval.elem (0, j) = prod;
	      if (j < nc - 1)
		prod *= elem (0, j+1);
	    }
	}
      else if (nc == 1)
	{
	  retval.resize (nr, 1);
	  Complex prod = elem (0, 0);
	  for (int i = 0; i < nr; i++)
	    {
	      retval.elem (i, 0) = prod;
	      if (i < nr - 1)
		prod *= elem (i+1, 0);
	    }
	}
      else
	{
	  retval.resize (nr, nc);
	  for (int j = 0; j < nc; j++)
	    {
	      Complex prod = elem (0, j);
	      for (int i = 0; i < nr; i++)
		{
		  retval.elem (i, j) = prod;
		  if (i < nr - 1)
		    prod *= elem (i+1, j);
		}
	    }
	}
    }
  return retval;
}

ComplexMatrix
ComplexMatrix::cumsum (void) const
{
  int nr = rows ();
  int nc = cols ();
  ComplexMatrix retval;
  if (nr > 0 && nc > 0)
    {
      if (nr == 1)
	{
	  retval.resize (1, nc);
	  Complex sum = elem (0, 0);
	  for (int j = 0; j < nc; j++)
	    {
	      retval.elem (0, j) = sum;
	      if (j < nc - 1)
		sum += elem (0, j+1);
	    }
	}
      else if (nc == 1)
	{
	  retval.resize (nr, 1);
	  Complex sum = elem (0, 0);
	  for (int i = 0; i < nr; i++)
	    {
	      retval.elem (i, 0) = sum;
	      if (i < nr - 1)
		sum += elem (i+1, 0);
	    }
	}
      else
	{
	  retval.resize (nr, nc);
	  for (int j = 0; j < nc; j++)
	    {
	      Complex sum = elem (0, j);
	      for (int i = 0; i < nr; i++)
		{
		  retval.elem (i, j) = sum;
		  if (i < nr - 1)
		    sum += elem (i+1, j);
		}
	    }
	}
    }
  return retval;
}

ComplexMatrix
ComplexMatrix::prod (void) const
{
  int nr = rows ();
  int nc = cols ();
  ComplexMatrix retval;
  if (nr > 0 && nc > 0)
    {
      if (nr == 1)
	{
	  retval.resize (1, 1);
	  retval.elem (0, 0) = 1.0;
	  for (int j = 0; j < nc; j++)
	    retval.elem (0, 0) *= elem (0, j);
	}
      else if (nc == 1)
	{
	  retval.resize (1, 1);
	  retval.elem (0, 0) = 1.0;
	  for (int i = 0; i < nr; i++)
	    retval.elem (0, 0) *= elem (i, 0);
	}
      else
	{
	  retval.resize (1, nc);
	  for (int j = 0; j < nc; j++)
	    {
	      retval.elem (0, j) = 1.0;
	      for (int i = 0; i < nr; i++)
		retval.elem (0, j) *= elem (i, j);
	    }
	}
    }
  return retval;
}

ComplexMatrix
ComplexMatrix::sum (void) const
{
  int nr = rows ();
  int nc = cols ();
  ComplexMatrix retval;
  if (nr > 0 && nc > 0)
    {
      if (nr == 1)
	{
	  retval.resize (1, 1);
	  retval.elem (0, 0) = 0.0;
	  for (int j = 0; j < nc; j++)
	    retval.elem (0, 0) += elem (0, j);
	}
      else if (nc == 1)
	{
	  retval.resize (1, 1);
	  retval.elem (0, 0) = 0.0;
	  for (int i = 0; i < nr; i++)
	    retval.elem (0, 0) += elem (i, 0);
	}
      else
	{
	  retval.resize (1, nc);
	  for (int j = 0; j < nc; j++)
	    {
	      retval.elem (0, j) = 0.0;
	      for (int i = 0; i < nr; i++)
		retval.elem (0, j) += elem (i, j);
	    }
	}
    }
  return retval;
}

ComplexMatrix
ComplexMatrix::sumsq (void) const
{
  int nr = rows ();
  int nc = cols ();
  ComplexMatrix retval;
  if (nr > 0 && nc > 0)
    {
      if (nr == 1)
	{
	  retval.resize (1, 1);
	  retval.elem (0, 0) = 0.0;
	  for (int j = 0; j < nc; j++)
	    {
	      Complex d = elem (0, j);
	      retval.elem (0, 0) += d * d;
	    }
	}
      else if (nc == 1)
	{
	  retval.resize (1, 1);
	  retval.elem (0, 0) = 0.0;
	  for (int i = 0; i < nr; i++)
	    {
	      Complex d = elem (i, 0);
	      retval.elem (0, 0) += d * d;
	    }
	}
      else
	{
	  retval.resize (1, nc);
	  for (int j = 0; j < nc; j++)
	    {
	      retval.elem (0, j) = 0.0;
	      for (int i = 0; i < nr; i++)
		{
		  Complex d = elem (i, j);
		  retval.elem (0, j) += d * d;
		}
	    }
	}
    }
  return retval;
}

ComplexColumnVector
ComplexMatrix::diag (void) const
{
  return diag (0);
}

ComplexColumnVector
ComplexMatrix::diag (int k) const
{
  int nnr = rows ();
  int nnc = cols ();
  if (k > 0)
    nnc -= k;
  else if (k < 0)
    nnr += k;

  ComplexColumnVector d;

  if (nnr > 0 && nnc > 0)
    {
      int ndiag = (nnr < nnc) ? nnr : nnc;

      d.resize (ndiag);

      if (k > 0)
	{
	  for (int i = 0; i < ndiag; i++)
	    d.elem (i) = elem (i, i+k);
	}
      else if ( k < 0)
	{
	  for (int i = 0; i < ndiag; i++)
	    d.elem (i) = elem (i-k, i);
	}
      else
	{
	  for (int i = 0; i < ndiag; i++)
	    d.elem (i) = elem (i, i);
	}
    }
  else
    cerr << "diag: requested diagonal out of range\n";

  return d;
}

bool
ComplexMatrix::row_is_real_only (int i) const
{
  bool retval = true;

  int nc = columns ();

  for (int j = 0; j < nc; j++)
    {
      if (imag (elem (i, j)) != 0.0)
	{
	  retval = false;
	  break;
	}
    }

  return retval;	      
}

bool
ComplexMatrix::column_is_real_only (int j) const
{
  bool retval = true;

  int nr = rows ();

  for (int i = 0; i < nr; i++)
    {
      if (imag (elem (i, j)) != 0.0)
	{
	  retval = false;
	  break;
	}
    }

  return retval;	      
}

ComplexColumnVector
ComplexMatrix::row_min (void) const
{
  Array<int> index;
  return row_min (index);
}

ComplexColumnVector
ComplexMatrix::row_min (Array<int>& index) const
{
  ComplexColumnVector result;

  int nr = rows ();
  int nc = cols ();

  if (nr > 0 && nc > 0)
    {
      result.resize (nr);
      index.resize (nr);

      for (int i = 0; i < nr; i++)
        {
	  int idx = 0;

	  Complex tmp_min = elem (i, idx);

	  bool real_only = row_is_real_only (i);

	  double abs_min = real_only ? real (tmp_min) : abs (tmp_min);

	  if (xisnan (tmp_min))
	    idx = -1;
	  else
	    {
	      for (int j = 1; j < nc; j++)
		{
		  Complex tmp = elem (i, j);

		  double abs_tmp = real_only ? real (tmp) : abs (tmp);

		  if (xisnan (tmp))
		    {
		      idx = -1;
		      break;
		    }
		  else if (abs_tmp < abs_min)
		    {
		      idx = j;
		      tmp_min = tmp;
		      abs_min = abs_tmp;
		    }
		}

	      result.elem (i) = (idx < 0) ? Complex_NaN_result : tmp_min;
	      index.elem (i) = idx;
	    }
        }
    }

  return result;
}

ComplexColumnVector
ComplexMatrix::row_max (void) const
{
  Array<int> index;
  return row_max (index);
}

ComplexColumnVector
ComplexMatrix::row_max (Array<int>& index) const
{
  ComplexColumnVector result;

  int nr = rows ();
  int nc = cols ();

  if (nr > 0 && nc > 0)
    {
      result.resize (nr);
      index.resize (nr);

      for (int i = 0; i < nr; i++)
        {
	  int idx = 0;

	  Complex tmp_max = elem (i, idx);

	  bool real_only = row_is_real_only (i);

	  double abs_max = real_only ? real (tmp_max) : abs (tmp_max);

	  if (xisnan (tmp_max))
	    idx = -1;
	  else
	    {
	      for (int j = 1; j < nc; j++)
		{
		  Complex tmp = elem (i, j);

		  double abs_tmp = real_only ? real (tmp) : abs (tmp);

		  if (xisnan (tmp))
		    {
		      idx = -1;
		      break;
		    }
		  else if (abs_tmp > abs_max)
		    {
		      idx = j;
		      tmp_max = tmp;
		      abs_max = abs_tmp;
		    }
		}

	      result.elem (i) = (idx < 0) ? Complex_NaN_result : tmp_max;
	      index.elem (i) = idx;
	    }
        }
    }

  return result;
}

ComplexRowVector
ComplexMatrix::column_min (void) const
{
  Array<int> index;
  return column_min (index);
}

ComplexRowVector
ComplexMatrix::column_min (Array<int>& index) const
{
  ComplexRowVector result;

  int nr = rows ();
  int nc = cols ();

  if (nr > 0 && nc > 0)
    {
      result.resize (nc);
      index.resize (nc);

      for (int j = 0; j < nc; j++)
        {
	  int idx = 0;

	  Complex tmp_min = elem (idx, j);

	  bool real_only = column_is_real_only (j);

	  double abs_min = real_only ? real (tmp_min) : abs (tmp_min);

	  if (xisnan (tmp_min))
	    idx = -1;
	  else
	    {
	      for (int i = 1; i < nr; i++)
		{
		  Complex tmp = elem (i, j);

		  double abs_tmp = real_only ? real (tmp) : abs (tmp);

		  if (xisnan (tmp))
		    {
		      idx = -1;
		      break;
		    }
		  else if (abs_tmp < abs_min)
		    {
		      idx = i;
		      tmp_min = tmp;
		      abs_min = abs_tmp;
		    }
		}

	      result.elem (j) = (idx < 0) ? Complex_NaN_result : tmp_min;
	      index.elem (j) = idx;
	    }
        }
    }

  return result;
}

ComplexRowVector
ComplexMatrix::column_max (void) const
{
  Array<int> index;
  return column_max (index);
}

ComplexRowVector
ComplexMatrix::column_max (Array<int>& index) const
{
  ComplexRowVector result;

  int nr = rows ();
  int nc = cols ();

  if (nr > 0 && nc > 0)
    {
      result.resize (nc);
      index.resize (nc);

      for (int j = 0; j < nc; j++)
        {
	  int idx = 0;

	  Complex tmp_max = elem (idx, j);

	  bool real_only = column_is_real_only (j);

	  double abs_max = real_only ? real (tmp_max) : abs (tmp_max);

	  if (xisnan (tmp_max))
	    idx = -1;
	  else
	    {
	      for (int i = 1; i < nr; i++)
		{
		  Complex tmp = elem (i, j);

		  double abs_tmp = real_only ? real (tmp) : abs (tmp);

		  if (xisnan (tmp))
		    {
		      idx = -1;
		      break;
		    }
		  else if (abs_tmp > abs_max)
		    {
		      idx = i;
		      tmp_max = tmp;
		      abs_max = abs_tmp;
		    }
		}

	      result.elem (j) = (idx < 0) ? Complex_NaN_result : tmp_max;
	      index.elem (j) = idx;
	    }
        }
    }

  return result;
}

// i/o

ostream&
operator << (ostream& os, const ComplexMatrix& a)
{
//  int field_width = os.precision () + 7;
  for (int i = 0; i < a.rows (); i++)
    {
      for (int j = 0; j < a.cols (); j++)
	os << " " /* setw (field_width) */ << a.elem (i, j);
      os << "\n";
    }
  return os;
}

istream&
operator >> (istream& is, ComplexMatrix& a)
{
  int nr = a.rows ();
  int nc = a.cols ();

  if (nr < 1 || nc < 1)
    is.clear (ios::badbit);
  else
    {
      Complex tmp;
      for (int i = 0; i < nr; i++)
	for (int j = 0; j < nc; j++)
	  {
	    is >> tmp;
	    if (is)
	      a.elem (i, j) = tmp;
	    else
	      goto done;
	  }
    }

done:

  return is;
}

ComplexMatrix
Givens (const Complex& x, const Complex& y)
{
  double cc;
  Complex cs, temp_r;
 
  F77_FCN (zlartg, ZLARTG) (x, y, cc, cs, temp_r);

  ComplexMatrix g (2, 2);

  g.elem (0, 0) = cc;
  g.elem (1, 1) = cc;
  g.elem (0, 1) = cs;
  g.elem (1, 0) = -conj (cs);

  return g;
}

ComplexMatrix
Sylvester (const ComplexMatrix& a, const ComplexMatrix& b,
	   const ComplexMatrix& c)
{
  ComplexMatrix retval;

  // XXX FIXME XXX -- need to check that a, b, and c are all the same
  // size.

  // Compute Schur decompositions

  ComplexSCHUR as (a, "U");
  ComplexSCHUR bs (b, "U");
  
  // Transform c to new coordinates.

  ComplexMatrix ua = as.unitary_matrix ();
  ComplexMatrix sch_a = as.schur_matrix ();

  ComplexMatrix ub = bs.unitary_matrix ();
  ComplexMatrix sch_b = bs.schur_matrix ();
  
  ComplexMatrix cx = ua.hermitian () * c * ub;

  // Solve the sylvester equation, back-transform, and return the
  // solution.

  int a_nr = a.rows ();
  int b_nr = b.rows ();

  double scale;
  int info;

  Complex *pa = sch_a.fortran_vec ();
  Complex *pb = sch_b.fortran_vec ();
  Complex *px = cx.fortran_vec ();
  
  F77_XFCN (ztrsyl, ZTRSYL, ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb,
			     b_nr, px, a_nr, scale,
			     info, 1L, 1L));

  if (f77_exception_encountered)
    (*current_liboctave_error_handler) ("unrecoverable error in ztrsyl");
  else
    {
      // XXX FIXME XXX -- check info?

      retval = -ua * cx * ub.hermitian ();
    }

  return retval;
}

/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/
----------
X-Sun-Data-Type: default
X-Sun-Data-Description: default
X-Sun-Data-Name: dMatrix.cc
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 3293

// Matrix manipulations.
/*

Copyright (C) 1996 John W. Eaton

This file is part of Octave.

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

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

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

*/

#if defined (__GNUG__)
#pragma implementation
#endif

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <cfloat>

#include <iostream.h>

#include "byte-swap.h"
#include "dbleAEPBAL.h"
#include "dbleDET.h"
#include "dbleSCHUR.h"
#include "dbleSVD.h"
#include "f77-fcn.h"
#include "lo-error.h"
#include "lo-ieee.h"
#include "lo-mappers.h"
#include "lo-utils.h"
#include "mx-base.h"
#include "mx-inlines.cc"
#include "oct-cmplx.h"

// Fortran functions we call.

extern "C"
{
  int F77_FCN (dgemm, DGEMM) (const char*, const char*, const int&,
			      const int&, const int&, const double&,
			      const double*, const int&,
			      const double*, const int&,
			      const double&, double*, const int&,
			      long, long);

  int F77_FCN (dgeco, DGECO) (double*, const int&, const int&, int*,
			      double&, double*);

  int F77_FCN (dgesl, DGESL) (const double*, const int&, const int&,
			      const int*, double*, const int&);

  int F77_FCN (dgedi, DGEDI) (double*, const int&, const int&,
			      const int*, double*, double*,
			      const int&);

  int F77_FCN (dgelss, DGELSS) (const int&, const int&, const int&,
				double*, const int&, double*,
				const int&, double*, double&, int&,
				double*, const int&, int&);

  // Note that the original complex fft routines were not written for
  // double complex arguments.  They have been modified by adding an
  // implicit double precision (a-h,o-z) statement at the beginning of
  // each subroutine.

  int F77_FCN (cffti, CFFTI) (const int&, Complex*);

  int F77_FCN (cfftf, CFFTF) (const int&, Complex*, Complex*);

  int F77_FCN (cfftb, CFFTB) (const int&, Complex*, Complex*);

  int F77_FCN (dlartg, DLARTG) (const double&, const double&, double&,
				double&, double&);

  int F77_FCN (dtrsyl, DTRSYL) (const char*, const char*, const int&,
				const int&, const int&, const double*,
				const int&, const double*, const int&,
				const double*, const int&, double&,
				int&, long, long);

  double F77_FCN (dlange, DLANGE) (const char*, const int&,
				   const int&, const double*,
				   const int&, double*); 

  int F77_FCN (qzhes, QZHES) (const int&, const int&, double*,
			      double*, const long&, double*);
 
  int F77_FCN (qzit, QZIT) (const int&, const int&, double*, double*,
			    const double&, const long&, double*,
			    int&);
 
  int F77_FCN (qzval, QZVAL) (const int&, const int&, double*,
			      double*, double*, double*, double*,
			      const long&, double*);
}

// Matrix class.

Matrix::Matrix (const RowVector& rv)
  : MArray2<double> (1, rv.length (), 0.0)
{
  for (int i = 0; i < rv.length (); i++)
    elem (0, i) = rv.elem (i);
}

Matrix::Matrix (const ColumnVector& cv)
  : MArray2<double> (cv.length (), 1, 0.0)
{
  for (int i = 0; i < cv.length (); i++)
    elem (i, 0) = cv.elem (i);
}

Matrix::Matrix (const DiagMatrix& a)
  : MArray2<double> (a.rows (), a.cols (), 0.0)
{
  for (int i = 0; i < a.length (); i++)
    elem (i, i) = a.elem (i, i);
}

// XXX FIXME XXX -- could we use a templated mixed-type copy function
// here?

Matrix::Matrix (const charMatrix& a)
  : MArray2<double> (a.rows (), a.cols ())
{
  for (int i = 0; i < a.rows (); i++)
    for (int j = 0; j < a.cols (); j++)
      elem (i, j) = a.elem (i, j);
}

bool
Matrix::operator == (const Matrix& a) const
{
  if (rows () != a.rows () || cols () != a.cols ())
    return false;

  return equal (data (), a.data (), length ());
}

bool
Matrix::operator != (const Matrix& a) const
{
  return !(*this == a);
}

Matrix&
Matrix::insert (const Matrix& a, int r, int c)
{
  Array2<double>::insert (a, r, c);
  return *this;
}

Matrix&
Matrix::insert (const RowVector& a, int r, int c)
{
  int a_len = a.length ();
  if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
    {
      (*current_liboctave_error_handler) ("range error for insert");
      return *this;
    }

  for (int i = 0; i < a_len; i++)
    elem (r, c+i) = a.elem (i);

  return *this;
}

Matrix&
Matrix::insert (const ColumnVector& a, int r, int c)
{
  int a_len = a.length ();
  if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
    {
      (*current_liboctave_error_handler) ("range error for insert");
      return *this;
    }

  for (int i = 0; i < a_len; i++)
    elem (r+i, c) = a.elem (i);

  return *this;
}

Matrix&
Matrix::insert (const DiagMatrix& a, int r, int c)
{
  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
    {
      (*current_liboctave_error_handler) ("range error for insert");
      return *this;
    }

  fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);

  for (int i = 0; i < a.length (); i++)
    elem (r+i, c+i) = a.elem (i, i);

  return *this;
}

Matrix&
Matrix::fill (double val)
{
  int nr = rows ();
  int nc = cols ();
  if (nr > 0 && nc > 0)
    for (int j = 0; j < nc; j++)
      for (int i = 0; i < nr; i++)
	elem (i, j) = val;

  return *this;
}

Matrix&
Matrix::fill (double val, int r1, int c1, int r2, int c2)
{
  int nr = rows ();
  int nc = cols ();
  if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
      || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
    {
      (*current_liboctave_error_handler) ("range error for fill");
      return *this;
    }

  if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }

  for (int j = c1; j <= c2; j++)
    for (int i = r1; i <= r2; i++)
      elem (i, j) = val;

  return *this;
}

Matrix
Matrix::append (const Matrix& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nr != a.rows ())
    {
      (*current_liboctave_error_handler) ("row dimension mismatch for append");
      return Matrix ();
    }

  int nc_insert = nc;
  Matrix retval (nr, nc + a.cols ());
  retval.insert (*this, 0, 0);
  retval.insert (a, 0, nc_insert);
  return retval;
}

Matrix
Matrix::append (const RowVector& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nr != 1)
    {
      (*current_liboctave_error_handler) ("row dimension mismatch for append");
      return Matrix ();
    }

  int nc_insert = nc;
  Matrix retval (nr, nc + a.length ());
  retval.insert (*this, 0, 0);
  retval.insert (a, 0, nc_insert);
  return retval;
}

Matrix
Matrix::append (const ColumnVector& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nr != a.length ())
    {
      (*current_liboctave_error_handler) ("row dimension mismatch for append");
      return Matrix ();
    }

  int nc_insert = nc;
  Matrix retval (nr, nc + 1);
  retval.insert (*this, 0, 0);
  retval.insert (a, 0, nc_insert);
  return retval;
}

Matrix
Matrix::append (const DiagMatrix& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nr != a.rows ())
    {
      (*current_liboctave_error_handler) ("row dimension mismatch for append");
      return *this;
    }

  int nc_insert = nc;
  Matrix retval (nr, nc + a.cols ());
  retval.insert (*this, 0, 0);
  retval.insert (a, 0, nc_insert);
  return retval;
}

Matrix
Matrix::stack (const Matrix& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nc != a.cols ())
    {
      (*current_liboctave_error_handler)
	("column dimension mismatch for stack");
      return Matrix ();
    }

  int nr_insert = nr;
  Matrix retval (nr + a.rows (), nc);
  retval.insert (*this, 0, 0);
  retval.insert (a, nr_insert, 0);
  return retval;
}

Matrix
Matrix::stack (const RowVector& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nc != a.length ())
    {
      (*current_liboctave_error_handler)
	("column dimension mismatch for stack");
      return Matrix ();
    }

  int nr_insert = nr;
  Matrix retval (nr + 1, nc);
  retval.insert (*this, 0, 0);
  retval.insert (a, nr_insert, 0);
  return retval;
}

Matrix
Matrix::stack (const ColumnVector& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nc != 1)
    {
      (*current_liboctave_error_handler)
	("column dimension mismatch for stack");
      return Matrix ();
    }

  int nr_insert = nr;
  Matrix retval (nr + a.length (), nc);
  retval.insert (*this, 0, 0);
  retval.insert (a, nr_insert, 0);
  return retval;
}

Matrix
Matrix::stack (const DiagMatrix& a) const
{
  int nr = rows ();
  int nc = cols ();
  if (nc != a.cols ())
    {
      (*current_liboctave_error_handler)
	("column dimension mismatch for stack");
      return Matrix ();
    }

  int nr_insert = nr;
  Matrix retval (nr + a.rows (), nc);
  retval.insert (*this, 0, 0);
  retval.insert (a, nr_insert, 0);
  return retval;
}

Matrix
Matrix::transpose (void) const
{
  int nr = rows ();
  int nc = cols ();
  Matrix result (nc, nr);
  if (length () > 0)
    {
      for (int j = 0; j < nc; j++)
	for (int i = 0; i < nr; i++)
	  result.elem (j, i) = elem (i, j);
    }
  return result;
}

Matrix
real (const ComplexMatrix& a)
{
  int a_len = a.length ();
  Matrix retval;
  if (a_len > 0)
    retval = Matrix (real_dup (a.data (), a_len), a.rows (), a.cols ());
  return retval;
}

Matrix
imag (const ComplexMatrix& a)
{
  int a_len = a.length ();
  Matrix retval;
  if (a_len > 0)
    retval = Matrix (imag_dup (a.data (), a_len), a.rows (), a.cols ());
  return retval;
}

Matrix
Matrix::extract (int r1, int c1, int r2, int c2) const
{
  if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }

  int new_r = r2 - r1 + 1;
  int new_c = c2 - c1 + 1;

  Matrix result (new_r, new_c);

  for (int j = 0; j < new_c; j++)
    for (int i = 0; i < new_r; i++)
      result.elem (i, j) = elem (r1+i, c1+j);

  return result;
}

// extract row or column i.

RowVector
Matrix::row (int i) const
{
  int nc = cols ();
  if (i < 0 || i >= rows ())
    {
      (*current_liboctave_error_handler) ("invalid row selection");
      return RowVector ();
    }

  RowVector retval (nc);
  for (int j = 0; j < nc; j++)
    retval.elem (j) = elem (i, j);

  return retval;
}

RowVector
Matrix::row (char *s) const
{
  if (! s)
    {
      (*current_liboctave_error_handler) ("invalid row selection");
      return RowVector ();
    }

  char c = *s;
  if (c == 'f' || c == 'F')
    return row (0);
  else if (c == 'l' || c == 'L')
    return row (rows () - 1);
  else
    {
      (*current_liboctave_error_handler) ("invalid row selection");
      return RowVector ();
    }
}

ColumnVector
Matrix::column (int i) const
{
  int nr = rows ();
  if (i < 0 || i >= cols ())
    {
      (*current_liboctave_error_handler) ("invalid column selection");
      return ColumnVector ();
    }

  ColumnVector retval (nr);
  for (int j = 0; j < nr; j++)
    retval.elem (j) = elem (j, i);

  return retval;
}

ColumnVector
Matrix::column (char *s) const
{
  if (! s)
    {
      (*current_liboctave_error_handler) ("invalid column selection");
      return ColumnVector ();
    }

  char c = *s;
  if (c == 'f' || c == 'F')
    return column (0);
  else if (c == 'l' || c == 'L')
    return column (cols () - 1);
  else
    {
      (*current_liboctave_error_handler) ("invalid column selection");
      return ColumnVector ();
    }
}

Matrix
Matrix::inverse (void) const
{
  int info;
  double rcond;
  return inverse (info, rcond);
}

Matrix
Matrix::inverse (int& info) const
{
  double rcond;
  return inverse (info, rcond);
}

Matrix
Matrix::inverse (int& info, double& rcond, int force) const
{
  Matrix retval;

  int nr = rows ();
  int nc = cols ();

  if (nr != nc || nr == 0 || nc == 0)
    (*current_liboctave_error_handler) ("inverse requires square matrix");
  else
    {
      info = 0;

      Array<int> ipvt (nr);
      int *pipvt = ipvt.fortran_vec ();

      Array<double> z (nr);
      double *pz = z.fortran_vec ();

      retval = *this;
      double *tmp_data = retval.fortran_vec ();

      F77_XFCN (dgeco, DGECO, (tmp_data, nr, nc, pipvt, rcond, pz));

      if (f77_exception_encountered)
	(*current_liboctave_error_handler) ("unrecoverable error in dgeco");
      else
	{
	  volatile double rcond_plus_one = rcond + 1.0;

	  if (rcond_plus_one == 1.0)
	    info = -1;

	  if (info == -1 && ! force)
	    retval = *this; // Restore matrix contents.
	  else
	    {
	      double *dummy = 0;

	      F77_XFCN (dgedi, DGEDI, (tmp_data, nr, nc, pipvt, dummy,
				       pz, 1));

	      if (f77_exception_encountered)
		(*current_liboctave_error_handler)
		  ("unrecoverable error in dgedi");
	    }
	}
    }

  return retval;
}

Matrix
Matrix::pseudo_inverse (double tol)
{
  SVD result (*this);

  DiagMatrix S = result.singular_values ();
  Matrix U = result.left_singular_matrix ();
  Matrix V = result.right_singular_matrix ();

  ColumnVector sigma = S.diag ();

  int r = sigma.length () - 1;
  int nr = rows ();
  int nc = cols ();

  if (tol <= 0.0)
    {
      if (nr > nc)
	tol = nr * sigma.elem (0) * DBL_EPSILON;
      else
	tol = nc * sigma.elem (0) * DBL_EPSILON;
    }

  while (r >= 0 && sigma.elem (r) < tol)
    r--;

  if (r < 0)
    return Matrix (nc, nr, 0.0);
  else
    {
      Matrix Ur = U.extract (0, 0, nr-1, r);
      DiagMatrix D = DiagMatrix (sigma.extract (0, r)) . inverse ();
      Matrix Vr = V.extract (0, 0, nc-1, r);
      return Vr * D * Ur.transpose ();
    }
}

ComplexMatrix
Matrix::fourier (void) const
{
  ComplexMatrix retval;

  int nr = rows ();
  int nc = cols ();

  int npts, nsamples;

  if (nr == 1 || nc == 1)
    {
      npts = nr > nc ? nr : nc;
      nsamples = 1;
    }
  else
    {
      npts = nr;
      nsamples = nc;
    }

  int nn = 4*npts+15;

  Array<Complex> wsave (nn);
  Complex *pwsave = wsave.fortran_vec ();

  retval = *this;
  Complex *tmp_data = retval.fortran_vec ();

  F77_FCN (cffti, CFFTI) (npts, pwsave);

  for (int j = 0; j < nsamples; j++)
    F77_FCN (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);

  return retval;
}

ComplexMatrix
Matrix::ifourier (void) const
{
  ComplexMatrix retval;

  int nr = rows ();
  int nc = cols ();

  int npts, nsamples;

  if (nr == 1 || nc == 1)
    {
      npts = nr > nc ? nr : nc;
      nsamples = 1;
    }
  else
    {
      npts = nr;
      nsamples = nc;
    }

  int nn = 4*npts+15;

  Array<Complex> wsave (nn);
  Complex *pwsave = wsave.fortran_vec ();

  retval = *this;
  Complex *tmp_data = retval.fortran_vec ();

  F77_FCN (cffti, CFFTI) (npts, pwsave);

  for (int j = 0; j < nsamples; j++)
    F77_FCN (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);

  for (int j = 0; j < npts*nsamples; j++)
    tmp_data[j] = tmp_data[j] / (double) npts;

  return retval;
}

ComplexMatrix
Matrix::fourier2d (void) const
{
  ComplexMatrix retval;

  int nr = rows ();
  int nc = cols ();

  int npts, nsamples;

  if (nr == 1 || nc == 1)
    {
      npts = nr > nc ? nr : nc;
      nsamples = 1;
    }
  else
    {
      npts = nr;
      nsamples = nc;
    }

  int nn = 4*npts+15;

  Array<Complex> wsave (nn);
  Complex *pwsave = wsave.fortran_vec ();

  retval = *this;
  Complex *tmp_data = retval.fortran_vec ();

  F77_FCN (cffti, CFFTI) (npts, pwsave);

  for (int j = 0; j < nsamples; j++)
    F77_FCN (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);

  npts = nc;
  nsamples = nr;
  nn = 4*npts+15;

  wsave.resize (nn);
  pwsave = wsave.fortran_vec ();

  Array<Complex> row (npts);
  Complex *prow = row.fortran_vec ();

  F77_FCN (cffti, CFFTI) (npts, pwsave);

  for (int j = 0; j < nsamples; j++)
    {
      for (int i = 0; i < npts; i++)
	prow[i] = tmp_data[i*nr + j];

      F77_FCN (cfftf, CFFTF) (npts, prow, pwsave);

      for (int i = 0; i < npts; i++)
	tmp_data[i*nr + j] = prow[i];
    }

  return retval;
}

ComplexMatrix
Matrix::ifourier2d (void) const
{
  ComplexMatrix retval;

  int nr = rows ();
  int nc = cols ();

  int npts, nsamples;

  if (nr == 1 || nc == 1)
    {
      npts = nr > nc ? nr : nc;
      nsamples = 1;
    }
  else
    {
      npts = nr;
      nsamples = nc;
    }

  int nn = 4*npts+15;

  Array<Complex> wsave (nn);
  Complex *pwsave = wsave.fortran_vec ();

  retval = *this;
  Complex *tmp_data = retval.fortran_vec ();

  F77_FCN (cffti, CFFTI) (npts, pwsave);

  for (int j = 0; j < nsamples; j++)
    F77_FCN (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);

  for (int j = 0; j < npts*nsamples; j++)
    tmp_data[j] = tmp_data[j] / (double) npts;

  npts = nc;
  nsamples = nr;
  nn = 4*npts+15;

  wsave.resize (nn);
  pwsave = wsave.fortran_vec ();

  Array<Complex> row (npts);
  Complex *prow = row.fortran_vec ();

  F77_FCN (cffti, CFFTI) (npts, pwsave);

  for (int j = 0; j < nsamples; j++)
    {
      for (int i = 0; i < npts; i++)
	prow[i] = tmp_data[i*nr + j];

      F77_FCN (cfftb, CFFTB) (npts, prow, pwsave);

      for (int i = 0; i < npts; i++)
	tmp_data[i*nr + j] = prow[i] / (double) npts;
    }

  return retval;
}

DET
Matrix::determinant (void) const
{
  int info;
  double rcond;
  return determinant (info, rcond);
}

DET
Matrix::determinant (int& info) const
{
  double rcond;
  return determinant (info, rcond);
}

DET
Matrix::determinant (int& info, double& rcond) const
{
  DET retval;

  int nr = rows ();
  int nc = cols ();

  if (nr == 0 || nc == 0)
    {
      double d[2];
      d[0] = 1.0;
      d[1] = 0.0;
      retval = DET (d);
    }
  else
    {
      info = 0;

      Array<int> ipvt (nr);
      int *pipvt = ipvt.fortran_vec ();

      Array<double> z (nr);
      double *pz = z.fortran_vec ();

      Matrix atmp = *this;
      double *tmp_data = atmp.fortran_vec ();

      F77_XFCN (dgeco, DGECO, (tmp_data, nr, nr, pipvt, rcond, pz));

      if (f77_exception_encountered)
	(*current_liboctave_error_handler) ("unrecoverable error in dgeco");
      else
	{
	  volatile double rcond_plus_one = rcond + 1.0;

	  if (rcond_plus_one == 1.0)
	    {
	      info = -1;
	      retval = DET ();
	    }
	  else
	    {
	      double d[2];

	      F77_XFCN (dgedi, DGEDI, (tmp_data, nr, nr, pipvt, d, pz, 10));

	      if (f77_exception_encountered)
		(*current_liboctave_error_handler)
		  ("unrecoverable error in dgedi");
	      else
		retval = DET (d);
	    }
	}
    }

  return retval;
}

Matrix
Matrix::solve (const Matrix& b) const
{
  int info;
  double rcond;
  return solve (b, info, rcond);
}

Matrix
Matrix::solve (const Matrix& b, int& info) const
{
  double rcond;
  return solve (b, info, rcond);
}

Matrix
Matrix::solve (const Matrix& b, int& info, double& rcond) const
{
  Matrix retval;

  int nr = rows ();
  int nc = cols ();

  if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
    (*current_liboctave_error_handler)
      ("matrix dimension mismatch solution of linear equations");
  else
    {
      info = 0;

      Array<int> ipvt (nr);
      int *pipvt = ipvt.fortran_vec ();

      Array<double> z (nr);
      double *pz = z.fortran_vec ();

      Matrix atmp = *this;
      double *tmp_data = atmp.fortran_vec ();

      F77_XFCN (dgeco, DGECO, (tmp_data, nr, nr, pipvt, rcond, pz));

      if (f77_exception_encountered)
	(*current_liboctave_error_handler) ("unrecoverable error in dgeco");
      else
	{
	  volatile double rcond_plus_one = rcond + 1.0;

	  if (rcond_plus_one == 1.0)
	    {
	      info = -2;
	    }
	  else
	    {
	      retval = b;
	      double *result = retval.fortran_vec ();

	      int b_nc = b.cols ();

	      for (volatile int j = 0; j < b_nc; j++)
		{
		  F77_XFCN (dgesl, DGESL, (tmp_data, nr, nr, pipvt,
					   &result[nr*j], 0)); 

		  if (f77_exception_encountered)
		    {
		      (*current_liboctave_error_handler)
			("unrecoverable error in dgesl");

		      break;
		    }
		}
	    }
	}
    }

  return retval;
}

ComplexMatrix
Matrix::solve (const ComplexMatrix& b) const
{
  ComplexMatrix tmp (*this);
  return tmp.solve (b);
}

ComplexMatrix
Matrix::solve (const ComplexMatrix& b, int& info) const
{
  ComplexMatrix tmp (*this);
  return tmp.solve (b, info);
}

ComplexMatrix
Matrix::solve (const ComplexMatrix& b, int& info, double& rcond) const
{
  ComplexMatrix tmp (*this);
  return tmp.solve (b, info, rcond);
}

ColumnVector
Matrix::solve (const ColumnVector& b) const
{
  int info; double rcond;
  return solve (b, info, rcond);
}

ColumnVector
Matrix::solve (const ColumnVector& b, int& info) const
{
  double rcond;
  return solve (b, info, rcond);
}

ColumnVector
Matrix::solve (const ColumnVector& b, int& info, double& rcond) const
{
  ColumnVector retval;

  int nr = rows ();
  int nc = cols ();

  if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
    (*current_liboctave_error_handler)
      ("matrix dimension mismatch solution of linear equations");
  else
    {
      info = 0;

      Array<int> ipvt (nr);
      int *pipvt = ipvt.fortran_vec ();

      Array<double> z (nr);
      double *pz = z.fortran_vec ();

      Matrix atmp = *this;
      double *tmp_data = atmp.fortran_vec ();

      F77_XFCN (dgeco, DGECO, (tmp_data, nr, nr, pipvt, rcond, pz));

      if (f77_exception_encountered)
	(*current_liboctave_error_handler)
	  ("unrecoverable error in dgeco");
      else
	{
	  volatile double rcond_plus_one = rcond + 1.0;

	  if (rcond_plus_one == 1.0)
	    {
	      info = -2;
	    }
	  else
	    {
	      retval = b;
	      double *result = retval.fortran_vec ();

	      F77_XFCN (dgesl, DGESL, (tmp_data, nr, nr, pipvt, result, 0));

	      if (f77_exception_encountered)
		(*current_liboctave_error_handler)
		  ("unrecoverable error in dgesl");
	    }
	}
    }

  return retval;
}

ComplexColumnVector
Matrix::solve (const ComplexColumnVector& b) const
{
  ComplexMatrix tmp (*this);
  return tmp.solve (b);
}

ComplexColumnVector
Matrix::solve (const ComplexColumnVector& b, int& info) const
{
  ComplexMatrix tmp (*this);
  return tmp.solve (b, info);
}

ComplexColumnVector
Matrix::solve (const ComplexColumnVector& b, int& info, double& rcond) const
{
  ComplexMatrix tmp (*this);
  return tmp.solve (b, info, rcond);
}

Matrix
Matrix::lssolve (const Matrix& b) const
{
  int info;
  int rank;
  return lssolve (b, info, rank);
}

Matrix
Matrix::lssolve (const Matrix& b, int& info) const
{
  int rank;
  return lssolve (b, info, rank);
}

Matrix
Matrix::lssolve (const Matrix& b, int& info, int& rank) const
{
  Matrix retval;

  int nrhs = b.cols ();

  int m = rows ();
  int n = cols ();

  if (m == 0 || n == 0 || m != b.rows ())
    (*current_liboctave_error_handler)
      ("matrix dimension mismatch in solution of least squares problem");
  else
    {
      Matrix atmp = *this;
      double *tmp_data = atmp.fortran_vec ();

      int nrr = m > n ? m : n;
      Matrix result (nrr, nrhs);

      for (int j = 0; j < nrhs; j++)
	for (int i = 0; i < m; i++)
	  result.elem (i, j) = b.elem (i, j);

      double *presult = result.fortran_vec ();

      int len_s = m < n ? m : n;
      Array<double> s (len_s);
      double *ps = s.fortran_vec ();

      double rcond = -1.0;

      int lwork;
      if (m < n)
	lwork = 3*m + (2*m > nrhs
		       ? (2*m > n ? 2*m : n)
		       : (nrhs > n ? nrhs : n));
      else
	lwork = 3*n + (2*n > nrhs
		       ? (2*n > m ? 2*n : m)
		       : (nrhs > m ? nrhs : m));

      Array<double> work (lwork);
      double *pwork = work.fortran_vec ();

      F77_XFCN (dgelss, DGELSS, (m, n, nrhs, tmp_data, m, presult, nrr, ps,
				 rcond, rank, pwork, lwork, info));

      if (f77_exception_encountered)
	(*current_liboctave_error_handler) ("unrecoverable error in dgelss");
      else
	{
	  retval.resize (n, nrhs);
	  for (int j = 0; j < nrhs; j++)
	    for (int i = 0; i < n; i++)
	      retval.elem (i, j) = result.elem (i, j);
	}
    }

  return retval;
}

ComplexMatrix
Matrix::lssolve (const ComplexMatrix& b) const
{
  ComplexMatrix tmp (*this);
  int info;
  int rank;
  return tmp.lssolve (b, info, rank);
}

ComplexMatrix
Matrix::lssolve (const ComplexMatrix& b, int& info) const
{
  ComplexMatrix tmp (*this);
  int rank;
  return tmp.lssolve (b, info, rank);
}

ComplexMatrix
Matrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const
{
  ComplexMatrix tmp (*this);
  return tmp.lssolve (b, info, rank);
}

ColumnVector
Matrix::lssolve (const ColumnVector& b) const
{
  int info;
  int rank;
  return lssolve (b, info, rank);
}

ColumnVector
Matrix::lssolve (const ColumnVector& b, int& info) const
{
  int rank;
  return lssolve (b, info, rank);
}

ColumnVector
Matrix::lssolve (const ColumnVector& b, int& info, int& rank) const
{
  ColumnVector retval;

  int nrhs = 1;

  int m = rows ();
  int n = cols ();

  if (m == 0 || n == 0 || m != b.length ())
    (*current_liboctave_error_handler)
      ("matrix dimension mismatch in solution of least squares problem");
  else
    {
      Matrix atmp = *this;
      double *tmp_data = atmp.fortran_vec ();

      int nrr = m > n ? m : n;
      ColumnVector result (nrr);

      for (int i = 0; i < m; i++)
	result.elem (i) = b.elem (i);

      double *presult = result.fortran_vec ();

      int len_s = m < n ? m : n;
      Array<double> s (len_s);
      double *ps = s.fortran_vec ();

      double rcond = -1.0;

      int lwork;
      if (m < n)
	lwork = 3*m + (2*m > nrhs
		       ? (2*m > n ? 2*m : n)
		       : (nrhs > n ? nrhs : n));
      else
	lwork = 3*n + (2*n > nrhs
		       ? (2*n > m ? 2*n : m)
		       : (nrhs > m ? nrhs : m));

      Array<double> work (lwork);
      double *pwork = work.fortran_vec ();

      F77_XFCN (dgelss, DGELSS, (m, n, nrhs, tmp_data, m, presult, nrr,
				 ps, rcond, rank, pwork, lwork, info));

      if (f77_exception_encountered)
	(*current_liboctave_error_handler) ("unrecoverable error in dgelss");
      else
	{
	  retval.resize (n);
	  for (int i = 0; i < n; i++)
	    retval.elem (i) = result.elem (i);
	}
    }

  return retval;
}

ComplexColumnVector
Matrix::lssolve (const ComplexColumnVector& b) const
{
  ComplexMatrix tmp (*this);
  return tmp.lssolve (b);
}

ComplexColumnVector
Matrix::lssolve (const ComplexColumnVector& b, int& info) const
{
  ComplexMatrix tmp (*this);
  return tmp.lssolve (b, info);
}

ComplexColumnVector
Matrix::lssolve (const ComplexColumnVector& b, int& info, int& rank) const
{
  ComplexMatrix tmp (*this);
  return tmp.lssolve (b, info, rank);
}

// Constants for matrix exponential calculation.

static double padec [] =
{
  5.0000000000000000e-1,
  1.1666666666666667e-1,
  1.6666666666666667e-2,
  1.6025641025641026e-3,
  1.0683760683760684e-4,
  4.8562548562548563e-6,
  1.3875013875013875e-7,
  1.9270852604185938e-9,
};

Matrix
Matrix::expm (void) const
{
  Matrix retval;

  Matrix m = *this;

  int nc = columns ();

  // trace shift value
  double trshift = 0;

  // Preconditioning step 1: trace normalization.

  for (int i = 0; i < nc; i++)
    trshift += m.elem (i, i);

  trshift /= nc;

  if(trshift < 0.0)
    trshift = 0.0;	// don't want stable eigenvalues to become unstable
  else
    for (int i = 0; i < nc; i++)
      m.elem (i, i) -= trshift;

  // Preconditioning step 2: balancing.

  AEPBALANCE mbal (m, "B");
  m = mbal.balanced_matrix ();
  Matrix d = mbal.balancing_matrix ();

  // Preconditioning step 3: scaling.

  ColumnVector work(nc);
  double inf_norm
    = F77_FCN (dlange, DLANGE) ("I", nc, nc, m.fortran_vec (),nc,
				work.fortran_vec ());

  int sqpow = (int) (inf_norm > 0.0
		     ? (1.0 + log (inf_norm) / log (2.0))
		     : 0.0);

  // Check whether we need to square at all.

  if (sqpow < 0)
    sqpow = 0;

  if (sqpow > 0)
    {
      double scale_factor = 1.0;
      for (int i = 0; i < sqpow; i++)
	scale_factor *= 2.0;

      m = m / scale_factor;
    }

  // npp, dpp: pade' approx polynomial matrices.

  Matrix npp (nc, nc, 0.0);
  Matrix dpp = npp;

  // Now powers a^8 ... a^1.

  int minus_one_j = -1;
  for (int j = 7; j >= 0; j--)
    {
      npp = m * npp + m * padec[j];
      dpp = m * dpp + m * (minus_one_j * padec[j]);
      minus_one_j *= -1;
    }

  // Zero power.

  dpp = -dpp;
  for(int j = 0; j < nc; j++)
    {
      npp.elem (j, j) += 1.0;
      dpp.elem (j, j) += 1.0;
    }

  // Compute pade approximation = inverse (dpp) * npp.

  retval = dpp.solve (npp);

  // Reverse preconditioning step 3: repeated squaring.

  while (sqpow)
    {
      retval = retval * retval;
      sqpow--;
    }

  // Reverse preconditioning step 2: inverse balancing.

  retval = retval.transpose();
  d = d.transpose ();
  retval = retval * d;
  retval = d.solve (retval);
  retval = retval.transpose ();

  // Reverse preconditioning step 1: fix trace normalization.

  return retval * exp (trshift);
}

Matrix&
Matrix::operator += (const Matrix& a)
{
  int nr = rows ();
  int nc = cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
      return *this;
    }

  if (nr == 0 || nc == 0)
    return *this;

  double *d = fortran_vec (); // Ensures only one reference to my privates!

  add2 (d, a.data (), length ());

  return *this;
}

Matrix&
Matrix::operator -= (const Matrix& a)
{
  int nr = rows ();
  int nc = cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
      return *this;
    }

  if (nr == 0 || nc == 0)
    return *this;

  double *d = fortran_vec (); // Ensures only one reference to my privates!

  subtract2 (d, a.data (), length ());

  return *this;
}

Matrix&
Matrix::operator += (const DiagMatrix& a)
{
  int nr = rows ();
  int nc = cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
      return *this;
    }

  for (int i = 0; i < a.length (); i++)
    elem (i, i) += a.elem (i, i);

  return *this;
}

Matrix&
Matrix::operator -= (const DiagMatrix& a)
{
  int nr = rows ();
  int nc = cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
      return *this;
    }

  for (int i = 0; i < a.length (); i++)
    elem (i, i) -= a.elem (i, i);

  return *this;
}

// unary operations

Matrix
Matrix::operator ! (void) const
{
  int nr = rows ();
  int nc = cols ();

  Matrix b (nr, nc);

  for (int j = 0; j < nc; j++)
    for (int i = 0; i < nr; i++)
      b.elem (i, j) = ! elem (i, j);

  return b;
}

// column vector by row vector -> matrix operations

Matrix
operator * (const ColumnVector& v, const RowVector& a)
{
  Matrix retval;

  int len = v.length ();
  int a_len = a.length ();

  if (len != a_len)
    gripe_nonconformant ("operator *", len, 1, 1, a_len);
  else
    {
      if (len != 0)
	{
	  retval.resize (len, a_len);
	  double *c = retval.fortran_vec ();

	  F77_XFCN (dgemm, DGEMM, ("N", "N", len, a_len, 1, 1.0,
				   v.data (), len, a.data (), 1, 0.0,
				   c, len, 1L, 1L));

	  if (f77_exception_encountered)
	    (*current_liboctave_error_handler)
	      ("unrecoverable error in dgemm");
	}
    }

  return retval;
}

// diagonal matrix by scalar -> matrix operations

Matrix
operator + (const DiagMatrix& a, double s)
{
  Matrix tmp (a.rows (), a.cols (), s);
  return a + tmp;
}

Matrix
operator - (const DiagMatrix& a, double s)
{
  Matrix tmp (a.rows (), a.cols (), -s);
  return a + tmp;
}

// scalar by diagonal matrix -> matrix operations

Matrix
operator + (double s, const DiagMatrix& a)
{
  Matrix tmp (a.rows (), a.cols (), s);
  return tmp + a;
}

Matrix
operator - (double s, const DiagMatrix& a)
{
  Matrix tmp (a.rows (), a.cols (), s);
  return tmp - a;
}

// matrix by diagonal matrix -> matrix operations

Matrix
operator + (const Matrix& m, const DiagMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
      return Matrix ();
    }

  if (nr == 0 || nc == 0)
    return Matrix (nr, nc);

  Matrix result (m);
  int a_len = a.length ();
  for (int i = 0; i < a_len; i++)
    result.elem (i, i) += a.elem (i, i);

  return result;
}

Matrix
operator - (const Matrix& m, const DiagMatrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
      return Matrix ();
    }

  if (nr == 0 || nc == 0)
    return Matrix (nr, nc);

  Matrix result (m);
  int a_len = a.length ();
  for (int i = 0; i < a_len; i++)
    result.elem (i, i) -= a.elem (i, i);

  return result;
}

Matrix
operator * (const Matrix& m, const DiagMatrix& a)
{
  Matrix retval;

  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nc != a_nr)
    gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  else
    {
      if (nr == 0 || nc == 0 || a_nc == 0)
	retval.resize (nr, a_nc, 0.0);
      else
	{
	  retval.resize (nr, a_nc);
	  double *c = retval.fortran_vec ();

	  double *ctmp = 0;

	  int a_len = a.length ();

	  for (int j = 0; j < a_len; j++)
	    {
	      int idx = j * nr;
	      ctmp = c + idx;

	      if (a.elem (j, j) == 1.0)
		{
		  for (int i = 0; i < nr; i++)
		    ctmp[i] = m.elem (i, j);
		}
	      else if (a.elem (j, j) == 0.0)
		{
		  for (int i = 0; i < nr; i++)
		    ctmp[i] = 0.0;
		}
	      else
		{
		  for (int i = 0; i < nr; i++)
		    ctmp[i] = a.elem (j, j) * m.elem (i, j);
		}
	    }

	  if (a_nr < a_nc)
	    {
	      for (int i = nr * nc; i < nr * a_nc; i++)
		ctmp[i] = 0.0;
	    }
	}
    }

  return retval;
}

// diagonal matrix by matrix -> matrix operations

Matrix
operator + (const DiagMatrix& m, const Matrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
      return Matrix ();
    }

  if (nr == 0 || nc == 0)
    return Matrix (nr, nc);

  Matrix result (a);
  for (int i = 0; i < m.length (); i++)
    result.elem (i, i) += m.elem (i, i);

  return result;
}

Matrix
operator - (const DiagMatrix& m, const Matrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nr != a_nr || nc != a_nc)
    {
      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
      return Matrix ();
    }

  if (nr == 0 || nc == 0)
    return Matrix (nr, nc);

  Matrix result (-a);
  for (int i = 0; i < m.length (); i++)
    result.elem (i, i) += m.elem (i, i);

  return result;
}

Matrix
operator * (const DiagMatrix& m, const Matrix& a)
{
  int nr = m.rows ();
  int nc = m.cols ();
  int a_nr = a.rows ();
  int a_nc = a.cols ();
  if (nc != a_nr)
    {
      gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
      return Matrix ();
    }

  if (nr == 0 || nc == 0 || a_nc == 0)
    return Matrix (nr, a_nc, 0.0);

  Matrix c (nr, a_nc);

  for (int i = 0; i < m.length (); i++)
    {
      if (m.elem (i, i) == 1.0)
	{
	  for (int j = 0; j < a_nc; j++)
	    c.elem (i, j) = a.elem (i, j);
	}
      else if (m.elem (i, i) == 0.0)
	{
	  for (int j = 0; j < a_nc; j++)
	    c.elem (i, j) = 0.0;
	}
      else
	{
	  for (int j = 0; j < a_nc; j++)
	    c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
	}
    }

  if (nr > nc)
    {
      for (int j = 0; j < a_nc; j++)
	for (int i = a_nr; i < nr; i++)
	  c.elem (i, j) = 0.0;
    }

  return c;
}

// matrix by matrix -> matrix operations

Matrix
operator * (const Matrix& m, const Matrix& a)
{
  Matrix retval;

  int nr = m.rows ();
  int nc = m.cols ();

  int a_nr = a.rows ();
  int a_nc = a.cols ();

  if (nc != a_nr)
    gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  else
    {
      if (nr == 0 || nc == 0 || a_nc == 0)
	retval.resize (nr, a_nc, 0.0);
      else
	{
	  int ld  = nr;
	  int lda = a_nr;

	  retval.resize (nr, a_nc);
	  double *c = retval.fortran_vec ();

	  F77_XFCN (dgemm, DGEMM, ("N", "N", nr, a_nc, nc, 1.0,
				   m.data (), ld, a.data (), lda, 0.0,
				   c, nr, 1L, 1L));

	  if (f77_exception_encountered)
	    (*current_liboctave_error_handler)
	      ("unrecoverable error in dgemm");
	}
    }

  return retval;
}

// other operations.

Matrix
Matrix::map (d_d_Mapper f) const
{
  Matrix b (*this);
  return b.apply (f);
}

Matrix&
Matrix::apply (d_d_Mapper f)
{
  double *d = fortran_vec (); // Ensures only one reference to my privates!

  for (int i = 0; i < length (); i++)
    d[i] = f (d[i]);

  return *this;
}

bool
Matrix::any_element_is_negative (void) const
{
  int nr = rows ();
  int nc = cols ();

  for (int j = 0; j < nc; j++)
    for (int i = 0; i < nr; i++)
      if (elem (i, j) < 0.0)
	return true;

  return false;
}


bool
Matrix::any_element_is_inf_or_nan (void) const
{
  int nr = rows ();
  int nc = cols ();

  for (int j = 0; j < nc; j++)
    for (int i = 0; i < nr; i++)
      {
	double val = elem (i, j);
	if (xisinf (val) || xisnan (val))
	  return 1;
      }
  return 0;
}

bool
Matrix::all_elements_are_int_or_inf_or_nan (void) const
{
  int nr = rows ();
  int nc = cols ();

  for (int j = 0; j < nc; j++)
    for (int i = 0; i < nr; i++)
      {
	double val = elem (i, j);
	if (xisnan (val) || D_NINT (val) == val)
	  continue;
	else
	  return false;
      }

  return true;
}

// Return nonzero if any element of M is not an integer.  Also extract
// the largest and smallest values and return them in MAX_VAL and MIN_VAL.

bool
Matrix::all_integers (double& max_val, double& min_val) const
{
  int nr = rows ();
  int nc = cols ();

  if (nr > 0 && nc > 0)
    {
      max_val = elem (0, 0);
      min_val = elem (0, 0);
    }
  else
    return false;

  for (int j = 0; j < nc; j++)
    for (int i = 0; i < nr; i++)
      {
	double val = elem (i, j);

	if (val > max_val)
	  max_val = val;

	if (val < min_val)
	  min_val = val;

	if (D_NINT (val) != val)
	  return false;
      }

  return true;
}

bool
Matrix::too_large_for_float (void) const
{
  int nr = rows ();
  int nc = cols ();

  for (int j = 0; j < nc; j++)
    for (int i = 0; i < nr; i++)
      {
	double val = elem (i, j);

	if (val > FLT_MAX || val < FLT_MIN)
	  return true;
      }

  return false;
}

// XXX FIXME XXX Do these really belong here?  They should maybe be
// cleaned up a bit, no?  What about corresponding functions for the
// Vectors?

Matrix
Matrix::all (void) const
{
  int nr = rows ();
  int nc = cols ();
  Matrix retval;
  if (nr > 0 && nc > 0)
    {
      if (nr == 1)
	{
	  retval.resize (1, 1);
	  retval.elem (0, 0) = 1.0;
	  for (int j = 0; j < nc; j++)
	    {
	      if (elem (0, j) == 0.0)
		{
		  retval.elem (0, 0) = 0.0;
		  break;
		}
	    }
	}
      else if (nc == 1)
	{
	  retval.resize (1, 1);
	  retval.elem (0, 0) = 1.0;
	  for (int i = 0; i < nr; i++)
	    {
	      if (elem (i, 0) == 0.0)
		{
		  retval.elem (0, 0) = 0.0;
		  break;
		}
	    }
	}
      else
	{
	  retval.resize (1, nc);
	  for (int j = 0; j < nc; j++)
	    {
	      retval.elem (0, j) = 1.0;
	      for (int i = 0; i < nr; i++)
		{
		  if (elem (i, j) == 0.0)
		    {
		      retval.elem (0, j) = 0.0;
		      break;
		    }
		}
	    }
	}
    }
  return retval;
}

Matrix
Matrix::any (void) const
{
  int nr = rows ();
  int nc = cols ();
  Matrix retval;
  if (nr > 0 && nc > 0)
    {
      if (nr == 1)
	{
	  retval.resize (1, 1);
	  retval.elem (0, 0) = 0.0;
	  for (int j = 0; j < nc; j++)
	    {
	      if (elem (0, j) != 0.0)
		{
		  retval.elem (0, 0) = 1.0;
		  break;
		}
	    }
	}
      else if (nc == 1)
	{
	  retval.resize (1, 1);
	  retval.elem (0, 0) = 0.0;
	  for (int i = 0; i < nr; i++)
	    {
	      if (elem (i, 0) != 0.0)
		{
		  retval.elem (0, 0) = 1.0;
		  break;
		}
	    }
	}
      else
	{
	  retval.resize (1, nc);
	  for (int j = 0; j < nc; j++)
	    {
	      retval.elem (0, j) = 0.0;
	      for (int i = 0; i < nr; i++)
		{
		  if (elem (i, j) != 0.0)
		    {
		      retval.elem (0, j) = 1.0;
		      break;
		    }
		}
	    }
	}
    }
  return retval;
}

Matrix
Matrix::cumprod (void) const
{
  Matrix retval;

  int nr = rows ();
  int nc = cols ();

  if (nr == 1)
    {
      retval.resize (1, nc);
      if (nc > 0)
	{
	  double prod = elem (0, 0);
	  for (int j = 0; j < nc; j++)
	    {
	      retval.elem (0, j) = prod;
	      if (j < nc - 1)
		prod *= elem (0, j+1);
	    }
	}
    }
  else if (nc == 1)
    {
      retval.resize (nr, 1);
      if (nr > 0)
	{
	  double prod = elem (0, 0);
	  for (int i = 0; i < nr; i++)
	    {
	      retval.elem (i, 0) = prod;
	      if (i < nr - 1)
		prod *= elem (i+1, 0);
	    }
	}
    }
  else
    {
      retval.resize (nr, nc);
      if (nr > 0 && nc > 0)
	{
	  for (int j = 0; j < nc; j++)
	    {
	      double prod = elem (0, j);
	      for (int i = 0; i < nr; i++)
		{
		  retval.elem (i, j) = prod;
		  if (i < nr - 1)
		    prod *= elem (i+1, j);
		}
	    }
	}
    }
  return retval;
}

Matrix
Matrix::cumsum (void) const
{
  Matrix retval;

  int nr = rows ();
  int nc = cols ();

  if (nr == 1)
    {
      retval.resize (1, nc);
      if (nc > 0)
	{
	  double sum = elem (0, 0);
	  for (int j = 0; j < nc; j++)
	    {
	      retval.elem (0, j) = sum;
	      if (j < nc - 1)
		sum += elem (0, j+1);
	    }
	}
    }
  else if (nc == 1)
    {
      retval.resize (nr, 1);
      if (nr > 0)
	{
	  double sum = elem (0, 0);
	  for (int i = 0; i < nr; i++)
	    {
	      retval.elem (i, 0) = sum;
	      if (i < nr - 1)
		sum += elem (i+1, 0);
	    }
	}
    }
  else
    {
      retval.resize (nr, nc);
      if (nr > 0 && nc > 0)
	{
	  for (int j = 0; j < nc; j++)
	    {
	      double sum = elem (0, j);
	      for (int i = 0; i < nr; i++)
		{
		  retval.elem (i, j) = sum;
		  if (i < nr - 1)
		    sum += elem (i+1, j);
		}
	    }
	}
    }
  return retval;
}

Matrix
Matrix::prod (void) const
{
  Matrix retval;

  int nr = rows ();
  int nc = cols ();

  if (nr == 1)
    {
      retval.resize (1, 1);
      retval.elem (0, 0) = 1.0;
      for (int j = 0; j < nc; j++)
	retval.elem (0, 0) *= elem (0, j);
    }
  else if (nc == 1)
    {
      retval.resize (1, 1);
      retval.elem (0, 0) = 1.0;
      for (int i = 0; i < nr; i++)
	retval.elem (0, 0) *= elem (i, 0);
    }
  else
    {
      if (nc == 0)
	{
	  retval.resize (1, 1);
	  retval.elem (0, 0) = 1.0;
	}
      else
	retval.resize (1, nc);

      for (int j = 0; j < nc; j++)
	{
	  retval.elem (0, j) = 1.0;
	  for (int i = 0; i < nr; i++)
	    retval.elem (0, j) *= elem (i, j);
	}
    }
  return retval;
}

Matrix
Matrix::sum (void) const
{
  Matrix retval;

  int nr = rows ();
  int nc = cols ();

  if (nr == 1)
    {
      retval.resize (1, 1);
      retval.elem (0, 0) = 0.0;
      for (int j = 0; j < nc; j++)
	retval.elem (0, 0) += elem (0, j);
    }
  else if (nc == 1)
    {
      retval.resize (1, 1);
      retval.elem (0, 0) = 0.0;
      for (int i = 0; i < nr; i++)
	retval.elem (0, 0) += elem (i, 0);
    }
  else
    {
      if (nc == 0)
	{
	  retval.resize (1, 1);
	  retval.elem (0, 0) = 0.0;
	}
      else
	retval.resize (1, nc);

      for (int j = 0; j < nc; j++)
	{
	  retval.elem (0, j) = 0.0;
	  for (int i = 0; i < nr; i++)
	    retval.elem (0, j) += elem (i, j);
	}
    }
  return retval;
}

Matrix
Matrix::sumsq (void) const
{
  Matrix retval;

  int nr = rows ();
  int nc = cols ();

  if (nr == 1)
    {
      retval.resize (1, 1);
      retval.elem (0, 0) = 0.0;
      for (int j = 0; j < nc; j++)
	{
	  double d = elem (0, j);
	  retval.elem (0, 0) += d * d;
	}
    }
  else if (nc == 1)
    {
      retval.resize (1, 1);
      retval.elem (0, 0) = 0.0;
      for (int i = 0; i < nr; i++)
	{
	  double d = elem (i, 0);
	  retval.elem (0, 0) += d * d;
	}
    }
  else
    {
      retval.resize (1, nc);
      for (int j = 0; j < nc; j++)
	{
	  retval.elem (0, j) = 0.0;
	  for (int i = 0; i < nr; i++)
	    {
	      double d = elem (i, j);
	      retval.elem (0, j) += d * d;
	    }
	}
    }
  return retval;
}

Matrix
Matrix::abs (void) const
{
  int nr = rows ();
  int nc = cols ();

  Matrix retval (nr, nc);

  for (int j = 0; j < nc; j++)
    for (int i = 0; i < nr; i++)
      retval (i, j) = fabs (elem (i, j));

  return retval;
}

ColumnVector
Matrix::diag (void) const
{
  return diag (0);
}

ColumnVector
Matrix::diag (int k) const
{
  int nnr = rows ();
  int nnc = cols ();
  if (k > 0)
    nnc -= k;
  else if (k < 0)
    nnr += k;

  ColumnVector d;

  if (nnr > 0 && nnc > 0)
    {
      int ndiag = (nnr < nnc) ? nnr : nnc;

      d.resize (ndiag);

      if (k > 0)
	{
	  for (int i = 0; i < ndiag; i++)
	    d.elem (i) = elem (i, i+k);
	}
      else if ( k < 0)
	{
	  for (int i = 0; i < ndiag; i++)
	    d.elem (i) = elem (i-k, i);
	}
      else
	{
	  for (int i = 0; i < ndiag; i++)
	    d.elem (i) = elem (i, i);
	}
    }
  else
    cerr << "diag: requested diagonal out of range\n";

  return d;
}

ColumnVector
Matrix::row_min (void) const
{
  Array<int> index;
  return row_min (index);
}

ColumnVector
Matrix::row_min (Array<int>& index) const
{
  ColumnVector result;

  int nr = rows ();
  int nc = cols ();

  if (nr > 0 && nc > 0)
    {
      result.resize (nr);
      index.resize (nr);

      for (int i = 0; i < nr; i++)
        {
	  int idx = 0;

	  double tmp_min = elem (i, idx);

	  if (xisnan (tmp_min))
	    idx = -1;
	  else
	    {
	      for (int j = 1; j < nc; j++)
		{
		  double tmp = elem (i, j);

		  if (xisnan (tmp))
		    {
		      idx = -1;
		      break;
		    }
		  else if (tmp < tmp_min)
		    {
		      idx = j;
		      tmp_min = tmp;
		    }
		}
	    }

	  result.elem (i) = (idx < 0) ? octave_NaN : tmp_min;
	  index.elem (i) = idx;
        }
    }

  return result;
}

ColumnVector
Matrix::row_max (void) const
{
  Array<int> index;
  return row_max (index);
}

ColumnVector
Matrix::row_max (Array<int>& index) const
{
  ColumnVector result;

  int nr = rows ();
  int nc = cols ();

  if (nr > 0 && nc > 0)
    {
      result.resize (nr);
      index.resize (nr);

      for (int i = 0; i < nr; i++)
        {
	  int idx = 0;

	  double tmp_max = elem (i, idx);

	  if (xisnan (tmp_max))
	    idx = -1;
	  else
	    {
	      for (int j = 1; j < nc; j++)
		{
		  double tmp = elem (i, j);

		  if (xisnan (tmp))
		    {
		      idx = -1;
		      break;
		    }
		  else if (tmp > tmp_max)
		    {
		      idx = j;
		      tmp_max = tmp;
		    }
		}
	    }

	  result.elem (i) = (idx < 0) ? octave_NaN : tmp_max;
	  index.elem (i) = idx;
        }
    }

  return result;
}

RowVector
Matrix::column_min (void) const
{
  Array<int> index;
  return column_min (index);
}

RowVector
Matrix::column_min (Array<int>& index) const
{
  RowVector result;

  int nr = rows ();
  int nc = cols ();

  if (nr > 0 && nc > 0)
    {
      result.resize (nc);
      index.resize (nc);

      for (int j = 0; j < nc; j++)
        {
	  int idx = 0;

	  double tmp_min = elem (idx, j);

	  if (xisnan (tmp_min))
	    idx = -1;
	  else
	    {
	      for (int i = 1; i < nr; i++)
		{
		  double tmp = elem (i, j);

		  if (xisnan (tmp))
		    {
		      idx = -1;
		      break;
		    }
		  else if (tmp < tmp_min)
		    {
		      idx = i;
		      tmp_min = tmp;
		    }
		}
	    }

	  result.elem (j) = (idx < 0) ? octave_NaN : tmp_min;
	  index.elem (j) = idx;
        }
    }

  return result;
}

RowVector
Matrix::column_max (void) const
{
  Array<int> index;
  return column_max (index);
}

RowVector
Matrix::column_max (Array<int>& index) const
{
  RowVector result;

  int nr = rows ();
  int nc = cols ();

  if (nr > 0 && nc > 0)
    {
      result.resize (nc);
      index.resize (nc);

      for (int j = 0; j < nc; j++)
        {
	  int idx = 0;

	  double tmp_max = elem (idx, j);

	  if (xisnan (tmp_max))
	    idx = -1;
	  else
	    {
	      for (int i = 1; i < nr; i++)
		{
		  double tmp = elem (i, j);

		  if (xisnan (tmp))
		    {
		      idx = -1;
		      break;
		    }
		  else if (tmp > tmp_max)
		    {
		      idx = i;
		      tmp_max = tmp;
		    }
		}
	    }

	  result.elem (j) = (idx < 0) ? octave_NaN : tmp_max;
	  index.elem (j) = idx;
        }
    }

  return result;
}

ostream&
operator << (ostream& os, const Matrix& a)
{
//  int field_width = os.precision () + 7;

  for (int i = 0; i < a.rows (); i++)
    {
      for (int j = 0; j < a.cols (); j++)
	os << " " /* setw (field_width) */ << a.elem (i, j);
      os << "\n";
    }
  return os;
}

istream&
operator >> (istream& is, Matrix& a)
{
  int nr = a.rows ();
  int nc = a.cols ();

  if (nr < 1 || nc < 1)
    is.clear (ios::badbit);
  else
    {
      double tmp;
      for (int i = 0; i < nr; i++)
	for (int j = 0; j < nc; j++)
	  {
	    is >> tmp;
	    if (is)
	      a.elem (i, j) = tmp;
	    else
	      goto done;
	  }
    }

done:

  return is;
}

template <class T>
static void
read_int (istream& is, bool swap_bytes, T& val)
{
  is.read ((char *) &val, sizeof (T));

  if (swap_bytes)
    {
      switch (sizeof (T))
	{
	case 1:
	  break;

	case 2:
	  swap_2_bytes ((char *) &val);
	  break;

	case 4:
	  swap_4_bytes ((char *) &val);
	  break;

	case 8:
	  swap_8_bytes ((char *) &val);
	  break;

	default:
	  (*current_liboctave_error_handler)
	    ("read_int: unrecognized data format!");
	}
    }
}

template void read_int (istream&, bool, char&);
template void read_int (istream&, bool, signed char&);
template void read_int (istream&, bool, unsigned char&);
template void read_int (istream&, bool, short&);
template void read_int (istream&, bool, unsigned short&);
template void read_int (istream&, bool, int&);
template void read_int (istream&, bool, unsigned int&);
template void read_int (istream&, bool, long&);
template void read_int (istream&, bool, unsigned long&);

static inline bool
do_read (istream& is, oct_data_conv::data_type dt, 
	 oct_mach_info::float_format flt_fmt, bool swap_bytes,
	 bool do_float_conversion, double& val)
{
  bool retval = true;

  switch (dt)
    {
    case oct_data_conv::dt_char:
      {
	char tmp;
	read_int (is, swap_bytes, tmp);
	val = tmp;
      }
      break;

    case oct_data_conv::dt_schar:
      {
	signed char tmp;
	read_int (is, swap_bytes, tmp);
	val = tmp;
      }
      break;

    case oct_data_conv::dt_uchar:
      {
	unsigned char tmp;
	read_int (is, swap_bytes, tmp);
	val = tmp;
      }
      break;

    case oct_data_conv::dt_short:
      {
	short tmp;
	read_int (is, swap_bytes, tmp);
	val = tmp;
      }
      break;

    case oct_data_conv::dt_ushort:
      {
	unsigned short tmp;
	read_int (is, swap_bytes, tmp);
	val = tmp;
      }
      break;

    case oct_data_conv::dt_int:
      {
	int tmp;
	read_int (is, swap_bytes, tmp);
	val = tmp;
      }
      break;

    case oct_data_conv::dt_uint:
      {
	unsigned int tmp;
	read_int (is, swap_bytes, tmp);
	val = tmp;
      }
      break;

    case oct_data_conv::dt_long:
      {
	long tmp;
	read_int (is, swap_bytes, tmp);
	val = tmp;
      }
      break;

    case oct_data_conv::dt_ulong:
      {
	unsigned long tmp;
	read_int (is, swap_bytes, tmp);
	val = tmp;
      }
      break;

    case oct_data_conv::dt_float:
      {
	float f;

	is.read ((char *) &f, sizeof (float));

	if (do_float_conversion)
	  do_float_format_conversion (&f, 1, flt_fmt);

	val = f;
      }
      break;

    case oct_data_conv::dt_double:
      {
	is.read ((char *) &val, sizeof (double));

	if (do_float_conversion)
	  do_double_format_conversion (&val, 1, flt_fmt);
      }
      break;

    default:
      retval = false;
      (*current_liboctave_error_handler)
	("read: invalid type specification");
      break;
    }

  return retval;
}

int
Matrix::read (istream& is, int nr, int nc,
	      oct_data_conv::data_type dt, int skip,
	      oct_mach_info::float_format flt_fmt)
{
  int retval = -1;

  bool ok = true;

  int count = 0;

  double *data = 0;
  int max_size = 0;

  int final_nr = 0;
  int final_nc = 0;

  if (nr > 0)
    {
      if (nc > 0)
	{
	  resize (nr, nc, 0.0);
	  data = fortran_vec ();
	  max_size = nr * nc;
	}
      else
	{
	  resize (nr, 32, 0.0);
	  data = fortran_vec ();
	  max_size = nr * 32;
	}
    }
  else
    {
      resize (32, 1, 0.0);
      data = fortran_vec ();
      max_size = 32;
    }

  oct_mach_info::float_format native_flt_fmt
    = oct_mach_info::float_format ();

  bool do_float_conversion = (flt_fmt != native_flt_fmt);

  // XXX FIXME XXX -- byte order for Cray?

  bool swap_bytes = false;

  if (oct_mach_info::words_big_endian ())
    swap_bytes = (flt_fmt == oct_mach_info::ieee_little_endian
		 || flt_fmt == oct_mach_info::vax_g
		 || flt_fmt == oct_mach_info::vax_g);
  else
    swap_bytes = (flt_fmt == oct_mach_info::ieee_big_endian);

  for (;;)
    {
      // XXX FIXME XXX -- maybe there should be a special case for
      // skip == 0.

      if (is)
	{
	  if (nr > 0 && nc > 0 && count == max_size)
	    {
	      final_nr = nr;
	      final_nc = nc;

	      break;
	    }

	  if (skip != 0)
	    is.seekg (skip, ios::cur);

	  if (is)
	    {
	      double tmp = 0.0;

	      ok = do_read (is, dt, flt_fmt, swap_bytes,
			    do_float_conversion, tmp);

	      if (ok)
		{
		  if (is)
		    {
		      if (count == max_size)
			{
			  max_size *= 2;

			  if (nr > 0)
			    resize (nr, max_size / nr, 0.0);
			  else
			    resize (max_size, 1, 0.0);

			  data = fortran_vec ();
			}

		      data[count++] = tmp;
		    }
		  else
		    {
		      if (is.eof ())
			{
			  if (nr > 0)
			    {
			      if (count > nr)
				{
				  final_nr = nr;
				  final_nc = (count - 1) / nr + 1;
				}
			      else
				{
				  final_nr = count;
				  final_nc = 1;
				}
			    }
			  else
			    {
			      final_nr = count;
			      final_nc = 1;
			    }
			}

		      break;
		    }
		}
	      else
		break;
	    }
	  else
	    {
	      ok = false;
	      break;
	    }
	}
      else
	{
	  ok = false;
	  break;
	}
    }

  if (ok)
    {
      resize (final_nr, final_nc, 0.0);

      retval = count;
    }

  return retval;
}

template <class T>
static void
write_int (ostream& os, bool swap_bytes, T val)
{
  if (swap_bytes)
    {
      switch (sizeof (T))
	{
	case 1:
	  break;

	case 2:
	  swap_2_bytes ((char *) &val);
	  break;

	case 4:
	  swap_4_bytes ((char *) &val);
	  break;

	case 8:
	  swap_8_bytes ((char *) &val);
	  break;

	default:
	  (*current_liboctave_error_handler)
	    ("write_int: unrecognized data format!");
	}
    }

  os.write ((char *) &val, sizeof (T));
}

template void write_int (ostream&, bool, char);
template void write_int (ostream&, bool, signed char);
template void write_int (ostream&, bool, unsigned char);
template void write_int (ostream&, bool, short);
template void write_int (ostream&, bool, unsigned short);
template void write_int (ostream&, bool, int);
template void write_int (ostream&, bool, unsigned int);
template void write_int (ostream&, bool, long);
template void write_int (ostream&, bool, unsigned long);

static inline bool
do_write (ostream& os, double d, oct_data_conv::data_type dt,
	  oct_mach_info::float_format flt_fmt, bool swap_bytes,
	  bool do_float_conversion)
{
  bool retval = true;

  switch (dt)
    {
    case oct_data_conv::dt_char:
      write_int (os, swap_bytes, (char) d);
      break;

    case oct_data_conv::dt_schar:
      write_int (os, swap_bytes, (signed char) d);
      break;

    case oct_data_conv::dt_uchar:
      write_int (os, swap_bytes, (unsigned char) d);
      break;

    case oct_data_conv::dt_short:
      write_int (os, swap_bytes, (short) d);
      break;

    case oct_data_conv::dt_ushort:
      write_int (os, swap_bytes, (unsigned short) d);
      break;

    case oct_data_conv::dt_int:
      write_int (os, swap_bytes, (int) d);
      break;

    case oct_data_conv::dt_uint:
      write_int (os, swap_bytes, (unsigned int) d);
      break;

    case oct_data_conv::dt_long:
      write_int (os, swap_bytes, (long) d);
      break;

    case oct_data_conv::dt_ulong:
      write_int (os, swap_bytes, (unsigned long) d);
      break;

    case oct_data_conv::dt_float:
      {
	float f = (float) d;

	if (do_float_conversion)
	  do_float_format_conversion (&f, 1, flt_fmt);

	os.write ((char *) &f, sizeof (float));
      }
      break;

    case oct_data_conv::dt_double:
      {
	if (do_float_conversion)
	  do_double_format_conversion (&d, 1, flt_fmt);

	os.write ((char *) &d, sizeof (double));
      }
      break;

    default:
      retval = false;
      (*current_liboctave_error_handler)
	("write: invalid type specification");
      break;
    }

  return retval;
}

int
Matrix::write (ostream& os, oct_data_conv::data_type dt, int skip,
	       oct_mach_info::float_format flt_fmt)
{
  int retval = -1;

  bool ok = true;

  int count = 0;

  const double *d = data ();

  int n = length ();

  oct_mach_info::float_format native_flt_fmt
    = oct_mach_info::float_format ();

  bool do_float_conversion = (flt_fmt != native_flt_fmt);

  // XXX FIXME XXX -- byte order for Cray?

  bool swap_bytes = false;

  if (oct_mach_info::words_big_endian ())
    swap_bytes = (flt_fmt == oct_mach_info::ieee_little_endian
		 || flt_fmt == oct_mach_info::vax_g
		 || flt_fmt == oct_mach_info::vax_g);
  else
    swap_bytes = (flt_fmt == oct_mach_info::ieee_big_endian);

  for (int i = 0; i < n; i++)
    {
      if (os)
	{
	  if (skip != 0)
	    os.seekp (skip, ios::cur);

	  if (os)
	    {
	      ok = do_write (os, d[i], dt, flt_fmt, swap_bytes,
			     do_float_conversion);

	      if (os && ok)
		count++;
	      else
		break;
	    }
	  else
	    {
	      ok = false;
	      break;
	    }
	}
      else
	{
	  ok = false;
	  break;
	}
    }

  if (ok)
    retval = count;

  return retval;
}



Matrix
Givens (double x, double y)
{
  double cc, s, temp_r;

  F77_FCN (dlartg, DLARTG) (x, y, cc, s, temp_r);

  Matrix g (2, 2);

  g.elem (0, 0) = cc;
  g.elem (1, 1) = cc;
  g.elem (0, 1) = s;
  g.elem (1, 0) = -s;

  return g;
}

Matrix
Sylvester (const Matrix& a, const Matrix& b, const Matrix& c)
{
  Matrix retval;

  // XXX FIXME XXX -- need to check that a, b, and c are all the same
  // size.

  // Compute Schur decompositions.

  SCHUR as (a, "U");
  SCHUR bs (b, "U");
  
  // Transform c to new coordinates.

  Matrix ua = as.unitary_matrix ();
  Matrix sch_a = as.schur_matrix ();

  Matrix ub = bs.unitary_matrix ();
  Matrix sch_b = bs.schur_matrix ();
  
  Matrix cx = ua.transpose () * c * ub;
  
  // Solve the sylvester equation, back-transform, and return the
  // solution.

  int a_nr = a.rows ();
  int b_nr = b.rows ();

  double scale;
  int info;

  double *pa = sch_a.fortran_vec ();
  double *pb = sch_b.fortran_vec ();
  double *px = cx.fortran_vec ();

  F77_XFCN (dtrsyl, DTRSYL, ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb,
			     b_nr, px, a_nr, scale, info, 1L, 1L));


  if (f77_exception_encountered)
    (*current_liboctave_error_handler) ("unrecoverable error in dtrsyl");
  else
    {
      // XXX FIXME XXX -- check info?
  
      retval = -ua*cx*ub.transpose ();
    }

  return retval;
}

ComplexColumnVector
Qzval (const Matrix& a, const Matrix& b)
{
  ComplexColumnVector retval;

  int a_nr = a.rows();
  int a_nc = a.cols();

  int b_nr = b.rows();
  int b_nc = b.cols();

  if (a_nr == a_nc)
    {
      if (a_nr == b_nr && a_nc == b_nc)
	{
	  if (a_nr != 0)
	    {
	      Matrix jnk (a_nr, a_nr, 0.0);
	      double *pjnk = jnk.fortran_vec ();

	      ColumnVector alfr (a_nr);
	      double *palfr = alfr.fortran_vec ();

	      ColumnVector alfi (a_nr);
	      double *palfi = alfi.fortran_vec ();

	      ColumnVector beta (a_nr);
	      double *pbeta = beta.fortran_vec ();

	      Matrix atmp = a;
	      double *pa = atmp.fortran_vec ();

	      Matrix btmp = b;
	      double *pb = btmp.fortran_vec ();

	      long matz = 0;
	      int info;

	      // XXX FIXME ??? XXX
	      double eps = DBL_EPSILON;

	      F77_FCN (qzhes, QZHES) (a_nr, a_nr, pa, pb, matz, pjnk);

	      F77_FCN (qzit, QZIT) (a_nr, a_nr, pa, pb, eps, matz, pjnk, info);

	      if (! info)
		{
		  F77_FCN (qzval, QZVAL) (a_nr, a_nr, pa, pb, palfr,
					  palfi, pbeta, matz, pjnk);

		  // Count and extract finite generalized eigenvalues.

		  int cnt = 0;

		  for (int i = 0; i < a_nr; i++)
		    if (beta(i) != 0)
		      cnt++;

		  ComplexColumnVector cx (cnt);

		  cnt = 0;

		  for (int i = 0; i < a_nr; i++)
		    {
		      if (beta(i) != 0)
			{
			  // Finite generalized eigenvalue.

			  cx(cnt++) = Complex (alfr(i), alfi(i)) / beta(i);
			}
		    }

		  retval = cx;
		}
	      else
		(*current_liboctave_error_handler)
		  ("qzval: trouble in qzit, info = %d", info);
	    }
	}
      else
	gripe_nonconformant ("qzval", a_nr, a_nc, b_nr, b_nc);
    }
  else
    (*current_liboctave_error_handler) ("qzval: square matrices required");

  return retval;
}

/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/
----------
X-Sun-Data-Type: default
X-Sun-Data-Description: default
X-Sun-Data-Name: RCS.tar.gz.uu
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 446

begin 600 RCS.tar.gz
M'XL(  L8Q30  ^P=^5/;1K>_6G_%ELZDMK&%9:XD#FD( at 91O"&& I$TS_3RR
MM,8B.EQ)YDB__N_?>[NKU6G9!G.D,3.)I3W?_?;MI>.=DY4?[OBOU5IK;;9:
M\-O6UM?;\-MJ:>UU]BO^?FAMK&UH;:W=:F,Y38 dot W'];O&C#\&P6A[D dot 7 at >&%
M(1U?CI;DX1]BT5Y;Q5] 8G6^0-[=WS'P?^>='OK6E6H8C8N[Z*.EH02,X7^[
MO=[:7(WYKZTA_U?AY8?670"3_?O dot ^3^ at NEG1U'9'T0V#!D%'":Z=GF?#@^T9
M7X(."4 VC+"C&)[C4#>LO%I9(:\ZBJ) +<740UIY]E1MP4M+;6NJMJZNMSL5
M?10./)]PHG8J0..0DMVK84?I^;IK#"BT[]*K$+K6.MB45M"4IJYMSM 4PF32
MP%!>O1+0V=Z9\JI>KQ/J#,-K J_$ 1SU,TH at 57FEA%!/082X!A!'=ZWAR-9#
MRW,#55FI*\J dot -[SVK;-!2*H[-:(]>[9!_N,-7/*;2G;UT',5Y71 at !:1OV93 
M[U#W0^+UR7LCU"^HJBC\ ;/Z/J4D\/KAI>[3#KGV1L307>)3TT(:]T: EQ42
MW357 %_',ZW^-20H(]>D/ at D'E(34=P)L'5_>'GX@;ZE+?=TF1Z dot >;1GDP#*H
M&U"B QR8$ at RH27K76%S9P]Y/1 dot ]DSX-6&9H=0BW(]\D%]0-X)^T&@>ZK>H@0
M^L0;8JD: at '6M &'B@BG<8A1,8KD,OH$W!(@'T X@=6G9-NE1,@IH?V0W")0D
MO^V?_OK^PZFR??B)_+9]?+Q]>/JI R6!W9!++RAOQW*&M at 7- MS ;V CX/]N
M]WCG5RB__7K_8/_T$P"L[.V?'NZ>G)"]]\=DFQQM'Y_N[WPXV#XF1Q^.C]Z?
M[*J$G% Z at 7)*GU$>"&324+=LD 'E$_ I )!LDPP06Y\:U+H @'1B@&Q,9H>B
MVYY[QA 38 at %*)4!A8K/S_NC3_N%;@'"_3UPO;)!+WP)A"#U6!CFG%'"N0=:?
MD5,0;6CBR-8-2IKD9(055U=;#?+:"T(L]&Z;$#"\FM;45EN;#?+A9!NPJJ\H
MRD]6'_#L6RX at 4^UV 8.WW6Y-^6GHZV>.S@A/4>59;\I/U 6)9+6@$OEU^^-N
M=^?]X=[^V^ZOD dot  at :]LBDY(7AN7WK3!V\3%20>7W;T\.7R20+H/2I[F %!13Q
M]]]_)WO[O[_;Q:=47R>?3KJGGXYV3U+=!=?!2G at ]I$%ACTL[@,/5]N[1Z^T#
M=;"4S7BS>UJ0>K+SZX?CHO2/;]*I_<W-9M]PTXFVUZ2^[_FY5(M2FDMT].$0
MU"F7/ at I1^%*ISE6SIP<TEVBY-G P .^=S/",L&D@U%@>";OG^2'H#^F/7(.9
M-W))P?[8- at @#6$'JNX#FDO*W0D"!0[*WN=G=VSDDU:]GU'$:Y(^WN^_>U4@5
MV!N$Q!CH?KU!"EZ at [I dot &4JE4"/M+)!8\[W@H8E?Y"B*CO-59"D&7A<5)7!Z5
MM,'^KX$GR5/!\! at 5=MX#%8K;2C[#__4$,*8')H$FH"CNP[18'V_V;])'7*&H
M:G&' at <TZ/#FX78>3NK&#@/5S<'(B96B<4$@V5HH!DJFL3#*'$[G>B*D]J2U9
M(P8=5 dot 70"X7G0O/K at ?>W7+#J!F^"]/LA\<%%H=:!$H%-!IO-3'8(+@L\"&^%
MMRUKZ?[9"(UI &;^=$"ON3/I4:C"7+W%7;5NFA8X"]WEC3#G9X #%:T-P?M8
MS$U7]>: at X36_U@@;$6'31$#<HP"PB\UX?=X,U8T!"48] ;::XY$!2('L[>SM
MG>YG dot %0BL5BKSVOMC:DUL7Z/UW]]D_I?;1AM at =+^ 7[^]*UL(=;Y C-3R:IB
MXJFHB] /KE%)3H]//AW,8 at  K9:(]P8*5Y9.)!2(,64&>E#5N0IR2E'3/**/D
MX=O=8C1CZ H->UHO,T6$HM4Z1/E'45!B87R4JA?]= at _UPZY/@Y$-(VZ/#9<P
M"0:E\AE1 *&.*HJQNV'K08#C=9;*$Y\_3[U&:/&W)T2O 2F>PP at )1I;7[1>B
M[$O0+=7W+ at -2K36(KAH0#<%CC?E&'"!644+ dot R19I=>#G!1$%X&5Y&5M,E+)X
M*0M*B1;A)2I%"(4!%A1LD/,:%-35Q'L'"34-,L?>Y4=JA)[_A/ at 7XQ#2&I"I
MVM0] at T$HHM526QF$DJ FBB;@Y=#!T-)":*&0 '=Z6'<\>^2X$;C&6'"-%*S:
M!'"-,G"!F"T$U[ at !N&\L_>Q&LC(!8+T<7BLC##,1F"5^6S+!$K\QT1")CTM"
MLO$3:3;! dot &'D>LE";PA90XP7=0S3'>N*FDV,G7 at <&P4&V,H 1C6_3$<(=!$)
M"HQ'2%I)*V<E$[94TBIA32?;R9[GV5EH/0BL=! FLK6585H,+G=!/.CIDZKH
MFORX%<-!_O>_R,#S=.D.$#"?AB.(G/JZ'5#F6$4"_6L$(\8J!.JZX+M\E(R=
M O0?IP)=]/ECM1[B5!3 at JTMY,"G$UR-@+,[3N 'UPQ43J!?2%9]Z/DXL\<XP
M(LPP_$D6+%X_YT39Z)GX_,>HR1!2[[H^25"R$Z<;)$%)/NA"\ at /_6TAOGRSS
MVB])D at M1ML&S#<A dot ,>-O(2W5NC'R?1@0=VVK)X8-+";O#G37M*D/XYLE'P<\
MA"4S0>3(+3$H$YQE).5I_RAC_#_",L'Y(S8%GM]?!ADVEHND6LF"D#5,$]B3
M< $E' )A9%W'UB;/BY=;)5S !AZ2#2D: S Y:^DC at =/V<@[437NK6Q(XIF,Q
MG8$!CYG 3(;G3>"48_V^+ Q.E5=A=-! C(U&#&L3ASL2-'CCB,TZ?A!&IV at ,
M<7NU2'NI(L9EAD=Q$SI'F-%G'E L+ dot #=6L""8?O"$,[=$!8%& at M[^&CM(8=6
M3*]=Z';,',::#&,86R13 dot $\8$UKDR1/"Z-T:'RNE1IW%>,M19R4=0 at %D-\<N
M/:OZ;\4RYJ%0-TWHF_CUV^*]/3OVOB85KRUU3*89+$W at BV4TM'S0-B\OGPV>
M;H at Z['D dot >HGH3]1*@068BW:-_,T0#9TA8J]U$&!X:'<0VBU,[T E7@=A1I*E
MZAA0Q\ Z!M0Q$G6R 1<6!)G at !8L#+NS?PB+8_]BIUCG+_T)"OE,)R0H(+N>[
M9GZ=(3G1-)N53,Y%W81UWB4Q+8>Z;-'0L0)'#XT!(PB'=3(?&6Q=,2I!"<+,
M]%P at 5+W ^2X7/+3+/+*<)A.C><A6HX$-ZZ%!6CC]6I"MLRS993(>X"5G94,R
M$+ at =)[1OB0/)B<8'YT$F3KBM0L2X?3,,T1X#&]+K%0NK]* :4;(8L>#&PW!C
MX2H>#2L6'N-1>(S"A>Z%J;I7C at 2A;GR97U!A%"QD3\\%I5)=,IANCN,$ W<Z
M1O at )1OAEC&!,D-LH7 dot -FC)#])8K<F!%S\17&+:W3@W%#>R1<F)>;,&[DM!]0
M&1*[H!X%(^;C(1:V:4Y:,:? at 8L&/>?)CX3(>$3,6GN-Q\6/A0![68 VH[UBA
MI;ND>N%9YDTXD$403W%$G)%2%R_W_BV)B 55^+&^4N0A(B-)7+HH7"E?#Y9M
M\\6=<[%U&C [)]7$ at D\MXEV":ASX253#0Y7!T OH?*F6HL($\MT!@5);4:>E
M#*?JF '(I*U 1;H1X1[M"FK)3=>H.%N9.MCS>=<<#7''?;SS&NNR'=B1 at D4'
M&'-!?4995E:(D$ at K8 dot ?JDGNIZ5\C"PKBJ3LD,-MX[;ET@K#0*Q 7(^3\*%DN
MSNY+OX]E5":9]+++!+9-FM at X^ dot Q dot (HO)+689(FN\\&)##5XILRLG*:28/6'G
M-&LIM4Z;$MAHN59LY='87AYM.;V#.B&Q>!A2, 'GDJ!#8>XM5?).CM*R[,,:
M'+Z<FA>X(BM:K;<2&PAO, dot =EN2!HELG at #0!/=DHCYX&RL M A#_*Y4KGPYW-
M=,?:A,\1IB)G)\I<3CE-\1 )J0=9N?^1! ]$+P80\K0>1.PT\&#%S_V?^7XX
M?-[[ dot 74(A&'"/2^UP1_$E>Q$I8 dot "2M&6NZ:8\L7J#X!WS*WDT#C+,*$O8_0 at 
MZ^YR>G#S_8,1? at * B2@FD2C4AE2!>#0V5B'0:T[4A?/LSL&\+DQ#W3O1B-M3
M[B9Z$:$TJVI(4HC35_/1CEO28-)P$+JA_KC!H.7V/6Q(;"WT(==,C#UD72S7
MX-FUR2/01+7P">LCV>]=]R6/F_,V^! &M,> at 22B*1W=3CY&CX(ZM'T5;V*9E
M.P/8QP$;C -(\-<(K]%Q&"1+Q0*%>#&U5T0".[/P N!Z2:SA12BM!"\<DOJ0
M)6^Q7+7/KUKI7E C at CYN)#KW0+ZF6A')T-)7: dot ;KN#;DT#<1<"8JPVBNRT:^
M6Y%]&M, dot 'OS_/7>S28-4HQ8P]L!0M$$8:HV(N\ dot OM; at 5Y$D?6J)7!F771'6I
M:W at C-Z0^-3&6GH9#(]<'"(!+ dot LHIWV1HN83!%2LGXU'E;PQ8+CR\H O*)D6[
M dot [1'01=&X( at Z)N  544.5L08.E-J"[-9 "3YW=3BXCP)TW"/\H]"HGGQ#!/8
M#1O'$!O at O5'02\AN^<"B'&@N694,J\R1XUP+(:L4,$5>!5/&%-:(O' !_X9?
M&V F:W&CDW at T!9,4;+V,4::%C*IP"UGY)W87,\Y # ,Z,KVN-#*"P:$G-IB/
MM2(BX^3C&QF!,-YP48TGEL at )TPP6/P26>S:R=;\+C8QH<2S^(2YNTW[8E74<
M$7(6U/D8UV'WQA568K42#B:P\,:M+7*BF at "LM'D8%4(JRXW#9?2$G5EM)Q 1
M][RR4^<I6Q=MQD>S6L%2;/-P770JSO?7(.'-ZX/N[M')_L'[P[16BDK&5)68
M7%P.4'^K[&05.P20K =6X07G.5.V9C-]7*5\#H#/FO#3]86V/<- at )-\'5<;E
M+;Y8[S?Q< at 1I?!("] ;*)UZK O!$?;]6(VKL7[,6/I(2[/ACKF,CW;'$\B-R
MY W\^^"KB7FZY-!L1EWK>R/?HGY^H#('5XV)PS  ? (=&PN2'APM+X[Y7#;H
MR\R;8S4N at "B0Y!=\>B[V:+!\T2"4T>*)J1R393,%]:+&XDD/G(U:JV.=96V=
M at 9IQU9<!WM54==VDNH dot KYNE;/+_(U>;=]6RN>MS]3)RZO/^2Z96(^G'4,N[J
M)M[ at DPBHS_A>/_\SW<>,$F8M1.R[%['$[5Y3B5AA'UBPH"/9TOF??")3OJU$
M0X<:JWP3X16RVS87TOO]2N^L!C*B*B- at DJ*,S 5DY'21ZUZ<?).IEJ$^FT/$
M9C/4QW0FC/="N=QJ5'(2'[50K# at A6)^ME,I:=;8&>OYG*CX=SP9L(@WC-!WG
MNX-R IJ;CZ.LA:586(I'Y.<6-NA.;%".P?.W07EV3K9);W9/LP;)I'A-O.7J
M ,5-9J%3]<?-#D_J=Y:IZ/EVF)F/+C#'T-1M9J"WMOAJ%C>ZF=F42/#-S^T_
M(X-I?FXAE]EL9)2B84HK3LE.9R",5;-68I<?[QQUVNOI?.?!Q/EJ+'?;V6K_
MNYBMEG/)B6GK2HD8):9F2Z>DN<S..AD=TYS1NT&TUKU./YNIZ><(P<I8A;K]
M+'7 at V3B$29_&Z<UH944CO0:9;0&NN'-^@?FT!O? dot ^IYL>X510)M0[=7R,$'&
MK:#*;#M[$*YD89 at 7<^8"P=0LFI^'9(]BY39^[J5/ 4YA!KA^5Y?$,DK!_ER+
MD6^$9 at 8_RX+?Q=!]=A<KN]_TFUOM77C2>_:D[0FN4GJ57B>WH"L6'L?'H+*O
MD/3XA6R]Q(5L(I<%$!+33,33BR^:K2!-*AF&RJ]YE#"4+19C9T\XP)]=GT6G
MPFM7IG#9 at B0"XJD<>*7,?P?L,B'98,^G^I>.> 5.5/Z9R6F7[:XJLI+I8R#S
M]QBWA&?>KGU^X#24S-=MQKJ4 at HUV=^E8,B>C%J[ENW(MRH1HX=_M76;V!V*O
MBG0 T[%A'F$;-_NWC\?LX$8A -M<H[M?$M93MA3;3R at PQ?A_ @AC;/C#0" >
M6;O3C_X'06:\PC.< N-=8+N=V%[+)^>NHH ;V dot GY6T- dot -?1M#JZ_D%_ at ]SEQ
MBS<#R5,S/E/-05 PH9S:ASX(RDYW.6,.=T5G97JY+PVD'4?"UK#:)2B"J^T&
M#,D7>22YG^+V]"4)V!FV;B"-KS"T]2$V$(SK)6F.V627VDI#<.GY7SH)[\% 
M >19!M1HUQVPVE4FPIP3[!' at S+B!N(*;J dot #$%9RQ+IC5K+(&"ARQ:!=_2FGI
MBX+5=5(7>$0DK>&FO[6H:5E2/+PD&A04+\^CY;PL WP!I9\"4W(A:M,O at S/E
M7A)?]0,'XX#H<NEMD-C9 dot +A2PCT,WZ;*9'P82 dot ^/= at C74:#3!F>F*,@!XF;P
M_F)/P*EP>"#<KEPFDWI:F:RC*3U dot + at S%>AJ[_XRNYA6X<_LPJ-A93!D(W<QK
MW1J<N?K0>4+#Q+7,I99'/P.Q #\WAWH7T4_*I5(=R,*/402X= at DJY#PROYHA
MN?2NY0NUTF^FE2[E,A?^<N$O%_YR=G^9\93CO&#:!^8=X$3OQ[X+" at 98=\ dot  
M=2>,&KT:0H#NAA9^R%>WC9'-  at )5?HPT^L"N;@('/_])MI@)7U=;F3_:U!J0
MH:G:1NIO4V9LY#+:(J/57M]8TZ+_-VASE6>T-IZN;F[(_]=H<PTSUM2GZQOM
M]37Y_RIM;O :JT\WUUM:_#]M;O*,9^W-UM/U]D9K37NZ_FSU*6T^:RC_=";>
M2S!T9MB^ELYP8IN>F;\#&RR]%C &CSI0$ at RL?DC8H9NX*<CCZ6)7!*]PA%+G
MFA9R"C]G'(1T2+3GHB'7\QW=MKY&C!PC7D;B\H"HEV6PP+DO?429*WSW%/.T
M59$( at JS;U5I\_B0%+YB_J)SEZ&=X9IX58E^!=G\ dot R27N40%!8]ICG5%7G#D*
M/=(#%!U*1B[/G at J-% dot RDN17U7DJW]O-$UZ2GV[IK0*::8 dot CV[M'K[8/MPYU=
MXO3P,#18HZ777+^1RYBH\IK4+#L#9:8+0S>9PT_C at %Q]3@)0SQBNA!_GYEA<
M6"#4%8Q=%^6 D66KY'O&2_M+_.0<_G,RYIJE,BM:8,OCR8[ at KR';S8:\J;&M
M2JQOO)A$;463X;^0JL8DPO; dot XC(ULL(3VFJK%JVA at  dot OAQY0X078&U/A"+@<T
M'% ?O]CI4FJBB( at CJWI(=-M6HR$@AT8*9 1<JY,ND+LM1E .R4R[?9V1-K4[
MJ4C\6%N1H4Y5K<, (#'$8&("R";+)#:R 9KN<- at @YG#XG-G:G_'>5M^[(D//
MOG8]!^TS$Q6#!FK>UD!E<<;+B,]X961OB(-,*!C1]1"H -!3/R#Z?Y\2557A
M5Y/?0'<LE\_>=L_%9II,1+>)$=U+'MDUF]E-MD.&<)T]+;,GYD(^G\?[OV09
M4Y:I)GN-JTA/FLJ.H(II^ ?U/8X30X/WT#01YRD^HI, /K[V!Z+-Y?0FM9+<
M" ZD_"CD7E,P$ at P@BQ2VXB-PT%*-4TA-[>3%'D2L!7GHVA5QEI97'(ZQ$#X=
M4O8)6Z88D:T09PJ9I&:OEQ(=BH=ZXGHAJ3IXS#")6BD08$LC[!*&E-7+?8 W
M&+ O\&* at ! ;^&O2/WQ9T:4&4=J /\2(T4!8[X/4M< 14-S'2TL&'HF2'T7%U
MT6GD[C)XJ8FKJ+B-1!KG4O/4,%/)9L03GE)4)]/F% 0#I]W'?3]C''=J' <0
MP7B$1*ZW%HWKQ,T-%]P;]*[9GF/QUGP9D6CLYVP5^57=>O'=L1>9+]FG+FL<
M&]"SK4P7R:U,F%(G>L&<?1$$!;,+> at 1( H)>"H+X]I;\3BH=6L=;G6[<_=1T
M*%XKX+=[760^]#?NXB]%WF[&[NUCEW4Q/3SSK2$%WPG!" at 0T(#"XRW<I1F&)
M?4F9?0%<BV[Y*EZ#C9O'E;FBZ(0U%#51282B1OF"6RJZ<QP6W+U[!['=TB& 
MQ_Z+F^: at JBUY7/\B]4WHU#>B^6?-94DCPO4 _]5P>!DM0):%=HC]=.NEK&1)
ME.>P"9Z)D1>>7??<R(-?H8[B.$"?43^7QWP$-B653_ FG#$["XL_N!ZDM61Y
MDI8T[P:,YJQP+).2>S(;<CQW7]28*Q3S)L9#B<B=0#.>."LKD6*!BF6U[F:Z
M%D,9P9V^C_46A,1RR^4.,4O(NP*C60[&LKQS),AP[YZ)<1]03*;%.)8\#$WN
M$YJ(-J!HL2^[ at :+EOB 9DW>K^.+\&6Z.E>.JHJ+BF[_%5VBQ*GPC'7_+?!5R
MPN!K>6LIOHH(VVKP-J* RI_A4\HE'^;ETUO+<_ at J;RQ&WS#5F_=)]>8\J1[+
M^GA=?>3$_S>(_#=+_,<N^05.,R*P,]E+2>([V>^V,XHZN4UW,W_F?3X:,!4/
MTBXW=4%5V?G9XMJ\NUK!FE.TIN_4.C,Q.KN=9H*R_;^]+_]NVU8:O;]:?P7B
M>UZN)5N*Y&RMEZ1IEM;O.&X;.[E=OCX=BJ1MVA*ID)*7].O__F8&"P$2I*C%
M6VKU-)9(;# dot 8#0/,H##5>=XN^?HFO/F53G at 9@T^>\,8,$S[A;.V54H*<^[BR
M6ZU at SJV dot -7M0!A\-][=E#P*()M/=&RY_"P\ dot I$ dot 8PC&7EN-'CL2)?=O dot B$$V
MZ7D][? ^4)EW 15PHP:S:( at 7HFF\'A#>;^HGYQUCXT0%!8B0K?(+&99DVSP_
MQ2!S]$\&0VF)E'.]M6?NC?;>,WW,UI(YJ$8A''\K3R:1]I804^8<&#U3<DK:
M#J?^T=7>5;OC2TM98' at GL[@NYW*GZ%(@*R'NICJXE?K?*50'@RGT?_XLR/3Z
M_^N;\%NI_YLW/N.-VSGC&3U?<8Y+]/T\<YPQ B9-^%IZ#L9\[QH&P%3SKELI
MYGDII8__*K0)M#AK4FN90_F.]3R^II)S';;GZY"K,KV;.1HS!]<H at $:?W1<L
M*SHF]6O,$Q)6_MXI*XQ9)>W.NA>D&/)NRMY[93OM=ML=G_![95OII,HMF?&[
MJVSUM?>]NKU7MXM6MW=ZQ7.O=6?6NG=ZWN^5[^S*]W[!>Z^#[W7PPG6P<?!H
ML>>-%K#W?]UZ>;;]_PI2 dot 8U;5%M&&#61V5/"\(:W83+& at /<H[%_B3<_0PJ$?
M^Z&+]Y dot QP24;QL&9,_*3!]B:XWGK;,4S#G6G at ?G:8;,YCG/<O7F<\1S'#<YC
M,N[1!5U7-)FY at U%W;T[O>;/DJ-7=F\Y_#HM"1^/0B2\-92I,X*()?I /DA?M
M*OLN G1:AK FYWQ-W8"=U_-S!<LH6BN*.[ ;I, ->(1 #5D[#X'GHBD&R79H
MNC[5J: at I!R7G=/$CRQW?F7)D U@S!4.@RX6/[-&<(_."L\#S%S^NK-M'&YX9
M]',3]#7;F*Z!O&8;V#50UVP#FYNXC(BE&4_6% 3%9)3L;&28ZED!2$4RG&=,
M*1DFAN+2QC8[&<XSL"LEPWD&)LEPT? at JB76J3EV+'E1)R-/-D]=<([L^,3?M
MR.8FL'D#QB;KW4%)8.0=\;??RGV6VQ-*\/5-^*W<8+D]H0133/A]*($1&G"5
M at 03W<01W,XY $O0-!1/,JLXMFZQW2\C?:_69M/J=G_=[Y3Z3<B^>]WL=/X>.
MO]?P]QK^ZL,%9SF;,84M<#<5P0T; "6NN8'NRU%?![I;AZFFYM#E=W/>;EB!
M3_)>+FCR)J;FN)N3-S?372/;?%V8_V>PS40[MM1XY3G7] ,8 at RG2TU:D%TN7
M [W+&1+RSM3SO:5>;*DC+'V/D4&VF3YP3) at 7E@"D0F;>]"0V_J\GYM4XIY\Y
M4N0YN<R\V(Y*S'L+\_*:%#^,(V^,XN$K,F0$3+=1'*<[G?.)X\_C:!3XX=<U
M;Q*HVSAQ<B-XOFG+<-M78?_\ [GMJYBWKY';'CUB$=V64^QT,7YM; P<L [=
MKMM][PRA#CLLOO"J!TH:3POKIF2O!=60P [Y $HZ\0HZ,2B(GX$6D!*"#%N2
MK>BGAL5U6KR<O(<N+C&$[)XXJ_N;^\! at :CWXEG%UVTV)W#%X@1D#N3I:%W3>
MNRI(7C%(Z2'P7A3U<W"$EUUT0:)Q%B1=O,()S++0"?.'ORL?Z)]X)T^IRY2S
M)MC1$I/\,A;33XI\>1$D,%P8I] dot O(W/";SYN^"WN at ^8KHGCL2X[7S<9#!RQV
MR5D?TK+(]&'$!%82= at SF*L dot KS@(Z0>_"L.AV/;RRR8K2?E^B-.DZL=_%Z]1N
M&)V(+P2!K6AXK-,U(2VQOM%04LO at SD 1B-HO>",3- G$(]&$5_B\?@^X2IB#
M19HP'O_(QPM) ? at HQ@92' Z=>-1B[%4_ 3Q?\#4\B#;6=^(C'P0&K I8,@!$
MXH\56Q-U)FZ5PZ)RL,< &JP5WK_ZM?OIU2Z]>K^SA]];K'BNQ$ at 3>5H35)]S
MT87FY6GAAWA#59?(;(;I4S&6;?;P(;[.75)F(_0V++UX?6)O<<]I5UY,1 0E
M[BO2" at 2B )]J7D"4$#!A[:YV'92 +'VL$PR]>2&K IFDC00%I;?TTJKM(.WR
M;WT1GB.[JQ<=M:429"Z5X9&$3IS#B06UO& at >?7K10"\:2]P).J.BN9E9,M!L
M+6JT^J:[M[-W(%HG9N?] at *R4KP+U*E"MF9 dot 2EYM*)M@8:A1%7>+B+LQ4][ ?
M.:-;K4D61 [O= at ^Z(':X&P/P&Q0\EQ--SW?V,N75\^G4E]TN ]$V*^8SAMEF
M)2F6FLJ= at ING\$HH[MF*]9MO2=0Q<>%>T0WSKGT7^E JM+9UE]?:C]C-9:P'
M<WUJ;+KF(_VYQ;\HB*Q;Q 40!;R114)4"(3T.$Z!^NQ 3G2 )[*F!,< UYA 
MT5E9;PIL_"C0!=PY%/QMM?(+> <,FSO%.^WY>>=!%4KK7!OOY"":GG<6"M&5
M\TY[7MYY, WO=!; .Z6^!W<\0.?9[&QDWVY:&#=E-UEHL!GK>P&SBLT:AZ!X
M ZQ)8UNB7AM:MR>K at BEF82G:G9D dot K,E<98 EE%(>+''@R I6L-I1'5?D*_UP
MVQ134 3U27TZULK"K$_ED at U@;-D"= dot 9,UUS<E(P'=X:9<*R+YR5HM825L,_5
M*^2D"3#-R$A9F$RR,F"Z?C;* at +P0+N( %S)1'N*%\="M5D>+61A9JC>R$WA3
M2YRRL06"L&_K8L52OY%U2<]+GK=9OB]FZ6&IOGJUU%E]$5$VMANASBF6 Y;Z
MJU= G<GGKYP^<_JOR(:T$PL8?\S+HOH&E\-Y* +-?I at !BMM$\]S"L(+(#09K
M%QI\TUD4KVE#^I/OCJ(XRQJ8^* P519_V:Y/UQC"?9KC,RNC93D-*>TT92-\
MW]QFI_+LHWB_E;ZGQ'RG>G29/C[$E>13C5&MG(K#H?%O\\);6*[ dot 7E(G&_ at C
MW=!*352L43=VDL3P,QE0=3K .CK->]+(K.N$$*R>6E4*$_#/V4'SE ?*E>9M
MG6'4^B%';:_,]6- Z19;QF8V at +X^C_UDY'MIXHUH3-NPL1,>^?\3+NO[N%[A
MA at W0$V[\X[9'E\XD\&'KQ(?5Y!F-;;$?D9'SXK1&M4T;7;27;4IKT@;/-O#]
M#LU35W at 25(B0(I#Y>"U0GU2%VF#$V@1A==T@5Y R..F#P'+&XU4<.Y=; ,D+
M -3S+[330*H dot O: at LT%2UM dot F'O&W+821#\LBPTVDV#"?;$R+V-5TAJSUT&I3Y
MHE9!%S'QH:E+PR)%;*744*/!D-" at ,[IW4>>%!*E)6B0PL^0IBHH=2:>7B.;2
M,B_%WJ7HJ at XB%XJEO]/SVO*DC'PCY!$?>+-#%,>%F:&S#<;N9-5[1B7S,-#,
M7JL) "]B!4 ?O*R9&7D][^G6 "CR;.L!HF((6Q*;10V>I.VETTAA'?)Q.AVB
M4=.9GAH]>MPU2OT5[ &5$0 O,-?=<_:Z(K)[0W:8^HB(2+46TB#;OQ4IVJ1$
M-H1[(L\Z%] dot +!ZPS at WC :O?B =&P2/% S16(!^?"% _.19%X@#?_;/'P0F*S
MJG at @O dot ?% SV^$O'@7"Q>/'R(SNT,*XRG::T'O5I>0E3JK8J$4 U=DWAPB\2#
M:Q$/1?9P1?& at 6P_>A>09BWBPVK<GM\F "+B$*%W3WQT),=& "*[$@#B9VX X
MN5(),:4!H5>;04)4M2&^9 at FA&Q"+D1 W9D-\71)BH@V1E1"+L2&FD! %-L1L
M$ at (/ZC^*:K4H&0%R!@_3T/ dot M+;8BG[*H-"TR-$+D?ACX?:]['GBC8QA+E+2&
ML>\&21"%F!1GE3W?+'# I)%F>3^,/>61Y'1ALD8)^=O at OT<-EOBC<[:B#:; dot 
M&H_P?>9Z)=Z^J$I> dot (DN at :>(!^,$.=R\> &-2-P$@)L)@7>Y2+KB.+HMUA'A
M:UNX 4&SG+3<ON\ at &J)D8Z/G>+U@5+>&R&O\MIG%G]4[/VD?0+![D"#,G)3I
M-SK'$A4+[F1.#J0%4\'!V%$TBH#C0^XG4V2)3S8TZK3=;5+[(3CST6^9R:)\
MD;N+X)*C7 at @6U]6VM)@+,S7R@8=B>%KCT?@\&+_OQ* dot C-?;[[JL/!S_40;*L
ML4MHVEW3ZM0M:=V.V,KZ&EOGKXXR&T&\<_F4=HUR3]OR:9(I2RTT ; at 3 #JI
MZQ[B(PM^]B_[9^ACCHOSZ5N?]]9H=JSOW G))(#M?_WU5_9NY]?W;^E;L\E"
MW_<P.LX]]MU3-CIV1MAU;XTB?%SFQ"#L^WV*($J< at <];037;$BUB;^ dot 1S_;=
MXW',/)]BN)) Q%"JT>R__O'C!^8 14#SRQ^7]:U+_ at [U12]]QUL_B)TPP=A3
M& N,,O3/ ?0H]H(0X_A:^0D>4R:(I#4 dot  at Y$37W8'>N2I631QC[NB=(*#U\OF
MV^U!R5[E=D7I7+L$EUG<1<4S=EK'?CP M#DD>QL ;P,ZE?.V'P&Y\%E0A.-_
M'E.H*DR7XYXV1Q)5:YG +3%G8)= at Z=:DH dot $>?]YKZ5YXF6S?=?J^+!B$AY&.
M*=88(CH)K;D84KU43Y3JE99"M+ at 7EB(U(RW'* :,@"0X^+#_VZZ1F*,C Y-[
M]._0D;^'/<J[P4T,_N1"OB,(T]<(HY&6 at TW*RH&U)B?E*,G&00 MV]6%C8,Y
MY^) 7TJC5^UR-($=&DAA2$L&A6F!V/F=V4>-VN;F)MN- !GLDQ,' dot ,1D at S4:
M_,4 at \OP-]GIU53UY&WK\=>-1[;M:K=9I=6K]Z*CVW4Z(/:*)?$:&!;P=^1>C
MVG=>Y^FS9^P)_+X___<5G/_[\'K_D<=_MD -G_WK"C[M3KO]I-W^5[N]WGGZ
M=!W^MMN=]:?TM]UY#B\?_ZO][,FSSGIGO;V dot Y3J/GSUY\J_V50PF^QDGH!> at 
MR\2-1B dot _N)Q?\ at X_",LZ @)_VNN/%SO(J_L<^XZWU&FM;]8<U_43L(Z2RT$/
MS.5-D /N:;+)P/P dot 7%C$ at (& 4<U+WX$P^VZ39,5ZS0-UOO3M-ZTV_&BWUD%\
M/&T]_69SR1F/CH$?.%(WEP#'8&J\O0!KM0=* 40?M!^"/(&N8=U(8B??5*?U
M]/D43>&8/#]Q:]]])T9'D at RD&P dot [<G3)X"<; (S.D4\R3XHT $CH]($3!L-Q
MGV>2:*$\!8X97L;!T?&(K;RNL\ZWWSYC_S<Z#ME_6^PMK%'"6NW at &*SVPP"#
M(A,*X\83!S^1] at "ES;_@J\/8!TT<'8[ dot P3K;9)?1F+D dot \JJ'2YR@A]980('>
MCP!>D-7!X24\J(UA\1F3#0&J:I! at Z_CCA[V/[ <_!#W49S^#E@]<MANX8+S[
M:*X-\4ER#$9B[Y*,B7?8^[[HG;T#Q><1F)O,#RB%!F at T6D&"G8V2# Q*&"%F
MUL!2=0QLKP%BTH(&;"D('BI#'-\QK.*X80I G0= at C/9\-D[\PS&H?2C)_KMS
M\.-/'P]JK_9^8_]]]>'#J[V#WS:A)$PWO/5A'4+M!""Q F at 6Q@WS#=,(\+]_
M^^'UCU#^U?<[NSL'OV$D_+N= at [VW^_OLW4\?V"OV,ZPO=EY_A&4&^_GCAY]_
MVG_;8FS?]R= at KG9(F <$>?[("?IHJ?X&\Y3 D/H>3WL VM^'19+''#!IAY>3
MIZ/F]*/PB  39 %,)89"9//ZIY]_V]G[ 4:X at VD61FOL/ Y&E/P"R^#,U2PS
MM\:>?LL dot ?!3H[ dot >^X_J at -/;'6/'Q8U!-WT?)" N]?\48"-Y.I]EYW'Z^QC[N
MOP*H0-_7_ at WZQO,/@Q" 6>EV 8(?NMUZ[=_#V#D:.(1X2F1 O=7^[8= D50+
M*K$?7WUZVWW]T]Z[G1^Z/\+#T.V//9]M8<Z9X*AU_$*KH-Y1\/,+_5$ at 7!]8
M(7V\W+L<^<WDW!FVCI>UQQ[8,Z_>_OS]J]W\\S=O#_(/:8EB>?SIC?D0#,/F
MH1N:#_M1D\R[W-/ ]_W<PP'E/$ERS\%Z[V>>#BZ:/2?Q<P^#L ^3D8 BUE^ 
M(=IT82XNL#SZ at -YQTYH=CD dot 7)!4[AQ4X+/=@7D&@^6  ++]>5IX1M?KV>"J\
M-Y at *KR[7K^ZQ$S?DHM7X 74?IA:U_M#R7:2%R)7GS\O;K%[FX=J$PLAI:_1O
M77FQ=?C=B dot !__5-=IK9HV*#AW^'?AM9V= at CV'I(^];"_JS \J1\+BAM6..T=
M> at %U^&:G D@3NI)?;)-N[[R?)-3[[OZ^@K>(1-0(EJS#U#M?LKU()V!"2^EP
M at 5GVHI%00RA+(U#E 1X,=(4]?W@X8C'H&^0[8",0L'B#*\K?$>@?4 >\%>GF
M$K6<^&A,R6Q 9A\<^Y=<,_1\J$)Z.^!ZU_&\ "2_$_)&2) dot Y at :0(IGEMG>;Q
M6M3\4F=DWE#V&#'BG at \##K&9Z) WXSON,=YZ*X;=RLV+"T !4;Q^]^Y@)S,K
MPOZWT2[6 dot N2UWA74FEB_Q^M_/TM]3_ at &WTC?8(;QK7+ ((/TJ:UYX7!X0PZ'
M:<3?4CE9%])MN7!;*A$.!FS\45:P"1+2L1<>^82]O1_>VL%+^[:*<P.03 DE
M\E at  dot L9^_@ F^QG[Y_<>W93(@+ULRH"-LAG!E^9Z"$7:T<U"AGUR'5G(R.U5E
MA02QC. ,<R#]\ONG5[NSP5K\)2MT<_C at NTC2^]AW$K1 dot Y4I>9M+C==6.[D,6
MGZ%O:P/,0-P)7M_B[;T at WW=\IM^Y2,=M_RK8 dot =*+:IM'RA%")VN@D#QGHV62
MR Y./[+VD+E%XW.-P74FC,\M&Y_T++F5QI>YD<4Z.-L%D!-&:+]W1K]P)G??
M#)_QO.>0%B/GM)J"50CNFN#*S&.#X,+WFJ/+H<^7)M) Q%9 at B>>_+ (9Q<0L
M()=!F]]IK++#*-UWNG\LLYTHC[!+4-1NX?9V[ at IM(V#D4*53I+R9ZA:1__U?
M>:"!/U? "1>E-0T<>O'[VE7O!1?/%P[V0?E at 12\/> I*!,W1:?:A:C" %6:<
MRU3JD!7$8OY';#"9,YI6A<)0D)_*R&5(+.]/$S36+OF&!2"$IC%E #D=N M/
M2:#8"W72A.9#OL!+<W at #+Y@Q+]*U7LEG3_$?PE=_2/YZ!&)9[4V;4"L/NY6P
M<3 Y!D;\K9HL;&PE5L6G*1OG1&F*.3MF7VS?:I2N!DB4<Z/4$.?%")TR$V\&
MS70 JIA^Z534]>,:G6XKF*Z>^'LM'6L3]:D:6I/"+&=06WR*5FVJJ\H\\?$)
M2Y92[4T3/UM\[JS"V;'RS1I3_<AM^*GA$<36$=0F_L;KXO?Z]/#&'45VZXK"
MU# dot 7G at D(L4P'69SGBH;RZKO+G[NB#GU? %4B^!-I4D !S+)>9W\1H.) 6F<3
M!PQ?UC=QM/QP"VUJ81T<,Z+,J dot -"'1?KN%#'U>K4,D3 at \K/WO* ]4V+,C];Q
M_H,B4V0"+2A20%] at Z)6;(]/1N6ZQS#)9T3GS@H$?DK]A$"0#9^0>$PKX6',S
M9\W]3>/K" at &+=(,O,RFK^0XK"A?M:G2U]2F%,Z$ dot 3&:^?YM[[= KU9=NF^13
MU=FQKMLF\R&^<V<0KIN>-X#RC/DR+[FGT-P-_'=N!NOF at O4 dot B9B<<KC=XB49
M.>[IXF2Z:UEM5D=X;6EEF9]9+T(Z#;<ZE<<:WF,KW at GGR@D0NK/A776D%9D"
M[PN1ZNZ<\N5FD-^Y,:0O2JZ[,VG3FR)US?UX0WA?C&2_%S05<4]'<8=1XB\ at 
MTS&_AIRNEXO55"B*R at =)E9]%G'2FTHQQ$=[D7!"&<4)0A:8()/ HH$GW&Q5X
MGJQ9K*3W*;T<B)]NE91!,4C>&&^(U.[]P4KD3[6XFTMGD2<1N4X L,<% B#)
M4%[_L3+!?9#U;5_'LIIHWC_O$BNLLR8V#III4WM%'(&O7/'*PA;8PAHOG7$]
MZ>2/KR?=;( M&0MV,V1<K-N%OZI##JO.ZHGAH3+#M"3VT9R%#H7H"UJU6AIP
MJ:5JL"0'*A+$ at 734!)JW>P9#.PB!< */1I@ 9+3'DY._:52G(8+U8$\A;[E\
MG> at U,\XNG^0$C)6R"U"&^TVLD60I^ %+KAL=-!*<JT8BIXERTOWG\#_<?8O?
MW_W'V LB$-KJI+]6J:]5VK54DA[BIE at JY@(%KAA@ dot @^N91XQ(X:+2#FK]G*D
M/+MG6X(B!E ,C9'GQ)C!3 84:4,4TW1< at :9/LEYLXT1]"0:OA++GP,XL]"UA
MF9;$%0[$9NIBJ'Q6X/-Z%5KTXR+CCD=%I?=*P5M/T]:J+H\LHM=UF_+6"HX>
M4JMZ3XMO/;V1BVIQXP#HWO7U?O-!A=4,6KF$(.<0WS$HNONQZJ32^#%57X!W
M_B6?QWBJF8?:%41/(9B,A_?S!WK*A.'9*).3:,0:0WJ\36]MUR*FC<CC#5^,
M1N3EBL,OT,B7HA:4/:CYL=*Z&$-.1F'9_8SX28/B] .J:SQ1 +:07MA)<*W)
MF1Y^J:>M3+IT>I[8-AI7RG=:<-59A,$24%8G[.ZP/TZZ>)'C-G^ 9F"KG69"
MR):BQ-,R'0*?;$PG((OS1_ at ,=_H>".K6DQ2K.6!T*5\RBA1-(1.,^-6$]OP*
M<KZ\\6!PR506B=S,J' dot V93-#C: at 3=/@9?L'(XWK:Z*2)JC!3%!A6-EM>0/>"
M<T%8?CUX1L ,$W_L15TE9P1V1I'8G=W_]$99\?)V5'B<.BK8/A$\V>!)$!Z-
M^P[=\(47$QI+K8]IN;Y_..JJPD;8K2C\*2U,P3C6TEEC( DPC&&;[;=X"MFZ
M$GPH]^AMN at A$;;4YK5P$M dot "6(9W[,B26W)O&C=4E+$5[KPW1J3@UA^'";[[?
M[;[]>7]G]Z<]D[U$);=2)9K;\V-DQ!4Z2D-[XGH]8.\M/HO$-<VF>9+!T.-R
M>1N*V^T1/*M at EC dot )Z/K84LO6-G>NQTT\8JBDAD8B;Z"\]G-%#%2K']?KK)6J
M1=6() ?L\5.N1]?L48#S"1'_!O[_&+<T#X]I%.M!DFH_/QK' 4;Z9^V%PDC]
MZL2##X>C! :=.-A8HNM:%(BI?LWX2[$:)R<D+TPE3)F$70FU;)!M<X+.Y,S-
M-6.I)QM+5_SH,7G2P#JKG:<T5*X^!29>L/,$3\"OA"'!G$:#\^?;_+U- ^:5
MJ*I<28,6G7KGV dot 7]E[ at 8)/;3)4'1@7C>X$,YJ#_P=^/D3[./2?&WRHZ[IZU_
M/&UIP1*5:,O:!Q:T=*1:.OF3N^_4KT=2I]>I\E14*XAVW7)GS#W9_E/(=EJ1
M*+%*"-0Q2FBVH)'C1>5KX.B;C+4,]LG?ALUFL(_/B1BO!7.YW1W==8WL)W9P
M<%C\OGJ%S*!!&V$G?QH+Q>)IP";,,5;I.-\=HVO=<#0:Q4^KU^Y%Q#]91-P>
MS78O?*Y$^.0F>/'")S^=)<+HS=L#)8(\']-D!*$3%ERS/L'3:]3/^6,+>YK&
MW3MO%QF?K]XC5)['RVOU[$HZ$M!X?ZS_J9R=?[1QSLB[)Y]T\$D[?:)D%8YM
MQ:N7B-?;ZND5.LOA^^.3W+U8;%YG;_R/</8J)ZSF]352 at TFJT;R994Y<3IK3
MNF]37!.>UUBG?8,.6PG?4IYO9G'H)I2!,'-HM3>E3!2-]-;8I TJ>W<\)T-5
M\;C WDHDY?Q;8AEA*3;'](VR!UH^QJJ[9(+O5I:%9]MRZ$[FA,3<0Y at QQDF3
M2B9W; _M7K)>DV1=GR!!E;3I;69WQL1^3_%Z0G5%64B16WI:P*2>'U$!FK%>
ML5KF6A1C/E56G9+YI%TW[.PA'_ ?84PK#13F+$UF7BK*!4K$D"L)]J6R^4XH
M(JXX,7H^A^,4BVU#^F7S#1>OO)'5U-Z=Z at TY2PK>^JR=%LCY:?OG@G\1HRB1
M_[,-RM!&UI- F:'I!^?M:G<&-3AUQPM3P//T7#H9N4-<5ZF2,X$?]TKY:U?*
MM0GV]U>MEZ?6I dot *H!5>=U[H dot XOIR&IU812#E"BY2/\XZ at *O3E0L8T=7JS8PR
M[R?3+5+I!(\3GFH=JB;2[J# 5+T5S,>5=R:^4A,5EJ?'2<:\YB\&%AUI49&#
M5"VJ;X.K6J8&H:D4?0=0P ]\)NCF!?H:%&C&Q6D at CC4T(0:X]<->PM\-%F8.
M%*G E) at $X7%B<6(;9\B/D[+0K$%!9)8,1^FULJ$3AG+6!#M5+@$-+)EN0L!M
MY8$S+8&$XL^Z2=X8P :2HEYTU4<NNE;;',%Y%)\J(X5H#(8"L-,+J/&X,0 -
MN;+>H D S*GE$(Q7/H:O^&V#J1OC&&8A6R&"Y^_I*Q;(:.JTGY#W$]K[P<<#
MZB>$9 at ;6?@9I/X-ZD4E%_:U0MWE<BL'@GTI&D99!%33R *B/$^ :2[7S #=8
MA$HF"ATFXHRIL)90>.#V"W2ZQB=$J(5K<VP &%8+*L[<OJ"XJ^CZ!8VS#.[3
MMY#TNR0*[E3(LUW95>'%"\V, )]YE3U!?Z$ *U,K,PUN% at OC!L96H "KVAGE
M0[-80KF135ZB3V\'S-#S8 at R0^3HNG(NRY?FQ dot #6Q,%/D*I;G=\$8L5SF329)
M^8ZZ,C;RER7WS,QY]T;&O9%1V<C at !@98&G?/R,B8%T6F@VDXY*V&JB9#-9&[
M*"=$*OPG+_HK#6) dot 1T16&2U^5%=H(-!%B-" at $XX2(ABA0/R+810"G>(5:*[3
M=\6=0*U:#1/T!ZY*X>]XP'9__,FV:4A/6^W,QV]VUN!%I]5Y9GR>JQ?/<B_6
MQ8OV^M-G3SKRWV=^\S%_T7[VS>/GS]2_3_SF$WSQI/7-TV?K3Y^H?Q_[S6>\
MQN-OGC]M=])__>9S_N+;]>?M;YZN/VL_Z7SS]-O'W_C-;]=J?V]:4E(,!_G3
M5'D?A7 at R2/5DQF$/4ZRL $ ^QN_X+#D.#D>,@L3238E1S!]+_SB4_AD%0^C1
M)95X(T(R\H>LLR%:":-XX/2#+W*F"B2 JV6+D'VL at F[+Y2^5+Q_Q\X-DMJS(
MAUMIT)<V4-!:2W1;1/B?$3O'0UM +B3% at B,_%%%PHXCU\+I-GXU#_EHW dot B: dot 
MF9E#9<UM dot 8)2-*UOL)[3=T(7GA!R^#TZK_9>OV6#'L;, at P)8_I[+5YP_?-CB
M57S/&I3GF:6 at X4PT7M%8'F_0?8UR)(;9@WI%IJ00E !,V\7))?"W2^X\6-Y9
MYM&8^/\ at HQG7X"'I-(O63/UHR><AG='$*:C3X3SJ&I/ M-I2F8.F!RX$.P!O
M%%-EZNP1?[ .E*$4_X8(G^/H>$V7/IX?^W3EUKFO[G$5<=!XC2O=X\-U*1^-
MBL^3 at VMOF@5RR83T:S^[APXAUCBE9Z,R:DOJ1Z-J8YNM:P8940< JY?1CF<"
MF.%PN,:\X7"#!.1_,/]?'%VP8=2_#*,!"E4B%%=<!2OH"6J)H$,W#3J4M#9$
MVQM*2$SN =PP7C] dot F//_OF&M5 at O^=M1E*(,@Y-LOW1-QKBSC;WB dot _H87W dot _0
M;&8/BP\)Q 9]6Z5O).G_.$E//JHRGBJSHO>:5E$FB_%:CBK%VN]^''&8" S>
M0],;#L7 at 2[ dot S:&-/TT"=U%&TZ:<S2][*8<B;@7'X8N:" 0E5W'A4@<'#89TC
MJ&4<2,<>A+Z%=XKP/_B\WK! 'L3^T*<[#( at 1I&00(:U$F2:8JC_QI:&E;E*L
M at E&N dot F2E@P !*8$S!66FJS2&E M#E()>+K T/S[/>.Q))/$GMCJ9-BN  *KP
M$$VA G5H&+ P(E#J3"HT??<BS6BM;BM8M=Q6,$5>M)K,O35EEG61E+3+\U:+
M1';=;+;JHS at 8^B""8:D1XD7*=%AZ61O[<AJGSV\+IC:T\\E%F:J+3D2;55)-
MU4!BR"S(<-;>ALD8_0Q1V+]DN"<<^X>P) +-BK)_< G3&9SAS= at /L#''\]:!
MP0KOE*AE!U V<\T[/'/-NS9SR;A'$>F+F[Z4\3+7TMR1*9R1^6:[0&=UIFL(
M+*QR-W$]([O,ANMF)5P#!XU#)P;6H4'B6:C<VC*]$F?6M)NIF=C3+OC.6GPV
MJZDTKR;'6B^SL_2 at * dot U=3T(L,E:=\05-[Y+"P<2OY@OI8K!@1&&B8??4G\EK
MQHPDR"7'!7BJR[/,+3%A019,08OXAL at /\UC6*M!=8YD$'%W2U9$), dot T'WM+F
M\4"4S7E'#<DFEE+Y[):?<S+<G>J6TS48YAZ,COY)6^8C;;55HIXS0U(;DIM?
M/)9>BB=!W<7_Y?&H27Y/A+W:(34J6>("'= F1;E;$CT0(,"B4*ZS+I (<;7F
M5"7 U8+;<^3"TB Z<KW9;T0SG'$.-(N91+5S,JGL6E"'S<(>'SV2* !D9/%3
M&2MR/)(3+9IB2JQ at N57F%&'E:CILB@X!*2F%S(P4T[X<E Q5R?-!5C621!_D
MSC7=D-5229&^UWUBU6U/64U75)EC2(-R05WM&JWL$8 at )!E(^(W2>.;^B*6Y^
MI5-<9I>53'%CFBDN.IMXI9,NISFN;),43*_5*K''#O#1<&,EN\DHFDP]E=Q"
M+MR43(=0W:I1Q?C1 WG&O-" *SA4)< at F/56EG4('^O(NH"@Z+#$K@KR#5J0"
M!_$.[\U#X([A0E3'VT7<5GFJ^B79.$\X,,B<T)(!45K:V5QO[9E[HZ#T3!^S
MM60 dot JE$(Q]_*/"2"WQ)RRIP$HV=*$T at ;/M0__'"Z:O]G:2D+# dot ]D%GMP-B6O
M2P,E(NZF!KAI)>\4BOF!=?5MU^?Y?<LI]/E7-)LWK<^;-S>=C:N9SFEGTZZG
M*TY?B;Z>:?HRVKM at +G4%GDZI:VCKJ295-RG,;7JE)O\JUM)I##1IFXSCR;&>
M9M8T9:[#]GP=< at VC=S-'8^;@&@70Z#/Z(G<MZ:1^C7E"*LI?E& dot %,:L[W?S:
M>"IM66Q,WQO2TQO2Y,CT&-N6]BE_X, at YKBW*VJ[@0TP'B__K+L2!YD+,[/UX
M3LZ'B.W<<A=B1$=44C)OY1SW V<(J.IZW?=XOV/,#BUGPWJ60WF]%I3O7P*\
M]EUG\=9H60]77\Q>G4VHVT2ZQ\UMF![\5K#/T8NB?CK^\+*+4 at 8G*DBZH7\$
M&#R;^5*QVF+V,9# =-DG3K M2:] at //9UV X=X$""K10X//P$]!<ZX0V#!Z(B
MO7F;9>XF6D+ at +X($1@O#=/IT1SO\YL dot &WR)4CR#O2!FD79?6MLQROR\1D72=
MV =4C! 5=PLC" at  dot (D3?=O9V] _%[>YO?R8Y*,PI'08CTL20CU$TJD?A*$<;)
MB4N2#_P1J) at O>,()>@8B8@)W&(SQG@4)O :EA5>AC/PC/VXQ]JJ?1/(V+FP&
MY!'K._&1GV!!CR4#F +\(<Y5XC/9^S& "M+N_:M?NY]>[=*K]SM[^+UEF4;1
M9R*][:"PG0N\M" - at QX$89?0,L-D%M]4+Q6GZ$Y-DKK5D%[RKFTOC21ZYI1<
M!Q%Q*L(7+R0(1"\I.-RL4<6V%!ZIF ),*V;0X(.4!BL3G#&[HRCJ$LUT =;N
M83]R"I)"W at ZNU!'Z;O>@"_2+?,E11P]V]G1L<*#SR$B%-_#-K[_^RM[M_/K^
M+7U[$R%WT!4\0/J7K dot ?WH_"( at :[W7S)V< dot Q?LN0X&H/!-7 N>SZ\IRWMON^$
MOL?&0^:P7H A(A$4_^\Q'ECM1>,1X DLDV2(9]* at N<-Q2)<U\4/UT"$VPK>L
MDY?YZZS[_07<N*G?W%C.;EK"X()=:-P?YK:C$21"K,?$F<FB"%;7[FE5ZK=M
M=61:^Q$ dot 2RUQ5#9AE+;\4WF/%P)1P46C5H at "WL@B(2H$0MKT4Z ^.Y 3'>")
M;"O!R=I/*;BBL[+>%-CXJ9 %K$)B1]2A=XIIVO,SS8,J)-:Y-J;)030]TRP4
MHBMGFO:\3/- at &J;I+)YIW/%@&$>6Y/0EGICJJ<<R.><+YD%C)\%F?-ID2!F.
M+V/K+6#JL-E-G0IY W1)(V"3>FUHW9ZL"LK/H#C/ dot 24 at DYO#!#F> >3);)7-
M =&V at BRVU*P@!ZL=U;$%Y$E09F<V)TJ+/:<%LUB FY/Z=!QHS8XA,+-D0PNV
M;$%-9EMS2J9+QI7"!V^*YW!XBV<Y:+6$X[#/U1MBN GPSLAO67A-NC+ at O:W<
MED',0IB-HZ60U_)XF8/5;E*YI717;,1D<&F=%TOM1F9.9N*-V89HG>ZR(08%
M[IF_3"&T77P$NM)2S;3["L1>)7S_9>OI9+I%H:5Z(^MQG(*,;U!=3""1]AQ4
MO'H]5&P;8E4J7KTV*K9M+E\-%5=?I5FJK\Y'Q<GGNT['NME at ,T3M1 dot 3Q@-=9
MS9 at %$KE]_(&R/J89?W5+\^II-6NY&+!I=HN=IB5\'+JIW$Z]9/XP+=[\XF dot U
M9+O2]#JDT6IHR1RITC+Y6/+4\%NCL["*ROQENSZI dot  at [X-(<K*[*R3CKDEM/4
M$8?OF]OL=%/GIM,T.P0VNDKO:YF<&EXJ/#13V>KDPW'0P+=YX2TL5\>$5?S"
M/)$)A0 at O-<FQ1BJHM7%G['!]YK" dot 3M">GH@I/2NT>FH]Z\0$X'-VT#SEY^Q*
M3SG-, dot JZ77:X? at PHW6++V PF/?@\]A/,>J#.Y^)62G2(&8J._/\)EW5R]8I)
M#2BI dot P at L&\YZ0OG0\R^T4QBJ#KTHH6-5,&WL(6^M)$4@/ZPX@VXK\4#+"^#5
M0C!=!.)HS!=5E ,3'W7,B!^%;QN'[S%#&H*O3Z]WH9W:D;OHHF!=4 9O2]QH
M94GL;LB[3E9%B 5E.H;\5F&^>Y5U)G6L:J,H<J1J7(7];$F(BQH[2=M*<4-1
M=A;G+-=TV;R$*] at 2\B\(%G& dot :<_9 _DB&A3'O'!2M5H4B"#GS'9\29V0+:9C
MYV)Z)L$ZE9 at $"_Z#F<2YJ,@DSL5=9Y(7$N*J3$*X6223 dot !?S,HD*(D\=M42>
M4VL2O9K&)^7ME_*)JGI-3.(6,4G5Y>=LFL2[D&2Z&$T2<":QNB9OHR8)YM(D
M)S-IDI-%,<F4FD2O5HE))BJ3KYA)=$U2RB0S:));QB03-4F62:;3)!68Q*))
MIF22*!D!E(.':00$K'%6Y%,6)=;PAT>/.!4<!G[?ZYX'WN at 8>HV2%N8V"RA[
M]TJ=K;+GA3$X:=!#/ at +'&BNB@B&$ 1$EM!J#_QXU6.*/SMF*-IHZ:SS"]YE0
M%=Z^J$IK-(D at @9F(GPD/<EAY\0(:D5@) "M%L5FY:([BR+DM=>G]EG2J!4D+
MS^LA_%&RL=%SO%XPLD=LI 1?S9LV:4=-L%N0(*R<2.DW[GPE*CK!0*BB9GR3
M,BYC1]$HPE2J/KU3!(A/-C0Z%"?P1_Y at V'=&/MMR^TZ2L(,7,E4O"ND:X-S#
M<[Z9"<##HBPY=X;=WN7(AP<'#_F)4YJ dot I(75V,J*>PS8;-390SH3C&(/UN at K
M!_5T)M(V3#),S at /,/K^BU2'W@NN #.AL(&1"0HAGZ_2,FEOG#6;ZYZX]L]*3
MM-*3RI6^22M]4Z62YQ\ZP/%4JVHDRLJR1/P&XQ$I1R% at PF,4'4/A1J,'1CB*
M/I4X<RP_<WS:0*; 6!]"W:KED^ (S[!.6VT<SE at Q dot 8[BT4P]35T3'L[4T93U
M\*CP3/W(BI(G at Q"O2*12P-!=SF8&9P)146+Z+FB dot LXT-^CJZ'/K,&ZTQ("Y\
M/W#<8XRWB#8VZ&AWEQ,4.^S#U\$HS]Y0CQY!E[P"-N['J&O2H_Z*_ZEH+)-R
MJ at  9R=+>2&=U8JGLH*$#H)F-FK)[B+VXN-.QMF;((#0MH(#H-A7-TD<I.;*D
MVR33KT;!U]'].-.]P4+7 C\RD X__KX6R+,]FSQ]'4- 89OVC[U<!]QFK[J 
MN8[N4;YHW5-(PW5 G>G7D'?7,0"28MH >(3+(462Y*V7P]1VH8)U%7)B$8AU
MOBCKZI)5>\U6L+G.FA2VU!0'Y7 &0+CPU>EV at O'%*U2"@!>U@D!-&D 4CUS9
M/_R54 at Q&2%(5PX@;191NVNF#_B2UE@QAT7,8N!1LF]YR(@; QU2S;5T"::5>
MZ+P>Y6L*\1<SX/-V2]4KEDU dot  at Z%1>H*R35<P"C%-?B<4J='H5%>A6(R"FYE,
M9J]">_DE2FVY\,%0,G(SR()\Q1 at Z?1[WK@J*9VY:4#O\R&F*PY(_!&J>+)1Q
MZ$ at [?##96/$EI@]+)NPI dot (J?-OUX?::F'Z]/ dot H? N\#F.V80?4$G_*76#W:A
M**QLKD,*9.Z**:<ZVZ45Q%JUR.A"WX1HC"[I-IKG57,!;,TF0_G)HMC#R'!8
MD+Z.G<N7JIM4PJ:\*< at A,])S:"+I]H*CKH_[NNAN$/<LZ$VH 6[G0 U\WP<^
M'XU NO FT"<#J_'B*F> ]:-JI<P5^U1C2F'2SE9L;IJ,8,,LC_G#R']?A at '"
M;X=+)[H+!Z0WM)8V at 7*"QT+IV_^TV/]+NNRR?D;\)MA_6]%A+BM6RN1:7C*-
MS>DH0L;+)JA8]$MC$[EF:61)*_']TR/^!N0<ND= 7FO>1>6ER)ZO$1G8Y#W?
M)!51IFFKES42GFKQ8:P[N*_ dot NN[ at 9H!^<#HZ-2-H:%!9Y^ A)AG)XY#7T^4(
MWMY!#LPE4RB*,>DB2M5Y1#I#2"OIL\S6D*4UN<.+% DW=%0J +#0'P3!ZNJ?
M>9>FZC(/-DPCF0!U'5P+9'_)+SJV\-8V*L$K+N4)37^(="8JXI%MP at M>]\Z+
M_9UV(/%C:Y7J6QK.MO)W'M4:"-8&<TS1V31:T_!=%$"E.C1YR.Y%)Y(70M7*
M>C8EF*DD*V0R_B#1Z^))4IF$>TV!J65LX at 6%O2'0HEK-VTH3G(3G<3#RN:UN
M.*_S7D+-27CO_"MV_BF$5O7^H?V1M*A:N=LUYR:TS)WN)\PYK4K*:[Z*::H9
M3HZI^D//P$P]35T3'L[4T93U</D[4S^B8H&;4-"%P9OR^&NYQY#$T(P>PT+K
M]6J\A ::3(<!,42=>=G5865/8%G;.O%.U\6X:A<FATP)A^Y7*X6#>&)*""JW
MGF&]Z;K1?&1E?2"S33G^BBT;_#Q=%[JGJZP+XN I1U^U;5-,3-=) at :>,+Z;0
M%T:GI!?I#+,I,IL'C at ]I+E_9S/XOK]*H/9O;K7C8"W24T5"NPE-F4R73^<2N
MP2G&#TWH6>]$4L'TJ \\,V]PN7???*7NF^R9C(*LLY'A><EZ0"+N 1F6>$"B
MG =$^CD4SZQ14L3)[ at ZKT91Q=23H# dot ( dot #^$,V+R=Z]+L,M,J8E2XU _ ":%*
M*,<NE*%Z::2R=&$QBSM!8"AW8VI3O]TV'AW1[;8?#GZHLQ5HXW+-K&$$5>'-
MLVMLG3\[RL3+N>ZF]I3BZW)/V^)IDBE*#303'>0C/31L_[)_AF$L<>ZFP>RY
MKU[V at 5N4!]<F-^1=N2Y=HCNBW&-KV":F]W,9W9_;[U dot >P,09^+P55%PMT:*\
MV73?/1['H*9<^!TE at 4QORMC^ZQ\_?F .S!LTO/R1*Q?^$ /)>NE#WN !WM")
MLA*ZAX&%_CG !V(I"#'-J'ZW[9C2U2:M<1B,G/C2>I\RF.I=42S!$69N4Y8M
M]:!(;W)+HEBN)1J[S'2-9PG'CG'3*&L , WH) at 5SGZXJ);RJJ?8_CTD!PP0X
M[FES)#&QELG#*&8AZH_5/:0HR8K2%O?X\U[Z/&46NF]8%D1Q:>P<#1%UA,)<
M<E^M4$\4ZI450JRX%_D<P4:*X%$,R #V//BP_]NND22X(S,E]^C?H2-_#WNI
M:!3O+N0[ DY>4*^E!A92J"0U<#5KJB0E,$&R;#\ at : dot -#SG\XTI><1 SQV!P[
M#?>B,>[E[L0M-,R0+_O^A1%8\\L7BN2<)%#J(JXFUT VOMH at N((TV9IA9="A
M09XN/=:+R]L]MFW7%NAOJ5%0=S(;=H]?!J*T-973M+52;X)=3\)37D at 2C=I dot 
MY,I6D"\6V\;"UBS7^#$0Y?0/8]ZLI2UZN4UE*K<6E+46\-:":JWU?+1W"UJC
ME]M4IK U at 3F';]\XEB'A8W@[J84>;Z%G&08^+FN!#H. !/["#7MI_:1";,G"
M:R]?OL1?F=[\(=J/;[[?[;[]>7]G]Z>]M+8R'#Y_.<;5\B^___AVOV[2"XHB
M$$(T&/ at &)%*W-H )/G_Y?>>@H#J,0F^#2ZV,8?> /U3[6%KKY,W]Y?=/KW8+
MVB>JX^G2H291#;Z ><X-?4DH=EPVH>(1^8)QKP"-U2,_!''7)Z>S'\"O,\H5
MW!+G^7'%I=9;19< B8L$5M4U0  at >#@8#K@2W"L!=;K_2;YM, GV[XH8C.?(9
M^C9WPK*CD)MA at ))W90AHB8TZ]V*%AHP&GA@OS @@'QI=(S:%+[CK)3H2>WAR
M#TG+>40:<].RJU1EG8]["LM$%QNP'.:4#HJ)TR&2$;3^?[QE16=VF]Y^D0$U
MNZS?72#U<B^]8,0X6%%)EXK1TFWU/K\3PH6U'49T!Z"5E^L%"X1'C=KFYB;;
MC4#=LT].'* N3C98H\%?#"+/WV"O5U?5D[>AQU\W'M6^ WN at T^K4^M%1[;L=
MF&#<DX_],XJ6 at +<C(/_:=U[G\>/'[&G-@;_/V7K1 M+55I#F[1Q-]$K0M?";
=T.C]^;^O^/S?O^X_]Y_[SS_V\_\![ 0RG0 H  at #1
 
end