conindex.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file conindex.h
///
/// Conley index computation routines.
/// This file contains the definition of routines for the computation
/// of the Conley index using the CHomP library.
///
/// @author Pawel Pilarczyk
///
/////////////////////////////////////////////////////////////////////////////
 
// Copyright (C) 1997-2014 by Pawel Pilarczyk.
//
// This file is part of my research software package. This is free software:
// you can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// This software 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 this software; see the file "license.txt". If not,
// please, see <https://www.gnu.org/licenses/>.
 
// Started on July 19, 2006. Last revision: February 2, 2011.
 
 
#ifndef _CMGRAPHS_CONINDEX_H_
#define _CMGRAPHS_CONINDEX_H_
 
 
// include some standard C++ header files
#include <iostream>
#include <algorithm>
#include <new>
 
// include selected header files from the CHomP library
#include "chomp/system/textfile.h"
#include "chomp/system/timeused.h"
#include "chomp/homology/homology.h"
#include "chomp/struct/autoarray.h"
 
// include local header files
#include "config.h"
#include "typedefs.h"
#include "mapcomp.h"
#include "indpair.h"
 
 
// --------------------------------------------------
// ------- the CLAPACK eigenvalues subroutine -------
// --------------------------------------------------
 
extern "C" {
int dgeev_ (char *jobvl, char *jobvr, int32 *n,
        double *a, int32 *lda, double *wr, double *wi, double *vl, 
        int32 *ldvl, double *vr, int32 *ldvr, double *work, 
        int32 *lwork, int32 *info);
}
 
 
// --------------------------------------------------
// ------------------ CONLEY INDEX ------------------
// --------------------------------------------------
 
template <class IndexPair, class euclidom>
class ConleyIndex;
 
template <class IndexPair, class euclidom>
std::istream &operator >> (std::istream &in,
        ConleyIndex<IndexPair,euclidom> &c);
 
template <class IndexPair, class euclidom>
std::ostream &operator << (std::ostream &out,
        const ConleyIndex<IndexPair,euclidom> &c);
 
/// The class that computes and returns properties of the Conley index.
/// The index pair object is used to get the cubical index pair.
template <class IndexPair = EmptyIndexPair,
        class euclidom = chomp::homology::integer>
class ConleyIndex
{
public:
        /// The default constructor.
        ConleyIndex (const IndexPair *c = 0);
 
        /// The copy constructor.
        ConleyIndex (const ConleyIndex<IndexPair,euclidom> &c);
 
        /// The destructor.
        ~ConleyIndex ();
 
        /// The assignment operator.
        ConleyIndex<IndexPair,euclidom> &operator =
                (const ConleyIndex<IndexPair,euclidom> &c);
 
        /// Sets the index pair to be used to compute the Conley index.
        void setIndexPair (const IndexPair *f);
 
        /// Computes the Conley index. If requested, uses the two-layer
        /// construction. Otherwise, first tries the flat model, and switches
        /// to the two-layer construction only in the case of failure.
        int compute (bool twolayer = false);
 
        /// Reduces the index map using some shift equivalences to certain
        /// extent. Shrinks the matrices and generators if requested to.
        int reduce (bool shrink = true);
        
        /// Truncates the Conley index to the new top level.
        void truncate (int newlevel);
 
        /// Computes the eigenvalues of the index map.
        /// The eigenvalues are appended to the given vectors of doubles
        /// with the "push_back" method.
        /// Returns the number of eigenvalues or -1 in case of failure.
        int eigenvalues (int level, std::vector<double> &Re,
                std::vector<double> &Im, bool nonzero = true,
                bool sorted = true) const;
 
        /// Verifies if the index is trivial. Returns true if yes,
        /// false if it is non-trivial. Note: This information is tabulated.
        bool trivial () const;
 
        /// Retrieves the nth Betti number.
        int BettiNumber (int n) const;
 
        /// Retrieves the 'i'th [torsion] coefficient at the nth level.
        /// Returns 0 if 'i' is beyond the available scope.
        euclidom Coefficient (int n, int i) const;
 
        /// Retrieves the top level, or the maximal dimension of the index.
        /// The numbers of available levels are from 0 to dim inclusive.
        /// Returns -1 if the index is undefined.
        int dim () const;
 
        /// Retrieves the matrix of the index map at the given level.
        /// Returns 0 if the map at this level is zero.
        const chomp::homology::mmatrix<euclidom> *Map (int n) const;
 
