UT/UT_Matrix.h

Go to the documentation of this file.

/*
 * PROPRIETARY INFORMATION.  This software is proprietary to
 * Side Effects Software Inc., and is not to be reproduced,
 * transmitted, or disclosed in any way without written permission.
 *
 * Produced by:
 *      David Wong
 *      Side Effects
 *      123 Front Street West, Suite 1401
 *      Toronto, Ontario
 *      Canada   M5J 2M2
 *      416-504-9876
 *
 * NAME:        Matrix of arbitrary size (C++)
 *
 * COMMENTS:    From Numerical Recipes
 *      See Matrix Computations, Gene H. Golub, Charles F. Van Loan
 *      Rewritten without row array by Jeff Lait
 *
 */

#ifndef __UT_Matrix_H__
#define __UT_Matrix_H__

#include "UT_API.h"
#include <stdlib.h>
#include <iostream.h>
#include <SYS/SYS_Types.h>
#include "UT_Assert.h"
#include "UT_Math.h"
#include "UT_Matrix2.h"
#include "UT_Vector.h"

class UT_DMatrix3;
class UT_DMatrix4;

class UT_API UT_Matrix {
public:
     UT_Matrix();
     
     // Input the row and column index range [nrl..nrh][ncl..nch].
     // The "standard" way to set up a matrix of R rows and C columns
     // would be:
     // UT_Matrix               foo(1, R, 1, C);
     
     UT_Matrix(int nrl, int nrh, int ncl, int nch);
     UT_Matrix(int nrl, int nrh, int ncl, int nch, float *array);
     UT_Matrix(int nrl, int ncl, UT_Matrix &mat);
     explicit UT_Matrix(const UT_Matrix &mat);  // deep copy constructor
    ~UT_Matrix();

    // Allocate space.  Resulting memory is not initialized.
    void         init(int nrl, int nrh, int ncl, int nch);

    // Enlarge/reduce matrix to have a new number of rows and
    // columns.  NRL & NCL are unchanged.  Afterwards, the
    // matrix always owns its data!
    void         resize(int nrows, int ncols);

    // Enlarge by one row/column and copy given data in (if given, otherwise
    // filled in with zeroes
    void         appendCol( float *new_col=0 );
    void         appendRow( float *new_row=0 );

    // insert new row of zeros (or nw_row, if given) at the specified row idx
    //
    // this has been tested and it works, we just don't use it anywhere.
    // uncomment it if you need it.
    // void      insertRowAt( int at, float *new_row=0 );

    int          isInit() const { return myMatrix ? 1 : 0; }

    // Steal from another matrix or array
    // Resulting matrix has origin 1, 1.
    void         submatrix(const UT_Matrix &A, 
                            int nrl, int nrh, int ncl, int nch);
    // &(*this)(nrl, ncl) == array.  If array is wider than
    // nch - ncl + 1, stride should be set.
    void         submatrix(float *array, 
                            int nrl, int nrh, int ncl, int nch, int stride=-1);

    // Sets to zeros.  The version with a range will set a sub-block
    // to zero.
    void         zero(int nrl, int nrh, int ncl, int nch);
    void         zero()
                    { zero(myNRL, myNRH, myNCL, myNCH); }

    /// These methods allow one to read out and write into blocks of
    /// the UT_Matrix using our other matrix classes.  This is very useful
    /// when working with block algorithms.
    /// Keep in mind that the usual UT_Matrix? methods are 0 based while
    /// this class is usually 1 based.
    /// Only double matrices are supported as they are the preferred
    /// matrix format.
    void         getSubmatrix2(UT_Matrix2 &mat, int row, int col) const;
    void         setSubmatrix2(int row, int col, const UT_Matrix2 &mat);
    void         addSubmatrix2(int row, int col, const UT_Matrix2 &mat);
    
    void         getSubmatrix3(UT_DMatrix3 &mat, int row, int col) const;
    void         setSubmatrix3(int row, int col, const UT_DMatrix3 &mat);
    void         addSubmatrix3(int row, int col, const UT_DMatrix3 &mat);
    
    void         getSubmatrix4(UT_DMatrix4 &mat, int row, int col) const;
    void         setSubmatrix4(int row, int col, const UT_DMatrix4 &mat);
    void         addSubmatrix4(int row, int col, const UT_DMatrix4 &mat);

    // Make identity:
    void         makeIdentity();

    // Negate the matrix.
    void         negate();

    // Get the row low index.
    int          getNRL() const { return myNRL; }

    // Get the row high index.
    int          getNRH() const { return myNRH; }

