MatrixCharacteristics.h

Go to the documentation of this file.

#ifndef SimTK_SIMMATRIX_MATRIX_CHARACTERISTICS_H_
#define SimTK_SIMMATRIX_MATRIX_CHARACTERISTICS_H_

/* -------------------------------------------------------------------------- *
 *                       Simbody(tm): SimTKcommon                             *
 * -------------------------------------------------------------------------- *
 * This is part of the SimTK biosimulation toolkit originating from           *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody.  *
 *                                                                            *
 * Portions copyright (c) 2005-13 Stanford University and the Authors.        *
 * Authors: Michael Sherman                                                   *
 * Contributors:                                                              *
 *                                                                            *
 * Licensed under the Apache License, Version 2.0 (the "License"); you may    *
 * not use this file except in compliance with the License. You may obtain a  *
 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.         *
 *                                                                            *
 * Unless required by applicable law or agreed to in writing, software        *
 * distributed under the License is distributed on an "AS IS" BASIS,          *
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   *
 * See the License for the specific language governing permissions and        *
 * limitations under the License.                                             *
 * -------------------------------------------------------------------------- */

#include "SimTKcommon/Scalar.h"

#include <iostream>
#include <cassert>
#include <complex>
#include <cstddef>
#include <utility> // for std::pair

namespace SimTK {


class MatrixStructure;
class MatrixStorage;
class MatrixOutline;
class MatrixCondition;
class MatrixCharacter;
class MatrixCommitment;


//  ------------------------------ MatrixStructure -----------------------------
//  ----------------------------------------------------------------------------
class SimTK_SimTKCOMMON_EXPORT MatrixStructure {
public:
    enum Structure {
        NoStructure      = 0x00000000,  
        Matrix1d         = 0x00000001,  
        Zero             = 0x00000002,  
        Identity         = 0x00000004,  
        Permutation      = 0x00000008,  
        RepeatedDiagonal = 0x00000010,  
        Diagonal         = 0x00000020,  
        BiDiagonal       = 0x00000040,  
        TriDiagonal      = 0x00000080,  
        BandedSymmetric  = 0x00000100,  
        BandedHermitian  = 0x00000200,  
        Banded           = 0x00000400,  
        Triangular       = 0x00000800,  
        QuasiTriangular  = 0x00001000,  
        Hessenberg       = 0x00002000,  
        Symmetric        = 0x00004000,  
        Hermitian        = 0x00008000,  
        SkewSymmetric    = 0x00010000,  
        SkewHermitian    = 0x00020000,  
        Full             = 0x00040000   
    };
    static const char* name(Structure);

    typedef unsigned int StructureMask; // 32 bits
    static const StructureMask AnyStructure         = 0x0007ffffU; // see above
    static const StructureMask UncommittedStructure = 0xffffffffU;
    static StructureMask calcStructureMask(Structure);

    enum Position {
        NoPosition = 0x0000,
        Lower      = 0x0001,  // matrix is lower triangular (default)
        Upper      = 0x0002   // matrix is upper triangular
    };
    static const char* name(Position);

    typedef unsigned short PositionMask; // 16 bits
    static const PositionMask AnyPosition         = 0x0003U;  // see above
    static const PositionMask UncommittedPosition = 0xffffU;
    static PositionMask calcPositionMask(Structure);

    enum DiagValue {
        NoDiagValue = 0x0000,
        StoredDiag  = 0x0001, // could be anything (default)
        ZeroDiag    = 0x0002, // zero (e.g. for skew matrices)
        UnitDiag    = 0x0004  // unit (one) diagonal is used frequently by Lapack
    };
    static const char* name(DiagValue);

    typedef unsigned short DiagValueMask; // 16 bits
    static const DiagValueMask AnyDiagValue         = 0x0003U;
    static const DiagValueMask UncommittedDiagValue = 0xffffU;
    static DiagValueMask calcDiagValueMask(Structure);

    MatrixStructure& setMissingAttributes() {
        if (structure == NoStructure)
            structure = Full;
        if (position == NoPosition)
            position = Lower;
        if (diagValue == NoDiagValue)
            diagValue = StoredDiag;
        return *this;
    }

    std::string name() const {
        return std::string(name(getStructure())) 
            + "|" + std::string(name(getPosition()))
            + "|" + std::string(name(getDiagValue()));
    }