        /// Retrieves a human-readable text representation of the map.
        std::string MapString (int n, const char *newline = "\n") const;
 
        /// Retrieves a human-readable text showing the eigenvalues.
        std::string EigenString (int n) const;
 
        /// Retrieves a human-readable text representation of the homology.
        std::string HomString () const;
 
        /// Defines the Conley index by the given Betti numbers
        /// and the given maps in homology.
        void define (const int *betti,
                const chomp::homology::mmatrix<euclidom> *matr,
                int _toplevel);
 
        /// Defines the Conley index by the given Betti numbers
        /// and adds the identity map.
        void define (const int *betti, int _toplevel);
 
private:
        /// An object that feeds this class with cubes.
        const IndexPair *cf;
 
        /// The number of the top homology level.
        int toplevel;
 
        /// The torsion coefficients of the corresponding homology
        /// generators. If delta () == 1, then no torsion.
        std::vector<std::vector<euclidom> > coef;
 
        /// The matrices of the index map at each level.
        chomp::homology::mmatrix<euclidom> *indmap;
 
        /// Converts the Conley index to a string (without the index pair).
        friend std::ostream &operator << <> (std::ostream &out,
                const ConleyIndex<IndexPair,euclidom> &c);
 
        /// Retrieves the Conley index from a string
        /// (without the index pair).
        friend std::istream &operator >> <> (std::istream &in,
                ConleyIndex<IndexPair,euclidom> &c);
 
        /// Rounds a real number to the given precision.
        /// If the absolute value of the number is below the given threshold
        /// then the number is truncated to 0.
        static double round (double x, double epsilon, double threshold);
 
        /// Variable that stores the tabulated information on whether
        /// the index is nontrivial. Values used: -1 = unknown, 0 = trivial,
        /// 1 = nontrivial.
        mutable int nontrivialindex;
 
}; /* class ConleyIndex */
 
// --------------------------------------------------
 
/// An auxiliary class that represents complex numbers.
/// It is used in the procedure for computing eigenvalues.
struct ppComplex
{
        double re, im;
        ppComplex (double _re, double _im): re (_re), im (_im) {}
 
}; /* struct ppComplex */
 
inline bool operator < (const ppComplex &x, const ppComplex &y)
{
        if (x. re < y. re)
                return true;
        else if (x. re > y. re)
                return false;
        return (x. im < y. im);
} /* operator < */
 
// --------------------------------------------------
 
template <class IndexPair, class euclidom>
inline ConleyIndex<IndexPair,euclidom>::ConleyIndex (const IndexPair *c)
{
        cf = c;
        toplevel = -1;
        indmap = 0;
        nontrivialindex = -1;
        return;
} /* ConleyIndex::ConleyIndex */
 
template <class IndexPair, class euclidom>
inline ConleyIndex<IndexPair,euclidom>::ConleyIndex
        (const ConleyIndex<IndexPair,euclidom> &c)
{
        cf = c. cf;
        toplevel = c. toplevel;
        coef = c. coef;
        if (toplevel >= 0)
                indmap = new chomp::homology::mmatrix<euclidom> [toplevel + 1];
        else
                indmap = 0;
        for (int i = 0; i <= toplevel; ++ i)
                indmap [i] = c. indmap [i];
        nontrivialindex = c. nontrivialindex;
        return;
} /* ConleyIndex::ConleyIndex */
 
template <class IndexPair, class euclidom>
inline ConleyIndex<IndexPair,euclidom> &
        ConleyIndex<IndexPair,euclidom>::operator =
        (const ConleyIndex<IndexPair,euclidom> &c)
{
        cf = c. cf;
        toplevel = c. toplevel;
        coef = c. coef;
        if (indmap)
                delete [] indmap;
        if (toplevel >= 0)
                indmap = new chomp::homology::mmatrix<euclidom> [toplevel + 1];
        else
                indmap = 0;
        for (int i = 0; i <= toplevel; ++ i)
                indmap [i] = c. indmap [i];
        nontrivialindex = c. nontrivialindex;
        return *this;
} /* ConleyIndex::operator = */
 
template <class IndexPair, class euclidom>
inline ConleyIndex<IndexPair,euclidom>::~ConleyIndex ()
{
        if (indmap)
                delete [] indmap;
        return;
} /* ConleyIndex::~ConleyIndex */
 
template <class IndexPair, class euclidom>
inline void ConleyIndex<IndexPair,euclidom>::setIndexPair
        (const IndexPair *f)
{
        cf = f;
        return;
} /* ConleyIndex::setIndexPair */
 
