---

Matrix.h

Go to the documentation of this file.

/* WARANTY NOTICE AND COPYRIGHT
This program 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
of the License, or (at your option) any later version.

This program 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 program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

Copyright (C) Michael J. Meyer

matmjm@mindspring.com
spyqqqdia@yahoo.com

*/

#ifndef martingale_matrix_h
#define martingale_matrix_h


#include <string>
#include <sstream>
#include <iostream>
#include <algorithm>                      // min,max
#include <cmath>
#include "TypedefsMacros.h"
//#include "tnt/tnt_array1d.h"
//#include "tnt/tnt_array2d.h"
#include "jama/jama_eig.h"
//#include "tnt_array1d_utils.h"
//#include "tnt_array2d_utils.h"
#include "Array.h"
//#include "Utils.h"                 

/*
 * DataStructures.h
 *
 * Created on February 19, 2003, 6:00 PM
 */

MTGL_BEGIN_NAMESPACE(Martingale)




/***************************************************************************************

                      Useful global functions

***************************************************************************************/




template<typename S>
S relativeError(S exact, S approx, S epsilon)
{
        S d=exact;
        if((-epsilon<d)&&(d<epsilon)) d=epsilon;
 
        return 100*(approx-exact)/d;
} 




/***************************************************************************************

                      C-Style vectors

***************************************************************************************/

// forward declarations
template<typename S,typename MatrixType> class Matrix;


template<class S>
class Vector {

protected:

    int dim;         // dimension
        int b;           // index base: indices i=b,b+1,...,b+dim-1
    S* dptr;         // pointer to data

public:
        
// ACCESSORS
        
   int getDimension() const { return dim; }
   void setDimension(int d) { dim=d; }
   int getIndexBase() const { return b; }
   void setIndexBase(int base) { b=base; }
   S* getData() const { return dptr; }
 
      
   S& operator[](int i)
   {
           #ifdef SUBSCRIPT_CHECK
              SubscriptCheck::checkSubscript(i,b,dim,"Vector");
           #endif          
           return dptr[i-b]; 
    }
   
   const S& operator[](int i) const 
   { 
           #ifdef SUBSCRIPT_CHECK
              SubscriptCheck::checkSubscript(i,b,dim,"Vector");
           #endif         
           return dptr[i-b]; 
   }
   
   
   // CONSTRUCTION - DESTRUCTION
   ~Vector(){ delete[] dptr; }
   
   explicit Vector(int d, int base=0) : 
   dim(d), b(base), dptr(new S[d]) 
   { 
           for(int i=0;i<dim;i++) dptr[i]=0; 
   }

   
   // shallow copy not good enough if we return from function
   Vector(const Vector& v) : 
   dim(v.getDimension()), b(v.getIndexBase()), dptr(new S[dim])
   {
           S* vdptr=v.getData(); 
           for(int i=0;i<dim;i++) dptr[i]=vdptr[i];
   }
   
   
   Vector(int n, S& v) : 
   dim(n), b(0), dptr(new S[n])
   {
            for(int i=0;i<dim;i++) dptr[i]=v;
   }

   
   Vector(int n, S* a) : dim(n), b(0), dptr(new S[n])
   {
            for(int i=0;i<dim;i++) dptr[i]=a[i];
   }
   
   
   template<int n>
   Vector(S a[n]) : dim(n), dptr(new S[n])
   {
            for(int i=0;i<dim;i++) dptr[i]=a[i];
   }

   
   Vector& operator=(const Vector& x)
   {
       dim=x.getDimension(); 
           b=x.getIndexBase();
           S* xdptr=x.getData();
           for(int i=0;i<dim;i++) dptr[i]=xdptr[i];
           return *this;
   } 
   
   
   
// ALGEBRAIC OPERATIONS NEEDED BY RANDOMOBJECT
   
   Vector& switchComponents(int i, int j)
   { 
       S tmp=dptr[i-b]; dptr[i-b]=dptr[j-b]; dptr[j-b]=tmp; 
   }
   
   
   Vector& operator +=(const Vector& x)
   { 
           S* xdptr=x.getData();
           for(int i=0;i<dim;i++) dptr[i]+=xdptr[i];
           return *this;
   }
   