    struct Mask {
        Mask() {setToUncommitted();}
        Mask(StructureMask sm, PositionMask pm, DiagValueMask dm)
        :   structure(sm), position(pm), diagValue(dm) {}
        Mask& setToUncommitted() 
        {   structure=UncommittedStructure; position=UncommittedPosition; 
            diagValue=UncommittedDiagValue; return *this; }
        bool isUncommitted() const
        {   return structure==UncommittedStructure && position==UncommittedPosition 
                && diagValue==UncommittedDiagValue; }
        bool isSatisfiedBy(Structure str, Position pos, DiagValue diag) const
        {   return ((StructureMask)str&structure)==(StructureMask)str 
                && ((PositionMask)pos&position)==(PositionMask)pos
                && ((DiagValueMask)diag&diagValue)==(DiagValueMask)diag; }
        bool isSatisfiedBy(const MatrixStructure& actual) const
        {  return isSatisfiedBy(actual.getStructure(), actual.getPosition(), 
                                actual.getDiagValue()); }

        StructureMask  structure;
        PositionMask   position;
        DiagValueMask  diagValue;
    };

    MatrixStructure() {setToNone();}

    MatrixStructure(Structure s, Position p=NoPosition, DiagValue d=NoDiagValue)
        :   structure(s), position(p), diagValue(d) {} 

    Mask mask() const;

    Structure   getStructure() const {return structure;}
    Position    getPosition()  const {return position;}
    DiagValue   getDiagValue() const {return diagValue;}

    MatrixStructure& setStructure(Structure s) {structure=s; return *this;}
    MatrixStructure& setPosition (Position  p) {position=p; return *this;}
    MatrixStructure& setDiagValue(DiagValue d) {diagValue=d; return *this;}

    MatrixStructure& set(Structure s, Position p, DiagValue d)
    {   structure=s; position=p; diagValue=d; return *this; }

    MatrixStructure& setToNone() 
    {   structure=NoStructure; position=NoPosition; 
        diagValue=NoDiagValue; return *this; }

private:
    Structure  structure:32;
    Position   position:16;
    DiagValue  diagValue:16;
};


//  ------------------------------ MatrixStorage -------------------------------
//  ----------------------------------------------------------------------------
class SimTK_SimTKCOMMON_EXPORT MatrixStorage {
public:
    enum Packing {
        NoPacking    = 0x0000,
        Full         = 0x0001,  // full storage layout
        TriInFull    = 0x0002,  // a triangular piece of a full storage layout
        TriPacked    = 0x0004,  // triangle packed into minimal storage, at performance cost
        Banded       = 0x0008,  // a packed, banded storage format
        Vector       = 0x0010,  // a possibly-strided or scattered vector
        Scalar       = 0x0020,  // a single scalar is stored
        Permutation  = 0x0040   // a permuted identity matrix
    };
    static const char* name(Packing);
    typedef unsigned short PackingMask;
    static const PackingMask AllPacking = 0x007fU; // see above
    static const PackingMask UncommittedPacking = 0xffffU;

    enum Placement {
        NoPlacement  = 0x0000,
        Lower        = 0x0001,  // stored in lower triangle of full matrix
        Upper        = 0x0002,  // stored in upper triangle of full matrix
    };
    static const char* name(Placement);
    typedef unsigned short PlacementMask;
    static const PlacementMask AllPlacement = 0x0003U; // see above
    static const PlacementMask UncommittedPlacement = 0xffffU;

    enum Order {
        NoOrder      = 0x0000,
        ColumnOrder  = 0x0001,  // matrix is stored by columns
        RowOrder     = 0x0002,  // matrix is stored by rows
    };
    static const char* name(Order);
    typedef unsigned short OrderMask;
    static const OrderMask AllOrder = 0x03U; // see above
    static const OrderMask UncommittedOrder = 0xffU;

    enum Diagonal {
        NoDiag       = 0x0000,
        StoredDiag   = 0x0001,  // matrix diagonal is stored
        AssumedDiag  = 0x0002   // matrix diagonal is not stored but has known value
    };
    static const char* name(Diagonal);
    typedef unsigned short DiagonalMask;
    static const DiagonalMask AllDiagonal = 0x0003U; // see above
    static const DiagonalMask UncommittedDiagonal = 0xffffU;