template <class IndexPair, class euclidom>
double ConleyIndex<IndexPair,euclidom>::round (double x, double epsilon,
        double threshold)
{
        if ((x < 0) && (x > -threshold))
                return 0;
        if ((x > 0) && (x < threshold))
                return 0;
        x /= epsilon;
        long xl = static_cast<long> (x);
        if (xl - x > 0.5)
                -- xl;
        else if (x - xl > 0.5)
                ++ xl;
        return (xl * epsilon);
} /* ConleyIndex::round */
 
template <class IndexPair, class euclidom>
inline int ConleyIndex<IndexPair,euclidom>::BettiNumber (int n) const
{
        if ((n < 0) || (n > toplevel))
                return 0;
        int size = coef [n]. size ();
        int count = 0;
        for (int i = 0; i < size; ++ i)
                if (coef [n] [i]. delta () == 1)
                        ++ count;
        return count;
} /* ConleyIndex::BettiNumber */
 
template <class IndexPair, class euclidom>
inline const chomp::homology::mmatrix<euclidom> *
        ConleyIndex<IndexPair,euclidom>::Map (int n) const
{
        if (!indmap || (n < 0) || (n > toplevel))
                return 0;
        return (indmap + n);
} /* ConleyIndex::Map */
 
template <class IndexPair, class euclidom>
inline euclidom ConleyIndex<IndexPair,euclidom>::Coefficient
        (int n, int i) const
{
        euclidom zero;
        zero = 0;
        if ((n < 0) || (n > toplevel))
                return zero;
        int size = coef [n]. size ();
        if ((i < 0) || (i >= size))
                return zero;
        return coef [n] [i];
} /* ConleyIndex::Coefficient */
 
template <class IndexPair, class euclidom>
inline int ConleyIndex<IndexPair,euclidom>::dim () const
{
        return toplevel;
} /* ConleyIndex::dim */
 
template <class IndexPair, class euclidom>
inline std::string ConleyIndex<IndexPair,euclidom>::MapString (int n,
        const char *newline) const
{
        if (!indmap || (n < 0) || (n > toplevel))
                return std::string ();
        std::ostringstream out;
        const chomp::homology::mmatrix<euclidom> &m = indmap [n];
        bool nonzero = false;
        for (int j = 0; j < m. getncols (); j ++)
                if (m. getcol (j). size ())
                {
                        if (!nonzero)
                        {
                                out << "Map " << n << ":" << newline;
                                nonzero = true;
                        }
                        out << "#" << (j + 1) << " = " << m. getcol (j) <<
                                newline;
                }
        if (!nonzero)
                out << "Map " << n << " = 0" << newline;
        return out. str ();
} /* ConleyIndex::MapString */
 
template <class IndexPair, class euclidom>
inline std::string ConleyIndex<IndexPair,euclidom>::HomString () const
{
        if (coef. size () <= 0)
                return std::string ("0");
        std::ostringstream out;
        int coefsize = coef. size ();
        for (int i = 0; i < coefsize; ++ i)
        {
                bool nonzero = false;
                if (i)
                        out << ", ";
                int counter = 0;
                int size = coef [i]. size ();
                for (int j = 0; j < size; ++ j)
                {
                        if (coef [i] [j]. delta () > 1)
                        {
                                if (nonzero)
                                        out << "+";
                                if (counter)
                                {
                                        out << "Z";
                                        if (counter > 1)
                                                out << "^" << counter;
                                        out << "+";
                                        counter = 0;
                                }
                                out << "Z_" << coef [i] [j];
                                nonzero = true;
                        }
                        else
                                ++ counter;
                }
                if (counter)
                {
                        if (nonzero)
                                out << "+";
                        out << "Z";
                        if (counter > 1)
                                out << "^" << counter;
                        nonzero = true;
                }
                if (!nonzero)
                        out << "0";
        }
        return out. str ();
} /* ConleyIndex::HomString */
 