   Vector& operator -=(const Vector& x)
   { 
           S* xdptr=x.getData();
           for(int i=0;i<dim;i++) dptr[i]-=xdptr[i];
           return *this;
   }
   
   Vector& operator /=(int N)
   { 
           Real f=1.0/N;
           for(int i=0;i<dim;i++) dptr[i]*=f;
           return *this;
   }
   
   
   Vector& operator *=(const S& lambda)
   { 
           for(int i=0;i<dim;i++) dptr[i]*=lambda;
           return *this;
   }

   
   template<typename MatrixBaseType>
   Vector& operator *=(const Matrix<S,MatrixBaseType>& A)
   { 
           int d=A.rows();        // new dimension
           S* Ax=new Real[d];         // new memory
           for(int i=0;i<d;i++){
                   
                   S sum=0.0;
                   for(int j=A.rowBegin(i);j<A.rowEnd(i);j++) sum+=A.entry(i,j)*dptr[j];
                   Ax[i]=sum;
           }

       // free the old memory and reset
           delete[] dptr; dptr=Ax; dim=d;
           return *this;
   }   
   
   
// MEAN, EUCLIDEAN NORM
   
   S mean() const
   {  
           S sum=0; 
           for(int i=0;i<dim;i++) sum+=dptr[i];
           return sum/dim;
   }
   
   Real norm() const
   {  
           Real sum=0, absx_i; 
           for(int i=0;i<dim;i++){ absx_i=fabs(dptr[i]); sum+=absx_i*absx_i; }
           return sqrt(sum);
   }
   
   

// TEST FOR EQUALITY


void testEquals(const Vector<S>& v, S precision, S epsilon, string test) const
{
        int vdim=v.getDimension(), vbase=v.getIndexBase(); 
        if((vdim!=dim)||(vbase!=b)){
                
                std::cout << "\n"+test+": failed, dimensions or index range not equal.";
                return;
        }
                
        for(int i=0;i<dim;i++)
        { 
                S err=relativeError(dptr[i],v[i+b],epsilon);
                if(err>precision){
                        
                     std::cout << "\n"+test+": failed, "
                               << "component v["<<i<<"], relative error " << err;
                     return;
                }
        } // end for i
        
    std::cout << "\n"+test+": passed.";
} // end testEquals


// PRINTING

std::ostream&  printSelf(std::ostream& os) const
{
        os << "\n\nVector of dimension " << dim << ":" << endl;
    for(int i=0;i<dim-1;i++) os << dptr[i] << ", ";
    os << dptr[dim-1];
    return os << endl;
} // end operator <<
        
              
}; // end Vector


template<class S>
std::ostream& operator << (std::ostream& os, const Vector<S>& v)
{
        return v.printSelf(os);
} 


typedef Vector<Real> RealVector;
        


/***************************************************************************************
 *
 *                             MATRIX CLASSES
 *
***************************************************************************************/

// IMPLEMENTATION NOTES: we use the Barton-Nachman trick.
// The type of matrix (triangular,symmetric,...) is determined by the base class
// from which the matrix class derives and this base class is a template parameter
// for the derived class template.
//
// Matrices are stored in row major order and usually traversed row by row.
// For each i we need to know where nonzero entries begin and end in this row.
// The functions rowBegin(i) and rowEnd(i) provide this information.
// The function rowEnd(i) returns the index one past the last nonzero entry
// in row i. The functions colBegin(i) and colEnd(i) do the same for
// the columns. The information returned by these functions is type specific
// (upper triangular, lower triangular,...) and not specific to each matrix 
// object and allows optimizations depending on the matrix type (matrix
// products).



/***************************************************************************************
 *
 *                      MATRIX BASE TYPES
 *
***************************************************************************************/

// forward declaration
template<typename S> class LowerTriangular;


// IMPLEMENTATION NOTES: the matrix base classes define the behaviour for the 
// various special matrix classes. They are put at the base of the derived class
// Matrix using the Barton-Nachman trick (base class is a template parameter for
// the derived template class). 
// 
// Consequently many functions will seem to have unnecessary general parameter
// signatures (the most general applying to any matrix type). This is necessary
// as the derived class Matrix has to be written in the most general case and
// has to call functions without any assumptions of the type of the base class.


template<typename S>
class UpperTriangular {
        
protected:
        
        int dim;      // number of rows = number of cols            
        S** dptr;     // pointer to data
        
public:
        