    struct Mask {
        Mask()
        :   packing(UncommittedPacking), placement(UncommittedPlacement), 
            order(UncommittedOrder), diagonal(UncommittedDiagonal) {}
        Mask(PackingMask pkm, PlacementMask plm, OrderMask om, DiagonalMask dm)
        :   packing(pkm), placement(plm), order(om), diagonal(dm) {}
        Mask& setToUncommitted()
        {   packing=UncommittedPacking; placement=UncommittedPlacement; 
            order=UncommittedOrder;     diagonal=UncommittedDiagonal; return *this; }
        bool isUncommitted() const 
        {   return packing==UncommittedPacking && placement==UncommittedPlacement 
                && order==UncommittedOrder     && diagonal==UncommittedDiagonal; }        
        bool isSatisfiedBy(Packing pack, Placement place, Order ord, Diagonal diag) const
        {   return ((PackingMask)pack    & packing)   == (PackingMask)  pack 
                && ((PlacementMask)place & placement) == (PlacementMask)place
                && ((OrderMask)ord       & order)     == (OrderMask)    ord     
                && ((DiagonalMask)diag   & diagonal)  == (DiagonalMask) diag; }
        bool isSatisfiedBy(const MatrixStorage& actual) const
        {   return isSatisfiedBy(actual.getPacking(), actual.getPlacement(), 
                                 actual.getOrder(),   actual.getDiagonal());}      

        PackingMask   packing;
        PlacementMask placement;
        OrderMask     order;
        DiagonalMask  diagonal;
    };

    static MatrixStorage calcDefaultStorage(const MatrixStructure&,
                                            const MatrixOutline&);

    std::string name() const {
        return std::string(name(getPacking()))
            + "|" + std::string(name(getPlacement()))
            + "|" + std::string(name(getOrder()))
            + "|" + std::string(name(getDiagonal()));
    }

    Mask mask() const {
        Mask ms; // initially uncommitted
        if (packing)   ms.packing   = (PackingMask)packing;
        if (placement) ms.placement = (PlacementMask)placement;
        if (order)     ms.order     = (OrderMask)order;
        if (diagonal)  ms.diagonal  = (DiagonalMask)diagonal;
        return ms;
    }

    MatrixStorage() 
    :   packing(NoPacking), placement(NoPlacement), order(NoOrder), diagonal(NoDiag) {}

    MatrixStorage(Packing pk, Placement pl=NoPlacement, Order o=NoOrder, Diagonal d=NoDiag)
    :   packing(pk), placement(pl), order(o), diagonal(d) {}

    MatrixStorage(Packing pk, Order o)
    :   packing(pk), placement(NoPlacement), order(o), diagonal(StoredDiag) {}

    MatrixStorage& setMissingAttributes() {
        if (packing==NoPacking) 
            packing = Full;
        if (placement==NoPlacement)
            placement = Lower;
        if (order==NoOrder)
            order = ColumnOrder;
        if (diagonal==NoDiag)
            diagonal = StoredDiag;
        return *this;
    }

    MatrixStorage& setToNone()
    {   packing=NoPacking; placement=NoPlacement; 
        order=NoOrder;     diagonal=NoDiag; return *this; }

    MatrixStorage& setPacking(Packing p)     {packing   = p; return *this;}
    MatrixStorage& setPlacement(Placement p) {placement = p; return *this;}
    MatrixStorage& setOrder(Order o)         {order     = o; return *this;}
    MatrixStorage& setDiagonal(Diagonal d)   {diagonal  = d; return *this;}

    Packing   getPacking()   const {return packing;}
    Placement getPlacement() const {return placement;}
    Order     getOrder()     const {return order;}
    Diagonal  getDiagonal()  const {return diagonal;}

private:
    Packing   packing:16;
    Placement placement:16;
    Order     order:16;
    Diagonal  diagonal:16;
};


//  ------------------------------- MatrixOutline ------------------------------
//  ----------------------------------------------------------------------------
class SimTK_SimTKCOMMON_EXPORT MatrixOutline {
public:
    enum Outline {
        NoOutline   = 0x0000,
        Scalar      = 0x0001,    // 1x1
        Column      = 0x0002,    // mx1, m != 1
        Row         = 0x0004,    // 1xn, n != 1
        Square      = 0x0008,    // mxn, m == n
        Wide        = 0x0010,    // mxn, m < n
        Tall        = 0x0020,    // mxn, m > n
        Rectangular = 0x0040     // mxn
    };
    static const char* name(Outline);