template <class IndexPair, class euclidom>
inline int ConleyIndex<IndexPair,euclidom>::compute (bool twolayer)
{
//      using chomp::homology::sbug;
//      using chomp::homology::chain;
//      using chomp::homology::chainmap;
//      using chomp::homology::timeused;
        using namespace chomp::homology;
 
        typedef typename IndexPair::CubeType CubeType;
        typedef typename IndexPair::CubSetType CubSetType;
        typedef typename CubeType::CoordType CoordType;
        typedef chomp::homology::mvmap<CubeType,CubeType> MvMapType;
 
        // measure the time of creating data structures
        sbug << "Preparing to compute the Conley index...\n";
        timeused prepTime;
 
        // create the cubical map for the homology computation
        CubSetType X, A, Y, B;
        int dim = cf -> dim ();
        int countS = cf -> countInv ();
        CoordType c [CubeType::MaxDim];
        for (int i = 0; i < countS; ++ i)
        {
                cf -> getcube (i, c);
                CubeType q (c, dim);
                X. add (q);
                Y. add (q);
        }
        int countExit = cf -> countExit ();
        for (int i = 0; i < countExit; ++ i)
        {
                cf -> getcube (countS + i, c);
                CubeType q (c, dim);
                A. add (q);
                B. add (q);
        }
        for (int i = 0; i < countExit; ++ i)
        {
                int countImg = cf -> imgcount (countS + i);
                for (int j = 0; j < countImg; ++ j)
                {
                        int n = cf -> getimgcube (countS + i, j);
                        cf -> getcube (n, c);
                        CubeType q (c, dim);
                        B. add (q);
                }
        }
 
        MvMapType F;
        for (int i = 0; i < countS + countExit; ++ i)
        {
                cf -> getcube (i, c);
                const CubeType q (c, dim);
                F [q];
                int countImg = cf -> imgcount (i);
                for (int j = 0; j < countImg; ++ j)
                {
                        int n = cf -> getimgcube (i, j);
                        cf -> getcube (n, c);
                        CubeType r (c, dim);
                        F [q]. add (r);
                }
        }
 
        // DEBUG !!!
        static char debugfilename [] = "_/_x.cub";
        if (debugfilename [1] == '9')
                debugfilename [1] = 'A';
        if (debugfilename [1] < 'Z')
                ++ (debugfilename [1]);
 
        // DEBUG !!!
        if (false)
        {
                { debugfilename [3] = 'x'; std::ofstream f (debugfilename); f << X; }
                { debugfilename [3] = 'a'; std::ofstream f (debugfilename); f << A; }
                { debugfilename [3] = 'y'; std::ofstream f (debugfilename); f << Y; }
                { debugfilename [3] = 'b'; std::ofstream f (debugfilename); f << B; }
                { debugfilename [3] = 'f'; std::ofstream f (debugfilename); f << F; }
        }
 
        // indicate the time used for the preparation
        sbug << "Preparation completed in " << prepTime << ".\n";
        prepTime = 0;
 
        // measure the time of creating data structures
        sbug << "Computing the Conley index map in homology...\n";
        timeused compTime;
 
        // compute the map in homology
        chain<euclidom> *homXA = 0, *homYB = 0;
        chainmap<euclidom> *homF = 0;
        int maxlevelXA = -1, maxlevelYB = -1, maxlevelF = -1;
        bool error = false;
        try
        {
        //      chomp::homology::outputstream::mute mcon (scon);
        //      chomp::homology::outputstream::mute mout (sout);
        //      chomp::homology::outputstream::mute mbug (sbug);
                if (twolayer)
                {
                        maxlevelF = chomp::homology::Homology2l
                                (F, "F", X, "X", A, "A",
                                homXA, maxlevelXA, homF, 0xFF);
                        maxlevelYB = maxlevelXA;
                }
                else
                {
                        maxlevelF = chomp::homology::Homology
                                (F, "F", X, "X", A, "A",
                                Y, "Y", B, "B", homXA, maxlevelXA,
                                homYB, maxlevelYB, homF, true, 0x03);
                }
        }
        catch (...)
        {
                // if failed for two-layer map, then this is a bad error
                if (twolayer)
                        throw;
 
                // otherwise the map must be computed using the new method
                else
                {
                        sout << "*** Homology computation failed. "
                                "Using a careful two-layer method. ***\n";
                        error = true;
                }
        }
        
        // if the levels don't agree, the map must be computed again
        if (!error && ((maxlevelF != maxlevelXA) ||
                (maxlevelXA != maxlevelYB)))
        {
                if (twolayer)
                        throw "Wrong homology levels in the Conley index.";
                sout << "*** Homology levels don't agree. "
                        "Using a careful two-layer method. ***\n";
                error = true;
        }
 
        // show the time used for the computation of the index map
        sbug << "The map computed in " << compTime << ".\n";
        compTime = 0;
 
        // if an error occurred for the flat model, use the two-layer one
        if (error)
        {
                delete [] homXA;
                delete [] homYB;
                delete homF;
                CubSetType emptyset;
                X = emptyset;
                A = emptyset;
                Y = emptyset;
                B = emptyset;
                MvMapType emptymap;
                F = emptymap;
                return this -> compute (true);
        }
 
        toplevel = dim; /* previously I had here "dim + 1" --- why? */
        for (int level = 0; level <= toplevel; ++ level)
        {
                std::vector<euclidom> coefficients;
                if (level <= maxlevelXA)
                {
                        for (int i = 0; i < homXA [level]. size (); ++ i)
                                coefficients. push_back
                                        (homXA [level]. coef (i));
                }
                coef. push_back (coefficients);
        }
        if (indmap)
                delete [] indmap;
        if (toplevel >= 0)
                indmap = new mmatrix<euclidom> [toplevel + 1];
        else
        {
                indmap = 0;
                nontrivialindex = 0;
        }
        for (int i = 0; i <= maxlevelXA; ++ i)
                indmap [i] = (*homF) [i];
 
        // DEBUG !!!
        if (false && homF)
        {
                debugfilename [3] = 'h';
                std::ofstream f (debugfilename);
                f << *homF;
        }
 
        delete [] homXA;
        delete [] homYB;
        delete homF;
        return 0;
} /* ConleyIndex::compute */
 