        typedef LowerTriangular<S> TransposeType;
        
int rows() const { return dim; }
int cols() const { return dim; }
void setCols(int) {   }

S** allocate(int nRows, int nCols)
{
        if(nRows!=nCols)
        { cout<<"\n\nUpper triangular matrix must be square. Terminating."; exit(1); }
        S** data=new S*[nRows];
        for(int i=0;i<nRows;i++){
                
        data[i]=new S[nRows-i];
        for(int j=i;j<nRows;j++)data[i][j-i]=0;
    }
        return data;
}


UpperTriangular(int nRows, int nCols) : dim(nRows), dptr(allocate(nRows,nCols)) {  } 
        
~UpperTriangular(){  for(int i=0;i<dim;i++) delete[] dptr[i]; delete[] dptr; }
        
        

S& entry(int i, int j){ return dptr[i][j-i]; }
const S& entry(int i, int j) const { return dptr[i][j-i]; }

                
int rowBegin(int i) const { return i; }

int rowEnd(int i) const { return dim; }

int colBegin(int j) const { return 0; }

int colEnd(int j) const { return j+1; }

        
}; // end UpperTriangular




template<typename S>
class LowerTriangular{
                
protected:
        
        int dim;      // number of rows = number of cols            
        S** dptr;     // pointer to data
        
public:
        
        typedef UpperTriangular<S> TransposeType;
        
int rows() const { return dim; }
int cols() const { return dim; }
void setCols(int) {   }


S** allocate(int nRows, int nCols)
{
        if(nRows!=nCols)
        { cout<<"\n\nLower triangular matrix must be square. Terminating."; exit(1); }
        S** data=new S*[nRows];
        for(int i=0;i<nRows;i++){
                
        data[i]=new S[i+1];
        for(int j=0;j<=i;j++)data[i][j]=0;
   }
   return data;
}

                
LowerTriangular(int nRows, int nCols) : dim(nRows), dptr(allocate(nRows,nCols)) {  } 

~LowerTriangular(){  for(int i=0;i<dim;i++) delete[] dptr[i]; delete[] dptr; }
        
        
S& entry(int i, int j){ return dptr[i][j]; }
const S& entry(int i, int j) const { return dptr[i][j]; }
        
                
int rowBegin(int i) const { return 0; }

int rowEnd(int i) const { return i+1; }

int colBegin(int j) const { return j; }

int colEnd(int j) const { return dim; }

                
}; // end LowerTriangular




template<typename S>
class Rectangular {
        
protected:
        
        int rows_,        // number of rows 
            cols_;        // number of columns
        S** dptr;         // pointer to data

        
public:
        
        typedef Rectangular<S> TransposeType;
        
int rows() const { return rows_; }
int cols() const { return cols_; }
void setCols(int nCols) { cols_=nCols;  }


S** allocate(int nRows, int nCols)
{
        S** data=new S*[nRows];
        for(int i=0;i<nRows;i++){
                
        data[i]=new S[nCols];
        for(int j=0;j<nCols;j++)data[i][j]=0;
   }
   return data;
}
        
                
Rectangular(int nRows, int nCols) : 
rows_(nRows), cols_(nCols), dptr(allocate(nRows,nCols))
{    } 

~Rectangular(){  for(int i=0;i<rows_;i++) delete[] dptr[i]; delete[] dptr; }
        
        
S& entry(int i, int j){ return dptr[i][j]; }
const S& entry(int i, int j) const { return dptr[i][j]; }

        
int rowBegin(int i) const { return 0; }

int rowEnd(int i) const { return cols_; }

int colBegin(int j) const { return 0; }

int colEnd(int j) const { return rows_; }

        
}; // end Rectangular



/***************************************************************************************
 *
 *                      MATRIX PRODUCT RETURN TYPE
 *
***************************************************************************************/


template<typename S, typename MatrixType1, typename MatrixType2>
struct ProductType { typedef Rectangular<S> type; };

template<typename S, typename MatrixType>
struct ProductType<S,MatrixType,MatrixType> { typedef MatrixType type; };
        



/***************************************************************************************
 *
 *                             MATRIX
 *
***************************************************************************************/


template<typename S, typename MatrixBaseType>
class Matrix : public MatrixBaseType {

protected:

        int a;           // row index base
        int b;           // column index base
        
        static S cache[SMALL][SMALL];          // workspace for small matrix optimization.


public:
         
        typedef typename MatrixBaseType::TransposeType TransposeType;


// FREE MEMORY
        
void deallocate(){ for(int i=0;i<rows();i++) delete[] dptr[i]; delete[] dptr; }
                
// ACCESSORS

S** getData() const { return dptr; }

int getRowIndexBase() const { return a; }
void setRowIndexBase(int r){ a=r; }

int getColIndexBase() const { return b; }
void setColIndexBase(int c){ b=c; }


S& operator()(int i, int j)
{ 
        #ifdef SUBSCRIPT_CHECK
          SubscriptCheck::checkSubscript(i,j,a,b,rows(),cols(),"Matrix");
        #endif          
        return entry(i-a,j-b);
}


const S& operator()(int i, int j) const 
{ 
        #ifdef SUBSCRIPT_CHECK
          SubscriptCheck::checkSubscript(i,j,a,b,rows(),cols(),"Matrix");
        #endif          
        return entry(i-a,j-b);
}       
        

// EQUALITY CHECK

void testEquals(const Matrix& B, S precision, S epsilon, string test) const
{
        int a_B=B.getRowIndexBase(), b_B=B.getColIndexBase();
        if((B.rows()!=rows())||(B.cols()!=cols())){
                
                std::cout << "\nMatrix::testEquals(): "+test+": failed, dimensions not equal.";
                return;
        }
        
        if((a_B!=a)||(b_B!=b)){
                
                std::cout << "\nMatrix::testEquals(): "+test+": failed, index bases not equal.";
                return;
        }
        
        for(int i=0;i<rows();i++)
        for(int j=rowBegin(i);j<rowEnd(i);j++) { 
                
                    S err=relativeError(entry(i,j),B.entry(i,j),epsilon);
                    if((err>precision)||(-err>precision)){
                        
                            std::cout << "\n"+test+": failed, "
                                      << "component B["<<i<<"]["<<j<<"], relative error " << err;
                        return;
                    } // endif
        } // end for i

    std::cout << "\n"+test+": passed.";
} // end testEquals


// CONSTRUCTION - DESTRUCTION


Matrix(int dim, int c=0) : MatrixBaseType(dim,dim), a(c), b(c) {    }


Matrix(int nRows, int nCols, int row_base, int col_base) : 
MatrixBaseType(nRows,nCols),
a(row_base), b(col_base) 
{    }



Matrix(const Matrix& A) : MatrixBaseType(A.rows(),A.cols()),
a(A.getRowIndexBase()), b(A.getColIndexBase()) 
{
        for(int i=0;i<rows();i++)
        for(int j=rowBegin(i);j<rowEnd(i);j++) entry(i,j)=A.entry(i,j);
}


template<int dim>
Matrix(S A[dim][dim], int row_base=0, int col_base=0) : 
MatrixBaseType(dim,dim),
a(row_base), b(col_base) 
{
        for(int i=0;i<rows();i++)
        for(int j=rowBegin(i);j<rowEnd(i);j++) entry(i,j)=A[i][j]; 
}


Matrix& operator =(const Matrix& B)
{
    if(this=&B) return *this;
        if((rows()!=B.rows())||(cols()!=B.cols())){
                
                cout << "\n\nMatrix assignement: matrix sizes different. Terminating.";
                exit(1);
        }

        a=B.getRowIndexBase(); b=B.getColIndexBase();
        for(int i=0;i<rows();i++)
        for(int j=rowBegin(i);j<rowEnd(i);j++) entry(i,j)=B.entry(i,j);
        
        return *this;
}



// NORMS

Real norm() const
{
        Real sum=0.0, absij;
        for(int i=0;i<rows();i++)
        for(int j=rowBegin(i);j<rowEnd(i);j++) { absij=std::abs(entry(i,j)); sum+=absij*absij; }
        
        return sqrt(sum);
} // end norm



Real rowNorm(int i) const
{
        Real sum=0, absij;
        for(int j=rowBegin(i);j<rowEnd(i);j++) { absij=std::abs(entry(i,j)); sum+=absij*absij; }
        
        return sqrt(sum);
} // end norm



Real colNorm(int j)
{
        Real sum=0, absij;
        for(int i=colBegin(j);i<colEnd(j);i++) { absij=std::abs(entry(i,j)); sum+=absij*absij; }
        
        return sqrt(sum);
} // end norm
        


Matrix& scaleRow(int i, Real f)
{
        for(int j=rowBegin(i);j<rowEnd(i);j++) entry(i-a,j)*=f; 
        return *this;
} // end scaleRow


Matrix& scaleCol(int j, Real f)
{
        for(int i=colBegin(j);i<colEnd(j);i++) entry(i,j-b)*=f;         
        return *this;
} // end scaleRow


        

// ALGEBRAIC OPERATIONS

S quadraticForm(const Vector<S>& x) const
{
        int dim_x=x.getDimension(), 
            b_x=x.getIndexBase();
        if((rows()!=dim_x)||(cols()!=dim_x)){
                
                cout << "\n\nMatrix::quadraticForm(): dimensional mismatch. Terminating.";
                exit(1);
        }
        S d=0, u=0, x_i,x_j;
        // diagonal
        for(int i=0;i<rows();i++) { x_i=x[i+b_x]; d+=x_i*x_i*entry(i,i); }
        // above the diagonal
        for(int i=0;i<rows();i++)
        for(int j=i+1;j<cols();j++) { x_i=x[i+b_x]; x_j=x[j+b_x]; u+=x_i*x_j*entry(i,j); }

        return d+2*u;
}



Matrix<S,TransposeType>& transpose() const
{
        Matrix<S,TransposeType>* 
        trnsps=new Matrix<S,TransposeType>(cols(),rows(),b,a);
    for(int i=0;i<rows();i++)
        for(int j=rowBegin(i);j<rowEnd(i);j++) trnsps->entry(j,i)=entry(i,j);
        
        return *trnsps;
}
        

Matrix& operator += (const Matrix& B)
{
        if((rows()!=B.rows())||(cols()!=B.cols())){
                
                cout << "\n\nMatrix +=: matrix sizes different. Terminating.";
                exit(1);
        }
        for(int i=0;i<rows();i++)
        for(int j=rowBegin(i);j<rowEnd(i);j++) entry(i,j)+=B.entry(i,j);
                
        return *this;
} // end operator +=


Matrix& operator *= (Real f)
{
        for(int i=0;i<rows();i++)
        for(int j=rowBegin(i);j<rowEnd(i);j++) entry(i,j)*=f;
                
        return *this;
} // end operator +=


// product will have the same type as this.
Matrix& operator *= (const Matrix& B)
{
        if(cols()!=B.rows()) {
                
                cout << "\n\nMatrix *=: matrix sizes not compatible. Terminating.";
                exit(1);
        }
        
        // if both matrices are small and square use the workspace, no memory allocation.
        // matrix size does not change and so the product matrix can use the old memory.
        if((rows()<SMALL)&&(rows()==cols())&&(cols()==B.cols())){
        
            for(int i=0;i<rows();i++)
            for(int j=rowBegin(i);j<rowEnd(i);j++){
                
            S sum=0.0;
                    int u=std::max(rowBegin(i),B.colBegin(j)),
                        v=std::min(rowEnd(i),B.colEnd(j));
                    for(int k=u;k<v;k++) sum+=entry(i,k)*B.entry(k,j);
                    cache[i][j]=sum;            
            }
            // copy the entries from the workspace
            for(int i=0;i<rows();i++)
            for(int j=rowBegin(i);j<rowEnd(i);j++) entry(i,j)=cache[i][j];
                return *this;
        } // end if
        
        // otherwise allocate new memory for the product
        Matrix<S,MatrixBaseType>* product = new Matrix<S,MatrixBaseType>(rows(),B.cols(),a,b);
    for(int i=0;i<rows();i++)
        for(int j=product->rowBegin(i);j<product->rowEnd(i);j++){
                
        S sum=0.0;
                int u=std::max(rowBegin(i),B.colBegin(j)),
                    v=std::min(rowEnd(i),B.colEnd(j));
                for(int k=u;k<v;k++) sum+=entry(i,k)*B.entry(k,j);
                product->entry(i,j)=sum;
        }
        deallocate();                // free old memory
        dptr=product->getData();
        setCols(B.cols());
        return *this;
                
} // end operator *=



// product will have the same type as this.
Matrix& operator ^= (const Matrix<S,TransposeType>& B)
{
         if(cols()!=B.cols()){
                   
                   cout << "\n\nMatrix ^=: dimensions not compatible. Terminating.";
               exit(1);
         }

        // if both matrices are small and square use the workspace, no memory allocation.
        // matrix size does not change and so the product matrix can use the memory of dptr
        if((rows()<SMALL)&&(rows()==cols())&&(rows()==B.rows())){
        
            for(int i=0;i<rows();i++)
            for(int j=rowBegin(i);j<rowEnd(i);j++){
                
            Real sum=0.0;
                    int u=std::max(rowBegin(i),B.rowBegin(j)),
                        v=std::min(rowEnd(i),B.rowEnd(j));
                    for(int k=u;k<v;k++) sum+=entry(i,k)*B.entry(j,k);
                    cache[i][j]=sum;            
            }
            // copy the entries from the workspace
            for(int i=0;i<rows();i++)
            for(int j=rowBegin(i);j<rowEnd(i);j++) entry(i,j)=cache[i][j];
                return *this;
        } // end if
        
        // otherwise allocate new memory for the product,
        Matrix<S,MatrixBaseType>* product = new Matrix<S,MatrixBaseType>(rows(),B.rows(),a,b);
    for(int i=0;i<rows();i++)
        for(int j=product->rowBegin(i);j<product->rowEnd(i);j++){
                
        S sum=0.0;
                int u=std::max(rowBegin(i),B.rowBegin(j)),
                    v=std::min(rowEnd(i),B.rowEnd(j));
                for(int k=u;k<v;k++) sum+=entry(i,k)*B.entry(j,k);
                product->entry(i,j)=sum;
        }
        deallocate();                     // free the old memory
        dptr=product->getData();
        setCols(B.rows());
        return *this;
                
} // end operator ^=




// MATRIX FUNCTIONS: some analytic functions f(A), where A=this.


Matrix< S,UpperTriangular<S> >& aat() const
{
    Matrix< S,UpperTriangular<S> >* 
        c = new Matrix< S,UpperTriangular<S> >(rows(),rows(),0,0); // aat
        for(int i=0;i<rows();i++)
        for(int j=i;j<rows();j++){
                
                S sum=0.0;                                      // row_i\cdot row_j
        int u=std::max(rowBegin(i),rowBegin(j)),
                    v=std::min(rowEnd(i),rowEnd(j));
                for(int k=u;k<v;k++) sum+=entry(i,k)*entry(j,k);
                c->entry(i,j)=sum;
        }
        return *c;
} // end aat



Matrix& exp() const
{
    if(rows()!=cols()){
                
                cout << "\n\nMatrix::exp(): matrix not square. Terminating.";
                exit(1);
        }
        
        int n=2;
        Real fn=norm(), f=4.0, g=1.0;
        while((f<fn)){ f*=2.0; n++; }                   // f=2^n
        g=1.0/f; g/=64.0;                               // g=1/2^{n+6}

        // temporary matrices
        Matrix<S,MatrixBaseType> U(rows());             // U=(A/2^{n+6})=g*A 
        Matrix<S,TransposeType>  Ut(rows());            // U', we use ^=U ie. *=U'
        Matrix<S,MatrixBaseType> SUk(rows());           // sum of U^k/k!, k=0,1,..,6

        for(int i=0;i<rows();i++)
        for(int j=rowBegin(i);j<rowEnd(i);j++){ 
                
                U(i,j)=Ut(j,i)=SUk(i,j)=g*entry(i,j);       // U=A/2^{n+6}
                if(i==j) SUk(i,j)+=1.0;                     // I+U
        }
        
        // Suk=I+U+U^2/2!+...+U^8/8! 
        for(int k=2;k<=8;k++){ U^=Ut; g=1.0/k; U*=g; SUk+=U; }
        
        // exp(A)=Suk^m, m=2^{n+6};
    Matrix<S,MatrixBaseType>& F=*(new Matrix<S,MatrixBaseType>(SUk));    
        for(int k=0;k<n+6;k++) F*=F;
        F.setRowIndexBase(a);   
        F.setColIndexBase(b);   
        return F;

} // end exp


// PRINTING

std::ostream& printSelf(std::ostream& os) const
{
        os << "\n\nMatrix, rows: " << rows() << endl << endl;
        for(int i=0;i<rows();i++){
                
                int u=rowBegin(i), v=rowEnd(i); 
                for(int j=0;j<u;j++) os << 0.0 << ", ";
                for(int j=u;j<v-1;j++) os << entry(i,j) << ", ";
                os << entry(i,v-1) << endl;
        }
    return os << endl << endl;
} 


}; // end Matrix