    typedef unsigned short OutlineMask;
    static const OutlineMask AnyOutline  = 0x007fU; // see above
    static const OutlineMask UncommittedOutline = 0xffffU;

    struct Mask {
        Mask() : outline(UncommittedOutline) {}
        explicit Mask(OutlineMask mask) : outline(mask) {}
        Mask& setToUncommitted() {outline=UncommittedOutline; return *this;}
        bool isUncommitted() const {return outline==UncommittedOutline;}
        bool isSatisfiedBy(const MatrixOutline& actual) const 
        {  return ((OutlineMask)actual.outline & outline) == (OutlineMask)actual.outline; }

        OutlineMask outline;
    };

    std::string name() const {return std::string(name(getOutline()));}

    MatrixOutline() : outline(NoOutline) {}

    MatrixOutline(Outline outline) : outline(outline) {}

    MatrixOutline& setToNone() {outline=NoOutline; return *this;}

    static OutlineMask calcMask(Outline);

    Mask mask() const {return Mask(calcMask(getOutline()));}

    bool isSizeOK(int m, int n) const;

    void getMinimumSize(int& m, int& n) const;

    static MatrixOutline calcFromSize(int m, int n);

    Outline getOutline() const {return outline;}

private:
    Outline outline:16;
};



//  ---------------------------- MatrixCondition -------------------------------
//  ----------------------------------------------------------------------------
class SimTK_SimTKCOMMON_EXPORT MatrixCondition {
public:
    enum Condition {
        UnknownCondition = 0x0000,
        Orthogonal       = 0x0001, // implies well conditioned
        PositiveDefinite = 0x0002, // implies well conditioned
        WellConditioned  = 0x0004, // implies full rank
        FullRank         = 0x0008, // but might have bad conditioning
        Singular         = 0x0010  // implies possible bad conditioning 
    };
    static const char* name(Condition);

    typedef unsigned short ConditionMask;   // 16 bits in mask
    static const ConditionMask AnyCondition          = 0x001fU;  // see above
    static const ConditionMask UncommittedCondition  = 0xffffU;

    enum Diagonal {
        UnknownDiagonal   = 0x0000,   
        ZeroDiagonal      = 0x0001,   
        OneDiagonal       = 0x0002,   
        RealDiagonal      = 0x0004,   
        ImaginaryDiagonal = 0x0008    
    };
    static const char* name(Diagonal);

    typedef unsigned short DiagonalMask;   // 16 bits in mask
    static const DiagonalMask AnyDiagonal          = 0x000fU;  // see above
    static const DiagonalMask UncommittedDiagonal  = 0xffffU;

    struct Mask {
        Mask() : condition(UncommittedCondition), diagonal(UncommittedDiagonal) {}
        Mask(ConditionMask cmask, DiagonalMask dmask) : condition(cmask), diagonal(dmask) {}
        Mask& setToUncommitted()    
        {   condition=UncommittedCondition; diagonal=UncommittedDiagonal; return *this;}
        bool isUncommitted() const 
        {   return condition==UncommittedCondition && diagonal==UncommittedDiagonal;}
        bool isSatisfiedBy(const MatrixCondition& actual) const 
        {   return ((ConditionMask)actual.condition & condition) == (ConditionMask)actual.condition
                && ((DiagonalMask) actual.diagonal  & diagonal)  == (DiagonalMask)actual.diagonal; }

        ConditionMask   condition;
        DiagonalMask    diagonal;
    };

    std::string name() const 
    {   return std::string(name(getCondition())) + "|" + std::string(name(getDiagonal()));}

    MatrixCondition() : condition(UnknownCondition), diagonal(UnknownDiagonal) {}

    MatrixCondition(Condition cond, Diagonal diag=UnknownDiagonal) 
    :   condition(cond), diagonal(diag) {}

    MatrixCondition& setToNone() {condition=UnknownCondition; diagonal=UnknownDiagonal; return *this;}

    static ConditionMask calcMask(Condition);

    static DiagonalMask calcMask(Diagonal);

    Mask mask() const 
    {   return Mask(calcMask(getCondition()), calcMask(getDiagonal())); }

    Condition getCondition() const {return condition;}
    Diagonal  getDiagonal()  const {return diagonal;}

    MatrixCondition& setCondition(Condition c) {condition=c; return *this;}
    MatrixCondition& setDiagonal (Diagonal d)  {diagonal=d; return *this;}

private:
    Condition condition:16;
    Diagonal  diagonal:16;
};



//  ------------------------------ MatrixCharacter -----------------------------
class SimTK_SimTKCOMMON_EXPORT MatrixCharacter {
public:
    MatrixCharacter() : nr(0), nc(0), lband(0), uband(0) {}