template <class IndexPair, class euclidom>
inline int ConleyIndex<IndexPair,euclidom>::reduce (bool shrink)
{
        // remove those generators whose image is zero or which do not
        // appear in any image, i.e., their preimage is zero
        const int twice_the_level = (toplevel + 1) << 1;
        for (int level = 0; level < twice_the_level; ++ level)
        {
                // retrieve the right matrix and determine: column or row?
                bool columns = (level > toplevel);
                chomp::homology::mmatrix<euclidom> &matr = indmap [columns ?
                        (level - toplevel - 1) : level];
 
                // get the lengths of columns/rows of the matrix
                int n = columns ? matr. getncols () : matr. getnrows ();
                if (n <= 0)
                        continue;
                int *lengths = new int [n];
                for (int i = 0; i < n; ++ i)
                        lengths [i] = columns ? matr. getcol (i). size () :
                                matr. getrow (i). size ();
 
                // remove from each column these elements which have
                // zero columns themselves
                for (int i = 0; i < n; ++ i)
                {
                        int length = lengths [i];
                        if (!length)
                                continue;
                        chomp::homology::chain<euclidom> col_row = columns ?
                                matr. getcol (i) : matr. getrow (i);
                        for (int j = 0; j < length; ++ j)
                        {
                                int row_col = col_row. num (j);
                                if (lengths [row_col] > 0)
                                        continue;
                                matr. add (columns ? row_col : i,
                                        columns ? i : row_col,
                                        -col_row. coef (j));
                        }
                }
 
                // forget the column lengths
                delete [] lengths;
        }
 
        // shrink the matrices and generators if requested to
        for (int level = toplevel; shrink && (level >= 0); -- level)
        {
                // retrieve the references to this level's data
                chomp::homology::mmatrix<euclidom> &matr = indmap [level];
                std::vector<euclidom> &coefs = coef [level];
 
                // if the matrix is empty
                if (!matr. getncols () || !coefs. size ())
                {
                        coefs. clear ();
                        chomp::homology::mmatrix<euclidom> emptymatr;
                        matr = emptymatr;
                //      if (level == toplevel)
                //              -- toplevel;
                        continue;
                }
 
                // select the numbers of generators/columns to keep
                std::vector<int> selected;
                int ncols = matr. getncols ();
                for (int i = 0; i < ncols; ++ i)
                        if (matr. getcol (i). size () >= 1)
                                selected. push_back (i);
                int length = selected. size ();
 
                // if the index at this level is trivial, skip the rest
                if (!length)
                {
                        coefs. clear ();
                        chomp::homology::mmatrix<euclidom> emptymatr;
                        matr = emptymatr;
                //      if (level == toplevel)
                //              -- toplevel;
                        continue;
                }
 
                // create a new matrix
                chomp::homology::mmatrix<euclidom> newmatr;
                newmatr. define (length, length);
                for (int i = 0; i < length; ++ i)
                {
                        for (int j = 0; j < length; ++ j)
                        {
                                newmatr. add (i, j, matr. get (selected [i],
                                        selected [j]));
                        }
                }
                matr = newmatr;
 
                // create a new sequence of coefficients
                std::vector<euclidom> newcoefs (length);
                for (int i = 0; i < length; ++ i)
                        newcoefs [i] = coefs [selected [i]];
                coefs = newcoefs;
        }
 
        // change the sign of negative matrix entries where possible
        // unless 1 == -1 in the Euclidean domain in question
        euclidom e1, minus1;
        e1 = 1;
        minus1 = -e1;
        if (e1 != minus1)
        {
                for (int level = 0; level <= toplevel; ++ level)
                {
                        // retrieve the reference to the matrix at this level
                        chomp::homology::mmatrix<euclidom> &matr =
                                indmap [level];
 
                        // find single -1's and change their sign
                        int n = matr. getncols ();
                        for (int i = 0; i < n; ++ i)
                        {
                                const chomp::homology::chain<euclidom> &column =
                                        matr. getcol (i);
                                if (column. size () != 1)
                                        continue;
                                if (column. coef (0) != minus1)
                                        continue;
                                if (column. num (0) <= i)
                                        continue;
                                matr. multiplyrow (i, minus1);
                                matr. multiplycol (i, minus1);
                        }
                }
        }
 
        return 0;
} /* ConleyIndex::reduce */
 