// cache must be initialized outside the class:
template<typename S,typename MatrixBaseType>
S Matrix<S,MatrixBaseType>::cache[SMALL][SMALL] = { };


template<typename S, typename MatrixBaseType>
std::ostream& operator << (std::ostream& os, const Matrix<S,MatrixBaseType>& A) 
{
    return A.printSelf(os);
} // end operator <<





// SOME  MATRIX TYPES

// since there are no templated typedefs:
template<typename S>
class UTRMatrix : public Matrix< S,UpperTriangular<S> >{
        
public:

UTRMatrix(int dim, int a=0) : Matrix< S,UpperTriangular<S> >(dim,dim,a,a) {  }
        
};
        
        
// upper triangular and symmetric real types are more elaborate below.
typedef Matrix< Real,Rectangular<Real> > RealMatrix;
typedef Matrix< Real,LowerTriangular<Real> > LTRRealMatrix;



/***************************************************************************************
 *
 *                             UPPER TRIANGULAR MATRICES
 *
***************************************************************************************/



// DIAGONALIZATION, RANK REDUCED FACTORIZATION, EIGEN ANALYSIS

// forward declarations
class UTRRealMatrix;


JAMA::Eigenvalue<Real>* eigenDecomposition(const UTRRealMatrix& C);


RealMatrix& rank_Reduced_Root(const UTRRealMatrix& C, int r);
                