    // Some handy predefined MatrixCharacters.
    class LapackFull;
    class Vector;
    class RowVector;

    MatrixCharacter& setToNone() {
        nr=nc=lband=uband=0;
        structure.setToNone(); outline.setToNone();
        storage.setToNone();   condition.setToNone();
        return *this;
    }

    int                nrow()       const {return nr;}
    int                ncol()       const {return nc;}
    std::pair<int,int> getSize()    const {return std::pair<int,int>(nrow(),ncol());}
    ptrdiff_t          nelt()       const {return (ptrdiff_t)nrow() * (ptrdiff_t)ncol();}

    int                getLowerBandwidth() const {return lband;}
    int                getUpperBandwidth() const {return uband;}
    std::pair<int,int> getBandwidth()      const 
    {   return std::pair<int,int>(getLowerBandwidth(), getUpperBandwidth()); }

    const MatrixStructure&  getStructure() const {return structure;}
    const MatrixStorage&    getStorage()   const {return storage;}
    const MatrixOutline&    getOutline()   const {return outline;}
    const MatrixCondition&  getCondition() const {return condition;}

    MatrixStructure&  updStructure() {return structure;}
    MatrixStorage&    updStorage()   {return storage;}
    MatrixOutline&    updOutline()   {return outline;}
    MatrixCondition&  updCondition() {return condition;}

    MatrixCharacter& setStructure(const MatrixStructure& sa)  {structure = sa; return *this;}
    MatrixCharacter& setStorage  (const MatrixStorage&   sa)  {storage   = sa; return *this;}
    MatrixCharacter& setOutline  (const MatrixOutline&   oa)  {outline   = oa; return *this;}
    MatrixCharacter& setCondition(const MatrixCondition& ca)  {condition = ca; return *this;}


    MatrixCharacter& setActualSize(int m, int n)
    {   setSize(m,n); outline = MatrixOutline::calcFromSize(m,n); return *this; }
    MatrixCharacter& setActualNumRows(int m)
    {   setNumRows(m); outline = MatrixOutline::calcFromSize(m,ncol()); return *this; }
    MatrixCharacter& setActualNumCols(int n)
    {   setNumCols(n); outline = MatrixOutline::calcFromSize(nrow(),n); return *this; }

    MatrixCharacter& setBandwidth(int lb, int ub) {
        assert(lb>=0 && lb>=0);
        lband = lb; uband = ub;
        return *this;
    }
    MatrixCharacter& setLowerBandwidth(int lb) {
        assert(lb>=0);
        lband = lb;
        return *this;
    }
    MatrixCharacter& setUpperBandwidth(int ub) {
        assert(ub>=0);
        uband = ub;
        return *this;
    }

    class Mask; // defined below

protected:
    MatrixCharacter(int m, int n,
                    int lb, int ub,
                    MatrixStructure structure,
                    MatrixStorage   storage,
                    MatrixCondition condition)
    :   nr(m), nc(n), lband(lb), uband(ub),
        structure(structure), storage(storage), 
        outline(MatrixOutline::calcFromSize(m,n)), 
        condition(condition) {}


    int              nr,            
                     nc;            
    int              lband,         
                     uband;         
    MatrixStructure  structure;
    MatrixStorage    storage;
    MatrixOutline    outline;
    MatrixCondition  condition;

private:
    // These are private because they don't set the outline as well.
    MatrixCharacter& setSize(int m, int n) 
    {   assert(m>=0 && n>=0); nr = m; nc = n; return *this; }
    MatrixCharacter& setNumRows(int m) 
    {   assert(m>=0); nr = m; return *this; }
    MatrixCharacter& setNumCols(int n) 
    {   assert(n>=0); nc = n; return *this; }
};

SimTK_SimTKCOMMON_EXPORT std::ostream& 
operator<<(std::ostream& o, const MatrixCharacter&);

class MatrixCharacter::LapackFull : public MatrixCharacter {
public:
    LapackFull(int m, int n)
    :   MatrixCharacter(m,n,0,0,
            MatrixStructure(MatrixStructure::Full),
            MatrixStorage(MatrixStorage::Full,MatrixStorage::ColumnOrder),
            MatrixCondition()) {}                   
};

class MatrixCharacter::Vector : public MatrixCharacter {
public:
    Vector(int m)
    :   MatrixCharacter(m,1,0,0,
            MatrixStructure(MatrixStructure::Matrix1d),
            MatrixStorage(MatrixStorage::Vector,MatrixStorage::ColumnOrder),
            MatrixCondition()) {}                   
};

class MatrixCharacter::RowVector : public MatrixCharacter {
public:
    RowVector(int n)
    :   MatrixCharacter(1,n,0,0,
            MatrixStructure(MatrixStructure::Matrix1d),
            MatrixStorage(MatrixStorage::Vector,MatrixStorage::RowOrder),
            MatrixCondition()) {}                   
};

//  -------------------------- MatrixCharacter::Mask ---------------------------
//  ----------------------------------------------------------------------------
class MatrixCharacter::Mask {
public:
    Mask() {setToUncommitted();}