template <class IndexPair, class euclidom>
inline void ConleyIndex<IndexPair,euclidom>::truncate (int newlevel)
{
        if ((newlevel >= -1) && (newlevel <= toplevel))
        {
                toplevel = newlevel;
                coef. erase (coef. begin () + toplevel + 1, coef. end ());
        }
        return;
} /* ConleyIndex<IndexPair,euclidom>::truncate */
 
template <class IndexPair, class euclidom>
inline int ConleyIndex<IndexPair,euclidom>::eigenvalues (int level,
        std::vector<double> &Re, std::vector<double> &Im,
        bool nonzero, bool sorted) const
{
        using chomp::homology::auto_array;
 
        // if the level is out of scope, return the error code
        if ((level < 0) || (level > toplevel))
                return -1;
 
        // retrieve the matrix and its size
        const chomp::homology::mmatrix<euclidom> &m = indmap [level];
        int fullsize = m. getncols ();
        if (fullsize <= 0)
                return 0;
 
        // retrieve the indices not related to torsion
        int size = fullsize;
        const std::vector<euclidom> &coefs = coef [level];
        std::vector<int> indices (size);
        int current = 0;
        for (int i = 0; i < fullsize; ++ i)
        {
                if (coefs [i]. delta () > 1)
                {
                        indices [i] = -1;
                        -- size;
                }
                else
                {
                        indices [i] = current;
                        ++ current;
                }
        }
        if (size <= 0)
                return 0;
 
        // prepare data for the LAPACK function to compute the eigenvalues
        char N = 'N';
        int32 dimension = size;
        auto_array<double> APtr (new double [size * size]);
        double *A = APtr. get ();
        auto_array<double> rePtr (new double [size]);
        double *re = rePtr. get ();
        auto_array<double> imPtr (new double [size]);
        double *im = imPtr. get ();
        int32 worksize = 3 * size + (2 * size + 128);
        auto_array<double> workPtr (new double [worksize]);
        double *work = workPtr. get ();
        int32 result = 0;
 
        // translate the chain matrix to the matrix of doubles
        double *ptr = A;
        for (int i = 0; i < size; ++ i)
        {
                for (int j = 0; j < size; ++ j)
                {
                        *(ptr ++) = 0;
                }
        }
        for (int col = 0; col < fullsize; ++ col)
        {
                int icol = indices [col];
                if (icol < 0)
                        continue;
                const chomp::homology::chain<euclidom> c = m. getcol (col);
                int len = c. size ();
                for (int i = 0; i < len; ++ i)
                {
                        int irow = indices [c. num (i)];
                        if (irow < 0)
                                continue;
                        euclidom e = c. coef (i);
                        int n = e. delta ();
                        euclidom e1;
                        e1 = n;
                        if (e1 != e)
                                n = -n;
                        A [size * irow + icol] = n;
                }
        }
 
        // call the LAPACK function to compute the eigenvalues of the map
        // (note: the eigenvalues of A^T are the same as of the eigenvalues
        // of A, so I don't care that the matrices are reprsented
        // in a different way in Fortran)
        dgeev_ (&N, &N, &dimension, A, &dimension, re, im,
                0, &dimension, 0, &dimension, work, &worksize, &result);
 
        // throw an error message in case of failure
        if (result < 0)
                throw "Conley Index eigenvalues: LAPACK error.";
        else if (result > 0)
                throw "The QR algorithm failed in the Conley Index "
                        "Eigenvalues computation.";
 
        // gather the output data
        std::vector<ppComplex> values;
        const double epsilon = 0.001;
        const double threshold = 0.01;
        for (int i = 0; i < size; ++ i)
        {
                re [i] = round (re [i], epsilon, threshold);
                im [i] = round (im [i], epsilon, threshold);
                if (nonzero && !re [i] && !im [i])
                        continue;
                values. push_back (ppComplex (re [i], im [i]));
        }
 
        // make a note of any nontrivial eigenvalues
        if (values. size ())
                nontrivialindex = 1;
 
        // sort the eigenvalues if requested to
        if (sorted)
                std::sort (values. begin (), values. end ());
 
        // write the results to the provided data structures
        for (std::vector<ppComplex>::const_iterator it = values. begin ();
                it != values. end (); ++ it)
        {
                Re. push_back (it -> re);
                Im. push_back (it -> im);
        }
 
        return values. size ();
} /* ConleyIndex<IndexPair,euclidom>::eigenvalues */
 
