Branch data Line data Source code
1 : : //- Class: CubitMatrix
2 : : //-
3 : : //- Description: This file defines the CubitMatrix class which is a
4 : : //- standard NxM Matrix. All relevant arithmetic operators are
5 : : //- overloaded so CubitMatrix's can be used similar to built-in types.
6 : : //-
7 : : //- Owner: Dan Goodrich
8 : : //- Checked by:
9 : :
10 : :
11 : : #ifndef CUBITMATRIX_HPP
12 : : #define CUBITMATRIX_HPP
13 : :
14 : : #include "CubitVector.hpp"
15 : : #include <vector>
16 : : #include <stdexcept>
17 : : #include <string>
18 : : #include "CGMUtilConfigure.h"
19 : :
20 : : class CUBIT_UTIL_EXPORT CubitMatrix
21 : : {
22 : : public:
23 : :
24 : : //- Heading: Constructors and Destructor
25 : : CubitMatrix();
26 : : //- Default constructor. 3x3 init to zeros.
27 : :
28 : : CubitMatrix( const int num_rows, const int num_cols);
29 : : //- Initializes NxM matrix to zeros.
30 : :
31 : : CubitMatrix( const CubitVector& vec1, const CubitVector& vec2,
32 : : const CubitVector& vec3);
33 : : //- Constructor: create 3x3 Matrix from three vector components
34 : : // where each vector forms one column of the matrix.
35 : :
36 : : CubitMatrix( const CubitVector& vec1, const CubitVector& vec2,
37 : : const CubitVector& vec3, const CubitVector& vec4 );
38 : : //- Constructor: create 3x4 Matrix from four vector components
39 : : // where each vector forms one column of the matrix.
40 : :
41 : : CubitMatrix( const CubitVector& vec1, const CubitVector& vec2 );
42 : : //- Constructor: create 3x3 Matrix from two vectors (otimes)
43 : : // ( vec1 otimes vec2 )_ij = vec1_i * vec2_j
44 : :
45 : : CubitMatrix( const int n );
46 : : //- Constructor: create n x n Identity matrix
47 : :
48 : : CubitMatrix( std::vector<int> &is, std::vector<int> &js, std::vector<double> &es, int n, int m );
49 : : // - uses vectors i, j, and s to generate an m-by-n
50 : : // - matrix, S, such that S(i(k),j(k)) = s(k).
51 : : // - Vectors i, j, and s are all the same length.
52 : : // - Any elements of s that have duplicate values of i and j are added
53 : : // - together.
54 : :
55 : : CubitMatrix( const CubitMatrix& matrix);
56 : : //- Copy Constructor
57 : :
58 : : virtual ~CubitMatrix();
59 : : //- destructor
60 : :
61 : : void reset_size( const int n, const int m, double default_value = 0.0 );
62 : : //- resize this matrix and set all values to a default value.
63 : :
64 : : //- Heading: Set and Inquire Functions
65 : 2802 : int num_rows() const { return numRows; }
66 : 8962 : int num_cols() const { return numCols; }
67 : :
68 : : void print_matrix() const;
69 : : void print_matrix( char *filename ) const;
70 : : //- Prints matrix
71 : :
72 : : void set_to_identity();
73 : :
74 : : bool is_identity( double tol = 1e-8 ) const;
75 : : bool is_equal( const CubitMatrix &other, double tol = 1e-8 ) const;
76 : :
77 : : //- Add an identity matrix into this matrix.
78 : : void plus_identity();
79 : :
80 : : //- Change Matrix component row n, column m to val.
81 : 12280 : void set(const int n, const int m, const double val)
82 : : {
83 [ + - ][ + - ]: 12280 : if(n < 0 || n >= numRows || m < 0 || m >= numCols)
[ + - ][ - + ]
84 [ # # ][ # # ]: 0 : throw std::invalid_argument (std::string("Index out of bounds"));
85 : : //assert (n >= 0 && n < numRows);
86 : : //assert (m >= 0 && m < numCols);
87 : 12280 : matrixPtr[n][m] = val;
88 : 12280 : }
89 : :
90 : : //- add the value to the component row n col m
91 : 352 : void add(const int n, const int m, const double val)
92 : : {
93 : 352 : matrixPtr[n][m] += val;
94 : 352 : }
95 : :
96 : : //- Gets the values of the matrix at position (n,m)
97 : 24243 : double get( int n, int m ) const
98 : : {
99 [ + - ][ + - ]: 24243 : if(n < 0 || n >= numRows || m < 0 || m >= numCols)
[ + - ][ - + ]
100 [ # # ][ # # ]: 0 : throw std::invalid_argument ("Index out of Bounds");
101 : : //assert (n >= 0 && n < numRows);
102 : : //assert (m >= 0 && m < numCols);
103 : 24243 : return matrixPtr[n][m];
104 : : }
105 : :
106 : : //- Heading: Matrix Algebra functions
107 : : // (assert if incompatible element sizes)
108 : : CubitMatrix operator= (const CubitMatrix& matrix);
109 : : CubitMatrix operator* (const CubitMatrix& matrix) const;
110 : : CubitVector operator* (const CubitVector& vector) const;
111 : : std::vector<double> operator* (const std::vector<double> & vector) const;
112 : : CubitMatrix operator* (double val ) const;
113 : : CubitMatrix operator/ (double val ) const;
114 : : CubitMatrix operator+ (const CubitMatrix& matrix) const;
115 : : CubitMatrix operator- (const CubitMatrix& matrix) const;
116 : :
117 : : CubitMatrix& operator+=(const CubitMatrix &matrix);
118 : : CubitMatrix& operator*=(const double multiplier);
119 : : //- Scale the matrix by a linear factor: {this = this * multiplier},
120 : :
121 : : CubitMatrix inverse(); // inverts 1x1, 2x2, or 3x3 matrix
122 : : // asserts if singular
123 : : CubitMatrix symm() const; // returns matrix^transpose * matrix
124 : :
125 : : CubitBoolean positive_definite() const; // returns true if matrix is positive definite
126 : :
127 : : double determinant () const;
128 : : double inf_norm () const; // infinity norm
129 : : double frobenius_norm_squared() const; // square of frobenius norm of A
130 : : double frobenius_norm_squared_symm() const; // square of frobenius norm of A^t A
131 : : double frobenius_norm_squared_adj() const; // square of frobenius norm of adjoint A
132 : : double frobenius_norm_squared_inv() const; // square of frobenius norm of A-inverse
133 : : double condition() const; // condition number of A using frobenius norm
134 : : int rank( void ) const;
135 : :
136 : : double cofactor (const int row, const int col) const;
137 : : CubitMatrix adjoint() const;
138 : : CubitMatrix transpose() const;
139 : :
140 : : // row and col indicate the row and column to leave out of the
141 : : // resulting matrix.
142 : : CubitMatrix sub_matrix( const int row, const int col ) const;
143 : :
144 : : // Create a matrix containing the rows and cols of this that are true in
145 : : // rows_to_include and cols_to_include.
146 : : void sub_matrix( const std::vector<bool> &rows_to_include,
147 : : const std::vector<bool> &cols_to_include,
148 : : CubitMatrix &submatrix );
149 : :
150 : : // routines to perform Gaussian elimination with pivoting and scaling
151 : : // on 3x3 matrix
152 : : int gauss_elim (CubitVector &b);
153 : : int factor (CubitVector &pivot);
154 : : void solve (CubitVector &b, const CubitVector& pivot);
155 : :
156 : : // routines to solve NxN system (from Numerical Recipes in C)
157 : : CubitStatus solveNxN( CubitMatrix& rhs, // must be NxM, where M is any
158 : : // number of colums.
159 : : CubitMatrix& coef ); // must be NxM, same as rhs.
160 : : CubitStatus solveNxN( const std::vector<double> &rhs,
161 : : std::vector<double> &coef );
162 : : CubitStatus ludcmp( double *indx, double& d );
163 : : CubitStatus lubksb( double *indx, double *b );
164 : :
165 : : // matrix must have been sized to be 3x3 before calling.
166 : : bool vector_outer_product(const CubitVector& vec1,
167 : : const CubitVector& vec2 );
168 : :
169 : : private:
170 : : double **matrixPtr;
171 : : double *matrixMem;
172 : : int numRows;
173 : : int numCols;
174 : : };
175 : :
176 : : #endif
177 : :
178 : :
179 : :
|