    typedef unsigned int SizeMask;
    static const SizeMask SizeUncommitted = 0xffffffffU;

    bool isResizeable()      const {return nr==SizeUncommitted || nc==SizeUncommitted;}
    bool isFullyResizeable() const {return nr==SizeUncommitted && nc==SizeUncommitted;}
    bool isNumRowsLocked()   const {return nr!=SizeUncommitted;}
    bool isNumColsLocked()   const {return nc!=SizeUncommitted;}

    unsigned int getNumRowsMask() const {return nr;}
    unsigned int getNumColsMask() const {return nc;}
    unsigned int getLowerBandwidthMask() const {return lband;}
    unsigned int getUpperBandwidthMask() const {return uband;}

    int getDefaultNumRows() const {return isNumRowsLocked() ? nr : 0;}
    int getDefaultNumCols() const {return isNumColsLocked() ? nc : 0;}

    bool isLowerBandwidthLocked()  const {return lband!=SizeUncommitted;}
    bool isUpperBandwidthLocked()  const {return uband!=SizeUncommitted;}
    int getDefaultLowerBandwidth() const {return isLowerBandwidthLocked() ? lband : 0;}
    int getDefaultUpperBandwidth() const {return isUpperBandwidthLocked() ? uband : 0;}

    Mask& setToUncommitted() {
        nr=nc=lband=uband=SizeUncommitted;
        structure.setToUncommitted(); storage.setToUncommitted();
        outline.setToUncommitted();   condition.setToUncommitted();
        return *this;
    }

    bool isUncommitted() const {
        return nr==SizeUncommitted       && nc==SizeUncommitted 
            && lband==SizeUncommitted    && uband==SizeUncommitted
            && structure.isUncommitted() && storage.isUncommitted()
            && outline.isUncommitted()   && condition.isUncommitted();
    }

    bool isSatisfiedBy(const MatrixCharacter& actual) const {
        return isSizeOK(actual.nr, actual.nc) 
            && isBandwidthOK(actual.lband, actual.uband)
            && structure.isSatisfiedBy(actual.getStructure())
            && storage.isSatisfiedBy(actual.getStorage())
            && outline.isSatisfiedBy(actual.getOutline())
            && condition.isSatisfiedBy(actual.getCondition());
    }

    bool isSizeOK(int m, int n) const 
    {   return ((SizeMask)m & nr)      == (SizeMask)m
            && ((SizeMask)n & nc)      == (SizeMask)n; }

    bool isBandwidthOK(int lower, int upper) const 
    {   return ((SizeMask)lower & lband) == (SizeMask)lower
            && ((SizeMask)upper & uband) == (SizeMask)upper; }

    SizeMask                nr,         
                            nc;         
    SizeMask                lband,      
                            uband;      
    MatrixStructure::Mask   structure;
    MatrixStorage::Mask     storage;
    MatrixOutline::Mask     outline;
    MatrixCondition::Mask   condition;

friend class MatrixCommitment;
};

//  ----------------------------- MatrixCommitment -----------------------------


//  ----------------------------------------------------------------------------
class SimTK_SimTKCOMMON_EXPORT MatrixCommitment {
public:
    MatrixCommitment() {} // set commitments to "none" and masks to "uncommitted"

    MatrixCommitment(const MatrixStructure& str)
    { new (this) MatrixCommitment(str, MatrixStorage(), MatrixOutline(), MatrixCondition());}

    class Vector;
    class RowVector;
    class Triangular;
    class Symmetric;
    class Hermitian;
    class SkewSymmetric;
    class SkewHermitian;