void factorAnalysis(const UTRRealMatrix& C, int r);
        


void factorizationTest(const UTRRealMatrix& C, int r, string message="");



class UTRRealMatrix : public Matrix< Real,UpperTriangular<Real> > {

public:
        
UTRRealMatrix(int d, int b=0) : Matrix< Real,UpperTriangular<Real> >(d,d,b,b) {  }

        
template<int n>
UTRRealMatrix(Real A[n][n], int b=0) : Matrix< Real,UpperTriangular<Real> >(A,b) {   }
        
        
RealMatrix& symmetricCompletion() const;


UTRRealMatrix& inverse() const;
        

LTRRealMatrix& ltrRoot() const;


UTRRealMatrix& utrRoot() const;


RealMatrix& rankReducedRoot(int r) const;


UTRRealMatrix& exp() const;


RealMatrix& matrixFunction(Real (*f)(Real));


void analyseFactors(int r) const;


void testFactorization(int r, string message="") const;


}; // end UTRRealMatrix




/***************************************************************************************
 *
 *            MATRIX ARRAYS
 *
***************************************************************************************/


class UTRMatrixSequence {
        
        int n;                         // number of matrices
        // array of pointers to const UTRMatrix: *matrix[t] is read only
        Array1D<const UTRRealMatrix*> matrices;    
        
public:
        
        UTRMatrixSequence(int m) : 
        n(m), matrices(m) 
    {  }  
        
        ~UTRMatrixSequence(){ for(int t=0;t<n;t++) delete matrices[t]; }

        
    const UTRRealMatrix& getMatrix(int t) const { return *(matrices[t]); }
        
    void setMatrix(int t, const UTRRealMatrix& matrix) { matrices[t]=&matrix; }
        
}; // end UTRMatrixSequence




class MatrixSequence {
        
        int n;                                    // number of matrices
        Array1D<const RealMatrix*> matrices;      // matrices[t]: pointer to matrix(t)
        
public:
        
        MatrixSequence(int m) : 
        n(m), matrices(m) 
    {  }  
        
        ~MatrixSequence(){ for(int t=0;t<n;t++) delete matrices[t]; }

        
    const RealMatrix& getMatrix(int t) const { return *(matrices[t]); }
        
    void setMatrix(int t, const RealMatrix& matrix) { matrices[t]=&matrix; }
        
}; // endMatrixSequence




MTGL_END_NAMESPACE(Martingale)

#endif

---