template <class IndexPair, class euclidom>
inline bool ConleyIndex<IndexPair,euclidom>::trivial () const
{
        // if the answer has already been computed before, use it now
        if (nontrivialindex >= 0)
                return !nontrivialindex;
 
        for (int level = 0; level <= toplevel; ++ level)
        {
                std::vector<double> Re, Im;
                int n = eigenvalues (level, Re, Im, true, false);
                if (n)
                {
                        nontrivialindex = 1;
                        return false;
                }
        //      if (coef [level]. size () && (indmap [level]. getncols () > 0))
        //              return false;
        }
        nontrivialindex = 0;
        return true;
} /* ConleyIndex::trivial */
 
template <class IndexPair, class euclidom>
inline std::string ConleyIndex<IndexPair,euclidom>::EigenString (int n) const
{
        // make sure the level is correct
        if (!indmap || (n < 0) || (n > toplevel))
                return std::string ();
 
        // prepare the output string stream
        std::ostringstream out;
 
        // compute the nonzero eigenvalues and sort them
        std::vector<double> Re, Im;
        int count = eigenvalues (n, Re, Im);
        if (count > 0)
                nontrivialindex = 1;
 
        // output the eigenvalues
        out << "Eigenvalues " << n << ": (";
        if (count >= 0)
        {
                for (int i = 0; i < count; ++ i)
                {
                        if (i)
                                out << ", ";
                        if (Re [i])
                        {
                                out << Re [i];
                                if (Im [i] > 0)
                                        out << "+";
                        }
                        if (Im [i])
                                out << Im [i] << "i";
                }
        }
        else
                out << "ERROR";
        out << ").";
 
        return out. str ();
} /* ConleyIndex::EigenString */
 
template <class IndexPair, class euclidom>
inline void ConleyIndex<IndexPair,euclidom>::define (const int *betti,
        const chomp::homology::mmatrix<euclidom> *matr, int _toplevel)
{
        // reset the potentially remembered nontriviality
        nontrivialindex = -1;
 
        // set the top level of homology
        toplevel = _toplevel;
 
        // create the coefficients vector
        coef. clear ();
        euclidom e;
        e = 1;
        for (int i = 0; i <= toplevel; ++ i)
        {
                std::vector<euclidom> b;
                for (int j = betti [i]; j > 0; -- j)
                        b. push_back (e);
                coef. push_back (b);
        }
 
        // set the index maps
        if (indmap)
        {
                delete [] indmap;
                indmap = 0;
        }
        if (toplevel >= 0)
        {
                indmap = new chomp::homology::mmatrix<euclidom>
                        [toplevel + 1];
        }
        for (int level = 0; level <= toplevel; ++ level)
        {
                indmap [level] = matr [level];
        }
        return;
} /* ConleyIndex::define */
 
template <class IndexPair, class euclidom>
inline void ConleyIndex<IndexPair,euclidom>::define (const int *betti,
        int _toplevel)
{
        // create an array of the identity maps
        chomp::homology::mmatrix<euclidom> *matr = (_toplevel < 0) ? 0 :
                new chomp::homology::mmatrix<euclidom> [_toplevel + 1];
        for (int level = 0; level <= _toplevel; ++ level)
        {
                if (betti [level] > 0)
                        matr [level]. identity (betti [level]);
        }
 
        // define the index
        this -> define (betti, matr, _toplevel);
 
        // clean up and finalize
        if (matr)
                delete [] matr;
        return;
} /* ConleyIndex::define */
 