    MatrixCommitment& commitSize(int m, int n) 
    {   commitNumRows(m); commitNumCols(n); return *this; }
    MatrixCommitment& commitNumRows(int m) 
    {   SimTK_SIZECHECK_NONNEG(m, "MatrixCommitment::commitNumRows()");
        masks.nr = m; return *this; }
    MatrixCommitment& commitNumCols(int n)  
    {   SimTK_SIZECHECK_NONNEG(n, "MatrixCommitment::commitNumCols()");
        masks.nc = n; return *this; }

    MatrixCommitment& commitBandwidth(int lb, int ub) 
    {  commitLowerBandwidth(lb); commitUpperBandwidth(ub); return *this;}
    MatrixCommitment& commitLowerBandwidth(int lb)
    {   SimTK_SIZECHECK_NONNEG(lb, "MatrixCommitment::commitLowerBandwidth()");
        masks.lband = lb; return *this; }
    MatrixCommitment& commitUpperBandwidth(int ub)
    {   SimTK_SIZECHECK_NONNEG(ub, "MatrixCommitment::commitUpperBandwidth()");
        masks.uband = ub; return *this; }

    MatrixCommitment& commitStructure(const MatrixStructure& s) 
    {   structure=s; masks.structure=s.mask(); return *this; }
    MatrixCommitment& commitStorage  (const MatrixStorage&   s) 
    {   storage=s;   masks.storage  =s.mask(); return *this; }
    MatrixCommitment& commitOutline  (const MatrixOutline&   o) 
    {   outline=o;   masks.outline  =o.mask(); return *this; }
    MatrixCommitment& commitCondition(const MatrixCondition& c) 
    {   condition=c; masks.condition=c.mask(); return *this; }

    MatrixCharacter calcDefaultCharacter(int minNumRows, int minNumCols) const;

    const MatrixStructure&  getStructureCommitment() const {return structure;}
    const MatrixStorage&    getStorageCommitment()   const {return storage;}
    const MatrixOutline&    getOutlineCommitment()   const {return outline;}
    const MatrixCondition&  getConditionCommitment() const {return condition;}

    const MatrixStructure::Mask&   getStructureMask() const {return masks.structure;}
    const MatrixStorage::Mask&     getStorageMask()   const {return masks.storage;}
    const MatrixOutline::Mask&     getOutlineMask()   const {return masks.outline;}
    const MatrixCondition::Mask&   getConditionMask() const {return masks.condition;}

    MatrixCharacter::Mask::SizeMask getNumRowsMask() const {return masks.nr;}
    MatrixCharacter::Mask::SizeMask getNumColsMask() const {return masks.nc;}
    MatrixCharacter::Mask::SizeMask getLowerBandwidthMask() const {return masks.lband;}
    MatrixCharacter::Mask::SizeMask getUpperBandwidthMask() const {return masks.uband;}

    int getDefaultNumRows() const {return masks.getDefaultNumRows();}
    int getDefaultNumCols() const {return masks.getDefaultNumCols();}

    bool isSizeOK(int m, int n) const {return masks.isSizeOK(m,n);} 
    bool isSizeOK(const std::pair<int,int>& mn) const
    {   return isSizeOK(mn.first, mn.second); }

    bool isBandwidthOK(int lower, int upper) const {return masks.isBandwidthOK(lower,upper);} 

    bool isSatisfiedBy(const MatrixCharacter& actual) const 
    {   return masks.isSatisfiedBy(actual); }
    bool isStructureOK(const MatrixStructure& s) const
    {   return getStructureMask().isSatisfiedBy(s); }
    bool isStorageOK(const MatrixStorage& s) const
    {   return getStorageMask().isSatisfiedBy(s); }
    bool isOutlineOK(const MatrixOutline& o) const
    {   return getOutlineMask().isSatisfiedBy(o); }
    bool isConditionOK(const MatrixCondition& c) const
    {   return getConditionMask().isSatisfiedBy(c); }

    bool isResizeable()      const {return masks.isResizeable();}
    bool isFullyResizeable() const {return masks.isFullyResizeable();;}
    bool isNumRowsLocked()  const {return masks.isNumRowsLocked();}
    bool isNumColsLocked()  const {return masks.isNumColsLocked();}

