1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
//- Class: CubitMatrix
//-
//- Description: This file defines the CubitMatrix class which is a
//- standard NxM Matrix. All relevant arithmetic operators are
//- overloaded so CubitMatrix's can be used similar to built-in types.
//-
//- Owner: Dan Goodrich
//- Checked by:


#ifndef CUBITMATRIX_HPP
#define CUBITMATRIX_HPP

#include "CubitVector.hpp"
#include <vector>
#include <stdexcept>
#include <string>
#include "CubitUtilConfigure.h"

class CUBIT_UTIL_EXPORT CubitMatrix
{
public:
  
    //- Heading: Constructors and Destructor
  CubitMatrix();
    //- Default constructor. 3x3 init to zeros.
  
  CubitMatrix( const int num_rows, const int num_cols);
    //- Initializes NxM matrix to zeros.
  
  CubitMatrix( const CubitVector& vec1, const CubitVector& vec2,
               const CubitVector& vec3);
    //- Constructor: create 3x3 Matrix from three vector components
    //  where each vector forms one column of the matrix.
  
  CubitMatrix( const CubitVector& vec1, const CubitVector& vec2,
               const CubitVector& vec3, const CubitVector& vec4 );
    //- Constructor: create 3x4 Matrix from four vector components
    //  where each vector forms one column of the matrix.

  CubitMatrix( const CubitVector& vec1, const CubitVector& vec2 );
    //- Constructor: create 3x3 Matrix from two vectors (otimes)
    //  ( vec1 otimes vec2 )_ij = vec1_i * vec2_j
  
  CubitMatrix( const int n );
    //- Constructor: create n x n Identity matrix

  CubitMatrix( std::vector<int> &is, std::vector<int> &js, std::vector<double> &es, int n, int m );
    // - uses vectors i, j, and s to generate an m-by-n
    // - matrix, S,  such that S(i(k),j(k)) = s(k).
    // - Vectors i, j, and s are all the same length.
    // - Any elements of s that have duplicate values of i and j are added
    // - together.

  CubitMatrix( const CubitMatrix& matrix);
    //- Copy Constructor
  
  virtual ~CubitMatrix();
    //- destructor
  
  void reset_size( const int n, const int m, double default_value = 0.0 );
    //- resize this matrix and set all values to a default value.

    //- Heading: Set and Inquire Functions
  int num_rows() const { return numRows; }
  int num_cols() const { return numCols; }
  
  void print_matrix() const;
  void print_matrix( char *filename ) const;
    //- Prints matrix

  void set_to_identity();
  
  bool is_identity( double tol = 1e-8 ) const;
  bool is_equal( const CubitMatrix &other, double tol = 1e-8 ) const;

   //- Add an identity matrix into this matrix.
  void plus_identity();
  
    //- Change Matrix component row n, column m to val.
  void set(const int n, const int m, const double val)
    {
      if(n < 0 || n >= numRows || m < 0 || m >= numCols)
        throw std::invalid_argument (std::string("Index out of bounds"));
      //assert (n >= 0 && n < numRows);
      //assert (m >= 0 && m < numCols);
      matrixPtr[n][m] = val;
    }

    //- add the value to the component row n col m
  void add(const int n, const int m, const double val)
    {
      matrixPtr[n][m] += val;
    }
  
    //- Gets the values of the matrix at position (n,m)
  double get( int n, int m ) const
    {
      if(n < 0 || n >= numRows || m < 0 || m >= numCols)
        throw std::invalid_argument ("Index out of Bounds");
      //assert (n >= 0 && n < numRows);
      //assert (m >= 0 && m < numCols);
      return matrixPtr[n][m];
    }
  
    //- Heading: Matrix Algebra functions 
    //  (assert if incompatible element sizes)
  CubitMatrix operator= (const CubitMatrix& matrix);<--- 'CubitMatrix::operator=' should return 'CubitMatrix &'.
  CubitMatrix operator* (const CubitMatrix& matrix) const;
  CubitVector operator* (const CubitVector& vector) const;
  std::vector<double> operator* (const std::vector<double> & vector) const;
  CubitMatrix operator* (double val ) const;
  CubitMatrix operator/ (double val ) const;
  CubitMatrix operator+ (const CubitMatrix& matrix) const;
  CubitMatrix operator- (const CubitMatrix& matrix) const;
  
  CubitMatrix& operator+=(const CubitMatrix &matrix);
  CubitMatrix& operator*=(const double multiplier);
    //- Scale the matrix by a linear factor: {this = this * multiplier},
  
  CubitMatrix inverse(); // inverts 1x1, 2x2, or 3x3 matrix
    // asserts if singular
  CubitMatrix symm() const; // returns matrix^transpose * matrix
 
  CubitBoolean positive_definite() const;  // returns true if matrix is positive definite
  
  double determinant () const;
  double inf_norm () const; // infinity norm
  double frobenius_norm_squared() const; // square of frobenius norm of A
  double frobenius_norm_squared_symm() const; // square of frobenius norm of A^t A
  double frobenius_norm_squared_adj() const; // square of frobenius norm of adjoint A
  double frobenius_norm_squared_inv() const; // square of frobenius norm of A-inverse
  double condition() const;  // condition number of A using frobenius norm
  int rank( void ) const;

  double cofactor (const int row, const int col) const;
  CubitMatrix adjoint() const;
  CubitMatrix transpose() const;

     // row and col indicate the row and column to leave out of the
    // resulting matrix. 
  CubitMatrix sub_matrix( const int row, const int col ) const;

  // Create a matrix containing the rows and cols of this that are true in
  // rows_to_include and cols_to_include.
  void sub_matrix( const std::vector<bool> &rows_to_include,
                   const std::vector<bool> &cols_to_include,
                   CubitMatrix &submatrix );
  
    // routines to perform Gaussian elimination with pivoting and scaling
    // on 3x3 matrix
  int  gauss_elim (CubitVector &b);
  int  factor (CubitVector &pivot);
  void solve (CubitVector &b, const CubitVector& pivot);

    // routines to solve NxN system (from Numerical Recipes in C)
  CubitStatus solveNxN( CubitMatrix& rhs,  // must be NxM, where M is any
                                           // number of colums.
                        CubitMatrix& coef ); // must be NxM, same as rhs.
  CubitStatus solveNxN( const std::vector<double> &rhs,
                        std::vector<double> &coef );
  CubitStatus ludcmp( double *indx, double& d );
  CubitStatus lubksb( double *indx, double *b );

  // matrix must have been sized to be 3x3 before calling.
  bool vector_outer_product(const CubitVector& vec1,
                            const CubitVector& vec2 );
  
private:
  double **matrixPtr;
  double *matrixMem;
  int numRows;
  int numCols;
};

#endif