// --------------------------------------------------
 
template <class IndexPair, class euclidom>
std::ostream &operator << (std::ostream &out,
        const ConleyIndex<IndexPair,euclidom> &c)
{
        // write the top level
        out << c. toplevel << " ";
        if (c. toplevel < 0)
                return out;
 
        // write the coefficients
        for (int i = 0; i <= c. toplevel; ++ i)
        {
                int size = c. coef [i]. size ();
                out << size << " ";
                for (int j = 0; j < size; ++ j)
                        out << c. coef [i] [j] << " ";
        }
 
        // write the maps
        for (int i = 0; i <= c. toplevel; ++ i)
        {
                const chomp::homology::mmatrix<euclidom> &m = c. indmap [i];
                out << m. getnrows () << " " << m. getncols () << " ";
                for (int j = 0; j < m. getncols (); ++ j)
                {
                        const chomp::homology::chain<euclidom> &col =
                                m. getcol (j);
                        out << col. size () << " ";
                        for (int k = 0; k < col. size (); ++ k)
                        {
                                out << col. num (k) << " ";
                                out << col. coef (k) << " ";
                        }
                }
        }
 
        return out;
} /* operator << */
 
template <class IndexPair, class euclidom>
std::istream &operator >> (std::istream &in,
        ConleyIndex<IndexPair,euclidom> &c)
{
        // reset the data
        c. nontrivialindex = -1;
        if (c. indmap)
        {
                delete [] c. indmap;
                c. indmap = 0;
        }
        c. coef. clear ();
 
        // read the top level
        in >> c. toplevel;
        if (c. toplevel < 0)
                return in;
 
        // read the coefficients
        for (int i = 0; i <= c. toplevel; ++ i)
        {
                int size = -1;
                in >> size;
                if (size < 0)
                        throw "Cannot decode the Conley index from input.";
                std::vector<euclidom> coef;
                for (int j = 0; j < size; ++ j)
                {
                        euclidom e;
                        e = 0;
                        in >> e;
                        if (e. delta () <= 0)
                                throw "Wrong coefficient in the "
                                        "Conley index from input.";
                        coef. push_back (e);
                }
                c. coef. push_back (coef);
        }
 
        // read the maps
        if (c. toplevel >= 0)
                c. indmap = new chomp::homology::mmatrix<euclidom>
                        [c. toplevel + 1];
        for (int i = 0; i <= c. toplevel; ++ i)
        {
                chomp::homology::mmatrix<euclidom> &m = c. indmap [i];
                int nrows = 0, ncols = 0;
                in >> nrows >> ncols;
                if ((nrows < 0) || (ncols < 0))
                        throw "Wrong matrix size in the Conley index "
                                "from input.";
                if ((nrows > 0) && (ncols > 0))
                        m. define (nrows, ncols);
                for (int j = 0; j < ncols; ++ j)
                {
                        int length = -1;
                        in >> length;
                        if (length < 0)
                                throw "Negative column length in the "
                                        "Conley index from input.";
                        for (int k = 0; k < length; ++ k)
                        {
                                int num = -1;
                                in >> num;
                                if ((num < 0) || (num >= nrows))
                                        throw "Wrong row number in the "
                                                "Conley index from input.";
                                euclidom e;
                                e = 0;
                                in >> e;
                                if (e. delta () <= 0)
                                        throw "Wrong matrix entry in the "
                                                "Conley index from input.";
                                m. add (num, j, e);
                        }
                }
        }
 
        return in;
} /* operator >> */
 
/// Writes the Conley index information in human-readable form
/// to an output stream.
template <class IndexPair, class euclidom>
std::ostream &writeIndexInfo (std::ostream &out,
        const ConleyIndex<IndexPair,euclidom> &c,
        const char *newline = "\n")
{
        int toplevel = c. dim ();
        out << "Homology: " << c. HomString () << newline;
        for (int n = 0; n <= toplevel; ++ n)
        {
                out << "- Level " << n << ":" << newline;
                out << "Nonzero eigenvalues: " << c. EigenString (n) <<
                        newline;
                out << c. MapString (n, newline);
        }
        return out;
} /* writeIndexInfo */
 
 
#endif // _CMGRAPHS_CONINDEX_H_