    bool isStructureCommitted() const 
    {   return !getStructureMask().isUncommitted(); }
    bool isStorageCommitted()   const 
    {   return !getStorageMask().isUncommitted();}
    bool isOutlineCommitted()   const 
    {   return !getOutlineMask().isUncommitted(); }
    bool isConditionCommitted() const 
    {   return !getConditionMask().isUncommitted();}

    void clear() {
        structure.setToNone(); 
        storage.setToNone(); 
        outline.setToNone(); 
        condition.setToNone();
        masks.setToUncommitted();
    }

protected:
    MatrixCommitment(const MatrixStructure& structure,
                     const MatrixStorage&   storage,                    
                     const MatrixOutline&   outline,
                     const MatrixCondition& condition)
    :   structure(structure), storage(storage), 
        outline(outline), condition(condition),
        masks() // set to all 1's 
    {
        if (outline.getOutline()==MatrixOutline::Scalar) commitSize(1,1);
        else if (outline.getOutline()==MatrixOutline::Column) commitNumCols(1);
        else if (outline.getOutline()==MatrixOutline::Row) commitNumRows(1);

        masks.structure = structure.mask();
        masks.storage   = storage.mask();
        masks.outline   = outline.mask();
        masks.condition = condition.mask(); 
    }
    
    MatrixStructure         structure;
    MatrixStorage           storage;
    MatrixOutline           outline;
    MatrixCondition         condition;

    MatrixCharacter::Mask   masks;
};


class MatrixCommitment::Vector : public MatrixCommitment {
public:
    Vector()
    :   MatrixCommitment
        (   MatrixStructure(MatrixStructure::Matrix1d), 
            MatrixStorage(),
            MatrixOutline(MatrixOutline::Column), 
            MatrixCondition())
    {
    }
    explicit Vector(int m)
    :   MatrixCommitment
        (   MatrixStructure(MatrixStructure::Matrix1d), 
            MatrixStorage(),
            MatrixOutline(MatrixOutline::Column), 
            MatrixCondition())
    {
        commitNumRows(m);
    }
};

class MatrixCommitment::RowVector : public MatrixCommitment {
public:
    RowVector()
    :   MatrixCommitment
        (   MatrixStructure(MatrixStructure::Matrix1d), 
            MatrixStorage(),
            MatrixOutline(MatrixOutline::Row), 
            MatrixCondition())
    {
    }
    explicit RowVector(int n)
    :   MatrixCommitment
        (   MatrixStructure(MatrixStructure::Matrix1d), 
            MatrixStorage(),
            MatrixOutline(MatrixOutline::Row), 
            MatrixCondition())
    {
        commitNumCols(n);
    }
};

class MatrixCommitment::Triangular : public MatrixCommitment {
public:
    Triangular()
    :   MatrixCommitment(MatrixStructure::Triangular, MatrixStorage(),
                         MatrixOutline(), MatrixCondition())
    {
    }
};

class MatrixCommitment::Symmetric : public MatrixCommitment {
public:
    Symmetric()
    :   MatrixCommitment(MatrixStructure::Symmetric, MatrixStorage(),
                         MatrixOutline(), MatrixCondition())
    {
    }
};

class MatrixCommitment::Hermitian : public MatrixCommitment {
public:
    Hermitian()
    :   MatrixCommitment
        (   MatrixStructure::Hermitian, 
            MatrixStorage(),
            MatrixOutline(), 
            MatrixCondition().setDiagonal(MatrixCondition::RealDiagonal))
    {
    }
};

class MatrixCommitment::SkewSymmetric : public MatrixCommitment {
public:
    SkewSymmetric()
    :   MatrixCommitment
        (   MatrixStructure::SkewSymmetric, 
            MatrixStorage(),
            MatrixOutline(), 
            MatrixCondition().setDiagonal(MatrixCondition::ZeroDiagonal))
    {
    }
};

class MatrixCommitment::SkewHermitian : public MatrixCommitment {
public:
    SkewHermitian()
    :   MatrixCommitment
        (   MatrixStructure::SkewHermitian, 
            MatrixStorage(),
            MatrixOutline(), 
            MatrixCondition().setDiagonal(MatrixCondition::ImaginaryDiagonal))
    {
    }
};

SimTK_SimTKCOMMON_EXPORT std::ostream& 
operator<<(std::ostream& o, const MatrixCommitment&);
     
} //namespace SimTK

#endif // SimTK_SIMMATRIX_MATRIX_CHARACTERISTICS_H_