    // Get the column low index.
    int          getNCL() const { return myNCL; }

    // Get the column high index.
    int          getNCH() const { return myNCH; }

    // Get number of rows
    int          rows() const { return myNRH - myNRL + 1; }

    // Get number of columns
    int          columns() const { return myNCH - myNCL + 1; }

    // Change the row and column low indices, 
    // and high indices will adjust themselves.
    void         changeNRLAndNCL(int nrl, int ncl);

    // Change the indices in a very shallow manner, i.e. simply change
    // the index without shifting the data. This is highly dangerous, yet
    // very useful if used with care.
    void         setShallowNRL(int nrl) { myNRL = nrl; }
    void         setShallowNRH(int nrh) { myNRH = nrh; }
    void         setShallowNCL(int ncl) { myNCL = ncl; }
    void         setShallowNCH(int nch) { myNCH = nch; }

    // Retrieve data
    float       &operator()(int row, int col)
                 {
                     UT_ASSERT_P(row >= myNRL && row <= myNRH);
                     UT_ASSERT_P(col >= myNCL && col <= myNCH);
                     return myMatrix[row*myStride + col];
                 }
    float        operator()(int row, int col) const
                 {
                     UT_ASSERT_P(row >= myNRL && row <= myNRH);
                     UT_ASSERT_P(col >= myNCL && col <= myNCH);
                     return myMatrix[row*myStride + col];
                 }

    UT_Matrix   &operator=(const UT_Matrix &m);

    // Component wise operators
    UT_Matrix   &operator+=(const UT_Matrix &m);
    UT_Matrix   &operator-=(const UT_Matrix &m);

    // left multiply this matrix[1..m][1..n] by the given matrix
    // A[1..p][1..m] and puts the product in 'result[1..p][1..n]'
    // NOTE: this, A, and result must all have the same lower index
    //       for rows and columns
    void         preMult( const UT_Matrix &A, UT_Matrix &result ) const;

    // right multiply this matrix[1..m][1..n] by the given matrix
    // A[1..n][1..p] and puts the product in 'result[1..m][1..p]'
    // NOTE: this, A, and result must all have the same lower index
    //       for rows and columns
    void         postMult( const UT_Matrix &A, UT_Matrix &result ) const;

    // Vector multiplication, before and after.
    // result = xT * this
    void         preMult(const UT_Vector &x, UT_Vector &result) const;
    // result = this * x
    void         postMult(const UT_Vector &x, UT_Vector &result) const
                                    { multVec(x, result); }

    // Fills the upper-triangular portion of this matrix with A^T . A
    void         upperNormalUpdate(const UT_Matrix &A);

    // Scaled outer product update.
    // this += b * x * yT
    void         outerproductUpdate(float b, 
                            const UT_Vector &x, const UT_Vector &y);

    // Add another matrix scaled appropriately:
    void         addScaledMatrix(const UT_Matrix &A, float scale);
    // this = A*scale
    void         setAndScale(const UT_Matrix &A, float scale);

    // Givens functions: do a pre or post multiplication with
    // the givens matrix.
    void         preMultGivensInPlace(float c, float s);
    void         postMultGivensInPlace(float c, float s);
    
    void         multVec(const UT_Vector &x, UT_Vector &result) const;

    // Calculate the L2-distance (Euclidean) between two rows for columns [cl..ch].
    float        rowsL2dist(int r1, int r2, int cl=-1, int ch=-1) const;

    // Norm calculations.
    float        normFrobenius() const;
    float        norm1() const;
    float        normInfinite() const;

    // place this matrix's transpose in the given result matrix
    void         transpose( UT_Matrix &result ) const;

    // Returns true if matrix is symmetric within given tolerance
    bool         isSymmetric(fpreal64 tolerance = UT_FTOLERANCE) const;

    // Get the pointer to a specific row.
    // Note real data starts at row[myNCL]
    float       *row(int i) const
                 {
                     UT_ASSERT_P(i >= myNRL && i <= myNRH);
                     return &myMatrix[i * myStride];
                 }

    void         clearAndDestroy();

    int          save(ostream &os, int binary) const;

    // I/O friends
    friend ostream      &operator<<(ostream &os, const UT_Matrix &m )
                        {
                            m.outTo(os);
                            return os;
                        }
protected:
    // I/O methods
    void         outTo( ostream &os ) const;

private:
    int          myNRL, myNRH, myNCL, myNCH;
    float       *myMatrix;
    int          myStride, myOffset;
    char         myOwnData;
};

#endif

---